Low-level user-defined functions for boundary conditions

For definitions using the legacy finite-volume scheme (i.e. not using CDO), the cs_user_boundary_conditions (in C) may be used. The Fortran equivalent, cs_f_user_boundary_conditions may still be used but is deprecated, and not described here (please refer to the documentation from previous versions if needed).

For more details about the treatment of boundary conditions, the user may refer to the theoretical and computer documentation [@theory] of the function cs_boundary_condition_set_coeffs (for wall conditions, see clptur) (to access this document on a workstation, use code_saturne info –guide theory).

From the user point of view, the boundary conditions are fully defined by the arrays specifying the boundary type and conditions assigned to each boundary face.

An array common to all variables defines the boundary condition type relative to the base flow:

bc_type[n_b_faces]

Other arrays are specific to each solved variable field, and accessible through its cs_field_t::bc_coeffs member:

icodcl[n_b_faces] defines the type of boundary condition for the variable.
rcodcl1[n_b_faces*dim], rcodcl2[n_b_faces*dim], and rcodcl3[n_b_faces*dim] contain the associated numerical values (value of the Dirichlet condition, of the flux ...).

In the case of standard boundary conditions, even without the GUI, it is sufficient to complete bc_type[face_id] and parts of rcodcl* arrays, as icodcl and most of rcodcl* are filled automatically based on the boundary condition type. For non-standard boundary conditions, those arrays must be fully completed.

icodcl	Description
1	Prescribed value (Dirichlet condition)
2	Convective boundary condition
3	Flux condition
4	Symmetry condition
5	Wall law (friction for velocity, wall value for scalars)
6	Rough wall law (deprecated))
9	Free outlet condition (for velocity)
11	Neumann on the normal component, Dirichlet on tangential components.
12	Imposed value for convection proportional to boundary cell value, imposed flux for diffusion (Robin condition)
13	Imposed value for convection, imposed flux for diffusion
14	Generalized symmetry condition for vectors (Marangoni effect for velocity for example)
15	User-defined wall law for scalars; for pressure, Neumann for advection and homogeneous Neumann for diffusion (walls with hydro. pressure for compressible module)

Note that in some cases, a negative value may be used, implying an automatic conversion:

For pressure, Dirichelt definition relative to the solved pressure (P*), (by default, definition relative to the hydrostatic pressure).
For enthalpy, conversion from given temperature

Coding of standard boundary conditions

The standard keywords used by the indicator bc_type are, in C: CS_INLET, CS_SMOOTHWALL, CS_ROUGHWALL, CS_SYMMETRY, CS_OUTLET, CS_FREE_INLET, CS_FREE_SURFACE, CS_CONVECTIVE_INLET and CS_INDEF.

If bc_type[face_id] == CS_INLET: inlet face.
- Zero-flux condition for pressure and Dirichlet condition for all other variables. The value of the Dirichlet condition must be given in f->bc_coeffs->rcodcl1[coo_id*n_b_faces + face_id] for each solved variable, except for the pressure (where coo_id is the coordinate id for vector or tensor fields, and 0 for scalars). The values of icodcl, rcodc2, and rcodcl3 and are filled automatically.
If bc_type[face_id] == CS_SMOOTHWALL: smooth solid wall face, impermeable and with friction.
- the eventual sliding wall velocity of the face is defined by CS_F_(vel)->bc_coeffs->rcodcl1[coo_id*n_b_faces + face_id] (with coo_id 0, 1, or 2). The initial rcodcl1 values of are zero for the three velocity components (and therefore need to be specified only if the velocity is not equal to zero).
  Warning
  The sliding wall velocity must belong to the boundary face plane. For safety, the code only uses the projection of this velocity on the face. As a consequence, if the velocity specified by the user does not belong to the face plane, the wall sliding velocity really taken into account will be different.
- For scalars, two kinds of boundary conditions can be defined:
  - Imposed value at the wall. The user must write:
    icodcl[face_id] = 5
    rcodcl1[face_id] = imposed value.
  - Imposed flux at the wall. The user must write:
    icodcl[face_id] = 3
    rcodcl3[face_id] = imposed flux value
    (depending on the variable, the user may refer to the case icodcl=3 of Coding of non-standard boundary conditions for the flux definition).
  - If the user does not fill these arrays, the default condition is zero flux.
If bc_type[face_id] == CS_ROUGHWALL: rough solid wall face, impermeable and with friction.
- The same definitions for scalars and for eventual sliding wall velocity apply as for a smooth wall.
- Roughness is also specified at each boundary face though associated fields: boundary_roughness for the dynamic roughness
  boundary_thermal_roughness for the thermal and scalar roughness.
