Effect of Cathode Cooling in Three-Dimensional Simulations of an Atmospheric Pressure Glow Discharge

Abstract

The Atmospheric Pressure Glow Discharge (APGD) is a relatively simple and versatile plasma source used in a wide range of applications. Active cooling of the cathode can effectively mitigate instabilities, leading to glow-to-arc transitions. This study investigates the effect of varying the degree of cathode cooling in APGD with a planar cathode in helium. The plasma flow model incorporates mass conservation, chemical species transport, momentum conservation, conservation of thermal energy of heavy species and of electrons, and electrostatics. The model is applied to time-dependent simulations through a three-dimensional computational domain describing the whole discharge, without geometric symmetry or steady-state assumptions. Simulations of an experimentally characterized APGD explore the effects of electric current and cathode cooling—ranging from thermally insulated to extreme convective cooling. Results show the formation of an annular region with high electric field over the cathode surface under conditions of high current and low cooling.

1. Introduction

Low-Temperature Plasma (LTP) sources are used in an increasing number of industries and applications ranging from materials processing and chemical synthesis to water treatment, agriculture, and medicine [1,2]. This interest is motivated by the ability of LTP sources to provide a highly reactive environment through the direct utilization of electric power. Among LTP sources, the Atmospheric Pressure Glow Discharge (APGD) is particularly relevant to a wide range of technological applications, such as lighting, spectroscopy, surface treatment, decontamination, and more recently, CO₂ conversion [3].

The APGD has been investigated using various working gases and geometric configurations. Stable APGDs have been demonstrated over spatial extents from tens of micrometers to tens of millimeters and for discharge currents ranging from micro-Ampere to milli-Ampere. The electric current range is generally dependent on the electrode set-up and operating conditions. Operation at high electric currents induces various instabilities, particularly thermal instabilities in the cathode fall region, which lead to the so-called glow-to-arc transition. An effective approach to avoid the development of fluctuations in the cathode fall region has been presented by Arkhipenko et al. [4] consisting of active cooling of the cathode, specifically by incorporating a planar cathode cooled by the impingement of a water jet. Their approach successfully extended the current range for stable cathodic APGD operation to over 1 A. In Figure 1, we schematically depict their experimental setup, where the plasma is generated between a plate cathode cooled by forced convection cooling and a slightly rounded pin anode. The figure also presents a snapshot of experimental results demonstrating a complex protruded circular pattern along the cathode. Different experimental studies of APGD have shown the formation of intricate patterns over the electrodes [5]. Pattern formation is a fascinating phenomenon with relevance to various applications. For example, pattern formation can be exploited for surface texturization but can be undesirable when it limits electrode life.

Computational modeling and simulation have become crucial tools for understanding, characterizing, and optimizing plasma sources and processes [6,7,8]. Computational simulations of APGD operation in a setup like the one presented in Figure 1 generally require a comprehensive mathematical model. Often, computational models are based on descriptions of the discharge relying on lower dimensionality (i.e., one- or two-dimensional, axisymmetric, and steady-state instead of three-dimensional and time-dependent) or limited physical fidelity (e.g., reduced set of chemical species and reactions, omission of thermal effects, simplified treatment of near-electrode regions, etc.). The large computational cost of 3D plasma flow simulations has been the primary limiting factor for their wider adoption. The recent work by Maerivoet et al. [9] on the study, an APGD for CO₂-CH₄ conversion, exemplifies the intricacies yet necessity of 3D simulations. Their model included gas flow, heat balance, and species transport, coupled with a comprehensive chemistry set, and was able to describe the major aspects of the chemical conversion process. Despite their relatively large computational cost compared to lower-dimensionality models, two- or three-dimensional models are generally required to unveil pattern formation phenomena [10,11]. Additionally, time-dependent models are essential to capture the occurrence of instabilities [12]. Our previous work, [13], addressed the modeling and simulation of the APGD in helium studied by Arkhipenko et al. [4]. That work revealed some important aspects of the discharge, particularly the predominance of thermal effects with increasing current; nevertheless, the formation of electrode patterns was not observed. In contrast, numerical simulations by Bienek and Hasan [14] of a micro-APGD in helium showed the occurrence of pattern formation over the cathode. Their findings demonstrated the importance of including thermal effects in captured pattern formation events. The state-of-the-art computational models for studying electrode pattern formation phenomena are based on stationary (steady-state) solvers, exemplarily depicted by the work of Benilov, Almeida, and collaborators [10,14,15]. In contrast, time-dependent (transient) solvers can generally capture a single pattern configuration [11]. Nevertheless, time-dependent solvers can describe the dynamic evolution of instabilities and can be better suited for the simulation of LTP applications.

This article presents a computational investigation of the effect of cathode cooling on the APGD studied by Arkhipenko et al. [4] depicted in Figure 1. The plasma flow model consistently describes chemical kinetics, thermal-fluid transport, and electrostatics. We considered two control parameters of the system: total discharge current and level of cathode cooling. The discharge current varies from 5 to 40 mA, well within the range experimentally investigated by Arkhipenko et al. [4]. Cathode cooling is modeled as convective cooling with a convective heat transfer coefficient ranging from 10⁻¹ to 10⁶ W·m⁻²·K⁻¹, describing from thermally insulated to extreme convective cooling. We conducted computational experiments using a complete 3D description of the spatial domain (i.e., no assumptions about solution symmetry) using a time-dependent numerical solution approach (i.e., no assumptions about steady-state). The model, together with the limited assumptions in the setup of simulations, is aimed at being a step towards comprehensive, predictive simulations capable of capturing the potential occurrence of instabilities and pattern formation events.

2. Materials and Methods

2.1. Mathematical Model

2.1.1. Transport Equations

The APGD is treated as a quasineural fluid in thermal and chemical nonequilibrium through the set of transport equations describing mass conservation, conservation of chemical species, momentum conservation, conservation of thermal energy of heavy species and of electrons, and electrostatics.

The Navier-Stokes equations, used to describe the fluid flow, are comprised of the equations of conservation total mass and of mass-averaged momentum, namely:

$${\partial }_t\rho +\mathit{u}\cdot \nabla \rho +\rho \nabla \cdot \mathit{u}=0, \mathrm{a}\mathrm{n}\mathrm{d}$$

(1)

