PDF

Dielectric Barrier Discharge
Introduction
This model simulates electrical breakdown in an atmospheric pressure gas. Because electrical breakdown is a complicated process, a 1D model is considered. To highlight the physics of the breakdown process, this example uses a simple argon chemistry, which keeps the number of species and reactions to a minimum.
Model Definition
The operating principle for a dielectric barrier discharge is as follows: there is a small gap filled with a gas between two dielectric plates. The gap between the two dielectric plates is typically less than one millimeter. On one of the dielectric plates, a sinusoidal voltage is applied. The other plate is electrically grounded. As the voltage applied to the top plate increases, a stronger electric field forms in the gap between the plates. Any free electrons in the gap1 are accelerated, and if the electric field is strong enough they may acquire enough energy to cause ionization. This can lead to a cascade effect where the number of electrons in the gap increases exponentially on a nanosecond time scale. Electrons created via electron impact ionization rush toward one of the dielectric plates, in the opposite direction to the electric field. An equal number of ions are also generated during electron impact ionization (electrons and ions must be created in equal pairs to preserve the overall charge balance). The ions rush toward the opposite dielectric plate in the same direction as the electric field. As a result, surface charge with opposite sign accumulates on both dielectric plates. This causes the electric field to become shielded from the gas filled gap. In fact, the electric field across the gap cannot exceed the breakdown electric field, which is highly dependent on the gas. The breakdown electric field is also a function of the surface properties of the dielectric material. Surface charge accumulation temporarily terminates the discharge until the field reverses direction and the process repeats in the opposing direction.
Modeling dielectric barrier discharges in more than one dimension is, of course possible, but the results can be difficult to interpret due to the amount of competing physics in the problem. In this simple model the problem is reduced to 1D by assuming the dielectric gap is much smaller than the diameter of the plates. It also makes it possible to quickly gain some insight into the characteristics of the discharge without excessive computation time.
The geometry for a typical dielectric barrier discharge is shown in Figure 1. The dielectric plates may be up to 15 cm in diameter; the dielectric and gap thickness are typically less than 1 millimeter.
Figure 1: Graphic illustration of a typical dielectric barrier discharge.
The electron density and mean electron energy are computed by solving a pair of drift-diffusion equations for the electron density and mean electron energy. Convection of electrons due to fluid motion is neglected. For detailed information on electron transport, see Theory for the Drift Diffusion Interface in the Plasma Module User’s Guide.
The electron source Re and the energy loss due to inelastic collisions Rε are defined later. The electron diffusivity, energy mobility, and energy diffusivity are computed from the electron mobility using:
The source coefficients in the above equations are determined by the plasma chemistry using rate coefficients. Suppose that there are M reactions that contribute to the growth or decay of electron density and P inelastic electron-neutral collisions. In general, P >> M. In the case of rate coefficients, the electron source term is given by:
where xj is the mole fraction of the target species for reaction j, kj is the rate coefficient for reaction j (SI unit: m3/s), and Nn is the total neutral number density (SI unit: 1/m3). The electron energy loss is obtained by summing the collisional energy loss over all reactions:
where Δεj is the energy loss from reaction j (SI unit: V). The rate coefficients can be computed from cross section data by the following integral:
where γ = (2q/me)1/2 (SI unit: C1/2/kg1/2), me is the electron mass (SI unit: kg), ε is energy (SI unit: V), σk is the collision cross section (SI unit: m2), and f is the electron energy distribution function. In this case, a Maxwellian EEDF is assumed.
For nonelectron 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 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:
Surface charge accumulation is added to the dielectric surfaces that are adjacent to the gap where the plasma forms by way of the following boundary condition:
where ρs is the surface charge density, which is computed by solving the following distributed ODE on the surfaces:
where n · Ji is the normal component of the total ion current density at the wall, and n · Je is the normal component of the total electron current density at the wall.
The discharge is driven by a sinusoidal electric potential applied to the exterior boundary of one of the dielectric plates:
where the applied peak voltage, V0 is 750 V and the angular frequency, the RF frequency being 50 kHz. The exterior boundary of the other dielectric plate is grounded.
plasma chemistry
Argon is an attractive gas to use in a benchmark problem since only a handful of reactions and a few species need to be considered. The list of chemical reactions considered is as follows (electron impact cross-section are obtained from Ref. 1):
Initially a small number of seed electrons are present. These are necessary in order to initiate the discharge on the first RF cycle. In addition to the volumetric reactions, the following surface reactions are implemented:
When the ions reach the wall they are assumed to change back to neutral argon atoms and donate their charge to the wall.
Results and Discussion
It is more convenient to analyze the results of this one-dimensional problem by extruding the solution into two dimensions. The extra dimension represents time. In COMSOL Multiphysics this is accomplished by adding a Parametric Extrusion 1D dataset. The surface plot is convenient because you can immediately see how the variables of interest evolve over time.
The mass fraction of excited argon atoms is plotted in Figure 2. The excited species have a much longer lifetime in the gap than the electrons or ions. This is because the primary mechanism for destruction of excited argon species is de-excitation upon contact with the wall. The excited argon atoms can only reach the wall via diffusion whereas the electrons and ions reach the wall very rapidly due to migration. It is also apparent from Figure 2 that the discharge reaches a periodic steady state solution after only two RF cycles.
Figure 2: Mass fraction of excited argon.
The electric potential is plotted in Figure 3. The voltage is relatively uniform across the discharge gap. This can be seen more clearly by examining the electric field in Figure 4. There is a much stronger electric field in the dielectric materials than in the gap. This is because the surface charge which accumulates on the dielectric surfaces tends to shield out the electric field.
Figure 3: Electric potential (x-axis) vs. time (y-axis).
Figure 4: Electric field across the gap (x-axis) vs. time (y-axis).
Figure 5: Extruded plot of the electron density.
Figure 6: Extruded plot of the mean electron energy.
Implicit in the equations solved for the number of charged particles and electrostatic potential is that the total electrical current is conserved. Mathematically, this means that:
where J is the total plasma current density (SI unit: A/m2), and ρ is the space charge density (SI unit: C/m3). Since electrons and ions are created in equal pairs the time derivative of the space charge density should be approximately zero. For this 1D model, this means the total current density must be constant across the gap at any snapshot in time. In Figure 7, the current density due to electrons (left) and ions (right) is plotted. The current density is not symmetric because of the different secondary emission coefficients used on the dielectric surfaces.
Figure 7: Plot of the electron current density (left) and the ion current density (right) in the discharge, excluding the first RF cycle.
The total plasma current density is plotted in Figure 8. As expected, the total current density is constant across the gap at any point in time.
Figure 8: Plot of the total plasma current density (sum of the electron and ion current density), excluding the first RF cycle. Conservation of charge requires that the total current density be constant across the gap at any point in time.
The total current at the grounded electrode is plotted in Figure 9. In the absence of the plasma, the current would be a perfect cosine wave. However, the presence of the plasma and flow of charged particles leads to a non-sinusoidal current waveform. The instantaneously absorbed power in the plasma is plotted in Figure 10. Time averaging this over 1 RF cycle yields the power absorbed by the plasma. The power is around 16.7 W on one half cycle and 17.7 W on the other half cycle. The difference is because the secondary emission coefficients are different on the upper and lower plates.
Figure 9: Plot of the total discharge current vs. time.
Figure 10: Plot of power vs. time for the dielectric barrier discharge.
Reference
1. Phelps database, www.lxcat.net, retrieved 2017.
Application Library path: Plasma_Module/Direct_Current_Discharges/argon_dbd_1d
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  1D.
2
In the Select Physics tree, select Plasma>Plasma (plas).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Time Dependent.
6
Geometry 1
You start by defining the geometry for the problem. This model has a simple 1D geometry consisting of 3 domains. Two dielectric domains and a gap where the plasma forms.
Interval 1 (i1)
1
In the Model Builder window, under Component 1 (comp1) right-click Geometry 1 and choose Interval.
2
In the Settings window for Interval, locate the Interval section.
3
From the Specify list, choose Interval lengths.
4
In the Left endpoint text field, type -1e-4.
5
6
Click  Build All Objects.
Add some parameters for the plate dimensions and excitation frequency.
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
Definitions
Variables 1
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Now add the charge conservation to the dielectric materials. This means that only a charge conservation equation is solved in the dielectric material and all the plasma components are solved for in the gap between the dielectrics.
Plasma (plas)
Charge Conservation 1
1
In the Model Builder window, under Component 1 (comp1) right-click Plasma (plas) and choose the domain setting Electrostatics>Charge Conservation.
2
Load in the argon cross sections from file. They form the basis of the plasma chemistry under investigation.
Cross Section Import 1
1
In the Physics toolbar, click  Global and choose Cross Section Import.
2
In the Settings window for Cross Section Import, locate the Cross Section Import section.
3
Click Browse.
4
5
Click Import.
6
In the Model Builder window, click Plasma (plas).
7
In the Settings window for Plasma, locate the Cross-Section Area section.
8
In the A text field, type As.
9
Locate the Plasma Properties section. Select the Use reduced electron transport properties check box.
Reaction 1
1
In the Physics toolbar, click  Domains and choose Reaction.
2
In the Settings window for Reaction, locate the Reaction Formula section.
3
In the Formula text field, type Ars+Ars=>e+Ar+Ar+.
4
Locate the Reaction Parameters section. In the kf text field, type 3.3734e8.
Reaction 2
1
In the Physics toolbar, click  Domains and choose Reaction.
2
In the Settings window for Reaction, locate the Reaction Formula section.
3
In the Formula text field, type Ars+Ar=>Ar+Ar.
4
Locate the Reaction Parameters section. In the kf text field, type 1807.
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 check box.
4
Locate the General Parameters section. From the Preset species data list, choose Ar.
Species: Ars
1
In the Model Builder window, click Species: Ars.
2
In the Settings window for Species, locate the General Parameters section.
3
In the x0 text field, type 1e-11.
4
From the Preset species data list, choose Ar.
Now let the initial number density of Argon ions be the same as the initial number of electrons. This forces the plasma to be initially charge neutral.
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 Initial value from electroneutrality constraint check box.
4
Locate the General Parameters section. From the Preset species data list, choose Ar.
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.
4
Locate the Model Inputs section. In the T text field, type 400[K].
The initial number density of seed electrons is very small, only one million free electrons per cubic meter. This corresponds to a near zero conductivity. So, the gap is truly acting as an insulator initially.
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 1e6.
4
In the ε0 text field, type 5.
Secondary emission of electrons is important when studying discharge curves from DBDs. In this example you add a higher secondary emission coefficient on the left wall.
Surface Reaction 1
1
In the Physics toolbar, click  Boundaries and choose Surface Reaction.
2
3
In the Settings window for Surface Reaction, locate the Reaction Formula section.
4
In the Formula text field, type Ar+=>Ar.
5
Locate the Secondary Emission Parameters section. In the γi text field, type 0.01.
6
In the εi text field, type 2.5.
Surface Reaction 2
1
In the Physics toolbar, click  Boundaries and choose Surface Reaction.
2
3
In the Settings window for Surface Reaction, locate the Reaction Formula section.
4
In the Formula text field, type Ar+=>Ar.
5
Locate the Secondary Emission Parameters section. In the γi text field, type 1E-6.
6
In the εi text field, type 2.5.
Surface Reaction 3
1
In the Physics toolbar, click  Boundaries and choose Surface Reaction.
2
3
In the Settings window for Surface Reaction, locate the Reaction Formula section.
4
In the Formula text field, type Ars=>Ar.
Wall 1
1
In the Physics toolbar, click  Boundaries and choose Wall.
2
Surface charge will begin to accumulate when the gas begins to break down. This will cause the electric field to be shielded in the gap. This is the phenomena responsible for terminating the discharge and also the reason why the breakdown voltage cannot be exceeded across the gap. COMSOL automatically computes the amount of surface charge accumulation when the feature is added to the model. The surface charge accumulation is computed by integration the electron and ion fluxes to the wall.
Surface Charge Accumulation 1
1
In the Physics toolbar, click  Boundaries and choose Surface Charge Accumulation.
2
Ground 1
1
In the Physics toolbar, click  Boundaries and choose Ground.
2
Terminal 1
1
In the Physics toolbar, click  Boundaries and choose Terminal.
2
3
In the Settings window for Terminal, locate the Terminal section.
4
In the Terminal name text field, type electrode.
5
In the V0 text field, type Vrf.
Now assign the relative permittivity to the dielectric material and the air gap where the plasma forms.
Materials
Dielectric 1
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
Right-click Material 1 (mat1) and choose Rename.
3
In the Rename Material dialog box, type Dielectric 1 in the New label text field.
4
5
6
In the Settings window for Material, locate the Material Contents section.
7
Mesh 1
There must be sufficient mesh density to resolve the sharp gradients in the electron and ion density in the gap. Therefore you specify that there are 200 elements across the width of the gap.
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Extremely fine.
Edge 1
In the Mesh toolbar, click  Edge.
Distribution 1
1
Right-click Edge 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
5
Locate the Distribution section. From the Distribution type list, choose Predefined.
6
In the Number of elements text field, type 200.
7
In the Element ratio text field, type 5.
8
From the Growth formula list, choose Geometric sequence.
9
Select the Symmetric distribution check box.
10
Click  Build All.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
Click  Range.
4
In the Range dialog box, choose Number of values from the Entry method list.
5
In the Stop text field, type 1e-4.
6
In the Number of values text field, type 201.
7
Click Replace.
8
In the Home toolbar, click  Compute.
Results
Parametric Extrusion 1D 1
1
In the Results toolbar, click  More Datasets and choose Parametric Extrusion 1D.
2
In the Settings window for Parametric Extrusion 1D, locate the Settings section.
3
Clear the Separate levels check box.
4
In the Level scale factor text field, type 50e3.
Now create a new parametric dataset, which ignores the first startup RF cycle so the current density can be visualized later.
Parametric Extrusion 1D 2
1
Right-click Parametric Extrusion 1D 1 and choose Duplicate.
2
In the Settings window for Parametric Extrusion 1D, locate the Data section.
3
From the Time selection list, choose From list. In the Times (s) list, choose 1E-5 through 1E-4.
Excited Argon Mass Fraction
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Excited Argon Mass Fraction in the Label text field.
Surface 1
1
Right-click Excited Argon Mass Fraction and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Mass fractions>plas.wArs - Mass fraction.
3
In the Excited Argon Mass Fraction toolbar, click  Plot.
4
Click the  Zoom Extents button in the Graphics toolbar.
Electric Potential
1
In the Model Builder window, right-click Excited Argon Mass Fraction and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Electric Potential in the Label text field.
Surface 1
1
In the Model Builder window, expand the Electric Potential node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Electric>V - Electric potential - V.
3
In the Electric Potential toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Electric Field
1
In the Model Builder window, right-click Electric Potential and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Electric Field in the Label text field.
Surface 1
1
In the Model Builder window, expand the Electric Field node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Electric>Electric field - V/m>plas.Ex - Electric field, x component.
3
In the Electric Field toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Electron Density
1
In the Model Builder window, right-click Electric Field and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Electron Density in the Label text field.
Surface 1
1
In the Model Builder window, expand the Electron Density node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Electron density>plas.ne - Electron density - 1/m³.
3
In the Electron Density toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Mean Electron Energy
1
In the Model Builder window, right-click Electron Density and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Mean Electron Energy in the Label text field.
Surface 1
1
In the Model Builder window, expand the Mean Electron Energy node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Electron energy density>plas.ebar - Mean electron energy - V.
3
In the Mean Electron Energy toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Electron Current Density
1
In the Model Builder window, right-click Mean Electron Energy and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Electron Current Density in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parametric Extrusion 1D 2.
Surface 1
1
In the Model Builder window, expand the Electron Current Density node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Current>Electron current density - A/m²>plas.Jelx - Electron current density, x component.
3
In the Electron Current Density toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Argon Ion Current Density
1
In the Model Builder window, right-click Electron Current Density and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Argon Ion Current Density in the Label text field.
Surface 1
1
In the Model Builder window, expand the Argon Ion Current Density node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Plasma>Species>Species wAr_1p>Ion current density - A/m²>plas.Jix_wAr_1p - Ion current density, x component.
3
In the Argon Ion Current Density toolbar, click  Plot.
4
5
Click the  Zoom Extents button in the Graphics toolbar.
Total Conduction Current Density
1
In the Model Builder window, right-click Argon Ion Current Density and choose Duplicate.
2
In the Settings window for 2D Plot Group, type Total Conduction Current Density in the Label text field.
Surface 1
1
In the Model Builder window, expand the Total Conduction Current Density node, then click Surface 1.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type plas.Jix_wAr_1p+plas.Jelx.
4
In the Total Conduction Current Density toolbar, click  Plot.
5
6
Click the  Zoom Extents button in the Graphics toolbar.
Terminal Current
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Terminal Current in the Label text field.
Global 1
1
Right-click Terminal Current and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Plasma>plas.I_electrode - Current, Terminal  electrode - A.
3
In the Terminal Current toolbar, click  Plot.
Total Power Deposition
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Total Power Deposition in the Label text field.
Global 1
1
Right-click Total Power Deposition and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Plasma>Power and collisions>plas.Pcap_tot - Total capacitive power deposition, electrons - W.
3
In the Total Power Deposition toolbar, click  Plot.
Global Evaluation 1
1
In the Results toolbar, click  Global Evaluation.
Use the timeavg operator to compute the time averaged power deposition for cycles 2-10.
2
In the Settings window for Global Evaluation, locate the Data section.
3
From the Time selection list, choose Last.
4
Locate the Expressions section. In the table, enter the following settings:
As you can see, the average power deposited to the plasma remains the same for each cycle after only 3 RF cycles. The power is around 16.7 W on one half cycle and 17.7 W on the other half cycle. The difference is due to the fact that the secondary emission coefficients are different on the upper and lower plates.
5
Click  Evaluate.
 

1
There are typically around 1,000,000 m-3 free electrons in air at sea level.