Calculation of Slow Invariant Manifolds for Reactive Systems Ashraf - - PowerPoint PPT Presentation

Calculation of Slow Invariant Manifolds for Reactive Systems

Ashraf N. Al-Khateeb Joseph M. Powers Samuel Paolucci

Department of Aerospace and Mechanical Engineering

Andrew J. Sommese Jeffery A. Diller

Department of Mathematics University of Notre Dame, Notre Dame, Indiana

47th AIAA Aerospace Science Meeting

Orlando, Florida 8 January 2009

Outline

• Introduction
• Slow Invariant Manifold (SIM)
• Method of Construction
• Illustration Using Model Problem
• Application to Hydrogen-Air Reactive System
• Summary

Introduction

Motivation and background

• Detailed kinetics are essential for accurate modeling of real

systems.

• Reactive flow systems admit multi-scale solutions.
• Severe stiffness arises in detailed gas-phase kinetics modeling.
• Computational cost for reactive flow simulations increases with

the spatio-temporal scales’ range, the number of species, and the number of reactions.

• Manifold methods provide a potential for computational saving.

Partial review of manifold construction in reactive systems

• ILDM, CSP

, and ICE-PIC are approximations of the system’s slow invariant manifold.

• MEPT, RCCE, and similar methods are based on minimizing a

thermodynamic potential function.

• Iterative methods may not converge.
• Davis and Skodje, 1999, present a technique to construct the 1-D

SIM based on global phase analysis.

• Creta et al.

and Giona et al., 2006, extend the technique to slightly higher dimensional reactive systems.

Long-term objective Create an efficient algorithm that reduces the computational cost for simulating reactive flows based on a reduction in the stiffness and dimension of the composition phase space. Immediate objective Construct 1-D SIMs for dynamical system arising from modeling unsteady spatially homogenous closed reactive systems.

Slow Invariant Manifold (SIM)

• The composition phase space for closed spatially homogeneous

reactive system:

dz dt = f (z) , z ∈ R3.

Phase space trajectory Phase space trajectory 1-D manifold 2-D manifold 0-D manifold (i.e. equilibrium point) Fast modes Slow modes

• An invariant manifold is defined as a subset S ⊂ RN−L−Q if

for any solution z(t), z(t0) ∈ S, implies that for any tf > t0,

z(t) ∈ S for all t ∈ [t0, tf].

• Not all invariant manifolds are attracting.
• SIMs describe the asymptotic structure of the invariant attracting

trajectories.

• Attractiveness of a SIM increases as the system’s stiffness in-

creases.

• On a SIM, only slow modes are active.
• SIMs can be constructed by identifying all critical points, finite and

infinite, and connecting relevant ones via heteroclinic orbits.

Mathematical Model

For a mixture of mass m confined in volume V containing N species composed of L elements that undergo J reversible reactions,

dni dt = V ˙ ωi, i = 1, . . . , N,

where,

˙ ωi =

J

X

j=1

νijkj N Y

i=1

“ni V ”ν′

ij − 1

Kc

j N

Y

i=1

“ni V ”ν′′

ij

! , i = 1, . . . , N, kj = Aj T βj exp „−Ej ℜT « , j = 1, . . . , J, Kc

j

= „ po ℜT «PN

i=1 νij

exp − PN

i=1 ¯

µo

i νij

ℜT ! , j = 1, . . . , J.

System reduction

• In chemical reactions, the total number of moles of each element

is conserved,

N

• i=1

φlin∗

i = N

• i=1

φlini, l = 1, . . . , L.

• Additional Q constraints can arise in special cases.
• The reactive system is recast as an autonomous dynamical

system,

dzi dt = fi (z1, . . . , zN−L−Q) , i = 1, . . . , N − L − Q,

where,

z = L (n)

• L :
• RN → RN−L−Q

.

Method of Construction

• For isothermal reactive systems, reaction speeds depend on

combinations of polynomials of z.

• The set of equilibria of the full reaction network is complex:

{ze ∈ CN−L−Q |f (ze) = 0}.

• The set consists of several different dimensional components and

contains finite and infinite equilibria.

• A 1-D SIM has a maximum of two branches that connect the

unique physical critical point (a sink) to two equilibria.

• These equilibria are identified by their special dynamical char-

acter: their eigenvalue spectrum typically contains only one unstable direction.

Sketch of SIM construction

R1 R2

SIM

R3

Projective space

• One-to-one mapping of the composition space, RN−L−Q →

RN−L−Q, Zk = 1 zk , k ∈ {1, . . . , N − L − Q}, Zi = zi zk , i = k, i = 1, . . . , N − L − Q.

• This maps equilibria located at infinity into a finite domain.
• To deal with the time singularity, we add the transformation

dt dτ = (Zk)d−1 ,

where d is the highest polynomial degree of f(z).

Computational strategy

• We use the Bertini software (based on a homotopy continu-

ation numerical technique) to compute the system’s equilibria up to any desired accuracy.

