PDF
Model of an Argon/Oxygen Inductively Coupled Plasma Reactor
Introduction
This tutorial models an argon/oxygen inductively coupled plasma reactor. The main goal is to show how to prepare a model with a mixture of different elements (in this case Ar and O2) in which one of the species can dissociate by electron impact (O2 dissociates into O) and where negative ions exist (the dissociative electron attachment of O2 creates O-).
A simplified plasma chemistry is used to discuss the main aspects of such discharges. It is important to keep in mind that a benchmark is not attempted and the idea is to provide a base case that can be used to develop more complex chemistries. In fact, it might be necessary to modify the data used and add more reactions to achieve experimental verification.
The plasma model is solved self-consistently with the Laminar Flow and Heat Transfer in Fluids interfaces. A Frequency-Stationary study step is used to do fast parameterization in the input power and oxygen mole fraction.
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.
For a nonmagnetized, nonpolarized plasma, the induction currents are computed in the frequency domain using the following equation:
The electromagnetic wave “sees” a plasma defined by the plasma conductivity in the cold plasma approximation that is set in the Plasma Conductivity Coupling 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 Electron Heat Source multiphysics feature.
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 following boundary condition for the electron flux:
and the electron energy flux:
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:
The walls of the reactor are grounded.
Plasma Chemistry
Negative ions are created in certain molecular gaseous discharges (like chlorine, oxygen, hydrogen, fluorocarbons, and so on) and these discharges tend to have complex plasma chemistries with many ions, dissociative products, and excited states. Here a simple plasma chemistry is used and no benchmark is attempted. In fact, it might be necessary to modify the data used and add more reactions to achieve experimental verification. Nevertheless, this plasma chemistry allows to show the main aspects of an electronegative discharge. The plasma chemistry for oxygen is based on the one presented in Ref. 1 (from the section “A Data Set for Oxygen,” page 270) but the electron impact reactions are mostly retrieved from Ref. 2 in the form of electron impact cross sections. A good discussion of the chemistry of an oxygen/argon plasma at low pressures can be found in Ref. 3.
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 presented in Table 1 (electron impact cross sections are obtained form Ref. 4).
Table 1: Argon reactions
reaction
formula
type
1
e+Ar=>e+Ar
Elastic
-
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
-
Oxygen has a much more rich reaction set that includes vibrational and rotational excitations, excitation of several electronic excited states, electron impact dissociation, dissociative attachment, and many others. Electron impact reactions with O2 are from Ref. 5 and electron impact reaction with O are from Ref. 6 except for e+O-=>O+e+e, which is from Ref. 1. The electron impact reactions used in this model are presented in Table 2. The following simplifications were made: 3-body attachment is not included, rotational and vibrational states are not treated explicitly but energy losses are considered, the dissociative excitation reaction at 14.7 eV is not included, polar dissociation is not included, reverse reaction by detailed balance are not included for O2 and O excited states. The only oxygen excited states that are solved explicitly are the singlet delta metastable of molecular oxygen O2(a1Δg) and the metastable state O(1D).
Table 2: oxygen Electron impact reactions
reaction
formula
type
1
e+O2=>O+O-
Dissociative attachment
-
2
e+O2=>e+O2
Elastic
-
3
e+O2=>e+O2
Rotational excitation
0.02
4
e+O2=>e+O2
Vibrational excitation
0.19
5
e+O2=>e+O2
Vibrational excitation
0.19
6
e+O2=>e+O2
Vibrational excitation
0.38
7
e+O2=>e+O2
Vibrational excitation
0.38
8
e+O2=>e+O2
Vibrational excitation
0.57
9
e+O2=>e+O2
Vibrational excitation
0.75
10
e+O2=>e+O2(a1Δg)
Excitation
0.977
11
e+O2=>e+O2
Excitation
1.627
12
e+O2=>e+O2
Excitation
4.5
13
e+O2=>e+O+O
Dissociative excitation
6.0
14
e+O2=>e+O+O(1D)
Excitation
8.4
15
e+O2=>e+O2
Excitation
9.97
16
e+O2=>e+O2+
Ionization
12.06
17
e+O=>e+O
Elastic
-
18
e+O=>e+O(1D)
Excitation
1.968
19
e+O=>e+O
Excitation
4.192
20
e+O=>e+c
Ionization
13.618
21
e+O-=>O+e+e
Electron impact detachment
12.00
On Table 3 are presented heavy species reaction involving ions. For reaction 6 it is used the same rate constant as for reaction 2.
Table 3: Heavy species reactions involving ions
Reaction
Formula
Type
1
O++O2=>O+O2+
Charge transfer
2
O++O-=>O+O
Mutual recombination
3
O-+O2+=>3O
Mutual recombination
4
O-+O2+=>O+O2
Mutual recombination
5
O-+O2=>O2+e
Detachment
6
O-+Ar+=>O+Ar
Mutual recombination
In addition to volumetric reactions, the following surface reactions are implemented.
Table 4: Surface reactions
reaction
formula
Sticking Coefficient
Secondary emission coefficient
Mean energy of secondary electrons (V)
1
Ars=>Ar
1
0
0
2
Ar+=>Ar
1
0.07
5.8
3
O=>0.5O2
0.2
0
0
4
O2+=>O2
1
0.05
4
5
O-=>O
1
0
0
6
O2(a1Δg)=>O
1
0
0
7
O(1D)=>0.5O2
0.2
0
0
8
O+=>O
1
0.05
4
It is by using surface reactions that boundary conditions for heavy species are introduced in the model. If no surface reactions that leads to the lost of a given species at a surface are introduced in the model, that species will not have losses by transport. This can lead to the unbounded growth of a given species and a steady state solution might not be possible.
Atomic recombination (reaction 3 in Table 4) at a surface is an important aspect of plasma discharges with molecular species since it influences the dissociation degree in the discharge. The sticking coefficient for atomic recombination is a function of the surface type and temperature.
Electronegative plasmas
Electronegative plasmas are plasmas that contain negative ions. Negative ions are mainly created by electron dissociative attachment (for example, e+O2=>O+O-). This reaction tends to be very effective at low electron energies and can reduce the electrons in a discharge to a point that an ion–ion discharge is obtained. The transport and volume creation/destruction mechanisms tend to be more complex than electropositive plasmas in many aspects. Here only a few are mentioned with emphasis on the numerical difficulties that they introduce. More information can be found in Ref. 1 section 10.3 and references therein.
In electronegative discharges negative ions are well confined by the ambipolar electric field and losses by transport are very small. This means that to achieve a steady state volume losses need to be included for negative ions. The mechanisms by which negative ions are lost depend on the gas mixture and pressure and they are: mutual recombination with positive ions (for example, O-+O+=>O+O or O-+Ar+=>O+Ar), detachment in collisions with excited or neutral atoms or molecules (for example, O-+O=>O2+e or O-+O2=>O+O2+e), and electron-impact detachment (for example, e+O-=>O+2e).
In electronegative discharges it is often possible to identify two spatial regions using the electronegativity (ratio of the negative ion density to the electron density): (i) one in the core of the discharge (the electronegative core) with high electronegativity where the dominant charge species are positive and negative ions; (ii) and the other close to the boundaries (electropositive edges) where the dominant charged species are electron and positive ions. In the transition between these two regions the negative ion density drops abruptly causing a chock-like phenomena. This transition needs to be well resolved spatially. If not, oscillations can be seen in the negative ion density and the model might not converge. Some strategies to deal with this are:
•
Increase the negative ion temperature of about 0.3 eV. An higher ion temperature makes the transport numerical easier. The ion temperature is defined in the section Mobility and Diffusivity Expressions in the species Settings. By default the ion temperature is the gas temperature.
•
Enable Isotropic diffusion for ions in the Inconsistent Stabilization section (the stabilization sections are visible when Stabilization is selected in Show More Options). This option adds artificial diffusion to all ions and helps smoothing the sharp transition of the negative ion density between the electropositive edge and the electronegative core, and also increase the density of the negative ions in the electropositive edge effectively increasing its losses by transport. This option should be used very carefully since completely wrong results can be obtained if too much diffusion is used (the tuning parameter for ions should not be larger than 0.1). A useful strategy is to start with a large Tuning parameter for ions (for example, 0.5) and ramp it down using an Auxiliary sweep.
Inlet and outflow
When solving for plasmas with chemistries that contain more than one element (for example, Ar and O2) with a stationary solver the mass fraction of each element is not conserved if no constraint is used. This problem is similar in nature to the one found when solving for Navier–Stokes equations in steady state without fixing the pressure somewhere. The Inlet boundary condition fixes the mass fraction or mole fraction of specified species and it is used in this model as a real inlet for O2 and as a strategy to fix its mass fraction. An important aspect to remember is that no Inlet is used for Ar since this species is being computed From mass constraint. It is assumed that ions are neutralized at the outlet boundary and as such no Outflow boundary condition is applied. For the species to which an Outflow boundary condition is applied no surface reaction is applied at the same boundary to represent a species that flows out of the system without interacting with a surface.
Results and Discussion
Figure 1 to Figure 4 show spatial distributions of the electron density, electron temperature, negative ion density, and power absorbed by the electrons at 1000 W. The power absorbed by the electrons and electron temperature profiles are typical of an ICP discharge. Most of the power deposition occurs below the coil and is shielded by the high dense plasma, and the electron temperature is relatively flat in the plasma region. The electron and negative ion density have a typical distribution of electronegative plasmas with an electronegative core (where the electrons have a flat profile) and electropositive edges (where the negative ion density drops fast toward the surface). These regions can be better observed in Figure 5 that shows the distribution of the charged species along the axis-of-symmetry.
Figure 6 presents the space average number density of charged species as a function of power. At 100 W the negative ions are the dominant negatively charged species. As power increases the discharge electronegativity decreases; at 1000 W, the averaged electronegativity is slightly smaller than one, with the electrons now being the dominant negative charged species. However, the electronegative core is always present, as shown in Figure 5.
Figure 7 presents the space average number density of charged species as a function of the mole fraction of oxygen. With low oxygen content the discharge has low electronegativity and O+ is the dominant positive ion. Increasing the oxygen content leads to higher electronegativity and a progressive change of the dominant positive ion to be O2+.
Figure 1: Electron density spatial distribution at1000 W.
Figure 2: Electron temperature spatial distribution at1000 W.
Figure 3: Negative ions density spatial distribution at1000 W.
Figure 4: Power absorbed by the electrons at 1000 W.
Figure 5: Charged species distribution along the axis-of-symmetry at 1000 W.
Figure 6: Spatial averaged number density of the charged species as a function of power.
Figure 7: Spatial averaged number density of the charged species as a function of oxygen mole fraction.
References
1. M.A. Lieberman and A.J. Lichtenberg, Principles of Plasma Discharges and Materials Processing, John Wiley & Sons, 2005.
2.  www.lxcat.net
3. J.T. Gudmundsson and E.G. Thorsteinsson, “Oxygen discharges diluted with argon: dissociation process,” Plasma Sources Sci. Technol., vol. 16, pp. 399–412, 2007.
4. Phelps database, www.lxcat.net, retrieved 2017.
5. Phelps database, www.lxcat.net, retrieved 2022.
6. Morgan database, www.lxcat.net, retrieved 2022.
Application Library path: Plasma_Module/Inductively_Coupled_Plasmas/icp_argon_oxygen
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Plasma > Nonisothermal Plasma Flow > Inductively Coupled Plasma.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select Preset Studies for Selected Multiphysics > Frequency–Stationary.
6
Click  Done.
Geometry 1
Create the geometry for an ICP reactor.
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose cm.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 30.
4
In the Height text field, type 30.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 30.
4
In the Height text field, type 3.
5
Locate the Position section. In the z text field, type 20.
Rectangle 3 (r3)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Position section.
3
In the r text field, type 3.
4
In the z text field, type 23.
Array 1 (arr1)
1
In the Geometry toolbar, click  Transforms and choose Array.
2
Select the object r3 only.
3
In the Settings window for Array, locate the Size section.
4
In the r size text field, type 4.
5
Locate the Displacement section. In the r text field, type 4.5.
Rectangle 4 (r4)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 15.
4
In the Height text field, type 5.
Rectangle 5 (r5)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 6.
4
In the Height text field, type 5.
5
Locate the Position section. In the r text field, type 15.
Rectangle 6 (r6)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 10.
4
In the Height text field, type 6.
5
Locate the Position section. In the r text field, type 20.
6
In the z text field, type 14.
Line Segment 1 (ls1)
1
In the Geometry toolbar, click  More Primitives and choose Line Segment.
2
In the Settings window for Line Segment, locate the Starting Point section.
3
From the Specify list, choose Coordinates.
4
In the r text field, type 20.
5
In the z text field, type 17.
6
Locate the Endpoint section. From the Specify list, choose Coordinates.
7
In the r text field, type 20.
8
In the z text field, type 18.
9
Click  Build All Objects.
Add some parameters to be used in the model.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
Pw
100[W]
100 W
 
