**3. Mathematical model**

inside the nozzle depends strongly on its density, which in turn depends on the temperature and, in the high velocity, compressible regime, on the velocity itself. As inside the nozzle the radial temperature profile is very sharp, small variations of its shape, for instance its peak width, lead to appreciable changes in the radial mass distribution that reflect directly on the velocity. All this is more clearly seen and quantified using simple 1-D, two zone models in which only axial dependence is considered of an inner hot zone and an outer cold region. By assuming sonic flow at the nozzle exit, these simple 1-D models are able to predict quite reasonably the measured temperatures (within the experimental uncertainty), but give very

The purpose of this work is to present a validation of the most frequently used plasma cutting torch models employing not only temperature but also velocity values as the experimental data to be confronted with. In order to do this, a 2-D model (similar to those proposed in the literature) was developed and applied to the same 30 A oxygen cutting torch that was used in a previous velocity measurement experiment. In this paper the plasma torch, model assump‐ tions, governing equations, boundary conditions and the physical details of the model are presented. The calculated distributions of temperature and velocity and its comparison with

The high energy density cutting torch used in this numerical study consisted of a cathode centered above an orifice in a converging-straight copper nozzle. The cathode was made of copper (7 mm in diameter) with a hafnium tip (1.5 mm in diameter) inserted at the cathode center. A flow of oxygen gas cooled the cathode and the nozzle and was also employed as the plasma gas. The gas passed through a swirl ring to provide arc stability. The nozzle consisted in a converging-straight bore (with a converging length of 1 mm, and a bore 1 mm in diameter, 4.5 mm length) in a copper holder surrounding the cathode (with a separation of 0.5 mm between the holder and the cathode surface). To avoid plasma contamination by metal vapors from the anode (usually the work piece to be cut), a rotating steel disk with 200 mm in diameter and 15 mm thickness was used as the anode (Freton et al., 2002). In this study, the disk upper surface was located at 6 mm from the nozzle exit. The arc was transferred to the edge of the disk, and the rotating frequency of the disk was equal to 24 Hz. At this velocity, a wellstabilized arc column was obtained, and the lateral surface of the anode disc was completely not melted. Thus, practically no metal vapors from the anode were present in the arc. A scheme

By performing a small orifice (1 mm in diameter) on the lateral of the cathode surface the pressure in the plenum chamber (*pch*) was measured by connecting a pressure meter at the upper head of the cathode. The gas mass flow (*dm/dt*) injected in the torch was also registered. In this experiment, the arc current, the plenum pressure and the gas mass flow, were fixed to

of the torch indicating several geometric dimensions is presented in Fig. 1.

values of 30 A for the arc current, *pch* = 0.7 MPa and *dm/dt* = 0.71 g s-1, respectively.

imprecise values (100% off) of the flow velocity (Kelly et al., 2004).

the experimental data are shown.

68 Computational and Numerical Simulations

**2. Arc cutting torch**

### **3.1. Computational domain**

The schematic of the modelled domain for the simulation is presented in Fig. 2. AC is the cathode part, including a 0.75 mm radius hafnium insert (AB). DE represents the copper nozzle. The edge of the domain EF is located at a radius of 10 *RN* (Gonzalez-Aguilar at al., 1999). FG represents the anode located at 6 mm from the nozzle exit. A mass flow rate of 0.71 g s-1 with a vortex injection that leads to a ratio of the azimuthal to the axial inlet velocity of tan (13º) ≈ 0.23 was used at the torch inlet CD.

**Figure 2.** Cutting torch computational domain.

### **3.2. Model assumptions**

The most frequently used cutting torch models use the LTE approximation, with the plasma flow in chemical equilibrium, and the internal energy of the fluid being characterized by a single temperature *T*. Other assumptions are:

**a.** The plasma flow is two-dimensional and axisymmetric. There are two effects that could break these assumptions. The first one is catastrophic: the double-arcing phenomenon in which a secondary arc appears from the cathode to anode passing across the nozzle (Prevosto et al., 2009c). The second one is a non-symmetric alignment of the torch. Both effects are not considered in this work.

Momentum conservation

Internal energy conservation

r

¶

r

¶

r

emission coefficient (NEC) (Naghizadeh-Kashani, et al., 2002).

and azimuthal –*uz*– components), *p* the pressure, *δ*¯

The first is the current continuity equation

and the second is one of Maxwell's equations

t

 m

permeability of free space.

In (3) the stress tensor is given by

where

( ) ( ) 0.

( ) ( ) <sup>4</sup> 0,

where *<sup>ρ</sup>* represents the total mass density, *u*¯ the fluid velocity (having, axial –*ux*–, radial –*uy*–

current density, *<sup>B</sup>*¯ the magnetic field (only the azimuthal component was considered), *<sup>e</sup>* the internal energy, *q*¯ the total heat flux, *E*¯ the electric field and *<sup>ε</sup>N* the plasma radiation net

Two further equations are required to describe the electromagnetic part of the plasma model.

*J* =- Ñs f

<sup>0</sup> Ñ´ = *B J* m

, ,

*ij e i j*

where *σ* is the electric conductivity, *ϕ* is the electrostatic potential and *μ*0 the magnetic

æ ¶ ¶ ¶ ö = +- <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> ¶¶ ¶ <sup>÷</sup> <sup>è</sup> <sup>ø</sup> *j i l*

*j i l u u u*

<sup>2</sup> , <sup>3</sup>

d

+ Ñ× + - × + Ñ× + <sup>=</sup> ¶ *<sup>N</sup>*

*<sup>e</sup> <sup>p</sup> ue q J E p u*

*<sup>u</sup> u u p JB*

 dt

*<sup>t</sup>* (3)

*ATM*

Ñ× = *J* 0, (5)

, (6)

, (7)

*xx x* (8)

¯ the stress tensor, *<sup>J</sup>*

Numerical Modelling of a Cutting Arc Torch http://dx.doi.org/10.5772/57045

¯ the

71

¯ the identity tensor, *τ*¯

p e

*<sup>t</sup> <sup>p</sup>* (4)

+ Ñ× Ä + - - ´ =

r

¶


### **3.3. Governing equations**

The fluid part of the thermal plasma model can be expressed as a set of general transport equations expressed in conservative form as a balance among accumulation, net flux and production, namely:

$$\frac{\partial \mathcal{\partial} \boldsymbol{\nu}}{\partial t} + \nabla \cdot \overline{\boldsymbol{f}\_{\boldsymbol{\nu}}} - S\_{\boldsymbol{\nu}} = \mathbf{0},\tag{1}$$

where Ψ is a conservative quantity, *t* represents time, *f* ¯ *<sup>ψ</sup>* is the total (i.e., convective plus diffusive) flux of Ψ, and *S*Ψ is the net production rate of Ψ. The set of conservation equations describing such a flow can be expressed as follows.

Total mass conservation

$$\frac{\partial \lrcorner \rho}{\partial t} + \nabla \cdot (\rho \,\overline{u}) = 0. \tag{2}$$

Momentum conservation

**a.** The plasma flow is two-dimensional and axisymmetric. There are two effects that could break these assumptions. The first one is catastrophic: the double-arcing phenomenon in which a secondary arc appears from the cathode to anode passing across the nozzle (Prevosto et al., 2009c). The second one is a non-symmetric alignment of the torch. Both

**b.** At atmospheric pressure or above, the plasma is generally collision dominated: the mean free path for all species is much smaller than the macroscopic characteristic lengths. Therefore, the continuum assumption is valid; and the plasma is considered as a Newto‐

**c.** The plasma gas is assumed to be pure oxygen in LTE with negligible concentration of metals in the plasma. The infiltration of metal atoms into the plasma can be due to the evaporation of copper and hafnium from the cathode and the nozzle. Atoms from the work-piece could also be diffused into the plasma, but generally in a negligible concen‐

**d.** Non turbulent fluctuations are considered in the electromagnetic parameters, which results in a considerable simplification to the problem. So, for the calculations of the

**h.** The anode was considered as a porous free boundary characterized by its electrostatic

The fluid part of the thermal plasma model can be expressed as a set of general transport equations expressed in conservative form as a balance among accumulation, net flux and

diffusive) flux of Ψ, and *S*Ψ is the net production rate of Ψ. The set of conservation equations

