COMSOL 6.4 - Thermal Plasma

Low pressure discharges are characterized by the fact that the electron temperature is much higher than the neutral gas temperature. As the gas pressure increases, the number of collisions between the electrons and neutrals increases. At high enough pressures the electron temperature becomes equal to the gas temperature. At this point the plasma is in local thermodynamic equilibrium and a much simpler MHD model can be used to model the plasma.

This model simulates a plasma at medium pressure (1 torr), where the gas temperature cannot be assumed to be constant but the plasma is still not in local thermodynamic equilibrium. In Figure 1 the electron (blue) and gas (black) temperatures are plotted as a function of pressure. At low pressures the two temperatures are decoupled, but as the pressure increases the temperatures tend toward the same limit. There are no axes on the plot since the exact temperature and pressure depend strongly on the gas in question.

This application requires the Plasma Module and AC/DC Module.

Figure 1: Plot of electron (blue) and gas (black) temperature versus pressure. At higher pressures the two temperatures become equal.

Model Definition

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.

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.

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.

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. 2):

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	-

Stepwise ionization (reaction 5) can play an important role in sustaining low pressure argon discharges. Excited argon atoms are consumed via superelastic collisions with electrons, quenching with neutral argon atoms, ionization or Penning ionization where two metastable argon atoms react to form a neutral argon atom, an argon ion and an electron. Reaction number 1, elastic collisions with electrons is primarily responsible for heating of the gas. 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).

electrical excitation

The reactor geometry is simply a cylindrical glass tube with a 4-turn coil wrapped around it. Gas flows in from the bottom and exits out of the top. The gas is heated through elastic and inelastic collisions. A fixed power of 700 W is applied to the coil.

Figure 2: Schematic of the ICP reactor. Flow enters from the base and leaves out the top.

Strategy to setup a model with inlets and outlets

In this model, an inlet or argon is added in the system. There is also an outlet responsible to fix the system pressure. After setting the and features in the interface the following reasoning should be used to set up the and subfeatures of the species in the interface:

•

If a species is set it does not need any and subfeature since its mass fraction is computed so that the mass in the system is conserved. There are no dependent variables being solved for the mass constraint species and as such does not need the boundary conditions set by the subfeatures.

•

If the feed into the system contains another species (the present model does not have) one of these species should be set as mass constraint and the other should use an subfeature that specifies the appropriate proportion. Consider as an example of a feed of Ar/O2 at 60/40 mole fraction. If Ar is set as mass constrain than O2 should have an subfeature with a mole fraction of 0.4.

•

If the species in created inside the reactor (as are Ar+ and Ars) no subfeature should be added.

•

If the species flows out of the system (as are Ar+ and Ars) an subfeature should be added.

Results

Results of the different quantities computed in this model are presented below. The electron density reaches high values while the electron temperature is kept very low meaning that electron losses by transport are small. The gas temperature in the system peaks to 1300 K showing that there is a considerable amount of energy transferred from the electrons to the heavy species that is not lost in transport. More that 90% of the heating of the gas is caused by electron impact elastic collisions defined in reaction 1.

Figure 3: Electron density inside the column.

Figure 4: Electron temperature inside the plasma source.

Figure 5: Plasma potential inside the plasma source.

Figure 6: Velocity field inside the plasma source.

Figure 7: Temperature inside the plasma source.

Figure 8: Mass fraction of ground state argon.

Figure 9: Mass fraction of electronically excited argon atoms.

References

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

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

Plasma_Module/Inductively_Coupled_Plasmas/thermal_plasma

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

Click the

button in the toolbar.

In the dialog, in the tree, select the checkbox for the node > .

In the tree, select the checkbox for the node > .

Click .

Import 1 (imp1)

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file thermal_plasma.mphbin.

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
p0	1[torr]	133.32 Pa	Initial and outlet pressure

Definitions

Plasma

In the toolbar, click

Select Domain 1 only.

In the window for , type Plasma in the text field.

Walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 9, 10, 34, and 35 only.

In the text field, type Walls.

Outlet

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 4 only.

In the text field, type Outlet.

Coil Walls

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundaries 12, 13, 15–17, 19, 21, 22, 24–26, and 28–32 only.

In the text field, type Coil Walls.

Start by importing the cross sections for argon and by activating the convection and thermodynamic property evaluation.

Plasma (plas)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Cross Section Import 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file Ar_xsecs.txt.

Click

In the window, click .

In the window for , locate the section.

Select the checkbox.

Reaction 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Ars+Ars=>e+Ar+Ar+.

Locate the section. In the kf text field, type 3.734E8.

Reaction 2

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Ars+Ar=>Ar+Ar.

Locate the section. In the kf text field, type 1807.

Ground 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Surface Reaction 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type Ar+=>Ar.

Surface Reaction 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type Ars=>Ar.

Wall 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

Species: Ar

In the window, click .

In the window for , locate the section.

Select the checkbox.

Locate the section. From the list, choose .

Species: Ars

In the window, click .

In the window for , locate the section.

In the x0 text field, type 1E-4.

From the list, choose .

Click to expand the section. In the Δh text field, type 11.5.

The thermodynamic properties for the electronically excited argon atoms can be the same as for the ground state species plus the threshold energy for the electron impact reaction. In this case this corresponds to an energy of 11.5eV. This is added in the text field .

Species: Ar+

In the window, click .

In the window for , locate the section.

Select the checkbox.

Locate the section. From the list, choose .

Locate the section. In the Δh text field, type 15.8.

The thermodynamic properties for the argon ions can be the same as for the ground state species plus the threshold energy for ionization. In this case this corresponds to an energy of 15.8eV. This is added in the text field .

You can set the gas temperature and pressure in the plasma model to the computed gas pressure and temperature from other physics interfaces. The velocity field is also set to the velocity field computed from the interface.

Plasma Model 1

In the window, click .

In the window for , locate the section.

In the μeNn text field, type 4E24.

Initial Values 1

In the window, click .

In the window for , locate the section.

In the ne,0 text field, type 1E15.

In the ε0 text field, type 3.

Magnetic Fields (mf)

In the window, under click .

Domain Coil 1

In the toolbar, click

and choose .

Select Domains 4–7 only.

In the window for , locate the section.

From the list, choose .

Select the checkbox.

In the Pcoil text field, type 700[W].

Laminar Flow (spf)

In the window, under click .

In the window for , click to expand the section.

Find the subsection. Clear the checkbox.

Define the pressure reference level to be 1 torr.

Locate the section. In the pref text field, type p0.

Locate the section. Click

Select Domain 1 only.

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

Inlet 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the Qsccm text field, type 100*tanh(1E5*t[1/s]).

Outlet 1

In the toolbar, click

and choose .

Select Boundary 4 only.

Heat Transfer in Fluids (ht)

In the window, under click .

In the window for , locate the section.

Click

Select Domain 1 only.

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the T text field, type 300.

Temperature 1

In the toolbar, click

and choose .

Select Boundaries 2, 9, 10, 34, and 35 only.

In the window for , locate the section.

In the T0 text field, type 300.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 4 only.

Materials

Dielectric

In the window, under right-click and choose .

Select Domain 2 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	4.5	1	Basic

In the text field, type Dielectric.

Air

Right-click and choose .

Select Domain 3 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	0	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type Air.

Copper coil

Right-click and choose .

Select Domains 4–7 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	6E7	S/m	Basic
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic