Modeling and Simulation of Low Current Atmospheric and High-Pressure Helium Plasma Discharges

A plasma discharge in a Helium gas reactor at different pressures and at low currents (0.25–0.45 A) has been investigated by Computational Fluid Dynamic modeling coupled with the Maxwell’s equations. The results show different discharge dynamics across the pressure range (0.1–8 MPa), with an arc discharge obtained at high pressure and a low current arc discharge observed at atmospheric pressure. A large density gradient at higher pressure causes a strong natural convection effect in the reactor. This density gradient affects drastically the discharge shape and the velocity field at high pressures while at atmospheric pressure, a lower density gradient was observed resulting in a low velocity magnitude. It has been observed that the velocity magnitude is not affected by the electric current. The discharge electric potential has been calculated by considering the electrical characterization of the electrodes and numerical results have been compared with experimental results. The comparison shows a good agreement between the measured and calculated discharge electric potential at lower pressures. These devices can be used as plasma sources for wastewater treatment.


INTRODUCTION
Freshwater supply is a global challenge, and it is estimated that approximately four billion people experience a water shortage once a year [1]. Pollution reduces water quality and has the effect of reducing the availability of water that can be used for potable or agricultural applications. Pollution refers to the introduction of harmful contaminants such as organic chemicals or microbes derived from human or animal waste into freshwater sources [2]. Increasing amounts of pollution have made the removal of hazardous organic pollutants in the water a critical topic in environmental research. Reactors using plasma-based technology have been proposed to remove micropollutants in wastewater. Plasma is known to contain chemically active species and charged particles, which can be used in a wastewater treatment application [3].
The treatment of wastewater has been done via a thermal [2] and non-thermal plasma technology [3,4]. In addition, the process could be carried out indirectly (i.e., plasma discharges generated above the water), directly (i.e., plasma discharges generated in water), and using a hybrid approach. These different techniques imply that there is a need to improve the design of the plasma reactors used in wastewater treatment applications to increase its effectiveness and efficiency [2][3][4][5][6][7][8][9][10][11][12]. An understanding of the physical phenomena describing the plasma dynamics and the calculation of the discharge temperature are vital in order to optimize the design of reactors and consequently the plasma treatment.
References [6], [9] and [2] have shown that an atmospheric pressure plasma discharge in gas above water can produce radicals and reactive species, which have been demonstrated to rapidly and efficiently degrade many organic compounds, including phenols. The main mechanism of the reactions in the plasma-liquid interface is the reactive species transfer.
Some techniques have been reported that involves the removal of pollutants in water by designing reactors which operate at high temperature and high pressure. Among these, the use of supercritical water oxidation (SCW) in the degradation of organic pollutants such as phenol have been discussed [13][14][15]. These studies indicated that SCW (T > 374.15°C; p > 22.1 MPa) possesses different properties such as density, dielectric constant, viscosity in comparison with normal water. Moreover, the use of this high temperature and pressure conditions in which the supercritical water is reactive, leads to high reaction rates and reaction times ranging from few seconds to few minutes [14]. Nonetheless, the underlying issue remains the reactor design and its construction since different hypothesis are under investigation. The same application of high-pressure range in reactors is documented for the well-established wet air oxidation (WAO), where the temperatures are between 127°C and 300°C and pressures from 0.5 to 20 MPa in the liquid phase [16,17]. However, the excessively long reaction times of the WAO process (normally up to several hours, with a high degradation efficiency of the organic material seldom achieved), makes the SCW a more preferred technology.
Plasma technology has been profiled in the last 2 decades as a potential approach to the degradation of persistent organic pollutants in wastewater and other categories of micropollutants for which the conventional wastewater treatment plants are ineffective. While reactor design in plasma treatment of wastewater has undergone several modifications including the use of low current and nanoseconds power supply systems to reduce the previously high energy consumption, there remain issues like reactor throughput, and removal efficiency depending on the class of the pollutants. Some of these issues are still unresolved because of inadequate information of plasma dynamics and plasma-liquid interactions. Even though there are no studies dedicated to the use of high-pressure plasma in wastewater treatment, the understanding of the stability and dynamics using an inert gas would be beneficial to the scientific community especially when mineralization of highly concentrated wastewater like the pharmaceutical and hospital wastewater is considered along with landfill leachates.
Our research group has previously used Computational Fluid Dynamics (CFD) to investigate a Local Thermodynamic Equilibrium (LTE) plasma discharge in a Helium atmosphere reactor, operating at 8 MPa pressure. The aim of this research is to investigate the stability and the discharge dynamics across the pressure range 0.1-8 MPa, in which typically the LTE model is valid. The originality of our research is based on two premises: 1) the high-pressure conditions, 0.1-8 MPa, which is usually studied at high current (I ≥ 1 A) and rarely studied at low current (I < 1 A), 2) the low-current conditions, at which the current is between the upper limit of the glow (10 -1 A and the start of the lowpressure arc at current of 1 A).
The generation of an arc discharge at current lower than 1 A and pressures higher than 1 MPa is technically challenging due to the constraints of the Paschen law over a wide range of the product of the pressure (p) and the gap between the electrodes (d). For example, according to the literature data of the Paschen curves presented for the breakdown voltage as a function of p × d values, the breakdown voltage will be around 10 3 V at 0.1 MPa for our system, and about 10 4 V at 2 MPa [18]. Therefore, to overcome the constraints of the Paschen law, a pin-to-plane electrode configuration was adopted. This implies that the electrodes (cathode and anode) are brought into contact before the high voltage power supply is turned on, and then the mobile electrode (in this case, the anode) is gradually moved apart to generate the arc discharge at current less than 0.5 A for which the power supply unit is rated for. From a simulation viewpoint, this can be achieved by introducing an artificial electrical conductivity between 1 and 100 S/m, and an initial hot channel greater than 1 eV, which correspond to a temperature of ∼11,600 K in order to rise the resistance of the gas and allow of the ignition of the discharge before stability can be achieved in the model. Furthermore, this paper will help to identify the velocity and temperature fields as the operating pressure increases, and how this could be key to the optimization of the process with a view to redesigning the reactor for plasma wastewater treatment. The numerical results will provide information's about the discharge characteristics, which will be used in the next phase of the research in which a plasma discharge above a liquid medium is considered.
Mathematical Model of this paper discusses the governing Maxwell equations to describe the electromagnetic behavior of the Helium gas, the simulation setup, and the boundary conditions. In Results of the paper, the results when 1) the pressure is varied from 2-8 MPa and the electric current is fixed at 0.35 A, and 2) atmospheric pressure (0.1 MPa) and varying electric current (0.25-0.45 A) are presented. Finally, in Discussion of the paper, conclusions are drawn based on the outcome and discussion of the results.

