Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
A theoretical investigation in multi-dimensional and non-thermal plasma effects on combustion processes and pollutant remediation
(USC Thesis Other)
A theoretical investigation in multi-dimensional and non-thermal plasma effects on combustion processes and pollutant remediation
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
A theoretical investigation in multi-dimensional and
non-thermal plasma effects on combustion processes and
pollutant remediation
by
Vyaas Gururajan
A dissertation presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In partial fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
May 2016
1
Acknowledgements
I’d like to thank my advisor, Fokion Egolfopoulos, for training me in the art of scientific inquiry,
for curating an intellectually nurturing environment, and for his exemplary strength of character.
My colleagues possess an infectious grade of curiosity that I am grateful for. My discussions with
Jagannath Jayachandran in particular have been fruitful; I’ve met no one else with such a staunch
discipline in scientific honesty and reason. Professor Ronney has been an inspiration from afar; the
content of his papers and lectures have provided a great deal of perspective, and his balance on
the tightrope that separates Combustion Science and Combustion Technology always serves as a
timely reminder. I am indebted to Professor Wittig’s Socratic lectures; they opened my eyes to the
beautiful mathematical structure of quantum mechanics. Last but not least, USC has served as a
warm and cozy shelter, and I am grateful to everyone in charge of running the place.
2
Dedicated to Libraries.
And Amma & Appa of course!
3
Contents
Contents 4
1 Introduction 5
2 A brief survey of computational methods 8
2.1 Dealing with Numerical Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Dealing with Incompressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Dealing with the Splitting of Operators . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 Direct numerical simulations of probe effects in low-pressure flame sampling 16
3.1 Modeling Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 Coupling Combustion and Non-thermal Plasma Kinetics 27
4.1 Electron Temperature and Relaxation Time . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Discharge Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Plasma Kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.4 The Boltzmann Transport Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5 Chemical Pathways in NOx remediation by non-thermal plasma 38
6 Transient plasma effects on the autoignition of DME/O2/Ar and C3H8/O2/Ar
mixtures 43
6.1 Modeling Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
7 Summary and Looking Ahead 58
References 62
4
1 Introduction
Combustion Science comprises the investigation of physical processes ranging over many orders of
magnitude in time scale. One may start with atoms and work their way up to forest fires as follows:
The electrons of atoms arrange themselves around protons in accordance with the Principles of
Quantum Mechanics
1
. Some arrangements aka molecules possess a low potential energy while
others possess a high potential energy
2
. The release of this energy from the latter type can be
brought about by disturbing the arrangement of resident electrons via a collision with another
particle (photon , electron, atom, or another molecule). The rate at which the totality of electrons
rearranges itself is determined by the potential energy landscape of their interaction which is harder
to imagine than to calculate
3
. Practical considerations require that the rate at which these electrons
rearrange about their parent protons be considered for a large number of molecules per unit volume.
Solvingfortheunitarytimetransformationofsuchsystemsisabandonedinfavorofapplyingresults
fromProbabilityTheory
4
; theratesofthevariouselectronrearrangementsshouldnowbethoughtof
asanaveragerate,weightedbythelikelihoodofthecollidingparticlepossessingtherequiredamount
of energy to trigger said rearrangement
5
. The release of energy from such rearrangements serves
to heat the surrounding molecules, i.e it is transformed into the kinetic energy of the surrounding
molecules. The presence of a large number of molecules gives rise to the phenomenon of diffusion;
the kinetic energy of molecules, their momenta, and their identities are transported to other regions
via a series of low energy collisions
6
. In the gas phase, these time scales are characterized by the
mean-free time, i.e. the time spent between collisions. The pressure of such a fluid is equilibrated
1
We owe much of the formalism to Dirac [1]. Studies in covalent bonding were pioneered by Pauling [2].
2
Contrary to most expositions, it is possible to think of the interaction between atoms in a molecule in terms of
force as well [3].
3
The Born-Oppenheimer approximation[4] (nuclear coordinates remain fixed while electronic coordinates are
solved for) makes these calculations feasible.
4
In the parlance of Quantum Mechanics, this is known as a mixed state [5], wherein the probability amplitudes
of various possible states are real numbers to be thought of as statistical weights. Regardless of the representation
of the states, these numbers codify the degree of likelihood to find the system in any one state. Once a particular
representation is chosen (energy, momentum, etc.), these numbers comprise the so called density matrix of that
representation. In the energy representation, the density matrix remains stationary at thermal equilibrium, a result
of Liouville’s Theorem [6]. A very useful simplification affords itself if we assume local and partial equilibrium,
i.e. the system is in thermal equilibrium but not in chemical equilibrium. This allows us to avoid solving for the
evolutionofthedensitymatrixandinsteadassumethatitisquasi-staticwhilesolvingfortheevolutionofmacroscopic
variables like Temperature, Pressure, Concentration, etc. This is the basis of the theory of linear non-equilibrium
thermodynamics, a regime we will avoid straying too far away from.
5
For the gas phase, expectation values for collision frequencies of such processes proceeds via straightforward
calculations once the assumption of thermal equilibrium is made [7]. These calculations must of course be supplied
with the dynamical information contained in the potential energy surfaces [8].
6
Arguably Boltzmann’s greatest contribution was to articulate the mechanism through which large groups of
atoms exhibited continuum behavior [9]. He was the founder of kinetic theory, a subject with the singular aim of
explaining the macroscopic properties of gases through microscopic mechanics.
5
first and foremost mechanically; an acoustic perturbation informs nearby fluid parcels to adjust
their pressures in anticipation of incoming fluid, a purely isentropic process [10]. The remainder
of pressure equilibration takes place due to the dissipation of the macroscopic velocity of the fluid
via viscous forces. Such an abstraction to the continuum aka gasdynamics, yields the necessary
simplifications in order to examine processes upward and downward of this hierarchy of scales; the
variationofthemacroscopicflow-fieldbringsaboutconditionsformicroscopicprocesseslikereaction
(accompanied by emergent diffusion) and vice-versa. The addition of long range forces like gravity
(and electricity if the fluid is ionized) can significantly alter flow-fields in the presence of density
gradients. The class of flows characterized by velocities lower than the speed of light is susceptible
to radiative heat transfer, wherein energy is transferred away from the fluid by photons belonging
to a spectrum that originates in the physical structure of the radiating atoms and molecules [11].
Larger scale relaxation processes can be added to extend the applicability of gasdynamics to cosmic
scales!
This is currently the theoretical treatment for chemically reacting flow in general, of which
Combustion is a special case; it is characterized by a high exothermicity and activation energies
much higher than the mean internal energies of the flow [12]. These two features alone are sufficient
to account for our everyday experience with combustion; it is restricted to thin layers spatially
in inhomogeneous mixtures (flames) and short durations in homogeneous mixtures (ignition delay
time). The theoretical foundations do not change, but our abstractions do: we are interested
in the speeds of chemically propagating fronts (flame speed), their stability, their topology, their
ignitabilityandsoon. AlthoughthesubjectmatterdiscussedinthisThesisisofimmediaterelevance
to Combustion Technology (the practice of taming the Science to make human related ends meet
7
),
we will invoke these abstractions rarely, since we wish to deal with fundamental additions to the
aforementioned theoretical treatment.
Controlling combustion is the foremost concern of Combustion Technology and is essentially
carried out in one of two ways: Altering flow fields and blending/engineering fuels. There is a third
way: altering the rates of chemical reactions by subjecting a fuel-air mixture to a swarm of hot
electrons, not for the purpose of increasing the temperature, but for creating radicals (atoms or
molecules with unpaired valence electrons; highly reactive) via electronic excitations and ioniza-
tions. This form of catalysis, and ultimately control (since the discharge parameters can be varied
externally), has recently been shown to hold promise [15]. Estimating these radical concentrations
would go a long way in aiding the design of combustion systems that utilize non-thermal plasma.
We wish to develop and discuss a computational framework that makes this possible.
An overarching aim of the current work is to integrate the solution procedure used to solve the
Navier-Stokes Equations
8
with the solution procedure used to solve non-thermal plasma dynamics,
7
Two representative perspectives on the subject are Weinberg [13] and Glassman [14].
8
The Navier-Stokes Equations [16], have sentenced us to Computational Fluid Dynamics (CFD), in large part due
6
which is described by the Boltzmann Transport Equation
9
. Each of these pieces is individually
treated in this work; putting them together is an ongoing effort including collaborations with many
other groups. Nonetheless, the specific reasons that make this problem intimidating are worth
discussing; every time a computational method is used, the domain of its applicability needs to be
stated, often in terms of the Physics it is unable to capture.
The outline of our presentation is as follows: We first briefly review key some of the difficulties
inherent in computationally tackling reacting flows. This will be followed by an application of a
chosen methodology to a particularly important problem in Combustion Science: the Molecular
Beam Mass Spectrometry of Burner Stabilized Flames [19]. After that, we briefly specify certain
results from Statistical Mechanics and Kinetic Theory inthe contextof Gas Discharge Physics. This
permits an introduction to the Boltzmann Transport Equation, the governing equation for electron
energy distribution which upon solving, allows one to calculate the reaction rates for various gas
discharge relevant processes. After summarizing the utilized numerical methods, we will conclude
by discussing a couple of applications of this framework: 1) NOx remediation and 2) Ignition delay
reduction for hydrocarbons (with an emphasis on low temperature chemistry).
to their inherent non-linearity (if they were linear, the principle of superposition would come to our rescue, making
solution procedures far simpler). CFD as a subject, with its own parlance and pedantry, began with Richardson
in 1922 [17], when he boldly attempted to predict the whether using finite differences at “points and instants”. His
groundwork is timeless; we have faster machines and cleverer algorithms, but CFD to a large extent is still confined
to the Standard Operating Procedure that he laid down. His introductory words carry a very familiar sentiment
felt by researchers in most fields today: “The scheme is complicated because the atmosphere is complicated. But it
has been reduced to a set of computing forms (what we call programs today). These are ready to assist anyone who
wishes to make partial experimental forecasts from such incomplete observational data as are now available. In such
a way it is thought that our knowledge of meteorology might be tested and widened, and concurrently the set of forms
might be revised and simplified. Perhaps some day in the dim future it will be possible to advance the computations
faster than weather advances and at a cost less than the saving to mankind due to the information gained. But this
is a dream.”
9
The Boltzmann Transport is an evolution equation for a single particle’s energy distribution function, which is
in fact a probability distribution. It is foundational to Kinetic Theory, and underpins the theoretical foundation of
the multicomponent Navier-Stokes equations for ideal gases. See Chapman & Cowling [18].
7
2 A brief survey of computational methods
The most robust combustion codes are those used to solve flow-field variables (concentration, tem-
perature, velocity, etc.) in canonical configurations. These include zero-dimensional batch reactors,
the freely propagating flame, and the strained flames that appear in stagnation flows (semi-infinite
and finite axial domains) [20]. Each of these configurations demand a numerical algorithm that
computes thermodynamic properties, chemical kinetic reaction rates, and detailed transport, while
addressing the numerical stiffness inherent in such systems.
2.1 Dealing with Numerical Stiffness
Numerical stiffness is a term used to refer to the disparity in time-scales inherent in a physical
phenomenon
10
. Combustion kinetics comprise reactants transforming to products via a number
of intermediate steps, each with its own characteristic timescale. These can vary from each other
by as much as four orders of magnitude. Another example is the phenomenon of diffusion; when
solving the diffusion equation
@T
@t
=r
2
T
for the scalar T, one observes that the initial distribution of T tends to spread outward in space as
time increases. This spreading in time can be written as the superposition of Fourier-components,
thus making it clear that the diffusion problem inherently contains a spectrum of time-scales;
slimmer the initial distribution in position space, wider is the distribution in frequency space. In
contrast, a linear advection equation
@T
@t
=urT
is not stiff; the initial condition simply displaces itself without change in shape in a direction
and speed imposed by the value of u.
Numerical stiffness raises the following question: “Do the smallest time-scales in the problem
need to be resolved for an accurate solution?” The answer is often no: in the combustion problem,
many of the highly reactive chemical species, essentially radicals, are in a quasi-steady state for
most of the duration of transience of the other variables; the concentrations of the radicals can
be estimated from algebraic relations that apply for the entire time-step. Similarly, a diffusion
problem may contain time-scales that contain negligible energy in its power spectrum that they can
be neglected without affecting the accuracy of the solution.
However, when an actual computation needs to be performed, the answer depends on the type
of numerical time integration adopted: the answer is yes for explicit integration and no for implicit
10
First articulated by Hirschfelder and Curtiss [21], it can be demonstrated with the toy equation
dy
dt
=
yG(t)
a(t)
.
In the interval t, the left hand side can become very small compared to the right hand side if the ratio
a
t
<< 1.
Then ytG(t); the faster timescale imposed by a(t) ceases to be present.
8
integration
11
. The system of partial differential equations that commonly governs a physical phe-
nomenon can be discretized in space to reduce the system to a set of ordinary differential equations
in time
12
. If they are explicitly integrated, we can write them as
dy
dt
n+1
=F (y
n
;t
n
)
wherey isanarrayofvariablesthataretobesolvedfor,F istheright-hand-sidefunction(usually
containing advective, diffusive, and source terms), t corresponds to the time, and the subscripts n
andn+1 correspond to current and future time-steps respectively. We can see that the rhs function
is readily available and all that needs to be done is extrapolate the value of the variables to the
next time-step using this known slope. The situation is different in implicit integration. Here,
dy
dt
n+1
=F (y
n+1
;t
n+1
)
Here, the rhs function is unknown, making it rightfully belong to the left hand side. This thus
requires a system of equations to be solved in order to arrive at y
n+1
. Computationally, this is
more intensive, but is also more stable. Implicit integration permits taking large time-steps, i.e. the
time-step size needed to be taken to capture transient behavior accurately is much smaller than the
time-step size required to guarantee stability, despite there being a multitude of timescales in the
problem. The time-step size required to guarantee stability in explicit integration is much smaller
than that required by accuracy. This makes computations in combustion kinetics, in general, very
expensive, due to the solution of large systems of equations using non-linear solvers
13
. There
have been some advancements in feasible explicit techniques in the past decade, but they are not
yet reliable. In any case, as far as the chemistry and diffusion is concerned, reacting flow codes
are in general limited to implicit integration. Thus, optimization comes mostly in the form of
advancements in linear algebra: preconditioners and iterative algorithms
14
.
Even if we do choose implicit integration, we are confronted with problems of efficiency. In a
2-dimensional grid, with a 1000 points in each direction, comprising not just velocities and temper-
ature, but also 50-60 species concentrations, the number of unknowns reaches 10
7
! This means
the storage of the accompanying matrix will need to accommodate roughly 10
14
numbers. This
becomes extremely prohibitive and is only really feasible for 1-dimensional problems. We are thus
forced to adopt a segregated approach [25]: each governing equation is solved for a single variable
over the entire grid while keeping all others constant, subsequently iterated over to minimize the
11
Hirsch [22] and Ascher & Petzold [23] present comprehensive discussions.
12
This approach is called the Method of Lines [22].
13
These require the Newton method, where an iterative search of the solution needs to be performed [20]. These
require Jacobian evaluations which can prove to be expensive.
14
For improvements in chemistry calculations for instance, see Perini [24]
9
difference between the lhs and rhs norms.
Having chosen implicit integration and a segregated approach, we must choose a method of
solution for the fluid dynamics and clarify just which steps are being solved explicitly/implicitly
and how the segregation is carried out.
2.2 Dealing with Incompressibility
Although we are discussing reacting flow problems in which the density changes by a factor of
5 7, a majority of our problems are confined to the so-called Low Mach Number regime
15
, where
thermodynamic pressure is considered a constant in space and time because the speed of sound
is much larger than the rate of any other phenomenon in the problem. This essentially means
that density changes in the problem are solely due to temperature changes; pressure is decoupled
from density. This introduces an inescapable difficulty in numerically solving the Navier-Stokes
equations. Let us consider the case in which density is constant everywhere.
The governing equations (momentum and continuity) are:
@U
@t
+UrU =
1
rp +r
2
U
rU =0
whereU is the velocity, is the density,p is the pressure and is the kinematic viscosity. We see
thatpressurelacksanevolutionequation
16
, implyingthatitacquiresitsfieldvaluesinstantaneously
compared to the flow time-scales
17
. Often we are faced with having an equation and solving for
its variable, but in this case we are presented with a variable (p) with no equation! We can write
one for it by taking the divergence of the momentum conservation equation and substituting the
continuity equation in it to get:
r
2
p =r (UrU)
This Pressure Poisson Equation will have to be solved (instead of the continuity equation),
15
If we consider deflagration waves for example, the fractional change in pressure between the burnt and uburnt
gas is given by
P
Pu
=
M
2
b
b
M
2
u
u
1+M
2
b
b
where the subscriptsu andb denote unburnt and burnt,M is the Mach number,
and
is the adiabatic index. The speed of sound in the burnt gas increases only with the square root of temperature;
M
b
=Mu
q
T
b
Tu
due to thermal expansion.
16
Such a system referred to as a Differential Algebraic Equation [26].
17
Since the pressure is not connected to density via an equation of state, it’s role as a physical variable is now
different. The momentum equation is to be thought of as being constrained by the continuity equation, and pressure
is the appropriate Lagrange multiplier [25].
10
alongsidethemomentumconservationequation, forthesimultaneousvaluesofpressureandvelocity.
Such a method is termed a Projection Method
18
; in this case, the momentum equation is projected
onto the space of velocity fields that are all solenoidal (their divergence is zero). An important
mathematical underpinning of this procedure is the fact that any vector field can be decomposed
into a solenoidal field and a conservative (curl-free) field; this result is codified in the Helmholtz
Decomposition Theorem [29].
This method of approach is still fraught with challenges: there is still no time-derivative for
the pressure, which means we cannot use standard ODE techniques to solve the combined set.
Practitioners of the Finite Element Method solve this set of equations by using consistent pairs of
element types [30], but this results in a large matrix that has to be inverted. This isn’t a problem
today, but given the profusion of finite volume and finite difference codes, formulating the problem
in one of these two ways seemed most natural.
In the finite volume and finite difference domains, another projection approach, is taken. Here,
instead of taking the divergence of the momentum equation, its discrete divergence is taken.
19
The method proceeds as follows
20
:
Discretize the momentum equation implicitly:
U
n+1
i
U
n
i
t
=H
i
(U
n+1
i
)G
i
(p
n+1
i
)
H
i
and G
i
correspond to spatial discretization operators that correspond to the advective &
diffusive terms for velocity and the gradient of pressure respectively. Now, by the Helmholtz
decomposition theorem, there must exist a vector field U
i
that can be decomposed as
U
i
=U
n+1
i
+G
i
(
n+1
i
)t
where is some scalar field, the gradient of which we are free to relate to the gradient of the
pressure field in a number of ways. If we adopt the Kim and Moin relationG
i
(
n+1
i
)t =G
i
(p
n+1
i
),
we have
U
i
U
n
i
t
=H
i
(U
n+1
i
)
18
A first application of such an idea was due to Harlow and Welch [27] and subsequently formalized by Chorin
[28].
19
Historically, the Projection method was first implemented in the finite volume context.
20
There are essentially three main types of projection methods of the kind we’re discussing now: One due to Kim
and Moin [31], another due to Bell, Collela and Glaz [32], and another due to Van Kan [33]. The methods differ
in boundary condition implementation, the type of discretization used, the type of grid used, and as will be made
clear, the degree of projection, which is decided by how much of the pressure gradient to retain while solving the
momentum equation. We adopt the Kim and Moin presentation here only for ease. The arguments that follow does
not rest on this choice.
11
which now lacks the pressure term. The right hand side is still unwieldy, in that it depends on
a future time. We can make the approximation
U
i
U
n
i
t
=H
i
(U
i
)
provided we obey the boundary conditions. Upon solving forU
i
, we can find the the scalar field
by solving
D
i
G
i
(
n+1
i
) =
1
t
D
i
(U
i
)
where D
i
is the discrete divergence operator and has been applied to the Helmholtz decompo-
sition relation. Upon solving for , we can substitute in back into the Helmholtz decomposition
relation to obtain the divergence free U
n+1
i
. Even though we did not specify the spatial discretiza-
tion, such a method stands out in two respects:
1. H
i
is not a linear operator, due to the presence of the non-linear advective term. This means
thatH
i
needs to be evaluated every time is calculated. Some get around this by evaluating
the advective term explicitly while solving the diffusive part implicitly. This introduces other
difficulties.
2. Even though we did not commit to a discretization scheme, upon closer examination, such a
projection method requires the use of a staggered grid
21
to assure consistent coupling between
velocity and pressure, which is not easily implementable for complex geometries. It is also
difficult to extend this procedure to transonic (compressible) flows.
We thus introduce a variant offered by Issa [34]. Issa commits to linearizing the operatorH
i
, making
it depend on past values of velocity. This, which to some may seem as too dear a price to pay,
actually comes with a couple of advantages: the need for staggered grids is not strictly required in
this case, and the method is easily generalized to compressible flows. Issa proceeds as follows:
1. Predictor step: First solve
U
i
U
n
i
t
=H
i
(U
i
)G
i
(p
n
i
)
Note that p is used from the previous time-step. The velocity U
i
need not satisfy the conti-
nuity equation.
2. First corrector step: Consider a new velocity fieldU
i
, which ought to satisfyD
i
U
i
= 0: To
accomplish this, use the predictions from the previous step to write
D
i
G
i
(p
i
) =D
i
H
i
(U
i
) +
D
i
(U
n
i
)
t
21
This means that velocity and pressure are solved on different points in the grid, as opposed to a collocated
arrangement (wherein velocities and pressures are solved for at the same point). This is required to avoid velocity-
pressure decoupling [25].
12
which upon solving yields p
i
which can then be substituted in
U
i
U
n
i
t
=H
i
(U
i
)G
i
(p
i
)
to solve for U
i
.
3. Second corrector step: Consider a newer velocity fieldU
i
, which ought to satisfyD
i
U
i
=
0: To accomplish this, use the predictions from the previous step to write
D
i
G
i
(p
i
) =D
i
H
i
(U
i
) +
D
i
(U
n
i
)
t
which upon solving yields p
i
which can then be substituted in
U
i
U
n
i
t
=H
i
(U
i
)G
i
(p
i
)
to solve for U
i
.
And so on and so forth! We can add as many corrector steps as accuracy requires. In general-
izing this to varying density (and compressible) flows, the velocities will now be weighted with
densities, the Discrete Pressure Poisson Equation will contain an additional term on the left hand
side denoting its coupling to density, and the energy equation needs to be solved iteratively with
pressure and velocity to ensure the desired coupling between thermal-expansion (and pressure work
for compressible flows) and velocity.
This procedure (called the PISO
22
method) is widely adopted in the OpenFOAM
23
suite of
CFD solvers, which is written in a form that applies to abstract tensor calculus [35]. Given the
linearization of operators proposed by Issa, such a framework can handle solving a wide variety of
flow problems. We utilize this very procedure in solving a laminar reacting flow problem (presented
in the next chapter). We must however discuss how the solution procedure for stiff chemistry is
handled alongside that of the fluid mechanics.
2.3 Dealing with the Splitting of Operators
The evolution of scalars, temperature and the concentrations of chemical species, must be coupled
with the fluid mechanics. For this, we can use Operator Splitting
24
, which was already used in the
previous section. The idea behind operator splitting is this: given an ODE that can be written as
dy
dt
= (T +R)y
22
Pressure-Implicit Splitting of Operators
23
Can be found at www.openfoam.org
24
For combustion problems, this is advocated in Oran and Boris [36]. Also see Toro [37].
13
where y is an array of scalars satisfying some initial condition y(0) = y
0
, and T and R are
transport and reaction operators respectively, the solution at a time t after t = 0 is given by
y =e
(T+R)t
y
0
if the operators can be considered/approximated as independent of y. Since we’re keeping with
the habit of linearizing, provided our time-step sizes are small enough to ensure accuracy, we can
attempt splitting the operators naively as follows:
1. Solve the transport part:
dy
T
dt
=Ty
T
y
T
(t) =e
Tt
y
0
2. Solve the reaction step with initial conditions taken from the transport step:
dy
R
dt
=Ry
R
y
R
(t) =e
Rt
y
T
(t)
y
R
(t) =e
Rt
e
Tt
y
0
which can be approximated as the solution of the solution of y at the end of the time-step.
Suchatime-steppingisexactlycorrectifT andRarelinearoperators. Wecanseethisbycomparing
the exact solution
e
(T+R)t
=I + (T +R)t +
T
2
+R
2
2
t
2
+
TR
2
t
2
+
RT
2
t
2
+O(t
3
)
with
e
Rt
e
Tt
=I + (T +R)t +
T
2
+R
2
2
t
2
+RT t
2
+O(t
3
)
which are exactly equal ifR andT commute (possible only ifR andT are linear operators). If this
is not true, we encounter a splitting error proportional to O(t
2
). This ends up being too costly.
Strang [38] suggested a splitting technique that produces second order accuracy:
1. Solve the transport part:
dy
T
dt
=Ty
T
y
T
(t) =e
Tt
y
0
2. Solve the reaction step assuming the initial condition is given by the previous solution evalu-
14
ated at
t
2
:
dy
R
dt
=Ry
R
y
R
=e
Rt
y
t
(
t
2
)
y
R
=e
Rt
e
T
t
2
y
0
3. Solve the transport part again, but this time using the initial condition from the previous
step evaluated at t:
dy
dt
=Ty
y =e
Tt
y
R
(t)
y =e
Tt
e
Rt
e
T
t
2
y
0
If this solution is evaluated at t, it reads as:
y(t) =e
T
t
2
e
Rt
e
T
t
2
y
0
which can be shown to be second-order accurate in time
25
.
The above considerations can be put together, which Cuoci et al. have done in their code
laminarSMOKE [40], which is built upon OpenFOAM’s abstract tensor classes. The scheme is
summarized in Figure 1.
25
An improvement to Strang splitting, where the steady-state errors that occure for large time-step sizes are
removed, can be found in a more recent paper [39]
15
Figure 1: Operator Splitting scheme used in laminarSMOKE (taken from [40]).
3 Directnumericalsimulationsofprobeeffectsinlow-pressure
flame sampling
Data obtained in laminar flames are key towards the development and validation of kinetic models.
Although experimental methods have advanced over the past few decades (e.g.,[41]), there are still
uncertainties associated with standard experiments (e.g., [42]) that could have a major effect on
kinetics given the relatively low sensitivity of flame properties to rate constants [43]. Minimizing
experimental uncertainties requires the in-depth understanding of the controlling physics of the
experimental approach as well as rigorous assessment of the assumptions of the hydrodynamic
model that is used in the numerical simulations.
Molecular beam mass spectrometry (MBMS) [44] is among the leading techniques in measuring
species concentrations in chemically reacting flows. The method mainly comprises probing a low-
pressure (40–100 mbar) burner-stabilized premixed laminar flame with a conical quartz nozzle at
various depths in order to characterize the chemical structure. Having been a staple in combustion
chemistry studies for more than fifty years and serving as an invaluable tool in kinetics studies [19],
the method nonetheless suffers from an inescapable drawback, namely the potential distortion by
the probe of the flame from its one-dimensional adiabatic structure that is assumed in numerical
simulations. This issue of probe-induced effects on flames is not limited to the MBMS technique.
For instance, gas chromatography measurements are just as prone to these distortions [45]. Efforts
16
in the past ([46], [47],and [48]) have focused on quantifying these systematic errors to be able to
“subtract” them out. Recently, Struckmeier et al. [49] published the results of a comprehensive
experimental investigation on the effects of probe angle, orifice diameter, and other experimental
parameters of the flame structure. Skovorodko et al. [50] developed a code to model the entire
experimental domain including the expansion section in the nozzle. However, this was done using
single-step chemistry and by spatially fixing the source terms in the energy and species conservation
equations.
A systematic and rigorous computational study of probe effects in low-pressure flame sampling
that delineates the effects of all the relevant boundary conditions while incorporating detailed
chemistry has yet to appear. One possible reason for this is the lack of freely available, open-
source numerical codes that could be modified to describe the exact laboratory conditions and used
by experimentalists. Recently, Cuoci et al.[40] have advanced the laminarSMOKE finite-volume
based code built using the OpenFOAM
26
suite of Computational Fluid Dynamics tools and the
OpenSMOKE [40] library of functions that handles detailed chemistry and transport while offering
an interface to various Ordinary Differential Equation solvers. The code also incorporates a robust
algorithm for time stepping, making it amenable to the stiff chemistry normally encountered in
combustion applications.
The goal of this investigation was to perform for the first time direct numerical simulations of
the probe effects on the structure of a low-pressure laminar flame and to assess the deviations from
results obtained in unperturbed flames; it should be noted that in direct numerical simulations all
pertinent spatial and temporal scales are resolved without any simplifying assumption in any type
of flow laminar or turbulent alike. Towards that goal, the laminarSMOKE [40] code was adopted to
conditions that are close to the experimental ones. The simulations are considered as direct given
that the fluid mechanics, chemical kinetics, and heat and mass transport processes were modeled
based on the current state-of-the-art knowledge and without invoking any simplifying assumptions.
This approach will lay the fundament for further study with the aim to provide correction rules
for existing data and result eventually in an optimized probe design that would cause a minimum
disturbance to the flame.
3.1 Modeling Approach
Numerical simulations of the three-dimensional problem coupling fluid dynamics and full chemistry
in this study are based on the laminarSMOKE code which addresses the stiff chemistry [21] compo-
nent of reacting flows by employing the Strang splitting [38], a procedure founded on the premise
that for suitably small time steps the processes of molecular transport and reaction can be separated
into a sequence of time integrations wherein the final solution of the transport step is the initial
26
Can be found at http://www.openfoam.org
17
Figure 2: Schematic of the axisymmetric computational domain. The orifice is placed at a constant
perpendicular distance of 5 cm from the flange so that the domain height is 6 cm when the orifice
is at 1 cm from the burner exit and 5.4 cm when the orifice is at 0.6 cm from the burner exit.
condition for the reaction step. The reaction step is thus reduced to a spatial array of homogeneous
reactors with appropriate initial conditions whose time integration is handled by standard ODE
solvers such as DVODE[51] that was used in the present study. The Strang splitting layer is super-
imposed on the well-established PISO algorithm (Pressure Implicit with Splitting of Operators)[34]
for the pressure-velocity coupling, to complete the algorithm for solving the discretized conservation
equations [52].
The study was performed for a 50 mbar propene/oxygen/argon mixture at an equivalence ratio
= 2:3 with mole fractions 0.255/0.495/0.25, respectively, a C/O ratio of 0.773, and a velocity
of 48 cm/s at the burner exit that is the inlet of the solution domain; matching the experiments
in [48]. The USC Mech-II kinetic model
27
was reduced from 117 species to 34 species using the
directed relation graph (DRG) method [53] to save computational cost while preserving the essential
aspects of the flame structure. Thermal radiation and Soret diffusion were not included, and the
thermodynamic, transport, and kinetic properties were evaluated using the OpenSMOKE library.
The Soret diffusion formulation is not available in the current version of laminarSMOKE. The
code will be modified to account for it in subsequent studies in which the aim will be to perform
quantitative comparisons against available experimental results.
The axisymmetric geometry of the solution domain as seen in Fig. 2 was enclosed by the surfaces
whose boundary conditions are defined in Table 1 for all cases considered. These cases were chosen
in order to isolate the various effects related to geometry, heat loss, and compressibility.
The problem was first solved on a series of coarse grids and subsequently mapped on to grids
of increasing cell density. In the last level of refinement, the number of cells in the incompressible
cases when the orifice is closed (UP, AFC, NAFC, AC, NAC, defined in Table 1) were approxi-
27
Can be found at http://ignis.usc.edu/USC_Mech_II.htm
18
Cases Orifice Probe Flange Burner exit/inlet Co-flow Outlet
Unperturbed
(UP)
N/A N/A N/A rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Adiabatic
Flange only;
Closed Orifice
(AFC)
N/A N/A rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Non-Adiabatic
Flange only;
Closed Orifice
(NAFC)
N/A N/A rp = 0;
T = 600K;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Adiabatic
Flange and
Probe; Closed
Orifice (AC)
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Non-Adiabatic
Flange and
Probe; Closed
Orifice (NAC)
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
T = 1600K;
rY
k
= 0;
Un = 0
rp = 0;
T = 600K;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Adiabatic
Flange and
Probe; Open
Orifice (AO)
p =
20mbar;
rT = 0;
rY
k
= 0;
rU = 0
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
rT = 0;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Non-Adiabatic
Flange and
Probe; Open
Orifice (NAO)
p =
20mbar;
rT = 0;
rY
k
= 0;
rU = 0
rp = 0;
T = 1600K;
rY
k
= 0;
Un = 0
rp = 0;
T = 600;
rY
k
= 0;
Un = 0
rp = 0;
T = 303K;
Y
k
+
D
k
U
rY
k
=e
k
;
Un = 48cm=s
rp = 0;
T = 303K;
Y
Ar
= 1:0;
Un = 48cm=s
p =
50mbar;
rT = 0;
rY
k
= 0;
rU = 0
Table 1: Boundary Conditions employed for all cases considered. e
k
is the flux fraction of the
species specified upstream of the burner exit. The velocity normal to the surface is denoted by U
n
.
19
mately 20,000 skewed by a geometric progression to achieve sub-micrometer cell sizes close to the
orifice and the reaction zones. For the compressible cases (AO, NAO, defined in Table 1) when
the orifice is open, approximately 80,000 cells were used to ensure reasonably mesh-independent
solutions in which species concentrations changed by less than 1% upon doubling the number of
cells. Interpolation schemes used to evaluate variables at the control volume faces were linear.
Implicit Euler time stepping was used to arrive at the steady-state solution. The geometric alge-
braic multi-grid (GAMG) solver was used for the pressure Poisson equation (PPE), preconditioned
with diagonal incomplete Cholesky factorization (DIC) whereas the compressible cases required the
diagonal incomplete LU (DILU) factorization owing to the asymmetry of the coefficient matrix in
PPE. For all other variables, the biconjugate gradient (BiCG) method with DILU preconditioning
was used. These linear algebra routines [54] and differencing schemes are all part of the Open-
FOAM toolbox. All grids were generated using Gmsh
28
. The modified PISO algorithm used in
OpenFOAM accompanied by an analysis of the discretization errors can be found in [55].
Neumann-type pressure boundary conditions were specified everywhere except at the outlet of
the domain in the UP, AFC, AC, NAFC, and NAC cases. The AO and NAO cases were simulated
using the Dirichlet-type boundary conditions to achieve transonic flow at the vicinity of the orifice.
For high Reynolds number flows, specifying the gradient of pressure, grad(p), as zero is a justified
approximation at the probe and flange owing to boundary layer theory [56]. At the main and
co-flow inlets,rp is not strictly zero and is proportional to the second spatial derivative of velocity
[57], but this coupling was neglected in the present study in order to separate out the complexities
that arise from such feedback effects.
The inlet boundary condition for species mass fractions was of the Robin type, wherein the
convective and diffusive flux contributions are required to add up to the net supply of species
specified upstream[58]. The temperature boundary conditions at the probe and flange were of the
Dirichlet type, thus asserting the temperature instead of solving for them via Robin boundary
conditions that require as input the heat transfer coefficient and an estimate of the amount of
heat transferred. This approximation is justified, as among the goals of the present study is to
characterize the effect of heat loss from the flame.
3.2 Results and Discussion
3.2.1 Unperturbed case (UP)
This reference case is void of probe and flange effects. The flame is attached to the burner with
the OH profile peaking at around 0.5 cm from the burner exit. The axial temperature profile of the
unperturbed case is slightly lower compared to the solution obtained using the PREMIX code[59],
due to radial temperature gradients caused by the cold argon co-flow as shown in Fig. 3. As will be
28
Can be found at http://www.geuz.org/gmsh
20
Figure 3: Streamlines superimposed on temperature contours for the unperturbed case (left) com-
pared to a perturbed case (right) with heat loss to the probe and flange along with suction through
the orifice.
shown in the following cases, the radial variations of species and temperature are highly coupled and
perturbations grow in severity in the presence of the probe as evident in Fig. 4, calling into question
the standard approach of simulating such flames by solving the species conservation equation for a
given measured temperature profile in one dimension.
3.2.2 Adiabatic Flange Only; Closed Orifice Case (AFC)
In the experiments, the probe and flange are moved towards and away from the flame as one unit.
To isolate the various effects, first simulations were carried out with an adiabatic flange only present
at 6 cm and 5.4 cm from the burner exit. These correspond to orifice locations of 1 cm (of which all
correspondingplotscanbefoundinSupplementarymaterial)and0.6cm, whentheprobeispresent.
A key difference in the AFC case compared to the UP case is the no-slip boundary condition in
the former. For the classic stagnation premixed flame, Egolfopoulos et al. [60] have shown that
the no-slip boundary condition decreases the strain rate throughout the flame, thereby increasing
the residence time of chemical species and correspondingly the maximum temperature with respect
to the UP case. This reasoning equally applies to the current burner-stabilized flames and was
confirmed in these simulations as shown in Fig. 5.
Thisstagnationflowconfigurationistheunderpinningforalargeportionoflaminarflamestudies
due to its simple quasi one-dimensional character. Flames studied in this setting are simulated
typically via the OPPDIFF code [61]. If the probe were thought to collapse onto the flange while its
orificeremainsonthecentralaxis, thedifficulttwo-dimensionalproblemof“probe-inducedeffects” is
21
Figure 4: OH mass fraction contours for the unperturbed case (left) compared to a perturbed case
(right) with heat loss to the probe and flange and suction through the orifice.
Figure 5: Centerline temperature profiles when the flange (in the absence of the probe) is placed
at 6 cm (left) and 5.4 cm (right) from the burner exit.
22
reducedtoareasonablyaccurateone-dimensionalproblemof“flange-inducedeffects”, thusamenable
to modeling via a slightly modified version of OPPDIFF, provided the open orifice effects are shown
via experiment to be negligible. This novel departure is the basis of the recent burner-stabilized
stagnation (BSS) flame approach used to measure soot particle size distributions in premixed flames
[62]. The flange-induced perturbations documented by that study agree qualitatively with the
present ones.
3.2.3 Non-Adiabatic Flange Only; Closed Orifice Case (NAFC)
As shown in Table 1, the flange temperature is set to 600 K, and as long it is sufficiently far away
from the flame, the heat loss does not perturb the flame significantly. The results of Fig. 5 indicate
that the maximum temperature slightly increases with respect to the UP case, similar to the BSS
flame studies.
3.2.4 Adiabatic Flange and Probe; Closed Orifice Case (AC)
Superposing an adiabatic probe with a closed orifice on the above configuration decreases the
temperature at the sampling location by about 40 K, as shown in Fig. 6. The OH concentration
has not reached its equilibrium value when the orifice is at 1 cm from the burner exit, as shown
in Fig. 7. Therefore, its consumption via CO +OH! CO
2
+H decreases due to the decreased
residence time caused by the presence of a stagnation surface. However, the low probe cone angle
of 20° mitigates this effect by introducing only a mild change in the direction of the streamlines
and hence even when the orifice is located at 0.6 cm from the burner exit, the perturbation is
relatively minor. Greater probe angles are expected to decrease the temperature further because of
the increased effective stagnation surface impinged upon by fluid upstream, and when conduction
losses are considered the situation could worsen. On the other hand, narrow probe angles cannot
be used because in the supersonic expansion section of the MBMS probe, satisfactory freezing of
the reactions is achieved only for wider probe angles [46].
3.2.5 Non-Adiabatic Flange and Probe; Closed Orifice Case (NAC)
Relaxing the assumption of adiabaticity, the probe temperature (1600 K) was chosen to mimic
the softening temperature of quartz, which is the maximum possible probe temperature attainable
in the actual experiment, thus introducing heat losses and decreasing the reactivity of the flame,
especially around the orifice where the temperature is the lowest. This effect is expected to worsen
when the orifice is brought closer to the burner exit, or the probe angle is increased as a larger
portion of the flame loses heat to the probe.
When the conduction losses to the probe become comparable to the conduction losses to the
burner surface, the flame can attach itself to the probe, a phenomenon observed during measure-
23
Figure 6: Centerline temperature profile when the orifice is at 0.6 cm from the burner exit.
Figure 7: Centerline OH mass fraction profile when the orifice is at 0.6 cm from the burner exit.
24
ments for large (around 65°) probe angles[49]. This tendency to attach to the probe can be seen
in the OH mass fraction contours for the perturbed flame in Fig. 4. Many chemical reactions are
thus left incomplete or quenched. The local quenching cascades upstream, especially via species
with large destruction time-scales, perturbing the entire chemical structure of the flame. Consider
for example, the consumption of the highly active ketynyl radical (HCCO) that depends on the
HCCO +H ! CH
2
+CO reaction. The temperature drop close to the probe reduces the H
concentration, thereby reducing the rate of consumption of HCCO and subsequently increasing the
concentration of CO. This elevated CO concentration at the orifice tends to increase the diffusive
flux upstream. At the same time, the Robin boundary condition imposed at the burner exit bal-
ances any diffusive flux with an equal and opposite convective flux, thus producing an increased
CO mass fraction at the burner exit.
3.2.6 Adiabatic Open Orifice Case (AO)
Before allowing for suction through the orifice for the NAC case, it is worthwhile to observe how the
temperature field changes solely due to suction in the absence of conduction losses. The dramatic
drop in temperature close to the probe shown in Fig. 6 is due to the rapid conversion of sensible
enthalpytokineticenergy. Thefirstorderextentofinfluenceinthetemperaturefieldisproportional
to the diameter of the orifice unlike the NAC case. The chemical reactions are nearly frozen in this
region causing the similar cascade effect mentioned above.
3.2.7 Non-Adiabatic Open Orifice Case (NAO)
This corresponds closely to the realistic case in which conduction losses and suction are simulta-
neously present. The thickness of the momentum boundary layer in the vicinity of the orifice is
decreased due to the increased local velocity. Since the Prandtl number is close to one, the thermal
boundary layer is also diminished by roughly the same extent, tending to reduce the overall heat
transfer to the probe.
When the orifice is at 1 cm from the burner exit, the temperature upstream of the compressible
flow region falls more or less between the NAC and the AO cases. Analysis showed that the
gas density, , at the tip is highest compared to both these cases. This results in a particularly
interesting change in the kinetic pathway of the formyl radical (HCO), which is mainly consumed
via HCO +M!H +CO +M , whose reverse rate is favored at high densities resulting thus in
decreased amounts of CO and H near the tip and increased amounts of HCO. Although the same
phenomenon is observed when the orifice is 0.6 cm from the burner exit (Fig. 8), the overwhelming
effects of increased strain rate and heat loss mentioned above nearly quench the central core of the
flame, causing the distorted OH field seen in Fig. 4.
Such effects have major influences on soot production as well. The propargyl radical (C
3
H
3
), a
25
Figure 8: Centerline HCO mass fraction profile when the orifice is at 0.6 cm from the burner exit.
major precursor to soot formation, is notably affected, as seen in Fig. 9. Similar effects were seen
for the acetylene (C
2
H
2
) and ethylene (C
2
H
4
) concentration profiles.
3.3 Concluding Remarks
The effects of probe and the supporting flange used in low-pressure flame sampling were investi-
gated through direct numerical simulations and the use of detailed descriptions of chemical kinetics
and molecular transport. A rich 50 mbar propene/oxygen/argon flame was modeled using lami-
narSMOKE, a fully compressible finite-volume open-source code. Given that there is vast variation
of scales in this problem resulting from both kinetics and compressibility, notable effort was ded-
icated to the various numerical algorithms involved as well as towards achieving sufficient spatial
resolution and solution convergence. Approximately 20,000 and 80,000 cells were used for the in-
compressible and compressible cases considered, which correspond to closed and open sampling
orifice, respectively.
Even under adiabatic conditions, it was determined that the presence of the flange with or with-
out the conical probce has a notable effect on the flame structure that becomes two-dimensional.
Introducing heat loss and compressibility by allowing for suction through the sampling orifice re-
sulted in further deviations of the flame structure compared to that of a one-dimensional unper-
turbed flame.
The results of this first study of its kind will guide future investigations towards developing
correction rules for existing data and eventually optimizing the probe design in order to minimize
flame structure perturbations.
26
Figure 9: Centerline C3H3 mass fraction profile when the orifice is at 0.6 cm from the burner exit.
4 Coupling Combustion and Non-thermal Plasma Kinetics
The following is a brief introduction to the concepts and vocabulary used in describing discharge
phenomena.
4.1 Electron Temperature and Relaxation Time
Fowler[7] first pointed out that the definition of Temperature needed to underlie the edifice of Ther-
modynamics in the justice of logical consistency and thus proposed the Zeroth law. The attempt to
explain the laws of Thermodynamics is taken up by Statistical Mechanics[6], wherein the quantity
Temperature is introduced in Gibbs’ canonical ensemble as the derivative of energy with respect to
entropy (the number of microstates; Boltzmann entropy) at constant volume, at equilibrium. As
functional as this definition is, it doesn’t shed light on how the temperature is partitioned when
one possesses more information about the system, for instance, the average vibrational energy of
a collection of particles (measurable via spectroscopy) or the average kinetic energy of a swarm
of electrons (measurable with a Langmuir probe). A concrete conception of Temperature can be
attained when one views the problem through the lens of Information Theory
29
: The trajectories
of a large number of particles are futile to compute; a large number of degrees of freedom ought to
be reduced to a few macroscopic variables (like Temperature, Pressure, Concentration, etc.); these
variables are parameters of probability distributions; the probability distributions themselves are
29
It was Jaynes [63] who unified Gibbs’s [64] ideas of phase space and Shannon’s Information Theory [65]. This
invaluablecontributiontostatisticalmechanicshasnotonlyclarifiedtheinterpretationofthermodynamicsquantities,
but has also inspired progress in non-equilibrium thermodynamics [66].
27
arrived at by maximizing uncertainty/ignorance/entropy subject to given constraints (energy and
its partitioning, volume, total number of atoms of particular kinds, etc.); the Lagrange multipliers
corresponding to each of these constraints are directly related to the sought after macroscopic vari-
ables, one of which, the temperature, is directly related to the Lagrange multiplier corresponding to
the measurable energy. This generalizes to any partitioning of energy; if the energies of particular
degrees of freedom are specified, then there exists a Lagrange multiplier for each partition and sub-
sequently, there exists a Temperature for each of these Lagrange multipliers. Therefore, under the
assumption of local and partial equilibrium, rotational temperature quantifies the internal energy
associated with the rotational degree of freedom, vibrational temperature quantifies the internal
energy associated with the vibrational degree of freedom and electron temperature quantifies the
internal energy.
In ideal gases, thermal equilibration occurs due to elastic collisions
30
between particles. The
timescale of this relaxation,
nn
is approximately equal to
1
nv
where n is the number of particles
per unit volume ( 10
26
m
3
), v is the average velocity of a particle ( 10
2
ms
1
) and is the
collision
31
cross-section
32
( 10
20
m
2
) . For standard conditions of temperature and pressure,
nn
t 10
8
s , the subscript nn indicates that the collision is between two neutral molecules (the
need for this distinction will soon become clear). This relaxation time-scale must not be confused
with the collision time-scale, which is the time spent by a particle during a collision. This is
approximately equal to
d
v
; where d is a typical length scale of interaction (for short range forces,
this is roughly 10
10
m, the square of which is roughly the collision cross-section). The collision
time scale ( 10
12
s) is thus smaller than the relaxation timescale by a factor ofnd
3
: When charged
particles are present in the gas, the Coulomb forces between them, unlike the short range Lennard-
Jones or Van der Waals forces between neutrals, does not present itself with a length scale in the
same sense of d, due to the relatively slow inverse square fall-off. The closest encounter between
charged particles occurs when their potential and kinetic energies are comparable;
Z
2
e
2
r0
t
3
2
kT,
where Z is the charge of the particle, e is value of a unit charge and r
0
is the encounter distance.
So the collision cross-section in such a case is given by
c
tr
2
0
t
4
9
Z
2
e
4
(kT )
2
Such collisions can be classified as strong
33
ones and are typically less frequent owing to the
30
Inelastic collisions engender chemical reactions and are considered to occur less frequently, which is another way
of saying that the local temperature has equilibrated faster than the species concentrations. Molecular energies follow
the Maxwell distribution, whose tail comprises the molecules that undergo chemical reaction.
31
By “collision” we mean an interaction between particles in which the energy transfer is capable of significantly
redirecting one of the particles.
32
The concept of a collision cross-section serves as a very useful abstraction for not just elastic collisions, but for
inelastic ones as well. The ratio of the number of particles in a particular state to the number of incident particles
per unit area is the cross-section for the particle to be in that particular state [67].
33
Actually, we have ignored the possibility of quantum mechanical interactions. This is a reasonable assumption
28
long range nature of Coulomb forces; these more frequent collisions will henceforth be regarded as
weak. Calculating the relaxation time for such charged particles requires taking into account all
such encounters for a large number of particles
34
. The velocity change v of a charged particle
flying by another is roughly
Ft
m
Z
2
e
2
mvr
. Since v can be either positive or negative, we can consider
the change of its square, (v)
2
Z
4
e
4
m
2
v
2
r
2
. The rate of change of this quantity for a flux of particles
is thus
d(v)
2
dt
nv
(v)
2
2rdr
2nvZ
4
e
4
m
2
v
2
dr
r
The lower and upper limits of this integral correspond to the closest (strongest) and farthest
(weakest) encounter respectively. The latter is given by the Debye radiusd, which is the length over
which the electric field of a charged particle has decayed by 63%(i.e 1
1
e
) due to the neutralizing
effect of surrounding oppositely charged particles[68]. The closest encounter is given by . Therefore,
dr
r
t
d
r0
dr
r
=ln(
d
r
0
) =ln()
=
3(kT )
3
2
2(4)
1
2
Z
2
e
3
n
1
2
We can now approximate the relaxation time-scale as
1
=
1
v
2
dv
2
dt
tn v
2
9
Z
4
e
4
(kT )
2
ln()
This applies for ion-ion collisions and electron-electron collisions and can be labeled
ii
and
ee
respectively. The collisions between electrons and ions (or neutrals) requires considering the great
discrepancy in their masses. If we consider only weak collisions, which are more frequent anyway,
the energy transfer between an electron and a proton is given by
E =
p
2
2m
i
t
m
e
m
i
v
e
E
where m
e
and m
i
are the masses of an electron and ion respectively, E =
1
2
m
e
v
2
e
is the kinetic
energy of the electron,is the angle of deflection (which is less than a radian for weak collisions). It
is clear that the greater the discrepancy in masses, the lesser is the energy transfer. Equivalently, it
takes many more collisions for electron-ion equilibration than electron-electron and ion-ion equlibra-
provided the time-averaged distance between electrons always exceeds the de-Broglie wavelength,h=p, each of whose
kinetic energy is kTe. This also implies that the electrons behave according to Boltzmann statistics and not Fermi-
Dirac statistics[6] (which electron gases exhibit in metals).
34
Zeldovich’s and Raizer’s book[11] are the source for most of the arguments that follow.
29
tion
35
. It is for this reason, that electrons in a plasma can be considered to possess a temperature
that is different from the rest of the gas.
The time for relaxation between electron and bath gas can be extended by imparting additional
energytotheelectronsviaanexternalelectricfield. Inthisway, theelectronscovergreaterdistances
without suffering elastic collisions with the bath gas while at the same time, the collisions that do
occur are of sufficient energy to excite or ionize the bath gas.
4.2 Discharge Physics
A gas discharge comprises a swarm of electrons originating in ionization and/or an external
source in the presence of an electric field. Various types of discharges exist, each with its own
characteristic electric field, electron density, and relaxation dynamics. As much as this remains
beyond the scope of this work, we will briefly go over the general ideas.
A simplified Langevin type kinetic equation for the motion of an electron
36
is
m
dv
dt
=eEm
m
m
=
c
(1cos)
Here, m is the mass of the electron, v its velocity, e its charge, Ethe ambient electric field,
c
the collision frequency ,nv
c
,
c
is the cross-section of elastic collisions
37
. The factor cos is to be
thought of as an average scattering angle, when it is equal to zero, we have isotropic scattering,
when it is 1, the scattering is mostly forward (in the direction of the field) and when it is -1, the
scattering is backward. Upon solving the above equation with an initial condition for velocity, we
note that the steady state solution gives usv
d
=
eE
mm
which is the so called drift velocity. Multiplied
by the charge density, this gives the current density in the field. If we write the current density as
j = E, is the conductivity of the ionized gas. It is immediately obvious that the conductivity
is inversely proportional to the collision frequency; this kinetic energy of the swarm of electrons is
dissipated in the gas as heat. The energy released by this current per second per unit volume is
given by the product jE. The kind of discharge we are interested in studying is one in which this
Joule heating is minimal, causing a temperature rise of about 10K. Our objective is to excite and
35
And for that matter ion-neutral and neutral-neutral equilibration. The relaxation time for electron-ion equi-
libration is approximately that of the ion-ion equilibration multiplied by the factor
m
i
me
. The temperature in this
expression is closer to that of the electrons than the ions since the contribution to the relative velocity between the
two particles is almost exclusively the electron’s.
36
This section is a summary of key arguments presented in Raizer[69].
37
We will be working with weakly ionized plasmas, where the degree of ionization
ne
n
<<1. If this is not true, the
Coulomb cross-section must be included in the calculation of collision frequency.
30
ionize molecules without heating them; i.e. to achieve a large frequency for inelastic rather than
elastic collisions.
If we equate the energy lost by an electron to an elastic collision with that gained by the electric
field, we find that the mean electron energy is roughly proportional to the reduced electric field
(E/N) and inversely to the ratio of electron mass to neutral mass. This is qualitatively correct. If
we are to take into account the distribution of electron energies within the swarm, we must solve a
detailed kinetic transport equation. We will discuss this in a later section.
In order to create an electron swarm in a gas, it must first be electrically broken down. The ratio
of the frequency of ionization
i
to the drift velocity v
d
is what is termed the Townsend ionization
coefficient . It is empirically given by
=Ape
Bp
E
It is directly proportional to the pressure since the fraction of ionizing collisions increases pro-
portionally. The exponential factor is a Boltzmann weighting; the fraction of electrons capable of
ionizing an atom is proportional toe
I=eEl
, wherel is the distance an electron travels without colli-
sion i.e. the mean free path, which is inversely proportional to pressure. If an anode and a cathode
are spearated by a distance d, and say the cathode is irradiated with a UV lamp and produces
photocurrent i
0
, then the current at the anode is given by i =i
0
e
d
. The number of positive ions
generated in the gap by an electron leaving the cathode is e
d
1. Therefore, the current at the
cathode isi
0
+i
0
(e
d
1). At the vicinity of the cathode, an electron produced by this mechanism
is capable of knocking out additional electrons from the cathode. The electronic part of this current
is thusi
1
=i
0
+
i
1
(e
d
1), where
is the number of electrons knocked out of the cathode per unit
electron present in the discharge. The current at the anode is thusi =i
1
e
d
=i
0
e
d
=[1
(e
d
1)].
Electrical breakdown is said to be achieved when the denominator of this expression goes to zero.
The voltage at which the breakdown occurs, when plotted as a function of the product of pressure
and electrode distance, exhibits a minimum. The left of this minimum, for which the pressure or
electrode distance (or both) are very small, the voltage required is high in order to compensate for
losses and lack of collisions. On the right, the rate of loss due to recombination and attachment are
appreciable. This curve is called the Paschen curve and is an important assessment in discharge
technology.
To obtain a non-equlibrium plasma discharge that excites molecular states, this resulting current
from breakdown must be kept to a minimum. Many techniques involve rapidly reducing the electric
field upon breakdown, causing inhomogeneous fields (glow discharges), which can be maintained
by flowing gas at a fixed flow rate between the electrodes. There is however another mechanism by
which electrical breakdown can take place: the so called streamer discharge. When the pressure and
/or electrode gap is high, and the voltage exceeds the breakdown voltage, the avalanche multipli-
31
cation by the aforementioned processes are too slow to produce the conducting channels that have
been experimentally observed. The secondary emission process does not seem to take place either,
as manifested in the independence of the formation of the conducting channel from the material
of the electrode. In a streamer discharge, an intense electron avalanche leaves behind a trail of
positively charged ions. The field caused by this separation causes secondary electron avalanches.
Photons that are emitted from these cause further photoionization. Once this resulting channel
reaches the anode, the completed circuit conducts a large amount of current and degenerates to an
arc discharge.
The electron density in the swarm that is headed towards the anode in a streamer discharge is
given by:
n
e
= (4D
e
t)
3=2
exp[
(xv
d
t)
2
+r
2
4D
e
t
+ (a)v
d
t]
whereD
e
is the diffusion coefficient of an electron, t is time, x is the spatial coordinate, r is the
center of the swarm, a is the attachment coefficient. The characteristic diffusion radius
38
is given
by
r
D
=
p
4D
e
t =
s
4D
e
x
0
e
E
0
=
r
8x
0
3eE
0
where
e
is the electron mobility, is the mean energy of the electron, and the subscript 0 refers
to values at the center of the cloud. However, the ions accumulate at each point and their total
density can be given by:
n
+
=
t
v
d
n
e
dt
0
For the density of negative ions, we may replace with a. This integral, in the limit of t!1,
we have
n
+
=
[r
D
(x)]
2
expfx
r
2
[r
D
(x)]
2
g
which says that the ion density increases away from the cathode and decreases away from the
central axis of the swarm. It is possible to estimate the current due to the motion of an electron
between two electrodes: if during its motion, an electron gains an energye(V=d)dx, wheredx =vdt
38
Using the diffusion coefficient of an electron here is justified provided the electron is far enough away from the
positively charged ions. Otherwise, the electrons are slowed down due to the attraction from the positive ions. The
positive ions experience a reactive force, resulting in the ions and electrons drifting together. This phenomenon is
called ambipolar diffusion and modifies the diffusion coefficient. On the other hand, if the electrons are too close
together, their mutual repulsion is what causes their motion, in contrast to diffusion. The argument here applies to
timescales after which this repulsion has taken place.
32
,V is the voltage, and d is the electrode spacing, and the work done by the power supply is iVdt;
then the current is given by i =ev=d.
39
For non-uniform discharges, the d needs to be replaced by
E=V, which is proportional only to the geometry of the electrode configuration and not the value
of the electric field itself.
The separation of charges creates an internal electric field. If this field is to make a difference
to the local electron density, it must be of the same order of magnitude as the external field. The
field due to the charges can be given by E
0
=eN
e
=R
2
, where R can be taken as the characteristic
diffusion radius of the electrons. If we equate this with the external field, knowing that N
e
=e
x
,
we can estimate a value for x that can be considered a criterion for streamer formation (Meek
criterion[69]).
Once the streamer completes the circuit, the highly conductive channel carries large amounts
of current. Designing streamer discharge systems that keep the current in this final phase to a
minimum is difficult but has largely advanced over the last two decades[70]. A simpler way to
attain such low currents is by using a corona discharge. Breakdown however is a function of not
only electric field distribution, but also of electrode material. Nevertheless, these have been used
to promising effect recently[71].
4.3 Plasma Kinetics
For the discharges we are interested in, the elastic scattering cross-sections of the electrons (the
process that is responsible for heating the gas) is typically orders of magnitude smaller than the
scatteringcross-sectionofinelasticcollisions. Ionizationsandexcitationscanbeconsideredthemost
crucial collisions as far as generating the required radical concentrations to initiate combustion
40
.
We will discuss these at an elementary level.
The ionization of an unexcited atom follows from its collision with an electron possessing suffi-
cient energyI. If it is possible to find an electron velocity distribution (which can be written as an
electron energy distribution) f
e
(v)dv, and given the the ionization cross-section
e
(v), the number
of ionization events per unit volume per unit time is given by
Z
e
ion
=n
a
n
e
1
v
k
e
(v)vf
e
(v)dv =n
a
n
e
e
where n
a
and n
e
are atom and electron densities respectively, and m
e
v
2
e
=2 =I. The integral
e
is contributed to only by the fraction of electrons which possess the required energy of ionization.
If the recombination rate is denoted by
e
, the total rate of production of electrons can is given by
39
One can find the the difference in the ionization and attachment coefficient by measuring the current since the
right hand side will now contain the electron density.
40
Actually, the radicals are produced as an indirect result of these processes; excited species can either de-excite
to radicals through collisions with neutrals, or ionized species can recombine with electrons to form radicals.
33
dn
e
dt
=
e
n
a
n
e
e
n
+
n
2
e
The recombination rate can be found from the principle of detailed balancing:
e
=
e
K(T
e
)
where the equilibrium constant, evaluated at the electron temperature, is given by the Saha
equation
41
:
K(T
e
) =
g
+
g
a
2(2m
e
kT
e
)
3=2
h
3
e
I=kTe
whereg
+
andg
a
are the degeneracies of the ion and atom, i.e. the number of states they exhibit
for a given energy andh is the Planck constant. As discussed in the previous section, the timescale
for electron avalanche formation can be given by
1
ena
. The electron energy distribution function
that occurs in the integral
e
is not always Maxwellian because of the large degree of disequilibrium
between the electrons and the surrounding atoms. This energy distribution must therefore be a
result of a computation that tracks the detailed transport of electrons through energy space. The
equation that describes this motion is the Boltzmann equation, which will be described in the next
section. The scattering cross-section is also an input, usually provided from experiments. It is
known to typically have a maximum, which can be explained by a formula derivable from classical
mechanics[67]; if the energy transfer between a particle of energy to a “target” particle is to exceed
the energy E, the the cross-section for such a transfer is given by
=
e
4
(
1
E
1
)
This formula gives the qualitiatively correct trend for the behavior of as a function of electron
energy; it typically exhibits a maximum (usually at 10
16
cm
2
) and is fairly linear nearE; a detailed
quantum mechanical explanation can be obtained from the Born approximation[11]. A large part
of the numerical solution of the electron energy distribution as a function of time requires frequent
lookup in a database comprising these cross-sections.
Similar reasoning applies in deriving the rates of excitation due to electronic collisions and
won’t be repeated here. For practical purposes, the distinction between many excited states, like
certain vibrational and rotational states, is unimportant for the purposes of estimating final radical
41
The appearance of Te here indicates that the ion and atom are at the same temperature. This is not strictly
true in a highly non-equilibrium setting, where the electron and gas temperatures vary significantly. Solving for the
electron temperature variation with time requires the Boltzmann equation.
34
concentrations upon discharge
42
. These can be “lumped”, i.e. their cross-sections for formation can
be suitably averaged over. Details can be found in[11].
4.4 The Boltzmann Transport Equation
TheBoltzmannTransportEquation(BTE),describestherateofchangeofasingle-particle-velocity-
distribution in position and velocity space[9, 6]. It is considered one of the triumphs of kinetic
theory, but has much left to be desired in its foundations
43
. We will assume that it describes the
energy transfer processes of electrons reasonably well for our purposes of estimating final radical
concentrations. The BTE resembles a typical transport partial differential equation, with the
exception of its source term which makes it an integro-differential equation. If f(v;x;t)is the
single-particle-velocity-distribution of an electron in the presence of an electric field E;the BTE
reads as:
@f
@t
+vrf
e
m
Er
v
f =C[f]
The left hand side describes the total time derivative represented as a sum of local time variation
andspatialvariationviaadvectioninpositionandvelocityspace. Therighthandsideisthecollision
integral that accounts for all possible collision processes (elastic and inelastic) and is to be thought
of as a functional rather than a function of f.
44
The integral contains the velocity distribution
weighted by the relative velocities between particles and the collision cross-section of events that
serve to change the energy distribution of the electron swarm; these are the ionization, excitation,
attachment and momentum cross-sections.
42
Such an argument needs to however be checked with such states considered in a test simulation. Molecular
nitrogen for instance possesses vibrational states that last as long as the time-scales of combustion. This added
complexity was ignored in our hydrocarbon ignition study (where argon was used), but was retained for our NO
remediation study.
43
The foundations of kinetic theory continue to remain controversial. This is directly related to a deep divide
that exists in the interpretation of probability; one camp believes that the probability of a certain event occurring
is a property of the system being examined, while the other believes that it is a function of the information about
the system one possess[72]. This second camp, in this author’s opinion, has articulated, clarified, and advanced
thermodynamics[73][74][75] more fruitfully than the first (see for example Keck[66]). The BTE is an equation that
describes the time evolution of a probability distribution, which arrives at a steady state, Maxwellian velocity
distribution. The statement of this fact, the H-theorem, has been shown to be false[76]. The emergence of the
Navier-Stokes equations from the BTE, upon introducing the integrals of motion (energy, mass, and momentum),
remains a puzzle, since the very same equations are obeyed by liquids, for which the BTE does not apply at all! This
being relegated to a “philosophical issue” is besides the point; the governing equations of macroscopic fluid mechanics
require stricter and more logical foundations if it is to make progress in areas such as super-critical flow and far-
from-equilibrium phenomena. For an example of a result from kinetic theory reinterpreted using sound probability
theory, see Jaynes[77].
44
A similar equation without the electric field term is considered foundational to the the description of transport
processes in multi-component gaseous mixtures, the solution of which is carried out by invoking perturbation theory;
the velocity distribution is perturbed about the thermal equilibrium solution (Maxwell distribution) as a series of
Sonine polynomials. Theapproachprovidesarigorousderivationoftransportcoefficientsforthedissipativetransport
of momentum, energy and species[18]. In the current case, where the electron swarm is far from thermal equilibrium,
we need to adopt a different strategy.
35
It is common practice to rewrite the BTE in spherical coordinates:
@f
@t
+vcos
@f
@x
e
m
E
cos
@f
@v
+
sin
2
v
@f
@cos
=C[f]
where v is the magnitude of the velocity, is the angle between the velocity and the field
direction, and x is the position along this direction. The velocity distribution can now be expanded
itintermsofLegendrepolynomials. Theaccuracyofthesolutionincreaseswiththenumberofterms
retained in this expansion. For a majority of discharge conditions, it is sufficient to expand to two
terms, wherein the first can be considered an isotropic contribution and the second is proportional
to the degree of anisotropy (caused by the presence of the electric field):
f(v;cos;x;t) =f
0
(v;x;t) +f
1
(v;x;t)cos
Substituting this into the BTE expressed in spherical coordinates, and rewriting velocity in
terms of energy ,
@f
0
@t
+
3
1=2
@f
1
@x
3
1=2
@
@
(Ef
1
) =C
0
@f
1
@t
+
1=2
@f
0
@x
E
1=2
@f
0
@
=N
m
1=2
f
1
where
= (2e=m)
1=2
is a constant, and = (v=
)
2
is the electron energy measured in electron
volts. The right hand side of the first equation denotes the change in f
0
due to collisions and the
m
refers to the total momentum scattering cross-section (which is the sum of the cross-section of
each species weighted by their respective mole fractions). A further simplification may be applied
to the distribution function: it can be separated into a product of terms, one which is dependent
on the electron energy solely and the other which varies only with space and time:
f
0;1
(;x;t) =
1
2
3
F
0;1
()n(x;t)
The numerical solution of F
0;1
is outlined in [78]. The input parameters to this code are the
reduced electric field, the electron density, and the electron collision cross-sections for various pro-
cesses. Upon solving the BTE, we can estimate the rate constant of a reaction i:
k
i
=
1
0
i
F
0
d
Once we are in possession of the rate constants, we can numerically integrate the species and
36
energy conservation equations, which, for the purposes of understanding only the time behavior of
such chemical systems (no spatial variation in any of the variables) are simply written as (for the
case of a constant pressure adiabatic batch reactor
45
):
dY
k
dt
=
_ !W
k
c
p
dT
dt
=
K
X
k=1
h
k
_ !W
k
where Y
k
is the mass fraction, _ !is the reaction rate, W
k
is the molecular weight, is the density,
c
P
is the specific heat at constant pressure, T is the temperature, h
k
is the enthalpy of formation.
It should be kept in mind that the electrons don’t contribute to any of the properties in the above
two equations since they are at a very different temperature throughout the simulation. Their
temperature is given by the velocity distribution that is a solution to the Boltzmann equation for
thatparticulartime-step. Theabovesetofequationsareverystiff, i.e. thereisgreatdisparityinthe
characteristic timescales of relaxation, making the system of ODEs numerically challenging to solve.
Thankfully, numerical integrators have matured in robustness; a careful consideration of tolerances
using the integrator CVODE[51] can yield reliable solutions. The thermodynamics and kinetic data
for combustion relevant reactions were computed using Cantera
46
while the rate computation of
ion and excited state chemistry were outputs from ZDPLASKIN
47
. These two codes, along with
BOLSIG+[78], were combined into a single code to handle a variety of discharge situations for a
range of fuel-air mixtures. The results of these computations will be discussed in the following
sections.
45
This can be easily generalized to other batch reactor systems (constant volume, plugged flow, continuously
stirred, etc.). See Kee[20] for the governing equations.
46
Can be found at http://www.cantera.org/docs/sphinx/html/index.html
47
Can be found at http://www.zdplaskin.laplace.univ-tlse.fr/
37
5 Chemical Pathways in NOx remediation by non-thermal
plasma
The NOx reduction/remediation is among the highest priorities in today’s environmentally con-
scious energy economy. Based on its properties, the deployment of non-thermal transient plasma
is also a viable method to eliminate NOx from exhaust gases. However, the elementary kinetic
pathways leading to NOx reduction have yet to be identified. Upon discharge, the resulting elec-
tron impact reactions give rise to pathways that depend on the oxidizing environment. According
to Penetrante et al.[79, 80] the plasma discharge leads to the formation of N radicals (neutral and
excited), following by
N +NO!N
2
+O (R1)
At high enough temperatures, this reduction is offset by the oxidation of NO to NO2 via
O +NO +M!NO
2
+M (R2)
O +O
2
+M!O
3
+M (R3)
O
3
+NO!NO
2
+O
2
(R4)
Thedominanceofthispathwaydependsontheamountofoxygeninthevicinityofthedischarge.
While lean-burn exhausts are conducive to this oxidation of NO to NO
2
, the OH radicals formed
via
O +H
2
O! 2OH (R5)
proceed to convert NO and NO2 to acids. In power plants, these acids can be treated with
ammonia to form salts [81]. On the other hand, in marine diesel engines the exhaust could be dis-
solved into the surrounding water after appropriately reducing pH levels. Recently, Manivannan et
al.[82] have proposed a technique involving electron beam microwave discharges for treating marine
diesel exhaust. Although this method is known to be most energy efficient, its design complexity
can be restrictive. Puchkarev and Gundersen[70] have shown that a simpler high repetition rate
corona discharge can achieve similar results at a comparable cost to the electron beam processing.
Studying parametrically the underlying plasma and chemical kinetics in a batch reactor, either of
these methods could be optimized. Solving for the full set of transport equations to account for
the spatio-temporal nature of streamer discharge can indeed reveal the coupling of all scales, but
38
is prohibitively expensive in most cases. On the other hand, a zero-dimensional analysis can assess
most pertinent details for a fraction of the cost.
The non-equilibrium nature of the plasma requires the solution of a Boltzmann transport equa-
tion for the electrons in addition to the chemical species conservation equations. In the present
study, the solution was carried out using the ZDPLASKIN code, which makes calls to BOLSIG+
for the evaluation of the reaction rates of electron impact reactions and CVODE for the time inte-
gration. The PUMPKIN
48
code was used for the reaction path analysis. The approach involved the
introduction of reduced electric field and electron densities as functions of time. A simple Gaussian
profile was used for the electron densities (n
e
) and a parabolic profile for the reduced electric field
(E), corresponding to representative experimental measurements of current and voltage. Three
different discharges were analyzed: 1. A 40 kHz repetition rate with a maximum electron density of
510
11
cc
1
2. A 4 kHz repetition rate with a maximum electron density of 510
10
cc
1
3. A 0.4 kHz
repetition rate with a maximum electron density of 5 10
9
cc
1
: The maximum reduced electric
field in all three cases was 70 Td and the discharge duration was 1 millisecond. The kinetic model
comprises reactions assembled from various sources. N
2
O
2
reactions were taken from Capitelli
et al. [83], Gordiets et al. [84], and Kossyi et al. [85], H
2
O reactions from Komuro et al. [86], and
CO
2
reactions from Dorai and Kushner[87]. Electron Scattering cross-sections for ionization, at-
tachment and excitation reactions were obtained from a central database [88]. Various proportions
of secondary diluents (O
2
, CO
2
and H
2
O, hereby abbreviated SD) in N2 (72% by volume in all
cases) were submitted into the zero-dimensional constant volume reactor, resulting in 6 cases that
were investigated. It is well known that unburned hydrocarbons (e.g., alkenes and aromatics) are
present in the exhaust in ppm levels, but due to the lack of available chemical kinetic models there
were not considered herein. The exhaust gas was fixed at a temperature of 500 K and a pressure
of 1 atm. The initial NO concentration was 2000 ppmv in all cases.
Case 1. N2:O2 = 72:28 molar ratio: The NO reduction efficiency is very poor for all
repetition rates (Figs. 10a-c). This indicates the dominance of R2, R3, and R4 leading to an NO
2
build-up, as can be seen in Figs. 10c-e. The oscillations seen in the NO and NO
2
profiles are
indicative of their direct coupling to the electron density.
Case 2. N2:O2:CO2 = 72:14:14 molar ratio: The reduction of NO is opposed by its
formation via
CO
2
+N!CO +NO (R6)
The slight net decrease in NO formation compared to Case 1 is only due to the presence of
lesser O
2
, as is apparent from the NO
2
build up.
Case 3. N2:CO2 = 72:28 molar ratio: Although the presence ofCO
2
retards NO reduction,
48
Can be found at http://www.pumpkin-tool.org/
39
R6 is still slower than R1. The absence ofO
2
however, implies that the oxidation pathways R2, R3,
and R4 are not effective due to a lack ofO radicals, as is also the case for theNONO
2
formation
cycle
NO
2
+O!NO +O
2
(R7)
Case 4. N2:O2:H2O = 72:14:14 molar ratio: When CO
2
is replaced by H
2
O, a dramatic
decrease in NO and NO
2
is observed. This is primarily due to R5, which is followed by the
formation of acids via three-body reactions:
OH +NO +NO
2
!HNO
2
+N
2
(R8)
OH +NO
2
+N
2
!HNO
3
+N
2
(R9)
R8 and R9 are the reactions that make repeated plasma discharge for marine diesel exhausts
feasible. The resulting acids can be suitably diluted and discharged into the surrounding water,
thus avoiding any bulky ammonia treatment. A noteworthy aspect about the 40 kHz discharge
rate case is that the efficiency of NO removal decreases after about half a millisecond (Fig. 10a).
Reaction path analysis reveals that the rapid electron density variation intensely couples with the
impact reaction
N
2
+e!N +N(2D) +e (R10)
These metastable nitrogen radicals, while rapidly accumulating, gradually stabilize into NO via
N(2D) +O
2
!NO +O (R11)
The rate of R10 is low at the lower repetition rates because not enough N(2D) has built up
over the course of a discharge cycle. As a result, NO
2
also gradually increases (Fig. 10d).
Case 5. N2:CO2:H2O = 72:14:14 molar ratio: OH radicals resulting from H
2
O improve
the NO removal efficiency compared to Case 4. The high frequency discharge case (Figs. 10a and
10d) show an enhancement compared to Case 4 due to the metastable nitrogen radicalsN(2D).
Case 6. N2:H2O = 72:28 molar ratio: Maximum NOx removal efficiency is reached via
R5, R8, and R9. NO concentrations decrease by two orders of magnitude compared to all other
cases.
It was shown that the optimization of non-thermal plasma treatment of exhaust gases relies
heavily on the kinetic pathways brought about by the mixture surrounding the discharge. These
40
pathways can vary significantly due to the coupling of the timescales of reactions and the repetition
rate of the plasma. The presence of water can greatly enhance NOx remediation, but only for
sufficiently low repetition rates. Further investigation into the role of hydrocarbons in the exhaust
is ongoing and will be the subject of a subsequent study.
41
Figure 10: NO and NO
2
concentrations as a function of time for various discharge frequencies and exhaust gas compositions.
42
6 TransientplasmaeffectsontheautoignitionofDME/O2/Ar
and C3H8/O2/Ar mixtures
The low to intermediate temperature chemistry, i.e. in the 500 K < T < 950 K range, can play
an important role in various combustion applications. A currently pertinent example is compres-
sion ignition engines, in which certain fuels exhibit notably better ignition propensity than others;
dimethyl ether (CH
3
OCH
3
or DME) for instance has remarkable autoignition characteristics, and
is considered a promising alternative fuel[89]. However, fuels with enhanced low temperature chem-
istry also exhibit the so-called Negative Temperature Coefficient (NTC) behavior, i.e. there exists
a temperature range within which an increase in temperature results in a decrease in reactivity [90].
The NTC results almost always in a two-stage ignition with consequences for engine design, which
are difficult to understate. For example, in homogeneous charge compression ignition (HCCI) en-
gines, thermal stratification is purposefully set up in the cylinder in order to control the heat release
rate[91] ; the central core ignites first followed by the peripheral gas. However, if the fuel-air mixture
is prone to ignite slower in the core region and faster in the periphery, which is a plausible situation
in the NTC regime, the performance of the HCCI engine could degrade appreciably. Last but not
least, in spark-ignition engines the events leading to knock are controlled largely by low tempera-
ture chemistry[92]. It is in this context that the control of such processes through non-equilibrium
transient plasma (TP) assistance is of growing importance [15]. The TP-assistance comes in the
form of radicals without relying on thermal energy, essentially decreasing the net global activation
energy. Two recent investigations have shed light on the coupling of plasma discharges with low
temperature chemistry. The study by Filimonova[93] was carried out at pressures ranging from
0.1 to 3 atm and the results showed that the TP-induced radicals greatly reduce and in certain
cases even eliminate the NTC behavior of mixtures of propane (C
3
H
8
) with air. Sun et al.[94]
performed experiments at 72 Torr and measured high formaldehyde (CH
2
O) concentrations prior
to the ignition of DME=O
2
=He non-premixed flames, indicating a distinct period of TP-assisted
low temperature chemistry. Capturing the effects of the NTC behavior on the ignition delay time,
ign
, through a combustion kinetics model, regardless of its response to TP discharges, is an area
of unceasing research. The kinetic underpinnings of NTC behavior have undergone many revisions,
the most recent of which, authored by Merchant et al.[95], is notable for its expository take on
the subject. In order to facilitate the explanation of TP effects on ignition, a brief review of the
interpretation of low temperature chemistry is presented next.
Fig. 11 summarizes the key steps that comprise typical low temperature combustion of hydro-
carbons. This regime of combustion is dictated by the rate of what has been considered to lie at
the heart of low temperature chemistry, the reaction
R +O
2
$RO
2
(RXN1)
43
Figure 11: Key pathways in low temperature hydrocarbon ignition (taken from Merchant et al. [7])
Blue arrows lead to chain branching and red arrows lead to chain termination. Chain propagation
is denoted by black arrows.
44
Figure 12: Equilibrium constants (at 10 atm) for RXN1 corresponding to CH4, C2H6, C3H8, and
DME. Equilibrium constants (at 10 atm) for RXN2 are also shown for C3H8 and DME. Note the
crossover temperatures.
where R is an alkyl radical. At sufficiently low temperatures, the forward rate of RXN1 is
favored while that of the isomerization
RO
2
$QOOH (RXN2)
is not; Q results from one H abstraction from R. Thus, the peroxy radical RO2 accumulates and
drivestheformationofthealkylhydroperoxideROOH,whichsubsequentlyleadstochainbranching.
A further increase in temperature can bring the system to a point where the backward rate
of RXN1 begins to be favored, when the equilibrium constant is less than unity. Note that this
happens only for certain fuels, namely for those that exhibit NTC behavior; in Fig. 12 it can be seen
45
that bothC
3
H
8
and DME possess this feature. The reversal of RXN1 has far reaching consequences
on the mixture’s ability to sustain chain branching. The ROOH formation gradually ceases with
increasing temperature, forcing the oxidation to proceed via RXN2, which is relatively much slower.
The mixture can autoignite only if it can sustain processes leading to OH formation through chain
branching, and this can happen only if a sufficient amount of QOOH accumulates as shown in Fig.
11. The OH radicals formed as a result are partly used to oxidize the fuel. Thus, the sufficient
condition for autoignition is the positive feedback of OH radicals. Two-stage ignition occurs mainly
because this feedback is stalled due to the formation of the ketohydroperoxide, OQ’OOH, as it will
be discussed in a later section. At sufficiently high temperatures, the chain branching relies almost
exclusively on the accumulation of hydrogen peroxide (HO
2
) radicals that are initially dormant.
The effects of TP on NTC chemistry can be understood by observing how the aforementioned
pathways are perturbed. In the recent studies[93, 94], partial answers have been given to the
question “Can plasma be used to fine-tune combustion systems at low temperatures?” The present
work aims to contribute a partial answer as well by delineating the effects of TP in a parametric
way, analyzed with an emphasis on the chemistry of neutral species. This is justified since the time
scales of TP chemistry are notably shorter compared to the time scales of autoignition, so that
one can reduce the effect of a TP discharge to simply the formation of a certain concentration of
radicals.
6.1 Modeling Approach
Zero dimensional isobaric batch-reactor simulations were carried out using CVODE for the time
integration of the species and energy conservation equations. Cantera was used for combustion
related rate calculations and ZDPLASKIN for TP related rate calculations. The combustion chem-
istry models used were those of Burke et al. [96] (for DME, CH
4
, and C
2
H
6
) and Merchant et
al. [95] (for C
3
H
8
). The TP kinetic model is identical to that used by Kosarev et al. [13] for
CH4, C2H6, and C3H8. In the present study, all alkyl radicals formed from electron excitation of
propane are exclusively n-propyl radicals. As it turns out, this is a generous estimate, as Merchant
et al.[95] have shown, the isopropyl radicals contribute to a pathway that diverts OH from chain
branching in the NTC regime. This will be made clear in the forthcoming discussion. For DME,
the electron scattering cross-sections of C
2
H
6
were used following Sun et al. [94]. It was assumed
that reactions involving electrons and ions do not contribute to any temperature increase in the
gas. The time integration of the Boltzmann equation for electron transport (called by ZDPLASKIN
via BOLSIG+) is carried out for every time-step in which the TP discharge is present, i.e. all TP
kinetic rates are evaluated simultaneously alongside combustion kinetics. In should be noted, that
the goal of the current study is not to assess the effects of different types of TP (dielectric barrier,
corona, streamer, etc.) and their various idiosyncratic effects on combustion chemistry. Since it is
46
Figure 13: Variation of electric field and electron density over one discharge.
desired to keep all the TP parameters the same across a variety of mixtures, a simple 50 ns pulse
of TP discharge with a reduced electric field (i.e. electric field per molecule) was adopted, which
varied in time as
E
E
max
=jsin(2ft)j
where t is time, f=20MHz, and E
max
= 100Td for all cases. The electron concetration (n
e
,
measured in particles per cc) varied in time as:
n
n
max
=exp
10sin
2
2f(t
2
)
As can be seen in Fig. 13, the electric field varies over a much longer time compared to
the electron density, in keeping with the basic assumptions of the two term expanded Boltzmann
Transport Equation.
n
max
was varied from 10
11
cc
1
to 10
13
cc
1
. These limits are referred henceforth to as “weak
TP” and “strong TP” respectively. This was done in order to quantify the combustion response (in
this case the decrease in
ign
) for a certain increase in n
max
. To that end, the term “plasma assist
effectiveness” (PAE) is introduced and is defined as the ratio of the change in
ign
when the weak
TP is applied to the change in
ign
when the strong TP is applied (see Fig. 14). In this way, we
can compare the improvements in
ign
for a fixed change in TP intensity. PAE can be interpreted
as the effectiveness of the weakest TP in comparison to the strongest one. Therefore, a high value
of PAE indicates the relative ease with which TP brings about ignition while a low value of PAE
47
Figure 14: The calculation of the plasma assist effectiveness at a given temperature is the ratio of
the vertical distances ‘a’ to ‘b’.
indicates a difficult enhancement.
All computations were carried out at a pressure of 10 atm that is relevant to those encountered
in engines, for DME and C
3
H
8
since both exhibit NTC behaviors with some distinct differences,
and Ar as the inert to circumvent complications stemming from rotational and vibrational states
of N
2
.
6.2 Results and Discussion
The results of calculations for stoichiometric mixtures are shown in Fig. 15. Each mixture’s
ign
is
shown as a region bounded by an upper curve and a lower curve. The upper curve corresponds to
the case of a purely thermal autoignition while the lower one corresponds to the case of autoignition
in the presence of the strong TP. The filled region in between represents the infinitude of curves
produced by varying TP intensities between the strong and weak cases. The
ign
values are seen
to decrease for most of the temperature range in the order C
3
H
8
< C
2
H
6
< CH
4
, as expected;
the dissociation energies of the C-H bond decreases in the same order for those molecules. DME
possesses the least
ign
due to the presence of the bond weakening O atom.
48
Figure 15: The ignition delays and plasma assists of the mixtures CH
4
:O
2
:Ar(1:2:7.52), C
2
H
6
:
O
2
:Ar (1:3.5:13.76),CH
3
OCH
3
:O
2
:Ar (1:3:11.28), andC
3
H
8
:O
2
:Ar (1:5:18.8). Filled regions
represent the infinitude of ignition delay time curves parametrized by increasing plasma intensities.
49
TP discharge is seen to decrease each mixture’s
ign
by an amount dictated by their chain
branching propensity. CH
4
is known to be an inhibitor due to the removal of H radicals via
CH
4
+H$CH
3
+H
2
(RXN3)
which eventually suppresses chain branching to O and OH.
Since TP results in alkyl radicals via electron collisions without forsaking the H radical via
RXN3, O and OH radical production proceeds via
CH
3
+O
2
$CH
3
O +O (RXN4)
CH
3
+O
2
$CH
2
O +OH (RXN5)
Note that the PAE is nonetheless poor for CH4 over a wide range of temperatures due to an
increase in the backward rate of
CH
4
+HO
2
$CH
3
+H
2
O
2
(RXN6)
Itiswellknownthathydrogenperoxide(H2O2)formationisthedominantprecursortoexplosive
chain branching at high temperatures for all alkanes via the high activation energy reaction
H
2
O
2
+M$OH +OH +M (RXN7)
This is why PAE for all mixtures decreases at high temperatures; TP-induced R radicals tend
to drive backwards RXN6, and its analogues for other alkanes. The fact that PAE increases before
decreasing close to 1000 K is evidence of the change in chain branching regimes. C
2
H
6
is a peculiar
molecule in that when C
2
H
6
=O
2
mixtures are subjected to a TP discharge, C
2
H
5
, C
2
H
3
, and O
radicals are formed, alongside C
2
H
4
. C
2
H
4
is susceptible to attack by O radicals, and can thus
enhance chain branching [2]:
C
2
H
4
+ (H;OH)$C
2
H
3
+ (H
2
;H
2
O) (RXN8)
The significant enhancements in autoignition due to TP discharge seems to be reserved for fuels
exhibiting multiple complex pathways at low temperatures. As mentioned in the Introduction,
at low enough temperatures, RXN1 proceeds forward and the resulting build-up of RO
2
leads to
the formation of ROOH, which is responsible for chain branching as discussed earlier. Therefore,
when C
3
H
8
and DME reacting mixtures are subjected to a TP discharge, the resulting R radicals
50
advance the system towards chain branching quicker. This can be seen in Figs. 16 and 17; when
the initial temperature is 550 K, the ROOH production is largely responsible for ignition. Hence
the significantly high PAE at very low temperatures. The reason the TP discharge is more effective
for DME lies in the fact that DME cannot undergo the step
RO
2
$alkene +HO
2
(RXN9)
while C3H8 can, and this results in an additional pathway that suppresses the branching propen-
sity of C3H8 reacting mixtures.
Around the temperature at which the equilibrium constant for RXN1 becomes less than unity,
the production of R radicals by TP results in the reduction of PAE for both C
3
H
8
and DME. For
slightly higher temperatures, OH chain branching [7] almost exclusively relies on the accumulation
of QOOH radicals. For this to happen, RO
2
must be available in sufficient quantities, but if the
temperatureissufficientlyhigh, theequilibriumofRXN1reverses, makingtheadditionofRradicals
by TP relatively inefficient. The fuel is essentially consumed by two channels via OH and HO
2
, to
form R radicals. In the absence of TP, it can be seen in Fig. 16 that for an initial temperature of
650 K the temperature profile exhibits two stages. On the other hand, the radical profiles exhibit
three stages, marked ast
1
,t
2
, andt
3
. Duringt
1
, fuel is consumed largely by OH compared toHO
2
.
This HO
2
builds up and at the beginning of t
2
, enough HO
2
exists to initiate
HO
2
+HO
2
$H
2
O
2
+O
2
(RXN10)
which manifests as a slower growth rate in its concentration profile. In this regard, at a compa-
rable temperature (775K), DME differs markedly in that it is mainly consumed by HO
2
and thus
the initiation of RXN10 is postponed, effectively eliminating the distinct stages t
1
and t
3
. But at
higher temperatures, DME is largely consumed by OH compared to HO
2
, and so HO
2
builds up
(similar to the C
3
H
8
=O
2
=Ar mixture at 650K) causing RXN10 to initiate abruptly, much before
autoignition (see Fig. 17 where T
0
= 950K).
For theC
3
H
8
case, in the presence of TP, it is seen that staget
1
is completely eliminated because
the introduction of R radicals frees up HO
2
to form H
2
O
2
. However, this does not shorten the
second ignition delay stage (t
2
+t
3
) as much due to a mechanism elucidated by Merchant et al.[95];
by shortening t1 with TP, heat has been released quicker and consequently, the reactions RXN1
and
QOOH +O
2
$O
2
QOOH (RXN11)
are promoted backwards, resulting thus in an OQ’OOH buildup. This buildup is directly re-
lated OH radicals being consumed via pathways that divert it from a positive feedback loop that
51
Figure 16: A comparison of key radical concentrations for stoichiometric C3H8/O2/Ar mixtures without and with a strong
plasma discharge at 10 atm. Note that the time has been scaled by the corresponding ignition delay times.
52
Figure 17: A comparison of key radical concentrations for stoichiometric DME=O
2
=Ar mixtures without and with a strong
plasma discharge at 10 atm. Note that the time has been scaled by the corresponding ignition delays.
53
contributes to fuel consumption and subsequent chain branching. It is important to note that these
pathways are a result of the formation of secondary products like CH
2
O and H
2
O
2
:
CH
2
O +OH$HCO +H
2
O (RXN12)
H
2
O
2
+OH$HO
2
+H
2
O (RXN13)
Eventually, the OQ’OOH accumulates to a point after which it rapidly releases OH, which upon
consumption and accompanying heat release, causes the first stage ignition. It can be seen from
Figs. 16 and 17 that for this regime (650 K and 775 K for C
3
H
8
and DME respectively), the TP
discharge does not influence the processes following the first ignition stage (t
1
+t
2
) and the reason
for this is now clear: the very acceleration in fuel consumption brought about by TP also causes
a rapid accumulation in OQ’OOH because secondary products like CH
2
O and HO
2
have been
formed that much faster. Thus, the barrier faced by TP-induced radicals is due to the OQ’OOH
accumulation.
At higher temperatures, OQ’OOH is no longer a bottleneck, and the PAE recovers (see Fig. 15)
albeit briefly, before succumbing to the reversal of reactions such as RXN6.
Athigherpressures, theNTCregionmovestowardshighertemperaturesbecausetheequilibrium
constant of RXN1 is pressure dependent. As a consequence, following from the foregoing analysis,
the range over which the PAE values are high is expected to increase (see Fig. 18).
Itshouldbementionedalsoinpassing,thatincreasingtheequivalenceratioforbothC3H8/O2/Ar
and DME/O2/Ar mixtures decreases
ign
, mainly because of an increase in the number of active
radicals at lower temperatures. PAE increases accordingly (see Figs. 19 and 20).
6.3 Concluding Remarks
The isobaric homogeneous batch reactor simulations indicate thatCH
4
=O
2
=Ar and C
2
H
6
=O
2
=Ar
mixtures, when subjected to a transient plasma (TP) discharge of varying intensities, appears to
possess a low Plasma Assist Effectiveness (PAE), in contradistinction to mixtures of DME=O
2
=Ar
andC
3
H
8
=O
2
=Ar, whose relatively higher PAEs for the same temperature range is a consequence of
favorable low temperature kinetics, specifically the ROOH pathway to chain branching. However,
the PAEs for these mixtures decrease as temperature increases since the Negative Temperature
Coefficient (NTC) regime contains chain branching pathways that lead to a two-staged ignition,
only the first stage of which is susceptible to a decrease in duration when subjected to TP. This is
due to the gradual buildup of OQ’OOH, a significant source of OH radicals and consequently chain
branching. Upon the release of these OH radicals (an event that marks the beginning of the second
ignition stage) the effects of TP induced radicals are effectively “lost”, drastically reducing the PAE.
54
Figure 18: The variation of ignition delay time for stoichiometric DME/O2/Ar mixtures and
C3H8/O2/Ar mixtures at two different pressures.
55
Figure 19: Variation of ignition delays with equivalence ratios for C3H8/O2/Ar mixtures at 10 atm.
Filled regions represent the infinitude of ignition delay time curves parameterized by increasing
plasma intensities.
56
Figure 20: Variation of ignition delays with equivalence ratios for DME/O2/Ar mixtures at 10 atm.
Filled regions represent the infinitude of ignition delay time curves parameterized by increasing
plasma intensities.
57
At temperatures higher than the NTC regime, the PAEs undergo a further decrease because the
TP induced radicals only slightly perturb the high activation energy production of OH radicals by
H
2
O
2
decomposition. The analysis reveals that the PAEs are high for a larger range of temperature
when the pressure is increased because the NTC regime is shifted towards higher temperatures. The
increase of equivalence ratio also increases the PAE, albeit marginally. Controlling the timing of
the plasma discharge could perhaps restore the PAE of a reacting mixture at NTC temperatures.
For instance, the OQ’OOH buildup that is characteristic of two staged ignition could be lessened
by bombarding the radical with electrons. This would be very effective since the dissociation
energy of the O-O bond in this radical is particularly small. If plasma were discharged during
such intermediate steps, the reduced chain branching due to thermally controlled bottlenecks may
lessen or even disappear. Reliable electron scattering cross-sections and the accompanying plasma
kinetics are currently unavailable for species that are formed during the intermediate steps, but
the foregoing results indicate that acquiring these will hold the key to controlling combustion in
compression ignition engines.
7 Summary and Looking Ahead
We have managed to use and modify a multidimensional reacting flow solver to examine an im-
portant problem in Combustion Science: the effects of a probe on a burner stabilized flame in the
context of Molecular Beam Mass Spectrometry. We quantified the extent to which different aspects
of the apparatus perturb the flame and hence the observed chemistry, indicating that conventional
one-dimensional codes are inadequate in modeling the phenomenon. This is one example of multi-
dimensional effects on flames which can only be explored using an infrastructure such as the one
used here.
With the confidence in its ability to capture essential physical phenomena, we have ventured
into trying to capture the effects of non-thermal plasma on combustion chemistry. Using theoretical
estimates on what the discharge parameters (reduced electric field and electron density) ought to
be for various types of discharges, we were able to model the non-thermal plasma kinetics in batch
reactors. We have shown how NOx can be remedied in exhaust gases by using particular combina-
tions of discharge parameters and the non-trivial effects of non-thermal plasma (i.e. its resultant
radicals) on low temperature chemistry, a matter of growing importance in the context of Homo-
geneous Charge Compression Ignition (HCCI) engines. Using the Operator Splitting techniques
discussed in Chapter 2, we are in a position to combine these solvers into a larger infrastructure.
The computational subroutines used in the preceding work have been written for modularity
with C, C++, and fortran codes. This allows any existing combustion code to interface with the
plasmasubroutinesinastraightforwardmanner. Workiscurrentlybeingdoneinthisarea; speeding
up the integration of the stiff set of equations is where most of the focus currently is. Although
58
Figure 21: The OH production rate for two different discharge energies in an opposed jet configu-
ration. The discharge is applied at the center of the domain, where the mixing between fuel and
oxidizer streams is maximum. It is seen that the 1J/cc case leads to ignition while the 0.65J/cc
does not.
a complete code that couples combustion and plasma processes spatially is yet to be written, as a
first cut, we can model the resulting radicals of a discharge as a source term in the energy & species
equations, neglecting the complexities of electron transport at timescales much smaller than the
kinetics that constitute combustion.
For instance, it is possible to add a source term that stimulates ignition in an opposed-jet
configuration, where the strain rate imposes a residence time for the active radicals formed during
the ignition of a diffusion flame
49
. In Figure 21, this is demonstrated for a H
2
=O
2
=N
2
diffusion
flame, where the fuel stream comprises 20%H
2
diluted in N
2
at 300 K, the oxidizer stream is air
at 300 K, and the imposed global strain rate is 300s
1
at a pressure of 1atm. It can be seen that
the minimum ignition energy is strain-rate dependent. The number of H and O radicals formed
are derived using simple estimates which use weighting coefficient that assume an average collision
cross-section for dissociation due to electron impact, and the bond dissociation energy. A detailed
plasma chemical model might not be required here, but the results of zero-dimensional simulations
will yield, approximately, the correct radical pools that can be used to accurately model the source
terms.
Another example is examining the aiding of spray combustion (typically relevant to diesel en-
gines) by a repeated plasma discharge in the vicinity of the fuel injector. Preliminary simulations
were run using the OpenFOAM
50
suite. Figure 22 shows a typical scenario realized in diesel engines.
Such discharges can be better quantified using tools similar to those used in the zero-dimensional
49
The method to simulate the transient opposed jet simulation is outlined in Egolfopoulos and Campbell[97].
50
Can be found at http://openfoam.org/
59
Figure 22: A mass of liquid heptane is sprayed into an environment (height=10cm and width=1cm)
that is at a pressure of 50atm and a temperature of 800K. The left-most column indicates the
position of the plasma discharge, which occupies a 2cc spherical region modeled as a reaction rate
enhancement to the Arrhenius collision frequency for the one step heptane combustion model.
60
calculations. The framework still lacks robustness due to different standards in coding and organiz-
ing used in different combustion codes. The intent of this work was to contribute in the direction
of integrability using free software. Many research groups are currently working towards this, and
it may soon be possible to design reliable and robust kinetically controllable combustion systems.
61
References
[1] Paul M Dirac. The Principles of Quantum Mechanics. Oxford university press, 1978.
[2] Linus Pauling and Bright E Wilson. Introduction to Quantum Mechanics with Appications to
Chemistry. Courier Corporation, 1935.
[3] Richard Phillips Feynman. Forces in molecules. Physical Review, 56(4):340, 1939.
[4] Max Born and Robert Oppenheimer. On the quantum theory of the molecules. 1924.
[5] Lev D Landau and Evgeny M Lifshitz. Chapter 01 - The Basic Concepts of Quantum Mechan-
ics. In A Course of Theoretical Physics Volume 3. Perganon Press, 1965.
[6] Mehran Kardar. Statistical Physics of Particles. Cambridge University Press, 2007.
[7] Ralph Howard Fowler and Edward Armand Guggenheim. Statistical thermodynamics. 1941.
[8] Keith J Laidler. Reaction kinetics: Homogeneous gas reactions, volume 1. Elsevier, 2013.
[9] Ludwig Boltzmann. Lectures on gas theory. Univ of California Press, 1964.
[10] Lev D Landau and Evgeny M Lifshitz. Chapter 01 - Ideal Fluids. In A Course of Theoretical
Physics Volume 6. Addison-Wesley, 1958.
[11] Ya B Zeldovich and Yu P Raizer. Physics of Shock Waves and High-Temperature Hydrody-
namic Phenomena: Volume 1, 1965.
[12] Forman A Williams. Combustion Theory. Addison-Wesley, 1985.
[13] Felix J Weinberg. The first half-million years of combustion research and todays burning
problems. Progress in Energy and Combustion Science, I:17–31, 1975.
[14] Irvin Glassman. Supersonic flight and cooking over wood-burning stoves: Challenges to the
combustion community. Proceedings of the Combustion Institute, 28:1–10, 2000.
[15] Yu Starikovskii, NB Anikin, IN Kosarev, EI Mintoussov, SM Starikovskaia, and VP Zhukov.
Plasma-assisted combustion. Pure and Applied Chemistry, 78(6):1265–1298, 2006.
[16] Lev D Landau and Evgeny M Lifshitz. Chapter 02 - Viscous Fluids. In A Course of Theoretical
Physics Volume 6. Addison-Wesley, 1958.
[17] Lewis Fry Richardson. Weather prediction by numerical process. Cambridge University Press,
2007.
62
[18] S Chapman and T G Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge
University Press, 1970.
[19] NilsHansen, Terrilla.Cool, PhillipR.Westmoreland, andKatharinaKohse-Hoinghaus. Recent
contributions of flame-sampling molecular-beam mass spectrometry to a fundamental under-
standing of combustion chemistry. Progress in Energy and Combustion Science, 35(2):168–191,
April 2009.
[20] Robert J Kee, Michael E Coltrin, and Peter Glarborg. Chemically Reacting Flow. John Wiley
& Sons, 2003.
[21] J O Hirschfelder and C F Curtiss. Integration of stiff equations. Proceedings of the National
Academy of Sciences, 38(1950):250–253, 1952.
[22] Charles Hirsch. Numerical Computation of Internal and External Flows. Butterworth-
Heinemann, 1991.
[23] Uri M Ascher and Linda R Petzold. Computer Methods for Ordinary Differential Equations
and Differential-Algebraic Equations, volume 61. Siam, 1998.
[24] Federico Perini, Emanuele Galligani, and RD Reitz. An analytical jacobian approach to sparse
reaction kinetics for computationally efficient combustion modeling with large reaction mech-
anisms. Energy & Fuels, 26:4804–4822, 2012.
[25] J H Ferziger and M Peric. Chapter 7: Solution of the Navier-Stokes Equations. In Computa-
tional Methods for Fluid Dynamics. Springer Science & Business Media, 2002.
[26] Linda R Petzold. Differential algebraic equations are not odes. SIAM Journal on Scientific
and Statistical Computation, 3(3):367–384, 1982.
[27] Francis H Harlow and Eddie J Welch. Numerical calculation of time-dependent viscous incom-
pressible flow of fluid with free surface. The Physics of Fluids, 8(12):2182–2189, 1965.
[28] A J Chorin. Numerical solution of the navier-stokes equations. Mathematics of computation,
22(104):745, October 1968.
[29] Rutherford Aris. Vectors, Tensors, and the Basic Equations of Fluid Mechanics. Courier
Corporation, 1962.
[30] Philip M Gresho and Robert L Sani. Incompressible Flow and the Finite Element Method.
John Wiley and Sons, Inc., New York, NY (United States), 1998.
[31] J Kim and P Moin. Application of a fractional-step method to incompressible navier-stokes
equations. Journal of computational physics, 323:308–323, 1985.
63
[32] J B Bell, Phillip Colella, and H M Glaz. A second-order projection method for the incompress-
ible navier-stokes equations. Journal of Computational Physics, 283:257–283, 1989.
[33] JJIM Van Kan. A second-order accurate pressure-correction scheme for viscous incompressible
flow. SIAM Journal on Scientific and Statistical Computing, 7(3):870–891, 1986.
[34] RIIssa. Solutionoftheimplicitlydiscretisedfluidflowequationsbyoperator-splitting. Journal
of Computational Physics, 62(1):40 – 65, 1986.
[35] HGWeller,GTabor,HJasak,andCFureby. Atensorialapproachtocomputationalcontinuum
mechanics using object-oriented techniques. Computers in Physics, 12(6):620, 1998.
[36] Elaine S Oran and Jay P Boris. Numerical Simulation of Reactive Flow, Second edition.
Cambridge University Press.
[37] Eleuterio F Toro. Splitting schemes for pdes with source terms. In Riemann Solvers and
Numerical Methods for Fluid Dynamics, pages 531–542. Springer, 2009.
[38] Gilbert Strang. On the construction and comparison of difference schemes. SIAM Journal on
Numerical Analysis, 5(3):506–517, 1968.
[39] L Raymond, William H Green, Shev Macnamara, and Gilbert Strang. Balanced splitting and
rebalanced splitting. Society of Industrial and Applied Mathematics, 51(6):3085–3105, 2014.
[40] Alberto Cuoci, Alessio Frassoldati, Tiziano Faravelli, and Eliseo Ranzi. Numerical modeling of
laminar flames with detailed kinetics based on the operator-splitting method. Energy & Fuels,
27(12):7730–7753, December 2013.
[41] E Ranzi, A Frassoldati, R Grana, A Cuoci, T Faravelli, A P Kelley, and C K Law. Hierarchical
and comparative kinetic modeling of laminar flame speeds of hydrocarbon and oxygenated
fuels. Progress in Energy and Combustion Science, 38(4):468 – 501, 2012.
[42] Chunsheng Ji, Enoch Dames, Yang L Wang, Hai Wang, and Fokion N Egolfopoulos. Propaga-
tion and extinction of premixed c5-c12 n-alkane flames. Combustion and Flame, 157(2):277 –
287, 2010.
[43] A T Holley, X Q You, Enoch Dames, Hai Wang, and Fokion N Egolfopoulos. Sensitivity of
propagation and extinction of large hydrocarbon flames to fuel diffusion. Proceedings of the
Combustion Institute, 32(1):1157 – 1163, 2009.
[44] J C Biordi, C P Lazzara, and J F Papp. Molecular beam mass spectrometry applied to deter-
mining the kinetics of reactions in flames. i. empirical characterization of flame perturbation
by molecular beam. Combustion and Flame, 23:73–82, 1974.
64
[45] E W Kaiser, T J Wallington, M D Hurley, J Platz, H J Curran, W J Pitz, and C K Westbrook.
Experimental and modeling study of premixed atmospheric pressure dimethyl ether air flames.
The Journal of Physical Chemistry A, 104(35):8194–8206, 2000.
[46] A N Hayhurst, D B Kittelson, and N R Telford. Mass spectrometric sampling of ions from
atmospheric pressure flames-ii.: Aerodynamic disturbance of a flame by the sampling system.
Combustion and Flame, 28:123–135, 1977.
[47] A C Yi and E L Knuth. Probe-induced concentration distortions in molecular-beam mass-
spectrometer sampling. Combustion and Flame, 63(3):369–379, March 1986.
[48] A T Hartlieb, B Atakan, and K Kohse-Hoinghaus. Effects of a sampling quartz nozzle on the
flame structure of a fuel-rich low-pressure propene flame. Combustion and Flame, 121(4):610–
624, June 2000.
[49] Ulf Struckmeier, Patrick Oß wald, Tina Kasper, Lena Böhling, Melanie Heusing, Markus Köh-
ler, Andreas Brockhinke, and Katharina Kohse-Höinghaus. Sampling probe influences on tem-
perature and species concentrations in molecular beam mass spectroscopic investigations of
flat premixed low-pressure flames. Zeitschrift für Physikalische Chemie, 223(4-5):503–537,
May 2009.
[50] P A Skovorodko, A G Tereshchenko, O P Korobeinichev, D A Knyazkov, and A G Shmakov.
Experimentalandnumericalstudyofprobe-inducedperturbationsoftheflamestructure. Com-
bustion Theory and Modelling, 17(April 2013):37–41, 2012.
[51] Peter N Brown, George D Byrne, and Alan C Hindmarsh. Vode: A variable-coefficient ode
solver. SIAM journal on scientific and statistical computing, 10(5):1038–1051, 1989.
[52] J H Ferziger and M Peric. Chapter 5: Solution of Linear Equation Systems. In Computational
Methods for Fluid Dynamics. Springer Science & Business Media, 2002.
[53] Tianfeng Lu and Chung K Law. A directed relation graph method for mechanism reduction.
Proceedings of the Combustion Institute, 30(1):1333–1341, January 2005.
[54] Yousef Saad. Iterative Methods for Sparse Linear Systems. Siam, 2003.
[55] Hrvoje Jasak. Error Analysis and Estimation for the Finite Volume Method with Applications
to Fluid Flows. PhD thesis, 1996.
[56] Lev D Landau and Evgeny M Lifshitz. Chapter 04 - Boundary Layers. In A Course of
Theoretical Physics Volume 6. Addison-Wesley, 1958.
65
[57] P M Gresho and R L Sani. On pressure boundary conditions for the incompressible navier-
stokes equations. International Journal for Numerical Methods in Fluids, 7(10):1111–1145,
October 1987.
[58] M D Smooke. Solution of burner-stabilized premixed laminar flames by boundary value meth-
ods. Journal of Computational Physics, 105:72–105, 1982.
[59] Robert J Kee, Joseph F Grcar, M D Smooke, J A Miller, and E Meeks. Premix: a fortran
program for modeling steady laminar one-dimensional premixed flames. Technical Report
SAND85-8249, 1985.
[60] Fokion N Egolfopoulos, H Zhang, and Z Zhang. Wall effects on the propagation and extinction
of steady, strained, laminar premixed flames. Combustion and Flame, 109(1-2):237–252, April
1997.
[61] Fran M Rupley, Andrew E Lutz, Robert J Kee, and J F Grcar. Oppdiff: A fortran program
for computing opposed-flow diffusion flames. Technical Report May, 1997.
[62] Aamir D Abid, Joaquin Camacho, David A Sheen, and Hai Wang. Quantitative measurement
of soot particle size distribution in premixed flames: The burner-stabilized stagnation flame
approach. Combustion and Flame, 156(10):1862–1870, October 2009.
[63] E T Jaynes. Information theory and statistical mechanics. Physical review, 1957.
[64] J Willard Gibbs. Elementary principles in statistical mechanics. Courier Corporation, 2014.
[65] Claude Shannon. A Mathematical Theory of Communication. PhD thesis, 1948.
[66] J C Keck. Rate-controlled constrained-equilibrium theory of chemical reactions in complex
systems. Progress in Energy and Combustion Science, 16:125–154, 1990.
[67] Lev D Landau and Evgeny M Lifshitz. Chapter 04 - Collisions Between Particles. In A Course
of Theoretical Physics Volume 1.
[68] Richard P Feynman, Robert B Leighton, and Matthew Sands. The Feynman Lectures on
Physics, Desktop Edition Volume II, volume 1. Basic Books, 2013.
[69] Yu P Raizer. Gas Discharge Physics. Springer Berlin, 1983.
[70] V Puchkarev and M Gundersen. Energy efficient plasma processing of gaseous emission using
a short pulse discharge. Appl. Phys. Lett., 71(23):3364–3366, Decemeber 1997.
[71] J B Liu, Paul D Ronney, and Martin A Gundersen. Premixed flame ignition by transient
plasma discharges. In Symposium (International) on Combustion, 2002.
66
[72] E T Jaynes. Probability theory: the logic of science. Cambridge University Press, 2003.
[73] E T Jaynes. Information Theory and Statistical Mechanics II. Physical Review, 1957.
[74] E T Jaynes. The Minimum Entropy Principle. Annual Review of Physical Chemistry, 1980.
[75] E T Jaynes. Macroscopic prediction. In Complex Systems: Operational Approaches in Neuro-
biology, Physics, and Computers, pages 254–269. Springer, 1985.
[76] E T Jaynes. Gibbs vs Boltzmann entropies. Am. J. Phys, 1965.
[77] E T Jaynes. Clearing up mysteries - the original goal. pages 1–27, 1989.
[78] GJMHagelaarandLCPitchford. Solvingtheboltzmannequationtoobtainelectrontransport
coefficients and rate coefficients for fluid models. Plasma Sources Science and Technology,
722(4):722–733, November 2005.
[79] BMPenetrante,MCHsiao,BTMerritt,GEVogtlin,andPHWallman. Comparisonofelectrical
discharge techniques for nonthermal plasma processing of no in n2. IEEE Transactions of
Plasma Science, 23(4):679–687, August 1995.
[80] B M Penetrante, R M Brusasco, B T Merritt, and G E Vogtlin. Environmental applications
of low-temperature plasmas. Pure Appl. Chem., 71(10):1829–1835, August 1999.
[81] Jiri Svachula, Natale Ferlazzo, Pio Forzatti, Enrico Tronconi, and Fiorenzo Bregani. Selective
reduction of nitrogen oxides (nox) by ammonia over honeycomb selective catalytic reduction
catalysts. Industrial & engineering chemistry research, 32(6):1053–1060, 1993.
[82] N Manivannan, W Balachandran, R Beleca, and M Abbod. Non-thermal plasma technology
for the abatement of nox and sox from the exhaust of marine diesel engine. Journal of Clean
Energy Technologies, 2(3):233–236, 2014.
[83] Mario Capitelli, Carlos M Ferreira, Boris F Gordiets, and Alexey I Osipov. Plasma kinetics in
atmospheric gases, volume 31. Springer Science & Business Media, 2013.
[84] Boris F Gordiets, Carlos M Ferreira, Vasco L Guerra, Jorge MAH Loureiro, Jacimar Nahorny,
Daniel Pagnon, Michel Touzeau, and Marinette Vialle. Kinetic model of a low-pressure n 2-o
2 flowing glow discharge. Plasma Science, IEEE Transactions on, 23(4):750–768, 1995.
[85] I A Kossyi, A Yu Kostinsky, A A Matveyev, and V P Silakov. Kinetic scheme of the non-
equilibrium discharge in nitrogen-oxygen mixtures. Plasma Sources Science and Technology,
1(3):207, 1992.
67
[86] Atsushi Komuro and Ryo Ono. Two-dimensional simulation of fast gas heating in an at-
mospheric pressure streamer discharge and humidity effects. Journal of Physics D: Applied
Physics, 47(15):155202, April 2014.
[87] R Dorai and MJ Kushner. Effect of multiple pulses on the plasma chemistry during the
remediation of nox using dielectric barrier discharges. J. Phys. D: Appl. Phys., 34:574–583,
2001.
[88] S Pancheshnyi, S Biagi, M C Bordage, G J M Hagelaar, W L Morgan, A V Phelps, and L C
Pitchford. The lxcat project: Electron scattering cross sections and swarm parameters for low
temperature plasma modeling. Chemical Physics, 398:148–153, 2012.
[89] Troy A Semelsberger, Rodney L Borup, and Howard L Greene. Dimethyl ether (dme) as an
alternative fuel. Journal of Power Sources, 156(2):497 – 511, 2006.
[90] MJ Piling. Low-Temperature Combustion and Autoignition. 1997.
[91] Yi Yang, John E Dec, Nicolas Dronniou, and Magnus Sjoberg. Tailoring hcci heat-release
rates with partial fuel stratification: Comparison of two-stage and single-stage-ignition fuels.
Proceedings of the Combustion Institute, 33(2):3047 – 3055, 2011.
[92] Kwang Min Chun, John B Heywood, and J C Keck. Prediction of knock occurrence in a spark-
ignition engine. In Symposium (International) on Combustion, volume 22, pages 455–463.
Elsevier, 1989.
[93] E A Filimonova. Discharge effect on the negative temperature coefficient behaviour and mul-
tistage ignition in c3h8-air mixture. Journal of Physics D: Applied Physics, 48(1):015201,
2015.
[94] Wenting Sun, Sang Hee Won, and Yiguang Ju. In situ plasma activated low temperature
chemistry and the s-curve transition in dme/oxygen/helium mixture. Combustion and Flame,
161(8):2054–2063, 2014.
[95] Shamel S Merchant, C Franklin Goldsmith, Aaron G Vandeputte, Michael P Burke, Stephen J
Klippenstein, and William H Green. Understanding low-temperature first-stage ignition delay:
Propane. Combustion and Flame, 2015.
[96] UltanBurke,KieranPSomers,PeterOToole,ChisMZinner,NicolasMarquet,GillesBourque,
Eric L. Petersen, Wayne K. Metcalfe, Zeynep Serinyel, and Henry J Curran. An ignition delay
and kinetic modeling study of methane, dimethyl ether, and their mixtures at high pressures.
Combustion and Flame, 162(2):315 – 330, 2015.
68
[97] Fokion N Egolfopoulos and C S Campbell. Unsteady counterflowing strained diffusion flames:
diffusion-limited frequency response. Journal of Fluid Mechanics, 318:1–29, 1996.
69
Abstract (if available)
Abstract
Numerical simulations are an indispensable tool in understanding the interplay between various mechanisms in complex physical phenomena. A computational infrastructure was developed that allows for the solution of conventional and novel reacting flow problems. We applied an algorithm that allows for the solution of the Boltzmann transport equation for electrons in order to address the effects of non-thermal plasma on combustion processes. The ignition of hydrocarbons that exhibit negative temperature coefficient (NTC) chemistry were found to couple very efficiently with non-thermal plasma for certain temperature ranges
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Modeling investigations of fundamental combustion phenomena in low-dimensional configurations
PDF
A theoretical study of normal alkane combustion
PDF
Mesoscale SOFC-based power generator system: modeling and experiments
PDF
Experimental studies of high pressure combustion using spherically expanding flames
PDF
Flame ignition studies of conventional and alternative jet fuels and surrogate components
PDF
Accuracy and feasibility of combustion studies under engine relevant conditions
PDF
Studies of combustion characteristics of heavy hydrocarbons in simple and complex flows
PDF
Studies of siloxane decomposition in biomethane combustion
PDF
Pressure effects on C₁-C₂ hydrocarbon laminar flames
PDF
Investigations of fuel effects on turbulent premixed jet flames
PDF
CFD design of jet-stirred chambers for turbulent flame and chemical kinetics experiments
PDF
Determination of laminar flame speeds under engine relevant conditions
PDF
Studies of Swiss-roll combustors for incineration and reforming applications
PDF
Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle
PDF
Studies on the flame dynamics and kinetics of alcohols and liquid hydrocarbon fuels
PDF
Investigations of fuel and hydrodynamic effects in highly turbulent premixed jet flames
PDF
An experimental investigation by optical methods of the physics and chemistry of transient plasma ignition
PDF
Studies of methane counterflow flames at low pressures
PDF
Experimental and kinetic modeling studies of flames of H₂, CO, and C₁-C₄ hydrocarbons
PDF
Physics and applications of energetic electrons resulting from nanosecond-scale transient plasma
Asset Metadata
Creator
Gururajan, Vyaas
(author)
Core Title
A theoretical investigation in multi-dimensional and non-thermal plasma effects on combustion processes and pollutant remediation
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
04/19/2016
Defense Date
03/25/2016
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
CFD,combustion,kinetics,OAI-PMH Harvest,plasma
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Egolfopoulos, Fokion (
committee chair
), Ronney, Paul (
committee member
), Wittig, Curt (
committee member
)
Creator Email
gururaja@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c40-230904
Unique identifier
UC11276883
Identifier
etd-GururajanV-4280.pdf (filename),usctheses-c40-230904 (legacy record id)
Legacy Identifier
etd-GururajanV-4280.pdf
Dmrecord
230904
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Gururajan, Vyaas
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
CFD
combustion
kinetics
plasma