Combustion and Flame 144 (2006) 2436

www.elsevier.com/locate/combustflame

Linear time reduction of large kinetic mechanisms with

directed relation graph: n-Heptane and iso-octane
Tianfeng Lu , Chung K. Law
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA
Received 7 September 2004; received in revised form 30 January 2005; accepted 10 February 2005
Available online 30 August 2005

Abstract
The algorithm of directed relation graph recently developed for skeletal mechanism reduction was extended to
overall linear time operation, thereby greatly facilitating the computational effort in mechanism reduction, partic-
ularly for those involving large mechanisms. Together with a two-stage reduction strategy and using the kinetic
responses of autoignition and perfectly stirred reactor (PSR) with extensive parametric variations as the crite-
ria in eliminating unimportant species, a detailed 561-species n-heptane mechanism and a detailed 857-species
iso-octane mechanism were successfully reduced to skeletal mechanisms consisting of 188 and 233 species,
respectively. These skeletal mechanisms were demonstrated to mimic well the performance of the detailed mech-
anisms, not only for the autoignition and PSR systems based on which the reduced mechanisms were developed
but also for the independent system of jet-stirred reactor. It was further observed that the accuracy of calculated
species concentrations was equivalently bounded by the user-specified error threshold value and that the reduction
time for a single reaction state is only about 50 ms for the large iso-octane mechanism.
2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

Keywords: Directed relation graph; Reduced mechanism; n-Heptane; Iso-octane

1. Introduction significant stiffness to the governing equations due

to the dramatic differences in the time scales of the
Chemical kinetics of fuel oxidation is an inte- species. Consequently, there exists the need to de-
gral component of combustion phenomena. However, velop from these detailed mechanisms corresponding
because of the complexity and strong nonlinearity reduced mechanisms of fewer variables and moder-
of the reaction pathways and the associated reaction ated stiffness, while maintaining the accuracy and
rates for the myriad of reactions and species involved, comprehensiveness of the detailed mechanism.
computational simulations with detailed mechanisms A major category of previous reduction methods
have been limited to the simplest fuels such as hy- was based on time scale analysis, through which both
drogen and the lower hydrocarbons in idealized flow the number of variables and the stiffness can be re-
fields. The simulation is further complicated by the duced by eliminating short time scales associated with
existence of highly reactive radicals which induces
quasi-steady-state species or partial equilibrium reac-
tions [18]. The method of intrinsic low-dimensional
* Corresponding author. Fax: +1 609 258 6233. manifold [9] performs eigenvalue analysis of the Ja-
E-mail address: tlu@princeton.edu (T. Lu). cobian matrix and assumes that the fast subspace van-
0010-2180/$ see front matter 2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
doi:10.1016/j.combustflame.2005.02.015
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 25

ishes quickly. The frequently reused reaction condi- species is formulated with the entries of the Jaco-
tions are handled with the method of in situ adap- bian matrix and the species that are strongly cou-
tive tabulation to reduce the overall simulation cost pled with the major species are identified iteratively.
[10,11]. The theory of computational singular pertur- This method is practically faster than those based
bation (CSP) considers the time dependence of the on sensitivity analysis, although formulation of the
Jacobian matrix, and higher-order accuracy can be ob- species coupling using Jacobian matrix is nontrivial
tained [1216]. and system-dependent knowledge is frequently re-
While time scale analyses can reduce the num- quired.
ber of species and the stiffness by eliminating species In a recent work [24], the method of directed re-
or reactions with short time scales, significant over- lation graph (DRG) was developed to identify unim-
head in central processing unit (CPU) cost is typically portant species by resolving species coupling with
involved in solving the algebraic equations obtained high efficiency and minimal requirement of system-
by assuming that the fast species or reactions are ex- dependent knowledge. In the present investigation
hausted, particularly when the size of the mechanism this method was extended to the reduction of large
is large. Therefore, it is frequently necessary to iden- mechanisms with hundreds of species, namely those
tify and eliminate the unimportant species and reac- for n-heptane and iso-octane [25,26], and the effi-
tions, i.e., to perform the skeletal reduction, before ciency and accuracy of this extended method were
reduction based on time scale analysis is conducted. demonstrated and analyzed.
Various methods for skeletal reduction have been
developed. One of the earliest developed methods
for skeletal reduction is sensitivity analysis [1719], 2. Methodology
which is simple to apply. However, it does not di-
rectly provide decoupled information about the reac- 2.1. DRG method
tions and species, and further postprocessing is re-
quired. The method of principal component analysis The skeletal reduction is to identify and elimi-
(PCA), based on sensitivity analysis, operates on the nate unimportant reactions or species. Typically, the
sensitivity matrices and systematically identifies the identification of unimportant species is more involved
redundant reactions [20]. A major restriction of algo- than that of reactions due to the complex coupling of
rithms based on sensitivity analysis is that they are the species. The method of DRG was developed to
time consuming to compute the sensitivity matrices resolve species coupling efficiently, by starting with
for large mechanisms. direct species coupling, which indicates that the re-
Reaction elimination can also be performed using moval of one species B from the mechanism induces
optimization approaches, such as integer program- immediate error to the production rate of another
ming, aiming to obtain an optimal set of reactions for species A. Such immediate error, noted as rAB , can
given constraints [21]. While the optimal set of re- be expressed as
actions can be identified, the optimal solution strictly