• Thermodynamic data is obtained from Chemkin-II.
• The SIM heteroclinic orbits are obtained by numerical integration
• f the species evolution equations using a computationally inex-

pensive scheme.

• Computation time is typically less than 1 minute on a 2.16 GHz

Mac Pro machine.

Simple Hydrogen-Oxygen Mechanism

• The kinetic model is adopted from Michael, 1992, Prog. Energy
• Combust. Sci. 18(4), p. 327.
• The mechanism consists of J = 8 bimolecular elementary

reactions involving N = 6 species {H, H2, O, O2, OH, H2O} and L = 2 elements {H, O}. In addition, since the total number

• f moles is constant, Q = 1. Subsequently, z ∈ R3.
• The system is spatially homogenous with isothermal and iso-

choric conditions, T = 1200 K, V = 103 cm3.

• Selected species are i = {1, 2, 3} = {H2, O, O2}.
• Initial number of moles of all species are n∗

i = 10−3 mol.

Reactive system evolution

10

• 8

10

• 6

10

• 4

5 10 15 10

• 10

×

mol/g

i

[ ]

t s

[ ]

10

• 3

20 25 O2 H2 H H O

2

O OH

n /m

Dynamical system

dz1 dt = 3.45 × 104 − 1.68 × 1011z1 − 3.47 × 1016z2

1

−1.35 × 1010z2 + 6.27 × 1016z1z2 + 1.40 × 1010z3 +1.11 × 1017z1z3 − 1.35 × 1016z2z3 − 2.04 × 1016z2

3,

dz2 dt = 7.69 × 105 + 2.66 × 1011z1 + 2.25 × 1016z2

1

−1.29 × 1012z2 − 2.47 × 1017z1z2 +3.91 × 1017z2

2 − 8.66 × 1011z3 − 1.51 × 1017z1z3

7.49 × 1017z2z3 + 2.46 × 1017z2

3,

dz3 dt = 6.84 × 1011z2 + 1.37 × 1017z1z2 − 2.74 × 1017z2

2

−4.10 × 1017z2z3 − 2.24 × 1015z3 ` 10−6 − z1 + z3 ´ ,

                                           ≡ f(z).

Finite equilibria

R1 ≡ (ze

1, ze 2, ze 3)

= ` −5.84 × 10−2, 6.85 × 10−4, −3.52 × 10−4´ mol/g, (λ1, λ2, λ3) = (5.93 × 106 ± i5.10 × 105, −1.18 × 106) 1/s, R2 ≡ (ze

1, ze 2, ze 3)

= ` 4.65 × 10−2, 0, 3.49 × 10−2´ mol/g, (λ1, λ2, λ3) = (−1.01 × 107, −3.35 × 106, 7.93 × 105) 1/s, R3 ≡ (ze

1, ze 2, ze 3)

= ` 3.73 × 10−3, 6.32 × 10−3, 1.61 × 10−2´ mol/g, (λ1, λ2, λ3) = (−1.02 × 107, −1.23 × 106, −4.30 × 105) 1/s, R4 ≡ (ze

1, ze 2, ze 3)

= ` 6.33 × 10−3, −1.86 × 10−3, 2.49 × 10−2´ mol/g, (λ1, λ2, λ3) = (6.88 × 106, 3.51 × 106, 1.57 × 106) 1/s, R5 ≡ (ze

1, ze 2, ze 3)

= ` 1.28 × 10−3, −5.98 × 10−2, 6.00 × 10−2´ mol/g, (λ1, λ2, λ3) = (5.65 × 107, 3.56 × 106, −1.06 × 104) 1/s, R6 ≡ (ze

1, ze 2, ze 3)

