Infinite channel inlet

A method of defining an inlet condition which converges towards an infinite channel profile is based simply on feedback from values computed downstream from the inlet.

Here, we define an inlet boundary condition for a very long channel or duct with a section matching the boundary faces of zone 'INLET'.

Initialization and finalization

Initialization and finalization is similar to that of the base examples, with the addition of the mapping object, described in a specific section hereafter

Base definition

Warning

We assume other boundary conditions are defined before this one (ideally, using the GUI).

The base definition given here ensures initialization of a (flat) inlet profile with the required mean velocity.

Note

We may skip the addition of the following block in the user subroutine if we define an equivalent inlet condition using the GUI, which will ensure the appropriate initialization before entering the user function.

  zn = cs_boundary_zone_by_name("inlet");

  const cs_real_t fmprsc = 1.; /* mean prescribed velocity */
  const cs_real_t xdh = 1.0;   /* Hydraulic diameter */

  for (cs_lnum_t e_idx = 0; e_idx < zn->n_elts; e_idx++) {

    const cs_lnum_t face_id = zn->elt_ids[e_idx];
    const cs_lnum_t c_id = b_face_cells[face_id];

    bc_type[face_id] = CS_INLET;

    vel_rcodcl1[n_b_faces*0 + face_id] = -fmprsc;

    cs_real_t uref2 = 0;
    for (cs_lnum_t ii = 0; ii< CS_F_(vel)->dim; ii++)
      uref2 += cs_math_pow2(vel_rcodcl1[n_b_faces*ii + face_id]);
    uref2 = cs_math_fmax(uref2, 1e-12);

    /* Turbulence example computed using equations valid for a pipe.

     * We will be careful to specify a hydraulic diameter adapted
     *   to the current inlet.

     * We will also be careful if necessary to use a more precise
     *   formula for the dynamic viscosity use in the calculation of
     *   the Reynolds number (especially if it is variable, it may be
     *   useful to take the law from 'cs_user_physical_properties'
     *   Here, we use by default the 'viscl0" value.
     *   Regarding the density, we have access to its value at boundary
     *   faces (b_rho) so this value is the one used here (specifically,
     *   it is consistent with the processing in 'cs_user_physical_properties',
     *   in case of variable density) */

    /* Calculation of turbulent inlet conditions using
       the turbulence intensity and standard laws for a circular pipe
       (their initialization is not needed here but is good practice) */

    cs_real_t b_rho = bpro_rho[face_id];

    cs_turbulence_bc_inlet_hyd_diam(face_id,
                                    uref2,
                                    xdh,
                                    b_rho,
                                    viscl0);

    for (int f_id = 0; f_id < n_fields; f_id++) {
      cs_field_t *f = cs_field_by_id(f_id);

      /* Here we only handle user scalar */
     if (cs_field_get_key_int(f, keysca) > 1)
       f->bc_coeffs->rcodcl1[face_id] = 1;
    }
  }

Mapping definition

To define the appropriate mapping object, the cs_boundary_conditions_map function is used. In this example, coordinates of the selected inlet face are shifted by a constant distance (5.95 meters along the x axis), and mapped to the corresponding mesh cells. Here, all cells are selected as point location candidates, but optionally, a restricted list of cells may be provided, which may accelerate location (or ensure it is done on a selected subset). In most cases, as in this example, a constant coordinate shift is used for all inlet points, but it is possible to define a specific shift for each face.

  ple_locator_t *inlet_l = NULL;
  zn = cs_boundary_zone_by_name("inlet");

  if (nt_cur == nt_prev) {

    cs_real_t coord_shift[1][3] = {{5.95, 0, 0}};

    cs_lnum_t *cells_ids;
    BFT_MALLOC(cells_ids, n_cells_ext, cs_lnum_t);
    for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
      cells_ids[c_id] = c_id;

    inlet_l = cs_boundary_conditions_map(CS_MESH_LOCATION_CELLS,
                                         n_cells,
                                         zn->n_elts,
                                         cells_ids,
                                         zn->elt_ids,
                                         coord_shift,
                                         0,
                                         0.10);
    BFT_FREE(cells_ids);
  }

In general, when defining a pseudo-periodic boundary condition, it is assumed that the mesh is not moving, so the mapping may be defined once and for all. Hence the test on nt_cur and nt_prev and the save attribute for the inlet_1 pointer.

Applying the map

At all time steps after the first (possibly even the first if the flow at the mapping locations is initialized to nonzero values), the values at points mapped to inlet face centers are applied to the field->bc_coeffs->rcodcl1 values of the corresponding faces, using the cs_boundary_conditions_mapped_set function. Optionally, a normalization by be applied, ensuring the mean values initially defined are preserved. Normalization is recommended for the velocity, and possibly scalars, but not for turbulent quantities, which should adapt to the velocity.

  if (nt_cur == 1) {

    const int interpolate = 0;

    for (int f_id = 0; f_id < n_fields; f_id++) {

      cs_field_t *f = cs_field_by_id(f_id);

      if (!(f->type & CS_FIELD_VARIABLE))
        continue;

      cs_real_t normalize = 0;
      const int sc_id = cs_field_get_key_int(f, keysca) - 1;
      if ((f == CS_F_(vel)) || (sc_id > -1))
        normalize = 1;

      cs_boundary_conditions_mapped_set(f,
                                        inlet_l,
                                        CS_MESH_LOCATION_CELLS,
                                        normalize,
                                        interpolate,
                                        zn->n_elts,
                                        zn->elt_ids,
                                        NULL);
    }
  }

Finalization

At the end of the computation, it is good practice to free the mapping structure:

  if (nt_cur == nt_max)
    inlet_l = ple_locator_destroy(inlet_l);