If bc_type[face_id] === CS_SYMMETRY: symmetry face (or wall without friction).
- Nothing to specify in icodcl and rcodcl arrays.
If bc_type[face_id] == CS_OUTLET: free outlet face (or more precisely free inlet/outlet with forced pressure)

The pressure is always treated with a Dirichlet condition, calculated with the constraint: $\frac{\partial }{\partial n}\left(\frac{ \partial P}{\partial \tau}\right) = 0$

The pressure is set to P₀ at the CS_OUTLET face closest to the reference point defined by cs_glob_fluid_properties->xyzp0. The pressure calibration is always done on a single face, even if there are several outlets.

If the mass flow is incoming, the velocity is set to zero and a Dirichlet condition for the scalars and the turbulent quantities is used (or zero-flux condition if no Dirichlet value has been specified).

If the mass flow is ougoing, zero-flux condition are set for the velocity, turbulent quantities, and scalars.

Nothing is set in icodcl or rcodcl for the pressure or the velocity. An optional Dirichlet condition can be specified for the scalars and turbulent quantities.
If bc_type[face_id] == CS_FREE_INLET: free outlet or inlet (based on Bernoulli relationship) face.
- If outlet, the equivalent to standard outlet. In case of ingoing flux, the Bernoulli relationship which links pressure and velocity is used (see the theory guide for more information). An additional head loss modelling the outside of the domain can be added by the user.
If bc_type[face_id] == CS_FREE_SURFACE: free-surface boundary condition.
If bc_type[face_id] == CS_CONVECTIVE_INLET: inlet with zero diffusive flux for all transported variables (species and velocity).
- This allows to exactly impose the ingoing flux, without adding a diffusive term.
If bc_type[face_id] == CS_INDEF: undefined type face (non-standard case).
- Coding is done in a non-standard way by filling icodcl, rcodcl1, rcodcl2, and rcodcl3 entries associated to that face (see Coding of non-standard boundary conditions).

Remarks

Whatever the value of the indicator bc_type[face_id], if the array icodcl[face_id] is modified for a given variable (i.e. filled with a non-zero value), the code will not use the default conditions for that variable at face face_id. It will take into account only the values of icodcl and rcodcl provided by the user (these arrays must then be fully completed, like in the non-standard case).

For instance, for a normal symmetry face where scalar 1 is associated with a Dirichlet condition equal to 23.8 (with an infinite exchange coefficient) for a given variable field with pointer f:

bc_type[face_id] = CS_SYMMETRY;
f->bc_coeffs->icodcl[face_id]) = 1;
f->bc_coeffs->rcodcl1[face_id] = 23.8;
 
// rcodcl2[face_id] = cs_math_infinite_r is already the default value.

The boundary conditions for the other variables are defined automatically.

The gradient boundary conditions in code_saturne boil down to determine a value for the current variable Y at the boundary face f_b, that is to say $\varia_\fib$ , value expressed as a function of $\varia_{\centip}$ , value of Y in I', projection of the center of the adjacent cell on the straight line perpendicular to the boundary face and crossing its center:
$\varia_\fib=A_{\fib}^g +B_{\fib}^g \varia_{\centip}.$

For a given face, the pair of coefficients $A_{\fib}^g , \, B_{\fib}^g$ may be accessed using the f->bc_coeffs->a[face_id] and f->bc_coeffs->b[face_id] arrays, where the f is a pointer to the variable's field structure.
- In the case of a vector or tensor, where d represents f->dim (3 or 6 respectively), f->bc_coeffs->a[face_id} is replaced by f->bc_coeffs->a[d*face_id + i] for coordinate i in the expressions above, and f->bc_coeffs->b[face_id} is replaced by f->bc_coeffs->b[d*d*face_id + d*i + j].
The flux boundary conditions in code_saturne boil down to determine the value of the diffusive flux of the current variable Y at the boundary face f_b, that is to say the $D_{\ib} \left(K_\fib, \, \varia \right)$ , value expressed as a function of $\varia_{\centip}$ , value of Y in I', projection of the center of the adjacent cell on the straight line perpendicular to the boundary face and crossing its center:

$D_{\ib} \left(K_\fib, \, \varia \right) = A_{\fib}^f +B_{\fib}^f \varia_{\centip}.$