depends on the selection of constraints, which is fre- i=1,I |A,i i Bi |
quently quite involved. Furthermore, the optimiza- rAB
,
i=1,I |A,i i |
tion approach is asymptotically slower than sensitiv-  1, if the ith reaction involves
ity analysis, since integer programming is in general
Bi = species B, (1)
an NP-hard problem. Another method, the detailed re-
0, otherwise,
duction [22], provides a fast way for the identification
of unimportant reactions by directly comparing the where the subscripts A and B indicate the species
reaction rates with a preselected critical value. The identity, i is the ith reaction, A,i the stoichiomet-
restriction of this method is the lack of consideration ric coefficient of species A in the ith reaction, and
of species or reaction coupling, because of which a the reaction rate. The terms in the denominator are the
slower reaction is not always unimportant, particu- contributions of the reactions to the production rate of
larly when it involves crucial radicals. Consequently, species A, and the terms in the numerator are those in
significant human interaction and extensive validation the denominator that involve species B. If the relative
may be required. error rAB is not small compared to a threshold value
While skeletal reduction based on reaction elim- , the removal of species B from the skeletal mecha-
ination has been extensively studied, those methods nism immediately induces a nonnegligible error in A.
involving direct species elimination are less estab- Consequently species B should be kept in the skele-
lished. A method was developed for the identification tal mechanism if species A is to be retained. Such an
of unimportant species by resolving species coupling immediate requirement of species A to species B is
using Jacobian analysis [23]. The coupling among denoted as A B.
26 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

Species can also be coupled indirectly through di-

rect coupling with intermediate species. For example,
if A B and BC, then A requires C indirectly.
The set of species required either directly or indirectly
by species A is defined as the dependent set of A,
and it is readily seen that if species A is to be re-
tained in the skeletal mechanism, so should be all the
species in its dependent set. Therefore, for a given set
of major species which must be retained in the skele-
tal mechanism, the species of the skeletal mechanism (a)
can be obtained from the union of all the dependent
sets of the major species. The dependent set of a major
species can be identified efficiently by mapping the
species coupling to a directed relation graph, which
can be constructed using the following rules:

(1) Each vertex in DRG is uniquely mapped to a

species in the detailed mechanism.
(2) There exists a directed edge A B if and only
if rAB
.
(3) The starting vertices of DRG correspond to the
major species in the mechanism. (b)

Fig. 1. Schematics showing typical configuration of the di-

A sample configuration of DRG is shown in
rected relation graph (DRG); the vertices correspond to the
Fig. 1a, in which vertices A, B, . . . , F correspond
species and the directed edges correspond to the immediate
to the respective species, and each directed edge in- requirement of one species to another. (a) Basic configura-
dicates an immediate requirement of one species of tion and identification of skeletal species. (b) Demonstration
another. The starting vertex A, enclosed in a bold of strongly coupled species groups.
circle, is a major species in the mechanism. The skele-
tal species in Fig. 1a correspond to all the vertices
reachable from the starting vertex. Graph searching species are removed and the skeletal mechanism is
methods such as depth first search (DFS) can thus be generated.
exploited to efficiently find all the vertices reachable An important observation of skeletal reduction
from the starting one. Species A, B, and D are there- with DRG is the existence of strongly coupled species
fore identified as the species of the skeletal mecha- groups, which are characterized by relatively stronger
nism in the current example. intragroup coupling than intergroup coupling. A sche-
It is noted that the set of starting vertices does not matic of strongly coupled species groups is shown
need to contain every major species in the mecha- in Fig. 1b. For example, the intragroup couplings in
nism. Instead, in most cases, it is sufficient to start the species groups {B, D} and {E, F} are stronger
from a single species, such as the fuel, to identify than the intergroup couplings. Such strongly coupled
all the species in the skeletal mechanism. This is be- species groups should be either retained or elimi-
cause the other major species such as the oxidizer, the nated together from the skeletal mechanism. Conse-
products, and the major radicals are required by the quently, the elimination of a strongly coupled species
fuel either directly or indirectly, and as such they can group induces a jump in the number of species in the
be retained together with the other important species skeletal mechanism as the error threshold value in-
automatically in the DFS. The inert species can be ei- creases, as was shown in the reduction of a detailed
ther added back to the skeletal species set after the ethylene mechanism [24]. Such jumps in the species
DRG reduction since they will be eliminated in the number can significantly reduce the possible number
DRG reduction or included as a starting species for of candidate skeletal mechanisms and thereby reduce
DRG reduction such that they can be retained conse- the amount of effort for skeletal mechanism valida-
quently. In cases in which multiple starting vertices tion.
are required, a DFS can be applied to each starting To obtain a skeletal mechanism over a sufficiently
vertex and, as mentioned above, the collection of all wide range of parameters such as pressure, temper-
the discovered vertices constitutes the species set of ature, equivalence ratio, and residence time, a set of
the skeletal mechanism. After the skeletal species are reaction states are sampled from typical applications
identified, reactions involving any of the eliminated such as perfectly stirred reactor (PSR) and autoigni-
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 27