MATHEMATICAL MODEL Assumptions
The time-dependent mathematical model is based on the following main assumptions: • The Helium plasma is considered as a single continuous fluid. • The plasma is optically thin and modelled using the Local Thermodynamic Equilibrium (LTE) model. The thermodynamic properties and transport coefficients are based on those provided by EquilTheTA [19][20][21][22][23][24]. • The gravitational force is included and acts in the negative x direction as per the coordinate system in Figure 1. • The Helium gas flow is turbulent.

Governing Equations
The fundamental laws which govern the mechanics of fluids are the continuity, momentum, and energy equations. When integrated over a finite control volume, the continuity, momentum, and energy equations are given by Equations 1-3 respectively.
⊗ Outer product, ρ Density (mass per unit volume), σ Stress tensor, a Area vector, E Total energy per unit mass, f b Resultant of the body forces per unit volume acting on the continuum, H Total enthalpy, I Moment of inertia, p Pressure, q Heat flux, S u and s u User-specified source term, T Viscous stress tensor, t Time, v Velocity, V Volume The fundamental laws which describe the electromagnetic behavior of a material are Maxwell's equations and the conservation of electric charge. Maxwell's equations are given by Eqs 4-7 while the continuity equation for the charge within a control volume is given by Eq. 8.
ρ e Electric charge density, B Magnetic flux density, D Electric flux density, E Electric field, H Magnetic field, J Electric current density The relationship between the electric current density (J) and the electric field (E) is described by the generalised Ohm's law and is given by Eq. (9). The electrical conductivity (σ) describes the behaviour of the material in response to electric currents and depends on the plasma pressure and temperature.
J σE Magnetohydrodynamics (MHD) describes the interaction between electrically conducting fluids (such as plasma) and electromagnetic fields. The plasma experiences a body force per unit volume known as the Lorentz force (f L ) and is described by Eq. 10. Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 748113

