Numerical Modelling of a Cutting Arc Torch

Mancinelli, Beatriz. Grupo de Descargas Electricas, Departamento Ing. Electromecanica, 
Universidad Tecnologica Nacional, Regional Venado Tuerto, Laprida 651, 
Venado Tuerto (2600), Santa Fe, Argentina


Thermal plasmas
Thermal plasmas are partially or strongly ionized gases, usually developed at atmospheric pressure. In fact, thermal plasmas can be generated by many methods, such as dc electrical discharges at current intensities higher than a few A and up to 10 5 A: free burning arcs, transferred arcs, or non-transferred plasma torches; ac or transient arcs (e.g., discharge lamps), circuit-breakers, or pulsed arcs; rf and microwave discharges at near-atmospheric pressure; and laser induced plasmas after the ignition and expansion phases (Gleizes et al., 2005;Boulos et al., 1994).
In this kind of plasma, the electric field remains rather weak (less than a few kVm −1 ) and the electron number density is rather high (more than 10 20 m −3 ) so that the velocity distribution functions of all the types of particles can be considered as Maxwellian. The energy delivered to the plasma, in general by the Joule effect, is first picked up by the electrons because of their high mobility. The electrons transfer part of this energy to the heavy particles by elastic collisions. Due to the high electron number density, elastic collision frequencies are very high, so energy transfer is important and leads to an almost even distribution of the energy: the mean heavy particle energy is the same as the electron energy and there is thermal equilibrium among all the kinds of particles, allowing a single temperature to be defined for the plasma at a given position. In the hottest regions of thermal plasmas this mean kinetic energy is of the order of 1 eV, which corresponds to a temperature of the order of 10 4 K.
In thermal plasmas, the electrons are also most responsible for inelastic collisions, such as ionization, recombination, excitation, de-excitation, attachment, and detachment. Due to the large value of the electron density, the inelastic collision frequencies are high, and they tend to establish a statistical equilibrium among all kinds of particles; i.e. the populations of the excited atoms, molecules and ions tend to obey equilibrium laws, such as the Boltzmann and Saha laws. This situation requires one condition: the influence of radiation on the populating mechanisms of the various species must be negligible in comparison to the inelastic electron collisions. This condition is, in general, verified in the hottest regions of thermal plasmas, corresponding to the regions with high electron density values, but not in the outer regions (Ghorui et al., 2007). In spite of the weak role of radiation in the population, the experimental evidence shows that thermal plasmas emit strong amounts of radiation, particularly in the UV and visible parts of the spectrum, and they cannot be considered as black body emitters. This radiation emission, together with the presence of temperature and number density gradients within the plasmas inducing diffusion effects, shows that thermal plasmas are not in a state of thermodynamic equilibrium. Nevertheless, locally and as a first approximation, thermal plasmas can be considered to be in a state of local thermodynamic equilibrium (LTE), meaning that the particle number densities are related by equilibrium laws (because of high collision frequencies), whereas radiation is not in equilibrium with the particle distribution (Planck's law is not valid). More specifically, for thermal plasma to be considered in LTE, it has to accomplish the following requirements: a. The different species of the plasma (atoms, ions, electrons, molecules) share a single Maxwellian distribution, characterized by a single temperature.
b. The ratio of the electrostatic energy density to pressure has to be small enough, and the temperature must be high enough, so charge carriers equilibrate through collisions the energy gained from the electric field.
c. Collisions (but not radiation) are the dominating mechanisms for ionization and excitation, and there must be micro-reversibility among collisions. Hence, Saha equilibrium and Boltzmann distribution laws are valid.
d. Spatial variations of the plasma properties are enough smooth, so a given particle that diffuses from one location to another has sufficient time to equilibrate.

Plasma cutting torch modelling
Plasma cutting is a process of metal cutting at atmospheric pressure by an arc plasma jet, where a transferred arc is generated between a cathode and a work-piece (the metal to be cut) acting as the anode (Ramakrishnan et al., 1997). Small nozzle bore, extremely high enthalpy and operation at relatively low arc current (≈ 10 ÷ 200) A are a few of the primary features of these torches (Nemchinsky & Severance, 2006). The physics involved in such arcs is very complicated. The conversion of electric energy into heat within small volumes causes high temperatures and steep gradients. Dissociation, ionization, large heat transfer rates (including losses by radiation), fluid turbulence and electromagnetic phenomena are involved. In addition, wide variations of physical properties, such as density, thermal conductivity, electric conductivity and viscosity have to be taken into account. These factors make hopeless the possibility of an analytical solution for such thermal plasmas. In the last years numerical plasma modelling has reached a state advanced enough to be of practical use in the study of cutting-arc processes (e.g., Gleizes et al., 2005).
Plasma modelling by numerical simulation in cutting torches is a powerful tool to predict the values of the fundamental physical quantities, namely the plasma temperature, the particles concentration and the fluid velocity. These numerical codes are employed to understand the relevant physical processes ruling the plasma behavior in order to interpret the experimental results of several plasma diagnostics, and ultimately to obtain optimized designs of such devices (Colombo et al., 2008). However, no complete predictive power has been possible as yet due to the complexity and variety of the processes involved. In particular, the practical use of cutting torch codes requires the introduction of some numerical coefficient whose value has to be obtained from a comparison between the model predictions and the experiment. For these reasons, the experimental validation of such models is of primary importance.
On the other hand, experimental data of cutting arcs are hard to obtain due to the harsh ambient conditions, and also because much of the plasma flow takes place in relatively small, bounded regions inaccessible to experimental probing (Nemchinsky & Severance, 2006). In practice, most of the available experimental data are related to spectroscopic Girard et al., 2006) and probe (Prevosto et al., 2008a(Prevosto et al., , 2008b(Prevosto et al., , 2009a  Because of the above quoted reasons, the experimental validation of the existing cutting torch models has been restricted to the temperature distribution in the nozzle-anode gap, together with other global parameters easy to obtain as the arc voltage and the arc chamber pressure. That experimental validation consists in obtaining a good matching with the experimental temperature values within the experimental uncertainty, which typically amounts to 10 % of the measured value (Peters et al., 2007). However, since in a cutting torch the plasma pressure is close to the atmospheric pressure, the plasma temperature is mainly determined by the energy equation (with the dominant energy loses being due to inelastic electron collisions that produce a sort of "anchorage" in the temperature), and so it is almost decoupled from the hydrodynamic flow. Hence, the calculated temperature value does not result in practice quite sensitive to changes in the numerical coefficients of the models at least within the temperature experimental uncertainty.
In high-energy density cutting torches it appears that the flow velocity is a more sensitive variable than the temperature to the modelling details. This can be understood considering the following argument. The acceleration of the flow is driven by a mostly axial pressure gradient along the torch nozzle, established between the pre-nozzle chamber and the nozzle exit. So, given the pressure profile, the velocity of each fluid element at different positions 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 imprecise values (100% off) of the flow velocity (Kelly et al., 2004).
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 assumptions, 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 experimental data are shown.

Arc cutting torch
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 of the torch indicating several geometric dimensions is presented in Fig. 1.
By performing a small orifice (1 mm in diameter) on the lateral of the cathode surface the pressure in the plenum chamber (p ch ) 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 values of 30 A for the arc current, p ch = 0.7 MPa and dm/dt = 0.71 g s -1 , respectively.

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 R N (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.

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.
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 Newtonian fluid following Navier-Stokes equation.
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 concentration considering the high mass flux rate of the working gas. h. The anode was considered as a porous free boundary characterized by its electrostatic potential.

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: where Ψ is a conservative quantity, t represents time, f ψ 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 Internal energy conservation Two further equations are required to describe the electromagnetic part of the plasma model. The first is the current continuity equation where , s f = -Ñ J (6) and the second is one of Maxwell's equations 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
where μ e is the effective viscosity and the 2/3 factor in the fluid dilatation term, ∂ u l / ∂ x l , 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 where κ e is the effective thermal conductivity and Γ e is the electron mass diffusion that can be approximated by where e ′ is the elementary electric charge and m is the electron mass. Equation (10) where C p , μ l and κ are the plasma specific heat at constant pressure, viscosity and thermal conductivity, respectively. The turbulent Prandtl number P r and the turbulent viscosity μ t are given in the section 3.5.

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 π ε N , where ε N was taken for a plasma radius of 0.

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 conductivity. 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: , l º m l c (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 R N . It has been found that for transferred arcs the turbulent Prandtl number can be approximated by unity (Fang et al., 1994).
1, º r P (14) thus only the parameter c in (13) 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 experiments. 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).

Boundary conditions
Finally, at the interface between the plasma and the anode, in order to maintain the conservation of the energy flux and current intensity at this boundary, the following relations (neglecting radiation) were used to calculate the local thermal and electric conductivities ( ) , here φ A , φ w and T anode are the anode voltage drop, the anode work function and the anode temperature; respectively 2005).

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 atmospheric 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 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.

Validation of the model
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 Laboratory for the same cutting torch.
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(Prevosto et al, , 2008b(Prevosto et al, , 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 technique were employed. One of the distribution corresponded to light emitted from the arc central core, while the other corresponded to an averaged emission over the whole emitting section of the arc (Prevosto et al., 2009b). The arc core velocity was obtained from spectrally filtered light fluctuations measurements using a band-pass filter (475 nm) to detect light emission fluctuations emitted only from the arc axis. Figure 3 presents the radial profiles of the calculated temperature at 3.5 mm from the nozzle exit for c = 0 (i.e., laminar flow), c = 0.08 and c = 0.20. As shown, all the temperature profiles are similar, their differences being smaller than the experimental temperature uncertainty. As an example, Fig. 4 shows the comparison among the theoretical profile corresponding to c = 0.08 and the experimentally derived temperature profiles. It can be seen from Fig. 4 that the model results are in good agreement with the experimental data. The measured plasma flow velocity obtained for different axial positions (z) with and without band-pass filter, together with a visible photograph of the arc is shown in Fig. 5. The theoretical distributions of the axial velocity on the axis for the same c values presented in Fig. 3 are shown in Fig. 6. For comparison purposes, the measured values of the axial velocity corresponding to light emitted from the arc central core are also included in Fig.  6. It can be seen that the theoretical profiles are close among them at the vicinities of the nozzle exit (reflecting the well known fact of the little importance of the turbulence inside the nozzle (Gleizes et al., 2005) but soon after the scatter in the u x values is larger than those found for the temperature, reaching about 100 % at the middle of the gap. Hence, it can be concluded that the fluid velocity strongly depends on the particular value of the turbulent parameter c. On the other hand, the theoretical profile presenting the best matching with the experimental data is that corresponding to c = 0.08.   Figure 7 shows the averaged theoretical axial plasma velocity over the emitting arc cross section, defined from the arc core (axis) to the ≈ 4000 K temperature line (Prevosto et al., 2009a) for the same c values presented in Fig. 3. As before, the experimental values of u x corresponding to an averaged emission over the whole emitting section of the arc are also included in this figure. It can be seen that the theoretical profile presenting the best matching with the experimental data is that corresponding to c = 0.08.

Conclusions
The modelling of dc arc plasma torches is quite challenging because the plasma flow is highly nonlinear, presents strong quantity gradients and is characterized by a wide range of time and length scales.
In the last years numerical plasma modelling has reached a state advanced enough to be of practical use in the study of cutting-arc processes. However, a self-consistent description of the plasma starting from only macroscopic parameters (such as the geometry, current intensity, nature of the gas and type of employed materials, mass flow rate and/or some boundary conditions) has not yet been possible because of the lack of precise knowledge of some phenomena (electrode phenomena, radiation, turbulence, wall ablation, etc), which impose simplifications on the models. In particular, the practical use of cutting torch codes requires the introduction of some numerical coefficient whose value has to be obtained from a comparison between the model predictions and the experiment. For these reasons, the experimental validation of such models is of primary importance.
Assuming LTE conditions, the properties of the plasma that must be computed are the temperature, pressure and velocity fields. Numerical models for the plasma generated in cutting torches published during the last ten years have been validated using temperature data derived from spectroscopic measurements in the nozzle-anode gap. It has been shown in this work that the plasma temperature is not the most appropriate quantity to validate numerical codes since it is not quite sensitive to changes in the model numerical parameters. Instead, it has been shown that the plasma velocity appears to be a more adequate quantity to perform such validation.
In order to realize this validation to such a sensitive variable as the plasma velocity, a 2-D model similar to those proposed in the literature was developed and applied to the same 30 A high-energy density cutting torch that was used in the velocity measurements recently published by some of the authors. Within the experimental uncertainties, it was found that a Prandtl mixing length turbulent parameter c = 0.08 allows to reproduce both the experimental data of velocity and temperature. However, this value has to be taken with caution, since that c value depends on the actual torch geometry, gas type and arc current. It can also be concluded that the simple Prandtl mixing length model is appropriated to predict the plasma characteristics in low-current cutting torches.