Settings of condensation mass source terms

Introduction

Source terms modelling condensation inside the fluid domain on internal metal structures and at the boundaries can be set respectively through the subroutines cs_user_metal_structures_source_terms and cs_user_wall_condensation.

Source terms for condensation on internal metal structures

This model can be enabled in the subroutine usppmo in the file cs_user_parameters.f90 as follows:

if ( ippmod(igmix).ge.0 ) then
  ! Specific condensation modelling
  !      if = -1 module not activated
  !      if =  0 condensation source terms with metal
  !                               structures activate
  icondv = -1
endif

The setting of the condensation source terms is then done in the subroutine cs_user_metal_structures_source_terms as follows below.

The following variables need to be declared:

integer          icmst
integer          ifac, iel, iscal
integer          ivarh
integer          izone
integer          f_id

double precision hvap, tk

type(gas_mix_species_prop) s_h2o_g

double precision, dimension(:), pointer :: cpro_cp
double precision, dimension(:), pointer :: cvar_h

Necessary species physical properties can be retrieved as follows:

call field_get_id_try("y_h2o_g", f_id)
if (f_id.ne.-1) &
  call field_get_key_struct_gas_mix_species_prop(f_id, s_h2o_g)

The zones on which the condensation mass source term will be imposed can be defined as follows:

!===============================================================================
! Select the cells which are associated to the metal structures volume
! with the subroutine getcel
!===============================================================================

izone = 0
met_znb = 1

do izmet = 1 , met_znb

  if (izmet.eq.1) then
    ! Cells associated to the geometric zone
    call getcel('z > -7.0d0 and z < 53.d0', ncmast, ltmast)

    izone = izone + 1

    do icmst = 1, ncmast
      iel = ltmast(icmst)
      izmast(iel) = izone
    enddo
  endif

enddo

