User boundary condition definition for the Lagrangian model

Boundary conditions

Lagrangian boundary conditions are based on boundary zone (cs_boundary_zone_t) definitions. Additional information may be provided for Lagrangian boundary types and injections.

As usual, definitions may be created using the GUI and extended with user functions.

Access to the Lagrangian boundary conditions structure, which is necessary to most of the following examples, may be done as follows:

  cs_lagr_zone_data_t *lagr_bcs = cs_lagr_get_boundary_conditions();

Boundary zones

In this example, we assign rebound conditions to all boundary zones, except for an inlet and outlet type to specified zones. The type assigned is an integer based on the cs_lagr_bc_type_t enumerator type.

  {
     const cs_zone_t  *z;
     int n_zones = cs_boundary_zone_n_zones();
 
     /* default: rebound for all types */
     for (int z_id = 0; z_id < n_zones; z_id++) {
       lagr_bcs->zone_type[z_id] = CS_LAGR_REBOUND;
     }
 
     /* inlet and outlet for specified zones */
     z = cs_boundary_zone_by_name("inlet");
     lagr_bcs->zone_type[z->id] = CS_LAGR_INLET;
 
     z = cs_boundary_zone_by_name("outlet");
     lagr_bcs->zone_type[z->id] = CS_LAGR_OUTLET;
  }

Injection sets

In the following example, a first injection set for an inlet zone is defined. Note that newly injected particles may also be modified using the cs_user_lagr_in function.

  {
    const cs_zone_t  *z = cs_boundary_zone_by_name("inlet");
    int set_id = 0;
    cs_lagr_injection_set_t *zis
      = cs_lagr_get_injection_set(lagr_bcs, z->id, set_id);
 
    /* Now define parameters for this class and set */
 
    zis->n_inject = 100;
    zis->injection_frequency = 1;
 
    /* Assign other attributes (could be done through the GUI) */
 
    if (cs_glob_lagr_model->n_stat_classes > 0)
      zis->cluster = set_id + 1;
 
    zis->velocity_profile = 0;
    zis->velocity_magnitude = 1.1;
 
    zis->stat_weight = 1.0;
    zis->flow_rate = 0.0;
 
    /* Mean value and standard deviation of the diameter */
    zis->diameter = 5e-05;
    zis->diameter_variance = 0.0;
 
    /* Density */
 
    zis->density = 2500.0;
 
    zis->fouling_index = 100.0;
 
    /* Temperature and Cp */
 
    if (cs_glob_lagr_specific_physics->solve_temperature == 1) {
      zis->temperature_profile = 1;
      zis->temperature = 20.0;
 
      zis->cp = 1400.;
      zis->emissivity = 0.7;
    }
 
  }

In the next example, a profile is assigned to the second injection set of an inlet zone (it is assumed this et was previously defined either through the GUI or user function).

This requires first defining a profile definition function, matching the profile of cs_lagr_injection_profile_compute_t. An example based on experimental profiles is given here:

static void
_injection_profile([[maybe_unused]] int          zone_id,
                   [[maybe_unused]] int          location_id,
                   [[maybe_unused]] const void  *input,
                   cs_lnum_t                     n_elts,
                   const cs_lnum_t               elt_ids[],
                   cs_real_t                     profile[])
{
  const cs_real_3_t  *b_face_coords
    = (const cs_real_3_t *)cs_glob_mesh_quantities->b_face_cog;
 
  const int itmx = 8;
 
  /* Data initializations with experimental measurements
     --------------------------------------------------- */
 
  /* transverse coordinate */
 
  cs_real_t zi[] = {0.e-3, 1.e-3, 1.5e-3, 2.0e-3, 2.5e-3, 3.0e-3, 3.5e-3,
                    4.0e-3, 4.5e-3, 5.0e-3};
 
  /* particle volume fraction */
 
  cs_real_t lvf[] = {0.377e-4, 2.236e-4, 3.014e-4, 4.306e-4, 5.689e-4,
                     8.567e-4, 7.099e-4, 4.520e-4, 2.184e-4, 0.377e-4};
 
  /* vertical mean velocity of the particles */
 
  cs_real_t ui[] = {5.544, 8.827, 9.068, 9.169, 8.923, 8.295, 7.151, 6.048,
                    4.785, 5.544};
 
  /* Loop en elements
     ---------------- */
 
  for (cs_lnum_t ei = 0; ei < n_elts; ei++) {
 
    /* Face center */
 
    const cs_lnum_t face_id = elt_ids[ei];
 
    const cs_real_t z = b_face_coords[face_id][2];
 
    /* Interpolation */
 
    int i = 0;
 
    if (z > zi[0]) {
      for (i = 0; i < itmx; i++) {
        if (z >= zi[i] && z < zi[i+1])
          break;
      }
    }
 
    /* Compute volume fraction and statistical weight */
 
    cs_real_t up = ui[i] +(z-zi[i])*(ui[i+1]-ui[i])/(zi[i+1]-zi[i]);
    cs_real_t lvfp = lvf[i] + (z-zi[i])*(lvf[i+1]-lvf[i])/(zi[i+1]-zi[i]);
 
    /* number of particles in the cell */
 
    profile[ei] = lvfp * up;
  }
}

