My Project
programmer's documentation
Infinite channel inlet with near-inlet feedback

A second method of defining an inlet condition which converges towards an infinite channel profile is based simply on feedback from values computed at inlet cell centers (combining the inlet boundary conditions and the effect of wall friction on the inlet walls). It assumes the mesh is very regular and orthogonal, at least on the first two inlet cell layers (using a gradient correction in the case of a less regular mesh might work, but has never been tested.

Local variables to be added

The following local variables need to be defined for the examples in this section:

integer ifac, iel, ii, ilelt, nlelt
double precision xustar2, xdh, d2s3, rhomoy
double precision acc(2), fmprsc, fmul, uref2, vnrm
integer, allocatable, dimension(:) :: lstelt, mrkcel
double precision, dimension(:), pointer :: bfpro_rom
double precision, dimension(:), pointer :: cvar_r11, cvar_r22, cvar_r33
double precision, dimension(:), pointer :: cvar_r12, cvar_r23, cvar_r13
double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_phi
double precision, dimension(:), pointer :: cvar_omg
double precision, dimension(:), pointer :: cvar_al, cvar_fb
double precision, dimension(:), pointer :: cvar_nusa
double precision, dimension(:), pointer :: cvar_scal
double precision, dimension(:,:), pointer :: cvar_vel
double precision, dimension(:,:), pointer :: cvar_rij

Initialization and finalization

Initialization and finalization is similar to that of the base examples

Body

Here, we define an inlet boundary condition for a very long channel or duct with a section matching the boundary faces of group 'INLET'.

We fix a mean inlet velocity, and use a feedback loop assuming a fixed-point type behavior will allow us to reach a state matching that of a very long inlet channel.

Warning
We assume other boundary conditions are defined before this one (ideally, using the GUI).
We also assume that the mesh is orthogonal at the inlet, and we are using a RANS (not LES) computation. to the current inlet.

For EBRSM of V2f models, to avoid laminarization at the inlet, the initial velocity (at the first time step) is divided by 10 on inlet faces adjacent to the boundary, so as to ensure a velocity gradient and turbulent production. Otherwise, the initialization may fail.

call getfbr('INLET', nlelt, lstelt)
!==========
fmprsc = 1.d0 ! mean prescribed velocity
if (ntcabs.eq.1) then
! For the Rij-EBRSM model (and possibly V2f), we need a non-flat profile,
! so as to ensure turbulent production, and avoid laminarization;
! here, we simply divide the initial velocity by 10 for inlet
! faces adjacent to the wall.
! The loop below assumes wall conditions have been defined first
! (in the GUI, or in this file, before the current test).
if (iturb.eq.32 .or. itytur.eq.5) then
allocate(mrkcel(ncelet))
do iel = 1, ncelet
mrkcel(iel) = 0
enddo
do ifac = 1, nfabor
if (itypfb(ifac) .eq. iparoi) then
iel = ifabor(ifac)
mrkcel(iel) = 1
endif
enddo
endif
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
iel = ifabor(ifac)
itypfb(ifac) = ientre
rcodcl(ifac,iu,1) = - fmprsc * surfbo(1,ifac) / surfbn(ifac)
rcodcl(ifac,iv,1) = - fmprsc * surfbo(2,ifac) / surfbn(ifac)
rcodcl(ifac,iw,1) = - fmprsc * surfbo(3,ifac) / surfbn(ifac)
if (iturb.eq.32 .or. itytur.eq.5) then
if (mrkcel(iel) .eq. 1) then
rcodcl(ifac,iu,1) = fmprsc/10.d0
endif
endif
uref2 = rcodcl(ifac,iu,1)**2 &
+ rcodcl(ifac,iv,1)**2 &
+ rcodcl(ifac,iw,1)**2
uref2 = max(uref2,1.d-12)
! Turbulence example computed using equations valid for a pipe.
! We will be careful to specify a hydraulic diameter adapted
! to the current inlet.
! We will also be careful if necessary to use a more precise
! formula for the dynamic viscosity use in the calculation of
! the Reynolds number (especially if it is variable, it may be
! useful to take the law from 'usphyv'. Here, we use by default
! the 'viscl0" value.
! Regarding the density, we have access to its value at boundary
! faces (romb) so this value is the one used here (specifically,
! it is consistent with the processing in 'usphyv', in case of
! variable density)
! Hydraulic diameter
xdh = 1.d0
! Calculation of turbulent inlet conditions using
! the turbulence intensity and standard laws for a circular pipe
! (their initialization is not needed here but is good practice)
rhomoy = bfpro_rom(ifac)
call turbulence_bc_inlet_hyd_diam(ifac, uref2, xdh, rhomoy, viscl0, &
rcodcl)
! Handle scalars
if (nscal.gt.0) then
do ii = 1, nscal
rcodcl(ifac,isca(ii),1) = 1.d0
enddo
endif
enddo
if (iturb.eq.32 .or. itytur.eq.5) then
deallocate(mrkcel)
endif
else
! Subsequent time steps
!----------------------
acc(1) = 0.d0
acc(2) = 0.d0
! Estimate multiplier
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
iel = ifabor(ifac)
vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2)
acc(1) = acc(1) + vnrm*surfbn(ifac)
acc(2) = acc(2) + surfbn(ifac)
enddo
if (irangp.ge.0) then
call parrsm(2, acc)
endif
if (acc(1).le.epzero) then
fmul = 0.d0 ! zero velocity in bulk domain
else
fmul = fmprsc/(acc(1)/acc(2)) ! 1 / estimate flow multiplier
endif
! Apply BC
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
iel = ifabor(ifac)
itypfb(ifac) = ientre
vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2)
rcodcl(ifac,iu,1) = - fmul * vnrm * surfbo(1,ifac) / surfbn(ifac)
rcodcl(ifac,iv,1) = - fmul * vnrm * surfbo(2,ifac) / surfbn(ifac)
rcodcl(ifac,iw,1) = - fmul * vnrm * surfbo(3,ifac) / surfbn(ifac)
if (itytur.eq.2) then
rcodcl(ifac,ik,1) = cvar_k(iel)
rcodcl(ifac,iep,1) = cvar_ep(iel)
elseif (itytur.eq.3) then
if (irijco.eq.1) then
rcodcl(ifac,ir11,1) = cvar_rij(1,iel)
rcodcl(ifac,ir22,1) = cvar_rij(2,iel)
rcodcl(ifac,ir33,1) = cvar_rij(3,iel)
rcodcl(ifac,ir12,1) = cvar_rij(4,iel)
rcodcl(ifac,ir13,1) = cvar_rij(6,iel)
rcodcl(ifac,ir23,1) = cvar_rij(5,iel)
else
rcodcl(ifac,ir11,1) = cvar_r11(iel)
rcodcl(ifac,ir22,1) = cvar_r22(iel)
rcodcl(ifac,ir33,1) = cvar_r33(iel)
rcodcl(ifac,ir12,1) = cvar_r12(iel)
rcodcl(ifac,ir13,1) = cvar_r13(iel)
rcodcl(ifac,ir23,1) = cvar_r23(iel)
endif
rcodcl(ifac,iep,1) = cvar_ep(iel)
if (iturb.eq.32) then
rcodcl(ifac,ial,1) = cvar_al(iel)
endif
elseif (itytur.eq.5) then
rcodcl(ifac,ik,1) = cvar_k(iel)
rcodcl(ifac,iep,1) = cvar_ep(iel)
rcodcl(ifac,iphi,1) = cvar_phi(iel)
if (iturb.eq.50) then
rcodcl(ifac,ifb,1) = cvar_fb(iel)
elseif (iturb.eq.51) then
rcodcl(ifac,ial,1) = cvar_al(iel)
endif
elseif (iturb.eq.60) then
rcodcl(ifac,ik,1) = cvar_k(iel)
rcodcl(ifac,iomg,1) = cvar_omg(iel)
elseif (iturb.eq.70) then
rcodcl(ifac,inusa,1) = cvar_nusa(iel)
endif
! Handle scalars (a correction similar to that of velocity is suggested
! rather than the simpler code below)
if (nscal.gt.0) then
do ii = 1, nscal
call field_get_val_s(ivarfl(isca(ii)), cvar_scal)
rcodcl(ifac,isca(ii),1) = cvar_scal(iel)
enddo
endif
enddo
endif
getfbr
subroutine getfbr(fstr, facnb, faces)
Build the list of boundary faces matching a criteria string.
Definition: cs_selector_f2c.f90:111