A Computational Magnetohydrodynamic Modelling Study on Plasma Arc Behaviour in Gasiﬁcation Applications

: The application of direct-current plasma arc furnace technology to the problem of coal gasiﬁcation is investigated using computational multiphysics models of the plasma arc inside such units. An integrated modelling workﬂow for the study of DC plasma arc discharges in synthesis gas atmospheres is presented. The thermodynamic and transport properties of the plasma are estimated using statistical mechanics calculations and are shown to have highly non-linear dependencies on the gas composition and temperature. A computational magnetohydrodynamic solver for electro-magnetically coupled ﬂows is developed and implemented in the OpenFOAM ® framework, and the behaviour of three-dimensional transient simulations of arc formation and dynamics is studied in response to different plasma gas compositions and furnace operating conditions. To demonstrate the utility of the methods presented, practical engineering results are obtained from an ensemble of simulation results for a pilot-scale furnace design. These include the stability of the arc under different operating conditions and the dependence of voltage–current relationships on the arc length, which are relevant in understanding the industrial operability of plasma arc furnaces used for waste coal gasiﬁcation.


Introduction
The adaptation of direct-current (DC) plasma arc furnace technology from the metallurgical industry for use in the valorisation of low-grade waste coal material is currently in development at Mintek in South Africa. Discarded coal from mining operations is currently stored in large surface mine dumps in South Africa and elsewhere and represents a significant environmental and health hazard via effects, such as mine-drainage contamination of the local ground water and soil. A means of processing large quantities of such wastes together with contaminated water into valuable intermediate energy storage products, such as synthesis gas (syngas, a mixture of carbon monoxide and hydrogen) is, therefore, of some interest as this can potentially make the remediation of waste coal dumps more economically viable.
Such processes may also be combined with carbon capture and sequestration technologies to produce hydrogen for fuel and electricity-generation applications. The advantages of using DC furnaces in the gasification process include scaling to much larger units than traditional plasma torch gasifiers and the ability to treat waste carbon sources. Waste carbon includes discarded coal containing large quantities of ash and other impurities and ultrafine coal fractions that are difficult to process in conventional applications. The concept of DC plasma furnace gasification is shown visually in Figure 1 and is described in more detail elsewhere [1]. An important element in improving the engineering of DC plasma arc furnaces for applications, such as gasification, is a deep understanding of the behaviour of the plasma arc itself. Arcs are high-temperature, high-velocity jets of ionised gas that are generated and sustained by the close coupling between fluid flow, energy and electromagnetic fields when large electric currents are passed through thermal plasmas [2,3]. Due to the extremely large driving forces and numerous instability modes, arcs may exhibit a broad range of dynamics from steady state through to oscillatory and chaotic behaviour depending on the operating parameters [3][4][5][6].
As temperatures in the arc jet can reach well in excess of 10,000 K, direct experimentation on plasma arcs is extremely challenging. Mathematical models of various types are of value as proxies to help build intuition and understanding of arc behaviour under different operating conditions. The modelling of arcs can be divided into two main tasks-estimation of the thermodynamic and transport properties of the plasma in use and the application of empirical or computational models based on the fundamental governing equations of magnetohydrodynamic (MHD) systems to calculate the arc behaviour.
Plasma properties are generally estimated using detailed statistical mechanics calculations and can have highly non-linear dependencies on temperature and pressure for even the simplest noble gases [2]. The complexity increases further when multiple gases are present, such as the CO and H 2 mixtures typical of syngas. An example of the plasma dissociation and ionisation reactions to be considered is shown below for the hydrogen species only.
In the full syngas mixture, the following species and the reactions between them are included: CO, H 2 , CH, OH, CO + , H + 2 , CH + , OH + , C, C + , C 2+ , H, H + , O, O + , O 2+ and e − . This system will be examined further as part of the present paper.
Computational modelling of arc behaviour by numerical solutions of the differential equations of fluid flow, heat transfer and electromagnetism is a challenging problem and has advanced in lockstep with increasing computer power and research into numerical methods. The earliest efforts used simplified geometries and pre-calculated electromagnetic fields in otherwise standard computational fluid dynamics (CFD) models to produce steady state solutions, such as those presented by Szekely and colleagues [7].
Soon after this, fully coupled MHD-CFD models became possible and then transient solutions, extended physics and other effects. These efforts have culminated in computational multiphysics models capable of simulating the transient dynamics of three-dimensional arcs under a wide variety of conditions [8][9][10]. Despite this progress, a number of important effects are often neglected in arc models and remain active areas of research and development today.
Among these are the formulation of accurate and mesh-independent boundary conditions able to represent non-equilibrium plasma sheath effects and the ablation of electrode and work material surfaces [11,12]. In addition, coupling between the electrode and plasma domains [13] and between the plasma and free-surface flows in the case of molten work materials [14] are significant in many cases but remain computationally expensive and algorithmically complex to implement. The convergence and computational cost of electromagnetic and thermal radiation solvers in coupled MHD-CFD models also remains a limiting step for performing high-fidelity transient simulations.
Despite the limitations on the state of the art, it was considered to be of interest to use modern plasma arc modelling tools to perform a preliminary examination of the behaviour of arcs in syngas mixtures. It is hoped that the methods presented during the course of this study may be of use for future work on arcs operating in unusual gas environments.

