- Email: [email protected]
Kinetic scheme reduction via geometric singular perturbation techniques
Pergamon PII:
Chemical Engineering Science, Vol. 51, No. 20, pp. 4661-4672, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0009-2509/96 $15.00 + 0.00
S001D-2509(96)00310-7
KINETIC SCHEME REDUCTION VIA GEOMETRIC SINGULAR PERTURBATION TECHNIQUES P. DUCHENE* Institut Francais du Prtrole, BP 3, 69390 Vernaison Cedex, France and P. ROUCHON l~cole des Mines de Paris, Centre Automatique et Syst~mes, 60 Bd Saint-Michel, 75272 Paris, Crdex 06, France (Received 19 M a y 1995; accepted in revised form 3 April 1996)
Abstract--An effective reduction method based on geometric singular perturbation and center manifold techniques is proposed. This method eliminates the fast and stable dynamics and gives the equations describing the slow ones. It is coordinate-free and extends the well-known quasi-steady-state method classically used for kinetic scheme reduction. This yields directly the slow dynamics even if the differential equations are not in standard two time-scale form (Tikhonov form). Application to combustion kinetics including 13 species and 67 reactions is presented and simulations are given. Copyright © 1996 Elsevier Science Ltd Keywords: Model reduction, geometric singular perturbations, center manifold, combustion, kinetics.
1. I N T R O D U C T I O N
Let us begin with a simple reaction scheme. It involves three chemical species Xx, X2 and X3 and three elementary independent reactions: XI k ~ x z ,
X2 k2~,X1, X I + X 2
4- k z x 2 -
2 2 = klx I -- k2X 2
3(1
=
--
klX1 + k2x2 - g k x l x 2
(2) 0 = klXl - k2x2
~k.... ' X 2 + X 3
the xis are compositions of the Xis; kl, k2 and ek are the kinetic constants. The small positive parameter indicates here that the third reaction is slow with respect to the first two reactions. From this reaction scheme, we derive an ordinary differential equation system. It corresponds to the conservation of each species in a closed homogeneous reactor (constant volume, molar volume independent of compositions, constant temperature). The species balances are: "?(1 = - - k l X l
The quasi-steady-state assumption applied to x 2 leads to the following reduced system:
or
Y¢1 = - e ( k k l / k 2 ) ( X l ) 2
(3)
x2 = ( k l / k 2 ) x l .
(4)
with
This system does not yield a correct approximation of the slow dynamics for small e. The reason is the following. From eq. (1), we have d(xl + x2)/dt = - e k x l x 2 . Thus, if eq. (4) was valid, one would have (1 + k l / k 2 ) X l = - e ( k k J k z ) ( X l ) 2.
ekXlX 2
(1)
-~3 = g k X l X 2
where the dots above the xs indicate differentiation with respect to time. Since x3 follows from the chemical invariant X l + X 2 + X s = X ° + x ° + x ° (the superscript 0 indicates initial compositions), we focus on the first two equations of (1).
This stands in contradiction with eq. (3). Here, the quasi-steady-state method cannot be directly applied. A change of coordinates is required first. Instead of using (x~, x2), consider the coordinates (~ = X l + X2, X2).
The first two equations of (1) are transformed to = -- ~k(¢ -- X2)X 2
(6) 3~2 = k l ( ~ -
* Corresponding author.
(5)
x2) -
k2x2.
Notice that ~ corresponds to an approximate lumping of xx and x2 [for more about lumping, see e.g. Kuo 4661
4662
P. DUCH~NEand P. ROUCHON
and Wet (1969); Li and Rabitz (1989, 1993)]. The application of the quasi-steady-state method to x2 leads then to the correct approximation: = -- 8k(~ -- X2)X2,
0 = kl(~. -- X2) -- k2x2,
that reads, in the original coordinates (x~, x2), = - e~-2 (Xl)
with x2 = ~ x l .
(7)
The difference with eq. (2) lies only in the coefficient (1 + kUk2) multiplying the derivative of x~. The algebraic equation, giving x2 with respect to x~ remains unchanged. For this simple system, the 'good coordinates' (~, x2) are easy to find. Van Breusegem and Bastin (1991, 1993) explained how to find such a linear change of coordinates when the reactions can be separated into fast and slow ones. For more complex systems, such as the combustion system of Section 5, this is no longer the case: there is no clear decomposition into slow and fast reactions and thus the change of coordinates is probably nonlinear. The need for model reduction is motivated by two purposes: a qualitative analysis of the model and decreasing the simulation times, in particular when kinetics is coupled with transport. In this paper, we propose a reduction method that overcomes the difficulty of finding the coordinate change. The obtained reduced model [eqs (8) and (9) of Section 2] is directly given in the modelling variables. A deeper analysis yields the following explanations. The quasi-steady-state method is not coordinate-free: it can be applied only to systems written with special coordinates, called Tikhonov coordinates in the sequel (see Section 3.2). Our extension provides a coordinate-free reduction method. This means that no change of coordinates is required. This extension admits a clear and rigorous justification within the theory of perturbations of dynamical systems (Arnold, 1980; Guckenheimer and Holmes, 1983; Wiggins, 1994). More precisely, we use here the theory of singular perturbations (Tikhonov et al., 1980), its geometric setting due to Fenichel (1979) and a special version of the center manifold approximation lemma (Carr, 1981); see also Ehrenstein (1990) and Ehrenstein and Nastoll (1991) for application of center manifold techniques to reduction of chemical kinetic schemes, in a slightly different context. Another approach to model reduction is proposed by Maas and Pope (1992). When the system admits two time-scales, one can prove that this method is equivalent to ours.
The content is the following. In Section 2, we describe the practical aspect of our method. The mathematical framework is the object of Section 3: we recall the classical Tikhonov form of singularly perturbed systems; we explain the geometric definition of slow/fast systems (Fenichel, 1979); the system reduction via invariant manifold techniques is derived on the basis of an approximation lemma that can be found in Carr (1981). Section 4 is devoted to a special
class of chemical systems where the two time-scales are due to fast and slow reactions. In Section 5, the application of our reduction method to a combustion scheme considered by Maas and Pope (1992) is presented. For this system, the quasi-steady-state method yields a rather poor approximation. 2. T H E R E D U C T I O N
METHOD
For readers not interested in technical developments, we just present here the formulae that are required for applying our reduction method to a dynamical system of the form = v(x)
with x = (xl . . . . . x~), v = (Vl. . . . . v~)
n being the dimension of x. Assume that this system admits two time-scales: • a slow time-scale that can be described by n~ < n variables. • a fast and (hyperbolically) stable time-scale that can be described by n I = n - ns variables. The integer ns will be the number of differential equations of the reduced slow system. The integer nf corresponds to the number of derivatives that are set to 0 in the quasi-steady-state method, i.e. the number of algebraic equations that must be solved. As in the quasi-steady-state method, we decompose the state x~R" into two sets of components, x = (xs, xs.) of dimension n~ for xs and nI for xf. Similarly, the system is decomposed into two parts
Yc, = v,(xs, xs),
~s = vs(x~, xs)
with v = (v,, vs). This decomposition is up to the choice of the user. We give below some guidelines to make this choice. In particular, the eigenvalue structure of the Jacobian matrix (cnvi/c3xs) provides an estimation of nf. Then, the reduced system is as follows:
2~ = C(x,, ho)" vs(x~, ho) 0 = Vs(X~, ho) xf=ho-
(8)
~[(c?vI~ 2 + ~l)f'~[)s]-xovf)7}
tLk~.U
'~, a,
ex.d"
Us(Xs,
ho)
where the matrix C(xs, ho) is given by
C = I . , - c,~xI L k ~ x #
+
~x,
axsJ ~
(9)
with 1,, the identity matrix of order n~. All the partial derivatives are computed at x = (xs, ho). Some algebraic manipulations provide another equivalent expression for C that might by preferable in practice:
C=[lo " + ax~,t~s)
77]x,j
-' "
Notice that formulae (8), applied to eq. (1) with x, = (x~, x3) and x I = x2, lead to the correct approximation, the correction factor (1 + kUk2) of eq. (7)
4663
Kinetic scheme reduction corresponding then to the matrix C - 1 (up to terms of order 1 in e). In fact, the only difference between eq. (8) and the equation resulting from the quasi-steady-state method lies in the correction matrix C. For systems with v, = el,, such as system (6) with x~ = ~, the matrix C is close to 1,, up to terms of order 1 in e. In this case, our method coincides, up to terms of order 2 in e, with the quasi-steady-state one. In order to avoid the computation and the inversion of the Jacobian matrices contained in eqs (8) and (9), it may be useful to approximate eq. (8) via the following differential-algebraic system: dx, dt
v~(x~, xs)
0 = vf(x.~, ho)
(10)
0 = vf(x, + 6vs(x~, x f), ho + 6vf(x,, x f))
where 6 is a very small parameter and (x~,ho, xy) are unknown. Some elementary computations show that if 6 ( ~ 0) is at least of order 2 in e then the variable x I in eq. (10) satisfies eq. (14) for p = 1. Such differentialalgebraic systems of index 1 can be easily solved numerically via Gear-like schemes and standard codes of public domain, such as DASSL or L S O D E / L S O D I [see e.g. Brenan et al. (1989)]. As far as we know, formulae (8) and its approximation (10) are new and have never been introduced elsewhere in such a simple formulation. Such equations are, in general, more complicated than the original ones. This difficulty can be overcome, if the formulae are approximated by off-line interpolations using only xs as independent variable. Then, the simulation just consists of integrating an explicit differential system ~s = P(x~), where function P results from such an interpolation. Note that, for the computation of the reduced slow dynamics, it is not required to exhibit e as a function of the model parameters. We just need to know that the system admits two time-scales, that the dimension of the fast part is n I and that it is asymptotically stable. This point and the simplicity of the formulae constitute the major interest of the method proposed here. In practice, we have to define n I (or equivalently n, = n - nl). Since our reduction method amounts to eliminating fast and stable time-scales, we may be guided by two criteria: 1. What is the fastest time-scale that cannot be eliminated? 2. The time-scales of the reduced model (slow timescales) must be significantly slower than those that have been eliminated. The answer to the first question is given by the following considerations: • The kinetic model may be part of a bigger model including, for example, transport. In this last case, one must keep in the reduced model the
time-scales that are of the same order of magnitude as those of transport phenomena (and of course the slower ones). We denote by ~ this typical time constant. • The time of evolution of the system, T, may be a priori known. In this case, one eliminates those time-scales that are significantly faster than T. The characteristic times of the system, (Ti) i 1 . . . . . may be estimated with the help of the eigenvalues of the Jacob]an Ov/Ox, (2i)i= 1...... : zi ~ 1/[Re()~i)[.
To satisfy the second criterion, one computes the eigenvalues along different trajectories of the system. Then one tries to separate the characteristic terms into two groups, (zl,r2 . . . . . r,s) and (r,~+l . . . . . T,), such that "gl ~ "C2 ~
"'" ~ ~'n! <<0~ ~ "Cnt+l ~
"'" ~ "(n.
The eigenvalues relative to the first group have a large, negative real part: they correspond to fast and stable phenomena. The eigenvalues of the second group correspond to the slow part of the system and thus are close to zero. After having determined ny, one has to parametrize the reduced model, i.e. choose a decomposition of x into (xs, xy). In view of the algebraic equations of (8), it seems reasonable to choose the x I variables among those that are at or near quasi-steady-state. A species is in quasi-steady-state when its creation rate is (almost) equal to its destruction rate. This may be evaluated by the ratio Icreation - destruction] creation + destruction for each species i. This ratio is based on the classical decomposition of vi v,(x) = ~,? (x) - v~- (x)
where the creation v [ ( x ) corresponds to reactions producing i, and the destruction v~(x) to reactions consuming i. The species for which this ratio is small (compared to 1) may be considered as being in quasisteady-state. This ratio can be numerically calculated along the trajectories of the system. We choose the x I variables among those for which this ratio is small. Other procedures for defining ns and (xs, x r) can be used. We just have sketched here a simple one that gives, for the combustion system of Section 5, an efficient reduced model of dimension one. Let us now explain how we have obtained eqs (8) and (9). 3. G E O M E T R I C SINGULAR PERTURBATIONS
3.1. Model reduction The following question underlies this paper: what is model reduction? A possible and reasonable answer to this question for dynamical systems of the form .;: = v(x), x ~ ~" is displayed in Fig. 1. Reduction can
4664
P. DUCHI~NEand P. ROUCHON
Fig. 1. Trajectories of ~ = v(x) are caught by an attractive invariant manifold Z.
The terminology 'singular perturbations' comes from the fact that the small parameter e multiplies the highest derivative (here d x l / d t ). More details on this classical standpoint can be found, e.g. in Tikhonov et al. (1980). In the sequel we always consider the time-scale t and approximations of trajectories for t~ [0, I/z] and with e > 0 but close to 0. Assume that the fast part is hyperbolically stable, i.e. the sub-system x.r = vf(x~,xr,e) with xs fixed admits (locally) an equilibrium with characteristic exponents (eigenvalues of Ovl/t3x I at this equilibrium) having a strictly negative real part. Then the slow approximation is obtained by the quasi-steady-state method (Tikhonov et al., 1980):
2~ = evs(xs, x f, ~) be defined via an attractive invariant manifold T,. A manifold E is invariant with respect to the vector field v, if v is tangent to E: if a trajectory starts on Y., it remains on Z for all time. E is called (locally) attractive if any trajectory, t ~ x(t) of ~ = v(x) (starting near Y.), tends to E as t tends to + ~ . In this case, the reduced system corresponds to the restriction of the vector field v to E. This restriction is meaningful since Z is invariant. It seems then natural to approximate trajectories of the complete system ~ = v(x) by trajectories on E. In fact, such an approximation is proved to be valid when, roughly speaking, the dynamics transverse to ~z (the dynamics that are neglected) are faster than the dynamics on E [see Arnold (1980) and Hirsch et al. (1977) for more details on normally hyperbolic invariant manifolds]. Reduction means then restriction of the dynamics to an invariant attractive manifold E. In practice the equations defining E are not known. Only the vector field, v, is explicitly available. The main difficulty is thus to obtain the equations of E or, at least, good approximations of them. 3.2. Tikhonov normal form of slow~fast dynamical system Here we restrict the study to a special class of dynamical systems, depending on a small parameter e and having an attractive invariant manifold E~. Such systems are characterized by dynamics transverse to E~ which is much faster than the dynamics on E~. This corresponds to systems having two time-scales: a fast and asymptotically stable one and a slow one (stable or unstable). The first kind of such systems is singularly perturbed systems of the so-called Tikhonov form,
dx~ d----t= e v~(x~,xy, e),
dxf dt - vI(x~'xI'e)
(11)
where 0 < 5<< 1. Very often, such systems are written with the time-scale z = st: dx~ d--~- = vs(x~, x I , e),
dx I e-d~t = vAx~, x l , 5).
0 = Vs(Xs, x s, e). The algebraic equations, vy = 0, correspond here to an approximation up to terms of order 1 in e, of E, equations. These coordinates (x,, xl) where the quasi-steady-state method applies, and where the vector field v is quasi-vertical (see Fig. 2) are clearly very specific. 3.3. Geometric settin9 Since we are interested in developing a reduction method that avoids changing coordinates, we need a coordinate-free point of view. The geometric definition of singularly perturbed systems due to Fenichel (1979) is as follows. Consider the dynamical system :~=v(x,e),
x~R", 0~
(12)
This system is said to have two time-scales, a fast and asymptotically stable one and a slow one, if, and only if, the following two assumptions are satisfied (Fenichel, 1979): (A1) for e = 0, eq. (12) admits an equilibrium manifold of dimension n~, 0 < n~ < n, denoted by ~0. (A2) for all Xo e Zo, the Jacobian matrix, c3v/Ox[txo.o~ admits n I = n - n~ eigenvalues with a strictly negative real part (the eigenvalues are counted with their multiplicities). This definition is illustrated in Fig. 3. Assumption (A1) implies that the velocity v(x, ~) is large everywhere except for x in a neighborhood of Y'o where v is small and of order unity in e. (A1) and (A2) imply that, for Xo~Y~o, the kernel, Uo(Xo), of the linear operator t~v/Oxllxo, O~ coincides with the tangent space of Eo at Xo. The linear space E~(xo) corresponding to the eigenvalues with real negative part satisfies
E~o(Xo) @ Eo(Xo) = R".
4665
Kinetic scheme reduction
a partition of the state x into two groups of components, x = (x,, x:), with dim(x~) = n~ and dim(x:) = n:, such that the projection of E~ on the x:coordinates is a local diffeomorphism around Xo. The equations defining E, around Xo admits the following approximation:
v( )
:(Xs, x:,
x , = ho(x,, e) + 0(~)
=
x d = ho(x~, e,)
(13)
x. tkVx,) + Fig. 2. In Tikhonov normal form, the trajectories of eq. (11) admit a vertical part before being caught by the attractive manifold E~.
axsj
&sl,.,.ho.,,=),,:)J
• v.,(x,, ho(x~, g), g) + 0(8 2)
where ho is the solution of vs(x~, ho(x~, ~:),~:) = 0. The trajectories of the perturbed system are captured by a trapping region around Eo and entering with a direction nearly parallel to E~o(Xo). Fenichel (1979, part of theorem 9.1) proves the following result• It asserts, for e small enough, the existence of a slow invariant attractive manifold Z~ for the perturbed system (12).
Theorem 1 (Feniehel, 1979): Consider eq. (12) satisfying (A1) and (A2). Then, for every open and bounded subset ~o of Y~o, there exists an open neighborhood Vo of f~o in ~", such that, for e positive and small enough, the perturbed system (12) admits an attractive invariant sub-manifold Z~ contained in Vo and close to Eo. We are interested in approximations, up to terms of order 1 in e, of slow trajectories for t e [0, l/e]. Thus we need an approximation up to terms of order 2 for the slow dynamics: errors like e 2 integrated over t e [0, I/e] will produce for t = I/e, distortions of magnitude less than or equal to ~: (e2 x(1/,:)= e). This means that an approximation, up to terms of order 2 in e, of £, equations is needed. 3.4. Approximation of Z~
Proposition 1. Consider eq. (12) satisfying (AI) and (A2). Take Xo ~ Zo. Then locally around Xo, there exists
v(x, e)
"" (""
!
In some generic sense, that can be easily defined, almost every partition of x in (Xs, x:) is admissible for eq. (13). But in practice, the coordinates x~ must be chosen such that the projection of Z~ onto the x:space is as well conditioned as possible•
Proof of Proposition 1. By the approximation lemma of invariant manifold, and its version for two timescale systems (Carr, 1981, theorem 5, p. 32), the equation x : = hp(xs, e) is an approximation up to terms of order p + 1 in e. of X~ equations, if hp satisfies the following equation stating the invariance of E,::
~h_ Vz(x..h,(x..e).e)=-" v.(x..h,(x..e). ~) + ~Xs (x. ~)
O(e'+b
(14) The zeroth order term, ho(x,,e), is given by the quasi-steady-state approximation
v:(xs, ho(x~, e), e) = 0 because v,(xs, ho(xs, e), e) = O(e). The first-order term, hi is written as hi = ho + epl where pl is a smooth function of xs that will be
x
v(x, O)
JX
ose to
1 _-:...y.o
Fig. 3. The geometric definition of singularly perturbed systems due to Fenichel (1979).
4666
P.
DUCH&NEaInd
wo
+ EPl)
4 + &PI, 4 = 8%
x u,(x,, h&s, E)+
&PI,
ROUCHON
Notice that approximation (8) applied to a system already in Tikhonov normal form yields an approximation of the dynamics on C,, up to terms of order 3 in E [use eq. (13) and replace o, in eq. (13) by EV, for systems like eq. (1 l)].
adjusted in order to verify (14) with p = 1, i.e. V/(X,, hhs,
P.
(I,, S)
4 + W).
(15) 4. CHEMICAL SYSTEMS WITH SLOW AND FAST
First, we develop the left side of the above equality. Noticing that u~(x,, ho(x,, a), E) = 0, we obtain
U&G,ho(&, E)+ &PI,E)=
au,
Eaxf
Pl + W).
(&.h&E),
@
Developing the rightmost term, we get r,(X,, ho@,, 8) + &PI, E) = Vs(xs, ho(xs, 81,E) +c!S ax, It
PI + ok+). (x,,h&.
REACTIONS
In this section we consider the case of kinetic systems where the reactions can be split into one group of slow reactions and one group of fast reactions. The reduction of such systems has been studied by Van Breusegem and Bastin (1991, 1993). We show here that the reduced models deduced from their calculations and the reduced model obtained from eq. (8) coincide, up to terms of order 2 in E. Consider the reaction system
E).c)
1 =
is a first-order term in &(0(s)) because %(&, ho@,, s), e) = O(E).
We now develop the right side of equality (15). d(ho + &Pl) r&G, ho(x,, c) + EPl, e) ax, (X,X E)
zdV(X,
E)
(16)
where x is the vector of concentrations, A is the stoichiometric matrix, v(x) is the vector of the chemical rates and E is a small number (0 < EK 1). Assume that v can be split into v = (EC,, a,), where ~6, (resp. 5,) is the vector of the chemical rates of the slow (resp. fast) kinetics. System (16) can be rewritten as k = A, E&(X) + A/d,(X).
=- aho
(17)
U&G, ho(xs, s) + &PI, ~1+ W2)
ax, (I,,E) =- aho a xs kc)
Let A,, be a square sub-matrix of A/ such that its rank is equal to the rank of A,. Let
u,(xs,
ho(xs> ~1,~1 A, = Pl + W2).
and Equating both sides of eq. (15) and using the relation
aho
-lavf av,_
(>
dx, = - ax, ax,
Then eq. (16) reads
we obtain J, = A,, at%,
XJ) + A,, cf(x,, x/).
The change of coordinates (~8, x/) H (5 = x, - &(AJixs>
x,)
hence eq. (13). The matrix leads to a quasi-vertical vector field. It yields the following Tikhonov form:
au,2 auf au,
ax, +ax’ax (-> s
f
g = (A,, - &(A,,)-‘A&G
is invertible: this results from (Al) (A2) and from some n elementary calculations. The reduced slow model (8) comes from eq. (13) and u,(x,,
ho +
EPI, ~1= t’stxs, ho, 6) + av, axr
+ O(E2).
&PI @..h,,e)
1, = Afs.$, + A,,Cf. Assuming that eigenvalues of A,,-diT,/dx, have strictly negative real parts, the quasi-steady-state method can be applied and leads to the following slow system: i = CA,, - &(A,,)-‘A~,l~~s
0 = A,,
EC3
+
A,,Cf
4667
Kinetic scheme reduction
The above formulae give the initial conditions (~o, y~) of the slow sub-system from the ones of the complete system (x°, x~).
Pulling back into the original coordinates ( x . x:) yields
l(~vs~-'~Vs 1 X, = II,, - a~s(af:)- \ ~ x : J -~XsJ 5. CASE-STUDY
× [Ass -- Asf(Aff)- XAfs]r.fs 0 = A:s~fs(xs, x:) + A:: fy(xs, xf). This approximation differs from eq. (8) only via terms of order 2 in e: at order 1 in e both approximations coincide. The main feature of eq. (8) in this case relies on the following point: one does not need to know precisely the decomposition between slow and fast kinetics. This can be an advantage if one is only interested in the slow dynamics. Notice that, even if the two timescale structure is not due to fast and slow kinetics, formulae (8) are still valid. However, for such special slow/fast kinetic schemes, one can, as displayed on Fig. 4, describe how and where trajectories starting far from the slow attractive manifold Z~ are caught by X~. The quantity = x., - Asf(Aff)- ~x: plays an important role here: it is a quasi-invariant during the fast transient, before arriving near Y~. As displayed on Fig. 4, the trajectory of eq. (18) with initial condition (x°, x}), becomes close, after a time interval of length of order 1, to the trajectory of the slow system here above, starting at the same time 0 with initial conditions (fro .~) solution of .~0 _ A s f ( a f f )
1 ~ = No __ A s f ( A f f )
~:Ix °, x~,
lx~,
~.) = 0.
In this section, we apply the above reduction method to the CO/H2/air combustion system also considered by Mass and Pope (1992). Thirteen chemical species are present: N2, CO, H2, 02, H20, C02, OH, H, O, HO2, CHO, H202, CH20. Sixty-seven reactions are taken into account. Their list is given in Table 1. Notice that N2 is inert. The system will be supposed closed, isobaric and adiabatic. It is modeled by a system of differential eel ,ations whose dimension is nine (the system is ist,_aric, adiabatic, and four atomic invariants are present: C, H, N and O). The determination of n: and of the partition (x~, x:) is based on a sample trajectory starting at t = 0 from the initial conditions given in Table 2. This corresponds to the beginning of the explosive phase that lasts about T = 10 -2 s, the time-scale of interest. The eigenvalues at different times along this trajectory are given in Table 3 (they have been sorted in increasing order). All eigenvalues, except one at the beginning of the trajectory, are real and negative. For the reduction purpose, one may notice that the eighth and ninth eigenvalues are well separated for t >/2 × 10 -4 s): the ratio of the latter to the former (when both are negative) ranges from 2 to more than 200 for t/> 2 × 10 4. This indicates that we should be able to build a one-dimensional reduced model. Thus n: = 8 and n~ = 1. For t close to 2 × 10 -4 s the separation is not very large (ratio of 2). This may result in poor accuracy of the reduced model. If so, the
i x' / . ' / / x o .................
......
.¢. ..........
............
.,
0
b
,,A--
-0 Xs
. ~
0 Xs
Fig. 4. F o r s y s t e m (18), fast t r a n s i e n t s a r e a l o n g x~ - A ~ : ( A : : ) - ~ x ~ = c o n s t a n t .
4668
P. DUCHI~NEand P. ROUCHON Table
1. Reaction
list
[mass
action
law
and
kinetics
constant
given
k = oT~exp(-Ea/RT)]
ko O2 + H --, O H + O O H + O ~ 02 + H H2 + O ~ O H + H O H + H ---, H2 + O H 2 + OH ~ H20 + H H 2 0 + H ~ H2 + O H OH + OH ~ H20 + O H20 + O ~ OH + OH H+H+M~Hz+M H2+M~H+H+M H+OH+M ~ H20+M HzO+M ~ H+OH+M O+O+M-*O2+M O2+M~O+O+M H + 02 + M ~ HO2 + M HO2 + M -.-* H + 02 + M HO2 + H -," OH + O H O H + O H ~ HO2 + H HO2 + H -~ H2 + 02 H2 + 0 2 -'¢ HO2 + H HO2 + H --, H 2 0 + H H 2 0 + O ~ HO2 + H HO 2 + O --* O H + 02 O H + O2 --' HO2 + O H O 2 + O H --' H 2 0 + 02 H 2 0 + 02 "-* HO2 + O H HO2 + HO2 ~ H202 + 0 2 O H + O H + M ~ H202 + i H2OE+M ~ OH+OH+M H 2 0 2 + H -'* H 2 + HO2 H2 + HO2 --* H202 + H H2O 2 + H -~ H 2 0 + OH H 2 0 + O H ~ H20 2 + H H 2 0 2 + O ~ O H + HO 2 O H + HO 2 ~ S 2 0 2 + O H 2 0 2 + O H ~ H 2 0 + HO2 H 2 0 + H 0 2 -'* H202 + O H CO + O H ~ CO2 + H C O 2 2r- H ~ CO + OH CO + HO2 ~ CO2 + O H CO2 + O H ~ CO + HO2 CO + O + M ~ C O 2 "~- M CO2 + M ~ CO + O + M CO + 02 --* CO2 + O CO 2 + O --* CO + 0 2 HCO + M -~ CO + H + M CO + H + M ~ H C O + M HCO + H ~ CO + H2 CO + H2 --' H C O + H HCO + O ~ CO + O H CO + OH ---, H C O + O HCO + O ~ CO2 + H CO2 + H --* H C O + O H C O + O H ~ CO + H 2 0 CO + H 2 0 ~ HCO + O H H C O + 02 ~ CO + HO2 CO + HO2 --' H C O + O2 C H 2 0 + M -~ H C O + H + M H C O + H + M -* C H 2 0 + M C H 2 0 + H --* HCO + H2 H C O + H2 --* C H 2 0 + H C H 2 0 + O ""¢ HCO + O H HCO + OH ~ C H 2 0 + O C H 2 0 + OH ~ HCO + H 2 0 HCO + H 2 0 -'-' C H 2 0 + O H C H 2 0 + H O 2 --* HCO + H 2 0 2 H C O + H 2 0 2 --~ C H 2 0 + HO2
2.00 x 10 TM 1.47x 1013 5.06 x 104 2.24 × 104 1.00 x 10 s 4.46× l0 s 1.50x 109 1.51 × 101° 1.80x 10 is 6.98 x l018 2.20x 1022 3.80×1023 2.90 x 1017 6.78 x 1018 2.30x 1018 2.66 x l018 1.50 × 101. 1.63 X 1013 2.50x 1013 8.39 x 1013 3.00x 1013 3.29 × 1013 1.80 × 1013 2.67 × 1013 6.00 × 1013 8.97 x 101'* 2.50× 1011 3.25 X 1022 2.11 x 1024 1.701 1012 9.35 x 1011 1.00x 1013 2.66 × 1012 2.80 × 1013 6.80x 1012 5.40x 1012 1.32× 1013 4.40× l06 6.12 x 10 s 1.50x 101'* 2.27 x 10 is 7.10× 1013 1.69 × 1016 2.50 x 1012 2.55 x 1013 7.10 x 1014 1.07 x 1015 2.00 x 1014 1.17 x 1015 3.00 x 1013 7.72 x 1013 3.00 x 1013 1.07 x 1016 1.00 x 1014 2.601 1015 3.00x 1012 5.21 × 1012 1.40x 1017 2.62 x 1015 2.50 × 1013 1.82 x 1012 3.50x 1013 1.12 × 1012 3.00x 1013 9.71 x 1012 1.00 x 1012 1.32x 1011
fl 0.00 0.00 2.67 2.67 1.60 1.60 1.14 1.14 -1.00 -1.00 -2.00 -2.00 -1.00 -1.00 -0.80 -0.80 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -2.00 -2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.50 1.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
E~ 16,791 502 6282 4395 3296 18,451 96 17,101 0 104,137 0 119,280 0 118,563 0 49,250 1003 37,738 693 55,603 1720 55,460 -406 52,618 0 70,030 --1242 0 49,393 3750 21,783 3583 73,469 6401 22,547 1003 34,155 -740 22,475 23,574 83,524 --4538 120,952 47,769 54,696 16,791 2054 0 89,424 0 87,537 0 110,753 0 104,543 0 34,513 76,431 -13,566 3989 18,152 3487 15,764 1194 30,477 8001 4132
by
Kinetic scheme reduction
4669
Table 2. Initial conditions of the complete system for the sample trajectory (mass fractions, temperature) CO Hz Oz H20
0.21 4.9 x 10 4 0.12 4.5x 10 -3
COz OH HO2 H202
0.11 3.7 x 10 -s 5.9 x 10 -5 7.5 x 10 --6
CH20 H O HCO
l a x 10 -n
N 2
6.0 X 10 -6
Temp.
0.55 1076K
1.4x 10 -4 2.3x 10 -7
Table 3. Eigenvalues of the CO/HE/air system at different times along the sample trajectory t=0s -
7.0
t=10-gs
-- 4.4 x 107
-
- 4 . 1 xl06 - 1.7x106 - 6 . 8 x 105
- 1.2x107 - 3 , 0 x 106 - 2 . 1 xl06
- 8.4x 106 - 8.0x 106 - 3.0x 106
-
3.9 x 104
-
1.4 x 1 0 5
-
2.6 x 105
-
1.6
x
-
1.9
x
-
1.2 x 1 0 4
-
4.7
-
1.6 x 1 0 5
-
5.5
× 10 5
--
1.0
X 10 6
8.8 × 1 0 3 - 1.2x103 1.3 x 104 --
-
x
10 4
1.8 x 1 0 4 --6.3x 103 1.5 x 104 --
1.9 x 1 0 7
t=5xl0-4s
- 1,2xl06 - 1,2x106 - 1.5x105
10 6
-
t=3xl0-4s
- 9 . 1 x 105 - 4 . 2 x l05 - 1.2x105
x
9.6 x 106
t=2xl0-4s
5.4 x 1 0 4 - 2 . 7 x 104 1.1 x 104 -
10 6
2.2 x 105 - 2.6x 104 - 5.6 x 1 0 3 --
6.1 x 1 0 7
10 6
- 4.3 x 10s - 6.0x 104 - 4.7 x 1 0 3
t=10-3s
t=10-2s
7 . 3 x 107 - 1.9x107 - 4 . 6 x 106 - 2 . 6 x 106
9.8 x 107 - 5.4x 107 - 2.4x 106 -2.3x106 8.8 × 1 0 5 5 . 4 x 105 - 2.6 x 105 - 6.8 x 10'* - 2.5 x 1 0 2
-
-
1.2
× 10 6
--
1.1
x
10 6
- 4.3 x 105 - 7 . 4 x 104 - - 1.6 x 1 0 3
Table 4. Creation/destruction ratios of the CO/H2/air system at different times along the sample trajectory t=0s HCO HO 2 OH H2Oz H CHzO H20 O 02 H2 CO CO2
0.00 0.00 0.00 0.51 0.01 0.08 0.93 0.08 0.58 0.74 0.91 0.99
t=10-4s 0.00 0.00 0.00 0.00 0.02 0.00 0.83 0.10 0.46 0.71 0.91 0.99
first-order reduced model of P r o p o s i t i o n 1 can be i m p r o v e d by adding higher-order terms (approximation l e m m a of i n v a r i a n t manifold for two time-scale systems used in the p r o o f of P r o p o s i t i o n 1). The ratios
t=2xl0-4s 0.00 0.00 0.00 0.01 0.01 0.01 0.03 0.05 0.19 0.19 0.88 0.94
(3)
t=3xl0-4s 0.00 0.00 0.00 0.00 0.00 0.00 0.06 0.01 0.04 0.00 0.48 0.50
t=5xl0-4s 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.03
the one-dimensional quasi-steady-state model
Yc~ = vs(x~, vy)
(19)
0 = Vs(X.~,x j-)
(20)
with n.~ = 1 a n d x~ = COz (dotted line). Icreation - destruction l creation + destruction at different times along this trajectory are given in Table 4 (they have been sorted in increasing order). The ratio for N2 is not given since this species is neither created n o r destroyed. Since the largest ratio c o r r e s p o n d s to CO2 (or CO), we can set x~ to the mass fraction of CO2 (or C O as well). The simulations of Figs 5-10 c o m p a r e three models: (1) the complete model of dimension nine (solid line), (2) the one-dimensional reduced slow model resulting from formulae (8), with ns = 1 a n d xs = CO2 (dashed line),
In these simulations, the starting value ofxs = C O 2 remains u n c h a n g e d between the three models. F o r xs = CO2, we take as initial value 0.16, c o r r e s p o n d i n g to the point of the sample trajectory reached at t = c~=2x 10-4s. F o r the reduced models, the initial conditions are entirely specified by xs = CO2, the remaining variables xy result from algebraic equations. Since these equations are different for the two reduced models, the starting values for xy m a y be different even if the initial value of xs = CO2 is the same (see Tables 5 a n d 6). F o r Figs 5-7 (resp. Figs 8-10), the initial condition (xs, x I) of the complete model coincides with the starting values of the reduced slow model (resp. quasisteady-state model) given in T a b l e 5 (resp. Table 6).
4670
P. DUCHgNE and P. ROUCHON
CO 2 0.35
...._......-.~ °°o°°°°°°0°"~°~
r" 0.3
0
..'f
--
co
le e
.-- 0.25
¢-- 0.3 ._o ¢n ¢/)
E 0.2 .._a.._
0.2
0.4
0.6
x 10.3
~.8 ~.~
"
0.4
0.6
0.8
time (S)
x 1 .3
Fig. 8. Mass fraction of CO2 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
H20
~
0.2
x 10.3
Fig. 5. Mass fraction of CO2 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
9
,/,
0.15
0.8
time (S)
..... / ~.
0.25 0.2
0'1-0
C02
0.35
9 x 10.3
H20
8~'~
~
._o
~6
......... "...~ . " ' ~ " Y
tIE
5
86 °°° I
~ - - c o le e
°°°°°°
F : ." ..~ 5 "" ~''"
--co
le~
"..!.."U" "\ /
4
~
0
0.2
0.4
0.6
0.8
time (S)
0
Fig. 6. Mass fraction of H20 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
0.4
0.6
0.8
time (s)
x 10.3
Fig. 9. Mass fraction of HzO calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
0
0.012 0.01
0.2
x 10.3
0
0.012/
,.
=,
0.01 ~ . . r
• °'7," \
t-
_o 0.008
,°
~ 0.006
0.006
\\
0.006 0'} (n
t~
~ 0.004 I
~ 0.004
/
0.002 0
0.002 0.2
0.4
0.6
time (s)
0.8
0.2 x 10.3
0.4
0.6
time (S)
0.8 x 10 .3
Fig. 7. Mass fraction of O calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
Fig. 10. Mass fraction of O calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
These choices of initial conditions for the complete model a d m i t the following explanation: reduction is a restriction to a n attractive i n v a r i a n t manifold. So the c o m p a r i s o n between the reduced models a n d the
complete one is meaningful only for trajectories starting near this manifold. W h e n the initial c o n d i t i o n of the complete model coincides with that of the reduced slow one, the
Kinetic scheme reduction
4671
Table 5. Initial conditions for the reduced slow model (mass fractions, temperature) CO H2 02 H20
0.18 1.6x10 -s 0.10 8.4 x 10 -3
CO 2
OH HO2 H202
0.16 2.7x10 4 2.4x 10 -5 1.6x 10 -6
CH20 H O HCO
2.6 x 10 -s 3.3x10 5 1.9x 10 -5 3.7x l0 -7
N2 Temp.
0.55 1297K
Table 6. Initial conditions for the quasi-steady-state model (mass fractions, temperature) CO H2 02 H20
0.18 4,2×10 -5 0.10 7.5 X 10 -3
CO 2
OH HO 2 H202
0.16 4.5 x 10 4 3.0× 10 -5 3.3x 10 6
CH~O H O HCO
trajectories remain very close. This is no longer true for the quasi-steady-state model.
e
6. CONCLUDING REMARKS
E
The reduction method presented here can be applied to any dynamical system having two time-scales (a slow one and a fast asymptotically stable one). As described by L6vine and Rouchon (1993), this method can be adapted for deriving slow control models where time derivatives of the control appear in the slow equation even if they are not present in the modeling equations. Notice that the explicit value of ~: is not required for computing the reduced slow system (8). This tends to indicate that such approximation techniques can be extended to more general systems than two time-scale ones, and, in particular, to systems admitting a normally hyperbolic and attractive invariant manifold (Hirsch et al., 1977).
Acknowledgements--The results presented here have been motivated by active discussions with Jean L6vine, especially during the Geometrical Colloquium in Moscow, May 93. This work benefits also from several discussions with Alain Bamberger, Zakia Benjelloun, Pierre Galtier, Anne Jaecker, Willy Nastoll and Arnaud Trouv6 from I.F.P. NOTATION
A C E h k n t T v
V x
stoichiometric matrix correction matrix linear subspace of R" approximation of the equation of the slow manifold kinetic constant dimension of the system to be reduced time evolution time of the system vector field (second member of a system of ordinary differential equations) or reaction rates (vector) neighborhood of R" concentration of a species, or vector of R"
Greek letters ct fastest time-scale that must be taken into account
f~
1.2 x 1.0x 3.8× 9.2x
10 -7 10 -4 10 3 10 -7
N2 Temp.
0.55 1297K
small nonnegative real number (0 ~< ~<< 1) new variable depending linearly on the old one (possibly a vector) invariant slow manifold time-scale open bounded subset of R"
Subscripts/superscripts f fast s slow
REFERENCES
Arnold, V. I., 1980, Chapitres Suppl~mentaires de la Th~orie des Equations Differentielles Ordinaires. Mir, Moscou. Brenan, K. E., Campbell, S. L. and Petzold, L. R., 1989, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. North-Holland, Amsterdam. Carr, J., 1981, Application of Center Manifold Theory. Springer, Berlin. Ehrenstein, U., 1990, Techniques de r6duction de seh+mas cin6tiques. Application & la mod61isation cin6tique de la pyrolyse de l'hexane. Technical Report 38397, Institut Franqais du P6trole. Ehrenstein, U, and Nastoll, W., 1991, Simplification de m6canismes r6actionnels par la th6orie de la vari6t6 centrale, application fi la mod61isation cin6tique de la production des NOx. Technical Report 39121, Institut Fran¢ais du P6trole. Fenichel, N., 1979, Geometric singular perturbation theory for ordinary differential equations. J. Differential Equations 31, 53-98. Guckenheimer, J. and Holmes, P., 1983, Nonlinear Oscillation.s, Dynamical Systems and Bifurcations of Vector Fields. Springer, New York. Hirsch, M. W., Pugh, C. C. and Shub, M., 1977, lnvariant Manifolds. Lecture Notes in Mathematics, Springer, Berlin. Kuo, J. C. and Wei, J., 1969, A lumping analysis in monomolecular reaction systems. Analysis of the approximately lumpable system. Ind. Enyn9 Chem. Fundam. 8, 124 133. L6vine, J. and Rouchon, P., 1993, An invariant manifold approach for robust control design and applications. Proceedinys MTNS'93, Regensburg, Germany. Li, G. and Rabitz, H., 1989, A general analysis of exact lumping in chemical kinetics. Chem. Engno Sci. 44, 1413--1430. Li, G. and Rabitz, H., 1993, Determination of approximate lumping schemes by a singular perturbation method. J. Chem. Phys. 99, 3562-3574.
4672
P. DUCHENEand P. ROUCHON
Maas, U. and Pope, S. B., 1992, Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combust. Flame 88, 239-264. Tikhonov, A., Vasil'eva, A. and Sveshnikov, A., 1980, Differential Equations. Springer, New York. van Breusegem, V. and Bastin, G., 1991, Reduced order dynamical modelling of reaction systems: a singular perturbation approach. Proceedings 30th Conference on Decision and Control, Brighton, U.K., pp 1049-1054.
van Breusegem, V. and Bastin, G. 1993, A singular perturbation approach to the reduced order dynamical modeling of reaction systems. Technical Report 92.14, CESAME, University of Louvain-la-Neuve, Belgium. Wiggins, S., 1994, Normally Hyperbolic lnvariant Manifolds in Dynamical Systems. Applied Mathematical Sciences, Vol. 105. Springer, Berlin.
Chemical Engineering Science, Vol. 51, No. 20, pp. 4661-4672, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0009-2509/96 $15.00 + 0.00
S001D-2509(96)00310-7
KINETIC SCHEME REDUCTION VIA GEOMETRIC SINGULAR PERTURBATION TECHNIQUES P. DUCHENE* Institut Francais du Prtrole, BP 3, 69390 Vernaison Cedex, France and P. ROUCHON l~cole des Mines de Paris, Centre Automatique et Syst~mes, 60 Bd Saint-Michel, 75272 Paris, Crdex 06, France (Received 19 M a y 1995; accepted in revised form 3 April 1996)
Abstract--An effective reduction method based on geometric singular perturbation and center manifold techniques is proposed. This method eliminates the fast and stable dynamics and gives the equations describing the slow ones. It is coordinate-free and extends the well-known quasi-steady-state method classically used for kinetic scheme reduction. This yields directly the slow dynamics even if the differential equations are not in standard two time-scale form (Tikhonov form). Application to combustion kinetics including 13 species and 67 reactions is presented and simulations are given. Copyright © 1996 Elsevier Science Ltd Keywords: Model reduction, geometric singular perturbations, center manifold, combustion, kinetics.
1. I N T R O D U C T I O N
Let us begin with a simple reaction scheme. It involves three chemical species Xx, X2 and X3 and three elementary independent reactions: XI k ~ x z ,
X2 k2~,X1, X I + X 2
4- k z x 2 -
2 2 = klx I -- k2X 2
3(1
=
--
klX1 + k2x2 - g k x l x 2
(2) 0 = klXl - k2x2
~k.... ' X 2 + X 3
the xis are compositions of the Xis; kl, k2 and ek are the kinetic constants. The small positive parameter indicates here that the third reaction is slow with respect to the first two reactions. From this reaction scheme, we derive an ordinary differential equation system. It corresponds to the conservation of each species in a closed homogeneous reactor (constant volume, molar volume independent of compositions, constant temperature). The species balances are: "?(1 = - - k l X l
The quasi-steady-state assumption applied to x 2 leads to the following reduced system:
or
Y¢1 = - e ( k k l / k 2 ) ( X l ) 2
(3)
x2 = ( k l / k 2 ) x l .
(4)
with
This system does not yield a correct approximation of the slow dynamics for small e. The reason is the following. From eq. (1), we have d(xl + x2)/dt = - e k x l x 2 . Thus, if eq. (4) was valid, one would have (1 + k l / k 2 ) X l = - e ( k k J k z ) ( X l ) 2.
ekXlX 2
(1)
-~3 = g k X l X 2
where the dots above the xs indicate differentiation with respect to time. Since x3 follows from the chemical invariant X l + X 2 + X s = X ° + x ° + x ° (the superscript 0 indicates initial compositions), we focus on the first two equations of (1).
This stands in contradiction with eq. (3). Here, the quasi-steady-state method cannot be directly applied. A change of coordinates is required first. Instead of using (x~, x2), consider the coordinates (~ = X l + X2, X2).
The first two equations of (1) are transformed to = -- ~k(¢ -- X2)X 2
(6) 3~2 = k l ( ~ -
* Corresponding author.
(5)
x2) -
k2x2.
Notice that ~ corresponds to an approximate lumping of xx and x2 [for more about lumping, see e.g. Kuo 4661
4662
P. DUCH~NEand P. ROUCHON
and Wet (1969); Li and Rabitz (1989, 1993)]. The application of the quasi-steady-state method to x2 leads then to the correct approximation: = -- 8k(~ -- X2)X2,
0 = kl(~. -- X2) -- k2x2,
that reads, in the original coordinates (x~, x2), = - e~-2 (Xl)
with x2 = ~ x l .
(7)
The difference with eq. (2) lies only in the coefficient (1 + kUk2) multiplying the derivative of x~. The algebraic equation, giving x2 with respect to x~ remains unchanged. For this simple system, the 'good coordinates' (~, x2) are easy to find. Van Breusegem and Bastin (1991, 1993) explained how to find such a linear change of coordinates when the reactions can be separated into fast and slow ones. For more complex systems, such as the combustion system of Section 5, this is no longer the case: there is no clear decomposition into slow and fast reactions and thus the change of coordinates is probably nonlinear. The need for model reduction is motivated by two purposes: a qualitative analysis of the model and decreasing the simulation times, in particular when kinetics is coupled with transport. In this paper, we propose a reduction method that overcomes the difficulty of finding the coordinate change. The obtained reduced model [eqs (8) and (9) of Section 2] is directly given in the modelling variables. A deeper analysis yields the following explanations. The quasi-steady-state method is not coordinate-free: it can be applied only to systems written with special coordinates, called Tikhonov coordinates in the sequel (see Section 3.2). Our extension provides a coordinate-free reduction method. This means that no change of coordinates is required. This extension admits a clear and rigorous justification within the theory of perturbations of dynamical systems (Arnold, 1980; Guckenheimer and Holmes, 1983; Wiggins, 1994). More precisely, we use here the theory of singular perturbations (Tikhonov et al., 1980), its geometric setting due to Fenichel (1979) and a special version of the center manifold approximation lemma (Carr, 1981); see also Ehrenstein (1990) and Ehrenstein and Nastoll (1991) for application of center manifold techniques to reduction of chemical kinetic schemes, in a slightly different context. Another approach to model reduction is proposed by Maas and Pope (1992). When the system admits two time-scales, one can prove that this method is equivalent to ours.
The content is the following. In Section 2, we describe the practical aspect of our method. The mathematical framework is the object of Section 3: we recall the classical Tikhonov form of singularly perturbed systems; we explain the geometric definition of slow/fast systems (Fenichel, 1979); the system reduction via invariant manifold techniques is derived on the basis of an approximation lemma that can be found in Carr (1981). Section 4 is devoted to a special
class of chemical systems where the two time-scales are due to fast and slow reactions. In Section 5, the application of our reduction method to a combustion scheme considered by Maas and Pope (1992) is presented. For this system, the quasi-steady-state method yields a rather poor approximation. 2. T H E R E D U C T I O N
METHOD
For readers not interested in technical developments, we just present here the formulae that are required for applying our reduction method to a dynamical system of the form = v(x)
with x = (xl . . . . . x~), v = (Vl. . . . . v~)
n being the dimension of x. Assume that this system admits two time-scales: • a slow time-scale that can be described by n~ < n variables. • a fast and (hyperbolically) stable time-scale that can be described by n I = n - ns variables. The integer ns will be the number of differential equations of the reduced slow system. The integer nf corresponds to the number of derivatives that are set to 0 in the quasi-steady-state method, i.e. the number of algebraic equations that must be solved. As in the quasi-steady-state method, we decompose the state x~R" into two sets of components, x = (xs, xs.) of dimension n~ for xs and nI for xf. Similarly, the system is decomposed into two parts
Yc, = v,(xs, xs),
~s = vs(x~, xs)
with v = (v,, vs). This decomposition is up to the choice of the user. We give below some guidelines to make this choice. In particular, the eigenvalue structure of the Jacobian matrix (cnvi/c3xs) provides an estimation of nf. Then, the reduced system is as follows:
2~ = C(x,, ho)" vs(x~, ho) 0 = Vs(X~, ho) xf=ho-
(8)
~[(c?vI~ 2 + ~l)f'~[)s]-xovf)7}
tLk~.U
'~, a,
ex.d"
Us(Xs,
ho)
where the matrix C(xs, ho) is given by
C = I . , - c,~xI L k ~ x #
+
~x,
axsJ ~
(9)
with 1,, the identity matrix of order n~. All the partial derivatives are computed at x = (xs, ho). Some algebraic manipulations provide another equivalent expression for C that might by preferable in practice:
C=[lo " + ax~,t~s)
77]x,j
-' "
Notice that formulae (8), applied to eq. (1) with x, = (x~, x3) and x I = x2, lead to the correct approximation, the correction factor (1 + kUk2) of eq. (7)
4663
Kinetic scheme reduction corresponding then to the matrix C - 1 (up to terms of order 1 in e). In fact, the only difference between eq. (8) and the equation resulting from the quasi-steady-state method lies in the correction matrix C. For systems with v, = el,, such as system (6) with x~ = ~, the matrix C is close to 1,, up to terms of order 1 in e. In this case, our method coincides, up to terms of order 2 in e, with the quasi-steady-state one. In order to avoid the computation and the inversion of the Jacobian matrices contained in eqs (8) and (9), it may be useful to approximate eq. (8) via the following differential-algebraic system: dx, dt
v~(x~, xs)
0 = vf(x.~, ho)
(10)
0 = vf(x, + 6vs(x~, x f), ho + 6vf(x,, x f))
where 6 is a very small parameter and (x~,ho, xy) are unknown. Some elementary computations show that if 6 ( ~ 0) is at least of order 2 in e then the variable x I in eq. (10) satisfies eq. (14) for p = 1. Such differentialalgebraic systems of index 1 can be easily solved numerically via Gear-like schemes and standard codes of public domain, such as DASSL or L S O D E / L S O D I [see e.g. Brenan et al. (1989)]. As far as we know, formulae (8) and its approximation (10) are new and have never been introduced elsewhere in such a simple formulation. Such equations are, in general, more complicated than the original ones. This difficulty can be overcome, if the formulae are approximated by off-line interpolations using only xs as independent variable. Then, the simulation just consists of integrating an explicit differential system ~s = P(x~), where function P results from such an interpolation. Note that, for the computation of the reduced slow dynamics, it is not required to exhibit e as a function of the model parameters. We just need to know that the system admits two time-scales, that the dimension of the fast part is n I and that it is asymptotically stable. This point and the simplicity of the formulae constitute the major interest of the method proposed here. In practice, we have to define n I (or equivalently n, = n - nl). Since our reduction method amounts to eliminating fast and stable time-scales, we may be guided by two criteria: 1. What is the fastest time-scale that cannot be eliminated? 2. The time-scales of the reduced model (slow timescales) must be significantly slower than those that have been eliminated. The answer to the first question is given by the following considerations: • The kinetic model may be part of a bigger model including, for example, transport. In this last case, one must keep in the reduced model the
time-scales that are of the same order of magnitude as those of transport phenomena (and of course the slower ones). We denote by ~ this typical time constant. • The time of evolution of the system, T, may be a priori known. In this case, one eliminates those time-scales that are significantly faster than T. The characteristic times of the system, (Ti) i 1 . . . . . may be estimated with the help of the eigenvalues of the Jacob]an Ov/Ox, (2i)i= 1...... : zi ~ 1/[Re()~i)[.
To satisfy the second criterion, one computes the eigenvalues along different trajectories of the system. Then one tries to separate the characteristic terms into two groups, (zl,r2 . . . . . r,s) and (r,~+l . . . . . T,), such that "gl ~ "C2 ~
"'" ~ ~'n! <<0~ ~ "Cnt+l ~
"'" ~ "(n.
The eigenvalues relative to the first group have a large, negative real part: they correspond to fast and stable phenomena. The eigenvalues of the second group correspond to the slow part of the system and thus are close to zero. After having determined ny, one has to parametrize the reduced model, i.e. choose a decomposition of x into (xs, xy). In view of the algebraic equations of (8), it seems reasonable to choose the x I variables among those that are at or near quasi-steady-state. A species is in quasi-steady-state when its creation rate is (almost) equal to its destruction rate. This may be evaluated by the ratio Icreation - destruction] creation + destruction for each species i. This ratio is based on the classical decomposition of vi v,(x) = ~,? (x) - v~- (x)
where the creation v [ ( x ) corresponds to reactions producing i, and the destruction v~(x) to reactions consuming i. The species for which this ratio is small (compared to 1) may be considered as being in quasisteady-state. This ratio can be numerically calculated along the trajectories of the system. We choose the x I variables among those for which this ratio is small. Other procedures for defining ns and (xs, x r) can be used. We just have sketched here a simple one that gives, for the combustion system of Section 5, an efficient reduced model of dimension one. Let us now explain how we have obtained eqs (8) and (9). 3. G E O M E T R I C SINGULAR PERTURBATIONS
3.1. Model reduction The following question underlies this paper: what is model reduction? A possible and reasonable answer to this question for dynamical systems of the form .;: = v(x), x ~ ~" is displayed in Fig. 1. Reduction can
4664
P. DUCHI~NEand P. ROUCHON
Fig. 1. Trajectories of ~ = v(x) are caught by an attractive invariant manifold Z.
The terminology 'singular perturbations' comes from the fact that the small parameter e multiplies the highest derivative (here d x l / d t ). More details on this classical standpoint can be found, e.g. in Tikhonov et al. (1980). In the sequel we always consider the time-scale t and approximations of trajectories for t~ [0, I/z] and with e > 0 but close to 0. Assume that the fast part is hyperbolically stable, i.e. the sub-system x.r = vf(x~,xr,e) with xs fixed admits (locally) an equilibrium with characteristic exponents (eigenvalues of Ovl/t3x I at this equilibrium) having a strictly negative real part. Then the slow approximation is obtained by the quasi-steady-state method (Tikhonov et al., 1980):
2~ = evs(xs, x f, ~) be defined via an attractive invariant manifold T,. A manifold E is invariant with respect to the vector field v, if v is tangent to E: if a trajectory starts on Y., it remains on Z for all time. E is called (locally) attractive if any trajectory, t ~ x(t) of ~ = v(x) (starting near Y.), tends to E as t tends to + ~ . In this case, the reduced system corresponds to the restriction of the vector field v to E. This restriction is meaningful since Z is invariant. It seems then natural to approximate trajectories of the complete system ~ = v(x) by trajectories on E. In fact, such an approximation is proved to be valid when, roughly speaking, the dynamics transverse to ~z (the dynamics that are neglected) are faster than the dynamics on E [see Arnold (1980) and Hirsch et al. (1977) for more details on normally hyperbolic invariant manifolds]. Reduction means then restriction of the dynamics to an invariant attractive manifold E. In practice the equations defining E are not known. Only the vector field, v, is explicitly available. The main difficulty is thus to obtain the equations of E or, at least, good approximations of them. 3.2. Tikhonov normal form of slow~fast dynamical system Here we restrict the study to a special class of dynamical systems, depending on a small parameter e and having an attractive invariant manifold E~. Such systems are characterized by dynamics transverse to E~ which is much faster than the dynamics on E~. This corresponds to systems having two time-scales: a fast and asymptotically stable one and a slow one (stable or unstable). The first kind of such systems is singularly perturbed systems of the so-called Tikhonov form,
dx~ d----t= e v~(x~,xy, e),
dxf dt - vI(x~'xI'e)
(11)
where 0 < 5<< 1. Very often, such systems are written with the time-scale z = st: dx~ d--~- = vs(x~, x I , e),
dx I e-d~t = vAx~, x l , 5).
0 = Vs(Xs, x s, e). The algebraic equations, vy = 0, correspond here to an approximation up to terms of order 1 in e, of E, equations. These coordinates (x,, xl) where the quasi-steady-state method applies, and where the vector field v is quasi-vertical (see Fig. 2) are clearly very specific. 3.3. Geometric settin9 Since we are interested in developing a reduction method that avoids changing coordinates, we need a coordinate-free point of view. The geometric definition of singularly perturbed systems due to Fenichel (1979) is as follows. Consider the dynamical system :~=v(x,e),
x~R", 0~
(12)
This system is said to have two time-scales, a fast and asymptotically stable one and a slow one, if, and only if, the following two assumptions are satisfied (Fenichel, 1979): (A1) for e = 0, eq. (12) admits an equilibrium manifold of dimension n~, 0 < n~ < n, denoted by ~0. (A2) for all Xo e Zo, the Jacobian matrix, c3v/Ox[txo.o~ admits n I = n - n~ eigenvalues with a strictly negative real part (the eigenvalues are counted with their multiplicities). This definition is illustrated in Fig. 3. Assumption (A1) implies that the velocity v(x, ~) is large everywhere except for x in a neighborhood of Y'o where v is small and of order unity in e. (A1) and (A2) imply that, for Xo~Y~o, the kernel, Uo(Xo), of the linear operator t~v/Oxllxo, O~ coincides with the tangent space of Eo at Xo. The linear space E~(xo) corresponding to the eigenvalues with real negative part satisfies
E~o(Xo) @ Eo(Xo) = R".
4665
Kinetic scheme reduction
a partition of the state x into two groups of components, x = (x,, x:), with dim(x~) = n~ and dim(x:) = n:, such that the projection of E~ on the x:coordinates is a local diffeomorphism around Xo. The equations defining E, around Xo admits the following approximation:
v( )
:(Xs, x:,
x , = ho(x,, e) + 0(~)
=
x d = ho(x~, e,)
(13)
x. tkVx,) + Fig. 2. In Tikhonov normal form, the trajectories of eq. (11) admit a vertical part before being caught by the attractive manifold E~.
axsj
&sl,.,.ho.,,=),,:)J
• v.,(x,, ho(x~, g), g) + 0(8 2)
where ho is the solution of vs(x~, ho(x~, ~:),~:) = 0. The trajectories of the perturbed system are captured by a trapping region around Eo and entering with a direction nearly parallel to E~o(Xo). Fenichel (1979, part of theorem 9.1) proves the following result• It asserts, for e small enough, the existence of a slow invariant attractive manifold Z~ for the perturbed system (12).
Theorem 1 (Feniehel, 1979): Consider eq. (12) satisfying (A1) and (A2). Then, for every open and bounded subset ~o of Y~o, there exists an open neighborhood Vo of f~o in ~", such that, for e positive and small enough, the perturbed system (12) admits an attractive invariant sub-manifold Z~ contained in Vo and close to Eo. We are interested in approximations, up to terms of order 1 in e, of slow trajectories for t e [0, l/e]. Thus we need an approximation up to terms of order 2 for the slow dynamics: errors like e 2 integrated over t e [0, I/e] will produce for t = I/e, distortions of magnitude less than or equal to ~: (e2 x(1/,:)= e). This means that an approximation, up to terms of order 2 in e, of £, equations is needed. 3.4. Approximation of Z~
Proposition 1. Consider eq. (12) satisfying (AI) and (A2). Take Xo ~ Zo. Then locally around Xo, there exists
v(x, e)
"" (""
!
In some generic sense, that can be easily defined, almost every partition of x in (Xs, x:) is admissible for eq. (13). But in practice, the coordinates x~ must be chosen such that the projection of Z~ onto the x:space is as well conditioned as possible•
Proof of Proposition 1. By the approximation lemma of invariant manifold, and its version for two timescale systems (Carr, 1981, theorem 5, p. 32), the equation x : = hp(xs, e) is an approximation up to terms of order p + 1 in e. of X~ equations, if hp satisfies the following equation stating the invariance of E,::
~h_ Vz(x..h,(x..e).e)=-" v.(x..h,(x..e). ~) + ~Xs (x. ~)
O(e'+b
(14) The zeroth order term, ho(x,,e), is given by the quasi-steady-state approximation
v:(xs, ho(x~, e), e) = 0 because v,(xs, ho(xs, e), e) = O(e). The first-order term, hi is written as hi = ho + epl where pl is a smooth function of xs that will be
x
v(x, O)
JX
ose to
1 _-:...y.o
Fig. 3. The geometric definition of singularly perturbed systems due to Fenichel (1979).
4666
P.
DUCH&NEaInd
wo
+ EPl)
4 + &PI, 4 = 8%
x u,(x,, h&s, E)+
&PI,
ROUCHON
Notice that approximation (8) applied to a system already in Tikhonov normal form yields an approximation of the dynamics on C,, up to terms of order 3 in E [use eq. (13) and replace o, in eq. (13) by EV, for systems like eq. (1 l)].
adjusted in order to verify (14) with p = 1, i.e. V/(X,, hhs,
P.
(I,, S)
4 + W).
(15) 4. CHEMICAL SYSTEMS WITH SLOW AND FAST
First, we develop the left side of the above equality. Noticing that u~(x,, ho(x,, a), E) = 0, we obtain
U&G,ho(&, E)+ &PI,E)=
au,
Eaxf
Pl + W).
(&.h&E),
@
Developing the rightmost term, we get r,(X,, ho@,, 8) + &PI, E) = Vs(xs, ho(xs, 81,E) +c!S ax, It
PI + ok+). (x,,h&.
REACTIONS
In this section we consider the case of kinetic systems where the reactions can be split into one group of slow reactions and one group of fast reactions. The reduction of such systems has been studied by Van Breusegem and Bastin (1991, 1993). We show here that the reduced models deduced from their calculations and the reduced model obtained from eq. (8) coincide, up to terms of order 2 in E. Consider the reaction system
E).c)
1 =
is a first-order term in &(0(s)) because %(&, ho@,, s), e) = O(E).
We now develop the right side of equality (15). d(ho + &Pl) r&G, ho(x,, c) + EPl, e) ax, (X,X E)
zdV(X,
E)
(16)
where x is the vector of concentrations, A is the stoichiometric matrix, v(x) is the vector of the chemical rates and E is a small number (0 < EK 1). Assume that v can be split into v = (EC,, a,), where ~6, (resp. 5,) is the vector of the chemical rates of the slow (resp. fast) kinetics. System (16) can be rewritten as k = A, E&(X) + A/d,(X).
=- aho
(17)
U&G, ho(xs, s) + &PI, ~1+ W2)
ax, (I,,E) =- aho a xs kc)
Let A,, be a square sub-matrix of A/ such that its rank is equal to the rank of A,. Let
u,(xs,
ho(xs> ~1,~1 A, = Pl + W2).
and Equating both sides of eq. (15) and using the relation
aho
-lavf av,_
(>
dx, = - ax, ax,
Then eq. (16) reads
we obtain J, = A,, at%,
XJ) + A,, cf(x,, x/).
The change of coordinates (~8, x/) H (5 = x, - &(AJixs>
x,)
hence eq. (13). The matrix leads to a quasi-vertical vector field. It yields the following Tikhonov form:
au,2 auf au,
ax, +ax’ax (-> s
f
g = (A,, - &(A,,)-‘A&G
is invertible: this results from (Al) (A2) and from some n elementary calculations. The reduced slow model (8) comes from eq. (13) and u,(x,,
ho +
EPI, ~1= t’stxs, ho, 6) + av, axr
+ O(E2).
&PI @..h,,e)
1, = Afs.$, + A,,Cf. Assuming that eigenvalues of A,,-diT,/dx, have strictly negative real parts, the quasi-steady-state method can be applied and leads to the following slow system: i = CA,, - &(A,,)-‘A~,l~~s
0 = A,,
EC3
+
A,,Cf
4667
Kinetic scheme reduction
The above formulae give the initial conditions (~o, y~) of the slow sub-system from the ones of the complete system (x°, x~).
Pulling back into the original coordinates ( x . x:) yields
l(~vs~-'~Vs 1 X, = II,, - a~s(af:)- \ ~ x : J -~XsJ 5. CASE-STUDY
× [Ass -- Asf(Aff)- XAfs]r.fs 0 = A:s~fs(xs, x:) + A:: fy(xs, xf). This approximation differs from eq. (8) only via terms of order 2 in e: at order 1 in e both approximations coincide. The main feature of eq. (8) in this case relies on the following point: one does not need to know precisely the decomposition between slow and fast kinetics. This can be an advantage if one is only interested in the slow dynamics. Notice that, even if the two timescale structure is not due to fast and slow kinetics, formulae (8) are still valid. However, for such special slow/fast kinetic schemes, one can, as displayed on Fig. 4, describe how and where trajectories starting far from the slow attractive manifold Z~ are caught by X~. The quantity = x., - Asf(Aff)- ~x: plays an important role here: it is a quasi-invariant during the fast transient, before arriving near Y~. As displayed on Fig. 4, the trajectory of eq. (18) with initial condition (x°, x}), becomes close, after a time interval of length of order 1, to the trajectory of the slow system here above, starting at the same time 0 with initial conditions (fro .~) solution of .~0 _ A s f ( a f f )
1 ~ = No __ A s f ( A f f )
~:Ix °, x~,
lx~,
~.) = 0.
In this section, we apply the above reduction method to the CO/H2/air combustion system also considered by Mass and Pope (1992). Thirteen chemical species are present: N2, CO, H2, 02, H20, C02, OH, H, O, HO2, CHO, H202, CH20. Sixty-seven reactions are taken into account. Their list is given in Table 1. Notice that N2 is inert. The system will be supposed closed, isobaric and adiabatic. It is modeled by a system of differential eel ,ations whose dimension is nine (the system is ist,_aric, adiabatic, and four atomic invariants are present: C, H, N and O). The determination of n: and of the partition (x~, x:) is based on a sample trajectory starting at t = 0 from the initial conditions given in Table 2. This corresponds to the beginning of the explosive phase that lasts about T = 10 -2 s, the time-scale of interest. The eigenvalues at different times along this trajectory are given in Table 3 (they have been sorted in increasing order). All eigenvalues, except one at the beginning of the trajectory, are real and negative. For the reduction purpose, one may notice that the eighth and ninth eigenvalues are well separated for t >/2 × 10 -4 s): the ratio of the latter to the former (when both are negative) ranges from 2 to more than 200 for t/> 2 × 10 4. This indicates that we should be able to build a one-dimensional reduced model. Thus n: = 8 and n~ = 1. For t close to 2 × 10 -4 s the separation is not very large (ratio of 2). This may result in poor accuracy of the reduced model. If so, the
i x' / . ' / / x o .................
......
.¢. ..........
............
.,
0
b
,,A--
-0 Xs
. ~
0 Xs
Fig. 4. F o r s y s t e m (18), fast t r a n s i e n t s a r e a l o n g x~ - A ~ : ( A : : ) - ~ x ~ = c o n s t a n t .
4668
P. DUCHI~NEand P. ROUCHON Table
1. Reaction
list
[mass
action
law
and
kinetics
constant
given
k = oT~exp(-Ea/RT)]
ko O2 + H --, O H + O O H + O ~ 02 + H H2 + O ~ O H + H O H + H ---, H2 + O H 2 + OH ~ H20 + H H 2 0 + H ~ H2 + O H OH + OH ~ H20 + O H20 + O ~ OH + OH H+H+M~Hz+M H2+M~H+H+M H+OH+M ~ H20+M HzO+M ~ H+OH+M O+O+M-*O2+M O2+M~O+O+M H + 02 + M ~ HO2 + M HO2 + M -.-* H + 02 + M HO2 + H -," OH + O H O H + O H ~ HO2 + H HO2 + H -~ H2 + 02 H2 + 0 2 -'¢ HO2 + H HO2 + H --, H 2 0 + H H 2 0 + O ~ HO2 + H HO 2 + O --* O H + 02 O H + O2 --' HO2 + O H O 2 + O H --' H 2 0 + 02 H 2 0 + 02 "-* HO2 + O H HO2 + HO2 ~ H202 + 0 2 O H + O H + M ~ H202 + i H2OE+M ~ OH+OH+M H 2 0 2 + H -'* H 2 + HO2 H2 + HO2 --* H202 + H H2O 2 + H -~ H 2 0 + OH H 2 0 + O H ~ H20 2 + H H 2 0 2 + O ~ O H + HO 2 O H + HO 2 ~ S 2 0 2 + O H 2 0 2 + O H ~ H 2 0 + HO2 H 2 0 + H 0 2 -'* H202 + O H CO + O H ~ CO2 + H C O 2 2r- H ~ CO + OH CO + HO2 ~ CO2 + O H CO2 + O H ~ CO + HO2 CO + O + M ~ C O 2 "~- M CO2 + M ~ CO + O + M CO + 02 --* CO2 + O CO 2 + O --* CO + 0 2 HCO + M -~ CO + H + M CO + H + M ~ H C O + M HCO + H ~ CO + H2 CO + H2 --' H C O + H HCO + O ~ CO + O H CO + OH ---, H C O + O HCO + O ~ CO2 + H CO2 + H --* H C O + O H C O + O H ~ CO + H 2 0 CO + H 2 0 ~ HCO + O H H C O + 02 ~ CO + HO2 CO + HO2 --' H C O + O2 C H 2 0 + M -~ H C O + H + M H C O + H + M -* C H 2 0 + M C H 2 0 + H --* HCO + H2 H C O + H2 --* C H 2 0 + H C H 2 0 + O ""¢ HCO + O H HCO + OH ~ C H 2 0 + O C H 2 0 + OH ~ HCO + H 2 0 HCO + H 2 0 -'-' C H 2 0 + O H C H 2 0 + H O 2 --* HCO + H 2 0 2 H C O + H 2 0 2 --~ C H 2 0 + HO2
2.00 x 10 TM 1.47x 1013 5.06 x 104 2.24 × 104 1.00 x 10 s 4.46× l0 s 1.50x 109 1.51 × 101° 1.80x 10 is 6.98 x l018 2.20x 1022 3.80×1023 2.90 x 1017 6.78 x 1018 2.30x 1018 2.66 x l018 1.50 × 101. 1.63 X 1013 2.50x 1013 8.39 x 1013 3.00x 1013 3.29 × 1013 1.80 × 1013 2.67 × 1013 6.00 × 1013 8.97 x 101'* 2.50× 1011 3.25 X 1022 2.11 x 1024 1.701 1012 9.35 x 1011 1.00x 1013 2.66 × 1012 2.80 × 1013 6.80x 1012 5.40x 1012 1.32× 1013 4.40× l06 6.12 x 10 s 1.50x 101'* 2.27 x 10 is 7.10× 1013 1.69 × 1016 2.50 x 1012 2.55 x 1013 7.10 x 1014 1.07 x 1015 2.00 x 1014 1.17 x 1015 3.00 x 1013 7.72 x 1013 3.00 x 1013 1.07 x 1016 1.00 x 1014 2.601 1015 3.00x 1012 5.21 × 1012 1.40x 1017 2.62 x 1015 2.50 × 1013 1.82 x 1012 3.50x 1013 1.12 × 1012 3.00x 1013 9.71 x 1012 1.00 x 1012 1.32x 1011
fl 0.00 0.00 2.67 2.67 1.60 1.60 1.14 1.14 -1.00 -1.00 -2.00 -2.00 -1.00 -1.00 -0.80 -0.80 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -2.00 -2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.50 1.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
E~ 16,791 502 6282 4395 3296 18,451 96 17,101 0 104,137 0 119,280 0 118,563 0 49,250 1003 37,738 693 55,603 1720 55,460 -406 52,618 0 70,030 --1242 0 49,393 3750 21,783 3583 73,469 6401 22,547 1003 34,155 -740 22,475 23,574 83,524 --4538 120,952 47,769 54,696 16,791 2054 0 89,424 0 87,537 0 110,753 0 104,543 0 34,513 76,431 -13,566 3989 18,152 3487 15,764 1194 30,477 8001 4132
by
Kinetic scheme reduction
4669
Table 2. Initial conditions of the complete system for the sample trajectory (mass fractions, temperature) CO Hz Oz H20
0.21 4.9 x 10 4 0.12 4.5x 10 -3
COz OH HO2 H202
0.11 3.7 x 10 -s 5.9 x 10 -5 7.5 x 10 --6
CH20 H O HCO
l a x 10 -n
N 2
6.0 X 10 -6
Temp.
0.55 1076K
1.4x 10 -4 2.3x 10 -7
Table 3. Eigenvalues of the CO/HE/air system at different times along the sample trajectory t=0s -
7.0
t=10-gs
-- 4.4 x 107
-
- 4 . 1 xl06 - 1.7x106 - 6 . 8 x 105
- 1.2x107 - 3 , 0 x 106 - 2 . 1 xl06
- 8.4x 106 - 8.0x 106 - 3.0x 106
-
3.9 x 104
-
1.4 x 1 0 5
-
2.6 x 105
-
1.6
x
-
1.9
x
-
1.2 x 1 0 4
-
4.7
-
1.6 x 1 0 5
-
5.5
× 10 5
--
1.0
X 10 6
8.8 × 1 0 3 - 1.2x103 1.3 x 104 --
-
x
10 4
1.8 x 1 0 4 --6.3x 103 1.5 x 104 --
1.9 x 1 0 7
t=5xl0-4s
- 1,2xl06 - 1,2x106 - 1.5x105
10 6
-
t=3xl0-4s
- 9 . 1 x 105 - 4 . 2 x l05 - 1.2x105
x
9.6 x 106
t=2xl0-4s
5.4 x 1 0 4 - 2 . 7 x 104 1.1 x 104 -
10 6
2.2 x 105 - 2.6x 104 - 5.6 x 1 0 3 --
6.1 x 1 0 7
10 6
- 4.3 x 10s - 6.0x 104 - 4.7 x 1 0 3
t=10-3s
t=10-2s
7 . 3 x 107 - 1.9x107 - 4 . 6 x 106 - 2 . 6 x 106
9.8 x 107 - 5.4x 107 - 2.4x 106 -2.3x106 8.8 × 1 0 5 5 . 4 x 105 - 2.6 x 105 - 6.8 x 10'* - 2.5 x 1 0 2
-
-
1.2
× 10 6
--
1.1
x
10 6
- 4.3 x 105 - 7 . 4 x 104 - - 1.6 x 1 0 3
Table 4. Creation/destruction ratios of the CO/H2/air system at different times along the sample trajectory t=0s HCO HO 2 OH H2Oz H CHzO H20 O 02 H2 CO CO2
0.00 0.00 0.00 0.51 0.01 0.08 0.93 0.08 0.58 0.74 0.91 0.99
t=10-4s 0.00 0.00 0.00 0.00 0.02 0.00 0.83 0.10 0.46 0.71 0.91 0.99
first-order reduced model of P r o p o s i t i o n 1 can be i m p r o v e d by adding higher-order terms (approximation l e m m a of i n v a r i a n t manifold for two time-scale systems used in the p r o o f of P r o p o s i t i o n 1). The ratios
t=2xl0-4s 0.00 0.00 0.00 0.01 0.01 0.01 0.03 0.05 0.19 0.19 0.88 0.94
(3)
t=3xl0-4s 0.00 0.00 0.00 0.00 0.00 0.00 0.06 0.01 0.04 0.00 0.48 0.50
t=5xl0-4s 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.03
the one-dimensional quasi-steady-state model
Yc~ = vs(x~, vy)
(19)
0 = Vs(X.~,x j-)
(20)
with n.~ = 1 a n d x~ = COz (dotted line). Icreation - destruction l creation + destruction at different times along this trajectory are given in Table 4 (they have been sorted in increasing order). The ratio for N2 is not given since this species is neither created n o r destroyed. Since the largest ratio c o r r e s p o n d s to CO2 (or CO), we can set x~ to the mass fraction of CO2 (or C O as well). The simulations of Figs 5-10 c o m p a r e three models: (1) the complete model of dimension nine (solid line), (2) the one-dimensional reduced slow model resulting from formulae (8), with ns = 1 a n d xs = CO2 (dashed line),
In these simulations, the starting value ofxs = C O 2 remains u n c h a n g e d between the three models. F o r xs = CO2, we take as initial value 0.16, c o r r e s p o n d i n g to the point of the sample trajectory reached at t = c~=2x 10-4s. F o r the reduced models, the initial conditions are entirely specified by xs = CO2, the remaining variables xy result from algebraic equations. Since these equations are different for the two reduced models, the starting values for xy m a y be different even if the initial value of xs = CO2 is the same (see Tables 5 a n d 6). F o r Figs 5-7 (resp. Figs 8-10), the initial condition (xs, x I) of the complete model coincides with the starting values of the reduced slow model (resp. quasisteady-state model) given in T a b l e 5 (resp. Table 6).
4670
P. DUCHgNE and P. ROUCHON
CO 2 0.35
...._......-.~ °°o°°°°°°0°"~°~
r" 0.3
0
..'f
--
co
le e
.-- 0.25
¢-- 0.3 ._o ¢n ¢/)
E 0.2 .._a.._
0.2
0.4
0.6
x 10.3
~.8 ~.~
"
0.4
0.6
0.8
time (S)
x 1 .3
Fig. 8. Mass fraction of CO2 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
H20
~
0.2
x 10.3
Fig. 5. Mass fraction of CO2 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
9
,/,
0.15
0.8
time (S)
..... / ~.
0.25 0.2
0'1-0
C02
0.35
9 x 10.3
H20
8~'~
~
._o
~6
......... "...~ . " ' ~ " Y
tIE
5
86 °°° I
~ - - c o le e
°°°°°°
F : ." ..~ 5 "" ~''"
--co
le~
"..!.."U" "\ /
4
~
0
0.2
0.4
0.6
0.8
time (S)
0
Fig. 6. Mass fraction of H20 calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
0.4
0.6
0.8
time (s)
x 10.3
Fig. 9. Mass fraction of HzO calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
0
0.012 0.01
0.2
x 10.3
0
0.012/
,.
=,
0.01 ~ . . r
• °'7," \
t-
_o 0.008
,°
~ 0.006
0.006
\\
0.006 0'} (n
t~
~ 0.004 I
~ 0.004
/
0.002 0
0.002 0.2
0.4
0.6
time (s)
0.8
0.2 x 10.3
0.4
0.6
time (S)
0.8 x 10 .3
Fig. 7. Mass fraction of O calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the reduced model).
Fig. 10. Mass fraction of O calculated by the complete model and the one-dimensional reduced models (the initial condition of the complete model is the same as that of the quasi-steady-state model).
These choices of initial conditions for the complete model a d m i t the following explanation: reduction is a restriction to a n attractive i n v a r i a n t manifold. So the c o m p a r i s o n between the reduced models a n d the
complete one is meaningful only for trajectories starting near this manifold. W h e n the initial c o n d i t i o n of the complete model coincides with that of the reduced slow one, the
Kinetic scheme reduction
4671
Table 5. Initial conditions for the reduced slow model (mass fractions, temperature) CO H2 02 H20
0.18 1.6x10 -s 0.10 8.4 x 10 -3
CO 2
OH HO2 H202
0.16 2.7x10 4 2.4x 10 -5 1.6x 10 -6
CH20 H O HCO
2.6 x 10 -s 3.3x10 5 1.9x 10 -5 3.7x l0 -7
N2 Temp.
0.55 1297K
Table 6. Initial conditions for the quasi-steady-state model (mass fractions, temperature) CO H2 02 H20
0.18 4,2×10 -5 0.10 7.5 X 10 -3
CO 2
OH HO 2 H202
0.16 4.5 x 10 4 3.0× 10 -5 3.3x 10 6
CH~O H O HCO
trajectories remain very close. This is no longer true for the quasi-steady-state model.
e
6. CONCLUDING REMARKS
E
The reduction method presented here can be applied to any dynamical system having two time-scales (a slow one and a fast asymptotically stable one). As described by L6vine and Rouchon (1993), this method can be adapted for deriving slow control models where time derivatives of the control appear in the slow equation even if they are not present in the modeling equations. Notice that the explicit value of ~: is not required for computing the reduced slow system (8). This tends to indicate that such approximation techniques can be extended to more general systems than two time-scale ones, and, in particular, to systems admitting a normally hyperbolic and attractive invariant manifold (Hirsch et al., 1977).
Acknowledgements--The results presented here have been motivated by active discussions with Jean L6vine, especially during the Geometrical Colloquium in Moscow, May 93. This work benefits also from several discussions with Alain Bamberger, Zakia Benjelloun, Pierre Galtier, Anne Jaecker, Willy Nastoll and Arnaud Trouv6 from I.F.P. NOTATION
A C E h k n t T v
V x
stoichiometric matrix correction matrix linear subspace of R" approximation of the equation of the slow manifold kinetic constant dimension of the system to be reduced time evolution time of the system vector field (second member of a system of ordinary differential equations) or reaction rates (vector) neighborhood of R" concentration of a species, or vector of R"
Greek letters ct fastest time-scale that must be taken into account
f~
1.2 x 1.0x 3.8× 9.2x
10 -7 10 -4 10 3 10 -7
N2 Temp.
0.55 1297K
small nonnegative real number (0 ~< ~<< 1) new variable depending linearly on the old one (possibly a vector) invariant slow manifold time-scale open bounded subset of R"
Subscripts/superscripts f fast s slow
REFERENCES
Arnold, V. I., 1980, Chapitres Suppl~mentaires de la Th~orie des Equations Differentielles Ordinaires. Mir, Moscou. Brenan, K. E., Campbell, S. L. and Petzold, L. R., 1989, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. North-Holland, Amsterdam. Carr, J., 1981, Application of Center Manifold Theory. Springer, Berlin. Ehrenstein, U., 1990, Techniques de r6duction de seh+mas cin6tiques. Application & la mod61isation cin6tique de la pyrolyse de l'hexane. Technical Report 38397, Institut Franqais du P6trole. Ehrenstein, U, and Nastoll, W., 1991, Simplification de m6canismes r6actionnels par la th6orie de la vari6t6 centrale, application fi la mod61isation cin6tique de la production des NOx. Technical Report 39121, Institut Fran¢ais du P6trole. Fenichel, N., 1979, Geometric singular perturbation theory for ordinary differential equations. J. Differential Equations 31, 53-98. Guckenheimer, J. and Holmes, P., 1983, Nonlinear Oscillation.s, Dynamical Systems and Bifurcations of Vector Fields. Springer, New York. Hirsch, M. W., Pugh, C. C. and Shub, M., 1977, lnvariant Manifolds. Lecture Notes in Mathematics, Springer, Berlin. Kuo, J. C. and Wei, J., 1969, A lumping analysis in monomolecular reaction systems. Analysis of the approximately lumpable system. Ind. Enyn9 Chem. Fundam. 8, 124 133. L6vine, J. and Rouchon, P., 1993, An invariant manifold approach for robust control design and applications. Proceedinys MTNS'93, Regensburg, Germany. Li, G. and Rabitz, H., 1989, A general analysis of exact lumping in chemical kinetics. Chem. Engno Sci. 44, 1413--1430. Li, G. and Rabitz, H., 1993, Determination of approximate lumping schemes by a singular perturbation method. J. Chem. Phys. 99, 3562-3574.
4672
P. DUCHENEand P. ROUCHON
Maas, U. and Pope, S. B., 1992, Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combust. Flame 88, 239-264. Tikhonov, A., Vasil'eva, A. and Sveshnikov, A., 1980, Differential Equations. Springer, New York. van Breusegem, V. and Bastin, G., 1991, Reduced order dynamical modelling of reaction systems: a singular perturbation approach. Proceedings 30th Conference on Decision and Control, Brighton, U.K., pp 1049-1054.
van Breusegem, V. and Bastin, G. 1993, A singular perturbation approach to the reduced order dynamical modeling of reaction systems. Technical Report 92.14, CESAME, University of Louvain-la-Neuve, Belgium. Wiggins, S., 1994, Normally Hyperbolic lnvariant Manifolds in Dynamical Systems. Applied Mathematical Sciences, Vol. 105. Springer, Berlin.