Assigning the profile to the injection set simply requires assigning the function to the pointer in the injection set structure:

  {
    const cs_zone_t  *z = cs_boundary_zone_by_name("inlet");
    int set_id = 1;
    cs_lagr_injection_set_t *zis
      = cs_lagr_get_injection_set(lagr_bcs, z->id, set_id);
 
    /* Assign injection profile function */
 
    zis->injection_profile_func = _injection_profile;
    zis->injection_profile_input = nullptr; /* default */
 
  }

An optional user-defined input function may also be associated.

Boundary-particle interactions

It is also possible to decide of the behavior of particle when they encounter a boundary (this boundary has to be of type CS_LAGR_BC_USER).

In the following example, the particle is simply deposited and marked for elimination:

 
# pragma omp atomic
  p_set.n_part_dep += 1;
 
# pragma omp atomic
  p_set.weight_dep += p_set.attr_real(
                                                      p_id,
                                                      CS_LAGR_STAT_WEIGHT);
 
  /* Mark particle as deposited and update its coordinates */
 
  p_set.set_flag(p_id, CS_LAGR_PART_DEPOSITED);
 
  cs_real_t *particle_coord
    = p_set.attr_real_ptr(p_id, CS_LAGR_COORDS);
  for (int k = 0; k < 3; k++)
    particle_coord[k] = c_intersect[k];
 
  /* Update event and particle state */
 
  *event_flag = *event_flag | (CS_EVENT_OUTFLOW | CS_EVENT_DEPOSITION);
  *tracking_state = CS_LAGR_PART_OUT;
 

Volume conditions

Lagrangian volume conditions are based on volume zone (cs_volume_zone_t) definitions. Additional information may be provided for Lagrangian injections.

As usual, definitions may be created using the GUI and extended with user functions.

Access to the Lagrangian volume conditions structure, which is necessary to most of the following examples, may be done as follows:

  cs_lagr_zone_data_t *lagr_vol_conds = cs_lagr_get_volume_conditions();

Injection sets

