Electric arcs examples

Local variables to be added

integer          ifac, ii, iel
integer          idim
integer          izone,iesp
integer          ilelt, nlelt

double precision uref2, d2s3
double precision rhomoy, dhy
double precision z1   , z2
integer          ipotr, ipoti, f_id, ipotva
character(len=80) :: f_name

integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: bfpro_rom
double precision, dimension(:,:), pointer :: cvara_potva
integer :: keyvar

Initialization and finalization

Initialization and finalization is similar to that of the base examples

Example 1

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

call field_get_val_s(ibrom, bfpro_rom)
!===================

call getfbr('1', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  itypfb(ifac) = ientre

  ! Zone Number (from 1 to n)
  izone = 1

  ! Zone localization for a given face
  izfppp(ifac) = izone

  rcodcl(ifac,iu,1) = 0.d0
  rcodcl(ifac,iv,1) = 0.d0
  rcodcl(ifac,iw,1) = 0.d0

  ! Turbulence

  if (itytur.eq.2 .or. itytur.eq.3                  &
       .or. iturb.eq.50 .or. iturb.eq.60            &
       .or. iturb.eq.70) then

    uref2 = rcodcl(ifac,iu,1)**2                    &
           +rcodcl(ifac,iv,1)**2                    &
           +rcodcl(ifac,iw,1)**2
    uref2 = max(uref2,1.d-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 acess 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
    dhy     = 0.075d0

    rhomoy = bfpro_rom(ifac)

    call turbulence_bc_inlet_hyd_diam &
  ( ifac, uref2, dhy, rhomoy, viscl0, rcodcl )

  endif

  ! --- Handle Scalars

  ! Enthalpy in J/kg (iscalt)
  ! 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)

  ii = iscalt
  icodcl(ifac,isca(ii))   = 1
  rcodcl(ifac,isca(ii),1) = 1.d6

  !  Electric potential  (ipotr)
  ! (could corresponds also to the real part of the electrical potential if
  ! Joule Effect by direct conduction)

  ! In the Cathode example (electric arc applications),
  ! we impose a constant value of the electrical potential which is zero,
  ! assuming that the potential is equal to "ipotr + an arbitrary constant"
  ! (What is important for electric arc is the difference between anode and
  ! cathode potentials)

  icodcl(ifac,ipotr)   = 1
  rcodcl(ifac,ipotr,1) = 0.d0

  ! Mass fraction of the (n-1) gas mixture components

  if (ngazge .gt. 1) then
    do iesp=1,ngazge-1
      write(f_name,'(a13,i2.2)') 'esl_fraction_',iesp
      call field_get_id(trim(f_name), f_id)
      call field_get_key_int(f_id, keyvar, ii)
      icodcl(ifac,ii)   = 1
      rcodcl(ifac,ii,1) = 0.d0
    enddo
  endif

  ! Specific model for Joule effect by direct conduction:
  ! Imaginary part of the potentiel (ipoti) is imposed to zero

  if (ippmod(ieljou).ge. 2) then
    icodcl(ifac,ipoti)   = 1
    rcodcl(ifac,ipoti,1) = 0.d0
  endif

  ! 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 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
  ! (well defined problem)
  ! 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 (see above)

  if (ippmod(ielarc).ge.2) then
      icodcl(ifac,ipotva    )   = 3
      rcodcl(ifac,ipotva    ,3) = 0.d0
      icodcl(ifac,ipotva + 1)   = 3
      rcodcl(ifac,ipotva + 1,3) = 0.d0
      icodcl(ifac,ipotva + 2)   = 3
      rcodcl(ifac,ipotva + 2,3) = 0.d0
  endif

enddo

Example 2

For boundary faces of color 5 assign an free outlet and example of electrode for Joule Effect by direct conduction.

call getfbr('5', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  itypfb(ifac)   = isolib

  ! Zone Number (from 1 to n)
  izone = 2

  ! Zone location for a given face

  izfppp(ifac) = izone

  ! --- Handle Scalars

  ! Enthalpy in J/kg  (By default zero flux with ISOLIB)
  !    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 (ippmod(ieljou).ge. 1) then
    icodcl(ifac,ipotr)   = 1
    if (ielcor.eq.1) then
      rcodcl(ifac,ipotr,1) = 500.d0*coejou
    else
      rcodcl(ifac,ipotr,1) = 500.d0
    endif
  endif

  if (ippmod(ieljou).ge. 2) then
    icodcl(ifac,ipoti)   = 1
    if (ielcor.eq.1) then
      rcodcl(ifac,ipoti,1) = sqrt(3.d0)*500.d0*coejou
    else
      rcodcl(ifac,ipoti,1) = sqrt(3.d0)*500.d0
    endif
  endif

enddo

Example 3

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

call getfbr('2', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  itypfb(ifac)   = isolib

  ! Zone number (from 1 to n)
  izone = 3

  ! Zone localization for a given face
  izfppp(ifac) = izone

  ! --- 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).

  icodcl(ifac,ipotr)   = 1

  if (ippmod(ielarc).ge.1 .and. ielcor .eq.1) then
    rcodcl(ifac,ipotr,1) = pot_diff
  else
    rcodcl(ifac,ipotr,1) = 1000.d0
  endif

  ! Mass fraction of the (n-1) gas mixture components
  !  zero flux by default with ISOLIB
  !  nothing to do

  ! vector Potential
  !  zero flux by default with ISOLIB
  !  nothing to do

enddo

Example 4

For boundary faces of color 3 assign a wall and example of potential vector Dirichlet condition

call getfbr('3', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  itypfb(ifac)   = iparoi

  ! Zone number (from 1 to n)
  izone = 4

  ! Zone localization for a given face
  izfppp(ifac) = izone

  ! 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 (ippmod(ielarc).ge.2) then
    if (cdgfbo(1,ifac) .le.  2.249d-2  .or.                      &
        cdgfbo(1,ifac) .ge.  2.249d-2  .or.                      &
        cdgfbo(3,ifac) .le. -2.249d-2  .or.                      &
        cdgfbo(3,ifac) .ge.  2.249d-2      ) then
      call field_get_val_prev_v_by_name('vec_potential', cvara_potva)
      iel = ifabor(ifac)
      icodcl(ifac,ipotva)   = 1
      rcodcl(ifac,ipotva,1) = cvara_potva(1, iel)

      icodcl(ifac,ipotva + 1)   = 1
      rcodcl(ifac,ipotva + 1,1) = cvara_potva(2, iel)

      icodcl(ifac,ipotva + 2)   = 1
      rcodcl(ifac,ipotva + 2,1) = cvara_potva(3, iel)
    endif
  endif

enddo

Example 5

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

call getfbr('51', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  itypfb(ifac)   = iparoi

  ! Zone number (from 1 to n)
  izone = 5

  ! Zone localization for a given face
  izfppp(ifac) = izone

  ! ---- Enthalpy (J/kg) :
  !      imposed heat transfer coefficient

  ii=iscalt
  icodcl(ifac,isca(ii))   = 1
  rcodcl(ifac,isca(ii),1) = 2.d4
  rcodcl(ifac,isca(ii),2) = 1.d5

  ! Real electrical potential: anode boundary condition;
  ! pot_diff calculated in cs_user_electric_scaling.c

  icodcl(ifac,ipotr)   = 1

  if (ippmod(ielarc).ge.1 .and. ielcor .eq.1) then
    rcodcl(ifac,ipotr,1) = pot_diff
  else
    rcodcl(ifac,ipotr,1) = 100.d0
  endif

  ! Restriking modeling:
  ! ===================
  !    example to fit depending on the case, the geometry etc...
  !     and also in agreement with cs_user_electric_scaling.c

  if (ippmod(ielarc).ge.1 .and. ielcor .eq.1) then
    if (irestrike.eq.1 .and. ntcabs.le.ntdcla+30) then
      z1 = restrike_point_z - 2.d-4
      if (z1.le.0.d0) z1 = 0.d0
      z2 = restrike_point_z + 2.d-4
      if (z2.ge.2.d-2) z2 = 2.d-2

      if (cdgfbo(3,ifac).ge.z1 .and. cdgfbo(3,ifac).le.z2) then
        icodcl(ifac,isca(ii))   = 1
        rcodcl(ifac,isca(ii),1) = pot_diff
      else
        icodcl(ifac,isca(ii))   = 3
        rcodcl(ifac,isca(ii),3) = 0.d0
      endif
    endif
  endif

  ! Vector potential : Zero flux

  if (ippmod(ielarc).ge.2) then
      icodcl(ifac,ipotva    )   = 3
      rcodcl(ifac,ipotva    ,3) = 0.d0
      icodcl(ifac,ipotva + 1)   = 3
      rcodcl(ifac,ipotva + 1,3) = 0.d0
      icodcl(ifac,ipotva + 2)   = 3
      rcodcl(ifac,ipotva + 2,3) = 0.d0
  endif

enddo

Example 6

For boundary faces of color 4 assign a symmetry.

call getfbr('4', nlelt, lstelt)
!==========

do ilelt = 1, nlelt

  ifac = lstelt(ilelt)

  ! Symmetries

  itypfb(ifac)   = isymet

  ! Zone number (from 1 to n)
  izone = 6

  ! Zone localization for a given face
  izfppp(ifac) = izone

  ! 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 (ippmod(ieljou).ge. 2) then
    icodcl(ifac,ipoti)   = 1
    rcodcl(ifac,ipoti,1) = 0.d0
  endif

enddo