Data setting for drift scalars

Introduction

This page provides an example of code blocks that may be used to perform a calculation with drift scalars.

Physical properties

Local variables to be added

The following local variables and associated initializations need to be defined for the examples in this section:

  const cs_lnum_t n_cells = domain->mesh->n_cells;
 
  auto cpro_viscl = CS_F_(mu)->get_val_s();
 
  /* Key id for drift scalar */
  const int keydri = cs_field_key_id("drift_scalar_model");
 
  /* Key id for diffusivity id */
  const int kivisl = cs_field_key_id("diffusivity_id");
 
  /* Number of fields */
  const int nfld = cs_field_n_fields();

Body

This example sets the scalar laminar dynamic diffusivity (for Brownian motion) to take thermophoresis into account.

Here is the corresponding code:

 
  /* Loop over fields which are scalar with a drift */
 
  for (int iflid = 0; iflid < nfld; iflid++) {
 
    cs_field_t  *f = cs_field(iflid);
 
    if (! (f->type & CS_FIELD_VARIABLE))
      continue;
 
    int drift_flag = f->get_key_int(keydri);
 
    if (drift_flag & CS_DRIFT_SCALAR_ADD_DRIFT_FLUX) {
 
      /* Position of variables, coefficients
       * ----------------------------------- */
 
      cs_span<cs_real_t> cpro_taup;
      cs_span<cs_real_t> cpro_taufpt;
      cs_span<cs_real_t> cpro_viscls;
 
      /* Scalar's diffusivity (Brownian motion) */
 
      int ifcvsl = f->get_key_int(kivisl);
      if (ifcvsl > -1)
        cpro_viscls = cs_field(ifcvsl)->get_val_s();
 
      /* Coefficients of drift scalar CHOSEN BY THE USER
         Values given here are fictitious */
 
      const cs_real_t diamp = 1.e-4;   /* particle diameter */
      const cs_real_t cuning = 1.;     /* Cuningham correction factor */
      const cs_real_t rhop = 1.e4;     /* particle density */
 
      /* Get corresponding relaxation time (cpro_taup) */
      cpro_taup = cs_field_by_composite_name(f->name, "drift_tau")->get_val_s();
 
      /* Corresponding interaction time particle--eddies */
 
      if (drift_flag & CS_DRIFT_SCALAR_TURBOPHORESIS) {
        cpro_taufpt = cs_field_by_composite_name(f->name, "drift_turb_tau")->get_val_s();
      }
 
      /* Computation of the relaxation time of the particles
       * --------------------------------------------------- */
 
      const cs_real_t diamp2 = cs_math_pow2(diamp);
 
      if (diamp <= 1.e-6) {
        /* Cuningham's correction for submicronic particles */
        for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
          cpro_taup[c_id] =   cuning*diamp2*rhop / (18.*cpro_viscl[c_id]);
        }
      }
      else {
        for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
          cpro_taup[c_id] = diamp2*rhop / (18.*cpro_viscl[c_id]);
        }
      }
 
      /* Compute the interaction time particle--eddies (tau_fpt)
       * ------------------------------------------------------- */
 
      if (drift_flag & CS_DRIFT_SCALAR_TURBOPHORESIS) {
 
        const cs_turb_model_t *t_mdl = cs_glob_turb_model;
        const int ksigmas = cs_field_key_id("turbulent_schmidt");
 
        /* k-epsilon or v2-f models */
        if (t_mdl->itytur == 2 || t_mdl->itytur == 5) {
          const auto cvar_k = CS_F_(k)->get_val_s();
          const auto cvar_eps = CS_F_(eps)->get_val_s();
          const cs_real_t turb_schmidt = f->get_key_double(ksigmas);
          for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
            cs_real_t xk = cvar_k[c_id];
            cs_real_t xeps = cvar_eps[c_id];
            cpro_taufpt[c_id] = (3./2.)*(cs_turb_cmu/turb_schmidt)*xk/xeps;
          }
        }
 
        /* Rij-epsilon models */
        else if (t_mdl->itytur == 3) {
          const auto cvar_rij = CS_F_(rij)->get_val_t();
          const auto cvar_eps = CS_F_(eps)->get_val_s();
          cs_real_t beta1 = 0.5 + 3.0/(4.*cs_turb_xkappa);
          for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
            cs_real_t xk = 0.5 * (  cvar_rij(c_id, 0)
                                  + cvar_rij(c_id, 1)
                                  + cvar_rij(c_id, 2));
            cs_real_t xeps = cvar_eps[c_id];
            cpro_taufpt[c_id] = xk/xeps/beta1;
          }
        }
 
        /* k-omega models */
        if (t_mdl->itytur == 6) {
          const auto cvar_omg = CS_F_(omg)->get_val_s();
          const cs_real_t turb_schmidt = f->get_key_double(ksigmas);
          for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
            cs_real_t xomg = cvar_omg[c_id];
            cpro_taufpt[c_id] = (3./2.)*(1./turb_schmidt)/xomg;
          }
        }
 
      }
 
      /* Brownian diffusion at cell centers
       * ---------------------------------- */
 
      /* Stop if the diffusivity is not variable */
      if (ifcvsl < 0)
        bft_error(__FILE__, __LINE__, 0,
                  "The diffusivity is uniform while a variable diffusivity\n"
                  "is computed.");
 
      /* Temperature and density */
      if (CS_F_(t) == nullptr)
        bft_error(__FILE__, __LINE__, 0,
                  "The temperature field on which physical properties depend\n"
                  "does not seem to be present.");
 
      const auto cvar_t = CS_F_(t)->get_val_s();
      const auto cpro_rom = CS_F_(rho)->get_val_s();
 
      /* Homogeneous to a dynamic viscosity */
      for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
        cs_real_t xvart = cvar_t[c_id];
        cs_real_t rho = cpro_rom[c_id];
        cs_real_t viscl = cpro_viscl[c_id];
        cpro_viscls[c_id] =   rho*cs_physical_constants_kb*xvart
                            * cuning/(3.*cs_math_pi*diamp*viscl);
      }
 
    } /* End of test on drift */
 
  } /* End of loop on fields */