Gas 3 PTCHEM example

The following sections illustrate subsets of a same example, combining high-level definitions (which are preferred), and lower-level setings where necessary.

This is assumed to be combined with the GUI or settings not specific to gas combustion, so wall and symmetry conditions are not shown here.

Definitions for the setup stage

Using a fixed mass-flow rate, the main inlet conditions can be defined in cs_user_boundary_conditions_setup, as follows for a zone named "inlet":

  const cs_zone_t *z = cs_boundary_zone_by_name("inlet");
 
  cs_real_t q_burner = cs_notebook_parameter_value_by_name("q_burner");
 
  cs_real_t q_inlet = q_burner / 360.;                                 /* W */
  cs_real_t m_inlet = q_inlet / cs_glob_combustion_gas_model->pcigas;  /* kg/s */
  cs_real_t r_inlet = 0.15;                                            /* m */
 
  cs_boundary_conditions_open_set_mass_flow_rate_by_value(z, m_inlet);
 
  cs_boundary_conditions_inlet_set_turbulence_hyd_diam(z, 2.*r_inlet);

Note that we have used a notebook variable here, allowing parametrization of the computation, but the mass flow should be fixed. For a variable mass flow rate, use cs_boundary_conditions_open_set_mass_flow_rate_by_func instead.

Definitions for the compute stage

Not all definitions can yet be handled at the setup stage, so additional settings are needed in cs_user_boundary_conditions.

Local variables to be added

The following aliases (for conciseness of the downstream code) and initializations are used here:

 
  const cs_lnum_t *b_face_cells = domain->mesh->b_face_cells;
  const cs_real_3_t *cdgfbo
    = (const cs_real_3_t *)(domain->mesh_quantities->b_face_cog);
  const cs_real_t *gxyz = cs_glob_physical_constants->gravity;
 
  const cs_fluid_properties_t *fp = cs_glob_fluid_properties;
  cs_combustion_gas_model_t  *cm = cs_glob_combustion_gas_model;
 
  const int kbmasf = cs_field_key_id("boundary_mass_flux_id");
  const int iflmab =  cs_field_get_key_int(CS_F_(vel), kbmasf);
  const cs_real_t *bmasfl = cs_field_by_id(iflmab)->val;
 
  /* Oxydant temperature needed by the combustion model
     when there is no oxydant inlet */
 
  cs_real_t coefg[CS_COMBUSTION_GAS_MAX_GLOBAL_SPECIES];
  for (int icg = 0; icg < CS_COMBUSTION_GAS_MAX_GLOBAL_SPECIES; icg++)
    coefg[icg] = 0;
  coefg[1] = 1;
 
  cm->tinoxy = fp->t0;
  cm->hinoxy = cs_gas_combustion_t_to_h(coefg, cm->tinoxy);
 

Inlet definitions

Based on the inlet defined at the setup stage, the following additional settings can be used. Note that the inlet temperature is global to the module, but matches an inlet, so we choose to (re)set it here. These settings will be handled in the setup stage in a future version.

  {
    const cs_zone_t *z = cs_boundary_zone_by_name("inlet");
 
    /* Fuel flow inlet type */
    cs_glob_bc_pm_info->ientfu[z->id] = 1;
 
    /* Inlet Temperature in K */
    cm->tinfue = fp->t0;
  }

Open boundary definition

Since it is not yet possible to handle open boundary conditions with backflow conditions with the high level API (i.e. at the setup stage), the following example illustrates have this can be managed, on an example zone named "open":

  {
    const cs_real_3_t *cvar_vel = (const cs_real_3_t *)CS_F_(vel)->val;
 
    cs_real_t *brom = CS_F_(rho_b)->val;
 
    int *p_icodcl = CS_F_(p)->bc_coeffs->icodcl;
    cs_real_t *p_rcodcl1 = CS_F_(p)->bc_coeffs->rcodcl1;
    cs_real_t *fm_rcodcl1 = CS_F_(fm)->bc_coeffs->rcodcl1;
    cs_real_t *fp2m_rcodcl1 = CS_F_(fp2m)->bc_coeffs->rcodcl1;
    cs_real_t *h_rcodcl1 = nullptr;
    if (cs_glob_physical_model_flag[CS_COMBUSTION_3PT] == 1) {
      h_rcodcl1 = CS_F_(h)->bc_coeffs->rcodcl1;
    }
 
    const cs_real_t xdh = 3.0;     /* acos(-1.d0) * 1 ( domain circumference) */
    const cs_real_t xitur = 0.02;  /* turbulence intensity */
 
    /* Loop on zone faces */
 
    const cs_zone_t *z = cs_boundary_zone_by_name("open");
 
    for (cs_lnum_t idx = 0; idx < z->n_elts; idx++) {
      const cs_lnum_t face_id = z->elt_ids[idx];
 
      bc_type[face_id] = CS_FREE_INLET;
 
      /* Imposed pressure */
 
      const cs_real_t pimp
        = fp->p0 + fp->ro0*cs_math_3_distance_dot_product(fp->xyzp0,
                                                          cdgfbo[face_id],
                                                          gxyz);
      p_icodcl[face_id] = 1;
      p_rcodcl1[face_id] = pimp;
 
      /* Outlet: Neumann conditions by defauls */
 
      if (bmasfl[face_id] > 0)
        continue;
 
      /* Inlet (backflow): Dirichlet conditions */
 
      const cs_lnum_t c_id = b_face_cells[face_id];
 
      cs_real_t uref2 = cs_math_3_square_norm(cvar_vel[c_id]);
      uref2 = cs_math_fmax(uref2, cs_math_epzero);
 
      cs_turbulence_bc_inlet_turb_intensity(face_id, uref2, xitur, xdh);
 
      /* Mixture fraction */
      fm_rcodcl1[face_id] = 0;
 
      /* Variance */
      fp2m_rcodcl1[face_id] = 0;
 
      /* Enthalpy */
      if (h_rcodcl1 != nullptr)
        h_rcodcl1[face_id] = cm->hinoxy;
 
      /* Soot */
      if (cm->isoot == 1) {
        CS_F_(fsm)->bc_coeffs->rcodcl1[face_id] = 0;
        CS_F_(npm)->bc_coeffs->rcodcl1[face_id] = 0;
      }
 
      /* Density */
      brom[face_id] = fp->p0/(  cs_physical_constants_r
                              * cm->tinoxy/cm->wmolg[1]);
    }
  }