Transition mechanism of negative DC corona modes in atmospheric air: from Trichel pulses to pulseless glow

A widely-used three-species hydrodynamic corona model is employed in which the dynamics of electrons, positive ions and negative ions are taken into consideration. The number densities of the charged species are calculated by solving the following set of continuity equations (1)–(3) coupled with a Poisson equation (4).

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {n}_{{\rm{e}}}}{\partial t}+{\rm{\nabla }}\cdot (-{\mu }_{{\rm{e}}}\vec{E}{n}_{{\rm{e}}}-{D}_{{\rm{e}}}{\rm{\nabla }}{n}_{{\rm{e}}})=\alpha {n}_{{\rm{e}}}\left|{\mu }_{{\rm{e}}}\vec{E}\right|-\eta {n}_{{\rm{e}}}\left|{\mu }_{{\rm{e}}}\vec{E}\right|\\ \,-{k}_{{\rm{e}}{\rm{p}}}{n}_{{\rm{e}}}{n}_{{\rm{p}}}+{k}_{{\rm{d}}{\rm{t}}}{n}_{{\rm{n}}}N\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {{\rm{n}}}_{{\rm{p}}}}{\partial t}+{\rm{\nabla }}\cdot ({\mu }_{{\rm{p}}}\vec{E}{n}_{{\rm{p}}}-{D}_{{\rm{p}}}{\rm{\nabla }}{n}_{{\rm{p}}})=\alpha {n}_{{\rm{e}}}\left|{\mu }_{{\rm{e}}}\vec{E}\right|-{k}_{{\rm{n}}{\rm{p}}}{n}_{{\rm{n}}}{n}_{{\rm{p}}}\\ \,-{k}_{{\rm{e}}{\rm{p}}}{n}_{{\rm{e}}}{n}_{{\rm{p}}}\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {n}_{{\rm{n}}}}{\partial t}+{\rm{\nabla }}\cdot (-{\mu }_{{\rm{n}}}\vec{E}{n}_{{\rm{n}}}-{D}_{{\rm{n}}}{\rm{\nabla }}{n}_{n})=\eta {n}_{{\rm{e}}}\left|{\mu }_{{\rm{e}}}\vec{E}\right|-{k}_{{\rm{n}}{\rm{p}}}{n}_{{\rm{n}}}{n}_{{\rm{p}}}\\ \,-{k}_{{\rm{d}}{\rm{t}}}{n}_{{\rm{n}}}N\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi =-\displaystyle \frac{{e}({{\rm{n}}}_{{\rm{p}}}-{{\rm{n}}}_{{\rm{e}}}-{{\rm{n}}}_{{\rm{n}}})}{\varepsilon }\end{eqnarray} \tag{ 4 }$$

where n_e, n_p and n_n are the electron, positive ion and negative ion number densities; μ_e, μ_p, μ_n, and D_e, D_p and D_n the mobility and diffusion coefficients for electrons, positive ions and negative ions respectively; α the Townsend ionization coefficient, η the attachment coefficient, k_dt the detachment coefficient, and k_ep and k_np the recombination coefficients of positive ions with electrons and negative ions, respectively. The detachment only accounts for the electrons from ${{{\rm{O}}}_{2}}^{-}$ ions and the coefficient was summarized in [33]. The charged particle mobility, ionization and attachment coefficients are supposed to depend on the local reduced field $\left|\vec{E}\right|/N,$ and these values are taken from Morrow's corona model [34]. Here N is the gas number density, which is related to the mass density ρ_g as ρ_g = m_gN, where m_g is the average mass of the air molecules. The electric field $\vec{E}$ is calculated according to $\vec{E}=-{\rm{\nabla }}\phi ,$ where ϕ is the electric potential; ε is the permittivity of air and e the electron charge.

The discharge current including the contributions of all charged species is proposed by Sato [35] and given by equation (5):