Computational Domain and Mesh
The computational domain is based on the geometry investigated by Ref. [25]. An annotated two-dimensional view of the computational domain is shown in Figure 1.
The mesh is refined in the inter-electrode gap since this is the physical space where the plasma arc is initially ignited, which shall result in high spatial gradients. The volume mesh was generated using polyhedral cells. Prism layers were generated next to wall surfaces and boundary conditions. The prism layers can capture near wall flow accurately, which is critical for determining the forces and heat transfer on walls.
A cell base size of 0.1 mm was used for the bulk Helium gas region. A volumetric control in the form of a cylinder was used in the inter-electrode gap to achieve a smaller cell base size in this region. The cylinder was created with a diameter of 2 mm and length of 1 mm. The volumetric control contained a cell size corresponding to 10% of the base size. The prism layers consisted of two layers with a total thickness corresponding to 10% of the cell base size. The volume mesh contained a total of 1763280 cells and can be viewed in Figure 2.

Simulation Setup and Boundary Conditions
In plasma arc simulations, the Ohmic heating represents a source term in the energy equation. However, modelling an electric current through an isolated gas at room temperature causes numerical instability due to the large scale of Ohmic heating and a low initial electrical conductivity. To overcome this obstacle, an artificial electrical conductivity (AEC) channel in the Helium region was defined. This channel allowed for the electrical conductivity to be artificially raised by specifying a minimum electrical conductivity which allowed for the ignition of an electrical spark. Similar approaches to ignite the arc have also been used by Refs. [25][26][27][28][29].
The AEC channel is defined as a cylinder with a diameter of 0.2 mm and length of 3 mm and is positioned in the interelectrode region. The AEC is defined using a field function and states that σ 1 S/m in the AEC channel.
The arc is further stabilized during the initial iterations by specifying an initial hot channel (IHC) as an initial condition for the Helium gas. The IHC is positioned in the inter-electrode region. A field function was used to define the initial temperature which states that T 16,000 K (σ 3686.39 S/m) in the IHC and 300 K (σ 7×10 -150 S/m) elsewhere in the computational domain.
These AEC and IHC conditions are specific to the p 0.1 MPa, I 0.35 A simulation. A variation in pressure and electric current requires an adjustment of the AEC and IHC conditions to ignite the plasma discharge.
The time step is set to 0.25 μs for the first 2 ms of the simulation. Thereafter the time step is set to 1 μs up to a physical time of 50 ms [25].
Net Emission Coefficients (NEC) are widely used in CFD simulations and are based on the assumption of Local Thermodynamic Equilibrium [30]. The use of NEC allows the simulation to pre-compute the energy balance within the arc among emitted, self-absorbed, and radiated energy. Therefore, it provides an estimation of the net radiation lost by the arc, by assuming a simple shape and currentdependent size for the arc. Refs. [31][32][33] have shown that the NEC is a good approximation to be used in numerical modelling in order to characterize the radiation losses in hot regions while taking into consideration the absorption. The data for the NEC values used in all the simulations are provided at atmospheric pressure [34].
The Lorentz force was specified as a user-defined momentum source. A field function was used to define the Lorentz force as specified in Eq. 10.
Reference [35] showed that an arc discharge can be split into three main domains: the near-cathode zone, the positive arc column and the near-anode zone. The plasma increasingly deviates from equilibrium conditions closer to the near-cathode and near-anode zones. The result is a strong difference in electrical behaviour between these zones. The discharge voltage (U) is expressed as the sum of the positive arc column voltage drop (U column ), the near-cathode (U near−cathode ) and the near-anode voltage drops (U near−anode ): An empirical formula proposed by Ref. [36] was used to define the near-cathode and near-anode voltage drop. The empirical formula is given in Eq. 12, where a, d are coefficients for a gas at atmospheric pressure.
All boundaries have been defined as walls. A wall boundary represents an impermeable surface which confines the solid (electrodes) and fluid (Helium plasma) regions. The cathode and anode boundaries were split as defined by AB and CD as shown in Figure 1. The boundary conditions for the simulation are given in Table 1.
A negative polarity was specified for the cathode and the electric current was imposed in the region of 0.25-0.45 A. At each time step (n), the cathode electric potential is fitted to maintain the target electric current.

