Physics based lumped element circuit model for nanosecond pulsed dielectric barrier

discharges
Thomas Underwood, Subrata Roy, and Bryan Glaz

Citation: Journal of Applied Physics 113, 083301 (2013); doi: 10.1063/1.4792665

View online: http://dx.doi.org/10.1063/1.4792665
View Table of Contents: http://scitation.aip.org/content/aip/journal/jap/113/8?ver=pdfcov
Published by the AIP Publishing

Articles you may be interested in

Simulations of nanosecond-pulsed dielectric barrier discharges in atmospheric pressure air
J. Appl. Phys. 113, 113301 (2013); 10.1063/1.4795269

Characterization of nanosecond pulse driven dielectric barrier discharge plasma actuators for aerodynamic flow
control
J. Appl. Phys. 113, 103302 (2013); 10.1063/1.4794507

Role of the electric waveform supplying a dielectric barrier discharge plasma actuator
Appl. Phys. Lett. 100, 193503 (2012); 10.1063/1.4712125

Modeling of dielectric barrier discharge plasma actuators driven by repetitive nanosecond pulses
Phys. Plasmas 14, 073501 (2007); 10.1063/1.2744227

Electrohydrodynamic force and aerodynamic flow acceleration in surface dielectric barrier discharge
J. Appl. Phys. 97, 103307 (2005); 10.1063/1.1901841

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 JOURNAL OF APPLIED PHYSICS 113, 083301 (2013)

Physics based lumped element circuit model for nanosecond pulsed

dielectric barrier discharges
Thomas Underwood,1,a) Subrata Roy,1,b) and Bryan Glaz2
1
Applied Physics Research Group, Department of Mechanical and Aerospace Engineering,
University of Florida, Gainesville, Florida 32611, USA
2
Vehicle Technology Directorate, U.S. Army Research Laboratory, Aberdeen Proving Ground, Maryland
21005, USA

(Received 25 December 2012; accepted 5 February 2013; published online 22 February 2013)
This work presents a physics based circuit model for calculating the total energy dissipated into
neutral species for nanosecond pulsed direct current (DC) dielectric barrier discharge (DBD) plasmas.
Based on experimental observations, it is assumed that the nanosecond pulsed DBD’s which have
been proposed for aerodynamic flow control can be approximated by two independent regions of
homogeneous electric field. An equivalent circuit model is developed for both homogeneous regions
based on a combination of a resistor, capacitors, and a zener diode. Instead of fitting the resistance to
an experimental data set, a formula is established for approximating the resistance by modeling
plasmas as a conductor with DC voltage applied to it. Various assumptions are then applied to the
governing Boltzmann equation to approximate electrical conductivity values for weakly ionized
plasmas. The developed model is then validated with experimental data of the total power dissipated
by plasmas. VC 2013 American Institute of Physics. [http://dx.doi.org/10.1063/1.4792665]

I. INTRODUCTION with relevant physics that could be implemented as an approx-

imation for the energy dissipated within a plasma for any
The need for improved control over aerodynamic flow
pulsed direct current (DC) dielectric barrier discharge (DBD)
separation has increased interest in the potential use of
configuration. Among the other goals in this paper is to probe
plasma actuators. The inherent advantages of plasma actuator
into the background processes that occur within plasmas and
flow control devices include: fast response time, surface com-
incorporate that knowledge into the model.
pliance, lack of moving parts, and inexpensiveness. However,
it has been established that the actuators which affect the
flow via directed momentum transfer are not effective at II. LUMPED ELEMENT CIRCUIT MODEL
Mach numbers associated with most subsonic aircraft appli-
cations. Recently, Roupassov et al.1 demonstrated that pulsed One of the primary assumptions in creating this model is
plasma actuators, in which energy imparted to the flow that nanosecond pulsed DBD’s can be approximated by two
appears to effectively control flow separation, seem to be independent regions of homogeneous electric field. One such
suitable at Mach numbers (M  0.3) beyond the capabilities region, dubbed the “hot spot” is the region adjacent to the
of the current plasma induced momentum based approaches. powered electrode. This region makes up a small portion of
Given the fundamental differences between the novel the total discharge area but was observed to be an important
pulsed discharge approach and the more conventional mo- component of the plasma discharge and necessary to obtain
mentum based approaches, there is a need to develop an effec- agreement with experimentally measured shock wave dynam-
tive and efficient model for the energy delivered to the flow ics by Roupassov et al.1 The other region, dubbed the “tail,”
by the plasma. Once calculated, that value can be input to a encompasses the rest of the plasma discharge and extends to
computational fluid dynamics solver as an energy source term the edge of the dielectric. As both regions are independent,
resulting in a coupled fluid/plasma dynamics model. the model presented in this paper consists of a single network
Multiphysics models of this type are required in order to study for each region containing a resistor, capacitors, and a diode.
detailed flow characteristics. However, detailed numerical As shown in Fig. 1, circuit elements that were used to
simulations involving plasma kinetics are computationally model the plasma include: an air capacitor Ca , a dielectric
prohibitive for a variety of coupled fluid/plasma design prob- capacitor Cd , a resistor Rf , and a zener diode Df . The air ca-
lems. To address this issue, efficient circuit element models pacitor represents the capacitance between the dielectric sur-
have been introduced to approximate the complex processes face and the exposed electrode. The dielectric capacitor
within plasmas. Circuit models such as those by Orlov2 rely represents the capacitance between the dielectric surface and
on empirical constants which may not be applicable to nano- insulated electrode and is proportional to the thickness of the
second pulsed discharges. To date, an approximate model of dielectric layer. Thus, the dielectric layer in the form of both
nanosecond pulsed plasma actuators has not been developed. its thickness and the value of its dielectric constant plays an
This paper deals primarily with establishing a flexible model important role in determining the effectiveness of the plasma
actuator. Finally, the zener diode, introduced by Orlov2 is
a)
Undergraduate Student. utilized in the model to enforce an energy threshold value
b)
Electronic mail: roy@ufl.edu. URL: http://aprg.mae.ufl.edu/roy/. below which plasma will not form.

