Infinite channel inlet

A method of defining an inlet condition which converges towards an infinite channel profile is based simply on feedback from values computed downstream from the inlet.

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

Local variables to be added

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

integer          ifac, iel, ii, ivar, iscal, ilelt, nlfac

integer          keyvar, keysca
integer          n_fields, f_id, normalize, interpolate
double precision xdh, rhomoy
double precision fmprsc, uref2

integer, allocatable, dimension(:) :: lstfac, lstcel
double precision, dimension(:), pointer :: brom

double precision, dimension(:,:), pointer :: vel
double precision, dimension(3,1) :: coord_shift
double precision, dimension(1) :: rvoid

type(c_ptr), save :: inlet_l = c_null_ptr

Note that the inlet_l pointer has the save attribute, as the matching structure usually needs to be built only once, then reused at each time step.

Initialization and finalization

Initialization and finalization is similar to that of the base examples, with the addition of the mapping object, described in a specific section hereafter

Base definition

Warning

We assume other boundary conditions are defined before this one (ideally, using the GUI).

The base definition given here ensures initialization of a (flat) inlet profile with the required mean velocity.

Note

We may skip the addition of the following block in the user subroutine if we define an equivalent inlet condition using the GUI, which will ensure the appropriate initialization before entering the user subroutine.

call getfbr('INLET', nlfac, lstfac)
!==========

do ilelt = 1, nlfac

  ifac = lstfac(ilelt)
  iel = ifabor(ifac)

  itypfb(ifac) = ientre

  rcodcl(ifac,iu,1) = fmprsc
  rcodcl(ifac,iv,1) = 0.d0
  rcodcl(ifac,iw,1) = 0.d0

  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 'cs_user_physical_properties'. 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 'cs_user_physical_properties', 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  = brom(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

Mapping definition

To define the appropriate mapping object, the boundary_conditions_map function is used. In this example, coordinates of the selected inlet face are shifted by a constant distance (5.95 meters along the x axis), and mapped to the corresponding mesh cells. Here, all cells are selected as point location candidates, but optionally, a restricted list of cells may be provided, which may accelerate location (or ensure it is done on a selected subset). In most cases, as in this example, a constant coordinate shift is used for all inlet points, but it is possible to define a specific shift for each face (defining a stride argument of 1 instead of 0 and coord_shift(:,ifac) instead of coord_shift(:,1)).

if (ntcabs.eq.ntpabs+1) then

  allocate(lstcel(ncel))

  do iel = 1, ncel
    lstcel(iel) = iel
  enddo

  coord_shift(1,1) = 5.95d0
  coord_shift(2,1) = 0.d0
  coord_shift(3,1) = 0.d0

  inlet_l = boundary_conditions_map(mesh_location_cells, ncel,                &
                                    nlfac, lstcel, lstfac,                    &
                                    coord_shift, 0,                           &
                                    0.1d0)

  deallocate(lstcel)

endif

In general, when defining a pseudo-periodic boundary condition, it is assumed that the mesh is not moving, so the mapping may be defined once and for all. Hence the test on ntcabs and ntpabs and the save attribute for the inlet_1 pointer.

Applying the map

At all time steps after the first (possibly even the first if the flow at the mapping locations is initialized to nonzero values), the values at points mapped to inlet face centers are applied to the rcodcl(ifac,1) values of the corresponding faces, using the boundary_conditions_mapped_set subroutine. Optionally, a normalization by be applied, ensuring the mean values initially defined are preserved. Normalization is recommended for the velocity, and possibly scalars, but not for turbulent quantities, which should adapt to the velocity.

if (ntcabs.gt.1) then

  call field_get_n_fields(n_fields)
  call field_get_key_id("variable_id", keyvar)
  call field_get_key_id("scalar_id", keysca)

  interpolate = 0

  do f_id = 0, n_fields-1
    call field_get_key_int(f_id, keyvar, ivar)
    if (ivar.ge.1) then
      call field_get_key_int(f_id, keysca, iscal)
      if (ivar.eq.iu .or. iscal.gt.0) then
        normalize = 1
      else
        normalize = 0
      endif
      call boundary_conditions_mapped_set(f_id, inlet_l, mesh_location_cells, &
                                          normalize, interpolate,             &
                                          nlfac, lstfac, rvoid,               &
                                          nvar, rcodcl)
    endif
  enddo

endif

Finalization

At the end of the computation, it is good practice to free the mapping structure:

if (ntcabs.eq.ntmabs) then
  call locator_destroy(inlet_l)
endif