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 cs_user_model function in the file cs_user_parameters.cpp as follows:

if (cs_glob_physical_model_flag[CS_GAS_MIX] > -1) {
  ! Specific condensation modelling
  !      if = -1 module not activated
  !      if =  0 condensation source terms with metal
  !                               structures activate
  icondv = -1;
}

The setting of the condensation source terms is then done in the cs_user_wall_consensation function, in addition to the other settings.

The cells which are associated with the metal structures volume are defined on second call as follows:

 
  /*--------------------------------------------------------------------------
   * Select the cells which are associated to the metal structures volume.
   *--------------------------------------------------------------------------*/
 
  if (iappel == 2) {
 
    const cs_zone_t *z = cs_volume_zone_by_name("metallic_structures");
    wall_cond->ncmast = z->n_elts;
 
    for (ieltcd = 0; ieltcd < z->n_elts; ieltcd++)
      wall_cond->ltmast[ieltcd] = z->elt_ids[ieltcd];
 
  }
 

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

 
  if (iappel == 2) {
    if (wall_cond->icondb == 0) {
 
      /* izcophc: model for the mass transfer (condensation) coefficient
       * ---------------------------------------------------------------
       * 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
       * ---------------------------------------------------
       * 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
       * --------------------------------------------
       * 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
         * ---------------------------------
         * 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 cell layers in the wall mesh
         * ---------------------------------------------
         */
        wall_thermal->znmur[iz] = 10;
 
        /* zepais: Total thickness of the solid wall [meters]
         * --------------------------------------------------
         */
        wall_thermal->zepais[iz] = 0.024;
 
        /* ztpar0: Initial temperature in the solid wall [celsius]
         * -------------------------------------------------------
         */
        wall_thermal->ztpar0[iz] = 26.57;
      }
    }
  }

