COMSOL 6.4 - Atmospheric Pressure Corona Discharge in Air

Atmospheric Pressure Corona Discharge in Air

This example presents a study of a coaxial DC corona discharge in dry air at atmospheric pressure. The dimensions and operation conditions are similar to the ones found in electrostatic precipitators with wire-to-plate configuration. The inner electrode has 100 μm radius and the gap between electrodes is 10 cm. The simulations presented are for steady-state regimes with the discharge sustained with 10 kV applied to the inner electrode while the exterior electrode is grounded. Emphasis is placed on the charged particles creation and transport and how that translates into the current-voltage characteristic of the discharge.

Model Definition

Figure 1 shows a cross section of the model geometry. The discharge is assumed to be diffuse and uniform in the radial direction. The model is one dimensional in the radial direction between the electrodes and describes the behavior of charged species using fluid-type equations. We assume that the gas temperature and air number density are constant. The gas temperature is kept at 600 K.

The model solves the electron and ions continuity and momentum equations, in the drift-diffusion approximation, self-consistently coupled with Poisson’s equation. The local field approximation is used, which means that transport and source coefficients are assumed to be well parameterized through the reduced electric field (E/N). In the local field approximation the fluid equation for the mean electron energy is not solved, which reduces significantly the complexity of the numerical problem.

The local field approximation is valid in a situation where the rate of electron energy gain from the electric field is locally balanced by the energy loss rate. When this condition is met the electrons are said to be in local equilibrium with the electric field and the electron mean properties can be expressed as a function the reduced electric field.

The model presented in the following section is used to simulate the ionization of the neutral gas as well as the flux of charged particles when the negative electric potential is applied at the inner conductor (cathode). The high electric field generated by the combination of high potential and small conductor curvature radius (inner conductor, ri) causes electron drift and ionization of the neutral gas surrounding the cathode. The resulting positive ions generate more electrons through secondary emission at the cathode surface. These electrons are accelerated through a small region away from the cathode, where they can acquire significant energy. This can lead to ionization which creates new electron–ion pairs. The secondary ions migrate toward the cathode where they eject more secondary electrons. This process is responsible for sustaining the discharge.

The model uses a Scharfetter–Gummel upwind scheme to eliminate numerical instabilities in the number density of the charged particles associated with the finite element method. This is needed, in particular close to the cathode, where the ion flux is particularly high.

Figure 1: Not-to-scale cross section of the coaxial configuration. The negative potential (-Vin) is applied at the inner conductor (cathode) and the outer electrode is grounded (anode). The shaded area represents the ionization region created by the positive space charge distribution generated in the vicinity of the cathode.

Domain Equations

Electron transport is modeled by solving the continuity equation and the momentum equation under the drift-diffusion approximation (for detailed information on electron transport, see Theory for the Drift Diffusion Interface in the Plasma Module User’s Guide)

The source coefficients in the above equations are determined by the plasma chemistry. The electron rate expression is defined as

where νe,j is the stoichiometric coefficient, and the reaction rate is defined as

where kjf is the forward rate constant and kjr is the reversed rate constant. Both the Electron Impact Reaction feature and Reaction feature can contribute to the electron rate expression. However, when using the Reaction feature it is important to note that the associated electron energy gain or loss is not included in the source term of the electron mean energy equation.

The rate constants can be computed from electron impact cross-section data

where γ = (2q/me)1/2 (SI unit: C1/2/kg1/2), me is the electron mass (SI unit: kg), ε is the electron energy (SI unit: V), σ is the electron impact collision cross section (SI unit: m2), and f is the electron energy distribution function.

When Townsend coefficients are used, the reaction rate is defined as

where αj/Nn is the reduced Townsend coefficient for reaction j (SI unit: m2) and Γe is the electron flux as defined above (SI unit: 1/(m2·s)). Townsend coefficients can increase the stability of the numerical scheme when the electron flux is field driven as is the case with DC discharges.

For heavy species, the following equation is solved for the mass fraction of each species (for detailed information on the transport of the nonelectron species, see Theory for the Heavy Species Transport Interface in the Plasma Module User’s Guide):