0021-8979/2013/113(8)/083301/11/$30.00 113, 083301-1 C 2013 American Institute of Physics

V

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-2 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

FIG. 1. Electric circuit model of a dielectric aerodynamic plasma actuator.

Since a uniform charge distribution along the top of the

dielectric is assumed, the typical asymmetric 2D plasma actu-
ator geometry featured in Fig. 1 can be simplified to a series
of homogeneous symmetric regions. As shown in Fig. 2, these
regions include: an anode sheath, “hot spot,” and “tail.” This
assumption results in a series of coupled 1D models that
account for the chordwise variation along the actuator. Fig. 3
shows the simplified circuit model within each homogeneous
region.

FIG. 3. Region of homogeneous potential, i.e., “hot spot.”

A. Circuit
As displayed in Fig. 1, the lumped element circuit is a 0 d Ad
function of the two capacitance values, Ca and Cd . In this Cd ¼ ; (2)
hd
model, the air is treated as both a conductor to generate a
physical relationship for the resistance Rf and a parallel plate where Ad is the cross-sectional area of the dielectric capaci-
capacitor to generate Ca . An advantage of modeling the tive element and hd is the height of the dielectric barrier
plasma as a conductor in addition to a parallel plate capacitor layer. As displayed in Fig. 5, Ad is the product of the span-
is that it generates a physical relationship for the resistance, wise length of the actuator za and dd , the width of the dielec-
Rf , a value that is traditionally empirically determined. The tric region.
air gap capacitor can be modeled as3 Treating the plasma as a conductor, the resistance for DC
voltage is proportional to rp , Aa , and ha : It can be given as3
0 a Aa
Ca ¼ ; (1) ha
ha Rf ¼ : (3)
rp Aa
where Aa is the cross-sectional area of the air and ha is the
approximate height of the plasma region of interest. The Starting from Kirchoff’s circuit laws,3 the governing dif-
height of the plasma has been shown by Roupassov et al.1 to ferential equation for the voltage drop experienced by the air
be approximately independent of applied voltage for nano- gap, DV, is given by
second pulsed DBD actuators. As displayed in Fig. 4, Aa is
 
the product of the spanwise length of the actuator za and la , dDVðtÞ dVapp Ca DVðtÞ
the chordwise distance from the exposed electrode to the end ¼ 1 j ;
dt dt Ca þ Cd Rf ðtÞðCa þ Cd Þ
of the dielectric region.
(4)
The capacitive element corresponding to the dielectric

can be modeled as3 1 if jEj > Ecrit
j¼ (5)
0 if jEj  Ecrit ;

where Vapp is the applied voltage and j is the contribution

from the zener diode. If the electric field magnitude, given as

FIG. 2. Plasma discharge regions analyzed in this study. FIG. 4. Sketch of the capacitive air element.

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-3 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

the physical evolution in the number of particles, f ðv; J; r; tÞ,

defined such that f ðv; J; r; tÞ dv is the number of particles in
a unit volume located at point r, time t, internal quantum
number J, and differential velocity range v þ dv. Using this
distribution function, the number of particles at point r and
FIG. 5. Sketch of the dielectric capacitive element. time t can be defined as

