REDUCTION OF DETAILED CHEMICAL REACTION
NETWORKS FOR DETONATION SIMULATIONS
Patrick Hung
Mechanical Engineering
California Institute of Technology
Pasadena, CA 91125
Joseph E. Shepherd
Graduate Aeronautical Laboratory
California Institute of Technology
Pasadena, CA 91125
While a detailed mechanism represents the state-of-the-
art of what
is known about a reaction network, its direct implementation in a
fully re
solved CFD simulation is all but impossible (except for the
simplest systems) with the computational power available today.
This paper discusses the concept of Intrinsic Low Dimensional
Manifold (ILDM), a technique that systematically reduces the
complexity of detailed mechanisms. The method, originally devel-
oped for combustion systems, has been successfully extended and
applied to gaseous detonation simulations
2,3,4
. Unfortunately, while
a one-d
imensional ILDM is reasonably easy to compute, manifolds
of higher dimensions are notoriously difficult. Moreover, the selec-
tion of the manifold dimension has been largely arbitrary,
with a
one-dimensional ILDM being the most popular if for no other rea-
son than that it is easiest to compute and store.
In this paper, we will present a technique that enables us to quanti-
tatively determine the dimensionality of the ILDM needed, as well
a
s a robust and embarrassingly parallel algorithm for computing
high-dimensional ILDMs. Finally, these techniques are demon-
strated in the context of a one-d
imensional ZND detonation with
detailed chemistry.
INTRODUCTION
Intrinsic Low Dimensional Manifold (ILDM) is
a technique for systematically identifying low-
dimensional attracting submanifolds of the
original state space of chemical reaction mecha-
nisms
1
. The method, originally developed for
low-speed combustion systems, has been suc-
cessfully
extended
and
applied
to
two-
dimensional gaseous detonation simulations
with the Hydrogen-Oxygen reaction mecha-
nism
2
, and gaseous detonation
3
and gas phase
RDX combustion
4
. While detailed reaction
mechanisms are now mature for many hydro-
carbon fuels and in a developmental stage for
nitramine explosives such as RDX and HMX, a
number of issues remain to be addressed before
they can be used in conjunction with the ILDM
method for detonation simulation.
First, the ILDM algorithm is computationally
expensive to apply, and the computed manifold
presents difficult tabulation, storage, and inter-
polation problems. While a one-dimensional
ILDM can be computed reasonably easily and
has been shown to work well for simple reaction
systems such as the Hydrogen-Oxygen reaction
mechanism, it is not reasonable to expect such
drastic amount of reduction to remain faithful to
even moderately complex hydrocarbon mecha-
nisms. We will present an algorithm that allows
us to determine, quantitatively, the number of
dimensions needed. The reason for using
ILDMs in detonation simulation is simple; we
want to extract as much information from the
detailed mechanism as we can afford, and as
little as we need. The ILDM technique allows us
to follow (if not reach) this goal systematically.
Unfortunately,
algorithms
for
computing
ILDMs, the most popular being continuation
methods, are far from robust. In this paper, a
new algorithm for the computation of ILDMs
that is more efficient, embarrassingly parallel,
and far more robust than continuation methods
is presented.
Finally, these techniques are applied to a one-
dimensional ZND detonation, giving us valuable
insights as well as demonstrating clearly a
“stepping-down” of dimensions as we move
away from the leading shock front. In other
words, the closer to the front we need to capture
the reaction dynamics, the higher the dimension
we need. Finally, remarks concerning the appli-
cation of this technique, as well as ramifications
of some of the underlying assumptions, are dis-
cussed.
INTRINSIC
LOW
DIMENSIONAL
MANIFOLDS
By using an operator split scheme to the reactive
Euler equations
2
, each of the finite volume in
the discretized domain during the reaction
source step is a constant-volume adiabatic com-
bustor. The governing equation is a system of
ODE, which can be written as,
( ;
, )
d
f
e
dt
φ
φ
ρ
=
(1)
where
1
2
1
(
,
,
,
) ,
2
,
,
,
T
k
n
k
k
n
y
W
f
φ
φ
φ
φ
φ
ω
ω
ω
ρ
ρ
ρ
=
=
=
k
φ
is the specific mole number of species
k
(with units [mol/kg]). The density
ρ
and spe-
cific internal energy
e
appear as parameters to
the governing ODE. The total number of spe-
cies in the reaction mechanism is denoted by
n
,
which is also the dimensionality of the chemical
state-space.
The definition of the ILDM is given below.
Given a vector field
f
, the Jacobian matrix is
defined below.
,
,
(
)
,
i j
i j
f
J
J
f
φ
φ
∂
=
=
∂
(2)
The Jacobian can then be transformed into a ba-
sis consisting of a direct sum of two subspaces
5
,
(
)
ˆ
0
(
)
ˆ
0
s
s
s
f
f
f
Z
J
Z Z
Z
φ
Λ
=
Λ
with,
(
)
1
ˆ
ˆ
s
s
f
f
Z
Z
Z
Z
−
=
where
s
Z
is a column partitioning of vectors
spanning the slow subspace, defined to be the
invariant eigenspace of
J
associated with the
least negative eigenvalues. A different basis, in
particular an orthonormal one, can be used in-
stead.
f
Z
is defined analogously, being spanned
by the remaining eigenvectors of
J
.
The equation that defines the ILDM is, finally,
(
)
ˆ
0
f
Z
f
φ
=
(3)
Or more formally, to seek a
k
-dimensional
ILDM, denoted by
k
,we have:
(
)
{
}
ˆ
:
0
k
f
Z
f
φ
φ
=
=
(4)
where
ˆ
:
n
n
k
f
Z
f
−
→
(5)
This definition is numerically awkward for the
following reasons. The ILDM is defined as the
zero level-set of a complicated nonlinear sys-
tem. One-dimensional level-sets are not difficult
to find by continuation methods, but higher di-
mensional level-surfaces are tricky, to put it
mildly.
Next, we have an additional complication from
mass conservation. Each (independent) elemen-
tal constraint increases the multiplicity of the
zero eigenvalue by one. Given
m
(independent)
elemental constraint and assume we are seeking
a
k-
dimensional ILDM, a remedy is to solve for
k
l
o
ij
j
ij
i
N
N
φ
φ
+
∩
=
(6)
where
ij
N
is the (
m
by
n
) species-element ma-
trix, and
o
i
φ
denotes some initial composition.
This is discussed in some detail by Eckett
2
.
Finally, and perhaps most importantly, the de-
fining function
ˆ
f
Z
f
is highly nonlinear. Being
a composite function where
ˆ
f
Z
comes possibly
from matrix inversion, the null space of
ˆ
f
Z
f
,
needed for continuation procedures, is difficult
to find accurately and needs to be approximated.
THE ILDM RECASTED
The disadvantages above notwithstanding, the
ILDM as formulated originally does have an
advantage: it suggests a direct method of solv-
ing for
k
, as long as we can compute level-
surfaces.
By examining Eq. (4), we see that a point
φ
in
chemical state space is in
k
when
(
)
f
φ
,
transformed to the new basis
(
)
s
f
Z Z
, has no
components in the fast subspace. In other words,
we have
{
}
span
k
s
Z
φ
φ
∈
⇔
∈
(7)
Eq. (7) has a very desirable property: only right
eigenvectors are needed to compute a basis for
s
Z
. Although not a
constructive
definition of
k
, Eq. (7) poses, as well as provides an an-
swer to, the important inverse question: What is
the ILDM-dimensionality of
φ
? We will de-
fine the ILDM-dimension of
φ
, denoted by
(
)
ILDM- dim
φ
, by
(
)
ILDM- dim
min( ) :
k
k
φ
φ
=
∈
(8)
The ILDM-dimension is well-defined because
of the (trivial) inclusion property:
0
1
n
⊆
⊆
⊆
(9)
We will see how this ILDM-dimension can be
computed in the next section.
ILDM-DIMENSIONALITY
AND
THE
GRAMMIAN PROCEDURE
We can use the original definition of the ILDM
to compute
(
)
ILDM- dim
φ
by writing
(
)
f
φ
in
a sorted eigenbasis and counting the number of
zeros, but this is numerically ill-posed, in part
because of numerical imprecision and round-off
errors, and more importantly, because
k
is
not an invariant manifold (see the concluding
remarks for more details).
Eq. (8) does provide us with a viable, and direct,
algorithm for estimating
(
)
ILDM- dim
φ
. We
define the
k
-dimensional Gram determinant (or
Grammian
6
) of
φ
, denoted by
(
)
k
φ
Γ
, by
(
)
det(
)
k
T
A
A
φ
Γ
=
(10)
where A is an
(
1)
n
k
×
+
matrix consisting of a
column partitioning of the k slowest eigenvec-
tor, augmented by an arc-length normalized
(
)
f
φ
,
(
)
(
)
1
2
,
,
,
,
n
A
v
v
v
g
φ
=
(11)
where
(
)
(
)
(
)
g
f
f
φ
φ
φ
=
(12)
This definition of
(
)
k
φ
Γ
satisfies the inclusion
relation of Eq. (9):
(
)
(
)
,
i
j
i
j
φ
φ
Γ
≤ Γ
∀
>
(13)
Additionally, the Gram determinant is non-
negative and bounded above by one,
(
)
1
0,
i
i
φ
≥ Γ
≥
∀
(14)
Furthermore,
(
)
k
φ
Γ
, viewed as a scalar valued
function on
(
1)
k
+
vectors, is continuous.
Finally, we have a method of computing
(
)
ILDM- dim
φ
, as follows
(
)
ILDM- dim
min( ) :
k
k
φ
ε
=
Γ
<
(15)
Note that the
(
)
ILDM- dim
φ
in Eq. (15) de-
pends on a parameter
ε
, exactly analogous to
the concept of the numerical rank
7
(
rank
ε
−
)
for matrices.
We will illustrate this technique by computing
the ILDM-dimension along a constant-volume
reaction trajectory. Given an ODE of the form
Eq. (1), subject to some initial conditions
o
φ
,
the reaction trajectory
( )
t
φ
satisfies,
(
)
(0)
( )
( );
,
o
d
t
f
t
e
dt
φ
φ
φ
φ
ρ
=
=
(16)
0
0.5
1
1.5
2
2.5
3
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Specific Mole Number: H
Specific Mole Number: OH
CV reaction trajectory, Hydrogen−Air Combustion
Figure 1. Constant-volume trajectory for Hydrogen-
Air Combustion.
Figure 1 shows a constant-volume trajectory for
the stoichiometric combustion of Hydrogen-Air
(2H
2
+O
2
+3.76N
2
) with a density of 4.58 kg/m
3
and a specific internal energy 1.27 MJ/kg. These
conditions correspond to an initial temperature
of 1543.4 K and an initial pressure of 2.8104
MPa. These conditions are taken from Eckett
2
and correspond approximately to the von Neu-
mann state of a CJ detonation of the mixture.
The first three Gram determinants along the tra-
jectory, i.e.
(
)
( ) ,
1, 2, 3
m
t
m
φ
Γ
=
are shown in
Figure 2. This figure has two interesting inter-
pretations.
0
0.5
1
1.5
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time (
μ
s)
Grammian
Γ
3
(t)
Γ
1
(t)
Γ
2
(t)
Figure 2: The first three Gram Determinants along the
CV trajectory of Figure 1 are plotted as a function of
time.
It can be seen, for example at 1 microsecond,
the only non-zero Gram determinant is
1
Γ
. This
follows from Eq. (15) that the ILDM-
dimensionality of the trajectory at that instant is
2. It means that when
2
ε
Γ
≤
and
1
ε
Γ
>
, the
two slowest eigenvectors are necessary and suf-
ficient to span
(
)
f
φ
.
The alternate point of view is the concept of the
time of arrival, introduced below. We will de-
fine the
time of arrival
k
t
of a trajectory
( )
t
φ
to
k
by
(
)
(
)
min
:
( )
k
k
k
t
t
t
t
φ
ε
Γ
≤
∀
>
(17)
It can be seen from Figure 2, with the definition
provided by Eq. (17), that the time of arrival to
3
2
1
,
,
is, approximately, 0.4, 0.5 and 1.5
microseconds, respectively.
Figure 3 decorates the reaction trajectory from
Figure 1 with the ILDM-dimension estimated by
the Grammian procedure. As will be seen
shortly, this ability to compute the ILDM-
dimension of any point
φ
in chemical state-
space forms the basis of our proposed embar-
rassingly parallel algorithm for ILDM computa-
tions.
0
0.5
1
1.5
2
2.5
3
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Specific Mole Number: H
Specific Mole Number: OH
Equilibrium
1D
2D
3D
4D and up
Figure 3: The ILDM dimension is shown along the
CV trajectory of Figure 1.
ILDM VIA CONGRUENCES
The definition of ILDM through Eq. (4) pro-
vides, by rewriting it in explicit form, a direct
method of computing
,
(
)
(
)
1
ˆ
0
k
f
Z
f
−
=
(18)
It is important to note what is meant by a (nu-
merical) solution of Eq. (18). Generally, a solu-
tion is comprised of a set of points
S
, perhaps
tabulated with some predetermined parameteri-
zations, in the original
n
-dimensional chemical
state space, with each point satisfying Eq. (3) to
within some tolerance:
(
)
ˆ
,
f
Z
f
S
φ
ε
φ
≤
∀
∈
(19)
Numerical schemes based on continuation
methods are typically used to solve Eq. (18),
which becomes very difficult for
1
k
>
.
On the other hand, armed with the tools from
the previous sections for computing the ILDM-
dimension, it is possible to take an indirect ap-
proach. We will redefine the set
S
as:
(
)
,
k
S
φ
ε
φ
Γ
≤
∀
∈
(20)
Given a trajectory
( )
t
φ
satisfying Eq. (1) sub-
ject to some initial condition
o
φ
, we compute
the time of arrival of
( )
t
φ
to
k
t
to
k
, as de-
fined in Eq. (17). Together with Eq. (20), we
see that
(
)
,
k
t
S
t
t
φ
∈
∀
≥
(21)
In other words, we can solve
k
(populate the
set
S
) by solving Eq. (1) subjected to different
initial conditions. This “filling of a manifold”
by curves is called a congruence
8
.
EXAMPLE: CONSTRUCTION OF A
ONE-DIMENSIONAL ILDM
We will use our procedure described above to
compute the one-dimensional ILDM for the
detonation problem studied in Rastigejev et al
9
.
We will represent the system by Eq. (1), re-
peated below,
(
)
;
,
d
f
e
dt
φ
φ
ρ
=
(22)
We will proceed as follows. A one-dimensional
ILDM
1
is, roughly speaking, a line through
the equilibrium point in chemical state-space.
Using any point
1
o
φ
∈
as initial data to Eq.
(22), the trajectory that results will be a part of
1
, starting at
o
φ
and ending at equilibrium.
The equilibrium point therefore divides
1
into
two pieces. Our procedure that follows will de-
termine, in the linear approximation, the two
initial data (one on each side of the equilibrium)
that maximize the extent of our solution to
1
.
We first compute the one-dimensional slow ei-
genspace of the system, which is an affine linear
space centered at the equilibrium point spanned
by the slowest eigenvector. This is shown in
Figure 4.
The meaning of this eigenspace, which we will
denote by
1
L
for convenience, is well known
and will not be elaborated. In the neighborhood
of the equilibrium point,
1
is tangent to
1
L
.
0
1
2
3
4
5
0
1
2
3
4
5
x 10
−4
Specific Mole Number: H2O
Specific Mole Number: H2O2
Reaction Trajectory
Equilibrium Point
Slowest Affine Eigenspace
Figure 4 CV Reaction trajectory, equilibrium point
and the slowest linear eigenspace.
Chemical state-space is compact as it is sub-
jected to the positivity constraint
0,
1..
i
i
n
φ
≥
∀
=
(23)
as well as elemental conservation
o
ij
j
ij
j
N
N
φ
φ
=
(24)
where, as alluded to earlier,
ij
N
is the species-
element matrix. Eq. (24) is automatically en-
forced in the solution of Eq. (22) because, when
each of the reactions in the detailed mechanism
conserves elements,
(
)
(
)
Null
ij
f
N
φ
∈
(25)
The maximum extent of
1
L
is a problem in lin-
ear programming. In one-dimension, it is akin to
extending
1
L
until one of the
n
inequalities in
Eq. (23) becomes an equality. The range of va-
lidity for the thermodynamic data imposes an
additional constraint: at a given value of density
and internal energy,
(
)
min
max
;
,
T
T
e
T
φ
ρ
≤
≤
(26)
The maximal linear eigenspace as well as its
thermodynamic boundary, determined by Eq.
(26), is illustrated in Figure 5.