An Investigation of the Accuracy of Manifold Methods and

Splitting Schemes in the Computational Implementation of

Combustion Chemistry
B. YANG* and S. B. P O P E
Sibley School of Mechanical Aerospace Engineering, Cornell University, Ithaca, New York 14853-7501

This paper is concerned with the efficient computational implementation of combustion chemistry, for use in
PDF methods and other applications. A new method of coupling reactions and mixing processes based on
manifold points with detailed chemistry is developed. Investigations are made of three different kinds of
methods: the direct numerical integration of the coupled reaction and mixingequations: the direct numerical
integration of the equations obtained by using operator-splitting to split reaction and mixing; and the new
method--solving the split system based on manifold points with detailed chemical kinetics. Errors between
the solutions are studied. It is found that chemical reactions have a significant influence on the accuracy of
operator-splitting methods. The solution of the split systems based on manifold points provides an accurate
representation for the solution of the full coupled equations. This means that tabulations can be made on
manifolds with no simplification made to the chemistry. © 1998 by The Combustion Institute

1. INTRODUCTION plifies turbulent combustion calculations, on

the other hand, it introduqes the infamous
In the last decade, great progress has been closure problems, especially the closure prob-
made in combustion research, especially in the lem with chemical reaction terms. Since in
computation of laminar flames [1, 2, 3, 4], and P D F calculations of turbulent combustion the
in the probability density function (PDF) averages of the chemical reaction terms can be
method for turbulent combustion [5, 6, 7]. For calculated, P D F methods overcome the closure
one-dimensional laminar flames, by consider- problerr~ with the reaction terms. It has been
ing the transport mechanism, the detailed shown that the P D F method is the most
chemical kinetic mechanism and the interac- promising metbod to calculate turbulent com-
tions between these two basic processes, today bustion [6]. P D F methods have been success-
it is a routine matter to calculate flame veloci- fully employed to calculate laboratory turbu-
ties, extinction, ignition, temperature and lent flames: tt~Ley can predict phenomena such
species distributions, from the governing equa- as super equilibrium radical levels, and local
tions. Results are in good agreement with those extinction [7]. Because of these advantages,
obtained from experiments [8, 9]. However, for P D F methods are becoming used increasingly
turbulent combustion, because of the complex- in industrial eombustor codes.
ities of turbulent flow, chemical reactions, and Although PDF methods have shown great
the interaction between them, in the foresee- promise in studies of turbulent combustion,
able future it is impossible to calculate the there is still a challenge to be o v e r c o m e - - c o u -
combustion flow field by directly integrating pling the detailed description of the turbulent
the basic governing equations. So averaging combustion flow field provided by PDF meth-
and modeling are necessary in turbulent com- ods with detailed chemical kinetic mechanisms.
bustion studies. Averaging, on one hand, sim- The problem can be stated thus: in a P D F
calculation of turbulent combustion, let ~b(t)
represent the composition of a particle at time
*Corresponding author. Address: 138 Upson Hall, Box
XX, Cornell University, Ithaca, New York 14853-7501. t, then what is the increment in composition
Telephone: 607-255-8270, Fax: 607-255-1222, e-mail A~b(t) due to reaction over a time step St? In
byang@mae.cornell.edu principle, this question can be answered by

COMBUSTION AND FLAME 112:16-32 (1998)

0010-2180/98/$19.00 © 1998by The Combustion Institute
pll soo 10-2180(97)00013-8 Published by Elsevier Science Inc.
COMPUTATIONS OF COMBUSTION CHEMISTRY 17

directly integrating the ordinary differential served, there is a wide range of time scales for
equations stemming from the detailed kinetic chemical reactions, from 10 -9 second to sec-
mechanism. But in practice, since a typical onds. Fast reactions, or reactions with small
combustion system involves dozens of chemical time scales, quickly bring composition points
species and hundreds of chemical reactions, down to attracting manifolds in the composi-
and it is needed to do such integrations on the tion space. Then composition points move
order of 10 9 times, the direct numerical inte- along on manifolds. By assuming that the
gration of the equations would require a huge movement of the composition point away from
amount of supercomputer time (several hun- manifolds to be zero, detailed chemistry can be
dred days) and thus make it impossible in simplified. The manifold method overcomes the
practical use. So simplifications of detailed ki- drawbacks of the reduced mechanism method.
netic mechanisms have been made in the past It requires no preliminary knowledge of which
in order to reduce the demand on computer chemical species are in steady-state and which
time. In practice there is a preprocessing stage chemical reactions are in partial equilibrium.
in which results from the calculations of sim- The only given assumption is the dimension of
plified chemistry are tabulated as functions of the manifold.
a few variables. Then these tables are used in The manifold method has been successfully
turbulent combustion calculations. used in both laminar flames and turbulent
There are basically two different ways of combustion studies [7, 15, 16]. In these studies,
doing the simplification of detailed chemistry: a manifold with fixed dimension, for example,
the reduced mechanism method [10, 11, 12] two dimensions, has been considered. The re-
and the intrinsic low-dimensional manifold sults from manifold calculations are tabulated
(ILDM) method [13]. For the reduced mecha- in a preprocessing stage. Then the method of
nism method, the simplification made to the tahle-!ook-up is used in PDF calculations.
detailed chemistry is achieved by introducing There are still some difficulties and inconve-
steady-state assumptions for some species, usu- niences with this approach to implementing
ally the intermediate species, and the partial manifold methods:
equilibrium assumptions for particular reac-
tions. The reduced mechanism method has • in general, it is not straightforward to
been employed in laminar flame calculations parametrize the manifold,
and in turbulent combustion calculations [11, • in different regions of the composition space,
14]. It has several disadvantages because of its manifolds of different dimension are appro-
fundamental philosophy. For the reduced priate,
mechanism method, one needs to know in ad- • the table generation (which is not fully auto-
vance which species are in steady-state, and mated) must be performed for each set of
which reactions are in partial equilibrium. Re- conditions of interest (fuel, pressure, equiva-
duced mechanism systems are derived manu- lence ratio, etc.),
ally from the given detailed chemistry. For • the whole of the manifold is (wastefully)
different fuel/oxidizer systems, or even for the tabulated since it is not known a priori which
same fuel/oxidizer system under different con- regions are needed.
ditions, different reduced mechanisms should
So wc have been motivated to investigate new
be used. Thus it requires a considerable amount
implementations of manifold methods which
of human time and labor to develop and test
can overcome these difficulties and inconve-
such systems. Assumptions of partial-equi- .
nlPnr,~,
librium and steady-state used in the reduced ........ s.
mechanism method are only valid in particular
reaction ranges. Also the accuracy cannot be 2. OUTLINE OF THE PAPER
given and controlled.
The manifold method is based on a more This paper addresses several issues related to
intrinsic study of the chemical reaction process the accurate numerical implementation of
happening in combustion [13]. As it is ob- combustion chemistry.
18 B. Y A N G A N D S. B. POPE