jEj ¼
jDVj
; (6) Xð
ha Nðr; tÞ ¼ f ðv; J; r; tÞdv: (9)
J
2
is greater than the breakdown electric field, Ecrit , required to
ionize air, then j takes on a value of one, otherwise it is zero This distribution function allows a mathematical
to signify that plasma has not formed. description to be developed for the temporal evolution in the
number of particles resulting from particle collisions within
a control volume. The time rate of change in the number of
B. Conductivity
particles due to externally applied fields can be described as5
To effectively calculate the resistance governed by Eq.
(3), an expression must first be developed for the electrical Df f ðv þ dv; J; r þ dr; t þ dtÞ  f ðv; J; r; tÞ
¼ : (10)
conductivity, rp , of the plasma. This value is one that tradi- Dt dt
tionally requires a numerical approach. To simplify the prob-
lem to a point where an analytic formulation can be used, A partial differential equation can be developed to
numerous simplifying assumptions were used and are describe Eq. (10) using formulas for the time rate of change of
described in the following paragraphs. v and r established using the equation of motion in the form of
For any plasma, the resulting electric current is com-
posed of two primary terms: the current from electrons and dv F
¼ ; (11)
that from ions. As the drift velocity of electrons, we , in a dt m
non-equilibrium plasma is significantly higher than ions, the
current density can be approximated as the portion from dr
¼ v; (12)
electrons if the number densities, Ne and Ni , are approxi- dt
mately the same. Using a form of the generalized Ohm’s
and the chain rule of calculus to obtain the Boltzmann ki-
Law, the current density vector, J and rp , respectively, can
netic energy equation given as
be written as

J  eNe we ¼ rp E; (7) Df @f @f F @f
¼ þv þ  : (13)
Dt @t @r m @v
rp ¼ eðNe le þ Ni li Þ; (8)
In order to approximate the total derivative, a relaxation
time can be introduced defined as the time taken for the sys-
where li and le represent the ion and electron mobilities,
tem to be reduced to an equilibrium distribution function. The
respectively. Much like Eq. (7), the electrical conductivity
tau approximation can be given as s  ðNa rea jvjÞ1 , where
relation can be simplified using the concept of quasineutral-
rea is the collision cross section between electrons and atoms,
ity which is defined as having approximate equal number
v is the average collision velocity, and Na is the number den-
densities for charged particles of opposite polarity. Thus, as
sity of atoms.5 If the number of particles within a control vol-
le is typically three orders of magnitude larger than li , it is a
ume is defined as a equilibrium distribution function, f0 , when
good assumption to approximate the electrical conductivity
each particle is at the same energy level as its neighbor, then
as only coming from electrons as long as Ne is at least of the
the particle evolution over time due to pairwise collisions af-
same order of magnitude as Ni .4 Quasineutrality is a typical
ter an external force has been applied can be given as5
assumption that is valid as long as the plasma being model-
ing is far enough away from the powered electrode to avoid Df f  f0
the boundary layer in plasma physics called the sheath. ¼ : (14)
Dt s
Since a pulsed DC voltage is assumed, the activation of
the external electric field will follow the voltage waveform As s only accounts for collisions between electrons and
as a step function. Thus, two expressions will be required for neutral atoms, it is only accurate in the event of weakly ion-
rp , where the first is valid for the period when an external ized plasma. Plasma actuators considered in this paper tradi-
electric field is applied, as shown in Fig. 9 from 0 to 60 ns, tionally feature a low degree of ionization, or simply the
and the second when the voltage drop over the air gap is amount of air that is ionized, and thus can be treated in a
zero. For the portions of the voltage waveform that DV is weakly ionized limit. In terms of the momentum-transfer
zero, the power is also zero according to Ohm’s Law and collision frequency which can be defined as the mass cor-
thus the conductivity during this time is of no importance. rected rate at which a particle of a specific species collides
To generate a analytic formula for the electrical conduc- with another, the criteria for a weakly ionized plasma can be
tivity, a distribution function must be introduced to describe given as4

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-4 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

 ei   en : (15) neglected. The electrical conductivity and electron mobility,

respectively, can be given as
This equation requires that the collision frequency ð t  
between electrons and ions be much less than those between  Ne ðt Þe2 t  t
rp ðt Þ ¼ EðtÞexp    dt; (19)
electrons and neutrals. Thus, if this requirement is met, the me Eðt Þ 0 Na rea 
collisional occurrences between electrons and charged par-
ticles can be effectively ignored and the relaxation time ð t  
e t  t
established is a good approximation of the total time rate of le ðt Þ  EðtÞexp    dt; (20)
me Eðt Þ 0 Na rea 
change in the number of particles within a control volume.
After introducing a relaxation time to approximate Eq. where rea is a function of electron energy and can be
(13), an equation of motion describing the average velocity obtained for various molecules found in air from Phelps.6
of electrons can be established by multiplying by me v and Numerical values used in the model are included in the
integrating over the electron velocity. An analytic formula- Appendix. The electron velocity can be obtained by assuming
tion for the average velocity an electron experiences due to a Maxwellian velocity profile. As a collection of electrons
an externally applied field, we , can be obtained by assuming within plasma have a range of velocities, the Maxwellian ve-
that the force term can be approximated as the Lorentz force locity profile represents the most probable distribution of
e these velocities. Thus, the distribution of velocities,
F ¼ eE  ðv  BÞ; (16) f ðu; v; wÞ, can be given by4
c
2 3
1
where B represents the magnetic field and E represents the mðu2 þ v2 þ w2 Þ
6 7
electric field. As plasma actuators have no applied magnetic f ðu; v; wÞ ¼ A3 exp4 2 5; (21)
kb T e
field, B can be set equal to zero. Therefore, the equation of
motion for an electron describing we can be given as
 12