The electrostatic field is computed using the following equation:

The space charge density ρ is automatically computed based on the plasma chemistry specified in the model using the formula

For detailed information about electrostatics see Theory for the Electrostatics Interface in the Plasma Module User’s Guide.

Boundary Conditions

Electrons are lost to the wall due to random motion within a few mean free paths of the wall and gained due to secondary emission effects, resulting in the boundary condition for the electron flux

(1)

The second term on the right-hand side of Equation 1 is the gain of electrons due to secondary emission effects, γp being the secondary emission coefficient. For the heavy species, ions are lost to the wall due to surface reactions and the fact that the electric field is directed toward the wall:

(2)

The discharge is driven by a DC electric potential (V0) applied to the inner conductor of the coaxial geometry (at coordinate r = ri). The other boundary (at coordinate r = ro) is grounded. To facilitate the beginning of the numerical simulation a step function is used to modulate V0 with the transient applied potential assuming the form

(3)

This numeric technique does not interfere with the results at steady state, which are the ones of interest in this work.

Plasma Chemistry

The chemistry of a plasma sustained in air can be very complex and a detailed study of the main excited states could easily have hundreds of reactions. The main goal of this model is to study the charge particle density profiles and currents. With that in mind it is used a simplified set of reactions that describes correctly the creation and destruction of charged species in a background of dry air.

Table 1 lists the chemical reactions considered Ref. 2. In the fluid equations nitrogen and oxygen are not treated separately as in a detailed chemistry. Instead a general species A is used for the background gas. A can be ionized forming positive ions p, and A can attach electrons forming negative ions n.

The creation and destruction of electrons in the volume is described by ionization and attachment Townsend coefficients, and by rate constants for a three-body attachment and electron–ion recombination. The Townsend coefficients are obtained as a function of the mean electron energy by suitably averaging over the electron energy distribution computed using a Boltzmann solver with a consistent set of electron scattering collisions cross sections of nitrogen and oxygen Ref. 3. It is used a mixture of 80% nitrogen and 20% oxygen. The relation between the mean electron energy and the reduced electric field is also obtain from the Boltzmann solver and is given in Figure 2.

Figure 2: Mean electron energy as a function of the reduced electric field for a mixture of 80% nitrogen and 20% oxygen.

For detailed information on how to compute source coefficients from a Boltzmann solver see the section The Boltzmann Equation, Two-Term Approximation Interface in the Plasma Module User’s Guide.

Table 1: Table of collisions and reactions modeled.
Reaction	Formula	Type	Δε (eV)	kf (m3/s)
1	e+A=>p+2e	Ionization	15	-
2	e+A=>n	Attachment	-	-
3	e+2A=>n+A	Attachment	-	-
4	e+p=>A	Reaction	-	5·10-14
5	n+p=>2A	Reaction	-	2·10-12

At steady state, the plasma main charged species are ions. For this reason, the initial conditions have an equal density of positive and negative ions and a small density of electrons. These initial conditions preserve charge neutrality as it is important for numerical reasons.

In addition to the volumetric reactions, the following surface reactions are implemented:

Table 2: Table of surface reactions.
Reaction	Formula	Sticking coefficient
1	p=>A	1
2	n=>A	1

When the ions reach the wall, they are assumed to change back to neutral atoms. Note that the secondary emission coefficient for positive ions is set to 0.05 on the cathode boundary (at coordinate r = ri) and to 0 at the outer electrode (at coordinate r = ro). The mean electron energy of the secondary electron is set to 4 eV. When using the local field approximation the mean energy of the secondary electron is only used in postprocessing.

Results and Discussion

Throughout most of this section, results are presented and discussed for a DC negative corona sustained with −45 kV applied to the inner electrode. At the end, the current-voltage characteristic of the discharge are presented. The background gas is kept at a constant density as obtained by the ideal gas law at atmospheric pressure and at a temperature of 600 K. All results presented and discussed correspond to a steady state operation.

When comparing with corona discharges sustained in noble gases, air corona discharges need much higher voltage to breakdown the background gas and to sustain a discharge. There are two main reasons for this: (i) the electron collision frequency in air is higher (in part due to rotational and vibrational interactions) making it more difficult to accelerate electrons; and (ii) oxygen is electronegative.