tion. A local DRG reduction can be performed for each elementary reaction is bounded by a small con-
each of these reaction states, and the union of the stant number Smax , say 68. Thus graph construction
species sets of all the local skeletal mechanisms yields can be performed in a more efficient way by evaluat-
the species set of the global skeletal mechanism. ing only the contributions of each reaction to the re-
It is to be noted that since the comprehensive lated edges. Since each reaction involves at most Smax
skeletal mechanism developed here is the union of a species, the maximum number of edges in the graph
set of local skeletal mechanisms each developed for a is bounded by a linear function of I and so is the
reaction state, many partially comprehensive skeletal time complexity of the graph construction step. The
mechanisms can be developed through the union of detailed procedure for graph construction is shown in
any number of these local skeletal mechanisms, with the coarse-grain pseudo code in Fig. 2.
different sizes, different parametric ranges, and dif- In the second step, each edge in the graph can
ferent extents of computational cost reduction. The be compared with the threshold value for trunca-
methodology is therefore applicable to the develop- tion. The resulting graph with a smaller number of
ment of mechanisms of various degrees of compre- edges can be searched with DFS to obtain the final
hensiveness. Indeed, each of the local skeletal mecha- species set of the skeletal mechanism. Since the num-
nisms being applicable to a given reaction state allows ber of edges is bounded by a linear function of I ,
the ultimate potential of transition between different and DFS features a linear searching time proportional
skeletal mechanisms during a simulation, from state to the number of edges in the graph, the time com-
to statea promising prospect that is worth further plexity of the second step is also O(I ). While this
pursuing. standard graph searching method is fast for reduction
with given accuracy requirement, it is not sufficient
2.2. Extended DRG algorithm for linear time if the number of species, instead of the accuracy, is
reduction given for the reduction. In such cases, the DFS algo-
rithm can be revised such that species can be sorted
The major advantage of DRG reduction is that the by the critical value of .
entire reduction can be completed in a time that is The revised DFS (RDFS) method first sorts the
linearly proportional to the number of reactions in edges by their respective values r. Then the edges
the mechanism. Linear time reduction is crucial for are inserted into the initially empty graph in decreas-
applications requiring real time or case-specific re- ing order of r, such that the strongly coupled species
duction. In CFD simulations using explicit integration groups, or kernels, can be formed first. The kernels
schemes, real time reduction based on Jacobian or grow larger as more edges are inserted into the graph,
sensitivity matrix analysis requires significant over- until a bridging edge is inserted and an otherwise iso-
head, which consequently diminishes the efficiency lated kernel becomes reachable from the starting ver-
of the reduction. In such cases, a reduction algorithm tices. The r value of this bridging edge then becomes
with linear reduction time is required to minimize the the critical value for each species in this kernel. The
overhead and provide the highest efficiency. DRG is procedure continues until all the edges are inserted,
such an algorithm in that the entire reduction time is as described in Fig. 2. While general purpose sort-
linearly proportional to the number of reactions in the ing algorithms take O(I log I ) time to sort the edges,
mechanism, as will be shown in the following sec- a simple linear time bucket sorting algorithm can be
tion. applied in the RDFS algorithm, considering that the
There are three major steps in DRG reduction, possible value of an edge is between 0 and 1, and
namely the graph construction step using Eq. (1), the the r value needs to be accurate only to 103 104 .
graph searching step to identify the skeletal species, The graph searching time for the RDFS algorithm is
and the skeletal mechanism generation step to create therefore still a linear function of I by the argument
the final mechanism files. that the total number of edges is O(I ) and no edge
The graph construction step involves evaluation needs to be searched more than once. Compared with
of the relative errors for each pair of species ac- the regular DRG reduction, the extended DRG pro-
cording to Eq. (1). There are totally K 2 terms to be vides the sorted species list in a single linear time
evaluated, where K is the number of species, and RDFS search. Therefore the skeletal mechanism of
each term contains O(I ) multiplications. Therefore any given accuracy or number of species can be gen-
a brute-force approach in graph construction would erated using this list.
require O(K 2 I ) evaluations, which can be asymp- Finally, the third step of the DRG reduction is
totically more expensive than the evaluation of the to eliminate the reactions that involve any removed
Jacobian matrix. However, the stoichiometric coeffi- species and subsequently generate the final skeletal
cient matrices for the detailed mechanisms are sparse mechanism. This step is straightforward and it is triv-
in general; i.e., the number of species involved in ial to show that the time complexity is also O(I ),
28 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

Fig. 2. Coarse-grain pseudo code of the extended DRG method with linear time graph construction and revised depth first search
(RDFS) algorithm.

provided that the number of species in each reaction formed on the skeletal mechanism. In such cases,
is bounded by a constant value, Smax . a two-stage DRG reduction can produce a skeletal
Since each of the three steps in DRG reduction can mechanism slightly smaller than that from a single-
be performed in linear time, the overall reduction time stage DRG reduction. The first stage of DRG reduc-
for DRG is O(I ), i.e., a linear function of the number tion is the major reduction stage, in which a signif-
of reactions in the detailed mechanism. icant number of species are eliminated from the de-
tailed mechanism, and the second stage is a minor
2.3. Two-stage DRG for reduction of large stage whose execution is optional.
mechanisms While it might be argued that more than two
DRG stages are needed to obtain a converged skele-
When the number of species is large in the de- tal species set for extremely large mechanisms, it is
tailed mechanism, the number of species eliminated reasonable to expect, and indeed was found for the
by a DRG reduction is expected to be correspondingly reduction of the largest mechanisms currently avail-
large. Consequently, the r value for each pair of the able such as those of n-heptane and iso-octane [25,
retained species can be changed due to the elimina- 26], that a two-stage DRG reduction is quite adequate.
tion of the unimportant species. Although the change It is further noted that, since the second stage of DRG
in the r values is a second-order effect, it is still pos- reduction can also be performed in linear time and
sible that some species in the skeletal mechanism can since the number of reactions in the skeletal mecha-
be further eliminated if a second-stage DRG is per- nism determined from the first DRG reduction stage is
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 29

typically much smaller than that of the detailed mech-

anism, the inclusion of an additional DRG stage does
not increase the overall reduction time significantly.
Consequently, a slightly smaller skeletal mechanism
can be obtained by implementing the second stage
with a slightly increased reduction time. Indeed, con-
sidering the minuscule amount of time involved in the
reduction, as will be shown later, the implementation
of the second stage of reduction is entirely appropri-
ate.
With the above extended DRG method, skele-
tal mechanisms with desired accuracy or size can
be generated. While such a skeletal mechanism al-
ready has reduced number of species and reactions,
there still could be some unimportant reactions that
contribute negligibly to every species in the skeletal
mechanism. To identify such unimportant reactions,
methods based on reaction elimination discussed in
the Introduction can be used. For example, a reliable Fig. 3. Comparison of ignition delay around NTC zone for
and efficient method using the CSP importance index, stoichiometric n-heptane mixture at atmospheric pressure
which measures the normalized contribution of a re- calculated with different versions of detailed mechanisms
action to the production rate of a species, can be found and a reduced mechanism, showing the large discrepancy of
in [14,15]. Versions 1 and 2 of the detailed mechanism and the necessity
for updating the skeletal mechanism.
In the next sections, the DRG method will be ap-
plied to derive skeletal mechanisms of two large de-
tailed mechanisms for n-heptane and iso-octane. The in previous works [3335]. However, the reduction
efficiency and accuracy control will be demonstrated processes using traditional reduction methods such as
and analyzed. sensitivity analysis are typically time consuming. Fur-
thermore, continuous discovery of new information
on the reaction pathways results in the frequent updat-
3. Skeletal mechanisms for n-heptane and ing of the detailed mechanisms, which in turn neces-
iso-octane sitate corresponding updating of the reduced mecha-
nisms. We demonstrate this point by showing in Fig. 3
3.1. Skeletal mechanism for n-heptane that while the skeletal mechanism in [34] agrees well
with Version 1 of the detailed mechanism [25], a large
n-Heptane is an important reference gasoline fuel, discrepancy in the computed ignition delay time is
and many detailed mechanisms for its oxidation have observed when the detailed mechanism is updated to
been developed by different groups [2531]. Due to Version 2 [26,27], hence indicating the need for the
the large molecule size, the oxidation process for n- development of a new skeletal mechanism for the up-
heptane is significantly more complex than that of dated detailed mechanism. In such a case, the DRG
lower hydrocarbons such as methane and ethylene. method can be applied for the fast reduction of large
Furthermore, n-heptane exhibits the negative temper- mechanisms whenever they are updated.
ature coefficient (NTC) and the two-stage or multi- We have thus applied the two-stage DRG method
stage ignition [32] behaviors in the low-temperature with the RDFS algorithm to Version 2 of the 2002
regime, whose descriptions require additional sets of detailed n-heptane mechanism [25,26], posted on the
reaction pathways. Therefore, the size of the detailed LLNL website [27]. This detailed mechanism consists
n-heptane mechanisms is typically much larger than of 561 species and 2539 reactions. Typical working
that of the smaller hydrocarbons. For example, the de- conditions were sampled from PSR and autoignition
tailed mechanism developed by Curran et al. [2527] with the equivalence ratio ranging from 0.5 to 1.5, and
contains more than 500 species and 2000 elementary the pressure ranging from 1 to 40 atm. The initial tem-
reactions. While detailed mechanisms of such a large perature for PSR was fixed at 300 K and the residence
size can still be used in the simulation of many homo- time range was from the state of extinction to that
geneous combustion phenomena, they can rarely be of near equilibrium. The initial temperature range for
applied to the simulation of diffusive systems with- autoignition was from 600 to 1800 K, which covers
out significant reduction in size. Indeed, several re- both the low-temperature NTC regime and the high-
duced mechanisms of n-heptane have been developed temperature ignition regime. More than 1000 points
30 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

Fig. 4. Dependence of number of species in the skeletal