m
A3 ¼ Ne : (22)
me dwe ðtÞ me we ðtÞ 2pkb Te
þ ¼ eEðtÞ: (17)
dt s
Using the non-relativistic definition of kinetic energy, a
Solving Eq. (17), a linear first-order differential equation, an relationship can be established for the kinetic energy of an
integral equation is obtained electron that is valid in the limit jvj  c, where c is the speed
of light. Using this approximation, the relativistic formula-
  ð t  
 e t t tion of the kinetic energy can be approximated as
we ðt Þ ¼  exp   EðtÞexp   dt: (18)
me s s
0 mc2 1
Ek ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  mc2  mv2 : (23)
2
Equation (18) can be solved in conjunction with Eq. (7) 1  jvj2 =c2
to obtain an expression for the time varying conductivity and
in conjunction with Eq. (8) to obtain an approximation for Averaging Eq. (21) and using the non-relativistic definition of
the time varying electron mobility if the ion mobility is kinetic energy, the mean kinetic energy of electrons becomes4

2 3
ð ð ð1 1 2 2 2
1 me ðu þ v þ w Þ
6 7
A3 me ðu2 þ v2 þ w2 Þexp4 2 5 du dv dw
1 2 k T
b e
Eav ¼ 2 3 ; (24)
ð ð ð1 1
me ðu2 þ v2 þ w2 Þ
6 7
A3 exp4 2 5 du dv dw
1 kb Te

where kb is Boltzmann’s constant. From the definition of ki- The required inputs for Eqs. (19) and (20) include: Ne , the
netic energy, the relationship between Eav and jvj with vector number density of electrons, Na , the number density of atoms,
components ðu; v; wÞ can be established and the average ther- and Te , the temperature of the electrons. Among these values,
mal velocity of electrons becomes Na can be assumed to be constant in time as the number den-
sity of atoms is significantly higher than that of free electrons.
sffiffiffiffiffiffiffiffiffiffiffi Many other models incorporate a constant electron
3kb Te temperature into their model.2,7 Using experimentally
jvj ¼ : (25)
me measured values of reduced electric field,8 E=Na vs.

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-5 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

electron temperature as detailed in Fig. 6, this model calcu- where A and B are empirical constants that have tabulated
lates a new electron temperature at each time step by com- values of 15 and 365, respectively, for air at atmospheric
paring the E=N value experienced by the plasma, produced pressure.9
from using Eqs. (4)–(6), with Fig. 6. It is assumed in this
paper that the effects of an applied electric field have an in-
C. Discharge development
stantaneous, on a time scale much faster than 109 s, effect
on electrons. The final remaining unknown required to close the equa-
A time-varying differential equation that governs Ne can tion system is to determine how the electric potential
be obtained from the drift-diffusion equations. The electron changes over the horizontal length of the actuator. To accu-
continuity equation that governs Ne can be given by rately model this, it is important to incorporate the wall
effects of the plasma actuator. In terms of potential variation,
@Ne these wall effects attract charged particles of opposite polar-
þ r  Ce ¼ ajCe j  bni ne ; (26)
@t ity and shield charged particles of the same polarity.
Therefore, for regions beside the anode (powered electrode
where Ce is called the charged species flux, a is the for positive pulses), an anode sheath is developed where an
Townsend coefficient of ionization, and b is the recombina- attraction of electrons occurs and a repulsion of positive ions
tion coefficient between electron and neutral atoms. By sim- occur. For regions above the dielectric region on top of the
plifying the plasma discharge into a combination of two grounded electrode, a cathode sheath is developed where
homogeneous regions plus an anode sheath and by invoking positive ions are collected and electrons are repelled. Fig. 7
the irrotational property of electric fields, Eq. (26) can be illustrates the two predominate sheaths that are developed in
simplified by ignoring any spatial variation in the number asymmetric plasma actuators.
density of electrons, i.e., Ce becomes It is important to consider the effects of charge collection
and repulsion in these regions as such phenomenon can have
Ce  Ne le jEj: (27) large effects on the variation of the electric potential over the
chordwise length and height of the plasma discharge.
Assuming the number densities of electrons and ions are
equal, i.e., in the quasineutral region, the electron continuity
equation for air, with a composition of 80% N2 and 20% O2, 1. Wall effects
can be written as As the model introduced in this paper is interested in
dNe solving for the amount of energy the plasma transfers to neu-
 aair jNe le Ej  0:80bN2 Ne2  0:2b02 Ne2 ; (28) tral species, the relative importance of both the anode and