$$\begin{eqnarray}&&\begin{array}{l}I=\displaystyle \frac{e}{{U}_{{\rm{a}}}}\displaystyle {\iint }_{D}2\pi r\left({\mu }_{{\rm{p}}}\vec{E}{n}_{{\rm{p}}}-{D}_{{\rm{p}}}{\rm{\nabla }}{n}_{{\rm{p}}}-{\mu }_{{\rm{e}}}\vec{E}{n}_{{\rm{e}}}-{D}_{{\rm{e}}}{\rm{\nabla }}{n}_{{\rm{e}}}\right.\\ \,\,-\left.{\mu }_{{\rm{n}}}\vec{E}{n}_{{\rm{n}}}-{D}_{{\rm{n}}}{\rm{\nabla }}{n}_{{\rm{n}}}\right)\cdot {\vec{E}}_{{\rm{L}}}{\rm{d}}r{\rm{d}}z\end{array}\end{eqnarray} \tag{ 5 }$$

where r and z are cylindrical coordinates, U_a the applied voltage, ${\vec{E}}_{{\rm{L}}}$ the Laplacian electric field and D the 2D axisymmetric domain.

The gas dynamics is influenced by the discharge through momentum and energy transfer. The compressible Navier–Stokes equations (6)–(8), including the conservation of mass, momentum and energy for neutral gas, are employed to describe the coupling between discharge and gas dynamics.

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}\vec{u})=0\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\rho }_{{\rm{g}}}\vec{u})}{\partial t}={\rm{\nabla }}\cdot \left[{\rho }_{{\rm{g}}}\vec{u}\otimes \vec{u}+p{\bf{I}}-\tau \right]=\mathop{{f}_{{\rm{e}}{\rm{h}}{\rm{d}}}}\limits^{\longrightarrow}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{{\rm{n}}}}{\partial t}+{\rm{\nabla }}\cdot \left[({E}_{{\rm{n}}}+p)\mathop{u}\limits^{\longrightarrow}-\kappa {\rm{\nabla }}T-\tau \cdot \mathop{u}\limits^{\longrightarrow}\right]={P}_{{\rm{heat}}}.\end{eqnarray} \tag{ 8 }$$

Here, ρ_g is the gas density, $\vec{u}$ the flow velocity, $\tau \,={\mu }_{v}\left[{\rm{\nabla }}\vec{u}+{({\rm{\nabla }}\vec{u})}^{{\rm{T}}}-2/3({\rm{\nabla }}\cdot \vec{u}){\bf{I}}\right]$ the stress tensor, μ_v the air dynamic viscosity, p the static pressure, I the unit tensor, T the gas temperature, κ the thermal conductivity of air, and ${E}_{{\rm{n}}}=p/(\gamma -1)+1/2{\rho }_{{\rm{g}}}\mathop{u}\limits^{\longrightarrow}\cdot \mathop{u}\limits^{\longrightarrow}$ is the total energy per unit volume, where γ = 1.4 is the ratio of specific heat for air.

The source term $\mathop{{f}_{{\rm{e}}{\rm{h}}{\rm{d}}}}\limits^{\longrightarrow}$ in the momentum equation (7) is the electrohydrodynamic (EHD) force [36], which is responsible for the momentum transfer from charged particles to neutral molecules and is written as:

$$\begin{eqnarray}&&{f}_{{\rm{e}}{\rm{h}}{\rm{d}}}={\rm{e}}({n}_{{\rm{p}}}-{n}_{{\rm{e}}}-{n}_{n})\vec{E}.\end{eqnarray} \tag{ 9 }$$

The source term P_heat in the energy equation (8) is the heating source of the gas. The energy released through collisions between ions and neutral species is assumed to be instantaneously transferred into gas heating. The thermal energy transferred from electrons to neutral molecules can be divided into several parts. Firstly, the energy transferred though elastic collisions and rotational excitation is immediately deposited into gas heating. Secondly, a fraction of the electron energy is expended on the excitation of the electronic degrees of freedom of N₂ molecules and instantaneously converted into gas heating [37]. Thirdly, the relaxation of vibrational energy stored in N₂ molecules occurs on a longer time scale. A relaxation time τ_vt is used for the process of vibration–translation (V–T) collisions. So the electrical power density deposited by the discharge is