Finally the source term type and values have to be set (in cs_user_wall_condensation.cpp) as follows:

    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_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) {
        assert(ik != -1 && iep != -1);
        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
            }
          }
        }
      }
    }
 
    /*------------------------------------------------------------------------
     * The user can specify here the values of the following arrays used by
     * the modelling of the metal structures condensation:
     *  svcond the scalar value to multiply by the sink term array
     *  of the metal structures condensation model.
     *  This array can be filled for each transported scalar.
     *------------------------------------------------------------------------*/
 
    const cs_lnum_t ncmast = wall_cond->ncmast;
 
    cs_gas_mix_species_prop_t s_h2o_g;
    const int key_mix = cs_field_key_id("gas_mix_species_prop");
    cs_field_get_key_struct(cs_field_by_name_try("y_h2o_g"),
                            key_mix, &s_h2o_g);
 
    // 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))
        continue;
      int iscal = cs_field_get_key_int(f, keysca);
      if (iscal <= 0)
        continue;
      int ivar = cs_field_get_key_int(f, var_id_key) - 1;
 
      int h_f_id = cs_thermal_model_field()->id;
 
      /* Loop over the cells associated to the metal structure
       * source terms zone */
 
      for (cs_lnum_t icmst = 0; icmst < ncmast; icmst++) {
        const cs_lnum_t c_id = wall_cond->ltmast[icmst];
 
        /* Compute the enthalpy value of vapor gas */
        cs_real_t tk = cvar_h[c_id]/cpro_cp[c_id];
        if (cs_glob_time_step->nt_cur < 2)
          tk = cs_glob_fluid_properties->t0;
        const cs_real_t hvap = s_h2o_g.cp*tk;
 
        // wall_cond->itypst[ivar*ncmast + icmst] = 1;
        wall_cond->svcond[ivar*ncmast + icmst] = 0.0;
        if (f_id == h_f_id)
          wall_cond->svcond[ivar*ncmast + icmst] = hvap;
      }
    }
 

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:     Uchida correlation
   * - CS_WALL_COND_MODEL_DEHBI:      Dehbi correlation
   */
 
  cs_wall_condensation_t *wall_cond = cs_get_glob_wall_condensation();
  wall_cond->icondb         = 0; /* activate wall condensation */
  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) {
 
    const cs_zone_t *zone = cs_boundary_zone_by_name_try("cold_wall");
 
    if (iappel == 1) {
      wall_cond->nfbpcd = zone->n_elts;
      wall_cond->nzones = izone;
    }
    else if (iappel == 2) {
      for (ieltcd = 0; ieltcd < zone->n_elts; ieltcd++) {
        cs_lnum_t ifac             = zone->elt_ids[ieltcd];
        wall_cond->ifbpcd[ieltcd]  = ifac;
        wall_cond->izzftcd[ieltcd] = izone;
      }
    }
 
  }
 

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.

 
  if (iappel == 2) {
    if (wall_cond->icondb == 0) {
 
      /* izcophc: model for the mass transfer (condensation) coefficient
       * ---------------------------------------------------------------
       * 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
       * ---------------------------------------------------
       * 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
       * --------------------------------------------
       * 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
         * ---------------------------------
         * 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 cell layers in the wall mesh
         * ---------------------------------------------
         */
        wall_thermal->znmur[iz] = 10;
 
        /* zepais: Total thickness of the solid wall [meters]
         * --------------------------------------------------
         */
        wall_thermal->zepais[iz] = 0.024;
 
        /* ztpar0: Initial temperature in the solid wall [celsius]
         * -------------------------------------------------------
         */
        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)]
         * ------------------------------------
         */
        wall_thermal->zhext[iz] = 1.e8;
 
        /* zhext: External temperature [celsius]
         * ---------------------------
         */
        wall_thermal->ztext[iz] = 26.57;
 
        /* zrob: Solid wall density [kilogram.meter^(-3)]
         * ------------------------
         */
        wall_thermal->zrob[iz] = 8000.0;
 
        /* zcondb: Solid wall thermal conductivity [watt.meter^(-1).celsius^(-1)]
         * ---------------------------------------
         */
        wall_thermal->zcondb[iz] = 12.8;
 
        /* zcpb: Solid wall specific heat [joule.kilogram^(-1).celsius^(-1)]
         * ------------------------------
         */
        wall_thermal->zcpb[iz] = 500.0;
      }
      else {
 
        /* ztpar: Constant wall temperature [celsius]
         * --------------------------------
         */
        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_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) {
        assert(ik != -1 && iep != -1);
        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
            }
          }
        }
      }
    }
 
    /*------------------------------------------------------------------------
     * The user can specify here the values of the following arrays used by
     * the modelling of the metal structures condensation:
     *  svcond the scalar value to multiply by the sink term array
     *  of the metal structures condensation model.
     *  This array can be filled for each transported scalar.
     *------------------------------------------------------------------------*/
 
    const cs_lnum_t ncmast = wall_cond->ncmast;
 
    cs_gas_mix_species_prop_t s_h2o_g;
    const int key_mix = cs_field_key_id("gas_mix_species_prop");
    cs_field_get_key_struct(cs_field_by_name_try("y_h2o_g"),
                            key_mix, &s_h2o_g);
 
    // 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))
        continue;
      int iscal = cs_field_get_key_int(f, keysca);
      if (iscal <= 0)
        continue;
      int ivar = cs_field_get_key_int(f, var_id_key) - 1;
 
      int h_f_id = cs_thermal_model_field()->id;
 
      /* Loop over the cells associated to the metal structure
       * source terms zone */
 
      for (cs_lnum_t icmst = 0; icmst < ncmast; icmst++) {
        const cs_lnum_t c_id = wall_cond->ltmast[icmst];
 
        /* Compute the enthalpy value of vapor gas */
        cs_real_t tk = cvar_h[c_id]/cpro_cp[c_id];
        if (cs_glob_time_step->nt_cur < 2)
          tk = cs_glob_fluid_properties->t0;
        const cs_real_t hvap = s_h2o_g.cp*tk;
 
        // wall_cond->itypst[ivar*ncmast + icmst] = 1;
        wall_cond->svcond[ivar*ncmast + icmst] = 0.0;
        if (f_id == h_f_id)
          wall_cond->svcond[ivar*ncmast + icmst] = hvap;
      }
    }