$$\rho {\partial }_t\mathit{u}+\rho \mathit{u}\cdot \nabla \mathit{u}+\nabla p+\nabla \cdot \mathit{\tau }=\mathbf{0},$$

(2)

where the operators ${\partial }_t$ and $\nabla$ represent the partial derivative with respect to time and the gradient operator, respectively; ρ is mass density, u velocity, p pressure, and τ the stress tensor. The Lorentz force due to the self-induced magnetic field is assumed negligible, and hence it is not included in Equation (2). Instead of pressure p, the model uses the gauge pressure $p^{\prime }$ as unknown, where $p^{\prime }=p-p_{ref}$ and $p_{ref}$ is set to 1 atm.

Given the state of thermal nonequilibrium of the plasma, heavy species and free electrons have different energy distribution functions and, therefore, dissimilar (equivalent) temperatures. The conservation equations for the heavy-species thermal energy and electrons thermal energy are expressed as:

$$\rho {\partial }_t{h}_h+\rho \mathit{u}\cdot \nabla h_h+\nabla \cdot {\mathit{q}}_h-D_t{p}_h+{\dot{Q}}_{eh}=0, \mathrm{and}$$

(3)

$$\rho {\partial }_t{h}_e+\rho \mathit{u}\cdot \nabla h_e+\nabla \cdot {\mathit{q}}_e-D_t{p}_e-{\dot{Q}}_{eh}-{\dot{Q}}_{ch}+{\mathit{J}}_q\cdot \mathit{E}=0$$

(4)

where $h_h$ and $h_e$ are the heavy species and electrons specific enthalpy, respectively; ${\mathit{q}}_h$ and ${\mathit{q}}_e$ are the total heat fluxes for heavy species and electrons, respectively; $D_t\equiv {\partial }_t+\mathit{u}\cdot \nabla$ represents the material derivative; and p_h and p_e are the partial pressure of heavy-species and of electrons, respectively. The term ${\dot{Q}}_{eh}$ describes the energy exchange between electrons and heavy species due to elastic collisions, whereas the term ${\dot{Q}}_{ch}$ the energy exchange due to inelastic collisions. These source terms are given by:

$${\dot{Q}}_{eh}=3k_B{n}_e{n}_{He}{\displaystyle \frac{m_e{m}_{He}}{{\left(m_e+m_{He}\right)}^2}}{k}_{r,1}\left(T_e-T_h\right) \mathrm{and}$$

(5)

$${\dot{Q}}_{ch}={\sum }_r\left(\mathsf{\Delta }{H}_r k_r{\prod }_s{n}_s^{{\beta }_{\mathrm{s}\mathrm{r}}}\right),$$

(6)

where k_B is the Boltzmann constant, n denotes number density and m the particle mass of specie s, $k_{r,1}$ represents the rate coefficient of elastic collisions between helium atoms and electrons, $\mathsf{\Delta }{H}_r$ is the change in enthalpy and $k_r$ the rate coefficient for reaction r, respectively, and ${\beta }_{sr}$ is the stochiometric coefficient for species s in reaction r. The chemical kinetics model, which includes the specification of the set of reactions and their parameters, is described in Section 2.1.2. Electron Joule heating is modeled through the term ${\mathit{J}}_q\cdot \mathit{E}$, where ${\mathit{J}}_q$ is the current density and $\mathit{E}$ is the electric field, while the Joule heating due to ion motion is neglected in Equation (3). Moreover, the terms describing enthalpy transport due to mass diffusion have been neglected (i.e., $h_e{\mathit{J}}_e$ in Equation (3) and ${\sum }_{s\ne e}{h}_s{\mathit{J}}_s$ Equation (4)) to reduce the computational cost of simulations.

As the plasma is considered quasi-neutral, therefore there is no charge accumulation, and the ions’ total charge balances the electron number density. This leads to the electric current continuity equation describing the electric potential $\varphi$:

$$\nabla ·(\sigma \nabla \varphi )=0,$$

(7)

where $\sigma$ represents the local electrical conductivity. The electric field is determined from the electric potential as $\mathit{E}=-\nabla \varphi$ and the current density is given by ${\mathit{J}}_q=\sigma \mathit{E}$. The quasineutrality assumption is generally valid in the bulk of the discharge. However, it is generally invalid in the vicinity of the electrodes. Nevertheless, our previous work showed that a quasineutral formulation throughout the domain can reasonably well represent a glow discharge’s overall thermal and electrical characteristics [13].

The helium plasma is modeled as composed of six species: electrons ($e^{-}$), helium atoms (He), atomic helium in an excited state ($H{e}^{*}$), singly ionized helium ion ($H{e}^{+}$), rotationally excited bimolecular helium ($H{e}_2^{*}$), and singly ionized molecular ion ($H{e}_2^{+}$). The excited helium atom $H{e}^{*}$ represents helium at three different electronic excitation levels, namely $He\left(2^{3}S\right)$, $He\left(2^{3}P\right)$, and $He\left(2^{1}S\right)$. The species $H{e}_2^{*}$ represents a molecular helium molecule in the rotational state $H{e}_2\left(a^3{\mathsf{\Sigma }}_u\right)$. The evolution of the plasma composition is determined by a set of conservation equations for species mass fraction of the form:

$$\rho {\partial }_t{y}_s+\rho \mathit{u}\cdot \nabla y_s+\nabla \cdot {\mathit{J}}_s-{\dot{R}}_s=0,$$

(8)

where $y_s$ represents the mass fraction associated with species $s$. The diffusive mass flux is denoted with ${\mathit{J}}_s$, and the source term ${\dot{R}}_s$ describes chemical reactions.

The species mass diffusion flux is given by:

$${\mathit{J}}_s=-\rho D_s\nabla y_s+{\displaystyle \frac{q_s}{\left|q_s\right|}}{\mu }_s{n}_s{\mathit{E}}_{amb},$$

(9)

where $D_s$ is the diffusion coefficient for species $s$, $q_s$ is the charge of species $s$ (equal to zero for neutral species), ${\mu }_s$ is the mobility coefficient of species $s$, and $n_s$ the number density of species $s$. The ambipolar electric field ${\mathit{E}}_{amb}$ is determined by:

$${\mathit{E}}_{amb}={\displaystyle \frac{\left(D_{{He}^{+}}-D_e\right)\nabla n_{{He}^{+}}+\left(D_{{He}_2^{+}}-D_e\right)\nabla n_{{He}_2^{+}}}{\left({\mu }_{{He}^{+}}+{\mu }_e\right)n_{{He}^{+}}+\left({\mu }_{{He}_2^{+}}+{\mu }_e\right)n_{{He}_2^{+}}}}.$$

(10)

The mass diffusion flux in Equation (9) does not account for multicomponent diffusion and is only valid for low degrees of ionization, which is a suitable approximation for APGDs. The ambipolar diffusion model in Equation (10) is similar to that adopted in Sun et al. [16] and in Zhang et al. [17], where it proved to be reasonably effective to represent the dynamics among equally charged ions of different masses [18].

Given the imposition of total mass conservation and quasineutrality, only four species conservation equations of the form of Equation (8) are solved in the model, namely, those for species $H{e}^{*}$, $H{e}^{+}$, $H{e}_2^{*}$, and $H{e}_2^{+}$ (i.e., s = {e⁻, He, $H{e}^{*}$, $H{e}^{+}$, $H{e}_2^{*}$, $H{e}_2^{+}$}). Moreover, while Equation (8) is written in terms of the species mass fraction $y_s$, our model adopts the logarithm of the mass fraction

$$z_S=\ln \left(y_s\right),$$

(11)

as the variable describing the population of species s. This formulation provides two numerical advantages: (i) $z_s$ is bounded to zero as its minimum value and (ii) while $y_s$ can vary by many orders of magnitude among different species and within a simulation, its logarithm $z_s$ is approximately of the same order of magnitude for every species.

Finally, the above species transport equations are complemented by the expressions describing the conservation of total mass and the quasineutrality condition, from which the density of helium atoms $n_{He}$ and of electrons $n_e$ are determined, respectively. These expressions are given by:

$${m_{He}n}_{He}=\rho -{\sum }_{s\ne He}{m_{s}n}_s \mathrm{and}$$

(12)

$$n_e=n_{He+}+n_{He2+}$$

(13)

respectively, where m_s is the mass of species s.

2.1.2. Chemical Kinetics

Our helium APGD model comprises 23 chemical reactions among the six species considered; these reactions are listed in

Table 1. Reactions for the APGD in the helium model. In the rate coefficient expressions, $k_r$ has units m³·s⁻¹ for binary reactions and m⁶·s⁻¹ for ternary reactions; $f(E/N)$ indicates that the rate coefficient is fitted from results from BOLSIG+ as a function of the magnitude of the reduced electric field $E/N$ assuming a Maxwellian velocity distribution; and $T_e$ and $T_h$ are Kelvin [19,23,24,25,26,27].

	Reaction			Rate Coefficient ${\mathit{k}}_{\mathit{r}}$	Energy $\mathsf{\Delta }{\mathit{H}}_{\mathit{r}}$
#				(m3·s−1, m6·s−1)	(eV)	Ref.
1	$e^{-}+He$	$\to$	$e^{-}+He$	$f(E/N)$	0	[19]
2	$e^{-}+He$	$\to$	$e^{-}+H{e}^{*}$	$f(E/N)$	19.82	[19]
3	$e^{-}+He$	$\to$	$2e^{-}+H{e}^{+}$	$f(E/N)$	24.58	[19]
4	$e^{-}+H{e}^{*}$	$\to$	$2e^{-}+H{e}^{+}$	$7.021\cdot {10}^{-14}{T}_e^{-0.1241}{e}^{-66400/T_e}$	4.78	[24]
5	$e^{-}+H{e}^{*}$	$\to$	$e^{-}+He$	$2.90\cdot {10}^{-15}$	−19.82	[25]
6	$e^{-}+H{e}_2^{*}$	$\to$	$e^{-}+2He$	$3.80\cdot {10}^{-15}$	−17.9	[25]
7	$2e^{-}+H{e}^{+}$	$\to$	$e^{-}+H{e}^{*}$	$7.80\cdot {10}^{-50}{\left(T_e/T_h\right)}^{-4.4}$	−4.78	[24]
8	$2e^{-}+H{e}_2^{+}$	$\to$	$e^{-}+He+H{e}^{*}$	$2.80\cdot {10}^{-32}$	0.0	[25]
9	$e^{-}+He+H{e}_2^{+}$	$\to$	$2He+H{e}^{*}$	$3.50\cdot {10}^{-39}$	0.0	[25]
10	$2e^{-}+H{e}_2^{+}$	$\to$	$e^{-}+H{e}_2^{*}$	$1.20\cdot {10}^{-33}$	0.0	[25]
11	$e^{-}+He+H{e}_2^{+}$	$\to$	$He+H{e}_2^{*}$	$1.50\cdot {10}^{-39}$	0.0	[25]
12	$2H{e}^{*}$	$\to$	$e^{-}+H{e}_2^{+}$	$2.03\cdot {10}^{-15}({T_h/300)}^{0.5}$	−18.2	[23]
13	$2H{e}^{*}$	$\to$	$e^{-}+He+H{e}^{+}$	$8.70\cdot {10}^{-16}({T_h/300)}^{0.5}$	−15.8	[23]
14	$2He+H{e}^{+}$	$\to$	$He+H{e}_2^{+}$	$1.4\cdot {10}^{-43}({T_h/300)}^{-0.6}$	0.0	[23]
15	$2He+H{e}^{*}$	$\to$	$He+H{e}_2^{*}$	$2.0\cdot {10}^{-46}$	0.0	[23]
16	$H{e}^{*}+H{e}_2^{*}$	$\to$	$e^{-}+2He+H{e}^{+}$	$5.0\cdot {{10}^{-16}\left(T_h/300\right)}^{0.5}$	−13.5	[25]
17	$H{e}^{*}+H{e}_2^{*}$	$\to$	$e^{-}+He+H{e}_2^{+}$	$2.0\cdot {10}^{-15}{\left(T_h/300\right)}^{0.5}$	−15.9	[25]
18	$2H{e}_2^{*}$	$\to$	$e^{-}+3He+H{e}^{+}$	$3.0\cdot {10}^{-16}({T_h/300)}^{0.5}$	−11.3	[25]
19	$2H{e}_2^{*}$	$\to$	$e^{-}+2He+H{e}_2^{+}$	$1.2\cdot {10}^{-15}({T_h/300)}^{0.5}$	−13.7	[25]
20	$e^{-}+H{e}^{+}$	$\to$	$H{e}^{*}$	$6.28\cdot {10}^{-21} T_e^{-0.5}$	0.0	[26]
21	$e^{-}+He+H{e}^{+}$	$\to$	$He+H{e}^{*}$	$7.4\cdot {10}^{-47}{\left(T_e/T_h\right)}^{-2}$	0.0	[24]
22	$e^{-}+H{e}_2^{+}$	$\to$	$He+H{e}^{*}$	$7.12\cdot {10}^{-21}{\left(T_e/T_h\right)}^{-1.5}$	0.0	[27]
23	$e^{-}+H{e}_2^{*}$	$\to$	$2e^{-}+H{e}_2^{+}$	$7.50\cdot {10}^{-13}{T}_e^{0.71}{e}^{-43400/T_e}$	3.4	[23]