For a given face, the pair of coefficients $A_{\fib}^f , \, B_{\fib}^f$ may be accessed using the f->bc_coeffs->af[face_id] and f->bc_coeffs->bf[face_id] arrays, where the f is a pointer to the variable's field structure.
- In the case of a vector or tensor, where d represents f->dim (3 or 6 respectively), f->bc_coeffs->af[face_id} is replaced by f->bc_coeffs->af[d*face_id + i] for coordinate i in the expressions above, and f->bc_coeffs->bf[face_id} is replaced by f->bc_coeffs->bf[d*d*face_id + d*i + j].
The divergence boundary conditions in code_saturne boil down to determining a value for the current variable Y (mainly the Reynolds stress components, the divergence $\divv \left(\tens{R} \right)$ used in the calculation of the momentum equation) at the boundary face f_b, that is to say $\varia_\fib$ , value expressed as a function of $\varia_{\centip}$ , value of Y in I', projection of the center of the adjacent cell on the straight line perpendicular to the boundary face and crossing its center:
$\varia_\fib=A_{\fib}^d +B_{\fib}^d \varia_{\centip}.$

For a given face, the pair of coefficients $A_{\fib}^d , \, B_{\fib}^d$ may be accessed using the f->bc_coeffs->ad[face_id] and f->bc_coeffs->bd[face_id] arrays, where the f is a pointer to the variable's field structure.
- In the case of a vector or tensor, where d represents f->dim (3 or 6 respectively), f->bc_coeffs->ad[face_id} is replaced by f->bc_coeffs->ad[d*face_id + i] for coordinate i in the expressions above, and f->bc_coeffs->bd[face_id} is replaced by f->bc_coeffs->bd[d*d*face_id + d*i + j].
Caution: to prescribe a flux (nonzero) to Rij, the viscosity to take into account is the molecular_viscosity field even if the turbulent viscosity field (ρ.C<sub>μ/ε).

Coding of non-standard boundary conditions

If a face does not correspond to a standard type, the user must completely fill the arrays bc_type, icodcl, rcodcl1, rcodcl2, and rcodcl3 fr each variable field. bc_type[face_id] is then equal to CS_INDEF or another value defined by the user (see note at the end of Coding of standard boundary conditions). The icodcl and rcodcl arrays must be filled as follows:

If f->bc_coeffs->icodcl[face_id] == 1: Dirichlet condition.
- f->bc_coeffs->rcodcl1[face_id] is the value of the variable at the given face.
  - For vectors and tensors, f->bc_coeffs->rcodcl1[face_id] is replaced by f->bc_coeffs->rcodcl1[n_b_faces*i + face_id] for coordinate i.
  - This value has the units of the variable:
    - m/s for the velocity
    - m²/s² for the Reynolds stress
    - m²/s³ for the dissipation
    - Pa for the pressure
    - °C for the temperature
    - J.kg ^-1 for the enthalpy
    - °C² for temperature fluctuations
    - J².kg^-2 for enthalpy fluctuations
- f->bc_coeffs->rcodcl2[face_id] is the value of the exchange coefficient between the outside and the fluid for the variable (usually referred to as "exterior" exchange coefficient)..
  - An "infinite" value (rcodcl2[face_id] == cs_math_infinite_r) indicates an ideal transfer between the outside and the fluid (default case).
  - It has the following units (defined in such way that when multiplying the exchange coefficient by the variable, the given flux has the same units as the flux defined below when icodcl == 3):
    - kg.m ^-2.s ^-1 for the velocity
    - kg.m ^-2.s ^-1 for the Reynolds stress
    - s.m ^-1 for the pressure
    - W.m ^-2.°C ^-1 for the temperature
    - kg.m ^-2.s ^-1 for the enthalpy
- f->bc_coeffs->rcodcl3[face_id] is not used.
If f->bc_coeffs->icodcl[face_id] == 2: radiative outlet.