In the following example, we inject 1 particle set at each time step:

 
  {
    /* The volume zone named "particle_injection" is created in the GUI */
 
    const cs_zone_t *z = cs_volume_zone_by_name("particle_injection");
 
    /* Inject 1 particle set every time step */
 
    int set_id = 0;
 
    cs_lagr_injection_set_t *zis
      = cs_lagr_get_injection_set(lagr_vol_conds, z->id, set_id);
 
    zis->n_inject = 1000;
    zis->injection_frequency = 1; /* if <= 0, injection at
                                     initialization only */
    zis->velocity_profile = -1; /* fluid velocity */
 
    zis->stat_weight = 1.0;
 
    zis->diameter = 5e-6; /* in meters */
    zis->diameter_variance = 1e-6 * 1e-6; /* in meters squared */
 
    zis->density = 2475.;
 
    zis->fouling_index = 100.0;
 
    /* Activation of agglomeration */
    if (cs_glob_lagr_model->agglomeration == 1 ||
        cs_glob_lagr_model->fragmentation == 1 ) {
      zis->aggregat_class_id = 1;
      zis->aggregat_fractal_dim = 3.;
      zis->diameter = cs_glob_lagr_agglomeration_model->base_diameter
                      * pow((cs_real_t)zis->aggregat_class_id,
                            1./zis->aggregat_fractal_dim);
    }
 
  }

 
  {
    /* The volume zone containing all cells always has id 0;
       a given zone may otherwise be selected using cs_volume_zone_by_name() */
 
    const cs_zone_t *z = cs_volume_zone_by_id(0);
 
    /* Inject 2 particle sets of different diameters */
 
    cs_gnum_t n_inject[] = {500, 500};
 
    cs_real_t diam[] = {1e-3, 1e-2};
    cs_real_t diam_dev[] = {1e-6, 1e-5};
    cs_real_t density[] = {2500., 1500.};
 
    for (int set_id = 0; set_id < 2; set_id++) {
 
      cs_lagr_injection_set_t *zis
        = cs_lagr_get_injection_set(lagr_vol_conds, z->id, set_id);
 
      zis->n_inject = n_inject[set_id];
      zis->injection_frequency = 0; /* if <= 0, injection at
                                       initialization only */
 
      if (cs_glob_lagr_model->n_stat_classes > 0)
        zis->cluster = set_id + 1;
 
      zis->velocity_profile = -1; /* fluid velocity */
 
      zis->stat_weight = 1.0;
      zis->flow_rate   = 0;
 
      zis->diameter = diam[set_id];
      zis->diameter_variance = diam_dev[set_id] * diam_dev[set_id];
 
      zis->density = density[set_id];
      zis->fouling_index = 100.0;
 
      zis->temperature_profile = 0; /* temperature seen */
 
      zis->cp = 1400.0;
      zis->emissivity = 0.7;
 
    }
  }
 

In the following example, we inject 2 particle sets at computation initialization (i.e. at the first time step of a computation sequence in which the Lagrangian module is activated). Note that newly injected particles may also be modified using the cs_user_lagr_in function.

 
  {
    /* The volume zone containing all cells always has id 0;
       a given zone may otherwise be selected using cs_volume_zone_by_name() */
 
    const cs_zone_t *z = cs_volume_zone_by_id(0);
 
    /* Inject 2 particle sets of different diameters */
 
    cs_gnum_t n_inject[] = {500, 500};
 
    cs_real_t diam[] = {1e-3, 1e-2};
    cs_real_t diam_dev[] = {1e-6, 1e-5};
    cs_real_t density[] = {2500., 1500.};
 
    for (int set_id = 0; set_id < 2; set_id++) {
 
      cs_lagr_injection_set_t *zis
        = cs_lagr_get_injection_set(lagr_vol_conds, z->id, set_id);
 
      zis->n_inject = n_inject[set_id];
      zis->injection_frequency = 0; /* if <= 0, injection at
                                       initialization only */
 
      if (cs_glob_lagr_model->n_stat_classes > 0)
        zis->cluster = set_id + 1;
 
      zis->velocity_profile = -1; /* fluid velocity */
 
      zis->stat_weight = 1.0;
      zis->flow_rate   = 0;
 
      zis->diameter = diam[set_id];
      zis->diameter_variance = diam_dev[set_id] * diam_dev[set_id];
 
      zis->density = density[set_id];
      zis->fouling_index = 100.0;
 
      zis->temperature_profile = 0; /* temperature seen */
 
      zis->cp = 1400.0;
      zis->emissivity = 0.7;
 
    }
  }
 

In the followin example, we define an injection for rain zone for cooling towers module:

  {
    const cs_mesh_t *m = cs_glob_mesh;
 
    /* Get all cells zone */
    const cs_zone_t *z = cs_volume_zone_by_id(0);
 
    /* Inject 1 particle set every time step */
    int set_id = 0;
 
    cs_lagr_injection_set_t *zis
      = cs_lagr_get_injection_set(lagr_vol_conds, z->id, set_id);
 
    zis->n_inject = 1000;
    zis->injection_frequency = 1; /* if <= 0, injection at
                                     initialization only */
    zis->velocity_profile = -1; /* fluid velocity */
 
    zis->stat_weight = 1.0;
 
    zis->diameter = 5e-3; /* 5mm */
    zis->diameter_variance = 0.;
 
    zis->density = 1000.;
    zis->cp = 4179.0;
 
    zis->fouling_index = 0.;
 
    /* Profile of injection below the packing */
    zis->injection_profile_func = _injection_profile;
    cs_real_t flow_rate = 0;
 
    cs_field_t *f_lagr_injection_profile = cs_field_try("lagr_injection_profile");
 
    if (f_lagr_injection_profile != nullptr) {
      auto lagr_injection_profile = f_lagr_injection_profile->get_val_s();
 
      for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
        flow_rate += lagr_injection_profile[cell_id];
    }
    else
      cs_log_warning(_("%s: Trying to compute lagrangian injetion flowrate "
                       "while field \"lagr_injection_profile\" does not "
                       "exist.\n"), __func__);
 
    cs::parall::sum(flow_rate);
 
    zis->flow_rate = flow_rate;
 
  }

