My Project
programmer's documentation
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
getfbr
subroutine getfbr(fstr, facnb, faces)
Build the list of boundary faces matching a criteria string.
Definition: cs_selector_f2c.f90:111