$$\begin{eqnarray}&&{P}_{{\rm{h}}{\rm{e}}{\rm{a}}{\rm{t}}}=({\vec{j}}_{{\rm{p}}}-{\vec{j}}_{{\rm{n}}})\cdot \vec{E}-{\eta }_{{\rm{e}}{\rm{l}}-{\rm{r}}}{\vec{j}}_{{\rm{e}}}\cdot \vec{E}-{\xi }_{{\rm{e}}{\rm{x}}}{\eta }_{{\rm{e}}{\rm{x}}}{\vec{j}}_{{\rm{e}}}\cdot \vec{E}+{P}_{{\rm{v}}{\rm{t}}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {P}_{{\rm{v}}{\rm{t}}}}{\partial {t}}+\displaystyle \frac{{P}_{{\rm{v}}{\rm{t}}}}{{\tau }_{{\rm{v}}{\rm{t}}}}=\displaystyle \frac{-{\eta }_{{v}}{\vec{j}}_{{\rm{e}}}\cdot \vec{E}}{{\tau }_{{\rm{v}}{\rm{t}}}}\end{eqnarray} \tag{ 11 }$$

where j_e, j_p and j_n are the current densities of electrons, positive ions and negative ions, η_el-r, η_ex and η_v are the fractional power deposited in elastic and rotational excitation of air molecules, electronic excitation and vibrational excitation of N₂ molecules respectively. The dependence of the three fractions on $\left|\vec{E}\right|/N$ is calculated using BOLSIG+ for discharges in air [38]. It is assumed that ξ_ex = 30% of the power expended on the electronic excitation is directly transferred to gas heating [39]. A phenomenological equation (11) of the V–T energy transfer is used in this paper [40]. The power density P_vt is first stored in vibrational states of the molecules and then released as heat with a relaxation time τ_vt = 20 μs [41]. The quantity τ_vt depends on the vibrational temperature [42, 43], but we treat it as a parameter and will later discuss it in section 4.1.

The above equations are closed with the equation (12) of state of perfect gases:

$$\begin{eqnarray}&&p=N{k}_{B}T=\displaystyle \frac{{\rho }_{g}}{{m}_{g}}{k}_{B}T.\end{eqnarray} \tag{ 12 }$$

The discharge model and the gas dynamic model are two-way coupled. The former model provides the EHD force and the heating power to the latter. In turn the changes in the gas density by gas heating have an impact on the transport coefficients of charged species, the ionization and attachment rates. In addition, the induced gas flow may modify the velocity of charged particles. However, the flow velocity is several meters per second and is much lower than the typical drift velocities of ion and electron, thereby hardly affecting the discharge dynamics in contrast to the gas heating.

In order to compare with experimental results, the negative corona model is investigated in a similar 2D axisymmetric geometry [31]. Figure 1 shows the schematic diagram of the needle-cylinder electrode arrangement. The needle has a hyperbolic tip with a radius of curvature of about 220 μm, and is placed 20 mm away from and collinearly to the cylinder electrode symmetry axis. The grounded electrode is a cylinder with an inner diameter of 34 mm and a height of 30 mm.

Zoom In Zoom Out Reset image size

The initial electrons and positive ion are uniformly distributed in the simulation domain and their number densities are equal to 10¹⁰ m⁻³. The density of negative ions is assumed to be zero. The initial gas velocity is 0 m s⁻¹, the initial pressure 1 bar and the initial temperature 293 K. The boundary conditions for the discharge model and the flow model can be found in [31, 32]. Secondary emitted electrons from the needle surface are required in the model to sustain the discharge. The secondary electron emission coefficient γ_se is set to 0.001, as we suggested in [31]. The boundary conditions for the heat transfer equation (8) are set as follows: we assume that the temperature on the needle electrode has a homogeneous Neumann boundary condition and that it is equal to the ambient temperature (293 K) on other boundaries.

The fully coupled corona discharge model is implemented in the commercial FEM-based package COMSOL Multiphysics. The discharge equations (1)–(4) are solved by the 'transport of diluted species' and 'electrostatics' modules, the gas dynamics equations (6)–(8) are solved by the 'Laminar Flow' and 'Heat Transfer in Fluids' modules, and the V–T transfer equation (11) is solved by the ordinary differential equation (ODE) module. The time step size is automatically determined with a maximum step of 10⁻⁸ s. In order to accelerate the calculation, the mesh is refined near the tip and is left coarser in the remaining region. A non-uniform triangular mesh was used with a maximum element size of 1.0 μm near the needle tip, which increases to 1 mm in the area far from the needle. In total around 80 000 triangular elements are used.

Firstly, we show the current waveforms calculated by the NCD-B model at voltages of −12 kV, −14 kV and −16 kV in figure 2. During the rising edge of the voltage, a large current pulse is seen for all voltages. This is because of rapid ionization and electron motion during the first several microseconds after the DC voltage is suddenly applied. The amplitude of the first current pulse is usually very high and reaches many hundreds of micro-amperes. After the first current pulse, the current gradually decreases. Since many electrons are generated and move away from the tip, they will be attached to oxygen molecules and negative ions drift away in the low fields region. At −12 kV the regular current pulses start to oscillate from t = 22 μs, indicating that the Trichel pulse mode of negative corona has started. The oscillating behavior appears later at −14 kV than at −12 kV. This implies that more negative charges are produced at −14 kV, and more time is required to clear the charges to recover the electric field near the cathode. The pulse repetition frequency at −14 kV is higher frequency than at −12 kV. When the voltage is increased to −16 kV, the current pulses do not appear until t = 88 μs, which is not shown in figure 2. The corona mode does not change to the pulseless glow stage in the voltage range of −12 to  −16 kV. Though we have not tested higher voltages to see the transition threshold voltage, the simulation results by NCD-B indicates that it does not describe the transition process of negative corona modes accurately.

Zoom In Zoom Out Reset image size

Figure 3 shows the current waveforms calculated by NCD-DT when the voltages are −12 kV and −14 kV. Similar to the calculated current by NCD-B, the current consists of a regular train of pulses at voltages of −12 kV and −14 kV. At −14 kV, the current pulses require more time to reach a steady state than in the NCD-B model results. The calculation results show that the mode transition does not occur at lower voltages when electron detachment is considered. These findings are not consistent with the results by Dordizadeh et al [30], who report that the pulsating mode disappears when detachment of electrons from O⁻ ions by collisions with unexcited N₂ is considered. The set of the coefficients of detachment reaction were taken from Soloviev et al [44]. However, this detachment practically does not occur when the gas temperature is not too high [33]. A different set of reactions is used in our model, where the detachment of electrons from ${{{\rm{O}}}_{2}}^{-}$ ions are considered. Moreover, the detachment rate coefficient from O⁻ ions is much higher than that from ${{{\rm{O}}}_{2}}^{-}$ ions. This will probably lead to more electrons and the formation of a glow. We will further discuss how the discharge mode changes for varying detachment rate in section 4.2.

Zoom In Zoom Out Reset image size

The simulated current by NCD-GH with voltages of −12 kV and −14 kV is shown in figure 4. The current at −12 kV shows pulsating behavior while at −14 kV the discharge becomes steady after the first large current pulse. This indicates that the transition from Trichel pulses to the pulseless glow occurs between −12 kV and −14 kV. This transition voltage is in good agreement with the experimental value [31]. The current pulses appear later at t = 156 μs, and the amplitudes of the pulses are tens of micro-amperes, lower than those by NCD-B and NCD-DT.

Zoom In Zoom Out Reset image size

In figure 5 the mean currents of the Trichel pulses are compared for the three models and the experimental data. The mean currents at −12 kV and −14 kV during the last current pulse are about 20 μA and 30 μA respectively, which is larger than the measured values of 15 μA and 24 μA. In spite of different corona modes at −14 kV, the differences in mean current calculated by the three models are small. This suggests that the ions in the drift region are the main contribution to the total current, while the ionization region, which is affected most by gas heating, is responsible for the corona modes. This will be clearly illustrated in the next section when the reduced electric fields are compared.