. This set of species and reactions has been demonstrated to provide an appropriate description of helium glow discharges [13]. The first three reactions represent elastic collisions between electrons and helium atoms, excitation of ground-state helium, and electron-impact ionization of ground-state helium. Their rate coefficients are given by fitting the results obtained via BOLSIG+ [19] as a function of the reduced electric field $E/N$, where $E=\left|\right|\mathit{E}\left|\right|$ is the magnitude of the electric field and $N={\sum }_s{n}_s$ the total number of particles. The collision cross-section data was gathered from the Morgan database [20]. The rate coefficients for the other twenty reactions were compiled and cross-checked from multiple sources are expressed as a function of the heavy-species temperature $T_h$ and the electron temperature $T_e$ [21,22].

2.1.3. Near-Electrode Region Model

In this study, the mathematical plasma flow model adopts the assumption of quasi-neutrality throughout the spatial domain. APGDs, particularly in configurations with larger gaps akin to the one under investigation, generally exhibit quasineutral behavior throughout most of the discharge, especially throughout the positive column, but also significant deviations from quasineutrality near the electrodes. The adoption of quasineutrality throughout the domain represents a trade-off between the computational cost of simulations (total compute time, including required mesh resolution and convergence of the numerical solution approach) and the degree of physical fidelity of the model (extent of phenomena described). Despite this limitation, our APGD model has demonstrated to provide reasonable agreement with experimental results [13], providing valuable insights and facilitating the analysis of atmospheric pressure glow discharge phenomena.

In order to account for the effect of dissimilar charge transport near electrode regions, we imposed the following mass fluxes of ions over the electrodes:

$${\Gamma }_{s,anode}={\displaystyle \frac{1}{4}}\rho y_s\sqrt{{\displaystyle \frac{8k_B}{\pi m_s}}},\mathrm{and}$$

(14)

$${\Gamma }_{s,cath}={\displaystyle \frac{1}{4}}\rho y_s\sqrt{{\displaystyle \frac{8k_B}{\pi m_s}}}+{\mu }_s{y}_s{E}_{cath}$$

(15)

where s = {He⁺, He₂⁺}, k_B is Boltzmann constant, and the term $E_{cath}$ depends on the characteristics of the cathode sheath region, and its value is set equal to 300 kV·m⁻¹. Moreover, to describe the near-cathode region, we employ a postprocessing approach to account for the voltage drop in the cathode sheath layer, fixed to 145 V based on experimental measurements [4,28,29]. Such an approach proved to be effective in reproducing experimentally observed APGD characteristics, as reported in [13].

2.1.4. Thermodynamic and Transport Properties

Closing the specification of transport equations in Section 2.1.1 requires the definition of thermodynamic and transport properties. The former relates intrinsic characteristics of the state of the plasma, i.e., enthalpy, pressure, number density, etc. The thermodynamic properties of the plasma are computed using standard results from kinetic theory of gases, e.g., ${\rho h}_e={\displaystyle \frac{5}{2}}{k}_B{n}_e{T}_e$, ${\rho h}_h={\sum }_{s\ne e}{\displaystyle \frac{5}{2}}{k}_B{n}_s{T}_h$, $p_e=k_B{n}_e{T}_e$, and $p_h={\sum }_{s\ne e}{k}_B{n}_s{T}_h$.

Transport properties are required to define diffusive fluxes τ, q_h, q_e, and J_s. The definitions of τ, q_h, and q_e can be considered standard in plasma flow models and are detailed in [13], whereas the specification of J_s used in the present work is given by Equation (9). Transport properties are obtained following two methodologies. A first-order Chapman–Enskog method [30] is used to calculate the species viscosity ${\mu }_s$ (used to compute the mass-averaged viscosity $\mu$, together with a mixing rule), and the thermal conductivities ${\kappa }_h$ and ${\kappa }_e$ associated with heavy species and electrons, respectively, based on a formulation presented in [31] using cross-sections from data in [32]. The diffusion coefficients and electrical mobilities are calculated from empirical expressions often used to simulate glow discharges [23,25]. The expressions for $\mu$_s and D_s used in the model are listed in

Table 2. Electrical mobility and diffusion coefficients. $T_h$ and $T_e$ are in units of eV, $p$ in Torr, and $N$ in m⁻³.

Species	Electrical Mobility ${\mathit{\mu }}_{\mathit{s}}$	Diffusion Coefficient ${\mathit{D}}_{\mathit{s}}$
	(m2·s−1·V−1)	(m2·s−1)	Ref.
$e^{-}$	$2.83\cdot {10}^{24}/N$	$(2.83\cdot {10}^{24}/N)T_e$	[23]
$H{e}^{*}$	-	$420p^{-1}{\left(T_h/0.025\right)}^{1.5}$	[25]
$H{e}^{+}$	$3.25\cdot {10}^{22}/N$	$3.25\cdot {10}^{22}{(T}_h/N)$	[23]
$H{e}_2^{*}$	-	$305p^{-1}{\left(T_h/0.025\right)}^{1.5}$	[25]
$H{e}_2^{+}$	$4.88\cdot {10}^{22}/N$	$4.88\cdot {10}^{22}(T_h/N)$	[23]

. The electrical conductivity $\sigma$ is computed based on the plasma composition and the charged species’ electrical mobilities. Detailed information about the definition of the thermodynamic and transport properties is found in [13].

2.1.5. Model Summary

The set of equations composing the APGD flow model is summarized in

Table 3. Plasma fluid model for APGD. Set of model equations expressed as a system of transport equations of the form Transient + Advective – Diffusive − Reactive = 0.