Modelling of the metal side (wall thermal model and metal properties can then be specified as follows:

!===============================================================================
! Parameters padding of the wall thermal model and condensation model
! ------------------------------------------------------------------
! The condensation model can be used alone with a constant temperature
! specified by the user at the cold wall (tmet=tmet0 in this case)
! or together with a 0-D thermal model. In the latter case, the two models are
! coupled.
!===============================================================================

if (icondv.eq.0) then

  ! Choice the way to impose the wall temperature (tmet)
  ! at the solid/fluid interface:
  !
  ! with the parameter itagms defined below:
  ! ----------------------------------------
  !  0 : A constant wall temperature imposed is given by the user
  !     ( tmet = tmet0 used as the wall temperature of the metal structures
  !      by the condensation model )
  !  1 : A variable wall temperature is imposed with a 0-D thermal model
  !     ( tmet = tmet(iel,1) computed by cs_metal_structures_tag.f90
  !      and used as the wall temperature of the metal structures
  !      by the condensation model )

  itagms = 1

  ! Wall temperature computed by a 0-D thermal model
  ! with an explicit scheme and variable over time.
  ! ------------------------------------------------
  ! Remark : the wall temperature unit is [°C].
  if(itagms.eq.1) then
    !----------------------------------------
    ! (xem) thickness of the metal structures
    !----------------------------------------
    xem = 0.024d0
    !-------------------------------------------
    ! Initial condition of the 0-D thermal model
    !-------------------------------------------
    tmet0 = 92.d0
    !--------------------------------------------
    ! Physical properties of the metal structures
    !--------------------------------------------
    ! (xro) density (kg.m-3)
    xro_m   = 8000.d0
    ! (xcond) thermal conductivity (W.m-1.C-1)
    xcond_m = 12.8d0
    ! (xcp)   Specific heat (J.kg-1.C-1)
    xcp_m   = 500.0d0
  else
    ! Wall temperature imposed as constant
    ! with a value specified by the user
    tmet = 92.d0
  endif

endif

Finally the source term type and values have to be set as follows:

!===============================================================================
! The user can specify here the values of the following arrays used by the
! modelling of the metal structures condensation:
!         - itypst(:,ivar) to specify the condensation source term type,
!         - svcond(:,ivar) the scalar value to multiply by the sink term array
!                          of the metal structures condensation model.
! These two arrays can be filled for each transported scalar.
!===============================================================================

!-- pointer to the specific heat
if (icp.ge.0) call field_get_val_s(icp, cpro_cp)

!-- pointer to the enthalpy value
ivarh = isca(iscalt)
call field_get_val_s(ivarfl(ivarh), cvar_h)

! loop over the cells associated to the metal structure
! source terms zone
do icmst = 1, ncmast
  iel = ltmast(icmst)

  ! Compute the enthalpy value of vapor gas
  if (ntcabs.le.1) then
    tk = t0
  else
    tk = cvar_h(iel)/cpro_cp(iel)
  endif
  hvap = s_h2o_g%cp*tk

  ! Condensation source terms associated
  ! to the metal structures imposed
  ! for each scalar.
  !----------------------------------------
  if (nscal.gt.0) then
    do iscal = 1, nscal
      if (iscal.eq.iscalt) then

        ! enthalpy value used for
        ! the explicit condensation term
        itypst(iel,isca(iscalt)) = 1
        svcond(iel,isca(iscalt)) = hvap
      else

        ! scalar values used for
        ! the explicit condensation term
        itypst(iel,isca(iscal)) = 1
        svcond(iel,isca(iscal)) = 0.d0
      endif
    enddo
  endif

enddo

Boundary source terms for condensation

The condensation of steam on surfaces can be activated by adding the following lines in function cs_user_model of file cs_user_parameters.c :

  /* Activated wall condensation model for natural convection
   *   CS_WALL_COND_MODEL_COPAIN : Legacy implementation of COPAIN correlation
   * (default) CS_WALL_COND_MODEL_COPAIN_BD : Update of COPAIN correlation from
   * Benteboula and Dabbene CS_WALL_COND_MODEL_UCHIDA : Correlation of Uchida
   *   CS_WALL_COND_MODEL_DEHBI  : Correlation of Dehbi
   */

  cs_wall_cond_t *wall_cond = cs_get_glob_wall_cond();
  wall_cond->icondb         = 0; // activate wall codnensation
  wall_cond->natural_conv_model
    = CS_WALL_COND_MODEL_DEHBI; // choose correlation

The subroutine cs_user_wall_condensation is called three times.

The first call computes the number of boundary faces and the number of zones on which a boundary mass source term is imposed, based on the selection criteria prescribed by the user.

In this example, all faces with tag "60" in the mesh are gathered in a single condensation zone.

  if (iappel == 1 || iappel == 2) {

    cs_zone_t *zone = cs_boundary_zone_by_name_try("cold_wall");

    for (ieltcd = 0; ieltcd < zone->n_elts; ieltcd++) {
      if (iappel == 2) {
        cs_lnum_t ifac             = zone->elt_ids[ieltcd];
        wall_cond->ifbpcd[ieltcd]  = ifac;
        wall_cond->izzftcd[ieltcd] = izone;
      }
    }
  }

  if (iappel == 1) {
    wall_cond->nfbpcd = ieltcd + 1;
    wall_cond->nzones = izone + 1;
  }

At the second call, connectivity tables are built between the global mesh numbering and the one dedicated to wall condensation (see snippet above). In addition, parameters related to the condensation model are set.

  /*
  ===============================================================================
   Parameters of the 1-D thermal model and condensation model
   ------------------------------------------------------------------
   Both models can be activated and coupled together or the condensation model
   can be used with a constant wall temperature specified by the user
   (at iappel=3 tpar=ztpar0(iz) in this case).
  ===============================================================================
  */

  if (iappel == 2) {
    if (wall_cond->icondb == 0) {

      /*
       *    izcophc = model for the mass transfer (condensation) coefficient
       *    ----------------------------------------------------------------
       *    Integer.
       *    1 : Turbulent wall law
       *    2 : Natural convection correlation
       *    3 : Maximum of the two previous (for mixed regime)
       *    */
      wall_cond->izcophc[iz] = 3;

      /*
       *    izcophg = model for the thermal exchange coefficient
       *    ----------------------------------------------------------------
       *    Integer.
       *    1 : Turbulent wall law
       *    2 : Natural convection correlation
       *    3 : Maximum of the two previous (for mixed regime)
       *    */
      wall_cond->izcophg[iz] = 3;

      /*
       *    iztag1d = on/off switch for 1D thermal module
       *    ----------------------------------------------------------------
       *    Integer.
       *    0 : Constant wall temperature (equal to ztpar0(iz))
       *    1 : Variable wall temperature computed with a 1D model
       *    */
      wall_cond->iztag1d[iz] = 1;

      if (wall_cond->iztag1d[iz] == 1) {

        /*
         *      ztheta = proportion of implicitation in the space discretization
         * scheme
         *      ----------------------------------------------------------------
         *      Float in the range [0, 1].
         *      Special values:
         *        0 : explicit scheme
         *        1 : implicit scheme
         *      */
        wall_thermal->ztheta[iz] = 1.0;

        /*
         *      zdxmin = Wall cell size parameters
         *      ----------------------------------------------------------------
         *      Float
         *      Special values:
         *        <=0 : Constant cell size
         *        > 0 : Variable cell size. In this case, the first cell size
         * (fluid side) is set to zdxmin [meters].
         *      */
        wall_thermal->zdxmin[iz] = 0.0;

        /*
         *      znmur = Number of cells in the wall mesh
         *      ----------------------------------------------------------------
         *      Positive integer
         *      */
        wall_thermal->znmur[iz] = 10;

        /*
         *      zepais = Total thickness of the solid wall [meters]
         *      ----------------------------------------------------------------
         *      Positive float
         *      */
        wall_thermal->zepais[iz] = 0.024;

        /*
         *      ztpar0 = Initial temperature in the solid wall [celsius]
         *      ----------------------------------------------------------------
         *      Float.
         *      */
        wall_thermal->ztpar0[iz] = 26.57;
      }
    }
  }

At the third call, properties related to the solid wall are set.

    // Fill in parameters for wall thermal mode
    // TODO : move to iappel == 2 ?
    if (wall_cond->icondb == 0) {
      if (wall_cond->iztag1d[iz] == 1) {

        /*
         *      zhext = External exchange coefficient
         * [watt.meter^(-2).kelvin^(-1)]
         *      ----------------------------------------------------------------
         *      Positive float.
         *      */
        wall_thermal->zhext[iz] = 1.e8;

        /*
         *      zhext = External temperature [celsius]
         *      ----------------------------------------------------------------
         *      Float.
         *      */
        wall_thermal->ztext[iz] = 26.57;

        /*
         *      zrob = Solid wall density [kilogram.meter^(-3)]
         *      ----------------------------------------------------------------
         *      Positive float.
         *      */
        wall_thermal->zrob[iz] = 8000.0;

        /*
         *      zcondb = Solid wall thermal conductivity
         * [watt.meter^(-1).celsius^(-1)]
         *      ----------------------------------------------------------------
         *      Positive float.
         *      */
        wall_thermal->zcondb[iz] = 12.8;

        /*
         *      zcpb = Solid wall specific heat
         * [joule.kilogram^(-1).celsius^(-1)]
         *      ----------------------------------------------------------------
         *      Positive float.
         *      */
        wall_thermal->zcpb[iz] = 500.0;
      }
      else {

        /*
         *      ztpar = Constant wall temperature [celsius]
         *      ----------------------------------------------------------------
         *      Float.
         *      */
        wall_cond->ztpar[iz] = 26.57;
      }
    }

Finally, the source terms associated with the condensation phenomenon are defined.

    const int nfbpcd = wall_cond->nfbpcd;
    for (ieltcd = 0; ieltcd < nfbpcd; ieltcd++) {

      // Enthalpy of steam
      cs_lnum_t ifac   = wall_cond->ifbpcd[ieltcd];
      cs_lnum_t iel    = ifabor[ifac];
      cs_real_t tk     = cs_glob_fluid_properties->t0;
      cs_lnum_t ntcabs = cs_get_glob_time_step()->nt_cur;
      if (ntcabs >= 1)
        tk = cvar_h[iel] / cpro_cp[iel];
      cs_real_t hvap = (cp_vap * tk);

      // Source term for momentum
      wall_cond->itypcd[iu * nfbpcd + ieltcd] = 0;
      wall_cond->itypcd[iv * nfbpcd + ieltcd] = 0;
      wall_cond->itypcd[iw * nfbpcd + ieltcd] = 0;

      wall_cond->spcond[iu * nfbpcd + ieltcd] = 0.0;
      wall_cond->spcond[iv * nfbpcd + ieltcd] = 0.0;
      wall_cond->spcond[iw * nfbpcd + ieltcd] = 0.0;

      // Source term for turbulence
      // TODO generalize to all turbulence models
      if (turb_mdl->itytur == 2) {
        wall_cond->itypcd[ik * nfbpcd + ieltcd]  = 0;
        wall_cond->itypcd[iep * nfbpcd + ieltcd] = 0;
        wall_cond->spcond[ik * nfbpcd + ieltcd]  = 0.0;
        wall_cond->spcond[iep * nfbpcd + ieltcd] = 0.0;
      }

      // Source term for scalars
      for (int f_id = 0; f_id < n_fields; f_id++) {
        f = cs_field_by_id(f_id);
        if (f->type & CS_FIELD_VARIABLE) {
          int iscal = cs_field_get_key_int(f, keysca);
          if (iscal > 0) {
            int ivar = cs_field_get_key_int(f, var_id_key) - 1;
            wall_cond->itypcd[ivar * nfbpcd + ieltcd] = 1;
            if (f == cs_thermal_model_field()) {
              wall_cond->spcond[ivar * nfbpcd + ieltcd] = hvap; // enthalpy
            }
            else {
              wall_cond->spcond[ivar * nfbpcd + ieltcd]
                = 0.0; // non-condensable species
            }
          }
        }
      }