Zoom In Zoom Out Reset image size

The 2D distributions of the electric field and net number density of charged species calculated by NCD-B were studied thoroughly in the Trichel pulse regime [29, 31]. It can be seen from the previous section that NCD-B and NCD-DT give similar current waveforms. This suggests that the distributions of the electric field and charged species are also quite similar for these two models. Figure 6 shows the time evolution of these distributions to illustrate the Trichel pulse formation. During a Trichel pulse, the steep rise of the current pulse is mainly due to rapid motion of positive ions in the high field near the needle tip, which looks like an ionization wave propagating towards the cathode. As the positive ions are neutralized by the cathode, the electric field adjacent to the cathode is strongly distorted by leftover space charges. A large amount of electrons at a further distance from the cathode is attached to molecules and negative ions are formed. This could reduce the field near the cathode and quench the ionization activity. As the negative ions drift away, the current decreases gradually to a low value due to the low mobility of ions. The maximum charged species density is not located on the symmetry axis. It means that a ring-like corona region is formed near the tip, which has also been reported in other simulations of Trichel pulses [45]. Moreover, a ring shape of the negative corona discharge was observed from a point electrode [46, 47]. The authors believed that it is due to the field distortion caused by negative space charges in the gap.

Zoom In Zoom Out Reset image size

Figure 7 presents results of the time evolution of the electric field and charged species density calculated by NCD-GH. The basic process of Trichel pulse formation from the NCD-GH model is similar to that from the NCD-DT model. There are still a few differences compared to figure 6 with respect to the Trichel pulse regime in figure 7(a). The positive ions are seen to move towards the needle more clearly in figure 6, while the negative ion clouds move away faster in figure 7. The maximum number density of charged species for NCD-GH is higher than that for NCD-DT. Figure 7(b) depicts the distributions of the electric field and charged species density of the glow mode at U = −14 kV. The distributions are stationary and no rapid motion of positive ions is seen as in figure 7(a).

Zoom In Zoom Out Reset image size

Figure 8 compares the distributions of gas density and gas temperature at the end of the simulation. The power density deposited in the cathode sheath region is quite large and leads to the gas heating. The maximum gas temperature in this region at U = −14 kV is 655 K, compared to 543 K at U = −12 kV. The calculated maximum temperatures are higher than the experimental results of 400–450 K [32], which were obtained from the spectroscopic measurements of the N₂ C–B transition. The difference is mainly because the temperature is measured at the region below the needle tip and is spatially averaged. Accordingly, the gas density at U = −14 kV is lower than at U = −12 kV. This will further enhance the reduced electric field near the tip and thus affect the ionization and attachment coefficients.

Zoom In Zoom Out Reset image size

Figure 9 shows the 1D distribution of the reduced field and ionization and attachment rates along the symmetry axis for NCD-B and NCD-GH model at U = −14 kV. The gas heating primarily influences these parameters in the region z < 100 μm. The reduced electric field calculated by NCD-GH near the needle can reach a maximum value of 1200 Td, which is triple that of NCD-B, while in the area far from the needle the reduced fields do not differ much. The electric fields do not have to be strong enough through the whole gap to reach a high net ionization rate, as discussed in [30], since the decreased gas density leads to higher reaction rates. In the ionization region at z = 100 μm, the ionization rate by NCD-GH is nearly two orders of magnitude greater than that by NCD-B, and so is the attachment rate.

Zoom In Zoom Out Reset image size

