Constructing Slow Invariant Manifolds for Reacting Systems with - - PowerPoint PPT Presentation

Constructing Slow Invariant Manifolds for Reacting Systems with Detailed Kinetics

• A. N. Al-Khateeb, J. M. Powers, S. Paolucci, J. D. Mengers

Department of Aerospace and Mechanical Engineering

• A. J. Sommese, J. D. Hauenstein, and J. Diller

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

61st APS Division of Fluid Dynamics Meeting

San Antonio, Texas 24 November 2008

Introduction

Motivation and background

• Detailed kinetics are essential for accurate modeling of real

reactive systems.

• Reactive systems contain many scales and subsequently severe

stiffness arises.

• 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 savings.

Slow Invariant Manifold (SIM)

• The composition phase space for closed spatially homogeneous

reactive system:

dz dt = f (z) , z ∈ RN−L−C.

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

Method of Construction

• For isothermal reactive systems, reactions speeds depend on

combinations of polynomials of species concentrations.

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

{ze ∈ CN−L−C |f (ze) = 0}; we focus on real equilibria.

• 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 0-D physical critical point (a sink) to two saddles.

• These saddles are identified by their special dynamical character:

their eigenvalue spectrum contains only one unstable direction.

Sketch of SIM Construction

R1 R2

S I M

R3

Projective Space for Equilibria at Infinity

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

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

• This transformation maps equilibria located at infinity into a finite

domain.

• To address the time singularity, we add the following transforma-

tion

dt dτ = (Zk)n−1 ,

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

Computational strategy