Conservation Equation	Transient	Advective	Diffusive	Reactive
Total mass	${\partial }_t\rho$	$\mathit{u}\cdot \nabla \rho +\rho \nabla \cdot \mathit{u}$	$0$	$0$
Momemtum	$\rho {\partial }_t\mathit{u}$	$\rho \mathit{u}\cdot \nabla \mathbf{u}+\nabla p$	$-\nabla \cdot \mathit{\tau }$	$0$
Heavy-species energy	$\rho {\partial }_t{h}_h$	$\rho \mathit{u}\cdot \nabla h_h$	$-\nabla \cdot {\mathit{q}}_h$	$D_t{p}_h+{\dot{Q}}_{eh}$
Electron energy	$\rho {\partial }_t{h}_e$	$\rho \mathit{u}\cdot \nabla h_e$	$-\nabla \cdot {\mathit{q}}_e$	$D_t{p}_e-{\dot{Q}}_{eh}+{\dot{Q}}_{ch}+{\displaystyle \frac{5k_B}{2e}}{\mathit{J}}_q\cdot \nabla T_e+{\mathit{J}}_q\cdot \mathit{E}$
Charge	$0$	$0$	$-\nabla \cdot {\mathit{J}}_q$	$0$
Species He*	${\rho }_{He*}{\partial }_t{z}_{He*}$	${\rho }_{He*}\mathit{u}\cdot \nabla z_{He*}$	$-\nabla \cdot {\mathit{J}}_{He*}$	${\dot{R}}_{He*}$
Species He+	${\rho }_{He+}{\partial }_t{z}_{He+}$	${\rho }_{He+}\mathit{u}\cdot \nabla z_{He+}$	$-\nabla \cdot {\mathit{J}}_{He+}$	${\dot{R}}_{He+}$
Species He2*	${\rho }_{He2*}{\partial }_t{z}_{He2*}$	${\rho }_{He2*}\mathit{u}\cdot \nabla z_{He2*}$	$-\nabla \cdot {\mathit{J}}_{He2*}$	${\dot{R}}_{He2*}$
Species He2+	${\rho }_{He2+}{\partial }_t{z}_{He2+}$	${\rho }_{He2+}\mathit{u}\cdot \nabla z_{He2+}$	$-\nabla \cdot {\mathit{J}}_{He2+}$	${\dot{R}}_{He2+}$

. The model describes the temporal evolution throughout the 3D domain of a total of 11 independent variables, namely {p′, u_x, u_y, u_z, T_h, T_e, ϕ, z_He*, z_He+, z_He2*, z_He2+}. The set of equations in Table 3 is expressed as a system of transient-advective-diffusive-reactive (TADR) transport equations of the form Transient + Advective – Diffusive − Reactive = 0, where each term represents a specific transport mechanism.

2.2. Computational Model

2.2.1. Solution Approach

The equations summarized in Table 3 are solved using our in-house TransPORT equation solver, TPORT. The discretization approach in TPORT is based on the Finite Element Method (FEM) using a nonlinear implementation of the Variational Multi-Scale (VMS) paradigm for the solution of transport problems [12]. The system of eleven transport equations is solved in a fully coupled, monolithic manner. This approach has shown its capacity to accurately solve a wide range of transport problems (incompressible, compressible, magnetohydrodynamic flows, radiative transport), including complex plasma flow problems [5,10].

State-of-the-art stationary solvers have demonstrated faster convergence and the ability to capture multiple solutions, as required to characterize pattern formation phenomena [14,15]. We use a second-order-accurate temporal solver, which generally requires longer computing time to achieve convergence to a steady-state solution (if it exists). However, this approach could be used to study the transient behavior of the discharge, including a potentially oscillatory steady-state. In this study, we implemented a segregated computational cycle to bridge the computational performance gap between stationary and time-dependent solvers. The computational cycle is composed of four stages, themselves composed of multiple timesteps, as follows:

• Stage 1: Two time steps during which only the total mass and momentum conservation equations are updated.
• Stage 2: Two time steps during which the heavy-species energy conservation equation is updated.
• Stage 3: Two time steps during which the electron energy conservation equation is updated.
• Stage 4: Four time steps during which the mass fraction conservation equations are updated alongside the charge conservation equation.

The timesteps of each of the stages are adjusted to achieve faster convergence toward steady-state. For example, the time step for Stage 1 can be multiple orders of magnitude larger than the time step used during Stage 4. This segregated solution approach was implemented by directly modifying the Jacobian matrix and the global residual vector resulting from the discretization of the system of equations in Table 3. In addition to achieving larger timesteps, this computational cycle provided better conditioned linear systems, requiring less iteration in our linear solver (scaled and flexible implementation of the Generalized Minimal RE Sidual (GMRES) method). Overall, we estimated that the convergence towards steady-state was fifty to a hundred times faster than using a single time-step (i.e., tightly coupled, monolithic) approach.

2.2.2. Computational Domain

The computational domain for the APGD corresponds to the experimental setup used by Arkhipenko et al. [4]. The rector geometry and the computational mesh used for its discretization are depicted in Figure 2. Figure 2a shows the geometry of the domain. The APGD reactor has a pin-to-plate configuration with a plate cathode undergoing controlled water cooling. The inter-electrode gap is h = 10 mm, and the total domain height is H = 15 mm. The circular plate describing the cathode has a diameter D = 24 mm, whereas the anode has a diameter d = 6 mm. The physical domain, as shown in Figure 2a, encompasses the complete cylindrical section (i.e., 360°), in contrast to the circular quadrant (i.e., 90° section) used in our previous work [13]. The latter implicitly assumes that the solution is symmetric along the axes defining the quadrant. Such an assumption is reasonable given the symmetry of the domain. Nevertheless, the model could depict solutions that are not symmetric (and can even be dynamic), as demonstrated in our previous work on the simulation of arc discharges [10].

Figure 2a also depicts the discretization of the spatial domain. The discretization grid comprises approximately 450 k hexahedron elements and over 435 k nodes and is conformal to the shape of the domain boundaries. Figure 2b represents a cut-view of the spatial domain passing through the main vertical axis together with the definition of the different domain boundaries.

2.2.3. Boundary and Initial Conditions

The boundary conditions used for the investigation of cathode pattern formation are listed in

