Electric arcs boundary conditions examples

Initialization and finalization

Initialization and finalization is similar to that of the base examples

Example 1

For boundary faces of zone "1" assign an inlet and assign a cathode for "electric" variables.

  {
    const cs_zone_t *zone = cs_boundary_zone_by_name("1");
    cs_field_t *th_f = cs_thermal_model_field();
    const auto bfpro_rom =  CS_F_(rho_b)->get_val_s();
 
    cs_real_t *rcodcl1_vel = CS_F_(vel)->bc_coeffs->rcodcl1;
 
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_INLET;
 
      /* Velocity */
      rcodcl1_vel[n_b_faces*0 + face_id] = 0;
      rcodcl1_vel[n_b_faces*1 + face_id] = 0;
      rcodcl1_vel[n_b_faces*2 + face_id] = 0;
 
      /* Turbulence */
 
      if (   (cs_glob_turb_model->itytur == 2)
          || (cs_glob_turb_model->itytur == 3)
          || (cs_glob_turb_model->model == CS_TURB_V2F_PHI)
          || (cs_glob_turb_model->model == CS_TURB_K_OMEGA)
          || (cs_glob_turb_model->model == CS_TURB_SPALART_ALLMARAS)) {
 
        cs_real_t uref2 =  0;
        for (cs_lnum_t ii = 0; ii < 3; ii++)
          uref2 += cs_math_pow2(rcodcl1_vel[n_b_faces*ii +face_id]);
        uref2 = cs::max(uref2, 1.e-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 (romb) 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)
 
         * Hydraulic diameter */
        const cs_real_t dhy = 0.075;
 
        const cs_real_t viscl0 = cs_glob_fluid_properties->viscl0;
 
        const cs_real_t rhomoy = bfpro_rom[face_id];
 
        cs_turbulence_bc_inlet_hyd_diam(face_id,
                                        uref2,
                                        dhy,
                                        rhomoy,
                                        viscl0);
 
        /* Handle Scalars
 
         * Enthalpy in J/kg
         * On this example we impose the value of the enthalpy
         * the arbitrary value of 1.d6 corresponds to a temperature
         * of 2200 Kelvin for argon at atmospheric pressure (see dp_ELE) */
 
        th_f->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
        th_f->bc_coeffs->rcodcl1[face_id] = 1e6;
 
        /* Specific model for Electric arc:
         * Vector Potential: Zero flux by default beacuse we don't a lot about
         * vector potential (what we know, is that A is equal to zero at
         * and infinite distance)
 
         * All the boundary conditions for A are zero flux, except on some
         * selected faces where we need to impose a value in order to have
         * a stable calculation (well defined problem).
         * These faces are selected where we are sure that the electrical current
         * density remains very low generally far from the center of the electric
         * arc and from the electrodes (see above) */
 
        if (potva != nullptr) {
          potva->bc_coeffs->icodcl[face_id] = CS_BC_NEUMANN;
          potva->bc_coeffs->rcodcl3[face_id] = 0;
          potva->bc_coeffs->icodcl[n_b_faces + face_id] = CS_BC_NEUMANN;
          potva->bc_coeffs->rcodcl3[n_b_faces + face_id] = 0;
          potva->bc_coeffs->icodcl[n_b_faces*2 + face_id] = CS_BC_NEUMANN;
          potva->bc_coeffs->rcodcl3[n_b_faces*2 + face_id] = 0;
        }
      }
    }

Example 2

For boundary faces of zone "2" assign an free outlet and example of electrode for Joule Effect by direct conduction.

  {
    const cs_zone_t *zone = cs_boundary_zone_by_name("2");
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_OUTLET;
 
      /* Handle Scalars
 
       * Enthalpy in J/kg  (By default zero flux with CS_OUTLET)
       *    Nothing to do
 
       * Mass fraction of the (n-1) gas mixture components
       * (Zero flux by defaut with ISOLIB)
       *    Nothing to do
 
       * Specific model for Joule Effect by direct conduction:
 
       * If you want to make a simulation with an imposed Power PUISIM
       * (you want to get PUISIM imposed in cs_user_parameters.c and
       * PUISIM = Amp x Volt) you need to impose IELCOR=1 in cs_user_parameters.c
       * The boundary conditions will be scaled by COEJOU coefficient
       * for example the electrical potential will be multiplied bu COEJOU
       * (both real and imaginary part of the electrical potential if needed)
 
       * COEJOU is automatically defined in order that the calculated dissipated
       * power by Joule effect (both real and imaginary part if needed) is equal
       * to PUISIM
 
       * At the beginning of the calculation, COEJOU ie equal to 1;
       * COEJOU is writing and reading in the result files.
 
       * If you don't want to calculate with by scaling,
       * you can impose directly the value. */
 
      if (cs_glob_physical_model_flag[CS_JOULE_EFFECT] > -1) {
        potr->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
        if (ielcor == 1) {
          const cs_real_t coejou = cs_glob_elec_option->coejou;
          potr->bc_coeffs->rcodcl1[face_id] = 500*coejou;
        }
        else {
          potr->bc_coeffs->rcodcl1[face_id] = 500.;
        }
      }
 
      if ( cs_glob_physical_model_flag[CS_JOULE_EFFECT] > 1) {
        poti->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
        if (ielcor == 1) {
          const cs_real_t coejou = cs_glob_elec_option->coejou;
          poti->bc_coeffs->rcodcl1[face_id] = sqrt(3)*500*coejou;
        }
        else {
          poti->bc_coeffs->rcodcl1[face_id] = sqrt(3)*500.;
        }
      }
 
    }