dt
cathode sheath regions needs to be established. Energy, Q, is
where aair and b, respectively, can be approximated as hav- a function of electrical conductivity and electric field
ing the form9 strength and can be given as
 
aair ¼ Ap exp 
Bp
; (29) Q ¼ rp jE2 jVvol : (32)
E
  Equation (32) shows that the energy transfer from both
7 300 sheath regions is proportional to the conductivity of the
bN2 ¼ 2:8  10 ; (30)
Te plasma within each region. As the cathode sheath region typ-
  ically has very few electrons due to repulsion effects, the
7 300 current within this region will be carried by positively
b02 ¼ 2  10 ; (31)
Te charged ions. The mobility of such ions are significantly less
than that of an electron, so as a first order approximation if
the local charged species number density and electric field
strength are assumed equal, the conductivity in the cathode
sheath will be significantly less than the conductivity in a
bulk plasma and the anode sheath, i.e.,

FIG. 7. Collective wall effects of the exposed powered electrode and virtual
FIG. 6. Plot of electron temperature vs. reduced electric field. electrode (dielectric region). (1) Anode sheath and (2) cathode sheath.12

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-6 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

rsc  rp  rsa : (33) equal number densities of both electrons and ions (ne ¼ ni ),
while the cathode and anode sheaths only have electron and
When combined with the fact that the volume of both ion densities, respectively. If this is assumed then the anode
sheath regions are orders of magnitude less than the bulk sheath can be approximated as having a net charge density
plasma, an order of magnitude approximation to the plasma equal to the quasineutral region’s calculated electron number
discharge can be obtained by ignoring the cathode sheath’s density. By assuming the anode sheath is devoid of any posi-
effects for nanosecond pulsed plasma discharges. Although tive ions and has a time-dependent number density of elec-
the cathode sheath is approximately negligible in terms of trons, Poisson’s equation for the sheath can be given by
energy transfer, the higher electric field strength and higher
conductivity present in the anode sheath contain important d2 / eNe ðtÞ
physics necessary for capturing the energy transferred by a ¼ : (41)
dx2 0
plasma discharge. Fig. 2 shows the updated “hot spot” and
“tail” regions with the anode sheath included in the hot spot If this form is assumed for a single Debye length (kD ), then
region adjacent to the powered electrode. the remaining discharge (rest of “hot spot” and “tail”) region
is part of the quasineutral bulk plasma and can be given by
2. Electric potential variation Laplace’s equation
To establish a variation in the electric potential, the gov-
erning Maxwell equations can be used which are written as3 d2 /
¼ 0: (42)
dx2
qf
rE¼ ; (34)
0 The ordinary differential equations for electric potential vari-
ation in the hot spot and tail regions, respectively, can be
r  B ¼ 0; (35) solved to provide an approximate 1D spatial variation
@B 8
rE¼ ; (36) < eNe x2 þ C x þ C if x  k
@t 1 2 D
/ðxÞ ¼ 20 (43)
@E :
r  B ¼ l0 J þ l0 0 ; (37) C3 x þ C4 if x > kD :
@t
Equation (43) requires a total of 4 boundary conditions.
in differential form where qf ¼ eðni  ne Þ is the net charge Those can be summarized as
density and l0 is the permeability of free space. In the ab-
sence of a time-varying magnetic field, Eq. (36), Faraday’s /1 ðx ¼ 0Þ ¼ Vapp ; (44)
Law of Induction, simplifies to
/2 ðx ¼ LÞ ¼ Vbreak ; (45)
r  E ¼ 0; (38)
/1 ðx ¼ kD Þ ¼ /2 ðx ¼ kD Þ; (46)
and since the curl of E is zero, the electric field can be solved
for as a potential function / and substituted into Gauss’ d/1 d/
Law, Eq. (34). The resulting equations, respectively, can be ðx ¼ kD Þ ¼ 2 ðx ¼ kD Þ; (47)
dx dx
given by
where /1 is the potential in the anode sheath, /2 is the poten-
E ¼ r/; (39) tial in the quasineutral region, L is length of the actuator, and
qf kD is a Debye length. The first boundary condition is Vapp at
r2 / ¼  : (40) the cathode and the second is based on the experimental obser-

vations by Roupassov et al.1 It was observed that a plasma dis-
The net charge density within the quasineutral region of a charge could be approximated as stopping at the edge of the
plasma is equal to zero as ne ¼ ni . For the anode sheath, the grounded electrode for asymmetric actuators independent of
net charge density can be approximated if the plasma is the applied voltage; so, the edge represents the absolute limit
assumed to uniformly distribute its charge density. Therefore, of ionization or the breakdown voltage of air. Using Eqs.
the quasineutral region of the “hot spot” and “tail” features (44)–(47) as boundary conditions, the potential becomes

8  
>
> eNe 2 eNe 2 eNe Vbreak  Vapp
< 20 x þ 20 L kD  0 kD þ
>
L
x þ Vapp if x  kD
/ðxÞ ¼   (48)
>
> eNe 2 Vbreak  Vapp
>
: kD þ ðx  LÞ þ Vbreak if x > kD :
20 L L

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-7 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

The Debye length provides an order of magnitude TABLE I. Experimental parameters.1

approximation for the extent of a plasma sheath by assuming
an exponential Boltzmann distribution in the charge density ha 0.4 mm
hd 0.3 mm
within the plasma discharge. Substituting this into Poisson’s
d 2.7
equation,
a 1
    An 30 mm2
d2 / e/
0 2 ¼ eN1 exp 1 ; (49) Ad 30 mm2
dx kb T e V 50 kV
DT 60 ns
where N1 is the charged particle density far away from the
electrode. Taking a first-order Taylor expansion, the Debye
length can be given as4
 12
 0 kb T e which allows a numerical ODE solver, such as those employ-
kD : (50) ing the Dormand-Prince method to minimize functional error
N 1 e2
by adjusting the step size after each time step.10 Numerical
Although at high voltages, a first-order approximation fails, integration required for Eqs. (19) and (20) and energy
Eq. (50) still provides an order of magnitude approximation derived from the circuit model, given as
of the extent of the anode sheath. ðt
DV 2
Q¼ dt; (60)
D. Numerical procedure 0 Rf ðtÞ

When solving for the energy imparted to neutral species, were performed via the Gauss-Kronrod quadrature method10
a coupled equation system results from Eqs. (4) and (28). at each time step.
The coupled terms include
III. RESULTS
Rf / Ne ; (51)
A. Validation against experiment
Ce / DV; (52)
To validate the accuracy of the model described in this
b / DV; (53) paper, comparisons with data presented in Roupassov
et al.1 are provided. The experimental parameters that were
/1 / Ne : (54) mentioned and used in the circuit model are given in Table
I. Reference 1 uses the electrode configuration detailed in
To solve the resulting equation system, the Dormand-Prince Fig. 8.
Runge-Kutta method was employed. This method provides Fig. 9 illustrates an approximation of the applied voltage
an efficient way to incorporate an adaptive step size that is square wave that was introduced in the experimental work of
important for computational efficiency in a problem that Roupassov et al.1 The slope that is introduced is to simulate
requires small time steps for convergence. The benefit of 0
a function that is differentiable. This is needed for Vapp ðtÞ in
such a procedure can be illustrated through a simplistic the governing differential equation as detailed by Eqs. (4)
example. If the error of each time step is defined as and (5). One could also generate a continuous function using
Fourier decomposition of a traditional square wave; how-
ik ¼ yðtÞ  yik ; (55) ever, this would not account for the minor rise and fall times
found in experimentation.
then if two step sizes are considered, h1 and h2 , the error of
each iteration and their relative error, respectively, can be
given as

yðtÞ  y1n1 ¼ 1n1 ¼ ah1 ; (56)

yðtÞ  y2n2 ¼ 2n2 ¼ ah2 ; (57)

y2n2  y1n1 ¼ aðh1  h2 Þ: (58)

Therefore, for a given error tolerance, , a sequence of step

sizes can be generated

ðhi  hiþ1 Þ FIG. 8. Experimental scheme used by Roupassov et al.1 for the discharge
hiþ2 ¼ q iþ1 ; (59) gap. (1) High-voltage electrode; (2) dielectric layer; (3) low-voltage elec-
jyn iþ1  yn ii j trode, (4) zone of discharge propagation, and (5) insulating plane.

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-8 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

FIG. 9. Plot of the input voltage, Vapp vs. time used in the model.
FIG. 10. Plot of electron number density vs. time for the “hot spot.”

B. “Hot spot” results

finer the relative error tolerances are required to be to guar-
For the small region of 0.4 mm  0.4 mm over all time
antee convergence of a solution.
outside of the anode sheath, the time variation in the number
Using the results displayed in Figs. 10 and 11, the total
density can be established. As an initial condition for Eq.
power dissipated to neutral species as a function of time can
(28), 1015 m3 electrons were assumed based on work by
be established. Using the result of Eq. (4) and the relation-
Ref. 1. As displayed in Fig. 10, there is a large gradient that
ship for the instantaneous power, the time-varying power
occurs during the rise time that peaks around 1:13  1019
imparted to the flow can be given as7
m3. It is also evident in Fig. 10 that the recombination of
electrons is largely negligible on the nanosecond time scale. DV 2 ðtÞ
If a frequency of 1 kHz is used, the recombination of elec- PðtÞ ¼ : (61)
Rf ðtÞ
trons with atoms allows a steady-state electron number den-
sity to be achieved on the nanosecond time scale. The As shown in Fig. 12, the instantaneous power is domi-
recombination of electrons becomes a significant quantity nate during the rise time for the “hot spot” region. Upon inte-
when exploring the dynamics of a plasma discharge on the grating the instantaneous power over time using Eq. (60),
microsecond time scale. this model produces an energy value of 2.1 mJ for this
Using Eq. (4) and Fig. 6, the model was able to produce region. When compared to the experimentally determined
a time-varying electron temperature. As displayed in Fig. 11, value of 4.2 mJ by Roupassov et al.,1 this model produces an
there is an initial spike in the electron temperature to 29 eV order of magnitude estimate for the energy imparted to neu-
(340 000 K) that coincides with the peak in electron num- tral species in this region.
ber density at 11 ns. Fig. 11 also suggests that the assumption
of a constant electron temperature is not an accurate assump-
C. “Hot spot” sheath results
tion for nanosecond pulsed DBD plasmas. The variability in
the electron temperature beyond 40 ns is due to the numerical Using the Debye length approximation provided by Eq.
error tolerances that are selected when solving Eqs. (4), (19), (50) and the solution from the quasineutral “hot spot” region
and (28). Fig. 11 shows that when reducing the relative error ðne 1019 m3, Te 10 eV), a Debye length of 7.43 lm is
in the Dormand-Prince method from 102 to 103 , the strong generated. The electron number density experienced in the
functional variability experiences a significant reduction. anode sheath is assumed to be equal to the values calculated
The higher the gradients are during the rise and fall time, the in the adjacent “hot spot” region. Fig. 13 shows the time

FIG. 11. Plot of electron temperature vs.

time for the “hot spot”: (a) 102 error
and (b) 103 error.

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-9 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

FIG. 12. Plot of power vs. time for the “hot spot.” FIG. 14. Plot of electron temperature vs. time for the “hot spot” sheath.

variation in the electron number density within the sheath

and also shows a peak number density of 1:13  1019 m3 .
Using the revised form of the electric potential within
the anode sheath given by Eq. (48) and using Fig. 6, the elec-
tric temperature variation over the duration of the pulse can
be established. As displayed in Fig. 14, there is an initial
spike in the electron temperature to (348 000 K) that coin-
cides with the peak in electron number density at 11 ns much
like the hot spot region.
As shown in Fig. 15, the instantaneous power is dominant
during the rise time for the anode sheath region. Upon inte-
grating the instantaneous power over time using Eq. (60), this
model produces an energy value of 0.045 mJ for this region.
This number is quite small compared to the 2.1 mJ experi-
FIG. 15. Plot of power vs. time for the “hot spot” sheath.
enced in the “hot spot” region. However, the anode sheath
does have a slightly higher linear energy density than the “hot D. “Tail” results
spot” region (6 J/m vs. 5.25 J/m). The reason that there is not
significant deviation predicted in the anode sheath and quasi- For the region 0.4 mm  4.6 mm, the time variation in
neutral “hot spot” regions is that explicit charge buildup is not the number density can be established. As an initial condi-
accounted for in the circuit model within the sheath region. tion for Eq. (28), 1015 m3 electrons were assumed, the same
As shown by Ref. 11, the electric field in the sheath and adja- number of electrons assumed for the “hot spot” region.
cent quasineutral regions are approximately equal until charge When comparing Figs. 10 and 16, the tail region experiences
buildup is allowed within the anode sheath during the length a lower growth rate in the number of electrons which is due
of the pulse or a series of pulses. to the lower electric field experienced by this region. As
described by Eq. (28), a lower electric field produces a lower

FIG. 13. Plot of electron number density vs. time for the “hot spot” sheath. FIG. 16. Plot of electron number density vs. time for the “tail.”

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-10 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

TABLE II. Comparison between calculated and experimentally measured

energy deposition.

Circuit model (mJ) Experimental1 (mJ) Abs. error (%)

Hot spot’ 2.1 4.2 50

Tail’ 6.6 8 17.5
Total 8.7 12.2 28.7

accuracy for the tail region can likely be attributed to the

larger distance from the cathode. This increase in distance
improves the assumptions of quasineutrality and that the spa-
tial diffusion of charged species is negligible. The region
close to the cathode features complicated ion and electron
buildup and as the relative distance from the cathode
FIG. 17. Plot of electron temperature vs. time for the “tail.” increases, its impact on the problem becomes negligible.
Table II summarizes the results obtained using the cir-
cuit model presented in this paper and associated experimen-
number of ionizations and therefore a more gradual rise and
tal measurements made by Roupassov et al.1 for both the
lower total peak in Ne , approximately 3:4  1017 m3.
“hot spot” and “tail” regions.
Using Eq. (4) and Fig. 6, the model was able to produce
a time-varying electron temperature. As displayed in Fig. 17, IV. CONCLUSION
the traditional assumption of 1 eV (11 600 K) does not
agree well with the results obtained in this model for the tail A new lumped element circuit model was presented that
region on the nanosecond time scale. Instead, Fig. 17 sug- is valid for pulsed DC DBD plasmas. The model approxi-
gests that the peak electron temperature is achieved during mates the total energy dissipated into neutral species using a
the rise time of the pulse, 220 000 K (19 eV) and then lumped element circuit while containing relevant plasma
trends downward during the plateau portion of the voltage physics in the form of a variable electron temperature and
waveform. Much like Fig. 11, the highest electron tempera- number density. An approximate expression was formulated
tures are achieved during the initial high gradient of the using the conductivity of the discharge to calculate the resist-
pulse. Fig. 17, unlike Fig. 11, also shows a rise in electron ance value for the air gap. Asymmetric wall effects were
temperature during the fall time of the pulse as well. This is also approximated in the model by including the effect of the
due to the lower electric potential and negative gradient anode sheath. Results of the model were verified against a
experienced in the tail region during the fall time. pulsed DC experiment conducted by Roupassov et al.1 and
As shown in Fig. 18, the instantaneous power is dominant order of magnitude agreement was obtained for the energy
during the plateau portion of the applied voltage pulse for the imparted into the plasma in both the homogeneous “hot spot”
“tail” region. Upon integrating the instantaneous power over region and “tail” region.
time using Eq. (60), this model produces an energy value of
6.6 mJ. When compared to the experimentally determined ACKNOWLEDGMENTS
value of 8 mJ by Roupassov et al.,1 this model produces an
The first author’s work was partially supported by an
absolute error of 17:5%. The significant improvement in
undergraduate internship at the US Army Research Laboratory.

APPENDIX: CROSS-SECTION DATA

The momentum cross sections used for the numerical
approximations in the model are shown in Table III. These
values were obtained by A. V. Phelps.6

TABLE III. Momentum cross sections.6

Kinetic energy (eV) N2 (1016 cm2) O2 (1016 cm2)

1 10 7.2
2.1 27 6.6
3 21.7 5.7
4 12.6 5.5
5 10.9 5.6
10 10.4 5
15 11 8.8
20 10.2 8.6
FIG. 18. Plot of power vs. time for the “tail.”

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26
 083301-11 Underwood, Roy, and Glaz J. Appl. Phys. 113, 083301 (2013)

1 7
A. Y. Starikovskii, A. A. Nikipelov, M. M. Nudnova, and D. V. D. Orlov, T. Corke, and M. Patel, “Electric circuit model for aerodynamic
Roupassov, “SDBD plasma actuator with nanosecond pulse-periodic dis- plasma actuator,” AIAA, paper no. 2006–1206, 2006.
8
charge,” Plasma Sources Sci. Technol. 18, 034015 (2009). L. C. Pitchford, J. P. Boeuf, and W. L. Morgan, “User-friendly boltzmann
2
D. Orlov, “Modelling and simulation of single dielectric barrier discharge code for electrons in weakly ionized gases,” in IEEE International Conference
plasma actuators,” In: Aerospace and Mechanical Engineering. University on Plasma Science, 1996 IEEE Conference Record-Abstracts (1996), p. 154.
9
of Notre Dame; 2006. Y. P. Raizer, Gas Discharge Physics, (Springer, 1991).
3 10
D. J. Griffiths, Introduction to Electrodynamics, 3rd ed. (Prentice Hall, G. W. Collins, “Fundamental numerical methods and data analysis,” NASA
Upper Saddle River, NJ, 1999). Astrophysics Data System, 1990, http://ads.harvard.edu/books/1990fnmd.book/.
4 11
M. Mitchner, Partially Ionized Gases (John Wiley & Sons, Inc., I. V. Adamovich, M. Nishihara, I. Choi, M. Uddi, and W. R. Lempert,
Indianapolis, IN, 1973). “Energy coupling to the plasma in repetitive nanosecond pulse dis-
5
B. M. Smirnov, Physics of Ionized Gases (Wiley, New York, 2001). charges,” Phys. Plasmas 16, 113505 (2009).
6 12
S. A. Lawton and A. V. Phelps, “Excitation of the b [sup 1] sigma [sup C. L. Enloe, T. E. McLaughlin, R. D. Van Dyken, K. D. Kachner, E. J. Jumper,
þ][sub g] state of o[sub 2] by low energy electrons,” J. Chem. Phys. 69, T. C. Corke, M. Post, and O. Haddad, “Mechanisms and responses of a dielec-
1055–1068 (1978). tric barrier plasma actuator: Geometric effects,” AIAA J. 42, 595–604 (2004).

[This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP:
128.252.67.66 On: Sun, 21 Dec 2014 19:33:26