Table 4. Boundary conditions for the simulation of a self-sustained APGD in helium. The simulations consider total current I_tot and cooling of the convective heat transfer coefficient over the cathode h_c as the controlling parameters.

Variable	Anode	Inlet	Outlet	Cathode
$p^{\prime }$	${\partial }_n{p}^{\prime }=0$	${\partial }_n{p}^{\prime }=0$	$p^{\prime }=0$	${\partial }_n{p}^{\prime }=0$
$\mathit{u}$	$\mathit{u}=\mathbf{0}$	$\mathit{u}={\mathit{u}}_{inlet}$	${\partial }_n\mathit{u}=\mathbf{0}$	$\mathit{u}=\mathbf{0}$
$T_h$	$-{\kappa }_h{\partial }_n{T}_h={-h}_a(T_h-T_{wa})$	$T_h=T_{inlet}$	${\partial }_n{T}_h=0$	${-{\kappa }_h\partial }_n{T}_h=-h_c\left(T_h-T_{wc}\right)$
$T_e$	${\partial }_n{T}_e=0$	${\partial }_n{T}_e=0$	${\partial }_n{T}_e=0$	${\partial }_n{T}_e=0$
$\varphi$	$\varphi =0$	${\partial }_n\varphi =0$	${\partial }_n\varphi =0$	$\varphi ={\varphi }_{cath} |{\int }_{cath}\sigma E_n=I_{tot}$
$z_s$	${\partial }_n{z}_s={\Gamma }_{s,anode}$	${\partial }_n{z}_s=0$	${\partial }_n{z}_s=0$	${\partial }_n{z}_s=$${\Gamma }_{s,cath}$

. In Table 4, n indicates the direction along the outer normal vector n to the given boundary and ${\partial }_n$ the derivative along that direction, i.e., ${\partial }_n=\mathit{n}·\nabla$. The inlet velocity profile is given by ${\mathit{u}}_{inlet}\left(r\right)={\mathit{u}}_{max}{\left(r/r_{u1}\right)}^{nu}\exp \left(-\left(r-d/2\right)/r_{u2}\right)$, where r is the radial coordinate from the main axis of the domain; nu, $r_{u1}$ and $r_{u2}$ are set equal to 4, 1.9 mm and 0.8 mm, respectively; and ${\mathit{u}}_{max}$ is the maximum velocity along the inlet equal to 0.164 $\mathrm{m}·{\mathrm{s}}^{-1}$. The velocity profile ${\mathit{u}}_{inlet}$ describes an induced flow rate of ~1 lpm. The inlet flow temperature is T_inlet = 500 K.

The anode and cathode surfaces experience convective heat transfer with reference temperature T_wa = T_wc = 500 K, respectively. The anode’s convective heat transfer coefficient is $h_a$= 10 W·m⁻²·K⁻¹, representative of free convection in ambient air, whereas the convective heat transfer coefficient for the cathode h_c is treated as a control parameter of the model with values ranging from thermal insulation to extreme convective cooling.

The total discharge current I_tot is imposed using a nonlinear boundary condition equivalent to an external electrical circuit. This boundary condition imposes a specified value of electric potential over the cathode ${\varphi }_{cath}$ such that the induced electric field over the cathode surface E_n, together with the distribution of electrical conductivity $\sigma$, lead to the satisfaction of continuity of total current, i.e., ${\int }_{cath}\sigma E_n=I_{tot}$.

Regarding species transport, for the ion species s = {$H{e}^{+}$, $H{e}_2^{+}$}, ${\Gamma }_{s,anode}$ and ${\Gamma }_{s,cath}$ are given by Equation (14) and Equation (15), respectively, whereas for the neutral species s = {$H{e}^{*}$, $H{e}_2^{*}$}, ${\Gamma }_{s,anode}={\Gamma }_{s,cath}=0$. A zero-gradient condition is imposed for all species over the inlet and outlet boundaries.

The computational experiments for the study are driven by two control parameters: the total discharge current $I_{tot}$ and the convective heat transfer coefficient for the cathode $h_c$. The study encompasses three levels of total current I_tot = {5, 20, 40} mA and four values of the convective heat transfer coefficient of the cathode h_c = {10⁻¹, 10¹, 10³, 10⁶} W·m⁻²·K⁻¹. The current levels were selected to be in the same range as those used in [13], allowing comparison with a different computational setup and with the experimental results from Arkhipenko et al. [4]. The cooling levels represent a well-insulated cathode (h_c = 10⁻¹), free convection cooling (h_c = 10¹), strong forced convection (h_c = 10³), and extreme convective cooling (h_c = 10⁶).

Given their time-dependent nature, the simulations require an appropriate set of initial conditions. In practice, the choice of initial conditions can have a significant effect on the simulations, from the total computing time to achieve a steady-state solution to the divergence of the numerical solution approach. As initial conditions, we used the steady-state results as a function of current from our previous study [13], interpolated to the new mesh (finer and spanning a complete cylinder) using Delaunay triangulation. This approach provided an initial condition well suited for the present numerical study, accelerating the convergence and robustness of the simulations.

3. Results and Discussion

3.1. Representative Characteristics of the APGD

Representative results of the helium APGD operating with a total discharge current I_tot = 20 mA and strong convective cooling of the cathode (h_c = 1000 W·m⁻²·K⁻¹) are shown in Figure 3. The results in Figure 3 depict the structure of the positive column in terms of the heavy-species temperature T_h, electron temperature T_e, and electron number density n_e. The distribution of T_h presents a globular shape with a maximum temperature of over a thousand Kelvin in the center of the discharge. The distribution of T_e shows a different structure from that of T_h; its shape is closer to a column, presenting a maximum of slightly over 3 eV, a value representative of APGDs. The distribution of n_e shows an hourglass shape, somewhat resembling the distribution of T_e, but depicting significant constriction near the center of the positive column.

Figure 4 presents the solution fields along a full radial slice across the center of the spatial domain. The figure’s first row shows the main fluid characteristic of the plasma, namely the distribution of gauge pressure $p^{\prime }$ as a contour plot, the fluid velocity $\mathit{u}$ as an overlay of streamlines and vector fields colored and sized following the local flow speed, and the distribution of heavy-species temperature T_h depicting the structure of the positive column. The pressure, velocity, and heavy-species temperature fields are similar to the results reported in [13].