= ` 1.43 × 10−3, −7.58 × 10−2, 7.08 × 10−2´ mol/g, (λ1, λ2, λ3) = (7.19 × 107, 4.47 × 106, 1.05 × 104) 1/s.

Infinite equilibria

• Employ the projective space mapping with d = 2 and k = 2:

d dτ B B B B B @ t Z1 Z2 Z3 1 C C C C C A = Z2

2 ·

B B B B B @ Z−1

2

f1 (Z1, Z2, Z3) − Z1 f2 (Z1, Z2, Z3) −Z2 f2 (Z1, Z2, Z3) f3 (Z1, Z2, Z3) − Z3 f2 (Z1, Z2, Z3) 1 C C C C C A ≡ F(Z), I1 ≡ (Ze

1, Ze 2, Ze 3)

= (−9.77, 0, −4.59) , (λ1, λ2, λ3) = (−5.74 × 1012 ± i7.83 × 1012, 6.10 × 1012), I2 ≡ (Ze

1, Ze 2, Ze 3)

= (0.60, 0, −0.48) , (λ1, λ2, λ3) = (−1.19 × 1013, 7.35 × 1011, 6.32 × 1011), I3 ≡ (Ze

1, Ze 2, Ze 3)

= (−0.01, 0, −0.67) , (λ1, λ2, λ3) = (−1.12 × 1013, −6.50 × 1011, 7.62 × 109).

The system’s 1-D SIM

0.5 1 2 2.5 3 3.5 4 4.5 0.01 0.02 0.03 0.01 0.02 0.03

S I M SIM ×10

• 2

R2 R3 I3 z mol/g

3

[ ]

z mol/g

2

[ ]

z mol/g

1

[ ]

1.5

Detailed Hydrogen-Air Mechanism

• A kinetic model is adopted from Miller et al., 1982, Proc. Com-
• bust. Ins. 19, p. 181.
• The mechanism consists of J = 19 reversible reactions involving

N = 9 species, L = 3 elements, and Q = 0, so that z ∈ R6.

• Closed and spatially homogenous system with isothermal and

isochoric conditions at T = 1500 K, and p∗ = 107 dyne/cm2.

• Stoichiometric mixture 2H2 + (O2 + 3.76N2).
• Selected species are

i = {1, 2, 3, 4, 5, 6} = {H2, O2, H, O, OH, H2O}.

Reactive system evolution

10

• 9

10

• 5

10

• 1

z mol/g 10

• 13

i

[ ]

10

• 17

10

• 5

10

• 8

10

• 2

10

• 11

t s

[ ]

10

1

H2 O2 H O

2

OH H O

System’s equilibria

• The system has 284 finite and 42 infinite equilibria.
• The set of finite equilibria contains 90 real and 186 complex 0-D,
• ne 1-D, one 2-D, and six 3-D equilibria.
• The set of infinite equilibria contains 18 real and 18 complex 0-D,

and six 1-D equilibria.

• Only 14 critical points have an eigenvalue spectrum that contains
• nly one unstable direction.
• Inside the physical domain there is a unique equilibrium:

R19 = ` 1.98 × 10−6, 9.00 × 10−7, 1.72 × 10−9, 2.67 × 10−10, 3.66 × 10−7, 1.44 × 10−2´ mol/g.

3-D projection of the system’s SIM

2 4 6 8 −1 1 2 3 4 5 −2 −1.5 −1 −0.5 0.5 1 SIM ×10

• 5

×10

• 5

×10

• 5

R19 R74 R79 z mol/g

5

[ ]

z mol/g

2

[ ]

z mol/g

1

[ ]

Summary and Conclusions

• Once the difficult task of identifying all equilbria is complete,

constructing the actual SIM is computationally efficient and al- gorithmically easy; thus, there is no need to identify it only approximately.

• Identifying all critical points, finite and infinite, plays a major role

in the construction of the SIM.

• The construction procedure can be systematically extended to

construct higher-dimensional SIMs.

Simple Reactive System

A + A ⇋ B kf = 1, kb = 10−5. B ⇋ C kf = 10, kb = 10−5.

• A reactive system adopted from D. Lebiedz, 2004, J. Chem. Phys.

120 (15), p. 6890.

• Model consists of J = 2 reversible reactions involving N = 3

species {cA, cB, cC}

• Conservation of mass, cA + cB + cC = 1, so that z ∈ R2.
• Major species are i = {1, 2} = {A, B},

The system’s global phase space

SIM R2 R1 Z2 Z1 I2 I1

The projective space.

R1 R2 SIM I

1

I

2

I

2

I

1

z

1

z

2 2

z

2 2 1

+ + z

1

z

1 2

z

2 2 1

+ +

Projection from Poincar´ e’s sphere.

The 1-D SIM vs. MEPT

0.1 0.2 0.5 1.0

z

2

z

1

0.5 1.0 0.02 0.04 0.06 0.08

z

2

R1

z

1

0.1

• 2
• 1

1 2 3 8 9 10 11 12 ×10

• 4

×10

• 4

z

1

z

2

R1 R2 SIM

Idealized Hydrogen-Oxygen

• Kinetic model adopted from Ren et al.a
• Model consists of J = 6 reversible reactions involving N =

6 species {H2, O, H2O, H, OH, N2} and L = 3 elements {H, O, N}, with Q = 0, so that z ∈ R3.

• Spatially homogenous with isothermal and isobaric conditions

with T = 3000 K, po = 1 atm.

• Major species are i = {1, 2, 3} = {H2, O, H2O},
• Initial conditions satisfying the element conservation constraints

are identical to those presented by Ren et al.

• aZ. Ren, S. Pope, A. Vladimirsky, J. Guckenheimer, 2006, J. Chem. Phys. 124, 114111.

The system’s 1-D SIM

−12 −10 −8 −6 −4 −2 2 4 6 0 2 4 6 8 2 4 R1 R6 ×10

• 3

×10

• 3

×10

• 3

z mol/g

3

[ ]

z mol/g

1

[ ]

z mol/g

2

[ ]

R7 SIM

The system’s 1-D SIM 7 −1 1 2 3 4 5 −1 1 2 3 ×10

• 3

6 5 4 1 2 3 4 5 SIM ×10

• 3

×10

• 3

z mol/g

3

[ ]

z mol/g

1

[ ]

z mol/g

2

[ ]

R1 R6 R7