mechanism as a function of the user-specified error thresh-
old value for the first and second stages of global DRG
reduction for n-heptane mixture, based on densely sampled
reactions states from PSR and autoignition for pressure from
1 to 40 atm, equivalence ratio from 0.5 to 1.5, initial temper-
ature range from 600 to 1800 K for autoignition, and 300 K
for PSR.

were sampled to ensure the coverage of every rep-

resentative reaction condition. The entire reduction
Fig. 5. Comparison of calculated ignition delay of n-
process using these sample reaction states took a few
heptane/air mixture autoignition as a function of initial tem-
seconds, which is negligible compared to the simula- perature under constant pressure and constant enthalpy with
tion time. detailed and 188-species skeletal mechanisms, for equiva-
The number of species in the skeletal mechanism lence ratio from 0.5 to 1.5, pressure from 1 to 40 atm, and
as a function of the critical relative error for the first initial temperature from 600 to 1800 K.
and second DRG reduction stages is shown in Fig. 4.
It is seen that, from the first stage of reduction, the ous reduction [34], which has 159 species and 770
number of species decreases rapidly as increases reactions.
when is small, say less than 0.2, while further in- The 188-species skeletal mechanism was vali-
creasing the value of leads only to slight reductions dated against the detailed one for PSR and autoigni-
in the species number. This trend implies that sig- tion within the parameter range of the reduction.
nificant numbers of species can be eliminated from Fig. 5 shows the calculated ignition delay time of
the detailed mechanism even with a small value, heptane/air mixtures as a function of the initial tem-
i.e., with high accuracy. By choosing an = 0.1, perature for = 0.51.5 and p = 140 atm. It is seen
290 species are retained after the first stage of reduc- that the agreement between the skeletal and the de-
tion. It may be noted that, because of the large number tailed mechanisms is good for both the NTC and
of total species involved, jumps in the species num- the high-temperature regimes. Fig. 6 compares the
ber representing the elimination of strongly coupled calculated temperature profiles as a function of the
groups consisting of 10 to 20 species are not sharply residence time for PSR over the same equivalence
contrasted on the graph. ratio and pressure ranges. Again, good agreement is
The second DRG reduction was then conducted observed. It is noted that, while the unstable middle
based on the 290-species first-stage skeletal mecha- branches of the extinction curve were difficult to ob-
nism. Facilitated by the species number jumps in the tain using the detailed mechanism, the convergence
reduction curve at this stage, 0.19 was selected can be readily achieved using the skeletal mechanism
and a 188-species skeletal mechanism was generated. because of the significantly smaller number of system
After further elimination of unimportant reactions us- variables. This facilitation in obtaining converged so-
ing the CSP importance index, a final skeletal mech- lutions using the skeletal mechanism for applications
anism with 842 reactions was obtained. The size of based on the Newton solver is an extra benefit in addi-
this skeletal mechanism is close to that in the previ- tion to the significantly reduced simulation time. The
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 31

Fig. 7. Dependence of number of species in the skeletal

mechanism as a function of the user-specified error thresh-
old value for the first and second stages of global DRG
reduction for iso-octane mixture, based on densely sampled
reactions states from PSR and autoignition for pressure from
1 to 40 atm, equivalence ratio from 0.5 to 1.5, initial temper-
ature range from 600 to 1800 K for autoignition, and 300 K
Fig. 6. Comparison of calculated temperature profile in PSR for PSR.
as a function of the residence time for n-heptane/air with the
detailed and the 188-species skeletal mechanisms, for equiv-
alence ratio from 0.5 to 1.5, pressure from 1 to 40 atm, and
to be retained as a function of the critical value of
initial temperature 300 K. for the first and second reduction stages. By select-
ing = 0.1 for the first stage, a skeletal mechanism
with 347 species was obtained for further reduction
good agreement in the comprehensive validation for
through the second stage, which subsequently yields
PSR and autoignition demonstrates the validity of the
a skeletal mechanism with 233 species and 959 re-
DRG reduction algorithm.
actions. The comprehensive validation of the skeletal
3.2. Skeletal mechanism for iso-octane mechanism for PSR and autoignition are shown in
Figs. 8 and 9, respectively. Close agreement is again
Iso-octane is another important component in fuel observed.
blends. It has also been extensively studied experi-
mentally, and many versions of mechanisms for iso-
octane oxidation have been developed [26]. While 4. Further validation
the ignition phenomena for iso-octane bear substan-
tial similarity to those of n-heptane, such as the cool While the accuracy of the skeletal mechanisms re-
flame and NTC behaviors, the difference in the mole- duced using DRG is bounded by the threshold error
cule structure renders the detailed mechanisms for value , the comprehensiveness is determined by the
these two large hydrocarbons significantly different parameter range and the phenomena based on which
from each other. Furthermore, the size of the detailed the mechanisms are developed. Thus in comprehen-
mechanism for iso-octane is expected to be larger sive reduction it is necessary to cover the ranges
than that of n-heptane. For example, the detailed 2002 of temperature, pressure, and reactant concentrations
mechanism by Curran et al. [26], downloadable from as extensively as possible, as has been done in the
[36], has 857 species and 3606 reactions. While a re- present study. Furthermore, the selected phenomena
duced mechanism for the iso-octane flame structure are expected to be relevant to more complex phenom-
has been developed based on a much smaller mecha- ena which consist of subregimes with reaction states
nism [37], the comprehensive reduction based on the similar to those covered in the reduction stage. Con-
large mechanism [26,36] is formidable using tradi- sequently, we have selected PSR and autoignition,
tional approaches. respectively, to capture the high-temperature burning
Using the same parameter range and sample ap- and extinction situations and the low-temperature ig-
plications as those for n-heptane, the detailed iso- nition situations. The reduced mechanisms generated
octane mechanism was reduced with the two-stage were shown to mimic the detailed mechanisms well
DRG method. Fig. 7 shows the number of species for these two phenomena.
32 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