It can be concluded from the above analysis that gas heating has effects on the negative mode transition. If gas heating was not taken into consideration, a higher voltage would be required to enhance the electric field in the whole gap and thus realizing the Trichel pulse—pulseless glow transition, as is discussed in [29]. The high electric field not only increases the rate of removal of negative ions into the low field gap, but also reduces the formation rate of negative ions. However, according to the results by NCD-B, this threshold voltage is so high that the discharge may quickly turn into the spark regime. That is why the removal of heat dissipated in the discharge is crucial to form glow mode in a short air gap at higher voltages [48]. The experiment results show that the negative corona can switch to a pulseless glow at a relative low voltage. This phenomenon can be simulated by NCD-GH. The gas heating decreases the gas density and thereby increases the reduced electric field. Therefore, the ionization and attachment rates are increased, which has a similar effect to raising the voltage. The difference is that only the cathode sheath region in the vicinity of the needle is affected instead of the entire discharge area. On the other hand, the detachment is not that important for the transition in our simulation, which is inconsistent with the findings by Dordizadeh et al [30]. The main reason is that the detachment rate is 3–4 orders of degree lower in our work. A higher detachment rate leads to suppression of negative ion density. Hence, the electric field will not be decreased to the same degree as the lower detachment rate. Besides, the extra secondary electrons in the vicinity of the tip will accelerate the ionization reaction and keep the charged particles in an equilibrium state. It can be concluded that the gas heating effect plays a greater role in the transition process than the detachment rate.

A phenomenological gas heating model is used in this paper. The energy relaxation time τ_vt is used to describe how fast the energy stored in the vibrational states can be transferred to the translational state. Hence, the gas temperature might be affected and so does the reduced electric field. The period of the Trichel pulses is several microseconds, which is lower than the τ_vt of 20 μs in the NCD-GH model. Here we investigate its effect by varying τ_vt from 2 μs to 200 μs. For τ_vt = 200 μs only a small fraction of the vibrational energy is transferred to translational energy, since τ_vt is now comparable to the total simulation time. Figure 10 shows the distributions of the reduced electric field and gas temperature on the symmetry axis for different values of τ_vt. The reduced electric fields are hardly affected by the value of τ_vt and they all decrease from about 950 Td to 10 Td for increasing distance to the tip. The reduced field is slightly higher near the needle tip for a lower τ_vt. The gas temperature distributions show similar trends and decrease from 540 K to the ambient temperature 293 K. The result shows that a smaller τ_vt leads to faster V–T energy transfer, whereas it enhances the reduced electric field only slightly.

Zoom In Zoom Out Reset image size

Next, figure 11 shows the time evolution of the discharge current for different values of τ_vt. For τ_vt = 200 μs the current goes into the Trichel pulse stage earlier than for τ_vt = 2 μs and 20 μs. This is because the energy is deposited on a shorter time scale and therefore the ionization activity is stronger, and a larger current pulse appears for smaller τ_vt, as can be seen in figure 11(a). It requires a longer time for negative ions to move from the tip. After the first large current pulse, the stable current pulses in figure 11(b) for different τ_vt are very similar. This can be explained by the similar reduced electric field and gas temperature in the vicinity of the tip. The simulation results demonstrate that τ_vt affects the time of Trichel–pulseless glow transition, but not the characteristics of the Trichel pulses.

Zoom In Zoom Out Reset image size

In section 3, it was shown that the addition of the detachment reaction in the NCD-DT model does not affect the transition of negative corona modes. So, we increase the detachment rate in order to better investigate its effect. Figure 12 shows the current waveforms calculated for different detachment rates of k_dt, 10² k_dt and 10⁴ k_dt. At −14 kV, the current calculated with 10² k_dt is close to that of the original k_dt. When the detachment rate is increased to 10⁴ times its original value, the current pulse disappears and the corona modes has changed into the pulseless glow mode. It should also be noted that there are small oscillations in the current waveform, which is not as smooth as that in figure 4. The calculation results illustrate that the higher detachment reaction increases the free electron density and thus reduces the negative ion density. The electrons can move away from the tip faster, so the electric field will not decrease to choke the ionization activity. However, the detachment rate used here has an unreasonably high value, while in the real case the effect of the detachment reaction on the negative corona modes transition cannot be compared to the gas heating.

Zoom In Zoom Out Reset image size