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