The second row of Figure 4 presents the electrical characteristics of the plasma, namely the distribution of electric potential $\varphi$, the magnitude of the current density $\left|\right|{\mathit{J}}_q\left|\right|$, and electron temperature $T_e$. The voltage drop represented here is postprocessed to account for the near-cathode voltage drop of 145 V, as experimentally determined by Arkhipenko et al. [4]. The current density magnitude is constricted towards the center of the positive column, with a maximum of 8650 A·m⁻². The distribution of $T_e$ is significantly different from that obtained in [13]. While the electron temperature range is similar, with a maximum of ~3 eV, the spatial distribution is significantly different. The center of the discharge column still presents the highest temperature; however, there is significant heating outside the discharge center, along both the anode and the cathode. The difference appears local to the anode and does not propagate through the discharge.

The final row of Figure 4 presents the distribution of representative chemical species through the APGD, namely the electron number density $n_e$, the helium ion number density $n_{He+}$, and the molecular helium ion $n_{He2+}$. Despite the use of the ambipolar diffusion model for ${\mathit{J}}_{He+}$ and ${\mathit{J}}_{He2+}$, together with a near electrode model (Section 2.1.3), the distributions of $n_e$, $n_{He+}$, and $n_{He2+}$ appear relatively unchanged from those reported in [13].

3.2. Effect of Current and Cathode Cooling on APGD Characteristics

The dependence of the APGD structure with total current I_tot and convective heat transfer coefficient for the cathode h_c is depicted in Figure 5 by the distribution of electron number density n_e and heavy-species temperature T_h across a radial slice through the center of the domain. The results in Figure 5 show that the discharge experiences significant warming with an increasing current. The maximum heavy-species temperature for I_tot = 40 mA is almost 1000 K higher than the maximum temperature for I_tot = 5 mA. The cathode cooling does not appear to drastically influence the maximum temperature of the discharge. However, the level of cathode cooling strongly affects the heavy-species temperature field near the cathode. For the lower cooling cases (i.e., $h_c\le 10$ W·m⁻²·K⁻¹), the temperature along the cathode is close to the temperature near the center of the discharge. For the higher current cases, the temperature field displays a high gradient towards the cathode, and the temperature near the cathode is close to the reference convective cooling temperature of 500 K. The impact of the cathode cooling on the discharge temperature appears to be localized and does not impact greatly the value of temperature in the center of the discharge.

The electron number density displays a slightly different dependence on the control parameters, I_tot and h_c. As the discharge current increases, the electron temperature increases, and the plasma becomes more ionized. The maximum electron number density increases by over three orders of magnitude, as I_tot is increased from 5 to 40 mA.

Figure 6 compares the discharge’s global characteristics against the experimental results from Arkhipenko et al. [4]. Figure 6a presents the current–voltage characteristic curve, Figure 6b presents the positive column’s temperature, and Figure 6c shows current density in the center of the discharge as a function of total current I_tot and convective heat transfer coefficient for the cathode h_c. The obtained current–voltage characteristic curve agrees with the results reported in [4], especially for higher values of current. The discrepancy between simulation and reference results could be explained by the insufficient physical fidelity of the model to describe near-electrode phenomena. The results in Figure 6b show that the model overestimates the temperature of the plasma while capturing the experimentally observed trend of increasing temperature with respect to current I_tot. This discrepancy can be explained by two main factors. First, the temperature is determined experimentally [4] via optical spectroscopy from optical emission through the positive column, whereas the temperature shown in Figure 6b corresponds to that at the center of the discharge, at a location 5 mm from the cathode. Second, given the adiabatic boundary condition (i.e., zero temperature gradient) over the outlet boundary (Table 4), heat losses from the plasma only occur through the anode and cathode surfaces. The heat flux through those surfaces might be slightly underestimated. Moreover, given that the electrical mobilities presented in Table 3 are inversely proportional to the total number density of the plasma N, if the chemical composition of the plasma remains constant, its electrical conductivity should increase as the temperature increases. Thus, the discrepancy in voltage drop presented in Figure 6a could be partially explained by the temperature overestimation presented in Figure 6b. Finally, Figure 6c compares the current density in the center of the discharge to the measurement from Arkhipenko et al. [4]. Overall, the simulation results are in good agreement with the experimental measurements.

3.3. Near-Cathode Characteristics

Complementary to the results in Section 3.2, here we present the characteristics near the cathode of the APGD as a function of total current I_tot and convective heat transfer coefficient for the cathode h_c. Figure 7a presents a similar arrangement of results as in Figure 5, but along a slice at 0.01 mm parallel from the cathode of the electron number density n_e and heavy-species temperature T_h. For clarity, the range of T_h is different for each value of current. At the same time, the electron number density is presented in a logarithmic scale with a single range for all values of current and cathode cooling level. As noted in Figure 5, the level of cathode cooling greatly influences the temperature distribution along the cathode. At the highest cooling level (i.e., $h_c$ = 10⁶ W·m⁻²·K⁻¹), the temperature appears to be uniform along the cathode, as expected from an extreme cooling condition. The cases with $h_c$ = 10³ W·m⁻²·K⁻¹ show a near-uniform temperature along the cathode. In contrast to the distribution of T_h, those of the electron number density n_e do not demonstrate a strong dependency on the level of cooling of the cathode. The relative insensitivity of the electron number density to cathode cooling levels, in contrast to the heavy-species temperature, suggests that the mechanisms controlling the distribution of electron density are less directly influenced by the surface cooling conditions.

The results in Figure 7b depict the discharge structure near the cathode for I_tot = 20 mA and two levels of cathode cooling, namely free convection cooling and strong forced convection (i.e., h_c = 10 and 1000 W·m⁻²·K⁻¹, respectively). These results show the drastic change in the structure of the distribution of T_h with increasing cathode cooling level, whereas the distribution of n_e remains nearly unchanged. As described next, this has important implications for the formation of patterns over the cathode.

