K M Tang

J D Yan

C Chapman

M T C Fang

THREE DIMENSIONAL MODELLING OF A DC ARC PLASMA IN A TWIN-TORCH SYSTEM

Three-dimensional modelling of a dc arc plasma in a twin-torch system

INTRODUCTION

Arc plasmas have been used for the production of fine powders of metallic and non-metallic materials in the past three decades [START_REF] Boulos | Thermal Plasmas -Fundamentals and Applications[END_REF]. In conventional systems, the raw material is normally introduced at the electrode tips, which often leads to problems such as nozzle blocking by agglomerates and excessive nozzle erosion [START_REF] Vardelle | Controlling Particle Injection in Plasma Spraying[END_REF]. To overcome these operational difficulties a revolutionary twin-torch configuration was developed in the early nineteen nineties [START_REF] Williams | Twin Torch Transferred Arcs with Remote Coupling for Materials Processing and Synthesis[END_REF] and subsequently used widely for forming ultra-fine particles [START_REF] Ageorges | Synthesis of Aluminum Nitride in Transferred Arc Plasma Furnace[END_REF], waste disposal [START_REF] Iwao | Development of Twin Torch Plasma Arc for Hazardous Waste Treatment ICOPS[END_REF] and spray deposition [START_REF] Huang | Ultrafast Thermal Plasma Physical Vapor Deposition of Yttria-Stabilized Zirconia for Novel Thermal Barrier Coatings[END_REF]. The former two applications take advantage of the fact that the design of the twin-torch system eliminates the requirement that the materials to be processed have to be electrically conducting to form part of the current path, just as in the case of a conventional transferred arc. The twin torch system, which has two-electrode assemblies oriented with an included angle of the torch axes between 90 and 120 degrees, has typically two distinguishable plasma columns in the regions close to the electrodes. The columns merge in a region commonly known as the coupling zone which is several centimetres away from the electrode tips. Despite the fact that the arc columns in the region very close to the electrode tips can be justifiably assumed to be axisymmetric (as will be shown in section 4.1), the bending and subsequent merging of the two columns makes the arc shape and the distribution of flow, temperature and other plasma parameters truly three dimensional.

It is generally accepted that the plasma environment in the coupling zone attains a high purity level in the sense that contamination due to electrode erosion is hardly detectable in this region. In the process of manufacturing fine powders, the feeder material is directly fed into this hot region to allow evaporation to take place before the evaporated material is rapidly condensed in a cold environment to achieve nucleation. It can therefore be seen that the thermal and flow environment in this coupling zone exerts significant influence on the quality and efficiency of the process. The two aspects that need to be first investigated are the formation of this coupling zone and its flow and thermodynamic characteristics. Unfortunately, there have been only very limited experimental [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF][START_REF] Iwao | Heating Efficiency of Twin Torch Plasma Arc[END_REF] or computational [START_REF] Barthelemy | Modeling a Pilot-Scale Combustion/Vitrification Furnace under Oxygen Plasma Arc Transferred between Twin Torches[END_REF][START_REF] You | Numerical Investigation of Flow Field of a Dual-Jet Plasma Generator[END_REF] investigations on the characteristics and behaviour of the plasma system due to the difficulties associated with modelling and diagnosis of the 3D arc. Barthelemy et al [START_REF] Barthelemy | Modeling a Pilot-Scale Combustion/Vitrification Furnace under Oxygen Plasma Arc Transferred between Twin Torches[END_REF] reported their work on modelling an arc drawn between two parallel electrodes. However radiation from the arc is neglected and a complete set of boundary conditions for the governing equations is not given. In addition, results are given only in the form of colour contours which cannot be quantitatively interpreted. There is no comparison between prediction and measurement. The work in [START_REF] You | Numerical Investigation of Flow Field of a Dual-Jet Plasma Generator[END_REF] neglected the important role that the Lorentz force plays and the temperature and velocity boundary conditions at the arc roots are specified according to experimental results.

The objective of the present work is to develop an arc model for such a twin torch plasma system and to numerically investigate the characteristics of the arcing environment. Section 2 gives a description of the mathematical model. Section 3 is devoted to the verification of the 3D model. In Section 4, the computational results are presented and discussed. Comparison is then given between computed results and available experimental data. Finally, conclusions are drawn in Section 5.

THE ARC MODEL