Fig. 9. Comparison of calculated temperature profile in PSR

as a function of the residence time for iso-octane/air with the
detailed and the 233-species skeletal mechanisms, for equiv-
alence ratio from 0.5 to 1.5, pressure from 1 to 40 atm, and
initial temperature 300 K.
Fig. 8. Comparison of calculated ignition delay for iso-
octane/air mixture autoignition as a function of initial tem-
perature under constant pressure and constant enthalpy with
with the detailed ones, with only slightly larger devi-
the detailed and the 233-species skeletal mechanisms, for
equivalence ratio from 0.5 to 1.5, pressure from 1 to 40 atm,
ations in a few situations with maximum errors less
and initial temperature from 600 to 1800 K. than 15%, which is consistent with the assigned ac-
curacy . The extension of the skeletal mechanisms
The final step in the development of the reduced developed based on PSR and autoignition to JSR can
mechanism is then to apply the reduced mechanisms therefore be considered successful.
to phenomena other than PSR and autoignition, for The fidelity of the skeletal mechanisms can be
an independent validation. In our previous studies further scrutinized by examining the species profiles.
on ethylene [24], we have extended the validation Fig. 12 shows the calculated mole fractions of sev-
of the performance of the reduced mechanisms to eral species as functions of initial temperature of the
the calculation of diffusive systems involving laminar iso-octane mixture using the 233-species skeletal and
flame speeds and nonpremixed counterflow ignition. the detailed mechanisms in a JSR, for = 1, 1.5,
Such a comparison, however, is formidable due to p = 10 atm, and residence time = 1 s. It is seen that
the heavy demand in CPU time for the present large the comparison is mostly very close, with the largest
mechanisms. As an alternative, we have selected the error for both the temperature and the species profiles
response of the jet-stirred reactor (JSR) for an in- occurring for the fuel-rich case within the temperature
dependent validation. This system covers both low- range from 900 to 1000 K. The largest error for the
and moderately high-temperature situations, with the species profiles is about 20% for CO2 around 900 K,
presence of combustion-generated radicals to initiate which is still comparable to the error in the tempera-
ignition as compared to the need for self-generation ture profile and the critical value of .
in autoignition. As a final example of validation comparison,
Figs. 10 and 11 show the temperature rise in JSR Figs. 13a and 13b, respectively, show the evolution
as a function of the initial temperature of the mix- of the mole fractions of selected major and radical
ture for stoichiometric 0.1% n-heptane and 0.1% iso- species for a stoichiometric 0.1% iso-octane mixture
octane mixtures, respectively, with = 0.51.5, p = in an atmospheric-pressure flow reactor, with initial
1040 atm, and residence time = 1 s. It is seen temperature 1080 K. Very close agreement is again
that the skeletal mechanisms mostly agree very well observed.
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 33

Fig. 10. Comparison of calculated temperature rise in JSR Fig. 11. Comparison of calculated temperature rise in JSR
as a function of the initial temperature for 0.1% n-heptane as a function of the initial temperature for 0.1% iso-octane
with the detailed and the 188-species skeletal mechanisms, with the detailed and the 233-species skeletal mechanisms,
for equivalence ratio from 0.5 to 1.5, pressure from 10 to 40 for equivalence ratio from 0.5 to 1.5, pressure from 10 to
atm, and residence time = 1 s. 40 atm, and residence time = 1 s.

5. Error and reduction time

based on reaction states sampled only from the flow
It is noticed from Figs. 1113 that an important reactor configuration. To demonstrate the differences
property of the DRG algorithm is that the relative between global and local reductions, Fig. 14 shows
error for each species and major system parameters the relative error in the temperature rise in a flow
such as the temperature is equally bounded by the reactor for a 0.1% stoichiometric iso-octane mixture
user-specified critical value ; i.e., every species re- under atmospheric pressure with initial temperature
tained in the skeletal mechanism developed by the 1080 K and residence time = 0.2 s. The circles and
DRG method should have an equivalent relative er- triangles, respectively, show the relative errors for lo-
ror bound, such that all the species concentrations in cal and global reductions. For reference, the fractions
the solution should be equivalently trustable. The rea- of eliminated species as functions of the threshold
son for this merit of the DRG method is that since value for both the local and the global reductions are
every species in the dependent set of a species A has shown by the solid curves. It is seen that the relative
to be retained in the skeletal mechanism if species A error for local reduction grows in a linear trend as
is kept, as discussed in Section 2.1, the accuracy for increases, and the line 10 bounds all the local rel-
every such species A in the skeletal mechanism is ative errors. However, the errors for global reduction
therefore bounded by . are much smaller than those of the local reduction,
In developing the reduced mechanism we have indicating that, by retaining many species that are not
required that all sampled states satisfy the critical important for the local conditions, the relative errors
value . The mechanism so derived is therefore com- for the global reduction are reduced. This effect can
prehensive in its applications, albeit at the expense be clearly seen from the two solid curves: the frac-
of keeping more species. It is therefore reasonable to tion of eliminated species for the global reduction is
expect that, at the expense of losing some degree of much smaller than that for the local reduction with the
comprehensiveness, a greater extent of reduction can same threshold value . Therefore, the accuracy con-
be achieved by performing the reduction locally, say trol scheme in DRG reduction is more clearly demon-
34 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