It reads $\dfrac{\partial \varia }{\partial t} + C \dfrac{\partial \varia}{\partial n} = 0$ , where C is a to be defined celerity of radiation.
- f->bc_coeffs->rcodcl1[face_id] is the value of the variable at I', projection of the center of the adjacent cell on the straight line perpendicular to the boundary face and crossing its center, at the previous time step.
- f->bc_coeffs->rcodcl2[face_id] is CFL number based on the parameter C, the distance to the boundary I'F and the time step: $CFL = \dfrac{C dt }{\centip \centf}$ .
- f->bc_coeffs->rcodcl3[face_id] is not used.
If f->bc_coeffs->icodcl[face_id] == 3: flux condition.
- f->bc_coeffs->rcodcl1[face_id] and f->bc_coeffs->rcodcl2[face_id] are not used.
- f->bc_coeffs->rcodcl3[face_id] is the flux value of of the variable at the wall. This flux is negative if it is a source for the fluid. It corresponds to:
  - $-(\lambda_T+C_p\frac{\mu_t}{\sigma_T})\grad T\cdot\vect{n}$ for a temperature (in )
  - $-(\frac{\lambda_T}{C_p}+\frac{\mu_t}{\sigma_h})\grad h\cdot\vect{n}$ for an enthalpy (in ).
  - $-(\lambda_\varphi+\frac{\mu_t}{\sigma_\varphi})\grad\varphi\cdot\vect{n}$ in the case of another scalar $\varphi$ (in $kg.m^{-2}.s^{-1}.[\varphi]$ , where $[\varphi]$ are the units of $\varphi$ ).
  - $-\Delta t\ \grad P\cdot\vect{n}$ for the pressure (in $kg.m^{-2}.s^{-1}$ ).
  - $-(\mu+\mu_t)\grad U_i\cdot\vect{n}$ for a velocity component (in $kg.m^{-1}.s^{-2}$ ).
  - $-\mu\grad R_{ij}\cdot\vect{n}$ for a R_ij-ε tensor component (in ).
If f->bc_coeffs->icodcl[face_id] == 4: symmetry condition, for symmetry faces or wall faces without friction.

This condition can only be( used for velocity components ( $\vect{U}\cdot\vect{n} = 0$ ) and the R_ij-ε tensor components (for other variables, a zero-flux condition type is usually used).
If f->bc_coeffs->icodcl[face_id] == 5: friction condition, for wall faces with friction. This condition can not be applied to the pressure.
- For the velocity and (if necessary) the turbulent variables, the values at the wall are calculated from theoretical profiles. In the case of a sliding wall, the three components of the sliding velocity are given by CS_F_(vel)->bc_coeffs->rcodcl1[coo_id*n_b_faces + face_id] (with coo_id 0, 1, or 2), the same as for a standard wall.
- For scalars, the condition icodcl[face_id] == 5 is similar to icodcl[face_id] == 1, but with a wall exchange coefficient calculated from a theoretical law. Therefore, the values of f->bc_coeffs->rcodcl1[face_id] and f->bc_coeffs->rcodcl2[face_id] must be specified: see [@theory].
If f->bc_coeffs->icodcl[face_id] == 5: friction condition, for rough wall faces with friction. This condition can not be applied to the pressure.
- The same rules apply as for a smooth wall, above.
- The dynamic roughness height is given by CS_F_(vel)->bc_coeffs->rcodcl3[face_id] only.
- For the other scalars, the thermal/scalar roughness height is given by f->bc_coeffs->rcodcl3[face_id] only.
If f->bc_coeffs->icodcl[face_id] == 9: free outlet condition for the velocity. This condition is only applicable to velocity components.
If the mass flow at the face is negative, this condition is equivalent to a zero-flux condition.
If the mass flow at the face is positive, the velocity at the face is set to zero (but not the mass flow).
rcodcl is not used.
If f->bc_coeffs->icodcl[face_id] == 14: generalized symmetry boundary condition for vectors (Marangoni effect for the velocity for instance).

This condition is only applicable to vectors and sets a Dirichlet boundary condition on the normal component and a Neumann condition on the tangential components.
For each of the 3 components coo_id, the required values are:

f->bc_coeffs->rcodcl1[coo_id*n_b_faces + face_id]: Dirichlet value for the coo_id coordinate.

f->bc_coeffs->rcodcl3[coo_id*n_b_faces + face_id]: flux value for the coo_id coordinate.

Therefore, the code automatically computes the boundary condition to impose to the normal and to the tangential components.

Consistency rules summary

In short, following consistency rules between icodcl codes for variables with non-standard boundary conditions:

Codes for vector or tensor components must be identical (the icodcl array is handled as a scalar array even for vectors and tensors).
If icodcl = 4 for velocity or Rij, it must be 4 for both.
If icodcl = 5 for a scalar or fluctuations, it must be 5 for the velocity.
- The same rule applies to icodcl = 6.
If icodcl = 5 for velocity or turbulence variables, it must be 5 for all such variables.
- The same rule applies to icodcl = 6.

Specific cases

Outlet faces

A standard CS_OUTLET outlet face amounts to a Dirichlet condition (icodcl == 1) for the pressure, a free outlet condition (icodcl == 9) for the velocity, and a Dirichlet condition (icodcl == 1) if the user has specified a Dirichlet value or a zero-flux condition (icodcl == 3) for the other variables.

Boundary condition types for enthalpy

For enthalpy, prescribed values in CS_F_(h)f->bc_coeffs->rcodcl3 may be defined using the temperature instead. In this case, CS_F_(h)f->bc_coeffs->icodcl must be replaced by - CS_F_(h)f->bc_coeffs->icodcl (i.e. its sign must bev inverted) to mark the face for automatic conversion.

Boundary condition types for compressible flows

For compressible flows, only one of the following boundary condition types be assigned:

CS_SMOOTHWALL (standard wall)
CS_SYMMETRY (standard symmetry)
CS_ESICF, CS_SSPCF, CS_SOPCF, CS_EPHCF, or CS_EQHCF (inlet/outlet)

Boundary conditions for ALE (Arbitrary Lagrangian Eulerian)

Fixed displacement on vertices

For a better precision concerning mesh displacement, one can also assign values of displacement to certain internal and/or boundary vertices. To do this one needs to fill the depale array and the vertex-based "mesh_displacement" field values. This property field tracks the total displacement of vertices relative to their initial position in the mesh. impale[vtx_id] = 1 indicates that the displacement of vertex vtx_id is imposed.

Note: Note that the impale array is initialized to 0; if its value is not modified, corresponding value in depale array will be ignored.

During the mesh deformation calculation at each time step, the position of the vertices whose displacement is fixed (i.e. impale=1) is not computed using the value of mesh velocity at the center of corresponding cell, but directly filled using the values of depale.

If the displacement is fixed for all vertices of a boundary face it is not necessary to prescribe boundary conditions at this face on mesh velocity. bc_coeffs->icodcl and bc_coeffs->rcodcl1 values of the "mesh_velocity" field will be overwritten:

icodcl is automatically set to 1 (Dirichlet)
rcodcl1 will be automatically set to face's mean mesh velocity value, that is calculated using the "mesh_displacement" field's values.

If a fixed boundary condition (ialtyb[face_id]) = CS_BOUNDARY_ALE_FIXED) is imposed to the face face_id, the displacement of each vertex vtx_id belonging to that face is considered to be fixed, meaning that impale[vtx_id] = 1 and displacment[vtx_id][*] = 0 (where mesh_displacement id a pointer to the "mesh_displacement's" values array).