A discussion on the plasma state is necessary before the governing equations are introduced. The design and typical operational conditions of such a twin-torch system expect the arc to operate in a statistically steady state regime, as is often assumed in the experimental investigation [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF][START_REF] Iwao | Heating Efficiency of Twin Torch Plasma Arc[END_REF] and modelling work [START_REF] Barthelemy | Modeling a Pilot-Scale Combustion/Vitrification Furnace under Oxygen Plasma Arc Transferred between Twin Torches[END_REF][START_REF] You | Numerical Investigation of Flow Field of a Dual-Jet Plasma Generator[END_REF]. The stability of the coupling zone is however dependent on the included angle. With a small included angle, the repulsive force between the arc columns is so strong that the shape and size of the coupling zone exhibits rapid fluctuation, as suggested in a recent paper [START_REF] Colombo | Three-Dimensional Time-Dependent Modeling of a DC Transferred Arc Twin-Torch System[END_REF] where no details of the arc model, electrode configuration, initial and boundary conditions are available. When the included angle increases, the repulsive force tends to weaken and the coupling zone becomes more stable [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. High speed video obtained for an industrial scale system shows that the main body of the coupling zone, with an included angle of 110 degrees, stays approximately in a definable space while its boundary oscillates rapidly in a manner similar to that observed at the tail of a turbulent plasma jet [START_REF] Trelles | Multiscale Finite Element Modeling of Arc Dynamics in a DC Plasma Torch[END_REF]. For convenience, this feature of oscillation of the discharge area is referred to as "wagging" in the paper. Considering the statistic behaviour of the coupling zone, it is proposed in this paper that the wagging effect be modelled by turbulence enhanced momentum and energy transfer introduced in a set of time averaged Navier-Stokes equations. Further discussion on this point will be given in the results and discussion section. The included angle chosen for this study is 90 degree [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF] for which experimental results are available.

As the plasma is generated in atmospheric environment with an arc current above 300A, local thermodynamic equilibrium (LTE) is assumed. Since both electrode assemblies are water cooled, the erosion rate of copper and tungsten were both found to be less than 0.01 parts per million (ppm) for a current range of 200A to 1000A [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. Therefore electrode erosion is neglected in the present work and the arc burns in pure argon. Gravitational effect is small in comparison with the Lorentz force or gas inertia and is neglected. For high temperature arcs laminar viscous effects are negligible [START_REF] Zhang | Theoretical Investigation of a 2ka arc in a supersonic nozzle[END_REF]. However to be consistent with the software (PHOENICS [START_REF]PHOENICS is provided[END_REF]) used to produce the results, the viscous terms are retained in the equation below without giving their detailed expression. The mass, momentum and energy equations in a Cartesian coordinate system can be written as:

( ) ( ) V S φ φ φ φ ∇ • -∇ • ∇ = v (1)
where φ is the solved-for variables as described in Table 1, ρ the density, t the time, V the velocity vector, Γ φ the diffusion coefficient, and S φ the source term. The symbols appearing in Table 1 together with those in the associated expressions are given in Table 2. The Lorentz forces in Table 1 can be expressed as

y z z y BX B J B J S - = , z x x z BY B J B J S - = and x y y x BZ B J B J S - =
. The calculation of J and B will be given later. The viscous terms are based on Newtonian fluid with Stokes hypothesis. The Ohmic heating is represented by σE 2 . The equation of state, thermodynamic, transport and electrical properties for argon are taken from [START_REF] Boulos | Thermal Plasmas -Fundamentals and Applications[END_REF]. They are all in a tabulated form and interpolation is used to calculate their values when the independent thermodynamic quantities do not fall on the data points. Radiation transfer from the arc is calculated as a volumetric energy source (q in Table 1), which is equal to the difference between emission and absorption per unit time and volume. The calculation of radiation absorption at one point in the arc involves the whole arc volume and over the frequency range of arc spectra ranging from infra-red to vacuum ultra violet. Detailed radiation transport calculation based on spectral radiation transport equation is computationally too costly. For high pressure arcs it has been found that radiation transport can be approximately accounted for by using the net emission coefficient (NEC) first introduced by Lowke [START_REF] Lowke | Radiation Emission Coefficients for sulphur hexafluoride arc plasmas[END_REF]. NEC is the difference between energy emitted (emission) and radiation energy absorbed due to radiation energy originated from other parts of the arc, in a unit volume per unit time. Thus, NEC represents the net radiation loss at a point. Emission coefficient is calculated on the axis of an infinitely long cylindrical arc column of radius R with a uniform column temperature T and pressure P. For such an arc column the radiation absorption can be converted to a one-dimensional problem [START_REF] Lowke | Radiation Emission Coefficients for sulphur hexafluoride arc plasmas[END_REF]. In our investigation q is equal to the NEC (multiplied by 4 ), (R,T,P), which has been tabulated as a function of R, T and P for argon [START_REF] Menart | Net Emission Coefficients for Argon-Iron Thermal Plasmas[END_REF]. Essoltani et. al [START_REF] Essoltani | Effect of the presence of iron vapours on the columetric emission of Ar/Fe and Ar/fe/H2 plasmas[END_REF]] show that radiation absorption in argon takes place in the first millimetre of the high temperature core of the arc. Such absorption can be accounted for by using the NEC for a cylindrical arc of radius of 1mm. Radiation re-absorption from the high temperature core outside the 1mm core where q is often negative [START_REF] Zhang | Theoretical Investigation of a 2ka arc in a supersonic nozzle[END_REF] is not significant for argon. NEC has also been used in the work of Blais et. Al [START_REF] Blais | Three-Dimensional Numerical Modelling of a Magnetically Deflected DC Transferred Arc in Argon[END_REF] to model a three dimensional d.c. transferred arc in argon. In the present investigation q is therefore computed using the data of [START_REF] Menart | Net Emission Coefficients for Argon-Iron Thermal Plasmas[END_REF] for 1mm arc radius at atmospheric pressure and the temperature of the point concerned. The computational domain used figure 1 which characteristically hereafter referred to as the comp A4, A1-C1-C4-A4, A4-A3-C3-C thermally insulating non-slip su distance (~50mm) between the temperature field inside the arc C3-B4 where the pressure is fix bigger than one. horiated tungsten electrode with a conic tip, as s d of plasma gas around the cathode and also se asma ignition. The radius of the flat cathode tip rod has a radius of 6.5mm and the thickness of quations is shown in in [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. This domain is oundaries A1-A2-A3-B4-C4 are treated as ions considering the ect on the flow and 2-C3-B3 and B1-C2-Peclet number much equations of the twin shown in figure 2. A erves as an auxiliary is 1mm and the full f the ring-shaped gas inlet is 2mm. The anode torch is similar to the cathode in structure. The anodic electrode tip is made of a high-purity copper button. The radius of the anode tip is 8mm and the thickness of the ring shaped gas inlet is 2mm. The copper nozzle around the anode is axially extended (downwards in figure 2) for 10mm from the tip of the electrode in order to better confine the anode attachment towards the centre of the tip. The torch separation, which is the distance along the straight line from the centre of the cathode tip to that of the anode, is 76mm. For such a device, the behaviour is dominated by the arc column and the bulk solution is not sensitive to the boundary conditions at the electrodes. With the volume of the arcing vessel being much greater than that of the arc, the pressure sufficiently far away from the arc is therefore considered to be constant. The background pressure is maintained at atmospheric. The boundary conditions are given in Table 2.

The continuity of current provides an equation for the calculation of electric potential:

( ) 0 = ∇ ⋅ ∇ ϕ σ (2)
where ϕ is the electric potential. The current density is calculated using Ohm's law

J σ ϕ = -∇ ur (3)
Argon at room temperature is an insulator, thus =0, which is a singularity of equation ( 2). However, if LTE is assumed for argon, argon has a very low (but not zero) electrical conductivity at room temperature [START_REF] Boulos | Thermal Plasmas -Fundamentals and Applications[END_REF]. When equation ( 2) is solved, σ is given a value of 10 -3 Ω -1 m -1 for computational convenience when the argon temperature is below 3900K. Numerical experiments show that electrical potential and electrical field within the arc plasma is not altered if a smaller value of σ is used. Such an approach has resulted in excellent agreement for an argon free burning arc [START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF]. It can also be shown that the electric field and current density in the arc column or in the region immediately surrounding it are not sensitive to the conditions specified on the boundary faces of the cubic domain because of the relatively large size of the cube. They are mainly determined by the distribution of the electrical conductivity in the conducting columns and the arc current. The electrical potential on the boundary faces is thus set to floating by specifying a zero electrical potential gradient in the normal direction and no electric current flows across these faces except where the electrodes intersect with two of the boundary faces. The nozzles surrounding the cathode and anode rods are metallic for the 90 degree configuration. In operation, they are used only for starting the arc by creating high frequency discharge in the gap between the electrode and the nozzle and do not shunt any of the arc current. It is to be noted that for electric potential equation the body of the electrodes and nozzles is part of the solution domain while it is not for the flow equations. Its electrical conductivity is set to a high value (1.0×10 5 Ω -1 m -1 ). A summary of the boundary conditions for equation ( 1) is given in Table 3.

With PHOENICS the unknowns such as pressure, enthalpy and velocities (flow variables) are iteratively computed. During the iteration for the computation of these variables electrical field is known and is equal to the value of the previous iteration. After the convergence for the flow variables is achieved, electrical conductivity is updated using the improved values of T and P. During the iteration to improve electrical potential the temperature field is frozen and the updated electrical conductivity is used. During the iteration for the electrical potential Equation ( 2) is therefore a linear equation as σ does not depend on the potential. For computational convenience, an arbitrary voltage of 5V is applied on the Cross section C in figure 2 and a zero potential on the Cross section A to obtain a current of I 5V and an electric potential field. The actual electric potential field is then obtained by multiplying a factor of I arc /I 5V to the potential field corresponding to 5V.

-A2-A3-A4, A1-C1-C4-A4, A4-A3-C3-C4, A1-A2-C2-C1, A2-B2-B3-A3, C1-B1-B4-C4 (Wall) N/A N.S. N.S. N.S. ** B2-C2-C3-B3, B1-C2-C3-B4 (Outlets) 1bar ** ** ** **
where N.S. means non-slip wall condition is used, T.I. implies a thermally insulating boundary face, ** means that the gradient of the variable concerned is zero in the normal direction of the boundary face, and N/A means boundary conditions are not necessary.

As explained in section 4.3, turbulence enhanced momentum and energy transfer is required to achieve better agreement between experimental and simulation results. Turbulence enhanced momentum and energy transfer is accounted for using the simple Prandtl mixing length model.

( )

2 t u u v v w w cD y z x z x y µ ρ ∂ ∂ ∂ ∂ ∂ ∂ = + + + + + ∂ ∂ ∂ ∂ ∂ ∂ (4)
where c is an adjustable turbulence parameter and D is a characteristic dimension of the arc column. D is given a typical value of 10mm in the present work representing the half width of the portion of the arc column (no longer cylindrical) near the coupling zone. The value of c is varied in the range from 0 to 0.2 and results are compared with the measured arc voltage in Section 4.

The magnetic vector potential, A v

, is employed to obtain the magnetic flux density through the relationship

2 0 A J µ ∇ = - ur ur (5)
B A = ∇ × ur ur (6)
where 0 µ is the permeability of free space. Equation ( 5) represents three differential equations for the three scalar components of the magnetic vector potential with the current density calculated by [START_REF] Williams | Twin Torch Transferred Arcs with Remote Coupling for Materials Processing and Synthesis[END_REF]. Because of the long range nature of the magnetic force, the domain for the calculation of the magnetic vector potential has to be larger than the CFD domain so that solution to Equation ( 5) is not sensitive to the conditions set at the extended domain boundary. The appropriate extension is shown in figure 3, where x, y and h are the dimensions of the CFD domain. In reality both electrode rods are longer than 100mm therefore the current density inside the cathode and anode rods is also extended in the computation halfway into the extended domain to better reflect the real situation. The boundary conditions for Equation ( 5) is that all its scalar components are set to zero on the bounding faces of the extended cube shown in figure 3.

Figure 3 Extended grid and current conducting area for the computation of magnetic vector potential. (x = y = h = 130mm)

VERIFICATION OF THE MODEL

Unlike the two dimensional axisymmetric arc models commonly used for modelling free burning arcs [START_REF] Chen | Heat Transfer and Fluid Flow in a High-Intensity Free-Burning Arc: an Improved Modelling Approach[END_REF] and switching arcs [START_REF] Yan | The Development of PC based CAD Tools for Auto-Expansion Circuit Breaker Design[END_REF], three dimensional arc models are more prone to numerical errors in the results. For example, it is well known that in a 2D model a free burning arc always stays axisymmetric in a cylindrical polar coordinate system. However in a 3D model there is no mechanism to automatically impose axisymmetry and the arc is free to move in any direction. Inside the arc region, Lorentz force can be very strong yet the plasma density very low. The velocity field of the arc under symmetric boundary conditions may still develop into a non-symmetric distribution if the numerical algorithms used to calculate the source terms for the momentum and energy equations are not carefully constructed. For example in PHOENICS [START_REF]PHOENICS is provided[END_REF], which is used to implement the arc models in the present work, a staggered grid system, with respect to the grid system for other scalar variables, is used for solving the velocity fields. The Lorentz force located at the cell centre of the staggered grid system has to be used in the momentum equations. It needs to be calculated from the solution of the electric and magnetic vector potential fields obtained using appropriate interpolation methods.

A free burning argon arc at 200A is modelled under identical experimental conditions as those stated in [START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF][START_REF] Haddad | Temperature Determinations in a Free-Burning Arc: I. Experimental Techniques and Results in Argon[END_REF]. The electrode and nozzle geometry is represented by a Cartesian grid system, as shown in figure 4. The cathode is made of tungsten with a rod diameter of 3.2mm and its tip has a full conical angle of 60 degree and a flat surface tip of 0.2mm. The inner and outer radius of the nozzle is 5mm and 7mm respectively. The argon flow rate is 0.5lpm. Copper is used as the anode material. At such a moderate current and the flow rate, electrode erosion is not severe and subsequently neglected in the modelling. The boundary conditions for the free burning arc are given in Table 4 with the boundaries as labelled in figure 4. Argon gas enters the nozzle at a temperature of 1123K, corresponds to the conditions used in [START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF][START_REF] Haddad | Temperature Determinations in a Free-Burning Arc: I. Experimental Techniques and Results in Argon[END_REF]. The nozzle is insulating and given a small electrical conductivity of 1x10 -3 Ω -1 m -1 for reasons stated in Section 2. The anode block (B1-B2-B3-B4 to C1-C2-C3-C4) is given an electrical conductivity of 10 5 Ω -1 m -1 . The anode surface is kept at a temperature of 1355K which is just below the boiling temperature of copper. The grid size near the anode and cathode tips is determined under the suggestions in [START_REF] Lowke | A Simplified Unified Theory of Arcs and their Electrodes[END_REF] .

x and yvelocity, u and v (m/s) z-velocity, w (m/s) Temperature,

T (K) or Enthalpy, h (J/kg) Electric Potential, (V) Cathode Top N/A N/A N/A N/A σ o j z = ∂ Φ ∂ Inlet Fix mass flux 0 W in 1123K ** C1-C2-C3-C4 (Anode) N/A N/A N/A N/A 0 = Φ B1-B2-B3-B4 (Anode) N/A N.S. Wall 1355K N/A B1-B2-C2-C1, B2-B3-C3-C2, B3-B4-C4-C3, B1-B4-C4-C1 (Anode) N/A N/A N/A N/A ** A1-A2-A3-A4, A1-A2-B2- B1, A2-A3-B3-B2, A3-A4-B4-B3, A1-A4-B4-B1 (Outlets) 1bar ** ** ** **
where J 0 is the current density obtained by dividing the arc current by the cross sectional area of the cathode top, W in is derived from the inlet conditions of the experiment, N.S. means non-slip wall condition is used, ** means that the gradient of the variable concerned is zero in the normal direction of the boundary face, and N/A indicates that boundary conditions are not necessary. In addition, the cathode and nozzle surface is treated as thermally insulating and non-slip condition applies. It is not a boundary for the electric potential equation.

A comparison of the axis temperature distribution and the radial temperature distribution between 2D and 3D modelling is given in figures 5 and 6. In the region with temperature above 10,000K, the overall level of agreement with experimental results is good, within 10% of the latter. Results from the two dimensional axisymmetric model in [START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF] have a similar level of agreement. Both models predict a bell-shaped arc column. The symmetry in radial temperature distribution is excellent while there is a very small level of asymmetry (5%) in the radial distribution of axial velocity, as indicated in figure 7.

The predicted arc voltage is 12.5V which is very close to the value predicted by others, 12.1V in [START_REF] Lowke | A Simplified Unified Theory of Arcs and their Electrodes[END_REF], 12.7 in [START_REF] Haddad | Temperature Determinations in a Free-Burning Arc: I. Experimental Techniques and Results in Argon[END_REF] and 12.9 in [START_REF] Lowke | Theory of Free-Burning Arc Columns including the influence of the Cathode[END_REF]. We thus conclude that the 3D model produces satisfactory results for known cases.

Figure 4 Computational domain for free burning arc and the grid system near the cathode. Typical cell size is 0.1mm x 0.1mm x 0.1mm.

RESULTS AND DISCUSSIONS

An included angle of 90 degree is used in the present work (figure 1) based on the consideration that at smaller angles the arc coupling zone tend to be highly unstable and the steady state approximation may not be a good representation of the physical situations. Experimental V-I characteristics are available for this configuration [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. The domain size and boundary conditions are already given in section 2. The pressure at the outlets is fixed to 1bar and the temperature at the inlets is 300K. The cubic CFD domain of 0.13m x 0.13m x 0.13m is divided into 129 x 79 x 129 cells in the x, y and z direction respectively. Another 10 cells are allocated in each of the 6 directions (positive and negative directions in each coordinate space) into which the CFD domain is extended for the computation of the magnetic vector potential. The typical cell size is 0.345mm x 0.345mm x 0.345mm in the arc region.

Since the coupling of the two arc columns takes place in the x-z plane (shown later in figure 9), the plasma regions only takes a fraction of the domain thickness in the y direction (as shown in figure 11 later on), and fewer cells are therefore needed in the y direction. To study the sensitivity of the results to the size of the domain which controls the location of the outlet, the cases were repeated with an enlarged CFD domain (0.16m x 0.13m x 0.16m in x, y and z direction). Results obtained are very similar to those with the CFD domain size of 0.13m x 0.13m x 0.13m. Thus results from the smaller domain are used in this section for discussion. Numerical results were obtained for the arc current of 300A, 500A, 700A and 900A.

Basic Features of the Coupled Torch under Laminar Flow Conditions

Previous studies suggest that laminar flow is a satisfactory assumption for free burning arcs in argon [START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF]. The laminar flow cases were therefore studied first. The results for the typical case of 500A with torch separation of 76mm are first discussed. The turbulence parameter is set to zero, corresponding to laminar flow. Discussion on the cases with different current and turbulence will be given in sections 4.2 and 4.3. As observed experimentally, the physical appearance of the arc resembles two plasma jets, directed at one another with an angle. Figure 8 shows a high speed video frame of a typical twin-torch plasma while figure 9 presents the predicted temperature contours of the 500A case under similar conditions. The latter represents the cross-sectional temperature distribution while the video frame is a result of the optical emission from the light of sight across the whole arc column within the given exposure time. Furthermore, the former only represents a snapshot of the arc column at a specific instant while the model predicts the time-averaged temperature distribution. Therefore it is likely that the arc column shape shown in figure 8 may be different from that in figure 9. Nevertheless, the basic features are observed both experimentally and numerically, as explained below.

The anode column is bent nearly 90 degrees with respect to the axis of the anode electrode body when it merges with the cathode jet although the latter is only slightly bent. The cathode column is thinner but longer. This is a direct consequence of the much stronger gas flow in the cathode jet resulting from magnetic pinch around the tip of the pointed cathode. In discussing the arc behaviour, the point where the two columns join needs to be referred to. However it is difficult to define such a point in the three dimensional space. In the following discussion, a point-of-separation is defined in the middle plane in the y direction of the cubic domain. It is on the line of equal electrical conductivity of 10 2 Ω -1 m -1 (approximately 6000K) along the inner arc edge where the curvature of the line is maximum. The whole plasma region can then be divided into three parts by a line passing through the point-ofseparation in a direction parallel to the line joining the centre of the two electrode tips, i.e. the cathode column, the anode column, and the coupled zone, as shown in figure 9. The value of 10 2 Ω -1 m -1 is chosen so that the point-of-separation is sufficiently close to the conducting column and meanwhile a reasonable part of the plasma region can be classified as the coupling zone.

The temperature at the cathode tip is highest, around 21,300K (figures 9 and 10), which is a result of the confinement of the arc root by magnetic pinch and cold gas entrainment (as clearly shown by the streamlines near the cathode in figure 10). The current density has also the highest value near the cathode tip (5×10 7 A/m 2 ), which is in contrast to the value of 5×10 6 A/m 2 near the anode where the arc root spreads over most of the tip surface. In the simulation and in the electrical characterisation experiment [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF], argon is used as the arcing gas. However in the temperature measurement part of [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF], argon is used for anode and nitrogen is used for cathode torch. The anode has a diffusive attachment with gas temperature around 12,500K which is close to the measurement reported in [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF], although a quantitative comparison is not possible due to the difference in operating conditions (current of 200A, included angle of 110 o and torch separation of 70mm).

In practice, the region where the two arc columns are coupled is expected to form an environment for the process of non-conducting materials. It is rather surprising in the present simulation to observe a "tissue" shaped linking layer (figure 11) between the two columns. This level of detail about the shape of the coupling zone has not been reported either experimentally or computationally. The layer forms a critical part of the current conducting loop. The Lorentz force creates pinch towards the middle part of the tissue and also towards the core of the two jets (figure 12). The pinch effect tends to bring cold gas towards the middle of the tissue thus adding cooling to it. On the other hand, radiation also contributes to the cooling of the tissue. On the other hand, the high current density inside the tissue produces strong ohmic heating, compensating the energy loss. The Lorentz force in the tissue-shaped linking layer is much stronger than that in the jets region (figure 12b) since the current density in the linking layer is much higher. The force as a vector at each point can be obtained by composing the force components in different directions in the three dimensional space. As long as the mass, momentum and energy can be balanced and a positive feedback mechanism does not exist, a steady arc column can be maintained.

The shape and size of the coupling zone and the velocity and temperature fields in it are important factors for heating the fine particles of raw materials. The temperature inside the coupling zone reached a highest value of 12634K for the given configuration and operating conditions (as shown in figure 9 and 10) and is in qualitative agreement to that measured in [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. A direct comparison with [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF] is not possible as the temperature measurements were performed on a mixture of argon and nitrogen at different current (200A). At the coupling zone, the predicted thickness of the linking layer with temperature above 10000K is only 2.7 mm (figure 11) while the diameter of the arc column immediately adjacent to the coupling zone is 20mm for the 500A case. Further examination of the results reveals that this is a consequence of the combined effect of "magnetic repelling" (figure 12), convective cooling (figure 10) and strong Ohmic heating (figure 13).

There are two regions where the Lorentz force is strongest. The first region is near the cathode tip and the second region is near the point-of-separation. For the arc columns near the tips of the two electrodes, it can be seen from figures 9 and 12 that the temperature and Lorentz force vector on the 10,500K isotherm has symmetric distribution with respect to the central axes of the electrodes. The interaction between the two jets is weak because of the distance involved (>50mm). Deviation from axisymmetry occurs at 20mm away from the cathode and anode tips where the distance between the arc columns originated from the cathode and anode is significantly less than 50mm. Further towards the coupling zone, it can be clearly seen from figure 12 that the repelling effect due to Lorentz force between the two arc columns starts to play an important role which renders the arc, also the magnetic and electrical fields, in the coupling zone three dimensional. The dashed line in figure 12a indicates a line in the middle plane of the y direction where the Lorentz force inside the plasma columns has its minimum magnitudes. Since the Lorentz force vectors in the x-z plane , as given in figure 12a at the temperature isotherm of 10500K (represents approximately the arc boundary as the temperature decays very fast to ambient away outside this isotherm) of the cathode and anode arc columns, are in opposite directions , the dashed line is a zero line of the magnetic field, B. It is clear that a major part of the coupling zone and the arc columns near it experiences an outward Lorentz force that tends to move the point-of-separation further away from the nozzle tips towards the outlets. At the point-ofseparation cold gas is driven into the coupling zone by momentum convection originated from the two gas jets from cathode and anode (figure 12). The two plasma jets manage to collide in spite of the repelling Lorentz force (figure 12) and form an interaction zone. The boundary of the interaction zone closer to the two electrodes has a small radius of curvature, thus producing a higher electric field than the far side (see figure 16). On reaching steady state a strong electric field is built up inside the linking layer to ensure current continuity. This is evident from the concentrated equi-potential lines shown in figure 16. This strong electric field produces intense Ohmic heating (figure 13) to balance the convective cooling driven by the Lorentz force and radiative energy loss near the point-of-separation and thus sustains the discharge. The gas near the anode is very stagnant (5.9m/s) and it only accelerates to a maximum of 52.0m/s at a location of 25mm away from the anode tip. In contrast, flow in the cathode jet reaches a maximum speed of 441m/s at 500A.

V-I Characteristics under Laminar Flow Conditions

Figure 14 shows a comparison of the predicted and measured arc voltage [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF] at different currents with the same included angle of 90 degree and torch separation of 76mm. In the electrical characterisation part of [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF] pure argon is used for both cathode and anode. Both experiment and prediction show an increasing arc voltage with current. As the current increases, the arc column attains higher temperature due to higher current density, producing a slightly broader arc (figure 15). The cathode plasma jet has considerably higher convection flow compared to the anode. The highest temperature located at the cathode tip increases from 21330K to 24766K. The maximum temperature in the coupling zone increases from 12634K to 14289K with a slightly larger size. A stronger convection flow from the cathode and stronger repulsive Lorentz force results in a more severe bent of the anode arc column at higher current. In terms of the point-of-separation location, there is no noticeable change in the z direction, but a minute gradual shift towards the negative x direction is observed. The stronger repulsive Lorentz force also produces a separated and distinct arc tails.

In figure 14, the laminar flow assumption leads to higher arc voltages (by 5V) than those from the measurements. This situation is not expected since in the model the voltage drop is only that over the arc columns and the coupling zone. The non-equilibrium layers in front of the cathode and anode, which are not included in our model, are not expected to substantially influence the overall behaviour of the twin-torch system, and the voltage drop across these two layers is in the range of 5-13V for thoriated tungsten cathode and copper anode [START_REF] Hsu | Study of the free-burning high-intensity argon arcs[END_REF][START_REF] Olsen | Thermal and electrical properties of an argon arc[END_REF]. Considering this, the predicted arc voltage should be at least 10V lower than it is under the laminar flow condition. For the laminar flow simulation at 300A, the convective flow of the arc column is too weak to overcome the repulsive Lorentz force. This has a similar effect as reducing the included angle which will reduce the stability as mentioned in [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF]. It has not been possible to obtain a converged solution for the laminar flow case at 300A, because the two columns are completely detached up to a region very close to the outlets.

The Effect of Turbulence

Further examination of the simu condition is associated with two point-of-separation, as shown in plasma system, the visual bound two arc columns near the coupl "wagging" effect. This type o momentum and energy transfer a exist. This would lead to broaden arc columns, thus shortening the indeed what the simulation resul of this type of "wagging".

Up to now, there is very little und development to turbulence. The c forced convection flow [START_REF] Jones | The physics of high power arcs[END_REF]. For to the arc instability. The waggin It enhances heat transfer to the su length model. However, Prandtl m r for current 500, 700 and 900A (all laminar e from 10000 to 25000K. ulation results reveals that the higher arc voltage arc columns which are separated by a narrow gap figure 10. In reality, as observed in a high speed dary of the coupling zone including its tail, and ling zone, are constantly subjected to a high fr of rapid spatial fluctuation has the statistical at the arc edge where large gradient of velocity, ned arc columns, and more importantly, an enhan e current conduction path and leading to a lower ts show when the concept of turbulence is used t derstanding and work on the cause of arc instabili consensus is that it is due to shear layer instability the twin torch system the magnetic driving force ng motion causes the arc voltage to fluctuate aroun urrounding region but may not be accurately mod mixing length model does give practically useful flow cases) using a e under laminar flow p of cold gas near the video of a twin-torch a small length of the requency, small scale effect of enhancing enthalpy and density nced merge of the two r arc voltage. This is to represent the effect ity and its y for arcs burning in e may also contribute nd an average value. delled by the mixing predication in arc jets burning in forced flow environment [START_REF] Zhang | Investigation of the effects of pressure ratios on arc behaviour in a supersonic nozzle[END_REF][START_REF] Yan | A comparative study of turbulence models for SF 6 arcs in a supersonic nozzle[END_REF]. In the present work we are concerned with the time averaged behaviour of the twin jet plasma system, not the cause of its instability and development into turbulence. With limited understanding and work on turbulence in arc discharges, Prandtl mixing length model is therefore adopted in the present work following previous work for arcs in convection flow [START_REF] Zhang | Investigation of the effects of pressure ratios on arc behaviour in a supersonic nozzle[END_REF].

The turbulence length scale is determined by multiplying the turbulence parameter, c, with the typical dimension of the system, D in equation [START_REF] Ageorges | Synthesis of Aluminum Nitride in Transferred Arc Plasma Furnace[END_REF]. A comparison of the temperature and electric potential distribution for the two cases with c = 0 and c = 0.2 at 500A is given in figure 16. With the introduction of turbulence, several changes can be observed. Firstly, the cathode arc column is broadened by the turbulence, which results in a small decrease of potential drop (<5V) over the section of the cathode arc column (from the cathode tip to the point-of-separation of c = 0.2). Secondly, because of the enhanced momentum and energy transfer, the small cold gas gap is no longer able to isolate the two arc columns as shown in figure 16 and the point-of-separation moves towards the two electrode tips. This results directly in two shorter arc columns. This enhanced coupling accounts for a substantial part of the arc voltage difference in the two cases. Thirdly, it is quite obvious from figure 11 that the introduction of turbulence results in a thicker (in the y direction) linking layer in the coupling zone which requires a much smaller potential drop (15V as compared with the 23V in the c=0 case) for the same current as a result of the larger cross section of the conducting linking layer. The difference in electric field strength along the anode arc column is insignificant because of the low gas flow speed. The overall shape of the twin-torch predicted by the turbulence case of c = 0.2 is closer to the image in figure 8, especially the tail of the coupling zone is narrowed in comparison with the laminar case.

The use of a turbulence parameter of 0.02 for c in equation ( 8) effectively brings the arc voltage in figure 14 to a level very close to the measurement. Considering the cathode and anode falls this parameter needs to be further increased. The predicted arc voltages corresponding to c = 0.2 are interesting in that the difference between predicted and measured arc voltages varies with current. At 300A the prediction is 10V lower than the measurement, which is qualitatively correct. The difference however increases to 35V at 900A which is far too larger than the expected cathode and anode falls and any experimental uncertainties.

By comparing the temperature distribution and arc column shape at different currents and values of c, it was found that c = 0.2 results in excessively "fat" cathode and anode columns and a much thick linking layer in the coupling zone at 700A and 900A. A possible explanation to the increased voltage difference (figure 14) is that in reality the discharge tends to be more stable at higher current as a result of the larger size of the arc columns. The use of c=0.2 is appropriate to account for the instability of the coupling zone at 300A which has two relatively thin arc columns, but would lead to excessive broadening of the plasma region at 900A. As discussed above, the asymmetric distribution of the Lorentz force, which is associated with the arc current density distribution, is one of the factors affecting the instability of the discharge. The turbulence parameter c in equation ( 4) does not need to be a constant for different currents. To better fit the measurement, the turbulence length scale need to be a decreasing function of the arc current. At the coupling zone, the conducting current path is relatively thin in comparison with the size of the arc column (figure 11), and its thickness is in the range of 2.7mm -3.1mm, for the current range of 300A to 900A and turbulence parameter of 0 to 0.02. For turbulence parameter of 0.2, the thickness is in the range of 6.0mm -8.7mm. Based on the discussion above, the thickness of the linking layer is in the range of 3-6mm. This raises the important issue that how the raw particles should be fed into the region to allow a sufficient length of the heating path within a required temperature range. In practice, additional gas inlet may be used to help feed the raw particles into the high temperature region. To track the motion of particles in the heating chamber, the particle dynamics in plasma flow needs to be considered. These are however out the scope of the present paper.

Cases with 180 Degree Included Angle

Figure 17 shows the simulation results of V-I characteristics where the two electrodes are coaxial, with a torch separation of 76mm. For all the simulations with 180 degree included angle cases, the metallic nozzle is set intentionally to a small electrical conductivity value of 10 -3 Ω -1 m -1 to avoid the situation where the nozzle becomes part of the current conducting path due to excessive broadening of the arc column. In Figure 17, the experimental results for an included angle of 175 degree [START_REF] Megy | Characterization of a Twin Torch Transferred DC Arc[END_REF] are also given for reference. For such a configuration, the addition of a small amount of turbulence (c = 0.02) brings down the arc voltage at 900A by 17% (13V) in comparison with the merely 2.5% decrease in the case with 90 degree included angle. The reduction in arc voltage is a direct consequence of the broadening of the middle section of the arc column and the velocity reduction of the flow due to turbulent viscosity. There is only a single arc column and no coupling zone can be identified.

The anode attachment remains distributed over most of the anode tip surface. At 500A, using c = 0.02 and c = 0.2 leads to negligible reduction in arc voltage, as shown in figure 18. The temperature distribution for c = 0.02 is similar to that of c = 0 (figure 18) and the two jets meet in a region near the anode and results in a ring-shaped outflow zone (although the axis velocity given in figure 19 is not negative, the shielding gas from the anode nozzle does flow in the negative axial direction). Further increasing c to 0.2 generates an anode jet (with negative velocity in figure 19 under the viscous effect of negative gas flow from the anode nozzle) and hence a smaller anode attachment (figure 18). Despite the broadening of the arc column, the arc voltage maintains approximately the same as that in the cases of c = 0 and c = 0.02 because the anode jet increases the potential drop near the anode, as indicated by more concentrated equi-potential lines at the anode. A strong anode jet is formed in the c = 0.2 case for all currents investigated from 300A to 900A with a similar shape indicated in figure 18.

The maximum axial velocity near the cathode is 840m/s at 900A and 400m/s at 500A in the laminar flow cases.

The 300A cases for c = 0 and c = 0.02 are not included in figure 19 as the anode root at this current deviates from axisymmetry and moves to the edge of the anode tip, and a discussion on the arc behaviour for this case is not directly relevant to the influence of turbulence. It is to be noted that no restrictions are imposed on the anode surface for the location of the anode attachment in the present work. At 300A the arc column is thin and the interaction between the cathode jet and gas flow from the anode nozzle may disturb the axisymmetry of the flow field and eventually lead to an asymmetric solution. A similar phenomenon on the cathode surface has been studied in [START_REF] Dabringhausen | Different Modes of Arc Attachment at HID Cathodes: Simulation and Comparison with Measurements[END_REF]. Figure 19 Axial velocity component on the geometric axis for the coaxial arrangement of electrodes at 500A (broken lines from top to bottom in increasing order of c) and 900A (solid lines in same order).

From the discussions above it can be seen that although the use of a simple Prandtl mixing length turbulence model can effectively lower the predicted arc voltage to a reasonable level at certain current, a single value of the turbulence parameter has however failed to properly model the "wagging" effect for the current range from 300A to 900A. Due to the difficulty in defining a proper turbulence length scale, D, in equation ( 4), it is fixed to a typical value of 10mm which is a representative radius of the arc columns approaching the coupling zone. It is an over-estimation for the arc columns near the electrode tips and under estimation for the dimension of the coupling zone.

Traditionally the turbulence length scale represents the dimension of the energy and momentum carrying eddies in a turbulent flow. In the present case, it should be an indication of the size of a spatial domain over which the "wagging" of the arc columns and the coupling zone takes place. In view of the fact that the assumption of laminar flow produces satisfactory results for free burning arcs in argon [START_REF] Haddad | Temperature Determinations in a Free-Burning Arc: I. Experimental Techniques and Results in Argon[END_REF][START_REF] Lowke | A Simplified Unified Theory of Arcs and their Electrodes[END_REF][START_REF] Lowke | Theory of Free-Burning Arc Columns including the influence of the Cathode[END_REF][START_REF] Zhu | A Unified Theory of Free Burning Arcs, Cathode Sheaths and Cathodes[END_REF], further improvement to the twin torch model may require the use of varying turbulence length scale depending on the arc size at different parts of the discharge region (cathode column, anode column and the coupling zone).

The non-LTE layers in front of the electrodes are not considered in the present work. The conical shape of the cathode and the converging nozzle together produces strong cold gas flow which helps contain the cathode arc root to the tip area. The maximum current density on the cathode surface is around 5.0×10 7 A/m 2 from the simulation results, a reasonable value for free burning arcs. It is thus expected that neglecting the sheath layers has limited effect on the flow field. Its possible influence on arc temperature is restricted to a small region near the cathode tip. It is possible that using orthogonal grids for a cone electrode tip can introduce local electric field enhancement near the edge of the tip. Since the majority of the current is supplied on the flat electrode tip in the cathode the electric potential and arc temperature on the axis is less affected by the steps approximating the shape of the cathode tip. For the typical case of the 180 o arrangement at 500A with c=0.02, halving the axial dimension of the cells around the cathode and anode tips leads to a change of 1000K in maximum temperature near the cathode (4%) and a 2V difference in the arc voltage (3% of the predicted 57V).

The shape of the arc column remains approximately the same. Therefore the uncertainty from using orthogonal grids to approximate the electrode shape does not significantly affect the features of the results. The discussions on the coupling of the arc jets and the conclusions based on the results in this paper will be valid.

CONCLUSIONS

The proposed three-dimensional arc model is able to reproduce the major characteristics of the twintorch plasma system observed in experiment. The predicted shape of the twin-torch is similar to that observed in high speed photography where the anode arc column is bent nearly 90 degrees under the influence of a strong cathode jet and asymmetric Lorentz force. The model correctly predicts a positive slope V-I characteristics and suggests that the coupling of the two columns is through a tissuelike linking layer which has not been reported previously. It is shown that the assumption of laminar flow fails to produce an arc column voltage that should be lower than the measurement for the 90 degree included angle cases. Turbulence has to be introduced in the model to account for the effect of the "wagging" style, small scale fluctuation of the arc columns and the coupling zone as observed in experiment. Although it is possible to lower the predicted arc voltage to a reasonable level at certain currents, a comparison of the predicted and measured arc voltage over the current range of 300A to 900A has however revealed that the use of a constant turbulence length scale is not satisfactory, as explained in section 4. Further improvement to the model may require the use of a varying length scale depending on the arc size at different parts of the discharge region.

E

  = Electric field (E x , E y and E z ) σ = electrical conductivity l k = laminar thermal conductivity t k = turbulent thermal conductivity p c = specific heat at constant pressure q = net radiation loss per unit volume

Figure 1 Figure 2

 12 Figure 1 Computational domain ( torch plasma system with an incl

Figure 5

 5 Figure 5 Predicted and measure axis temperature for a 200A free burning arc. Solid line with dimonds -present 3D model; Solid line -2D model [19]; * -measurement [22].

Figure 6

 6 Figure 6 Predicted and measure radial temperature distribution for a 200A free burning arc. Solid line with rectangles -present 3D model; Solid line -2D model[START_REF] Fang | On the Use of Langmuir Probes for the Diagnosis of Atmospheric Thermal Plasmas[END_REF]; * -measurement[START_REF] Haddad | Temperature Determinations in a Free-Burning Arc: I. Experimental Techniques and Results in Argon[END_REF]

Figure 7 (

 7 Figure 7(b) Temperature at x and y cross section located 3mm from the cathode tip. Data overlaps showing good symmetry in the x and y direction

Figure 8

 8 Figure 8 Diagram showing the image of a typical twin-torch plasma. The thinner jet originates from the cathode and the bent jet is from the anode.

Figure 9

 9 Figure 9 Temperature contour ranging from 10,000K to 21,500K (highest temperature is 21330K) for the 500A case in the x-z plane located at the middle of the domain in the y direction. The line passing through the point-of-separation is parallel to the line connecting the tip centres of the two electrodes.

Figure 10 Flow

 10 Figure 10 Flow streamlines with temperature contour in the range from 1,000K up to 21,000K and rectangle indicates the coupling zone. The point-of-separation is at the centre of the broken circle. The x-z plane is the same as that used in figure 9.

Figure 11 3D

 11 Figure 11 3D temperature contours from 10,000K to 12,500K at 500K intervals of the cases c = 0.0 (left, laminar flow) and c = 0.2 (right, turbulent flow) at 500A.

Figure 12

 12 Diagrams showing the vector and contour of magnitude of the Lorentz force. Fig. 12a: Lorentz force vector is plotted along the 10,500K iso-temperature line. The black broken line indicates the region of zero magnitude of the Lorentz force. The x-z plane is the same as that used in figure 9. Fig.12b: It shows the Lorentz force vector in the y-parallel cross sectional plane passing through the white broken line in the upper diagram.

Figure 13

 13 Figure 13 Ohmic heating in the coupling zone with maximum value of 3.50x10 10 W/m 3 . Inset shows position of the displayed area in the domain. The maximum Ohmic heating at the cathode tip is 1.38×10 11 W/m 3 .

Figure 14

 14 Figure 14 Voltage-Current characteristics for included angle of 90 degree and torch separation of 76mm. Measurement is from [7].

Figure 15

 15 Figure 15 Temperature contour common temperature scale range

Figure 16

 16 Figure 16 Temperature and electric potential distribution of the cases c = 0.0 (left, laminar flow with 105V arc voltage) and c = 0.2 (right, turbulent flow with 84V arc voltage) at 500A. The circle joined by the two white lines represents the location of the point-of-separation in the other diagram. The equal-potential lines have an interval of 1V. The thick white line divides the plasma region into three parts.

Figure 17 Figure 18

 1718 Figure 17 Voltage-Current characteristics for included angle of 180 degree and torch separation of 76mm

Table 1

 1 Transport properties and source terms for the governing equations Equation

		φ	Γ	S φ
	Continuity	1 0	0
	x-momentum	u µ l + µ t	-∂P/∂z +S BX + viscous terms
	y-momentum	v µ l + µ t	-∂P/∂y+S BY + viscous terms
	z-momentum	w µ l + µ t	-∂P/∂z +S BZ + viscous terms
	Energy (enthalpy)	h (k l +k t )/c p	σE 2 -q +dP/dt + viscous dissipation
	Table 2 Symbols for Physical Properties
	P = Pressure			h = enthalpy
	J = current density (J x , J y and J z )	B = magnetic flux density (B x , B y and B z )

Table 3

 3 Boundary Conditions

	Boundary	Pressure, P (bar)	x and y-velocity, u and v (m/s)	z-velocity, w (m/s)	Temperature, (J/kg) T (K) or Enthalpy, h	Electric Potential, (V)
	Electrode (Cathode)	N/A	N.S.	N.S.	T.I.

Value on Cross section C in figure

2

adjusted to produce the specified current Electrode N/A N.S. N.S. T.I. Value on Cross section A set

Table 4

 4 Boundary Conditions used for 3D free burning arc modelling

	Boundary	Pressure, P

Confidential: not for distribution. Submitted to IOP Publishing for peer review 24 June 2010