Example 3

For boundary faces of zone "2" assign a free outlet and example of anode for electric arc.

  {
    const cs_zone_t *zone = cs_boundary_zone_by_name("2");
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_OUTLET;
 
      /* Handle scalars
 
       * Enthalpy in J/kg  (Zero flux by default with ISOLIB)
       *    Nothing to do
 
       * Real component of the electrical potential */
 
      /*  For electric arc model,
       * ======================
       * We generally calculate the "electric variables" assuming that the total
       * intensity of the electrical current is imposed (COUIMP is the value of
       * the imposed total current).
 
       * In that case, you need to impose IELCOR=1 in cs_user_parameters.c
       * The "electrical variables" will be scaled by COEPOT coefficient :
       * for example the electrical potential will be multiplied by COEPOT,
       * Joule effect will be multipied by COEPOT * COEPOT and so on (see
       * cs_user_electric_scaling.c).
 
       * COEJOU is defined in cs_user_electric_scaling.c: different possibilities
       * are described in cs_user_electric_scaling.c, depending on the different
       * physics you want to simulate (scaling from current, from power, special
       * model for restriking ...).
 
       * The variable DPOT is defined: it corresponds to the electrical potential
       * difference between the electrodes (Anode potential - cathode Potential).
       * DPOT is calculated in cs_user_electric_scaling.c. DPOT is saved at each
       * time step, and for a following calculation
 
       * DPOT is the value of the boundary condition on anode assuming that
       * the cathode potential is equel to zero.
 
       * It is also possible to fix the value of the potential on the anode.
       * (for example, 1000 Volts).*/
 
      potr->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
 
      if ((cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] > -1) && (ielcor == 1))
        potr->bc_coeffs->rcodcl1[face_id] = cs_glob_elec_option->pot_diff;
      else
        potr->bc_coeffs->rcodcl1[face_id] = 1000.;
 
      /* Mass fraction of the (n-1) gas mixture components
       *  zero flux by default with CS_OUTLET
       *  nothing to do
 
       * vector Potential
       *  zero flux by default with CS_OUTLET
       *  nothing to do */
    }

Example 4

For boundary faces of zone "2" assign a wall and example of potential vector Dirichlet condition

  {
    cs_span_2d<cs_real_t> cvara_potva = cs_field("vec_potential")->get_val_v(1);
 
    const cs_zone_t *zone = cs_boundary_zone_by_name("2");
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_SMOOTHWALL;
 
      /* Wall: zero flow (zero flux for pressure)
       *       friction for velocities (+ turbulent variables)
       *       zero flux for scalars
 
       * Handle scalars
       * Enthalpy in J/kg  (Zero flux by default)
       *   Nothing to do
 
       * Real component of the electrical potential
       *   Zero flux by default
       *   Nothing to do
 
       * Specific model for Electric arc:
       * ===============================
 
       * Vector potential  A (Ax, Ay, Az)
 
       * Zero flux by default because we don't a lot about vector potential
       * (what we know, is that A is equal to zero at the infinite)
 
       * All the boundary conditions for A are zero flux, except on some chosen
       * faces where we need to impose a value in order to have a stable
       * calculation. These faces are chosen where we are sure that the
       * electrical current density remains very low generally far from the
       * center of the electric arc and from the electrodes:
 
       * On the following example, we choose to impose a "dirichlet" value for
       * the 3 components of A on a small zone of the boundary located near the
       * vertical free outlet of the computation domain.
 
       * In this example, the electric arc is at the center of the computational
       * domain, located on z axis (near x = 0 and y = 0).
       * The x (1st) and y (the 3rd) coordinates are contained between
       * -2.5 cm nd 2.5 cm:
 
       *    Ax(t, x,y,z) = Ax(t-dt, x=2.5cm, y=2.5cm, z)
       *    Ay(t, x,y,z) = Ay(t-dt, x=2.5cm, y=2.5cm, z)
       *    Az(t, x,y,z) = Az(t-dt, x=2.5cm, y=2.5cm, z) */
 
      if (cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] > 1) {
        if (   (cdgfbo[face_id][0] <= 249e-2)
            || (cdgfbo[face_id][0] >= 249e-2)
            || (cdgfbo[face_id][2] <= 249e-2)
            || (cdgfbo[face_id][2] >= 249e-2)) {
 
          const cs_lnum_t c_id = b_face_cells[face_id];
          potva->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
          potva->bc_coeffs->rcodcl1[face_id] = cvara_potva(c_id, 0);
 
          potva->bc_coeffs->icodcl[n_b_faces + face_id] = CS_BC_DIRICHLET;
          potva->bc_coeffs->rcodcl1[n_b_faces + face_id] = cvara_potva(c_id, 1);
 
          potva->bc_coeffs->icodcl[n_b_faces*2 + face_id] = CS_BC_DIRICHLET;
          potva->bc_coeffs->rcodcl1[n_b_faces*2 + face_id]
            = cvara_potva(c_id, 2);
 
        }
      }
    }