2.1. Splitting Strategies sought based on the particle's composition ~.

The table is searched for an entry 4, ~" suffi-
Most numerical implementations of table-look- ciently close to 4,. If one exists, then A~b"~ is
up chemistry correspond to the zero-order used to approximate A~h. If a sufficiently close
splitting strategy. That is, the evolution equa- table entry does not exist, then A~p is com-
tions for the fluid composition with the omis- puted by the direct integration of equations
sion of the reaction term are advanced for a from detailed chemistry, and the result is added
small time step 6t: and then the increment in to the table. This tabulation scheme is not
composition over the time 6t is added to the further described or used here. Rather we ad-
result. In PDF methods, the first part corre- dress the question: how can information about
sponds to mixing, and so this strategy is re- the manifold be used to determine the evolu-
ferred to as mix-then-react. tion of composition ~b(t) due to mixing and
Important questions addressed here are: reaction? The result, described in Section 5.3,
How large are the numerical errors introduced is an exact solution to the equations linearized
by this splitting? And, do higher-order splitting about a manifold point.
strategies offer advantages? Splitting strategies The accuracy of this solution compared to
are described in Section 4, and numerical re- the detailed kinetics is investigated in Section
sults are given in Section 6. 6.

2.2. Variable-Dimension Manifold 2.4. Pairwise Mixing Stirred Reactor (PMSR)

In previous implementations of the intrinsic A simple test case that has been used previ-
low dimensional manifold, a manifold of fixed ously is the partially stirred reactor (PaSR)
dimension is considered. It is preferable to [17]. For a PaSR, the system can either be
define the manifold instead in terms of a time premixed or diffusion. The mixing model used
scale r*, which is specified to be somewhat is the interaction-by-exchange-with-the-mean
smaller than the smallest fluid-mechanical time (IEM) model [17]. In the statistically stationary
scale in the problem (usually the mixing time state of a PaSR calculation, the composition
scale Zm~x). The definition of this manifold is ~b(t) of a particle is a unique function of the
given in Section 5.1. particle's age (i.e., time since it entered the
In order to implement the manifold method reactor). As a consequence, all particle compo-
it is necessary to solve the "closest manifold sitions lie on a one-dimensional manifold since
point" problem. That is, given ~'* and a compo- chemical reactions always pull particle compo-
sition ~b°', determine the closest point to ~b~°~ sition toward the equilibrium state which is a
that is on the manifold. The solution to this zero-dimensional manifold; and so, in fact, the
problem is given in Section 5.2. PaSR is not a strenuous test of simplified
chemistry schemes.
To provide a better test case we have devel-
2.3. Solution of the Mixing and Reaction oped a different implementation of the PaSR,
Equation referred to as the "pairwise mixing stirred re-
actor" (PMSR). This is described in Section 3
In the overall method currently being devel-
and used for all of the tests reported in Section
oped, properties of the manifold are stored
6.
using "in situ adaptive tabulation" (ISAT). As
the PDF 'calculation is performed, an unstruc-
tured table is generated, containing N pairs of 3. THE PAIRWISE MIXING STIRRED
compositions and their corresponding incre- REACTOR
ments, {~"~, A ~bt"), n = 1, 2 . . . . . N}. The table
is stored in a structure that is initially empty The PMSR is characterized by a residence
( N = 0). For each particle on each time step in time rrc~, a pairing time ~'pa~r,and a mixing time
the PDF calculation, the increment A~b is rmix. In numerical implementations the reactor
COMPUTATIONS OF COMBUSTION CHEMISTRY 19

consists of M particles, which arc arranged as Operator splitting techniques, also known as
L = M / 2 pairs of particles. the method of fractional steps, were intro-
At the initial time t = 0, the particle compo- duced in fifties [18, 19]. The mathematical
sitions 4/"~(0) are initialized to the chemical fe,undations are summarized in [20]. Operator
equilbrium state. splitting methods offer a wide variety for split-
With At being a time step, "reactor events" ting a system. Here we investigated ,several
take place at the discrete times, At, 2At . . . . . kinds of operator splitting methods: the zero-
These events consist of outflow, inflow, and order splitting method, the first-order splitting
pairing. At each event, mou t = LAt/rrc ~ parti- method, predictor-eorrector methods, and
cAe pairs (selected at random) flow out of the Strang's sequential splitting method. Splitting
reactor (i.e., they are discarded). A "'candidate errors depend on fractional time steps. To
pile" is created consisting of 2 r a i n particles investigate the errors, we divide the time step
with the specified inlet composition (where mi, in the PMSR, At, into several subtime steps,
= m,ut) , and mpair = LAt/rpair pairs ran- 6t, 6t = At/nsu b, where n~ub is the number of
domly selected from the reactor. Then the subtime steps.
particles in the candidate pile arc paired ran- For the zero-order splitting method, Eqs. l
domly, and are then added to thc reactor. and 2 are split into a pure mixing system:
After this event, the reactor again consists of
L pairs of particles. dcb"~/dt = (tb ~/' - tb"')/rmi ~, (3)
Between reactor events, the particle compo-
sitions evolve continuously. With i and j being d~b'J~/tlt = - (tb 'j' - ~b"0/'/'mix, (4)
the indices of a pair of particles, their composi-
tions evolve according to the equations and a pure chemical reaction system:

ddpO~/dt = S(~b ti~) + (4) ~i~ - ~(i))/7"mi x, (1) dd)"~/dt = St&m), (5)

dd)~J~/dt = S t + ~ ) + (4' "~ - 4;')/r,,,~,, (2) dd)tJ~/dt = S(~btJ~). (6)