This requires in preambule the definition of the function pointer:

static void
_injection_profile([[maybe_unused]] int          zone_id,
                   [[maybe_unused]] int          location_id,
                   [[maybe_unused]] const void  *input,
                   cs_lnum_t                     n_elts,
                   const cs_lnum_t               elt_ids[],
                   cs_real_t                     profile[])
{
  /* Loop on elements
     ---------------- */
 
  cs_field_t *f = cs_field_try("lagr_injection_profile");
 
  if (f != nullptr) {
    auto lagr_injection_profile = f->get_val_s();
 
    for (cs_lnum_t ei = 0; ei < n_elts; ei++) {
 
      const cs_lnum_t cell_id = elt_ids[ei];
 
      /* number of particles in the cell proportional to
       * the lagr_injection_profile */
 
      profile[ei] = lagr_injection_profile[cell_id];
 
      assert(profile[ei] >= 0.);
 
    }
  }
  else
    bft_error(__FILE__, __LINE__, 0,
              "%s: Called while field \"lagr_injection_profile\" doesn't exist.\n",
              __func__);
}

External force

User definition of an external force field acting on the particles.

It must be prescribed in every cell and be homogeneous to gravity (m/s^2) By default gravity and drag force are the only forces acting on the particles (the gravity components gx gy gz are assigned in the GUI or in cs_user_parameters)

void
cs_user_lagr_ef([[maybe_unused]] cs_real_t            dt_p,
                [[maybe_unused]] const cs_lnum_t      p_id,
                [[maybe_unused]] const cs_real_t     *taup,
                [[maybe_unused]] const cs_real_3_t   *tlag,
                [[maybe_unused]] const cs_real_3_t   *piil,
                [[maybe_unused]] const cs_real_33_t  *bx,
                [[maybe_unused]] const cs_real_t      tsfext,
                [[maybe_unused]] const cs_real_3_t    gradpr,
                [[maybe_unused]] const cs_real_33_t   gradvf,
                [[maybe_unused]] cs_real_t            rho_p,
                [[maybe_unused]] cs_real_3_t          fextla)
{
  fextla[0] = 0;
  fextla[1] = 0;
  fextla[2] = 0;
}

Impose motion on a particle

Impose the motion of a particle flagged CS_LAGR_PART_IMPOSED_MOTION by modifying the particle position and its velocity.

#pragma weak cs_user_lagr_imposed_motion
void
cs_user_lagr_imposed_motion
(
  [[maybe_unused]] cs_lagr_particle_set_t       &p_set,
  [[maybe_unused]] cs_lnum_t                     p_id,
  [[maybe_unused]] const cs_real_t               coords[3],
  [[maybe_unused]] const cs_real_t               dt,
  [[maybe_unused]] cs_real_t                     disp[3]
)
{
  /* Angular velocity */
  cs_real_t omega = 1.0;
 
  /* Here we impose the particle to move arround a cylinder with
   * the axis  is (*, 0, 1) */
  cs_real_t rcost = (coords[1] - 0.0);
  cs_real_t rsint = (coords[2] - 1.0);
 
  /* Imposed displacement */
  disp[0] = 0.;
  disp[1] = rcost * (cos(omega*dt) - 1.0 ) - rsint * sin(omega*dt);
  disp[2] = rsint * (cos(omega*dt) - 1.0 ) + rcost * sin(omega*dt);
 
  /* If there is multiple class specify specific behaviour for a given class */
  if (cs_glob_lagr_model->n_stat_classes > 0) {
    if (p_set.attr_lnum(p_id, CS_LAGR_STAT_CLASS) == 1) {
      for ( int i = 0; i < 3; i++)
        disp[i] *= -2;
    }
  }
}