Example 5

For boundary faces of zone "51" assign a wall and restriking model for electric arc (anode boundaray condition).

  {
    cs_field_t *th_f = cs_thermal_model_field();
 
    auto th_f_val_ext = th_f->bc_coeffs->get_val_ext();
    auto th_f_h_ext = th_f->bc_coeffs->get_h_ext();
 
    const cs_zone_t *zone = cs_boundary_zone_by_name("51");
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_SMOOTHWALL;
 
      /* Enthalpy (J/kg) :
       *   imposed heat transfer coefficient */
 
      th_f->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
      th_f_val_ext[face_id] = 2.e4;
      th_f_h_ext[face_id] = 1.e5;
 
      /* Real electrical potential: anode boundary condition;
       * pot_diff calculated in cs_user_electric_scaling.c */
 
      potr->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
 
      if ((cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] > 1) && (ielcor == 1)) {
        if (   (cs_glob_elec_option->irestrike == 1)
            && (cs_glob_time_step->nt_cur <= cs_glob_elec_option->ntdcla+30)) {
          cs_real_t z1 = cs_glob_elec_option->restrike_point[2] - 2.e-4;
          if (z1 <= 0.0)
            z1 = 0.0;
          cs_real_t z2 = cs_glob_elec_option->restrike_point[2] + 2.e-4;
          if (z2 >= 2e-2)
            z2 = 2e-2;
 
          if ((cdgfbo[face_id][2] >= z1) && (cdgfbo[face_id][2] <= z2)) {
            th_f->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
            th_f->bc_coeffs->rcodcl1[face_id] = cs_glob_elec_option->pot_diff;
          }
          else {
            th_f->bc_coeffs->icodcl[face_id] = CS_BC_NEUMANN;
            th_f->bc_coeffs->rcodcl3[face_id] = 0;
          }
        }
      }
 
      /* Vector potential : Zero flux */
 
      if (potva != nullptr) {
        potva->bc_coeffs->icodcl[face_id] = CS_BC_NEUMANN;
        potva->bc_coeffs->rcodcl3[n_b_faces*0 + face_id] = 0;
        potva->bc_coeffs->rcodcl3[n_b_faces*1 + face_id] = 0;
        potva->bc_coeffs->rcodcl3[n_b_faces*2 + face_id] = 0;
      }
    }

Example 6

For boundary faces of zone "4" assign a symmetry.

  {
    const cs_zone_t *zone = cs_boundary_zone_by_name("4");
    for (cs_lnum_t ilelt = 0; ilelt < zone->n_elts; ilelt++) {
 
      const cs_lnum_t face_id = zone->elt_ids[ilelt];
 
      bc_type[face_id] = CS_SYMMETRY;
      /* For all scalars, by default a zero flux condition is assumed
       * (for potentials also)
 
       * In Joule effect direct conduction,
       * we can use an anti-symetry condition for the imaginary component of the
       * electrical potential depending on the electrode configuration: */
 
      if (poti != nullptr) {
        poti->bc_coeffs->icodcl[face_id] = CS_BC_DIRICHLET;
        poti->bc_coeffs->rcodcl1[face_id] = 0.;
      }
    }