RESULTS
A parametric study has been performed to determine the effect of pressure. This study was performed in the pressure range of 2-8 MPa with the electric current fixed at 0.35 A. Thereafter, the simulation results at atmospheric pressure (p 0.1 MPa) and at different electric currents (I 0.25-0.45 A) are presented. The model has been verified comparing numerical results with experimental measurements obtained at atmospheric pressure.

Influence of Pressure When Electric Current Is Fixed (I = 0.35 A)
The pressure was varied between 2 and 8 MPa to determine the influence of pressure on electric potential (ϕ), electric current density (J), temperature (T), and velocity magnitude (|v|). The electric current was fixed at 0.35 A. The results are extracted at the center of the inter-electrode gap in the x axis direction (see Figure 1).
The temperature and density variation in the x axis direction at the center of the inter-electrode gap is given in Figure 3 and Figure 4 respectively. There is a small deviation between maximum and minimum temperatures for all cases. The high temperature in the discharge region as opposed to the lower temperature at the walls of the reactor can be clearly observed.
This variation in temperature results in high density gradients as can be seen in Figure 4. The density variation at the cooler regions of the reactor (i.e., −6 and +6 mm distance) is much larger at 8 MPa when compared to 2 MPa. At +6 mm distance, the increase in density between 2 and 8 MPa is 311%.
The large density gradient in the x axis direction, combined with the gravitational force acting in the negative x direction as per Figure 1, becomes the driving force of the natural convection in the reactor at higher pressures. The result of the increased natural convection effect at higher pressure can be seen in the velocity magnitude variation ( Figure 5). A larger velocity magnitude variation is observed at higher pressures.
The electric current density variation in the x axis direction at the center of the inter-electrode gap is shown in Figure 6. A higher pressure causes a displacement of the arc core and anode arc root towards the top of the reactor. This is due to a higher pressure causing an increased natural convection effect. Figures 7, 8 compares the velocity magnitude and temperature profiles respectively, across the entire pressure range. In the 8 MPa case, the natural convection effect causes the Helium gas with the highest velocity magnitude to move towards the top of the reactor. The movement of the gas towards the top of the reactor becomes less pronounced as the pressure decreases. In the 0.1 MPa case, recirculation zones form alongside the cathode. The natural convection effect is greatly reduced, and the maximum velocity magnitude is 66% lower than the 8 MPa case.
The velocity profiles from the 8 and 6 MPa simulations considering a laminar flow model are shown alongside the turbulent results in Figure 7. This comparison was performed to evaluate the effect of turbulent and laminar viscous regimes on the Helium gas velocity field in the reactor. A similar velocity field in the reactor is observed for both the turbulent and laminar models at a physical time of 50 ms. The velocity magnitude calculated by the turbulent model is 14% lower than that from the laminar flow model for the 8 MPa case, while there are no appreciable differences for the 6 MPa case.
The natural convection effect causes a slight bending of the discharge as seen in the 8 MPa temperature profile. A decrease in pressure causes the bending effect to diminish, and a low current arc type of discharge is observed at lower pressures.
Atmospheric Pressure Discharge ( p = 0.1 MPa) Plasma discharges in contact with liquids are often investigated at atmospheric pressure (p 0.1 MPa) [10,[37][38][39]. It is for this reason that the simulated discharge will be further investigated at atmospheric pressure. The investigation has been performed at an electric current range of 0.25-0.45 A. A low current arc discharge can be expected under these operating conditions [40,41].  The Experimental Setup and the Evolution of the Electric Potential The reactor used in the experiment and the one simulated in this study have the following similarities: the dimensions are equivalent and therefore the cylindrical discharge chamber is approximately 2.5 cm 3 , the reactor is aligned in the horizontal direction, the cathode has a conical tip and is connected to the high voltage port of the power generator which is operated under negative polarity and DC voltage, the reactor operates with Helium, and the inter-electrode gap is 1 mm and the operating pressure is 0.1 MPa, while varying the electric current. The experimental electrical characterization was conducted in a  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 748113 6 tip-to-plane electrode configuration in a reactor chamber with capability of up to 20 MPa. The reactor is equipped with a K-type thermocouple to measure the temperature in the reactor and a water-cooling jacket. The reactor is powered by a Technix high voltage power supply, which is current controlled and set at an ignition voltage of 9 kV. The operating current is set on the power supply unit and can be measured using a Chauvin hall effect current clamp, while the voltage can be measured using an Elditest 30-kV high voltage probe, connected to the Wave Jet 3540 A digital oscilloscope via BNC connect cables from the current and voltage probes. The Helium gas is filled into the reactor from the Helium  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 748113 7 baseline 5.0 cylinder to achieve the operating pressure. The Helium gas cylinder was purchased from Afrox, with a minimum purity of 99.999%.
In the experiment, the two electrodes are initially in contact, and the discharge chamber was filled with Helium at atmospheric pressure. A voltmeter was also connected in parallel to verify the consistency of the recorded voltage. The same procedure was repeated for the range of electric currents investigated.
The experimental results are compared with the simulation results in Figure 9. Both the simulation results and experiment measurements show a general trend of the electric potential decreasing as the electric current increased. The experimental and simulation results are in closer agreement at lower electric currents.
Discharge Behavior at Fixed Conditions (p 0.1 MPa and I 0.35 A) The temperature, density variation and velocity magnitude in the x axis direction at the center of the inter-electrode gap for the atmospheric pressure discharge when I 0.35 A is given in Figures 10-12 respectively.
The density variation between maximum and minimum values for the atmospheric pressure discharge is not as large as that observed for higher pressures (see Figure 4). This observation indicates that the natural convection effect will not be as dominant as is the case with higher pressures. These plots for the atmospheric pressure discharge show a symmetric shape around the 0 mm point which is expected from a low current arc type discharge.

Discharge Behavior Under Varying Electric Current (p 0.1 MPa)
The temperature in the x axis direction at the center of the inter-electrode gap for the atmospheric pressure discharge for varying electric current is given in Figure 13. Panel (a) shows a general trend of increasing temperature with increasing electric current. Panel (b) shows that an almost identical discharge shape is obtained for the electric currents in the range considered.
The electric current density in the x axis direction at the center of the inter-electrode gap at atmospheric pressure discharge for  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 748113 8 different electric currents is given in Figure 14. The trend shown in panel (a) is similar to the temperature plot (see Figure 13) showing that the electric current density tends to increase with electric current. As the electric current increases, the ions bombarding the cathode also increase, inducing an increase in electric current density. Panel (b) shows that an almost identical electric current density profile is obtained across the electric current range which was investigated.
The velocity magnitude in the x axis direction at the center of the inter-electrode gap for the atmospheric pressure discharge for varying electric current is given in Figure 15.  The trend in panel (a) shows a very slight fluctuation in velocity magnitude as the electric current is varied. This observation shows that electric current does not have a large impact on velocity magnitude. Panel (b) again shows that an almost identical velocity magnitude profile is obtained across the electric current range which was investigated.

DISCUSSION
It is important to mention that given the operating conditions of the experiment (high pressure ≥ 0.1 MPa and low current ( < 1 A), we could safely state that there is a tendency for our discharge to be at the regime of "glow-to-arc transition" rather than a normal glow discharge or an abnormal glow discharge.  At a relatively low current (less than a microampere), the electric field strength is not sufficiently large enough to generate a self-sustained discharge. The glow discharge is categorized into two types: the normal glow, and the abnormal glow. As the voltage and current continue to increase, a glow-to-arc discharge transition begins to occur due to cathode ions bombarding which eventually leads to an arc discharge generated at current <1 A.
The main highlight in this research is the fact that according to Ref. [18], arc discharges are characterized by high current (I ∼ 1-10 5 A) while glow discharges by low current (I ∼ 10 -6 to  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 748113 11 10 -1 A). However, to the best of the knowledge of the authors, it is technically challenging to generate an arc discharge at a current < 0.1 A and pressure ≥ 1 MPa. Indeed, according to Ref. [18], high pressure refers to pressure, p ∼ 0.1-0.5 atm at which the plasma of the positive column is typically in equilibrium. The column of the atmospheric pressure arc is usually considered as a dense low-temperature equilibrium plasma with temperature, T ∼ 6000-12,000 K.
Values of pressure higher than, 1 MPa are considered as very high pressure, while for arc generated at p ∼10 -3 -10 0 Torr, the plasma in the positive column is strongly in nonequilibrium and does not differ in principle from the glow discharge. However, given that the pressure in our study is from 0.1 to 8 MPa, it can be referred to as both the high pressure and very high-pressure regime, and the plasma column could be assumed to be in a state of local thermal equilibrium. Moreover, our voltage-current curve ( Figure 9) follows that of a "glow-to-arc transition" from the current range of 0.25-0.45 A which is not previously reported in literature for the condition investigated. Furthermore, the large density gradient at the higher operating pressure causes a strong natural convection in the reactor, which drastically affected the discharge shape and the velocity field at higher pressures. A lower density gradient was observed at atmospheric pressure resulting in a reduced velocity magnitude.

CONCLUSION
The pressure was varied between 2 and 8 MPa to determine the influence of pressure on discharge parameters. The large density gradient in the x axis direction is present at higher pressure and has a strong influence on the velocity magnitude in the reactor.
The simulation and experimental results at atmospheric pressure showed the general trend of decreasing electric potential values as the current increased. Both sets of results were in closer agreement at the lower range of electric current. Any discrepancy between the experimental and simulated electric potential is mainly due to the electrode characterization, a phenomenon which is not fully understood within the context of low current and high pressure as supported by Ref. [42].
The natural convection effect is not dominant at atmospheric pressure. This can be observed in the symmetry of the temperature, density and velocity magnitude plots which were generated at an electric current of 0.35 A.
At atmospheric pressure, a larger temperature value is calculated as the electric current increases. An almost identical low current arc discharge shape was observed across the electric current range which was investigated. An increase in electric current resulted in the calculation of larger electric current density values. This was expected as more energy is being added to the reactor system.
Increasing electric current does not have a large impact on velocity magnitude at atmospheric pressure. A lower velocity magnitude was observed at atmospheric pressure compared to a high pressure. The velocity plot at atmospheric pressure is different to that at high pressure due to the less pronounced natural convection effect.
The study has shown that the CFD results can provide useful trends to design a high pressure, low current plasma reactor. This data will be used to influence the design of a reactor. Future works will investigate the interaction of the plasma discharge with a liquid at atmospheric pressure for applications of wastewater treatment. Future developments can include investigation on streamers and ionization waves [43]. Such a simulation can present the challenges of low pressures and super short timescales.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
AM was the main author and responsible for the first draft of the manuscript. All authors provided review and comment on the subsequent versions of the manuscript. AM was responsible for CFD simulations; experimental measurements; analysis and interpretation of results. AD and GC provided suggestions to improve the mathematical model and the manuscript. They also performed the calculation of thermodynamic and transport properties of Helium plasma at different temperature and pressures using the EquilTheTA tool. SI was the main academic and research supervisor. He has provided support with the experimental measurements, model setup, analysis, and interpretation of results. He also provided input, review and editing towards the manuscript. All authors read and approved the final manuscript.

FUNDING
The corresponding author, SI, is funded by the Royal Society as a FLAIR Fellow (FLR\R1\201683).