Influence on boundary conditions related to fluid velocity

The effect of fluid velocity and ALE modeling on boundary faces that are declared as walls really depends on the physical nature of this interface.

Indeed when studying an immersed structure the motion of corresponding boundary faces is the one of the structure, meaning that it leads to fluid motion. On the other hand when studying a piston the motion of vertices belonging to lateral boundaries has no physical meaning, therefore it has no influence on fluid motion.

Whatever the case, the mesh velocity component that is normal to the boundary face is always taken into account ( $\vect{u}_{fluid} \cdot \vect{n} = \vect{w}_{mesh} \cdot \vect{n}$ ).

The modeling of tangential mesh velocity component differs from one case to another.

The influence of mesh velocity on boundary conditions for fluid modeling is managed and modeled in code_saturne as follows (using array names from cs_user_boundary_conditions_ale):

If ale_bc_type[face_id] = CS_BOUNDARY_ALE_FIXED: the face is motionless; mesh velocity equals 0.
If ale_bc_type[face_id] = CS_BOUNDARY_ALE_IMPOSED_VEL: tangential mesh velocity is modeled as a sliding wall velocity in fluid boundary conditions unless a value for fluid sliding wall velocity has been specified.
If ale_bc_type[face_id] = CS_BOUNDARY_ALE_SLIDING: tangential mesh velocity is not taken into account in fluid boundary conditions.
If `impale[vtx_id] = 1 for all vertices of a boundary face: tangential mesh velocity value that has been derived from vertices displacement is modeled as a sliding wall velocity in fluid boundary conditions unless a value for fluid sliding wall velocity has been specified by USER in code_saturne Interface or in 'cs_user_boundary_conditions' function.

Note that mesh velocity has no influence on modeling of boundary faces with "inlet" or "free outlet" fluid boundary conditions.

For "non standard" conditions, the user has to manage the influence of ALE boundary conditions (i.e. mesh velocity) on the ones for Navier Stokes equations(i.e. fluid velocity). (Note that fluid boundary conditions can be specified in this function.)

Combining approaches

Definitions may be based on standard boundary conditions and extended though non-stand conditions. For example, bc_type[face_id] can be set to CS_SMOOTHWALL (whether through the GUI or cs_user_boundary_conditions), and for a specific variable, the associated icodcl and rcodcl* arrays may be modified.

As always, it is recommended to specify only the values which need to be modified relative to the GUI definitions and default values, so as to keep user-defined functions concise, readable, and maintainable.