The subroutine Examples of data settings for mass source terms (cs_user_mass_source_terms.f90) is called at three different stages in the code (iappel = 1, 2 or 3):

iappel = 1: Calculation of the number of cells where a mass source is imposed: ncesmp. Called once at the beginning of the calculation.
iappel = 2: Identification of the cells where a mass source term is imposed: array icesmp(ncesmp). Called once at the beginning of the calculation.
iappel = 3: Calculation of the values of the mass source terms. Called at each time step.

The equation for mass conservation becomes: $\der{\rho}{t}+\divs{(\rho \vect{u})}=\Gamma$

The equation for a variable $\varia$ becomes: $\frac{\Delta(\rho\varia)}{\Delta t} = ... + \Gamma(\varia^{in} - \varia)$ discretized as $\rho \dfrac{\varia^{(n+1)} - \varia^{(n)}}{\Delta t} = ... + \Gamma(\varia^{in} - \varia^{(n+1)})$

$\varia^{in}$ is the value of $\varia$ associated to the injected mass.

Two options are available:

the mass flux is injected with the local value of variable $\varia$ : $\varia^{in} = \varia^{(n+1)}$ (the equation for $\varia$ is therefore not modified)
the mass flux is injected with a specific value for $\varia$ : $\varia^{in}$ is specified by the user

Variables to be specified by the user

ncesmp: number of cells where a mass source term is imposed
icetsm(ieltsm): identification of the cells where a mass source term is imposed. For each cell where a mass source term is imposed (ielstm in [1;ncesmp]), icetsm(ieltsm) is the global index number of the corresponding cell (icestm(ieltsm) in [1;ncel])
smacel(ieltsm,ipr): value of the injection mass rate gamma (in $kg \cdot m^3 \cdot s^{-1}$ ) in the ieltsm cell with mass source term
itypsm(ieltsm,ivar): type of treatment for variable ivar in the ieltsm cell with mass source term. itypsm = 0 --> injection of ivar at local value itypsm = 1 --> injection of ivar at user specified value
smacel(ieltsm,ivar): specified value for variable ivar associated to the injected mass in the ieltsm cell with a mass source term except for ivar=ipr

Remarks

if itypsm(ieltsm,ivar)=0, smacel(ieltsm,ivar) is not used
if smacel(ieltsm,ipr)<0, mass is removed from the system, therefore Code_Saturne automatically considers $\varia^{in}=\varia^{(n+1)}$ , whatever the values of itypsm or smacel specified by the user
if a value ivar is not linked for a mass source term is imposed, no source term will be taken into account.
if a scalar doesn't evolve following the standard equation $\dfrac{\partial{(\rho \varia)}}{\partial{dt}} + \divs{(\rho \vect{u} \varia)} = ...$ (alternate convective field for instance), the source term set by this routine will not be correct (except in case of injection at the local value of the variable). The proper source term should be added directly in ustssc.

Idenfication of the cells:
The selection of cells where to apply the source term 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          ieltsm
integer          ifac, ii
integer          ilelt, nlelt
integer          izone
 
double precision wind, wind2
double precision dh, ustar2
double precision xkent, xeent
double precision flucel
double precision vtot  , gamma
 
integer, allocatable, dimension(:) :: lstelt
 
type(var_cal_opt) :: vcopt
 

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)
 
return
end subroutine cs_user_mass_source_terms

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.

First or second call

First call: iappel = 1: ncesmp: calculation of the number of cells with mass source term
Second call (if ncesmp>0): iappel = 2: icetsm: index number of cells with mass source terms

Warning

Do not use smacel in this section (it is set on the third call, iappel = 3 )
Do not use icetsm in this section on the first call (iappel=1)

This section (iappel = 1 or 2) is only accessed at the beginning of a calculation. Should the localization of the mass source terms evolve in time, the user must identify at the beginning all cells that can potentially become a mass source term.

if (iappel.eq.1.or.iappel.eq.2) then

Example 1

No mass source term (default)

ieltsm = 0

Example 2

Mass source term in the cells that have boundary face of color 3 and the cells with a coordinate X between 2.5 and 5.

In this test in two parts, one must pay attention not to count the cells twice (a cell with a boundary face of color 3 can also have a coordinate X between 2.5 and 5). One should also pay attention that, on the first call, the array icetsm doesn't exist yet. It mustn't be used outside of tests (iappel.eq.2).

  izone = 0
  ieltsm = 0
 
  ! Cells with coordinate X between 2.5 and 5.
 
  call getcel('X > 2.5 and X < 5.0',nlelt,lstelt)
 
  izone = izone + 1
 
  do ilelt = 1, nlelt
    ii = lstelt(ilelt)
    izctsm(ii) = izone
    ieltsm = ieltsm + 1
    if (iappel.eq.2) icetsm(ieltsm) = ii
  enddo
 
 
  ! Cells with a boundary face of color 3
 
  call getfbr('3',nlelt,lstelt)
 
  izone = izone + 1
 
  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    ii   = ifabor(ifac)
    ! The cells that have already been counted above are not
    ! counted again.
    if (.not.(xyzcen(1,ii).lt.500.d0.and.                     &
         xyzcen(1,ii).gt.250.d0)    )then
      ieltsm = ieltsm + 1
      izctsm(ii) = izone
      if (iappel.eq.2) icetsm(ieltsm) = ii
    endif
  enddo