where S is the chemical reaction term, and rmix In the first fractional step, the pure mixing
is the mixing time scale. system is integrated over a time step St, to get
Notice that the mixing--i.e., the second ~b~nli~x and -~m~x,d~' J)"then, in the second fractional
terms on the right-hand sides of Eqs. 1 and 2 step, the pure reaction system is integrated
- - i s between pairs of particles, rather than (from initial conditions t~mixti~ and "rmix-th(J)
) over a
with the mean (as in the PaSR). Hence we time step 6t to get 4,"~(t + 6t) and thm(t +
refer to this test case as a pairwise mixing 6t). The traditional table-look-up method is
stirred reactor (PMSR). Compared to the basically this zero-order splitting method.
PaSR, the particles access a greater region of For the first-order splitting method, the first
composition space. fractional time step is the same as that in the
In the calculations reported in Section 6, the zero-order splitting method. The values of t~mix "~
PMSR parameters are specified as M = 100, -~m~x
a n d tfh{ J) are obtained after the first fractional
rr~~ = 10 -2 (s), yp~i~ = l0 3 (s), 7mix = 10 3 (s), step. In the second fractional step, the system
A t = 6.0 × 10 -4 (s). is split as

dd~"~/dt = S(~b")) + F, (7)

4. OPERATOR SPLITFING SOLUTIONS
F O R THE FULL COUPLING EQUATIONS d~btJ)/dt = S(~btj~) - F, (8)
In order to use the table-look-up method, it is
where F is the constant mixing vector defined
necessary to split Eqs. I and 2 so that for a
as:
given composition point, the increment of it
only depends on the given value of this compo-
sition point and the information about mixing. r = -[&'J",mix- d?~x)/6t'. (9)
20 B. Y A N G A N D S. B. POPE

This system is integrated over a time step ~t to For a chemically reacting system, the time
get the solutions ~b")(t + t~t) and th(/)(t + ~t). scales of reaction processes cover a wide
Two different predictor-corrector methods range. For a given point in composition space,
were investigated. For both methods, the first because of the existence of fast reaction
fractional step is the same as that in the first- components, the point moves quickly to a low-
order splitting method. In the second frac- dimensional manifold; most of the slow reac-
tional step, for the first predictor-corrector tion processes take place in the low-dimen-
method, the split systems Eqs. 7 and 8 are sional manifold. The low-dimensional manifold
integrated over a time step 8t to get the pre- is identified by studying the Jacobian matrix of
dict values: ,hi(,i) and ebb,J); then the mixing the chemical reaction source term, S. This is
vector F is calculated based on [qbti)(t)+ shown in the following.
¢b),i)]/2 and [~b(J~(t) + ~b~,J)]/2; in the corrector Let the composition of a particle be ~b = {Yt,
step, the split systems are solved over a time Y2 . . . . . Y,, h} T, where Y/, i = 1, 2 . . . . . n~, is the
step 8t to get the finaJ solutions ~p(i)(t + St) mass fraction of species i, n~ is the number of
and ~bt/)(t + ~St). species, and h is the specific enthalpy. Let
For the second predictor-corrector method in n = n s + l be the dimension of ~b. For a ho-
the second fractional step, the split systems mogeneous, adiabatic, isobaric system, the evo-
Eqs. 7, 8 are integrated over half the time step lution of th is determined by the equation:
~ t / 2 to get the predict values: ~b~,i~ and d~(J).
-rp ,
then the mixing vector F is calculated based on d4,
-- = S, (10)
,1ti~ and d~(J'
-rp In the corrector step, the split
. dt
systems are integrated over a full time step ,St
to get the solutions ~b(°(t + ~;t) and ~b(J)(t + where S is the source term due to chemical
~t). reactions (with enthalpy conservation giving S n
Strang's sequential splitting method [21] is = 0). For a given composition point, ~, the
performed as follows: advance the pure mixing Jacobian matrix J is
system by half the time step, 8t/2, to get
)"' I/2) and ~b~i~/2>;"then starting from &ese tgSi
conditions, integrate the pure reaction system
over a full time step, ,St, to get "wrcd,"~c and -t~'rcc,
,~(J)"
finally, from these conditions, integrate the The Jacobian matrix has n eigenvalues, Ai,
pure mixing system by half the time step, ~t/2, i = 1, 2 . . . . . n, and the corresponding n eigen-
to get 4~")(t + ~t) and 4~(J'(t + ~t). This vectors. These eigenvalues identify n timescales
method can also be thought of as mix-react-mix. of the movement of a particle in composition
The numerical errors incurred in these dif- space. The corresponding eigenvectors define
ferent methods have been quantified, and are the directions associated with the eigenvalues
presented in Section 6. [13]. A manifold point in the composition space
is the point at whieh the movements in direc-
tions corresponding to small timeseales vanish.
5. S O L U T I O N S OF T H E SPLIT SYSTEMS It is common sense that for a chemical reac-
BASED ON MANIFOLD P O I N T S tion system, in some regions of the composi-
5.1. M a n i f o l d P o i n t s tion space, the movement goes very fast, as
those near equilibrium; while in other regions
In order to restrict the points to be tabulated of the composition space, reactions progress
in a relatively small region in the composition slowly. If we use a fixed dimension for the
space, we represent our solutions by informa- manifold in the whole composition space, large
tion of manifold points. The manifold method errors can be expected in slow reaction re-
is described in detail in [13]. Manifold methods gions. So it is preferable to change the dimen-
have been used in the calculations of laminar sion of a manifold according to the reaction
flames [16, 15] and turbulent flames [7]. rate. This can be done by introducing a time
COMPUTATIONS OF COMBUSTION CHEMISTRY 21

scale ~-*. The time scale ¢* should be smaller in the order of descending real parts, then the
than the smallest fluid-mechanical time scale subspace spanned by the first n m Schur vec-
in the problem, here the mixing time scale ~'mix" tors, q],q~ . . . . . q~,, is the slow invariant sub-
At every point in the composition space a space. If in the matrix D, the eigenvalues are
f a s t s u b s p a c e and a s l o w s u b s p a c e are defined organized in a different way such that t h e n I
based on 7* and the eigenvalues and eigenvec- eigenvalues with the smallest real parts appear
tots of the Jacobian J at that point. If the first in D, then the corresponding subspace
eigenvalue Ai has real part greater than - l / r * spanned by the first n I Schur vectors q{,
(i.e., Re(A i) > - l / r * ) , then the correspond- q( . . . . .. q ln !, is the fast invariant subspace [22].
ing eigenvector is in the slow subspace. Con- Taking the reaeuon-rate vector S as an ex-
versely, if Re(Ai)< - l / r * , then the corre- ample, any n-vector can be expressed as the
sponding eigenvector is in the fast subspace. sum of slow and fast components
We denote by nm the number of eigenvalues
with Re(Ai)> - l / r * ) : so that nm is the di- S = S ~ + S?. (16)
mensionality of the slow subspace. The dimen-
sionality of the fast subspace is n t = n - n , , . Here, S t and S I are in the slow and fast
Put another way, the fast and slow subspaces subspaces respectively. Consequently, S r can
are the invariant subspaces of J corresponding be written
to eigenvalues with real part less than, and,
respectively, greater than, - l / r * . S f = V/S:, (17)
The eigenvectors can be taken as basis vec- where S: is an nr-vector giving the compo-
tors in the two subspaces, but computations nents of S / in the VLbasis. These components
are more stable if orthonormai bases are em- are determined as follows. Let the n × n ma-
ployed. Let the orthonormal basis for the slow trix V be
subspace be:
V - [WW], (18)
V'= [q~,q~ . . . . . q~,~], (12)
and let its inverse be
and let
I I
vf = [ql'q' ..... q~t]' (13) V -~ = W = Wf , (19)..

be the orthonormal basis for the fast subspace. where W ~ is an n,,, × n matrix, and W / is an
V s and V / can be readily computed by the n: × n matrix. Then the components S/ are
Schur decomposition method [22]. It is briefly given by
described below.
The Jacobian matrix J can be decomposed g: = WfS. (20)
by the Schur decomposition into the following
form: The manifold is defined as the points in the
composition space where the fast reactions are
Q r j Q = D + N, (14) equilibrated. In terms of equations, a point on
the manifold satisfies S / = 0, or equivalently
where, Q is a n × n unitary matrix, D = S / = 0, or
diag(A I, A2,..., A,), and N is a n x n strictly
upper triangular matrix. If Q is partitioned W/S = 0. (21)
like:
To summarize, at every point in the compo-
Q = [ql,q2 . . . . . q . ] , (15) sition space, in terms of the local Jacobian J
and the time scale r* we define the slow and
then qi are referred to as Schur vectors [22]. If fast subspaces and the corresponding matrices
in the matrix D, the eigenvalues are organized W, V r, W s, and W f. A point is part of the
22 B. YANG AND S. B. POPE

manifold if, and only if, Eq. 21 is satisfied at A~b"~, Eqs. 22 and 23 become
that point.
W"A,b"1 = 0, (26)
5.2. The Closest Manifold Point [W/J]A~b"' = W / { J ( ~ b " ' - 4,~''') - S(d,"D}.
The composition d' of a particle in a PDF or (27)
PMSR calculation corresponds to a point in
composition space. This is usually close to but Equations 26 and 27 are underdetermined.
not exactly on the manifold with the dimension Generally, they are solved by the singular value
determined by the given time scale z*. As decomposition method to get a minimum norm
mentioned in Section 2, we want to describe solution [22]. If this solution has negative mass
the evolution of &"~)(t) in terms of quantities fractions (in violation of conditions Eq. 24)
on the manifold (which will be tabulated). then instead we use the quadratic program-
Hence we need to find the closest manifold ming method [23, 24] to get the minimum
point. The closest manifold point &tin) of the norm solution for Eqs. 26 and 27 under con-
given point &(0) is a point which satisfies the straints (24).
following conditions:

(a) &~,n~ has the same enthalpy and element 5.3. Coupling Mixing and Reaction in
compositions as ~b~°). Manifolds
(b) It is on the manifold. This section describes how the coupling of
(c) It is realizable, which means that mass mixing and reaction are treated to get the
fractions of species are non-negative. value of ¢k'°'(6t) for a given particle, q~,0),
(d) It is as close as possible to ~0~ while based on its closest manifold point, ~btin).
satisfying conditions (a) to (c). The sy~,em obtained from operator splitting
Let W s and W / be the matrices defined by methods can be written in the general form:
Eq. 19 and let W e be the matrix determining
the element mass and enthalpy conservation, dd~ S + F. (28)
then the closest manifold point to tht°~ is de- dt
termined as the solution q~t'~ to the following
problem. For the zero-order splitting system, and the
Minimize the 2-norm of ~bt ' ' ' - th°'', subject system from Strang's algorithm, F is zero; for
to the conditions: the first order split system, F is the mixing
vector defined by Eq. 9.
W"(4, t''' - 4, ~°)) = 0, (22)
With 8 ~ - ~ b - . ~ " , the linearization of
Eq. 28 about ~ t " is
W / S ( d , t ' ) ) = 0, (23)
dS~b
~b~ > 0, i = 1,2 . . . . . n~. (24) T = S ( ~ " ' ) + J6O + F. (29)

These equations are solved by an iterative The exact solution of this equation is:
method. Let Oto be the estimate of th''') after
the i-th iteration: the next iteration is 6~b(6t) = A6d~(0) + BS(~b~m)) + BF, (30)

~b(,n)= ~(i+I)~ ~(i)+ ~f~(i)~_~b(o)+ A~(i) where the matrices A and B are defined by
(25)
A = e TM, (31)

With W evaluated at ~b(i), linearize Eq. 23

around 4,(i), write the equations in terms of
B = [~' e J~st-
"0
,) dT. (32)
C O M P U T A T I O N S OF C O M B U S T I O N C H E M I S T R Y 23

Since chemical reactions take place in a fold Point_s." For the divergent subspacc ~ ,
wide range on timescale, some processes the corresponding eigenvalues are R e ( h ) > _
quickly, it is preferable to represent the second l / r 't. i = 1,2 . . . . . n a, where na is the dimen-
term in the above equation by an exact solu- sion of the divergent subspace, r a = et6t is a
:ion: time scale of the divergent subspace, a is a
small positive value, we chose it to be 0.1, the
corresponding matrix is V a. For the the slow
61~ R = --10~t S tit, (33)
subspacc ,'/'. the corresponding eigenvalues,
-I/r* < R e ( A / ) < l / r a, i = !,2 . . . . . n,i ....
integrating from the starting point d,~"L Thus here n,l,,,, is the dimension of the slow sub-
Eq. 30 is modified to: space, the corresponding matrix is V ~. For the
fast subspacc :7, the corresponding eigenval-
34,(3t) = A6~b(0) + 6d) n + BF. (34) ues, Re(A,) _< - l / r * , i = 1,2 . . . . . n.t, here nf
is the dimcnsion of the fast subspace, the cor-
Due to the wide range of eigenvalues of the
rcsponding matrix is V f. Let
Jacobian matrix, difficulties may arise in the
calculation of matrices A and B. The numerical V = [VaV'Vf]. (37)
technique used in this paper is the method
proposed in [25]. Here is a brief description of The inverse of the matrix V is W , W = V i
the basic idea of the method. According to which can bc written as
[25], if

w = . (38)
tw'l
where 1 is an identity matrix of the same order For ;,lily vector X,
of the matrix J, then, lor 6t >_ Ik

Wx = . (39)
[¢,'J
where A and B are the matrices defined by where .~a .~, and ~,f, are the components of x
Eqs. 31 and 32. So if the matrix e c~' is calcu- in the bases of the divergent, slow, and fast
lated, then A and B are known. The Padd subspaces. Thus premultiplying the iinearized
approximation method [22] is used to calculate
coupling equation Eq. 29 by W, we get
e c6t.
The procedure mentioned above can be used
to calculate the matrices A and B when all the
eigenvalues of the Jacob(an matrix J are not
large compared with the inverse of the time
step, 1 / 6 t . If there exist positive eigcnvalues
that are large compared to l / 6 t , then the
system is not stable over the time step. This
may correspond to the region near ignition
+
[:::l , (40)

where reaction rates increase rapidly. In this t~Jj

situation, a special treatment is needed to solve where
the linearized system Eq. 29. The composition
space is divided into three subspaces: the di- [D o] = w Jr,
vergent slow subspace 3 , the regular slow [o E
subspace ,Sa and the fast subspace ~ They are
determined by Sehur decomposition in the D is a n a x n a matrix, E is a (nq,,,,. + nf) ×
same manner mentioned in the section "Man(- (n~t,,,v + nf) matrix.
24 B. Y A N G A N D S. B. P O P E

The equation for the components in the 6. RESULTS AND DISCUSSIONS

divergent subspace is
The chemical kinetic mechanism used in the
d8 ~,l calculations is shown in Table 1 which was
- - = ga(4~ ..... ) + D,~& d + F " . (41) used by S. M. Correa at General Electric Re-
dt
search Center. This is a mechanism for
The Eulcr explicit scheme is used to get m e t h a n e / a i r combustion, with sixteen species
and forty-one reactior, s. In our PMSR calcula-
ag'~(a,) tions, the incoming particles to the PMSR are
stoichiometric m e t h a n e / a i r at 3(}0 K. The
= a~,'(0) + a,[~,'(~,,,',) + nag,,' + ~,'1 pressure of the PMSR is one atmosphere. The
(42) particles in the PMSR are initialized to the
equilibrium condition. The mixing time scale is
For the slow and fast components, the solu- Tmix = 10 3 s, the pairing time scale is "/'pair :
tion is 10 3 s, and the residence time scale is rr¢~ =
10 2 s. The time scale used to determine the
dimension of manifolds is r* = lt)- 4 s. There
=,- +f,,e F, . . . . , a , are 100 particles in the reactor, the time step is
At = 6 x 1(1-4 s.
From direct integrations of Eqs. I and 2, the,
changes of the average mass fractions of all the
particles in the PMSR with time are shown in
Figs. 1 and 2. Figure 1 shows changes of the
Replacing the S(~bI''')) terms in Eqs. 42 and 43 average mass fractions of major stable species.
by the exact solution Figure 2 shows changes of the average mass
fractions of major radicals. It can be seen that
3g, r = W&b r, (44) after about ().6 residence times, the average
mass fractions reach a statistically steady state.
then the increment of the composition vector There is some fuel, CH 4, and oxidizer, O 2 in
in the bases of V is the reactor, while initially none of these reac-
tants exists since the reactor is initialized to
the equilibrium conditon. So the statistically
steady-state condition of the reactor is signifi-
cantly different from the equilibrium condi-
e~:~'ag,g0) ]
tion.
Figure 3 shows (kb m. - cbt, ol) which is the
2-norm of vector 4~t,<, - d,,.,
+ .[,~' e ~:~' ~ ) d r F ".~ .
J,? e ~:''~' "' dT~'J
(45)