Figure 3 presents the spatial profiles of the charged species. The discharge can be separated in two regions: (a) corresponding to a region of less than 1 mm near the cathode where most of the ionization occurs; and (b) the rest of the volume that reaches to the ground electrode.

The strongly negative potential at the inner electrode accelerates positive ions toward it and repeals negative charged particles. The result is a region of positive charge separation (region (a)) where strong electric fields exist and the electrons are accelerated to energies capable of ionizing the background gas. The electron temperature, the electric potential, and the reduced electric field are represented in Figure 4, Figure 5, and Figure 6. As can be seen, region (b) has weak electric field and consequent electron temperatures. In this region, electrons have not enough energy to ionize and are efficiently attached forming negative ions. The result is a long spatial portion of the discharge dominated by negative ions that drift toward the ground electrode. Note also that in region (b) the charge separation barely deforms the applied potential.

Figure 7, Figure 8, and Figure 9 are 2D representations (obtained by a revolution of the 1D solution) of the charged species number density. A representation of the charge distribution with distances in linear scale helps to build a more realistic image of this type of discharges. Notice how small the inner electrode and the ionization region are, and how most of the volume is filled with negative ions drifting in the directions of the anode. This observation does not mean that region (a) should be neglected. In fact, it is in region (a) that the mechanisms that sustain the discharge occur.

Comparing with corona discharges in noble gases, this type of corona discharges have a much higher resistivity because the charged particle number density is much lower, and ions, much less mobile than electrons, are the main charge carriers.

Figure 10 shows the absolute value of the total ion current at the ground electrode as a function of the absolute value of the applied voltage at the inner electrode. The current-voltage characteristic follows a quadratic law as expected. The obtained values of the current density are also coherent with the ones found in this type of discharges.

Figure 3: Spatial profiles of the charged species number density at steady state: electrons (blue), positive ions (green), and negative ions (red).

Figure 4: Spatial profile of the electron temperature.

Figure 5: Spatial profile of the electric potential.

Figure 6: Spatial profile of the reduced electric field.

Figure 7: 2D representation of the electron density (the scale is in log base 10).

Figure 8: 2D representation of the negative ion density (the scale is in log base 10).

Figure 9: 2D representation of the positive ion density (the scale is in log base 10).

Figure 10: Total ion current density (absolute value) at the ground electrode as a function of the applied voltage (absolute value) at the inner electrode.

References

1. M.A. Lieberman and A.J. Lichtenberg, Principles of Plasma Discharges and Materials Processing, John Wiley & Sons, 2005.

2. A.A. Kulikovsky, “Positive Streamer Between Parallel Plate Electrode in Atmospheric Pressure Air,” J. Phys. D: Appl. Phys., vol. 30, pp. 441–450, 1997.

3. Phelps database, www.lxcat.net, retrieved 2016.

Plasma_Module/Corona_Discharges/corona_discharge_air_1d

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Follow the steps below to create the model geometry: a simple 1D geometry consisting of a single domain bounded by the cathode (left, inner conductor) and the anode (right, outer conductor).

In the window, under click .

In the window for , locate the section.

From the list, choose .

Interval 1 (i1)

Right-click > and choose .

In the window for , locate the section.

In the table, enter the following settings:


Coordinates (cm)
10
0.01

Click

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
mueN	3.74e24(plas.Erd1e21)^-0.22[1/(Vms)]		Reduced electron mobility
muiN	6e21[1/(Vsm)]	1/(V·m·s)	Reduced ion mobility
rnp	2e-6[cm^3/s]	m³/s	ion-ion recombination
rei	5e-8[cm^3/s]	m³/s	electron-ion recombination
Vapp	-V0*ramp		Applied Voltage
ramp	tanh(1e5*t)
p0	760[torr]	Pa	Gas pressure
t0	600[K]	K	Gas temperature
ni0	1e17[m^-3]	1/m³	Initial ion number density
ne0	1e10[m^-3]	1/m³	Initial electron number density