Modification of newly injected particles

User modification of newly injected particles.

This function is called after the initialization of the new particles in order to modify them according to new particle profiles (injection profiles, position of the injection point, statistical weights, correction of the diameter if the standard-deviation option is activated).

Warning

The modification of particles position must be made in cs_user_lagr_in_force_pos to get an initialization of particle properties coherent with the local fields.

This function is called for each injection zone and set. Particles with ids between pset->n_particles and n_elts are initialized but may be modified by this function.

void
cs_user_lagr_in
(
  [[maybe_unused]] cs_lagr_particle_set_t         &p_set,
  [[maybe_unused]] const cs_lagr_injection_set_t  *zis,
  [[maybe_unused]] const cs_lnum_t                 particle_range[2],
  [[maybe_unused]] const cs_lnum_t                 particle_face_id[],
  [[maybe_unused]] const cs_real_t                 visc_length[]
)
{
  /* Simple changes to selected attributes
     ------------------------------------- */
 
  for (cs_lnum_t p_id = particle_range[0]; p_id < particle_range[1]; p_id++) {
 
    /* velocity */
 
    auto part_vel = p_set.attr_n_real_ptr(p_id, 1, CS_LAGR_VELOCITY);
    part_vel[0] = 1.0;
    part_vel[1] = 0.0;
    part_vel[2] = 0.0;
 
    /* diameter */
 
    p_set.attr_real(p_id, CS_LAGR_DIAMETER) = 5e-05;
 
    /* Temperature profile */
 
    p_set.attr_real(p_id, CS_LAGR_TEMPERATURE) = 20.0;
 
    /* Statistical weight profile */
 
    p_set.attr_real(p_id, CS_LAGR_STAT_WEIGHT) = 0.01;
 
    /* User variables */
 
    for (int attr_id = CS_LAGR_USER;
         attr_id < CS_LAGR_USER + cs_glob_lagr_model->n_user_variables;
         attr_id++) {
      cs_real_t *user_var = p_set.attr_real_ptr(p_id, (cs_lagr_attribute_t) attr_id);
      *user_var = 0.;
    }
 
  }
 
  if (particle_range[1] <= particle_range[0])
    return;
 
  /* Modifications occur after all the initializations related to
     the particle injection. */
 
  /* if new particles have entered the domain  */
 
  for (cs_lnum_t p_id = particle_range[0]; p_id < particle_range[1]; p_id++) {
    _inlet2(p_set, p_id);
  }
 
  /*
   * Trick to average the statistics at iteration nstist
   * starting from an unsteady two-coupling calculation
   */
 
  if (cs_glob_time_step->nt_cur > cs_glob_lagr_stat_options->nstist) {
    cs_glob_lagr_source_terms->nstits = cs_glob_lagr_stat_options->nstist;
    cs_glob_lagr_time_scheme->isttio = 1;
  }
 
}

To modify the particle position use:

#pragma weak cs_user_lagr_in_force_coords
void
cs_user_lagr_in_force_coords
(
  [[maybe_unused]] cs_lagr_particle_set_t         &p_set,
  [[maybe_unused]] const cs_lagr_injection_set_t  *zis,
  [[maybe_unused]] const cs_lnum_t                 particle_range[2],
  [[maybe_unused]] const cs_lnum_t                 particle_face_id[],
  [[maybe_unused]] const cs_real_t                 visc_length[]
)
{
  cs_real_t dest[3]  = {0., 0., 0.};
  for (cs_lnum_t p_id = particle_range[0]; p_id < particle_range[1]; p_id++) {
    cs_real_t *part_coord = p_set.attr_real_ptr(p_id, CS_LAGR_COORDS);
    for (int i = 0; i < 3; i++) {
      part_coord[i] = dest[i];
    }
  }
}