( ) 0.

r

¶ + Ñ× =

y y

 + Ñ× - = ¶ *f S*

0,

¯

*<sup>t</sup>* (1)

¶ *<sup>u</sup> <sup>t</sup>* (2)

*<sup>ψ</sup>* is the total (i.e., convective plus

electromagnetic parameters the mean plasma state parameters are considered.

effects are not considered in this work.

70 Computational and Numerical Simulations

potential.

production, namely:

Total mass conservation

**3.3. Governing equations**

nian fluid following Navier-Stokes equation.

tration considering the high mass flux rate of the working gas.

**e.** The electrodes sheath phenomena are not included in the modelling.

**g.** In the energy equation the viscous dissipation term is considered negligible.

**f.** Hall currents and gravitational effects are considered negligible.

¶y

> r

where Ψ is a conservative quantity, *t* represents time, *f*

describing such a flow can be expressed as follows.

$$
\frac{
\partial \left(
\rho \,\overline{u}
\right)
}{
\partial t
} + \nabla \cdot \left(
\rho \,\overline{u} \otimes \overline{u} + p \,\overline{\overline{\delta}} - \overline{\overline{\tau}} \right) - \overline{J} \times \overline{B} = 0.
\tag{3}
$$

Internal energy conservation

$$\frac{\partial(\rho e)}{\partial t} + \nabla \cdot \left(\rho \overline{u} \, e + \overline{q}\right) - \overline{J} \cdot \overline{E} + p \nabla \cdot \overline{u} + 4\pi \, \varepsilon\_N \frac{p}{p\_{ATM}} = 0,\tag{4}$$

where *<sup>ρ</sup>* represents the total mass density, *u*¯ the fluid velocity (having, axial –*ux*–, radial –*uy*– and azimuthal –*uz*– components), *p* the pressure, *δ*¯ ¯ the identity tensor, *τ*¯ ¯ the stress tensor, *<sup>J</sup>* ¯ the current density, *<sup>B</sup>*¯ the magnetic field (only the azimuthal component was considered), *<sup>e</sup>* the internal energy, *q*¯ the total heat flux, *E*¯ the electric field and *<sup>ε</sup>N* the plasma radiation net emission coefficient (NEC) (Naghizadeh-Kashani, et al., 2002).

Two further equations are required to describe the electromagnetic part of the plasma model. The first is the current continuity equation

$$
\nabla \cdot \overline{J} = 0,\tag{5}
$$

where

$$
\overline{J} = -\sigma \nabla \phi,\tag{6}
$$

and the second is one of Maxwell's equations

$$
\nabla \times \overline{B} = \mu\_0 \,\overline{J},
\tag{7}
$$

where *σ* is the electric conductivity, *ϕ* is the electrostatic potential and *μ*0 the magnetic permeability of free space.

In (3) the stress tensor is given by

$$
\pi\_{i,j} = \mu\_e \left( \frac{\partial \boldsymbol{u}\_i}{\partial \boldsymbol{\alpha}\_j} + \frac{\partial \boldsymbol{u}\_j}{\partial \boldsymbol{\alpha}\_i} - \frac{2}{3} \boldsymbol{\delta}\_{i,j} \frac{\partial \boldsymbol{u}\_l}{\partial \boldsymbol{\alpha}\_l} \right), \tag{8}
$$

where *μe* is the effective viscosity and the 2/3 factor in the fluid dilatation term, ∂ *ul* / ∂ *xl* , comes from the Stokes hypothesis for the dilatation viscosity. The total heat flux in (4) describes the heat transported by conduction and the enthalpy transport by mass diffusion, and is defined as

$$
\overline{q} \equiv -\kappa\_{\mathfrak{e}} \nabla T + \overline{\Gamma}\_{\mathfrak{e}} h\_{\mathfrak{e}},\tag{9}
$$

**3.5. Turbulent model**

The closure of the system equations requires extra relationships (which are commonly known as the turbulence model) to calculate the turbulent enhanced viscosity and thermal conduc‐ tivity. There are in the literature turbulence models of varying degrees of sophistication: the Reynolds stress model (RSM), the *k*- ε model and its variants, the renormalization group (RNG) *k*-ε model, the RNG *k*- ε model taking into account the low Reynolds number effect, the realizable *k*- ε model and the Prandtl mixing length model (Zhou et al., 2009). The two most

Following previous published models (Freton et al. 2002; Freton et al., 2003), the Prandtl mixing

where *c* is an adjustable parameter and λ is a local thermal radius which characterizes the boundary of the high velocity arc core, defined as the radial distance from the axis to the point at 2000 K (Yan et al., 1999). As inside the nozzle the radial temperature profile is very sharp (Prevosto et al., 2009c), λ results very close to the nozzle radius *RN*. It has been found that for transferred arcs the turbulent Prandtl number can be approximated by unity (Fang et al., 1994).

thus only the parameter *c* in (13) needs to be adjusted by comparing the numerical results with the experiment. The turbulent viscosity for isotropic turbulence was calculated taking into

> . <sup>æ</sup> <sup>ö</sup> æ ö ¶ æ ö ¶ æ ö <sup>=</sup> <sup>ç</sup> <sup>+</sup> <sup>÷</sup> ç ÷ ç ÷ ç ÷ <sup>ç</sup> ¶ ¶ <sup>÷</sup> è ø è ø è ø <sup>è</sup> <sup>ø</sup> *<sup>x</sup> <sup>z</sup> t m <sup>u</sup> <sup>u</sup> <sup>l</sup> <sup>y</sup>*

*y yy*

Table 1 summarizes the prescribed values of the physical quantities (or its spatial derivatives) on the boundaries shown in Fig. 2. In addition, the voltage drop between the cathode AC and the anode FG was adjusted in order that the integrated value of the axial current density on a given section corresponds to the value of the electric current of the torch. An external source term to increase the temperature was applied at the axis of the torch AG to initiate the current. A current value of 30 A was used in this study, in according to the value used in the experi‐ ments. Also, at the hafnium insert AB the maximum value of the axial current density on the

1/ 2 <sup>2</sup> <sup>2</sup>

, *ml c* (13)

Numerical Modelling of a Cutting Arc Torch http://dx.doi.org/10.5772/57045 73

º1, *Pr* (14)

(Freton et al., 2003). Besides, the electrostatic

(15)

º l

usual models are the Prandtl mixing length model and the *k*- ε model.

length model was chosen. Such length is given as:

account the effect of the vortex injection (Freton et al., 2003)

m r

axis of the geometry was limited to ≤ 170 A mm−2

**3.6. Boundary conditions**

2

where *κe* is the effective thermal conductivity and *Γ*¯ *<sup>e</sup>* is the electron mass diffusion that can be approximated by

$$
\overline{\Gamma}\_e \approx -\frac{m}{e'} \overline{J},
\tag{10}
$$

where *e*′ is the elementary electric charge and *m* is the electron mass. Equation (10) neglects the charge transported by ions. In (9) *he* = 5 *kB T* / (2 *m*) represents the specific electron enthalpy (*kB* is the Boltzmann's constant).

The effective viscosity is

$$
\mu\_{\mathfrak{e}} = \mu\_l + \mu\_t,\tag{11}
$$

and the effective thermal conductivity is

$$
\kappa\_c = \kappa + \frac{\mu\_t C\_p}{P\_r},
\tag{12}
$$

where *Cp*, *μl* and *κ* are the plasma specific heat at constant pressure, viscosity and thermal conductivity, respectively. The turbulent Prandtl number *Pr* and the turbulent viscosity *μt* are given in the section 3.5.

#### **3.4. Thermodynamic properties and transport coefficients**

The thermodynamic properties and transport coefficients data for a pure oxygen plasma in the temperature range 300 ÷ 30000 K (with a 100 K temperature intervals) and for nine different pressure values (0.1 ÷ 10 atm) were taken from (Murphy & Arundell, 1994). The source terms in (4) account for the Joule effect, the compression work, and the radiation loses 4 *π ε<sup>N</sup>* , where *εN* was taken for a plasma radius of 0.5 mm (Freton et al., 2003). The NEC of pure oxygen for one atmosphere, has been multiplied by the factor *p* / *pATM* for other pressures (Zhou et al., 2009).

### **3.5. Turbulent model**

, comes

are

where *μe* is the effective viscosity and the 2/3 factor in the fluid dilatation term, ∂ *ul* / ∂ *xl*

º- Ñ +G k

G »- , ¢ *<sup>e</sup>*

mmm

is the elementary electric charge and *m* is the electron mass. Equation (10) neglects

,

and *κ* are the plasma specific heat at constant pressure, viscosity and thermal

m

conductivity, respectively. The turbulent Prandtl number *Pr* and the turbulent viscosity *μt*

The thermodynamic properties and transport coefficients data for a pure oxygen plasma in the temperature range 300 ÷ 30000 K (with a 100 K temperature intervals) and for nine different pressure values (0.1 ÷ 10 atm) were taken from (Murphy & Arundell, 1994). The source terms in (4) account for the Joule effect, the compression work, and the radiation loses 4 *π ε<sup>N</sup>* , where *εN* was taken for a plasma radius of 0.5 mm (Freton et al., 2003). The NEC of pure oxygen for one atmosphere, has been multiplied by the factor *p* / *pATM* for

*r C*

k k= + *t p*

**3.4. Thermodynamic properties and transport coefficients**

*e*

the charge transported by ions. In (9) *he* = 5 *kB T* / (2 *m*) represents the specific electron enthalpy

where *κe* is the effective thermal conductivity and *Γ*¯

as

approximated by

72 Computational and Numerical Simulations

(*kB* is the Boltzmann's constant).

and the effective thermal conductivity is

The effective viscosity is

where *e*′

where *Cp*, *μl*

given in the section 3.5.

other pressures (Zhou et al., 2009).

from the Stokes hypothesis for the dilatation viscosity. The total heat flux in (4) describes the heat transported by conduction and the enthalpy transport by mass diffusion, and is defined

, *e ee q Th* (9)

*<sup>e</sup>* is the electron mass diffusion that can be

*<sup>m</sup> <sup>J</sup> <sup>e</sup>* (10)

= + , *elt* (11)

*<sup>P</sup>* (12)

The closure of the system equations requires extra relationships (which are commonly known as the turbulence model) to calculate the turbulent enhanced viscosity and thermal conduc‐ tivity. There are in the literature turbulence models of varying degrees of sophistication: the Reynolds stress model (RSM), the *k*- ε model and its variants, the renormalization group (RNG) *k*-ε model, the RNG *k*- ε model taking into account the low Reynolds number effect, the realizable *k*- ε model and the Prandtl mixing length model (Zhou et al., 2009). The two most usual models are the Prandtl mixing length model and the *k*- ε model.

Following previous published models (Freton et al. 2002; Freton et al., 2003), the Prandtl mixing length model was chosen. Such length is given as:

$$J\_m \equiv \mathcal{C} \mathcal{A},\tag{13}$$

where *c* is an adjustable parameter and λ is a local thermal radius which characterizes the boundary of the high velocity arc core, defined as the radial distance from the axis to the point at 2000 K (Yan et al., 1999). As inside the nozzle the radial temperature profile is very sharp (Prevosto et al., 2009c), λ results very close to the nozzle radius *RN*. It has been found that for transferred arcs the turbulent Prandtl number can be approximated by unity (Fang et al., 1994).

$$P\_r \equiv 1,\tag{14}$$

thus only the parameter *c* in (13) needs to be adjusted by comparing the numerical results with the experiment. The turbulent viscosity for isotropic turbulence was calculated taking into account the effect of the vortex injection (Freton et al., 2003)

$$
\mu\_t = \rho J\_m^{-2} \left( \left( \frac{\partial u\_x}{\partial \, \mathcal{Y}} \right)^2 + \left( \mathcal{Y} \frac{\partial}{\partial \, \mathcal{Y}} \left( \frac{u\_z}{\mathcal{Y}} \right) \right)^2 \right)^{1/2}. \tag{15}
$$

### **3.6. Boundary conditions**

Table 1 summarizes the prescribed values of the physical quantities (or its spatial derivatives) on the boundaries shown in Fig. 2. In addition, the voltage drop between the cathode AC and the anode FG was adjusted in order that the integrated value of the axial current density on a given section corresponds to the value of the electric current of the torch. An external source term to increase the temperature was applied at the axis of the torch AG to initiate the current. A current value of 30 A was used in this study, in according to the value used in the experi‐ ments. Also, at the hafnium insert AB the maximum value of the axial current density on the axis of the geometry was limited to ≤ 170 A mm−2 (Freton et al., 2003). Besides, the electrostatic potential value of the nozzle DE was calculated so as to preserve the zero current balance at its surface (i.e., the nozzle is electrically floating).

Finally, at the interface between the plasma and the anode, in order to maintain the conserva‐ tion of the energy flux and current intensity at this boundary, the following relations (neglect‐ ing radiation) were used to calculate the local thermal and electric conductivities

$$
\left[ -\kappa \left( \frac{\partial T}{\partial \mathbf{x}} \right) \right]\_{anode} = \left[ -\kappa \left( \frac{\partial T}{\partial \mathbf{x}} \right) \right]\_{plasma} + J\_x \left( \frac{\mathsf{S}}{2} \frac{k\_B}{e} (T - T\_{anode}) + \varphi\_A + \varphi\_w \right),
\tag{16}
$$

$$
\left[ -\sigma \left( \frac{\partial \phi}{\partial \mathbf{x}} \right) \right]\_{anode} = \left[ -\sigma \left( \frac{\partial \phi}{\partial \mathbf{x}} \right) \right]\_{planar},\tag{17}
$$

mass density at each position was then calculated from these initial pressure and temperature distributions using the equation of state. The electrostatic potential was set to zero initially. No initial distribution was required for the magnetic field as it is governed by an elliptic

For the external flow, the pressure within the domain was set to the ambient value of one atmosphere initially. The initial temperature was set to the ambient value throughout the domain. The initial value of the density was then evaluated from the temperature and pressure using the equation of the state. The electrostatic potential was set to zero initially. The specific

The set of governing equation was written in conservation form and discretized in time using a Taylor series first-order accuracy. These equations were then discretized in space using the finite volume method and solved with the given boundary conditions on a 81 15 non uniform internal grid and 39 47 non uniform external grid, by using the predictor-corrector algorithm (Ferziger &Períc, 2002; Fletcher, 1991). The time-step used in the time-marching algorithm was chosen so that the Courant-Friedrich-Levy criterion was fulfilled (Ferziger &Períc, 2002; Fletcher, 1991). The calculation was stopped when the following condition was achieved

values used for these initial guesses did not impact the final converged results.

( ) ( ) ( )

*t t*

 c

c

c

generating the results to be presented in the following section.

1


sufficient. Use of lower values for the convergence criterion resulted in negligible differences

The accuracy of the calculations was tested by repeating them with a 38 × 15 internal grid and 19 × 47 external grid. The change in the plasma temperature was everywhere less than 15 %, while the changes in the axial velocity were less than 20 %. The finer grid was then used for

This section is devoted to the torch model validation (for the temperature and the velocity) by comparing the model results with experimental data on these quantities obtained in our

For the temperature, its radial profile at 3.5 mm from the nozzle exit was used. This profile was derived from electrostatic probe measurements (Prevosto et al, 2008a, 2008b, 2009a) by using a rotating Langmuir probe system; and from a Schlieren technique (Prevosto et al, 2010) by using a Z-type optical configuration with a laser as light source. Both techniques give an experimental uncertainty of ≈ 10 % in the temperature values. For the axial velocity, two axial distributions (with an experimental uncertainty of ≈ 10 %) derived from a time-of-flight

<sup>3</sup> 10 ,

the value of variable *χ* at the time *t*. This convergence criterion was found to be

*<sup>t</sup>* (18)

Numerical Modelling of a Cutting Arc Torch http://dx.doi.org/10.5772/57045 75

equation (equation (7)).

being *χ* (*t*)

in the final results.

**4. Validation of the model**

Laboratory for the same cutting torch.

here *φA*, *φw* and *Tanode* are the anode voltage drop, the anode work function and the anode temperature; respectively (Gleizes et al; 2005).


**Table 1.** Model boundary conditions

### **3.7. Numerical aspects**

The unsteady form of (1) was solved using a time-marching method (Ferziger &Períc, 2002; Fletcher, 1991). Consequently, initial conditions had to be supplied in order to complete the formulation of the problem.

For the internal flow calculations the initial pressure distribution was prescribed as a linearly decreasing function of the axial position, from a specific value at the torch inlet to the atmos‐ pheric one at the exit. The fluid temperature was set to the ambient value everywhere, while the initial fluid velocity was given taking into account the mass flow conservation. The total mass density at each position was then calculated from these initial pressure and temperature distributions using the equation of state. The electrostatic potential was set to zero initially. No initial distribution was required for the magnetic field as it is governed by an elliptic equation (equation (7)).

For the external flow, the pressure within the domain was set to the ambient value of one atmosphere initially. The initial temperature was set to the ambient value throughout the domain. The initial value of the density was then evaluated from the temperature and pressure using the equation of the state. The electrostatic potential was set to zero initially. The specific values used for these initial guesses did not impact the final converged results.

The set of governing equation was written in conservation form and discretized in time using a Taylor series first-order accuracy. These equations were then discretized in space using the finite volume method and solved with the given boundary conditions on a 81 15 non uniform internal grid and 39 47 non uniform external grid, by using the predictor-corrector algorithm (Ferziger &Períc, 2002; Fletcher, 1991). The time-step used in the time-marching algorithm was chosen so that the Courant-Friedrich-Levy criterion was fulfilled (Ferziger &Períc, 2002; Fletcher, 1991). The calculation was stopped when the following condition was achieved

$$\frac{\mathbf{x}^{(t)} - \mathbf{x}^{(t-1)}}{\mathbf{x}^{(t)}} \le 10^{-3},\tag{18}$$

being *χ* (*t*) the value of variable *χ* at the time *t*. This convergence criterion was found to be sufficient. Use of lower values for the convergence criterion resulted in negligible differences in the final results.

The accuracy of the calculations was tested by repeating them with a 38 × 15 internal grid and 19 × 47 external grid. The change in the plasma temperature was everywhere less than 15 %, while the changes in the axial velocity were less than 20 %. The finer grid was then used for generating the results to be presented in the following section.