• We use the Bertini software (based on a homotopy continuation

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.

Zel’dovich Mechanism for NO Formation

• The mechanism (see Baulch et al., 2005) consists of J

= 2 reversible bimolecular reactions involving N = 5 species {NO, N, O, N2, O2} and L = 2 elements {N, O}.

In addition, since the total number of moles is constant, C = 1. Subsequently, z ∈ R2.

• Spatially homogenous with isothermal and isochoric conditions,

T = 4000 K, p0 = 1.65 atm.

• We find three 0-D finite equilibria (R1: source, R2: saddle, R3:

sink, physical) and three 0-D infinite equilibria (I1: saddle-node,

I2: source, I3: source)

The system’s 1-D SIM

×10

• 2

R SIM ×10

• 2

N mol/g

[ ]

NO mol/g

[ ]

2

R1 R3 I1 2.5

• 0.5

1 2 −2 −1 1 2 1.5 0.5

Equilibrium Thermodynamics and SIM

Within the physically accessible domain,

σ = − 1 T (∇G · f) ≥ 0,

at equilibrium

Hσ = − 2 T (HG · Jf) .

5 ×10

• 5

3 5 ×10

• 3

SIM I1 R2 R3 NO mol/g

[ ]

N mol/g

[ ]

4 3 2 1 6 7 8 4 2 1 3 SIM R3 ×10

• 7
• 10

N N [mol/g]

_

e

NO NO [mol/g]

_

e

R2 I1 5 10

• 5
• 10

×10 4 2 1

• 1
• 2
• 3
• 4
• The major/minor axes are aligned

with the Hessian eigenvectors.

• Eigenvectors of equilibrium thermo-

dynamic potentials do not coincide with system’s SIM, even at the physical equilibrium point!

Hydrogen-Air System

• The mechanism (Miller et al., 1982) consists of J = 19 reversible

reactions involving N = 9 species, L = 3 elements, and C =

0, so that z ∈ R6.

• Closed and spatially homogenous system of 2H2 + (O2 +

3.76N2) with isothermal and isochoric conditions at T = 1500 K,

and p0 = 107 dyne/cm2.

• The system has 284 finite (90 0-D real) and 42 infinite (18 0-D

real) equilibria.

• Only 14 critical points have an eigenvalue spectrum that contains
• nly one unstable direction.
• There is a unique physical equilibrium, R19.

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 mol/g

[ ]

O mol/g

[ ]

m

• l

/ g

[ ]

2

OH H2

Summary

• Constructing the actual SIM is computationally efficient and algo-

rithmically easy, thus there is no need to identify it only approxi- mately.

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

in the construction of the SIM.

• Irreversibility production rate and equilibrium thermodynamic po-

tentials do not provide information on the dynamics towards physical equilibrium.

The 2nd International Workshop on Model Reduction in Reacting Flow

March 30—April 2, 2009

HenryCurran,NationalUniversityofIreland YannisKevrekidis,PrincetonUniversity MarcMassot,CNRS—UniversiteClaudeBernard LindaPetzold,UniversityofCalifornia—SantaBarbara JamesRawlings,UniversityofWisconsin JamesRobinson,UniversityofWarwick

Center for Applied Mathematics

ACCEPTEDINVITEDSPEAKERS

ScientificCommittee: M.Giona,UniversityofRome“LaSapienza” D.Goussis,UniversityofAthens H.Najm,SandiaNationalLaboratories,Livermore S.Paolucci,UniversityofNotreDame J.M.Powers,UniversityofNotreDame A.J.Sommese,UniversityofNotreDame M.Valorani,UniversityofRome“LaSapienza” For up-to-date information please go to the following web site: http://cam.nd.edu/upcoming-conferences/spring2009 LocalOrganizingCommittee: S.Paolucci,UniversityofNotreDame J.M.Powers,UniversityofNotreDame A.J.Sommese,UniversityofNotreDame

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 C = 0, so that z ∈ R3.

• Spatially homogenous with isothermal and isobaric conditions

with T = 3000 K, p0 = 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

1-D SIM vs. 2-D ICE manifold

z mol/g

3

[ ]

z mol/g

1

[ ]

z mol/g

2

[ ]

1 2 3 4 2 4 4 2 6 ×10

• 3

SIM

1.5 2 2.5 0.2 3.4 3.8 4.2

S I M

×10

• 3

×10

• 3

R6 R1 R7

Outline

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

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 The construction of 1-D Slow Invariant Manifolds (SIMs) for dynam- ical system arising from modeling unsteady spatially homogenous closed reactive systems.

Partial review of manifold construction in reactive systems

• ILDM, CSP

, and ICE-PIC are approximations of the reaction slow invariant manifold.

• MEPT and similar methods are based on minimizing a thermody-

namics potential function.

• Iterative methods require “reasonable” initial conditions.
• 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.

• An invariant manifold is defined as an open subset S ⊂ RN−L−C

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.

Reactive system evolution

NO N 10

• 6

10

• 8

10

• 4

10

• 10

t s 10

• 4

10

• 3

10

• 2

z mol/g 10

• 5

[ ]

i

[ ]

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

Finite equilibria

dz1 dt = 2.51 × 102 + 1.16 × 107z2 + 6.99 × 108z2

2

−9.98 × 104z1 − 3.22 × 109z2z1, dz2 dt = 2.51 × 102 − 1.17 × 107z2 − 6.98 × 108z2

2

+8.47 × 104z1 − 1.84 × 109z2z1,

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

R1 ≡ (ze

1, ze 2)

= ` −1.78 × 10−5, −1.67 × 10−2´ , (λ1, λ2) = (4.18 × 107, 2.35 × 107) source, R2 ≡ (ze

1, ze 2)

= ` −4.20 × 10−3, −2.66 × 10−5´ , (λ1, λ2) = (−4.64 × 106, 7.11 × 105) saddle, R3 ≡ (ze

1, ze 2)

= ` 3.05 × 10−3, 2.94 × 10−5´ , (λ1, λ2) = ` −1.73 × 107, −1.91 × 105´ sink.

R3 is the physical equilibrium. Stiffness = |λ1/λ2| = 90.5

Infinite equilibria

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

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

1

B B @ Z−1

1

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

1, Ze 2)

= (0, 0) , (λ1, λ2) = ` −1.53 × 1013, 0 ´ saddle − node, I2 ≡ (Ze

1, Ze 2)

= (0, 1.01) , (λ1, λ2) = ` 2.12 × 1013, 9.36 × 1012´ source, I3 ≡ (Ze

1, Ze 2)

= (0, 2.60) , (λ1, λ2) = ` 3.04 × 1013, 2.41 × 1013´ source.

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 and 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 S I M