Here is another example of the modification of newly injected particles:

 
  const cs_coal_model_t  *cm = cs_glob_coal_model;
 
  /* Changes to selected attributes to define coal
     --------------------------------------------- */
 
  constexpr int n_max_layers = 5;
 
  const int  n_layers = cs_glob_lagr_model->n_temperature_layers;
  assert(n_layers <= n_max_layers);
 
  int        coal_id = 0;
  cs_real_t  cp = zis->cp;
  cs_real_t  water_mass_f = 0.0, density = 0;
  cs_real_t  initial_diameter = 0, shrinking_diameter = 0;
 
  cs_real_t  temperature[n_max_layers];
  cs_real_t  coal_mass_fraction[n_max_layers];
  cs_real_t  coke_density[n_max_layers];
  cs_real_t  coke_mass_fraction[n_max_layers];
 
  for (int l_id = 0; l_id < n_layers; l_id++) {
    temperature[l_id] = 0;
    coal_mass_fraction[l_id] = 0;
    coke_density[l_id] = 0;
    coke_mass_fraction[l_id] = 0;
  }
 
  /* Determine coal properties based on injection zone and set
     --------------------------------------------------------- */
 
  if (zis->zone_id == 1 && zis->set_id == 0) {
 
    coal_id = 0;
 
    density
      =   cm->xashch[coal_id] * cm->rho0ch[coal_id]
        +    (1.0 - cm->xwatch[coal_id] - cm->xashch[coal_id])
           * cm->rho0ch[coal_id]
           * (1.0 - (cm->y1ch[coal_id] + cm->y2ch[coal_id]) / 2.0);
 
    cp = cm->cp2ch[coal_id]; /* specific heat */
    water_mass_f = 0.0;      /* water mass fraction */
 
    for (int l_id = 0; l_id < n_layers; l_id++) {
 
      /* temperature profile (in degrees C) */
      temperature[l_id] = 800 - 273.15;
 
      /* reactive coal mass fraction */
      coal_mass_fraction[l_id] = 0.;
 
      /* coke density after pyrolysis */
      coke_density[l_id]
        =   (1.0 - cm->xwatch[coal_id] - cm->xashch[coal_id])
          * cm->rho0ch[coal_id]
        * (1.0 - (cm->y1ch[coal_id] + cm->y2ch[coal_id]) / 2.0);
 
      /* coke mass fraction */
      coke_mass_fraction[l_id] = coke_density[l_id] / density;
 
    }
 
    /* coke diameter */
    shrinking_diameter = zis->diameter;
 
    /* initial particle diameter */
    initial_diameter = zis->diameter;
 
  }
 
  /* Now set newly injected particle values
     -------------------------------------- */
 
  for (cs_lnum_t p_id = particle_range[0]; p_id < particle_range[1]; p_id++) {
 
    /* specific heat */
    p_set.attr_real(p_id, CS_LAGR_CP) = cp;
 
    /* water mass fraction in the particle */
    p_set.attr_real(p_id, CS_LAGR_WATER_MASS) = water_mass_f;
 
    cs_real_t diam = p_set.attr_real( p_id, CS_LAGR_DIAMETER);
    cs_real_t mass = density * cs_math_pi/6 * (diam*diam*diam);
 
    p_set.attr_real(p_id, CS_LAGR_MASS) = mass;
    p_set.attr_real(p_id, CS_LAGR_WATER_MASS) = water_mass_f * mass;
 
    cs_real_t *particle_temperature
      = p_set.attr_real_ptr(p_id, CS_LAGR_TEMPERATURE);
    cs_real_t *particle_coal_mass
      = p_set.attr_real_ptr(p_id, CS_LAGR_COAL_MASS);
    cs_real_t *particle_coke_mass
      = p_set.attr_real_ptr(p_id, CS_LAGR_COKE_MASS);
    cs_real_t *particle_coal_density
      = p_set.attr_real_ptr(p_id, CS_LAGR_COAL_DENSITY);
 
    for (int l_id = 0; l_id < n_layers; l_id++) {
      particle_temperature[l_id] = temperature[l_id];
      particle_coal_mass[l_id] = coal_mass_fraction[l_id] * mass / n_layers;
      particle_coke_mass[l_id] = coke_mass_fraction[l_id] * mass / n_layers;
      particle_coal_density[l_id] = coke_density[l_id];
    }
 
    p_set.attr_real(p_id, CS_LAGR_SHRINKING_DIAMETER) = shrinking_diameter;
    p_set.attr_real(p_id, CS_LAGR_INITIAL_DIAMETER) = initial_diameter;
  }