Fig. 13. Comparison of calculated species profiles in a

flow reactor as a function of the residence time for 0.1%
iso-octane with the detailed and the 233-species skeletal
Fig. 12. Comparison of calculated species profiles in JSR mechanisms, for = 1.0, p = 1 atm, and initial tempera-
as a function of the initial temperature for 0.1% iso-octane ture 1080 K.
with the detailed and the 188-species skeletal mechanisms,
for = 1.0, p = 10 atm, and = 1 s. discussions in Section 2.2. Furthermore, the reduc-
tion time for the large iso-octane mechanism with 857
strated in the local reduction than in the global reduc- species takes only about 50 ms, which is negligible
tion. As a result, the value describes relative errors compared to the simulation time and thereby demon-
better in local reduction than in global reduction. Con- strates the efficiency of the DRG reduction. Fig. 15b
sequently, the value of should be selected conser- shows that the reduction time also varies quite linearly
vatively for the local reduction and aggressively for with the number of species in the mechanism. This is
the global reduction. For example, to obtain a skeletal possibly a consequence of the practical observation
mechanism using local DRG reduction with a rela- that the number of reactions in large detailed mech-
tive error tolerance of 10%, a conservative value of anisms is approximately linearly proportional to the
1% should be used. On the other hand, by using number of species.
the global reduction with the same error tolerance, an
aggressive 1020% can be selected.
We next investigate the actual time involved in 6. Conclusion
processing the reduction. Fig. 15 shows the measured
average local reduction time for a single reaction state In the present study we have significantly im-
as a function of the size of the mechanism, i.e., the proved the capability of the DRG method for skeletal
number of reactions or number of species in the mech- mechanism reduction. Specifically, overall linear time
anism. Each data point was averaged from more than reduction is achieved by employing a linear time al-
500 local reductions using the reaction states sam- gorithm for graph construction and a revised deep first
pled from PSR and autoignition over a parametric graph search algorithm. The RDFS algorithm gen-
range similar to that for the comprehensive reduction erates all possible skeletal mechanisms ordered by
for n-heptane and iso-octane in Section 3. Detailed accuracy by a single linear time searching and can
and skeletal mechanisms of various sizes for differ- be applied for reductions requiring fixed numbers of
ent fuels were taken as the starting mechanisms for species in the skeletal mechanism.
DRG reduction. It is seen from Fig. 15a that the re- The reduction time for large mechanisms was
duction time varies fairly linearly with the number of shown theoretically and experimentally to be lin-
reactions in the mechanism, as anticipated from the early proportional to the number of reactions in the
 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436 35

Fig. 14. Calculated relative error in temperature rise for 0.1% iso-octane, as function of the user-specified error threshold value ,
in a flow reactor with = 1.0, p = 1 atm, initial temperature 1080 K, and residence time = 0.2 s, for both local and global
reduction.

efficiency and fully automated reduction render DRG

particularly suitable for case-specific or real time re-
duction, in which the reduction time is required to be
much shorter than or at least comparable to the simu-
lation time.
The extent of reduction can be further increased,
albeit to a much smaller degree, through a second-
stage application of DRG and the subsequent applica-
tion of the CSP importance index.
The usefulness of the present reduction algorithm
and strategy clearly increases rapidly with increas-
ing mechanism size. Indeed, this is probably the only
approach that can expeditiously and automatically
generate skeletal mechanisms from extremely large
detailed mechanisms. Specifically, a 188-species and
a 233-species skeletal mechanism, respectively, were
generated from a detailed 561-species n-heptane
mechanism and a detailed 857-species mechanism.
Not only do these skeletal mechanisms mimic well
the responses of autoignition and PSR over extensive
parametric ranges, from which they were developed,
but the agreement was almost equally satisfactory for
Fig. 15. Measured average reduction time of DRG for a the independent system of JSR.
single reaction state as a linear function of (a) number of The relative errors for local and global DRG re-
reactions and (b) number of species of the detailed mecha- ductions were found to be bounded by a linear func-
nism. tion of the accuracy threshold. Furthermore, while the
error for the global reduction is much smaller than
mechanism. Furthermore, the DRG method exploits a that of the local reduction, the skeletal mechanism de-
normalized error threshold value which is very easy rived with global reduction is much larger. As such,
to specify. Consequently the reduction requires only the threshold value employed in mechanism reduction
minimal system-dependent knowledge and the entire should be selected conservatively for local reduction
reduction process can be fully automated. The high- but aggressively for global reduction.
36 T. Lu, C.K. Law / Combustion and Flame 144 (2006) 2436