Figure 8 presents contour plots of the magnitude of the electric field E = ||E|| as a function of total current I_tot and cathode convective heat transfer coefficient h_c. The distribution of the magnitude of the electric field in Figure 8a displays a clear ring-like pattern along the surface of the cathode for the higher current (I_tot $\ge$ 20 mA) and lower cooling cases (h_c $\le$ 10 W·m⁻²·K⁻¹). This behavior was not observed in the APGD simulations reported in [13], which did not consider the full geometry or diffusion due to the ambipolar electric field. The patterns appeared spontaneously with a marked dependency on the level of cathode cooling and the magnitude of total current. The 5-mA discharges do not display pattern formation independently of the cooling of the cathode. The 20 and 40 mA cases display pattern formation for the four cases with $h_c$ < 10² W·m⁻²·K⁻¹. The structure of a ring pattern is more clearly depicted in the results in Figure 8b for I_tot = 20 mA and for free convection cooling and strong forced convection (i.e., h_c = 10 and 1000 W·m⁻²·K⁻¹, respectively), corresponding to the distributions of n_e and T_h shown in Figure 7b. The ring pattern is localized to the near-cathode region, depicting a toroidal shape centered at the axis of the discharge. This localized structure of the electric field is in contrast to the columnar structure of the electron number density observed in Figure 7b.

The emergence of ring-like patterns in the electric field at higher currents and lower cooling levels can be attributed to the increased electrical energy transport together with reduced heat dissipation. The occurrence of ring patterns or any other type of pattern has not been observed in previous simulations [13]. However, the patterns obtained here are not observed in the distributions of electron number density n_e, electron temperature T_e, or any other property (Figure 7). Moreover, the ring patterns emerging in our simulation differ from those observed experimentally (e.g., experimental image in Figure 1) [4]. Nevertheless, the spontaneous formation of patterns in our simulations provides insights into the broader dynamics affecting the discharge, urging further refinement and expansion of APGD models to encompass the nuanced intricacies of the near-cathode layer and enhance the overall accuracy of the simulations. The occurrence of patterns is generally associated with the existence of multiple solutions of the model equations (i.e., Table 3), which can also manifest via instabilities. In the context of the present study, instabilities are of primary importance given their predominance in APGDs operating at high currents. Particularly, APGDs are prone to thermal instabilities in the cathode region and tend to stabilize with strong cathode cooling, as demonstrated by Arkhipenko et al. [4]. Similarly, the ring patterns in Figure 8 appear to be mitigated by strong cathode cooling. The above facts underscore the complex interplay among the fidelity of the mathematical model, geometric representation of the computational domain, and numerical accuracy of the simulations to describe pattern formation in APGDs.

To assess the influence of cathode cooling on the overall discharge more clearly, Figure 9 presents the distribution of electron number density n_e and heavy-species temperature T_e along the axis of the domain. The results in Figure 9a show the distribution of n_e for the different values of current I_tot and levels of cathode cooling h_c. The three plots on the right side of Figure 9a are magnified views of the distribution of n_e over the first millimeter of the discharge. In our previous study [13], we noted that the discharge had a regime change between the low-current and the high-current cases, with a threshold of approximately 15 mA. The electron number density for the lowest current cases showed a monotonically increasing trend with the distance from the cathode. In contrast, the higher current cases displayed a local minimum in the center of the discharge and a higher number density at the anode and cathode. In this study, the trend is slightly different. Specifically, for the low-current cases of I_tot = 5 mA, the results show a local minimum for high cathode cooling. Nevertheless, this minimum is not in the center of the discharge, and the global trend of the discharge is still monotonically increasing with the distance from the cathode. The electron number density distribution for I_tot = 20 and 40 mA presents a local minimum similar to that reported in [13]. The slight difference in this trend may be explained by the change in boundary condition at the anode. In [13], the temperature profile was imposed along the anode, while in this study, the heat flux is modeled using a convective cooling model, as shown in Table 4. We associated the trend change in plasma properties as being a marker of competing phenomena during the progressive thermalization of the discharge.

The distribution of heavy-species temperature T_h along the center line of the discharge for all the studied cases is shown in Figure 9b. As previously noted in relation to the results in Figure 5 and Figure 6, the cooling of the cathode has a greater effect with increasing I_tot. While the 5-mA discharge shows a difference in T_h only within the first few millimeters from the cathode, as I_tot increases to 20 and 40 mA; the distribution of T_h drastically changes along the whole discharge. However, all the results for a given current value converge to a single temperature near the anode. The thermalization of the plasma for the higher currents and lower cooling levels seems to lead not only to the change in regime reported in [13], but also to the distribution of a ring-shaped concentrated electric field over the cathode (Figure 8).

4. Conclusions

We developed a self-consistent computational model for the three-dimensional time-dependent simulation of an Atmospheric Pressure Glow Discharge (APGD) with a planar cathode in helium to study the effect of cathode cooling. The mathematical model describes the plasma as a quasi-neutral fluid in thermal and chemical nonequilibrium and is comprised of the conservation equations for total mass, species mass, momentum, heavy-species energy, electron energy, and electric charge. We utilized the model to simulate the 10 mm gap pin-to-plate APGD experimental setup studied by Arkhipenko et al. [4]. The model is used in time-dependent three-dimensional computational simulations of the APGD that do not rely on assumptions of spatial symmetry or steady-state. We perform simulations of APGD operation as a function of two control parameters: total electric current, from 5 to 40 mA, and level of cooling of the cathode, the latter spanning from thermally insulated to extreme cooling.

Simulation results show overall good agreement with experimental measurements of voltage drop, positive column temperature, and current density as a function of discharge current. The results under the highest level of cooling exhibit the closest alignment with the experimental data. The distributions of electron temperature and of electron number density show columnar structures throughout the interelectrode gap, regardless of the cathode cooling level or total current. This behavior can be attributed to the requirement for continuity of electric current through the electrodes. In contrast, the distribution of heavy-species temperature reveals an increasingly pronounced negative gradient towards the cathode with higher cooling levels. This behavior is consistent with the enhanced transfer of thermal energy from the plasma to the cathode as the cathode cooling level increases. A consequence of these behaviors is the formation of a ring pattern in the distribution of the electric field over the cathode for the higher electric current and cathode cooling levels. However, no pattern was observed in simulations employing a similar model that used only one quadrant of the domain in [13]. Importantly, these patterns are absent in the distributions of other properties. This finding highlights the interdependence between the fidelity of the mathematical model, the geometric representation of the computational domain, and the numerical accuracy of the simulations to describe pattern formation in APGDs. Future work will focus on integrating electrode sheath models, incorporating the effects of charged species transport on cathode heat transfer, adopting mesh refinement strategies—particularly in the near-cathode region—and running simulations aimed at capturing pattern formation events and the occurrence of striations (stratification phenomena) along the plasma column.