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)
xdh = 1.d0
rhomoy = brom(ifac)
call turbulence_bc_inlet_hyd_diam(ifac, uref2, xdh, rhomoy, viscl0, &
rcodcl)
if (nscal.gt.0) then
do ii = 1, nscal
rcodcl(ifac,isca(ii),1) = 1.d0
enddo
endif
enddo
subroutine getfbr(fstr, facnb, faces)
Build the list of boundary faces matching a criteria string.
Definition: cs_selector_f2c.f90:111
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