Qf
250
250
 
xO2
0.9
0.9
 
xO
1e-4
1E-4
 
p0
0.02[torr]
2.6664 Pa
 
Create explicit selections to be used later in the model when defining selections of different features.
Definitions
Walls
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type Walls in the Label text field.
3
Select Domain 2 only.
4
Locate the Output Entities section. From the Output entities list, choose Adjacent boundaries.
Coils
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type Coils in the Label text field.
3
Select Domains 5–7 and 9 only.
Coil Boundaries
1
Right-click Coils and choose Duplicate.
2
In the Settings window for Explicit, type Coil Boundaries in the Label text field.
3
Locate the Output Entities section. From the Output entities list, choose Adjacent boundaries.
Average 1 (aveop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
Add an average operator to evaluate space averages of several quantities in the results section.
2
Select Domain 2 only.
Assign material properties to the different regions of the modeling domain. The domain where the plasma exists is automatically assigned a relative permittivity of 1.
Definitions
In the Model Builder window, collapse the Component 1 (comp1) > Definitions node.
Geometry 1
In the Model Builder window, collapse the Component 1 (comp1) > Geometry 1 node.
Add Material
1
In the Materials toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in > Air.
4
Click the Add to Component button in the window toolbar.
5
In the tree, select Built-in > Glass (quartz).
6
Click the Add to Component button in the window toolbar.
7
In the tree, select Built-in > Copper.
8
Click the Add to Component button in the window toolbar.
9
In the Materials toolbar, click  Add Material to close the Add Material window.
Materials
Air (mat1)
1
In the Settings window for Material, locate the Geometric Entity Selection section.
2
Click  Clear Selection.
3
Select Domain 4 only.
Glass (quartz) (mat2)
1
In the Model Builder window, click Glass (quartz) (mat2).
2
Select Domain 3 only.
Copper (mat3)
1
In the Model Builder window, click Copper (mat3).
2
In the Settings window for Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Coils.
The Plasma Chemistry Import feature
The next steps have instructions to use the Plasma Chemistry Import feature to import a file that automatically creates the argon-oxygen plasma chemistry.
The following is set or created automatically:
a
Species properties using Preset species data
b
Reaction group features for argon and oxygen
c
Surface reaction group features
The documentation accompanying the Plasma Chemistry Import feature contains more information about the file structure and what can be set automatically.
Plasma (plas)
Plasma Chemistry Import 1
1
In the Physics toolbar, click  Global and choose Plasma Chemistry Import.
2
In the Settings window for Plasma Chemistry Import, locate the Plasma Chemistry Import section.
3
Click  Browse.
4
Browse to the model’s Application Libraries folder and double-click the file Ar_O2_plasma_chemistry.txt.
5
Click  Import.
6
In the Model Builder window, click Plasma (plas).
7
In the Settings window for Plasma, locate the Domain Selection section.
8
Click  Clear Selection.
9
Select Domain 2 only.
10
Click the  Show More Options button in the Model Builder toolbar.
11
In the Show More Options dialog, select Physics > Stabilization in the tree.
12
In the tree, select the checkbox for the node Physics > Stabilization.
13
Click OK.
This model needs stabilization because the density of the negative ions can drop sharply when approaching the reactor edges.
14
In the Settings window for Plasma, click to expand the Inconsistent Stabilization section.
15
Select the Isotropic diffusion for ions checkbox.
16
In the δid,i text field, type 0.1.
17
Locate the Transport Settings section. Select the Mixture diffusion correction checkbox.
The surface reactions used in the model were created automatically but it is still necessary to specify the boundaries where they take place.
On the outlet boundary, it is assumed that ions are neutralized but nothing happens to neutrals. For this reason, neutral species do not have a Surface Reaction defined at the outlet.
Surface Reactions - Ions
1
In the Model Builder window, click Surface Reactions - Ions.
2
In the Settings window for Surface Reaction Group, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
Surface Reactions - Neutrals
1
In the Model Builder window, click Surface Reactions - Neutrals.
2
In the Settings window for Surface Reaction Group, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
4
Select Boundaries 3, 4, 6, 27, 33–36, 38, and 40 only.
In the following, the initial mole fraction for O2 and the initial number density for ions are specified. The mass fraction of Ar is found from a mass constraint and the initial density of Ar+ is found by requiring electroneutrality.
The mole fraction of O2 is set at the inlet. Since the mass fraction of Ar is found from mass constraint nothing needs to be done to it.
Only neutral species have an Outflow boundary condition. It is assumed that ions are neutralized at the outlet.
Species: O2
1
In the Model Builder window, expand the Component 1 (comp1) > Plasma (plas) > Group - Species node, then click Species: O2.
2
In the Settings window for Species, locate the General Parameters section.
3
In the x0 text field, type xO2.
Species: O
1
In the Model Builder window, click Species: O.
2
In the Settings window for Species, locate the General Parameters section.
3
In the x0 text field, type xO.
Species: O-
1
In the Model Builder window, click Species: O-.
2
In the Settings window for Species, locate the General Parameters section.
3
In the n0 text field, type 1E10[1/m^3].
Species: O2+
1
In the Model Builder window, click Species: O2+.
2
In the Settings window for Species, locate the Species Formula section.
3
Select the Initial value from electroneutrality constraint checkbox.
Species: O+
1
In the Model Builder window, click Species: O+.
2
In the Settings window for Species, locate the General Parameters section.
3
In the n0 text field, type 1E10[1/m^3].
Species: Ar
1
In the Model Builder window, click Species: Ar.
2
In the Settings window for Species, locate the Species Formula section.
3
Select the From mass constraint checkbox.
Inflow 1
1
In the Physics toolbar, click  Boundaries and choose Inflow.
2
Select Boundary 35 only.
3
In the Settings window for Inflow, locate the Inflow section.
4
Click  Add.
5
In the table, enter the following settings:
Species names
Mole fraction (1)
O2
xO2
Outflow 1
1
In the Physics toolbar, click  Boundaries and choose Outflow.
2
Select Boundary 39 only.
Species: Ar+
1
In the Model Builder window, under Component 1 (comp1) > Plasma (plas) > Group - Species click Species: Ar+.
2
In the Settings window for Species, locate the General Parameters section.
3
In the n0 text field, type 1E10[1/m^3].
Set the Additional enthalpy contribution field for species with internal energy so that the reaction enthalpy is correctly computed.
Species: O2a1Dg
1
In the Model Builder window, click Species: O2a1Dg.
2
In the Settings window for Species, click to expand the Species Thermodynamic Parameters section.
3
In the Δh text field, type 9.77e-1.
Species: O1D
1
In the Model Builder window, click Species: O1D.
2
In the Settings window for Species, locate the Species Thermodynamic Parameters section.
3
In the Δh text field, type 1.968.
Species: O2+
1
In the Model Builder window, click Species: O2+.
2
In the Settings window for Species, locate the Species Thermodynamic Parameters section.
3
In the Δh text field, type 12.06.
Species: O+
1
In the Model Builder window, click Species: O+.
2
In the Settings window for Species, locate the Species Thermodynamic Parameters section.
3
In the Δh text field, type 13.618.
Species: Ars
1
In the Model Builder window, click Species: Ars.
2
In the Settings window for Species, locate the Species Thermodynamic Parameters section.
3
In the Δh text field, type 11.5.
Species: Ar+
1
In the Model Builder window, click Species: Ar+.
2
In the Settings window for Species, locate the Species Thermodynamic Parameters section.
3
In the Δh text field, type 15.8.
Group - Species
In the Model Builder window, collapse the Component 1 (comp1) > Plasma (plas) > Group - Species node.
Set the plasma model to use the temperature, fluid velocity, and pressure computed by the Heat Transfer in Fluids and Laminar Flow interfaces.
Plasma Model 1
1
In the Model Builder window, click Plasma Model 1.
2
In the Settings window for Plasma Model, locate the Electron Density and Energy section.
3
From the Electron transport properties list, choose From electron impact reactions.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the ne,0 text field, type 1E15[1/m^3].
4
In the ε0 text field, type 2[V].
Ground 1
1
In the Physics toolbar, click  Boundaries and choose Ground.
2
In the Settings window for Ground, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
4
Select Boundaries 4, 6, 27, 33–36, and 38–40 only.
The Wall node sets boundary conditions for the electron transport equations.
Wall 1
1
In the Physics toolbar, click  Boundaries and choose Wall.
2
In the Settings window for Wall, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
Add the coil that is responsible for exciting the plasma.
Plasma (plas)
In the Model Builder window, collapse the Component 1 (comp1) > Plasma (plas) node.
Magnetic Fields (mf)
1
In the Model Builder window, under Component 1 (comp1) click Magnetic Fields (mf).
2
Select Domains 2–7 and 9 only.
Domain Coil 1
1
In the Physics toolbar, click  Domains and choose Domain Coil.
2
In the Settings window for Domain Coil, locate the Domain Selection section.
3
From the Selection list, choose Coils.
4
Locate the Coil section. Select the Coil group checkbox.
5
From the Coil excitation list, choose Power.
6
In the Pcoil text field, type Pw.
Set the thermodynamic properties and a heat source to use the values computed in the Plasma interface.
Set the wall temperature to 300 K.
Heat Transfer in Fluids (ht)
1
In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Fluids (ht).
2
In the Settings window for Heat Transfer in Fluids, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domain 2 only.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
In the Settings window for Temperature, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
4
Locate the Temperature section. In the T0 text field, type 300[K].
Set the fluid properties to use values computed in the Plasma interface. And set the flow inlet and outlet in the system.
Laminar Flow (spf)
1
In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
2
In the Settings window for Laminar Flow, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domain 2 only.
5
Locate the Physical Model section. In the pref text field, type p0.
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
Select Boundary 35 only.
3
In the Settings window for Inlet, locate the Boundary Condition section.
4
From the list, choose Mass flow.
5
Locate the Mass Flow section. From the Mass flow type list, choose Standard flow rate (SCCM).
6
In the Qsccm text field, type Qf.
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Select Boundary 39 only.
Mesh 1
Size 1
1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 2 only.
5
Locate the Element Size section. From the Calibrate for list, choose Plasma.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Finer.
Size 2
1
In the Model Builder window, right-click Mesh 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Walls.
5
Select Boundaries 4, 6, 27, 33–36, and 38–40 only.
6
Locate the Element Size section. From the Calibrate for list, choose Plasma.
7
From the Predefined list, choose Finer.
Mapped 1
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Coils.
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Boundary Selection section.
3
From the Selection list, choose Coil Boundaries.
4
Locate the Distribution section. From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type 25.
6
In the Element ratio text field, type 20.
7
Select the Symmetric distribution checkbox.
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Boundary Layers 1
1
In the Mesh toolbar, click  Boundary Layers.
2
In the Settings window for Boundary Layers, click to expand the Corner Settings section.
3
From the Handling of sharp corners list, choose No special handling.
4
Click to expand the Transition section. Clear the Smooth transition to interior mesh checkbox.
5
Locate the Domain Selection section. From the Geometric entity level list, choose Domain.
6
Select Domain 2 only.
Boundary Layer Properties
1
In the Model Builder window, click Boundary Layer Properties.
2
In the Settings window for Boundary Layer Properties, locate the Boundary Selection section.
3
From the Selection list, choose Walls.
4
Select Boundaries 4, 6, 27, 33, 34, 36, and 38–40 only.
5
Click  Build All.
The first study solves for 100 W and xO2=0.9. The solution of this study will be used in the subsequent study.
Base Case
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Base Case in the Label text field.
3
In the Study toolbar, click  Get Initial Value.
Organize the results from this study with a Group node.
Results
Electric Potential (plas), Electron Density (plas), Electron Temperature (plas), Magnetic Flux Density (mf), Magnetic Flux Density, Revolved Geometry (mf), Pressure (spf), Temperature (ht), Velocity (spf), Velocity, 3D (spf)
Right-click and choose Group.
Base Case
In the Settings window for Group, type Base Case in the Label text field.
Base Case
Step 1: Frequency–Stationary
1
In the Model Builder window, expand the Base Case > Solver Configurations node, then click Base Case > Step 1: Frequency–Stationary.
2
In the Settings window for Frequency–Stationary, locate the Study Settings section.
3
In the Frequency text field, type 13.56[MHz].
Solution 1 (sol1)
1
In the Model Builder window, expand the Base Case > Solver Configurations > Solution 1 (sol1) > Stationary Solver 1 node, then click Fully Coupled 1.
2
In the Settings window for Fully Coupled, click to expand the Results While Solving section.
3
Select the Plot checkbox.
4
Click  Run.
Results
Electron Density (plas)
Add a study to do a power sweep. This study will use the solution of the previous study as initial conditions to save computation time.
The Initial damping factor is set to 1 because a solution is used as initial condition and because the next parameter will use the previous solution.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select Preset Studies for Selected Multiphysics > Frequency–Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Frequency–Stationary
1
In the Settings window for Frequency–Stationary, locate the Study Settings section.
2
In the Frequency text field, type 13.56[MHz].
3
Click to expand the Values of Dependent Variables section. Find the Initial values of variables solved for subsection. From the Settings list, choose User controlled.
4
From the Method list, choose Solution.
5
From the Study list, choose Base Case, Frequency–Stationary.
6
Click to expand the Study Extensions section. Select the Auxiliary sweep checkbox.
7
Click  Add.
8
In the table, click to select the cell at row number 1 and column number 2.
9
Click  Range.
10
In the Range dialog, type 100 in the Start text field.
11
In the Step text field, type 50.
12
In the Stop text field, type 1000.
13
Click Replace.
14
In the Model Builder window, click Study 2.
15
In the Settings window for Study, type Power Sweep in the Label text field.
16
In the Study toolbar, click  Get Initial Value.
Results
Electric Potential (plas) 1, Electron Density (plas) 1, Electron Temperature (plas) 1, Magnetic Flux Density (mf) 1, Magnetic Flux Density, Revolved Geometry (mf) 1, Pressure (spf) 1, Temperature (ht) 1, Velocity (spf) 1, Velocity, 3D (spf) 1
Right-click and choose Group.
Results
Power Sweep
1
In the Model Builder window, expand the Power Sweep > Solver Configurations node, then click Results > Group 2.
2
In the Settings window for Group, type Power Sweep in the Label text field.
Power Sweep
Solution 2 (sol2)
1
In the Model Builder window, expand the Power Sweep > Solver Configurations > Solution 2 (sol2) > Stationary Solver 1 node, then click Fully Coupled 1.
2
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
3
In the Initial damping factor text field, type 1.
4
Locate the Results While Solving section. Select the Plot checkbox.
5
In the table, enter the following settings:
Plot group
Plot window
Electron Density (plas) 1
Graphics
6
In the Study toolbar, click  Compute.
Results
Electron Density (plas) 1
Add a study to do an O2 mole fraction sweep. This study will use a solution of the previous study as initial conditions to save computation time.
The Initial damping factor is again set to 1 because a solution is used as initial condition and because the next parameter will use the previous solution.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select Preset Studies for Selected Multiphysics > Frequency–Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 3
Step 1: Frequency–Stationary
1
In the Settings window for Frequency–Stationary, locate the Study Settings section.
2
In the Frequency text field, type 13.56[MHz].
3
Locate the Values of Dependent Variables section. Find the Initial values of variables solved for subsection. From the Settings list, choose User controlled.
4
From the Method list, choose Solution.
5
From the Study list, choose Power Sweep, Frequency–Stationary.
6
From the Parameter value (Pw (W)) list, choose 500 W.
7
Locate the Study Extensions section. Select the Auxiliary sweep checkbox.
8
From the Sweep type list, choose All combinations.
9
Click  Add.
10
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
xO2
 
 
11
In the table, click to select the cell at row number 1 and column number 2.
12
Click  Range.
13
In the Range dialog, type 0.9 in the Start text field.
14
In the Step text field, type -0.1.
15
In the Stop text field, type 0.1.
16
Click Replace.
17
In the Settings window for Frequency–Stationary, locate the Study Extensions section.
18
Click  Add.
19
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
Pw
500
W
20
In the Model Builder window, click Study 3.
21
In the Settings window for Study, type xO2 Sweep in the Label text field.
22
In the Study toolbar, click  Get Initial Value.
Results
Electric Potential (plas) 2, Electron Density (plas) 2, Electron Temperature (plas) 2, Magnetic Flux Density (mf) 2, Magnetic Flux Density, Revolved Geometry (mf) 2, Pressure (spf) 2, Temperature (ht) 2, Velocity (spf) 2, Velocity, 3D (spf) 2
Right-click and choose Group.
Results
xO2 Sweep
1
In the Model Builder window, expand the xO2 Sweep > Solver Configurations node, then click Results > Group 3.
2
In the Settings window for Group, type xO2 Sweep in the Label text field.
xO2 Sweep
Solution 3 (sol3)
1
In the Model Builder window, expand the xO2 Sweep > Solver Configurations > Solution 3 (sol3) > Stationary Solver 1 node, then click Fully Coupled 1.
2
In the Settings window for Fully Coupled, locate the Method and Termination section.
3
In the Initial damping factor text field, type 1.
4
Locate the Results While Solving section. Select the Plot checkbox.
5
In the table, enter the following settings:
Plot group
Plot window
Electron Density (plas) 2
Graphics
6
In the Study toolbar, click  Compute.
Results
Negative Ion Density
1
In the Results toolbar, click  2D Plot Group.
Create plots for the negative ion density, the power absorbed by the electrons, the space distribution of the charged species along the symmetry axis, and for averaged densities of the charged species as functions of power and oxygen mole fraction.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Dataset list, choose Power Sweep/Solution 2 (sol2).
4
In the Label text field, type Negative Ion Density.
Surface 1
1
Right-click Negative Ion Density and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type plas.n_wO_1m.
4
In the Negative Ion Density toolbar, click  Plot.
Absorbed Power
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Dataset list, choose Power Sweep/Solution 2 (sol2).
4
In the Label text field, type Absorbed Power.
Surface 1
1
Right-click Absorbed Power and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type mf.Qrh.
4
Locate the Coloring and Style section. From the Color table list, choose ThermalWave.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
Select Domain 2 only.
3
In the Absorbed Power toolbar, click  Plot.
Charged Species Along Axis-of-Symmetry
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Charged Species Along Axis-of-Symmetry in the Label text field.
3
Locate the Data section. From the Dataset list, choose Power Sweep/Solution 2 (sol2).
4
From the Parameter selection (Pw) list, choose Last.
5
Click to expand the Title section. From the Title type list, choose None.
6
Locate the Plot Settings section.
7
Select the y-axis label checkbox. In the associated text field, type Number density (1/m<sup>3</sup>).
8
Locate the Legend section. From the Position list, choose Lower middle.
Line Graph 1
1
Right-click Charged Species Along Axis-of-Symmetry and choose Line Graph.
2
Select Boundary 3 only.
3
In the Settings window for Line Graph, locate the x-Axis Data section.
4
From the Parameter list, choose Expression.
5
In the Expression text field, type z.
6
Click to expand the Legends section. Select the Show legends checkbox.
7
Find the Include subsection. Clear the Solution checkbox.
8
Select the Expression checkbox.
Line Graph 2
1
Right-click Line Graph 1 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type plas.n_wO2_1p.
Line Graph 3
1
Right-click Line Graph 2 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type plas.n_wO_1p.
Line Graph 4
1
Right-click Line Graph 3 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type plas.n_wAr_1p.
Line Graph 5
1
Right-click Line Graph 4 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type plas.n_wO_1m.
4
Click the  y-Axis Log Scale button in the Graphics toolbar.
Space Averaged Charged Species vs. Power
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Space Averaged Charged Species vs. Power in the Label text field.
3
Locate the Data section. From the Dataset list, choose Power Sweep/Solution 2 (sol2).
4
Locate the Title section. From the Title type list, choose None.
5
Locate the Plot Settings section.
6
Select the x-axis label checkbox. In the associated text field, type Power (W).
7
Select the y-axis label checkbox. In the associated text field, type Number density (1/m<sup>3</sup>).
8
Locate the Legend section. From the Position list, choose Lower right.
Global 1
1
Right-click Space Averaged Charged Species vs. Power and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
aveop1(plas.ne)
1/m^3
ne
aveop1(plas.n_wAr_1p)
1/m^3
Ar+
aveop1(plas.n_wO2_1p)
1/m^3
O2+
aveop1(plas.n_wO_1p)
1/m^3
O+
aveop1(plas.n_wO_1m)
1/m^3
O-
4
Click to expand the Legends section. Find the Include subsection. Clear the Solution checkbox.
5
In the Space Averaged Charged Species vs. Power toolbar, click  Plot.
6
Click the  y-Axis Log Scale button in the Graphics toolbar.
Space Averaged Charged Species vs. xO2
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Space Averaged Charged Species vs. xO2 in the Label text field.
3
Locate the Data section. From the Dataset list, choose xO2 Sweep/Solution 3 (sol3).
4
From the Parameter selection (Pw) list, choose First.
5
Locate the Title section. From the Title type list, choose None.
6
Locate the Plot Settings section.
7
Select the x-axis label checkbox. In the associated text field, type xO2.
8
Select the y-axis label checkbox. In the associated text field, type Number density (1/m<sup>3</sup>).
9
Locate the Legend section. From the Position list, choose Lower middle.
Global 1
1
Right-click Space Averaged Charged Species vs. xO2 and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
aveop1(plas.ne)
1/m^3
ne
aveop1(plas.n_wAr_1p)
1/m^3
Ar+
aveop1(plas.n_wO2_1p)
1/m^3
O2+
aveop1(plas.n_wAr_1p)
1/m^3
O+
aveop1(plas.n_wO_1m)
1/m^3
O-
4
Locate the x-Axis Data section. From the Axis source data list, choose xO2.
5
Locate the Legends section. Find the Include subsection. Clear the Solution checkbox.
6
In the Space Averaged Charged Species vs. xO2 toolbar, click  Plot.
7
Click the  y-Axis Log Scale button in the Graphics toolbar.