We close this exposition by emphasizing the high [15] T.F. Lu, Y. Ju, C.K. Law, Combust. Flame 126 (12)
efficiency of the present linear time algorithm, which (2001) 14451455.
is responsible for the extremely short time in mech- [16] M. Valorani, H.N. Najm, D.A. Goussis, Combust.
anism generation, as for the case of 50 ms for iso- Flame 134 (12) (2003) 3553.
[17] H. Rabitz, M. Kramer, D. Dacol, Annu. Rev. Phys.
octane.
Chem. 1134 (1983) 419461.
[18] T. Turanyi, J. Math. Chem. 5 (1990) 203248.
[19] A.S. Tomlin, T. Turanyi, M.J. Pilling, Comprehensive
Acknowledgment Chemical Kinetics, Elsevier, 1997, pp. 293437.
[20] S. Vajada, P. Valko, T. Turanyi, Int. J. Chem. Kinet. 17
This work was supported by the Air Force Office (1985) 5581.
of Scientific Research under the technical monitoring [21] B. Bhattacharjee, D.A. Schwer, P.I. Barton, W.H.
of Dr. Julian M. Tishkoff. Green, Combust. Flame 135 (2003) 191208.
[22] H. Wang, M. Frenklach, Combust. Flame 87 (34)
(1991) 365370.
[23] T. Turanyi, New J. Chem. 14 (1990) 795803.
References [24] T.F. Lu, C.K. Law, Proc. Combust. Inst. 30 (2005)
13331341.
[1] N. Peters, F.A. Williams, Combust. Flame 68 (2) (1987) [25] H.J. Curran, P. Gaffuri, W.J. Pitz, C.K. Westbrook,
185207. Combust. Flame 114 (12) (1998) 149177.
[2] K. Seshadri, N. Peters, Combust. Flame 73 (1) (1988) [26] H.J. Curran, P. Gaffuri, W.J. Pitz, C.K. Westbrook,
2344. Combust. Flame 129 (3) (2002) 253280.
[3] T. Turanyi, A.S. Tomlin, M.J. Pilling, J. Phys. Chem. 97 [27] http://www-cms.llnl.gov/combustion/combustion2.
(1993) 163172. html#n-C7H16, August 2004.
[4] J.Y. Chen, Combust. Sci. Technol. 57 (1988) 8994. [28] J.D. Blouch, C.K. Law, Proc. Combust. Inst. 28 (2000)
[5] C.J. Sung, C.K. Law, J.Y. Chen, Proc. Combust. 16791686.
Inst. 27 (1998) 295304. [29] A. Chakir, M. Belliman, J.C. Boettner, M. Cathonnet,
[6] C.J. Sung, C.K. Law, J.Y. Chen, Combust. Flame 125 Int. J. Chem. Kinet. 24 (1992) 385410.
(2001) 906919. [30] R.P. Lindstedt, L.Q. Maurice, Combust. Sci. Tech-
[7] T. Turanyi, J. Toth, Acta Chim. Hung. 129 (1992) 903 nol. 107 (1995) 317353.
914. [31] T.J. Held, A.J. Marchese, F.L. Dryer, Combust. Sci.
[8] H.S. Soyhan, F. Mauss, C. Sorusbay, Combust. Sci. Technol. 123 (1997) 107146.
Technol. 174 (1112) (2002) 7391. [32] N. Peters, G. Paczko, R. Seiser, K. Seshadri, Combust.
[9] U. Maas, S.B. Pope, Combust. Flame 88 (34) (1992) Flame 128 (12) (2002) 3859.
239264. [33] M. Bollig, H. Pitsch, J. Hewson, K. Seshadri, Proc.
[10] S.B. Pope, Combust. Theory Modeling 1 (1997) 4163. Combust. Inst. 26 (1996) 729737.
[11] B. Yang, S.B. Pope, Combust. Flame 112 (1998) 85 [34] R. Seiser, H. Pitsch, K. Seshadri, W.J. Pitz, H.J. Curran,
112. Proc. Combust. Inst. 28 (2000) 20292037.
[12] S.H. Lam, D.A. Goussis, Proc. Combust. Inst. 22 [35] S.L. Liu, J.C. Hewson, J.H. Chen, H. Pitsch, Combust.
(1988) 931941. Flame 137 (3) (2004) 320339.
[13] S.H. Lam, D.A. Goussis, Int. J. Chem. Kinet. 26 (1994) [36] http://www-cms.llnl.gov/combustion/combustion2.
461486. html#i-C8H18, August 2004.
[14] A. Massias, D. Diamantis, E. Mastorakos, D.A. Gous- [37] H. Pitsch, N. Peters, Proc. Combust. Inst. 26 (1996)
sis, Combust. Flame 117 (4) (1999) 685708. 763771.