Model Description
As described above, the development of a practical model for plasma arcs operating in syngas mixtures requires two steps-calculation of the thermophysical properties of the plasma and numerical solutions of the governing equations.

Calculation of Plasma Properties
The basis for calculation of the material properties is a statistical mechanics description of the plasma and the assumption of a certain degree of thermodynamic equilibrium. Local thermodynamic equilibrium (LTE), in which it is assumed that particle collision processes happen very fast and a single temperature may be used to represent both heavy and light particles (molecules, atoms, ions and electrons), is frequently assumed. This is generally accepted to be a reasonable assumption in the case of atmospheric-pressure thermal plasma systems, such as plasma arcs [2], with the possible exception of the coldest outlying regions and the plasma sheaths near to the boundaries.
Under the assumption of LTE, the fundamental quantum mechanical properties of the particles present, such as their vibrational, rotational and electronic energy levels, determine their statistical partition functions, and the partition functions, in turn, define the thermodynamic quantities of a mixture, such as the chemical potential and Gibbs free energy.
Minimisation of the Gibbs free energy with respect to the particle concentrations using appropriate mass and charge balance constraints then yields an equilibrium plasma composition, from which thermodynamic material properties, such as the heat capacity C P and density ρ, may be determined [2]. For the present work, these calculations were performed using a module written for the Python programming language [15] combined with species data from the NIST Atomic Spectra and Chemistry WebBook databases [16,17].
In order to calculate transport properties, such as the viscosity µ, thermal conductivity κ and electrical conductivity σ, additional information about the particle collision processes in the plasma mixture is required. In particular, expressions for the collision cross sections and resulting collision integrals must be obtained for each particle pair and collision type. Once the collision integrals are known, Chapman-Enskog collision theory may be applied to determine the transport of various properties of interest. The procedure used in this work follows that presented by Devoto [18]. Generalised empirical and theoretical collision integrals were obtained from a variety of sources [19][20][21] and implemented for the various species present in CO-H 2 plasma mixtures.