Generic subsection: do not modify

For iappel = 1: specification of ncesmp. This block is valid for both examples.

  if (iappel.eq.1) then
    ncesmp = ieltsm
  endif

Third call (for ncesmp > 0)

iappel = 3: - itypsm: type of mass source term

smacel : mass source term

Remarks: If itypsm(ieltsm,var) is set to 1, smaccel(ieltsm,var) must be set.

Example 1

Simulation of an inlet condition by mass source terms and printing of the total mass rate.

  wind = 0.1d0
  wind2 = wind**2
  dh     = 0.5d0

Calculation of the inlet conditions for k and epsilon with standard laws in a circular pipe

  ustar2 = 0.d0
  xkent  = epzero
  xeent  = epzero
 
  call turbulence_bc_ke_hyd_diam(wind2, dh, ro0, viscl0,  &
                                 ustar2, xkent, xeent )
 
  flucel = 0.d0
  do ieltsm = 1, ncesmp
    smacel(ieltsm,ipr) = 30000.d0
    itypsm(ieltsm,iv) = 1
    smacel(ieltsm,iv) = wind
    if (itytur.eq.2) then
      itypsm(ieltsm,ik) = 1
      smacel(ieltsm,ik) = xkent
      itypsm(ieltsm,iep) = 1
      smacel(ieltsm,iep) = xeent
    else if (itytur.eq.3) then
      itypsm(ieltsm,ir11) = 1
      itypsm(ieltsm,ir12) = 1
      itypsm(ieltsm,ir13) = 1
      itypsm(ieltsm,ir22) = 1
      itypsm(ieltsm,ir23) = 1
      itypsm(ieltsm,ir33) = 1
      smacel(ieltsm,ir11) = 2.d0/3.d0*xkent
      smacel(ieltsm,ir12) = 0.d0
      smacel(ieltsm,ir13) = 0.d0
      smacel(ieltsm,ir22) = 2.d0/3.d0*xkent
      smacel(ieltsm,ir23) = 0.d0
      smacel(ieltsm,ir33) = 2.d0/3.d0*xkent
      itypsm(ieltsm,iep) = 1
      smacel(ieltsm,iep) = xeent
    else if (iturb.eq.50) then
      itypsm(ieltsm,ik) = 1
      smacel(ieltsm,ik) = xkent
      itypsm(ieltsm,iep) = 1
      smacel(ieltsm,iep) = xeent
      itypsm(ieltsm,iphi) = 1
      smacel(ieltsm,iphi) = 2.d0/3.d0
      ! There is no mass source term in the equation for f_bar
    else if (iturb.eq.60) then
      itypsm(ieltsm,ik) = 1
      smacel(ieltsm,ik) = xkent
      itypsm(ieltsm,iomg)= 1
      smacel(ieltsm,iomg)= xeent/cmu/xkent
    endif
    if (nscal.gt.0) then
      do ii = 1, nscal
        itypsm(ieltsm,isca(ii)) = 1
        smacel(ieltsm,isca(ii)) = 1.d0
      enddo
    endif
    flucel = flucel+                                            &
         volume(icetsm(ieltsm))*smacel(ieltsm,ipr)
  enddo
 
  if (irangp.ge.0) then
    call parsom (flucel)
  endif
 
  call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
 
  if (vcopt%iwarni.ge.1) then
    write(nfecra,1000) flucel
  endif

Example 2

Simulation of a suction (by a pump for instance) with a total rate of $80 000 \: kg \cdot s^{-1}$ . The suction rate is supposed to be uniformly distributed on all the cells selected above.

Calculation of the total volume of the area where the mass source term is imposed (the case of parallel computing is taken into account with the call to parsom).

  vtot = 0.d0
  do ieltsm = 1, ncesmp
    vtot = vtot + volume(icetsm(ieltsm))
  enddo
  if (irangp.ge.0) then
    call parsom (vtot)
  endif

The mass suction rate is $\Gamma = - \dfrac{80000}{v_{tot}}$ (in $kg \cdot m^{-3} \cdot s^{-1}$ ). It is set below, with a test for cases where $v_{tot} = 0$ . The total mass rate is calculated for verification.

  if (vtot.gt.0.d0) then
    gamma = -80000.d0/vtot
  else
    write(nfecra,9000) vtot
    call csexit (1)
  endif
 
  flucel = 0.d0
  do ieltsm = 1, ncesmp
    smacel(ieltsm,ipr) = gamma
    flucel = flucel+                                          &
         volume(icetsm(ieltsm))*smacel(ieltsm,ipr)
  enddo
 
  if (irangp.ge.0) then
    call parsom (flucel)
  endif
 
  call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt)
 
  if (vcopt%iwarni.ge.1) then
    write(nfecra,2000) flucel, vtot
  endif

End third call

endif

Formats

 1000 format(/,'Mass rate generated in the domain: ',e14.5,/)
 
 2000 format(/,'Mass flux rate generated in the domain: ',e14.5,/,         &
               '                         distributed on the volume: ',e14.5)
 
 9000 format(/,'Error in cs_user_mass_source_terms',/,                     &
               '   the volume of the mass suction area is = ',e14.5,/)