COMSOL 6.4 - Microwave Cavity Plasma Reactor

This tutorial shows how to prepare a model of a hydrogen plasma sustained in a microwave cavity at moderate pressures. The model is based on the work in Ref. 1 and solves the plasma transport equations fully coupled with Maxwell’s equations, fluid flow, and heat transfer. A microwave cylindrical chamber contains a bell jar were a hydrogen plasma is created. The reactor is carefully designed so that the electric field has its maximum intensity above a substrate and much lower intensity at the bell jar boundary. This is important in order for the plasma to keep its maximum density in the substrate region even at high power.

The model requires the Plasma Module and the RF Module.

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.

In a microwave reactor, the high-frequency electric field is computed in the frequency domain through the equation

The electromagnetic wave “sees” a plasma defined by the plasma conductivity in the cold plasma approximation that is set in the multiphysics feature:

where ne is the electron density, q is the electron charge, me is the electron mass, νe is the collision frequency, and ω is the angular frequency. The Joule heating term that is responsible to heat the electrons is set in the multiphysics feature.

The fluid flow and heat transfer in the fluid are also solved for. The gas temperature plays an important role in the reactor operation since it can significantly influence the reduced electric field E/N.

Plasma Chemistry

The plasma chemistry is based on the simplified chemistry presented in table 4 of Ref. 1 and is presented in Table 1. The electron impact reactions are from Ref. 2 and retrieved from Ref. 3, and the rates are from Ref. 4 and Ref. 5. The model includes eight species: electrons, H2, H, H2+, H+, H3+, and two excited states of hydrogen corresponding to the levels n = 2 and n = 3 that are represented by Hn2 and Hn3.

Table 1: Modeled collisions and reactions.
Reaction	Formula	Type
1	e+H2=>e+H2	Elastic	-
2	e+H2=>e+H2	Vibrational excitation	0.516-1.5
3	e+H2=>e+2H	Dissociation	7.93-11.72 and17.22-17.53
4	e+H2=>e+H2	Excitation	12.4-14.6
5	e+H2=>e+H+Hn2	Dissociative excitation	14.68
6	e+H2=>e+H+Hn3	Dissociative excitation	16.57
7	e+H2=>2e+H2+	Ionization	15.4
8	e+H2=>2e+H+H+	Ionization	19
9	e+H=>e+H	Elastic	-
10	e+H=>e+Hn2	Excitation	10.2043
11	e+H=>e+Hn3	Charge transfer	12.1
12	e+H=>e+H	Excitation	12.755 and 13.0615
13	e+H=>2e+H+	Ionization	13.6057
14	e+H3+=>3H	Recombination	0
15	e+H2+=>H+Hn2	Recombination	0.01
16	e+H2+=>H+Hn3	Recombination	0.01
17	e+H+=>Hn2	Recombination	0
18	e+H+=>Hn3	Recombination	0
19	Hn2+H2=>H3++e	Ionization	-
20	Hn3+H2=>H3++e	Ionization	-
21	H2+H2+=>H3++H	Ionization	-
22	H2+H2=>2H+H2	Dissociation	-
23	2H+H2=>H2+H2	Association	-
24	H2+H=>3H	Dissociation	-
26	3H=>H2+H	Association	-

In addition to the volume reactions, the surface reaction listed in Table 2.

Table 2: surface reactions.
Reaction	Formula	Sticking coefficient
1	H=>0.5H2	0.02
2	Hn2=>H	1
3	Hn3=>H	1
4	H+=>H	1
5	H2+=>H2	1
6	H3+=>H2+H	1

Results and Discussion

The figures in this section present model results for a hydrogen plasma sustained at 25 kPa with 5 kW. In general, the results agree well with results presented in figure 21 and figure 22 of Ref. 1. A detailed comparison is not attempted because many aspects of the models are different. The plasma chemistry is not exactly the same. Even if the reactions are mostly the same, the rates and cross sections were obtained independently. The reactor dimensions and configurations are only approximate. Additionally, there are differences in the model equations.

Using a global model fully coupled with a Boltzmann equation in the two-term approximation, the study performed in the model Hydrogen Global Model Coupled with the Two-Term Boltzmann Equation — also in the Plasma Module Application Library — showed that for main quantities like electron density, H-atom density, and gas temperature there is not a significant difference when using a computed or a Maxwellian EEDF. This justified doing the space-dependent simulations using a Maxwellian EEDF. However, in this first attempt the results were wrong in many ways. An important aspect not captured was the plasma localization at the subtract. With a Maxwellian EEDF it is very easy to ionize in regions of low electron energy and the maximum of the plasma easily moves somewhere else. Further investigations with the global model showed significant differences in results for low mean electron energies. The EEDF was than adjusted to use a generalized distribution with a power law of 1.2.

The reactor dimensions were adjusted to have high electric field intensity at the substrate and to have a minimum in reflected power without plasma. This part of the investigation is done with the Electromagnetic Waves, Frequency Domain interface only and it is not shown here. With the plasma present, the electromagnetic characteristics of the reactor change considerably. This is possible to observe in Figure 5 and Figure 6, which show the S-parameter and Smith plot when increasing the power from 500 to 5000 W. Having a high-dense plasma in a considerable part of the reactor strongly increases the reflected power.

Figure 1: Electron density for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 2: Electron temperature for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 3: RF electric field norm for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 4: RF power absorbed by electrons for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 5: S-parameter as a function of the input power.

Figure 6: Smith plot as a function of the input power. For higher powers the curve moves away from the center.

Figure 7: Fluid velocity for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 8: Gas temperature for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

Figure 9: Number density of atomic hydrogen for a hydrogen plasma operating at 25 kPa with 5 kW of input power.

References

1. K. Hassouni, F. Silva, and A. Gicquel, “Modelling of diamond deposition microwave cavity generated plasmas,” J. Phys. D: Appl. Phys., vol. 43, p. 153001, 2010.

2. L. Marques, J. Jolly, and L.L. Alves, “Capacitively Coupled Radio-Frequency Hydrogen Discharges: The Role of Kinetics,” J. Appl. Phys., vol. 102, p. 063305, 2007.

3. IST-Lisbon database, www.lxcat.net, retrieved 2023.

4. M. Capitelli, C.M. Ferreira, B.F. Gordiets and A.I. Osipov, Plasma Kinetics in Atmospheric Gases, Springer, 2000.

5. R.K. Janev; W.D. Langer; K. Evans, Jr.; and D.E. Post, Jr., Elementary Processes in Hydrogen-Helium Plasmas, Springer-Verlag, 1987.

Plasma_Module/Wave-Heated_Discharges/microwave_cavity_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

Import a file with parameters to be used in the model. Most of them are used to create the reactor geometry.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

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

Create a profile to be used as initial conditions for the gas temperature.

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
r0	0[cm]	m	Maximum r-coordinate
z0	14[cm]	m	Maximum z-coordinate
prf	exp(-((r0-r)/4[cm])^2 -((z0-z)/2[cm])^2)		Initial temperature profile
Tinit	2500[K]*prf+1200[K]	K	Initial temperature