COMSOL 6.4 - DC Glow Discharge Coupled with the Two-Term Boltzmann Equation

DC Glow Discharge Coupled with the Two-Term Boltzmann Equation

DC glow discharges in the low pressure regime have long been used for gas lasers and fluorescent lamps. DC discharges are attractive to study because the solution is time independent. This model shows how to use the Plasma interface to set up an analysis of a positive column. The discharge is sustained by emission of secondary electrons at the cathode.

This tutorial models a DC discharge by solving plasma fluid equations fully coupled with the homogeneous and time-independent electron Boltzmann equation in the classical two-term approximation. The approximated Boltzmann equation is solved for each position of space and is coupled with the fluid equations by way of the electron mean energy. The rate constants of electron impact reactions and the electron transport parameters are obtained by suitable integration of the computed electron energy distribution function over electron scattering cross sections. Simulation results from this model are compared with results from the model DC Glow Discharge, 1D where the electron energy distribution function (EEDF) is assumed to be Maxwellian.

Model Definition

The DC discharge consists of two electrodes, one powered (the anode) and one grounded (the cathode):

Figure 1: Schematic of the DC discharge. The voltage applied across the electrodes leads to formation of a plasma.

domain equations

Electron transport is modeled by solving the continuity equation, the momentum equation under the drift-diffusion approximation, and the mean electron energy equation (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.

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.

The total electron energy loss or gained is calculated by summing the collisional energy changes from all reactions defined with the Electron Impact Reaction feature as

where Δεj is the energy loss from reaction j (SI unit: V) and F is the Faraday constant (SI unit: C/mol). For excitation and ionization collisions Δεj corresponds to the energy of the excited state being excited/deexcited or ionized, for attachment Δεj is set to zero, and for elastic collisions

where me and mk are the electron and heavy species mass in kg, Te is the electron temperature in eV, and Tgas is the gas temperature in K. 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.

Boltzmann Equation in the Two-Term Approximation

The EEDF used in this model is obtained from the solution of the homogeneous and time-independent electron Boltzmann equation in the classical two-term approximation

where F0 is the isotropic part of an EEDF constant in time and space that verifies the following normalization:

The different terms are presented below:

The following definitions apply:

Here:

•

γ = (2q/me)1/2 (SI unit: C1/2/kg1/2)

•

me is the electron mass (SI unit: kg)

•

ε = (v/γ)2 is energy (SI unit: V)

•

σε is the total elastic collision cross section (SI unit: m2)

•

σm is the total momentum collision cross section (SI unit: m2)

•

is the normalized total momentum collision cross section (SI unit: m2)

•

q is the electron charge (SI unit: C)

•

ε0 is the permittivity of free space (SI unit: F/m)

•

T is the temperature of the background gas (SI unit: K)

•

kb is the Boltzmann constant (SI unit: J/K)

•

ne is the electron density (SI unit: 1/m3)

•

Nn is the background gas density (SI unit: 1/m3)

•

λ is a scalar-valued renormalization factor that ensures that the EEDF is normalized to 1as explained above. An ODE is implemented to solve for the value of λ.

•

M is the mass of the target species (SI unit: kg).

The source term, S represents energy loss due to inelastic collisions. Because the energy loss due to an inelastic collision is quantized, the source term is nonlocal in the energy space. The source term can be decomposed into four parts where the following definitions apply:

where xk is the mole fraction of the target species for reaction k, σk is the collision cross section for reaction k, Δεk is the energy loss from collision k, and δ is the delta function at ε = 0.

The mean electron energy is defined by the integral

The external excitation in the Boltzmann equation comes from an electric field. If the Boltzmann equation is to be solved for a given mean electron energy (as in the present model) a Lagrange multiplier is introduced to solve for the reduced electric field, such that the following equation, presented in the weak form of the constraint, is satisfied:

where the tilde denotes a test function. For detailed information about electrostatics see Theory for the Boltzmann Equation, Two-Term Approximation Interface in the Plasma Module User’s Guide.

Coupling Between the Fluid Model and the Boltzmann Equation

The fluid model and the Boltzmann equation are solved fully coupled. The Boltzmann equation is solved at every point in space using as input the space dependent mean electron energy and the mole fraction of the different heavy species obtained from the fluid model. The EEDFs obtained from the solution of the Boltzmann equation are used to compute macroscopic rate constants and transport coefficients intervening in the fluid model equations thus closing the loop.

Electron Rate Coefficients, Townsend Coefficients, and Transport Parameters

Rate coefficient for electron impact reaction are computed using cross section data and computed EEDFs by the integral

The Townsend coefficients are computed using

where νk is the frequency associated with the reaction rate kk

and P is the power absorbed by the electrons from the DC electric field. The reduced electron mobility is computed using

boundary conditions

Unlike RF discharges, the mechanism for sustaining the discharge is emission of secondary electrons from the cathode. An electron is emitted from the cathode surface with a specified probability when struck by an ion. These electrons are then accelerated by the strong electric field close to the cathode where they acquire enough energy to initiate ionization. The net result is a rapid increase in the electron density close to the cathode in a region often known as the cathode fall or Crookes dark space.

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 following boundary condition for the electron flux:

(1)

and the electron energy flux:

(2)

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. The second term in Equation 2 is the secondary emission energy flux, εp being the mean energy of the secondary electrons. 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:

plasma chemistry

Argon is one of the simplest mechanisms to implement at low pressures. The electronically excited states can be lumped into a single species which results in a chemical mechanism consisting of only 3 species and 7 reactions (electron impact cross sections are obtained from Ref. 3):

Table 1: table of collisions and reactions modeled.
Reaction	Formula	Type
1	e+Ar=>e+Ar	Elastic	0
2	e+Ar=>e+Ars	Excitation	11.5
3	e+Ars=>e+Ar	Superelastic	-11.5
4	e+Ar=>2e+Ar+	Ionization	15.8
5	e+Ars=>2e+Ar+	Ionization	4.24
6	Ars+Ars=>e+Ar+Ar+	Penning ionization	-
7	Ars+Ar=>Ar+Ar	Metastable quenching	-

In this discharge, the electron density and density of excited species is relatively low so stepwise ionization is not as important as in high density discharges. In addition to volumetric reactions, the following surface reactions are implemented:

Table 2: table of surface reactions.
Reaction	Formula	Sticking coefficient
1	Ars=>Ar	1
2	Ar+=>Ar	1

When a metastable argon atom makes contact with the wall, it reverts to the ground state argon atom with some probability (the sticking coefficient).

Results and Discussion

In this section, results are presented from the model. For reference, results of the fluid model when assuming a Maxwellian EEDF are also presented.

Figure 2 through Figure 7 present the electron density, the electron temperature, the potential, the mass fraction of the argon excited state, the current densities of the charged species, and the electron transport parameters along the discharge. All quantities presented change significantly when the EEDF is computed. In summary: the electron density is more than 50% smaller and the maximum deviates toward the cathode; the electron temperature profile is similar but it is 1 eV higher corresponding to an increase of 50%; the mass fraction of the argon excited is considerable different with a much more pronounced maximum near the cathode and a much smaller population in the positive column; the different current densities maintain roughly the same profile and magnitude with a difference of only 10%.

Perhaps the most important difference between the two models is the fact that the creation of excited species and ions in the positive column is much smaller for the case when the EEDF is computed. This can be verified from the plots of the source terms of these species (not presented). For the ions, the source term can be verified to be close to zero since the ion current above 20 cm is practically constant meaning that the source term is zero (the divergence of the ion flux is zero).

At a first look this might seem incoherent with the fact that the electron temperature increases for the case of the computed EEDF. However, this can be understandable by looking at the shape of the computed EEDFs in the region of the positive column presented in Figure 8. The EEDFs deviates considerable from a Maxwellian with the high energy region above the excitation and ionization thresholds being strongly unpopulated. EEDFs with these shapes when integrated over the excitation or ionization cross section give small rate coefficients.

When modeling a DC discharge one might find that using a computed EEDF gives better agreement with experiments. By inspecting the EEDFs the use of an analytic EEDF with a depleted tail, representative of the EEDFs in positive column, might not be adequate since the EEDF in the cathode region has a different shape closer to a Maxwellian. For this case, a good option to avoid large computational times is to give to the fluid model a previously computed EEDF in the form of a lookup table as a function of the mean electron energy.

Figure 2: Electron density obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 3: Electron temperature obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 4: Potential obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 5: Mass fraction of the excited argon obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 6: Electron current density (blue), the ion current density (green), and the total current density (red) obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 7: Reduced electron mobility (blue) and diffusivity (green) electron current density (blue), the ion current density (green) and the total current density (red) obtained with the present model that computes the EEDF and from the same model where the EEDF was assumed to be Maxwellian.

Figure 8: Computed EEDFs in the positive column region (1.9 cm) and in the cathode fall region (18 cm). For reference, it is presented a Maxwellian EEDF with the mean energy obtained from the model with the computed EEDF.

References

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

2. G.J.M. Hagelaar and L.C. Pitchford, “Solving the Boltzmann Equation to Obtain Electron Transport Coefficients and Rate Coefficients for Fluid Models,” Plasma Sources Science and Technology, vol. 14, pp. 722–733, 2005.

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

Plasma_Module/Space-Dependent_EEDF_Modeling/positive_column_1d_boltzmann

Modeling Instructions

Root

The following instructions show how to create a model that solves plasma–fluid-type equations fully coupled with the Boltzmann equation in the two-term approximation.

Application Libraries

The DC positive column application-library model is used as an example

From the menu, choose .

In the window, select > > in the tree.

Click

Choose to solve the Boltzmann equation in the two-term approximation, set the number of elements in the extra dimension to 50, and choose to compute the maximum energy automatically.

Component 1 (comp1)

In the window, expand the node.

Plasma (plas)

In the window, expand the > node, then click .

In the window for , locate the section.

From the list, choose .

In the N text field, type 50.

Select the checkbox.

Click to expand the section. From the list, choose .

Change the two reactions that use Townsend coefficients from lookup tables to use cross-section data and select . This way, the Townsend coefficients are computed using the electron impact cross section and the computed EEDF.

2: e+Ar=>e+Ars

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

4: e+Ar=>2e+Ar+

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

Plasma Model 1

In the window, click .

In the window for , locate the section.

From the list, choose . This way, the electron mobility is obtained from suitable integration of the EEDF over the cross section for momentum transfer. The remaining transport parameters are computed using the electron mobility.

Maxwellian EEDF

In the window, click .

In the window for , type Maxwellian EEDF in the text field.

Add Study

Add an study to solve for the EEDF only. Later, the solution of this study will serve as the initial condition for the study that solves the fluid-type equations coupled with the Boltzmann equation.

In the toolbar, click