MHD-CFD Multiphysics Model
The governing equations of a plasma arc system include the compressible Navier-Stokes and continuity equations for velocity u and pressure P, the energy transport equation for enthalpy h (related monotonically to plasma temperature T) and Maxwell's equations for electric potential φ and magnetic vector potential A.
In (4),τ is the viscous stress tensor, j is the current density vector (here representing the direction of flow of electrons, not positive charge carriers), and B is the magnetic field. In (6), k B is the Boltzmann constant, and e is the fundamental charge. Q m is the mechanical source term due to heating by pressure and kinetic energy changes in compressible flows (in OpenFOAM's formulation, this is ∂(ρK)/∂t + ∇ · (ρUK), where K is the kinetic energy of the flow field). Q r is the thermal radiation source term obtained by solution of the radiative transport Equation (the formulation for this varies depending on the radiation model in use; however, for example, in the P 1 model, it is equal to αG − 4 σ SB T 4 , where α and are the material absorptivity and emissivity, respectively, G is the radiation intensity field, and σ SB is the Stefan-Boltzmann constant). In (8), µ 0 is the magnetic permeability of the plasma medium, taken as equal to the vacuum permeability.
Coupling exists in (4)-(8) due to the presence of the Lorentz momentum source term j × B, the Ohmic heating term j · j/σ, the induced current term u × B and the strong dependence of all material properties on temperature. In the present model, the only magnetic fields present are those that are self-induced by the current flow in the arc, although external fields could easily be applied in future studies. The action of the arc current, the induced magnetic field and the resulting Lorentz force term are shown in Figure 2. The solution of the coupled set of MHD-CFD field equations is performed using the finite volume method in which the region to be modelled is decomposed into a mesh of many small control volumes, each of which enforces a problem-specific set of conservation laws with respect to its immediate neighbours. In the present study, the OpenFOAM ® open source computational mechanics framework [22] was chosen to implement a solver for the plasma arc model.
The solver uses a segregated solution of the u, P and h fields with a combined pressureimplicit splitting of operators (PISO) and semi-implicit method for pressure-linked equations (SIMPLE) algorithm to resolve the pressure-velocity coupling. An iterative segregated algorithm was implemented for the solution of the A and φ fields, combined with fast lookup tables for the material properties. Full details of the solver algorithm and implementation are available in previous work [23]. Adaptive time stepping with a Courant number limit of 0.95 was used for all simulations unless otherwise stated.

Results and Discussion
Calculations of the thermophysical properties of syngas plasmas are presented, followed by the results from MHD-CFD simulations of a small pilot-scale DC gasification furnace.

Material Properties for Syngas Plasmas
The gasification of a carbon source with water produces a mixture consisting primarily of CO and H 2 with some CO 2 and CH 4 present as contaminants. "Ideal" syngas contains only carbon monoxide and hydrogen in varying ratios. A range of 50% to 100% CO was chosen for the present study; pure CO is of interest as a reference as it is the typical gas environment in metallurgical reductive smelting processes that DC plasma arc furnaces are often used for.
At elevated temperatures, molecules become less theromdynamically stable compared with their constituent atoms, ions and electrons. A CO/H 2 plasma at moderate temperatures can, therefore, be expected to contain multiple intermediate species and decomposition products. As mentioned earlier, in these calculations, the following were included: CO, H 2 , CH, OH, CO + , H + 2 , CH + , OH + , C, C + , C 2+ , H, H + , O, O + , O 2+ and e − . As the system is assumed to be at LTE, no transport or chemical reaction kinetics are considered here. At a specified temperature, pressure and initial composition, the equilibrium mole fractions of each species (and, hence, the thermodynamic and transport properties) were calculated using the statistical mechanics methods described earlier. The results obtained for various syngas mixtures at one atmosphere pressure are shown in Figure 3. Several effects are immediately obvious and have important implications for the behaviour of arcs in syngas environments. The first is that the properties are all very strong functions of temperature, which can be expected to reinforce the coupling between the thermal energy and other fields in the MHD-CFD model. Many of the properties also vary non-linearly and non-monotonically with temperature, especially C P and κ. The degree of non-linearity increases with the increasing hydrogen content in the mixture with a high peak developing at temperatures below 5000 K; this can be expected to introduce some additional sources of instability in the cooler outlying areas of the arc.
It is interesting to observe that the changing gas composition has relatively little effect on σ except at extremely high temperatures. This suggests that the electrical behaviour of the arc would be similar over a range of syngas mixtures for a given temperature field (although, given the large differences in the thermal properties, it is likely that the temperatures will vary widely).
The total radiation emission curves show higher volumetric emission from H 2 -rich syngas plasmas at intermediate temperatures, but this trend reverses at very high temperatures. Thermal excursions in the arc may, therefore, be somewhat less likely to be brought under control by negative feedback from the radiation source term in cases where more hydrogen is present.

MHD-CFD Simulations
In order to assess the impact of plasma properties on the behaviour of arcs in syngas mixtures, a MHD-CFD model of one of Mintek's DC furnace pilot plants was developed. This particular plant will be used for future testing of the DC furnace gasification concept and is supplied by a transformer and rectifier rated at 100 kVA nominal power (in practice, this delivers a maximum usable power of between 40 and 50 kW). The furnace unit consists of a graphite-lined cylindrical crucible connected to a bottom anode with a single graphite cathode electrode mounted vertically through the refractory-lined roof. The electrode is attached to a hydraulic mount and can be moved vertically. For context, some photographs of the pilot facility are shown in Figure 4.  Table 1. Unless otherwise specified, all models were initialised with a constant temperature of 10,000 K and a velocity field of zero and were ran for a total of 200 ms simulated time. For each simulation at a given arc length L a (defined as the clearance between the electrode tip and the anode surface below) and syngas CO mole fraction x CO , the DC current was stepped down from 1000 to 200 A and back up again in steps of 200 A every 20 ms. Three-dimensional computational meshes consisting primarily of hexahedral elements were constructed for each arc length investigated. The meshes were successively refined in the central region between the electrode tip and the anode surface, where the arc jet is usually located. An example model geometry and mesh for the 0.05 m arc length case is shown in Figure 5. Boundary conditions for the various fields in the model are given in Table 2, where n and t are normal and tangent vectors at the boundary surface, j k is the cathode spot current density (taken as 2 kA/cm 2 in this work), and T lim is a numerical limit placed on the surface temperature (typically the melting or vaporisation temperature of the boundary material). The boundary conditions for A represent a magnetically-insulating situation where the modelled region is surrounded by a highly conductive material, such as metal or graphite. The gradient boundary conditions for φ describe the current density on those surfaces, while the fixed value of zero represents the electrical ground or earth potential.
The conditional boundary conditions for φ and h are described in [23] and are switched according to which portion of the cathode surface is acting as the arc attachment spot at any given time. This, in turn, is determined dynamically from the local temperature field in the plasma at each time step in the simulation. Table 2. Boundary conditions for MHD-CFD model (see Figure 5).

Field
Cathode Anode Walls Atmosphere A limited mesh-dependence study was conducted to evaluate the MHD-CFD model's sensitivity to mesh resolution. For these tests, the simulations were run for only 10 ms at a fixed current and arc length of 1000 A and 0.05 m, respectively, in order to compare the initial transient behaviour. The arc voltage as a function of time was measured as the maximum of the φ field. The results are shown in Figure 6.
The finest two mesh resolutions behave quite similarly during the early phase of the simulations and capture the transition to transient dynamics shortly before 1 ms well. The dynamics of the finest three resolutions are similar during the latter stages of the simulations and show the system settling toward irregular oscillatory behaviour, whereas the coarsest resolution produces lower and steadier voltages. This is borne out by the quantitative results in Table 3, which show that, while the models are still not entirely mesh-independent, the finest three resolutions are well within one standard deviation interval of each other, and the results may, therefore, be considered as mesh-insensitive if not mesh-independent.
It is important to note that the MHD arc system is a chaotic transient flow, and only statistical comparisons are of value here-the exact trajectory of each case is likely to diverge exponentially after the initial conditions have decayed. The 1 mm resolution was chosen for the remaining simulations as a reasonable compromise of model performance and accuracy.  Although a rigorous model validation exercise must await data from future experiments on waste coal gasification using the 100 kVA test furnace or other facilities, a preliminary comparison was conducted against electrical models of DC plasma arcs originally developed by Bowman [3]. From observations of arcs in the 1-10 kA range, Bowman proposed a scale-invariant shape for the conducting core of the arc, which includes constriction effects near to the cathode attachment. By assuming a representative average plasma conductivity σ avg over the conducting core, the empirical shape function may then be integrated to obtain a relationship between the electrical variables (9).
Here, f , the integrated cell constant of the empirical arc shape, consists of several exponential decay terms and a linear term that depends only on L a [24].
It is important to note that the Bowman model applies only to steady, verticallyoriented arc jets. If any instability in the arc or asymmetry in its direction is present, the results of the two models cannot be easily compared. Validation was, therefore, performed using the parameters most likely to satisfy these conditions-the shortest arc length, 0.01 m, and the most stable plasma gas mixture, pure CO. The results of the analysis, using the average voltage in each current period, are shown in Figure 7.
The qualitative agreement is good, and the non-linear shape of the Bowman curve is captured well by the MHD-CFD model. Further work on validation against experimental data is, however, strongly advised to confirm that the computational model is working as expected. In order to explore the parameter space of the problem, a set of 30 simulations was then performed covering a range of different arc lengths and syngas compositions as indicated in Table 1. As described earlier, the simulations also varied the current to the arc in discrete steps every 20 ms. Two such current sweeps were performed during each simulationdecreasing from 1000 to 200 A and then increasing back up to 1000 A. Quantitative data in the form of the arc voltage were sampled from the model during each run. An example of this is shown in Figure 8. It is interesting to observe the noticeable changes in the arc dynamics during each current interval, in particular the transition from regular oscillations to a chaotic regime as the current increases. In order to further quantify the results sections of the voltage time series, representative intervals of fixed current were sampled to obtain a representative average voltage and were also processed by Fourier transform. The power spectrum thus obtained was then analysed to find its peaks.
The peaks' magnitudes were taken to be indicative of the strength of the dynamics in the system, in particular if the arc was exhibiting steady or unsteady behaviour. An example of this is shown in Figure 9. Note that, since each current level is visited twice during a simulation, the pair of results from each level were averaged. It is interesting to observe that, in this case, the arc has not chosen the shortest path between anode and cathode (vertically downward from the electrode tip) but is, instead, offset to one side. The strong turbulent jet creates very asymmetric flow and heat transfer patterns in the furnace vessel-at least over the short time periods modelled here.
The distribution of plasma species in the furnace vessel shows that, due to the high operating temperatures of the arc and the low dissociation energy of H 2 , the bulk of the hydrogen is, in fact, present as individual atoms. CO only dissociates in the core of the arc jet, producing significant quantities of free C and O atoms in the same region. The most easily formed ion (C + ) is only present at levels of a few percent even in the hottest parts of the arc, suggesting that the plasmas in syngas arcs are not heavily ionised.  The average voltages calculated from the MHD-CFD results at each x CO , L a and I value were then converted into furnace power by multiplying with the current to assess the operating regimes in which the model furnace is practically operable. The results are shown for all simulations in Figure 12. The axes of each graph refer to the operational variables of the furnace, and a separate graph is drawn for each combination of syngas chemistry.
Using a power limit of 70 kW, it can be seen that arcs in syngas mixtures will generally be limited to somewhat shorter arc lengths and/or lower currents when compared to the pure CO case (x CO = 1.0); however, the effect is not particularly strong and is unlikely to prohibit operation entirely. A furnace that is operable with a CO atmosphere should, therefore, be operable with syngas mixtures, at least from a power-input perspective. In addition to the power input, the stability of the arc was examined using spectral analysis for each condition as described earlier. The results were then categorised in terms of "steady state" (no appreciable peaks in spectrum), "fully transient" (peaks observed during both the falling and rising current stages of each simulation) and "partially transient" (peaks observed during only one of the stages).
Transient behaviour can result in spikes of high arc voltage and may cause arc extinction in an operating furnace, thereby limiting its practical use. Figure 12 shows that the steady state window, which covers most of the parameter space for pure CO, retreats consistently as increasing amounts of H 2 are introduced into the system. For x CO ≤ 0.7, almost none of the simulations show steady state behaviour at any combination of current and arc length. This agrees with earlier observations on the plasma properties, which suggested that the additional non-linearities from the hydrogen reactions could lead to more instabilities in plasma arcs.

Conclusions
A modelling workflow for studying DC plasma arcs operating in syngas mixtures was successfully developed and demonstrated using simple test cases. The calculation of the plasma properties showed that syngas mixtures become increasingly non-linearly dependent on temperature as the fraction of hydrogen increases, in particular the thermal properties, such as C P and κ. This is anticipated to create an additional source of instability in the tightly coupled MHD system of plasma arcs.
A computational MHD-CFD model was then applied to the problem of arc behaviour in a small pilot-scale DC plasma arc furnace. Preliminary numerical testing showed that the model was in reasonable agreement with the empirical descriptions of plasma arcs from the literature. The parameter space of the problem was explored using simulations at a wide range of arc currents, arc lengths and syngas compositions. The results showed that, while increasing the fraction of hydrogen did not have a large impact on the electrical power delivered by the furnace, it did cause the arc behaviour to become more transient and unstable. This effect potentially limits the operability window of DC furnaces used for syngas generation as such furnaces may be restricted to undesirable short-arc, low-current, highvoltage operation for a given power level; this can potentially make furnace control more sensitive and difficult.
Naturally, much work still remains to be performed in this area. In particular, validation of the model remains at an early stage and should be revisited once additional measured data from pilot experiments are available. In particular, accurate measurements of the relationship between voltage, current and arc length should be taken at various syngas compositions and compared with the model predictions in terms of both the dynamic behaviour and quantitative values. The design and implementation of more physicallyrealistic boundary conditions, especially at the anode and cathode conducting surfaces, would be of value in improving the accuracy of the model and permitting true mesh independence.
A deeper exploration of the plasma arc parameter space with computational simulations in order to develop reduced-order pragmatic or data-driven models of arcs in syngas mixtures will also be of interest. Ultimately, it is hoped that a combination of tools, such as those presented here, can provide practical and accessible engineering guidance in the design, operation and optimisation of arc furnaces and other units using similar technology in the future.