programmer's documentation
For turbulence model: cs_user_turbulence_source_terms subroutine

Additional right_hand side source terms for turbulence models

Usage

The additional source term is decomposed into an explicit part $(crvexp)$ and an implic it part $(crvimp)$ that must be provided here. The resulting equations solved by the code are:

\[ \rho \norm{\vol{\celli}} \DP{\varia} + .... = crvimp \cdot \varia + crvexp \]

where $ \varia $ is the turbulence field of index $ f_{id} $.

Note
$ crvexp \text{, } crvimp$ are defined after the Finite Volume integration over the cells, so they include the $\norm{\vol{\celli}}$ term. More precisely:
  • $crvexp$ is expressed in $ kg\cdot m^{-2} \cdot s^{-2} $
  • $crvimp$ is expressed in $ kg \cdot s^{-1} $

The $crvexp$, $crvimp$ arrays are already initialized to 0 before entering the routine. It is not needed to do it in the routine (waste of CPU time).

For stability reasons, Code_Saturne will not add -crvimp directly to the diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is treated implicitely only if it strengthens the diagonal of the matrix. However, when using the second-order in time scheme, this limitation cannot be done anymore and -crvimp is added directly. The user should therefore test the negativity of crvimp by himself.

When using the second-order in time scheme, one should supply:

The selection of cells where to apply the source terms is based on a getcel command. For more info on the syntax of the getcel command, refer to the user manual or to the comments on the similar command getfbr in the routine cs_user_boundary_conditions.

Local variables

! Local variables
integer iel
double precision ff, tau
type(var_cal_opt) :: vcopt
character(len=80) :: fname
integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: cpro_rom
double precision, dimension(:), pointer :: cvar_var

Initialization and finalization

The following initialization block needs to be added for the following examples:

! Allocate a temporary array for cells selection
allocate(lstelt(ncel))

At the end of the subroutine, it is recommended to deallocate the work array:

! Deallocate the temporary array
deallocate(lstelt)

In theory Fortran 95 deallocates locally-allocated arrays automatically, but deallocating arrays in a symetric manner to their alloacation is good pratice, and avoids using a different logic C and Fortran.

Remaining initialization

Get the density array in cpro_rom

call field_get_val_s(icrom, cpro_rom)

Get the array of the current turbulent variable and its name

call field_get_val_s(f_id, cvar_var)
call field_get_name(f_id, fname)
! --- Get variable calculation options
call field_get_key_struct_var_cal_opt(f_id, vcopt)
if (vcopt%iwarni.ge.1) then
write(nfecra,1000)
endif

Example

Example of arbitrary additional source term for turbulence models (Source term on the TKE "k" here).

Source term for $\varia:$

\[ \rho \norm{\vol{\celli}} \frac{d(\varia)}{dt} = ... - \rho \norm{\vol{\celli}} \cdot ff - \rho \frac{ \varia}{\tau}\]

with $ ff $ = 3.d0 an $ \tau $ = 4.d0

Body

Note
The turbulence variable names are:
  • 'k' and 'epsilon' for the k-epsilon models
  • 'r11', 'r22', 'r33', 'r12', 'r13', 'r23' and 'epsilon' for the Rij-epsilon LRR and S SG
  • 'r11', 'r22', 'r33', 'r12', 'r13', 'r23', 'epsilon' and 'alpha' for the EBRSM
  • 'k', 'epsilon', 'phi' and 'f_bar' for the phi-model
  • 'k', 'epsilon', 'phi' and 'alpha' for the Bl-v2-k model
  • 'k' and 'omega' for the k-omega turbulence model
  • 'nu_tilda' for the Spalart Allmaras model

Calculation of the explicit and implicit source terms

if (trim(fname).eq.'k') then
ff = 3.d0
tau = 4.d0
! --- Explicit source terms
do iel = 1, ncel
crvexp(iel) = -cpro_rom(iel)*cell_f_vol(iel)*ff
enddo
! --- Implicit source terms
! crvimp is already initialized to 0, no need to set it here
do iel = 1, ncel
crvimp(iel) = -cpro_rom(iel)*cell_f_vol(iel)/tau
enddo
endif

Format

1000 format(' User source terms for turbulence model',/)