[(/',,(- 4'..I =
~i=1
('b.(.,- 4',,,,.,)-"
)- •

So the increment of the composition vector is

averaged over the particles in the PMSR as a
function of time. Here the subscript, D means
8qb(St) = v/ 6g,-,/ (46) direct integrations, C means the original full
coupling equations (Eqs. 1 and 2), and 0 means
La '] zero-order splitting system. The average value,
(I4't~c - tbt>0l) is thus defined by
In summary, Eq. 46 provides a stable and
accurate solution to the coupled mixing and 1 M
rcactic.n equation linearized about !he mani- (14,,,c - 4,o,,I) = ~ i=El
: "*~i~
" ( - ,/,g~,l,
fold point &(").
COMPUTATIONS OF COMBUSTION CHEMISTRY 25

TABLE I

T h e Chemical Kinetic Mechanism ~

N o. React ion An n En

I tl + O , = Ot1 + O 1.5t)E + 17 -0.927 16874.

2 O + i l , = O I t + ti 3.87E -~ 04 2.70 6262.
3 0 t t + H , = t t , 0 + tt 2.16E + 08 1.51 3430.
4 Oti + Ott = O + II,O 2.1lIE + 08 1.40 -397.
5 tl + tl + M = H , + - l i d 6.40E + 17 -1.0 0.
6 II + O i l + M = H , O + M 8.40E + 21 -2.04;) 0.
7 t l 4- O , + M = ttO. + M 7 . 0 0 E + 17 -0.80 O.
8 t t O , + II = O t l + O I I 1.50E + 14 0.0 1004.
9 I I O , + tl = t t , + O , 2.50E + 13 0.0 693.
10 I I O , + O = O , 4- O11 2.llOE + 13 0.0 0.
II IIO, + Ott = ti,O + O, 6.02E + 13 0.0 0.
t2 11202 + M - 0 i i -r 0 1 t + M i.l)tiE + i7 0.11 45411.
13 C O + O t t = C O , + II 1.51E + 07 1.3 -758.
14 C O + O + M = CO, + M 3.01E + 14 0.0 3011.
15 t l C O + tt = H , + C O 7.23E + 13 0.0 O.
16 l t C O + O = OI1 + C O 3.00E + 13 0.0 0.
17 HCO + OH = ti,O + CO 1.00E + 14 0.0 0.
18 tlCO+ O, = tlO, + CO 4.20E + 12 0.0 0.
19 HCO + M = H + CO + M 1.86E + 17 - 1.0 16993.
211 C t i , O + H = H C O + H., 1.26E + 08 1.62 2175.
22 CH20 + O = HCO +'OH 3.50E + 13 0.0 3513.
23 CH20 + OH = HCO + H,O 7.23E + 05 2.46 -970.
24 CH20 + 0 2 = HCO + [IO, I.(X)E + 14 0.0 39914.
25 C H , O + C H 3 = H C O + CI-I a 8.91E - 13 7.40 -956.
26 CH,O + M = HCO + H + M 5.00E + 16 0.0 76482.
27 CH 3 + O = CH20 + H 8.43E + 13 0.0 0.
28 C H 3 + O H = C H , O + 1-1, 8.(XiE + 12 0.0 0.
29 CH 3 + 02 = C H 3 0 + O 4.30E + 13 0.0 30808.
30 CH 3 + 02 = CH20 + OH 5.20E + 13 0.0 34895.
31 CH 3 + HO 2 = CH30 + OH 2.28E + 13 0.0 0.
32 C H 3 + H C O = CH4 + C O 3.20E + II 0.50 0.
33 CH 4 + H = CH 3 + H 2 7.80E + t)6 2.11 7744.
34 CH 4 + O = CH 3 + O H 1.90E + 09 1.44 8676,
35 C H 4 + 0 2 = C H 3 + i-IO 2 5.60E + 12 0.0 55999.
36 CH a + OH = CH 3 + H20 1.50E + 06 2.13 2438.
37 CH 4 + HO 2 = CH 3 + n20 2 4.60E + 12 0.0 17997.
38 CH30 + H = CH20 + H, 2.00E + 13 0.0 0.
39 C H 3 0 + O11 = C H . , O + H , O 5.00E + 12 0.0 0.
40 C H 3 0 + O 2 = C H 2 0 + HO., 4.28E - 13 7.60 -3528.
41 CH30 + M = CH20 + H + M I.---~,E + 14 0.0 25096.

a R a t e constants are in the form k , = AnT n e x p [ - E J ( R T ) ] , here R is the universal gas constant. Units are moles,
cubic centimeters, seconds, Kelvins and c a l o r i e s / m o l e .
26 B. YANG AND S. B. POPE

0.14 ' ' ' ' I , , , i , , , , I , , , , I , , t ,

0.12

0.10

.2
¢J H20
~ 0.o8- Jk CO2

~1 ~ 02
0.06- v CH4

"~ "t -- co

0 . 0 ~

0.0

o.oo . . . . . T. ,~.-/-~ : : " ~

10 20 30 40

Time step
Fig. 1. Average mass fractions of 1420, CO 2, O~, CH~. and CO as functions of time.

where M is the number of particles in the can see that when the ratio of the subtime step
reactor. This is a measure of the error between ~t to the mixing time scale rmix decreases, the
the solution from the direct integration of the error uniformly decreases. Figure 4 is a plot of
original full coupled equations and the solu- the time average value (from time step 1 to the
tion from the direct integration of the zero- given time step) of the average error ( 1 4 J o e -
order splitting system. As shown in the plot, we $o(~l). After some time, this value becomes a

O . O 0 4 1 1 ' l l l i l ' l i ' l , l i l i , l i l l ,

U) 0.003
e-
.2
¢J

0.002 -----~---- o
M

°°°1~,~c~ ~
oo o o ~ - : : ~ ~ ~ o ~ 30 4o-~ 50
Time step
Fig. 2. Average mass fractions of H, OH, and O as functions of time.
COMPUTATIONS OF COMBUSTION CHEMISTRY 27

0.012 . . . . ~ . . . . ; ' ' ' ' i , , i . . . .

^
0.010

o 0.008
/~j/

/; ' ~ ~ ~

,,
~ I

#
i

.e- ~ o.oo6 /
V - - 6t I x=~. = 0.6
....... 6 t l x=, = 0.12
0.004-
............ ~t I x . u = 0.06
.......... 6t I ~.k = 0.03
0.002- ',

. . . . 1'0 . . . . 2"0 . . . . 3'0 . . . . 4'0 ' " "

Time step
Fig. 3. The average error. (l(h;,, - d,;).l), as a function of time.

constant one. Errors of other methods also In these figures, the first subscript of g)
behave in the same manner as those of the indicates D is the direct numerical integration
zero-order splitting system. So we can use the of the equations, M is the solution based on
time average value of the errors at time step 51) manifold points (i.e., Eq. 46). The remaining
as a measure of errors of different methods. subscripts indicate the splitting scheme. C is
These are plotted in Figs. 5-7. the full coupled mixing and reaction (no split-

0.008 a , , i ~ , , , J I i , , l I t l i i ' . . . . l

0.007
A
o
o 0.006
t

g 0.005
! ....... 6tlx~ = 0.12
V ............ 6t I -emil = 0.06
... 0.004-
O ...... 6t I x,.~. = 0 . 0 3

0.003-
~ -
a 0-0022
o

E
1,7. 0 . 0 0 1 "
._~--::: :::,: .....................................................................
0.000
10 . . . . 20 . . . . 3~ . . . . 4'0 . . . .
Time step
Fig. 4. T h e t i m e a v e r a g e o f t h e e r r o r . (l&/)(. - d,o.I), as a f u n c t i o n o f time.
28 B. YANG AND S. B. POPE

, , , i i i i i I i i i i i

1 0 .2
< 1 % 0 " % s q I;. //~]~
" ~ - - < I Ooc'¢os, l > slope 3/2 J ~ ' J
- < > - - < '~oo'~,.'> ,-~J/
+ < I eoc ~oPca I > ~. " " .." slope I

10 "=

10 -4.

., , . . . . . ,~.,

8t I "Cmlx
Fig. 5. Variations of splitting errors with the change of the ratio, (~//7"mi x.

ting), Strg is the system obtained from Strang's the system obtained from the first predictor-
sequential splitting method (mix-react-mix), S1 corrector method, PC2 is the system obtained
is the system obtained from first-order splitting from the second predictor-corrector method.
method, SO is the system obtained from zero- There are basically three different types of
order splitting method (mix-then-react), PC1 is errors: the splitting errors shown in Fig. 5,

I i i i i I I 1 1 i i I i I

10 -a.
<l¢ostrI - Ous~l>
;" < I ~osl " OMsl I >
< I eos* " ~).so t >
slope I

10"3

10 .4.

8t I ~ u
Fig. 6, Variations o f chemistry errors with the change o f the ratio, Bt/Tmi x.
C O M P U T A T I O N S OF C O M B U S T I O N C H E M I S T R Y 29

i , , i i i i J ',

< I Ooc - Ousmu I>

1 0 2-
~" < 1 % c " Ousl I >
~ - < 1 % c "Ouso I >

/ / . 4 ~ ~ .-''" ?lope 1

1 0 .3

10"4

~-2 10 "~
~t I ~m~=
Fig.7. V:tri~tion,,~l o~cr-;dlerror'.~ith the ch:tngcof the ratio. ~t/r,.,~.

indicated by ([qS/~u - ,ht~,,,]), here ,~t refers to of the real parts of eigcnvalues of the Jacobian
different splitting schcmc; the chcmistr 3, errors matrix J. A typical value of T,ma, is about
in Fig. 6, showing how wcll thc coupling be- 5 x Ill s (s). The formal order of accuracy is
tween chemistry and mixing based on manifl~ld to bc cxpcctcd only when (6t/r~m~,) is small
points is doing, they arc indicated by (l(h~,,, - compared to unity, which is not the case here.
4)Mml); the overall errors in Fig. 7, showing thc It is c,fident from Fig. 5 that the splitting error
errors between the solutions from the coupling is smallest for the simple zero-order splitting
method based on manifold points and those (SO), i.e.. mix-then-react. For the smallest time
from direct integration of the original full cou- step shown, these errors are about 0.01% of
pled equations, they are den(~ted by ([(b~,(.- the major species concentrations.
,t,~,,.l>. For selected splitting schemes, Fig. 6 shows
It may be seen from Fig. 5 that all of the thc error incurred in using the manifold
splitting schemes convergc. That is, the errors method (M) compared to direct integration
incurred tend to zero as the time step (St tends (D). In this case there is no reason to suppose
to zero. The slope of the curve corresponding that crror tends to zero with 6t. For, even in
to Strang's algorithm is about !.5, the others the limit 61---, 0, the compositions lie some
are about unity, which means that Strang's distance from the manifold, and hence the
method is of order one-and-a-half accurate, linearization about manifold points involves
the other methods, arc first-order accurate. some error. (This error does converge to zero
Theoretically, Strang's algorithm and the prc- as r* tends to zero.)
dictor-corrector methods arc of second-order From Fig. 6 we can see that for the zero-
accuracy. The drop in the ordcr of accuracy order and first-order splitting methods, errors
can be attributed to the fact that the time step decrease gradually as the ratio of the sub-time
used in the calculations, ~t, is large compared step to that of the mixing time scale decreases.
with the smallest combustion time scales ~,.,~,, For Strang's algorithm, in the region of large
which are represented by negative reciprocals 8t/r.~.~, the error decreases with the decreas-
30 B. Y A N G A N D S. B. P O P E

ing of fit/~,,,~ X, as ~t/~',,~x dccrcascs further, Sometimes, the dimension jumps to 14, then
the error increases, then decreases. decreases to 6 within one time step. The jump
The over-all errors arc shown in Fig. 7. of the dimension to high dimension corre-
These are the errors that we arc ultimately sponds to the increase in the mass fraction of
interested in as they represent errors of the methane as seen in the figure. In this situation,
method compared with direct integrations of what happens in the PMSR is that this particle
the original coupled equations. Errors of the is ejected from the reactor and replaced by an
zero-order and first-order splitting methods de- incoming particle which is a mixture of
crease a s t ~ l / ' r m i x decreases, the slope of the methane and air at stoichoimetric condition
curves is about unity. The behavior of the and room temperature. At these conditions,
ovcr-aU error of Strang's algorithm is the same reactions progress very slowly so the dimension
as that (~f the chemistry error of Strang's algo- of the manifold is high. As reaction and mixing
rithm. The chemist," 3' error is the controlling go on, the composition of the particle moves
factor for Strang's algorithm. ~towards high temperature. So correspondingly,
To give some idea about the change of di- reactions happen faster and the dimension of
mension of the manifold, we.plot the manifold manifold decreases.
dimension n,, versus subtime step fi~r one
particle in the PMSR. It is shown in Fig. 8. 7. C O N C L U S I O N S
Plotted together in this figure is the mass frac-
tion of methane for the same particle. The In this paper, different errors of different solu-
tion schemes for the original coupled equa-
case selected is the first-order splitting method
with five sub-time steps. One can see from this tions of a pair of particles in a pairwise mixing
figure that most of the time, the dimension of stirred reactor have been investigated. The fol-
the manifold is 6. The dimension of the con- lowing conclusions can be drawn:
served subspace is five (four due to element (1) The pairwise mixing stirred reactor (PMSR)
conservation and one due to enthalpy conser- has been formulated as a test case for
vation). Thus n,,, = 6 corresponds to ~ one-di- simplified chemistry schemes. This is a sig-
mensional manifold in the reactive subspace. nificantly richer and more strenuous test

........ , ........ ,~,,,,,,,..J..t,, .... s ......... i~).O S

cE 16-" ....... n-

-0.04 >."
,~0 14- ' t I

E 1,-
~ ,' ,
it ',"
i -11.03
|
"6 ! :'I "6
lo- :: ,; o.o2

E • ,
:'0 fi "0.01

-- . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . ;g'O0

sub-time step
Fig. 8. Changes of the dimension of manifold and methane mass fraction with subtime step.
COMPUTATIONS OF COMBUSTION CHEMISTRY 31

than the PaSR. (In the statistically-sta- formulation and sobaion o f the closest manifold
tionary state, the compositions in a PaSR point problem. This paper was prepared with the
lie on a one-dimensional manifold.) support o f the U.S. Department o f Energy, Mor-
(2) A variant of the intrinsic low dimensional gantown Energy Technology Center, Cooperative
manifold (ILDM) method has been formu- Agreement No. DE-FC21-92MC2906L
lated and used in which the manifold is
defined by a time scale ~'*. Rather than
having a specified, fixed dimension, the REFERENCES
manifold has different dimensions in dif- i. Smookc,M. D., J. Comput. Phy. 48:72 (1982).
ferent regions of the composition space. 2. Smcokc,M. D., J. Opt. TheoryAppl. 39(4):.489(1983).
This results in the manifold method provid- 3. Smooke, M. D., Miller, J. A., and gee, R. J., in
ing an accurate description even in regions Numerical Boundary Value ODEs (U. M. Ascher and
where the chemistry is slow. R. D. Russell, Eds.), Birkhauser, Basel, 1985, p. 303.
4. Dixon-Lewis,G., David, T., Gaskell, P. H., Fukutani,
(3) The "closest manifold point" problem that S., Jinno, H., Miller, J. A., Kee, R. J., Smooke, M. D.,
arises in the implementation of manifold Peters, N., Effelsberg, E., Wamatz, J., and Behrendt,
methods has been formulated as a mini- F., Twentieth Symposium (International) on Combus-
mization problem, and a reliable iteration tion, The Combustion Institute, Pittsburgh, 1984, p.
1893.
algorithm for its solution has been pre-
5. Pope,S. B., Prog. EnergyCombust. Sci. 11:119(1985).
sented and used. 6. Pope, S. B., Twenty-Third Symposium (International)
(4) An exact solution (Eq. 30) is given to the on Combustion, The Combustion Institute, ~tts-
coupled reaction-mixing equation, with burgh, 1990, p. 591.
frozen mixing vector and linearization 7. Norris, A. T., and Pope, S. B., Combust. Flame 100:211
about the closest manifold point (Eq. 29). (1995).
8. Smooke, M. D., Clump, J., Seshadri, K., and Gio-
A numerically stable solution methodology vangigli, V., Twenty-Third Symposium (International)
for this equation has been developed. on Combustion, The Combustion Institute, Pitts-
(5) Different schemes for splitting reaction and burgh, 1990,p. 463.
mixing have been investigated. The time 9. Smooke,M. D., Puri, I. K., and Seshadri, K., Twenty.
step sizes ~t of interest are large com- First Symposium (International) on Combustion, The
Combustion Institute, Pittsburgh, 1986, p. 1783.
pared to the smallest time scale of the 10. Peters, N., in Numerical Simulation of Combustion
chemistry TsmalI. As a consequence, schemes Phenomena (R. Glowinski, B. Larrouturou, Temam,
that are formally second-order accurate ex- Roger, Eds.), Lecture Notes in Physics, Springer-
hibit lower-order accuracy. All schemes ex- Verlag, Vol. 2,d, 1985, p. 90.
hibit first-order accuracy except for Strang's 11. Smooke, M. D., (Ed.), Reduced Kinetic Mechanisims
and Asymptotic Approximations for Methane-Air
algorithm (mix-react-mix) which is of order Flames, Lecture Notes in Physics, Springer-Verlag,
3/2. 1991, Vol. 384.
(6) Over the range of time step size investi- 12. Peters, N., and Rogg, B., (Eds.), Reduced Kinetic
gated, the simplest zero-order splitting- Mechanisms for Applications in Combustion Systems,
mix-then-react has the smallest splitting er- Lecture Notes in Physics, Springer-Verlag, 1993.
13. Maas, U., and Pope, S. B., Combust. Flame 88:239
rors, which is as low as 0.01%. This is an
(1992).
important conclusion because it justifies 14. Chen, J. Y., Dibble, R. W., and Bilger, R. W.,
the current practice in table-look-up Twenty-third Symposium (International) on Combus-
methodologies. tion, The Combustion Institute, Pittsburgh, 1990, p.
(7) The combinations of zero-order splitting 775.
15. Maas, U., and Pope, S. B., Twenty-Fifth Symposium
and the manifold method produce accurate
(International) on Combustion, The Combustion In-
solutions. stitute, Pittsburgh, 1994, p. 1349.
16. Maas, U., and Pope, S. B., Twenty-Fourth Symposium
(International) on Combustion, The Combustion In-
We are grateful to Professors U. A. Maas, T. F. stitute, Pittsburgh, 1992, p. 103.
Coleman and C. F. Van Loan f o r useful sugges- 17. Correa, S. M., and Braaten, M. E., Combust. Flame
tions. In particular Dr. Maas was involved in the 94:469 (1993).
32 B. Y A N G A N D S, B. P C ' P E

18. Douglas, J., J. Soc. Ind. Appl. Math. 3:42 (1955). 23. Coleman, T. F., and Y. Li, Cornell Theory, Ceater
19. Pcaceman, D., and Rachford, H., J. Soc. Ind. Appl. report CTC92-111, Cornell University, 1992.
Math. 3:28 (1995). 24. Land, A. H., FORTRAN Codes for Mathematical Pro-
20. Yanenko, N. N., The Method of Fractional Steps, gramming: Linear, Quadratic and Discrete, John Wiley
Springer-Verlag, Berlin, 1971. & Sons, New York, 1972.
21. Strang, G., SIAMJ. Numer. Anal. 5(3):506 (1968). 25. Van Loan, C. F., IEEE Trans. Auto. Con. AC-
22. Golub, G. t-!., and Van Loan, C. F., Matrir Computa- 23(3):395 (1978).
tions, 2nd ed., Johns Hopkins University Press, Balti-
more, 1989. Receit'ed April 30, 1996; accepted Febn~ary 12, 1997