Supercomputer Research in Chemistry and Chemical Engineering

Chapter 21 Extinction of Counterflow Premixed Laminar Flames 2

M. D.Smooke1and V. Giovangigli 1Department of Mechanical Engineering, Yale University, New Haven, CT 06520 Laboratoire de Mécanique Théorique, Université Paris VI, 4 Place Jussieu, 75230 Paris Cedex 05, France and Laboratoire d'Energétique, Ecole Centrale des Arts et Manufactures, 92290 Chatenay-Malabry, France Downloaded by RUTGERS UNIV on November 30, 2016 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch021

2

Problems in combustion and heat and mass transfer often depend upon one or more physical/chemical parameters. In many cases the combustion scientist is interested in knowing how the solution will behave if one or more of these parameters is varied. For some parameter regimes the governing equations can produce multiple solutions and the branches of the solution curve are linked via singular points. It is at these singular points, however, that the system exhibits special behavior. To be able to predict the solution structure in the neighborhood of these points, we employ a phase-space, pseudo arclength, continuation method that utilizes Newton-like iterations and adaptive gridding techniques. We apply the method in the solution of counterflowpremixed laminar flames. The modeling of steady-state problems in combustion and heat and mass transfer can often be reduced to the solution of a system of ordinary or partial differential equations. In many of these systems the governing equations are highly nonlinear and one must employ numerical methods to obtain approximate solutions. The solutions of these problems can also depend upon one or more physical/chemical parameters. For example, the parameters may include the strain rate or the equivalence ratio in a counterflow premixed laminarflame(1-2). In some cases the combustion scientist is interested in knowing how the system will behave if one or more of these parameters is varied. This information can be obtained by applying afirst-ordersensitivity analysis to the physical system (3). In other cases, the researcher may want to know how the system actually behaves as the parameters are adjusted. As an example, in the counterflow premixed laminar flame problem, a solution could be obtained for a specified value of the strain

0097-6156/87/0353-0404$06.00/0 © 1987 American Chemical Society

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

21.

SMOOKE & GIOVANGIGLI

Counterflow

Premixed

Laminar

Fiâmes

405

rate; this s o l u t i o n c o u l d t h e n be used as a s t a r t i n g estimate for a predictorcorrector c o n t i n u a t i o n process i n the same parameter.

I n m a n y applications

this procedure works well. However, w h e n there are m u l t i p l e solutions a n d the s o l u t i o n curve, w h e n p l o t t e d versus the parameter t h a t is v a r i e d , t u r n s back o n itself by passing t h r o u g h a singular p o i n t , this procedure encounters difficulty. It is at these singular points, however, t h a t the system often experiences special behavior. F o r e x a m p l e , i n the s o l u t i o n of counterflow p r e m i x e d l a m i n a r

flames,

as the s t r a i n rate is increased past a c r i t i c a l value, the flame w i l l e x t i n g u i s h . S i m i l a r l y , for a given s t r a i n rate, as the equivalence r a t i o is either increased or Downloaded by RUTGERS UNIV on November 30, 2016 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch021

decreased past a c r i t i c a l value, the flame w i l l also e x t i n g u i s h . Hence, accurate knowledge of these points provides i n f o r m a t i o n a b o u t the

flame's

flammability

limits. T o solve problems of p h y s i c a l interest n u m e r i c a l l y one often discretizes the s p a t i a l operators of the governing o r d i n a r y or p a r t i a l differential equations w i t h the result t h a t the s o l u t i o n of the o r i g i n a l continuous p r o b l e m is reduced to finding a n a p p r o x i m a t i o n to this s o l u t i o n at each point of a discrete mesh. T h e discrete nonlinear equations c a n be solved by N e w t o n ' s m e t h o d , or a variant thereof. I n a d d i t i o n , i f the s o l u t i o n contains regions of h i g h s p a t i a l activity, t h e n it is essential for reasons of b o t h efficiency a n d accuracy t h a t adaptive spat i a l grids be employed. However, i n the n e i g h b o r h o o d of singular p o i n t s , i.e., t u r n i n g points, where the derivative of the s o l u t i o n curve w i t h respect to the p a r a m e t e r possesses a v e r t i c a l tangent, N e w t o n ' s m e t h o d encounters difficulties. T o be able to follow the s o l u t i o n i n the n e i g h b o r h o o d of t u r n i n g points, we m o d ify the basic s o l u t i o n a l g o r i t h m a n d employ a p a t h t r a c i n g c o n t i n u a t i o n m e t h o d i n w h i c h points o n the arc of solutions are c o m p u t e d . I n this paper we a p p l y a powerfull c o m b i n a t i o n of c o n t i n u a t i o n techniques, Newton-like iterations a n d adaptive g r i d d i n g i n the study of the e x t i n c t i o n of strained p r e m i x e d l a m i n a r flames. I n p a r t i c u l a r , we investigate the e x t i n c t i o n of hydrogen-air flames by a c o m b i n a t i o n of s t r a i n rate a n d equivalence r a t i o v a r i a t i o n . T h e governing conservation equations are m o d e l e d w i t h c o m p l e x t r a n s p o r t a n d detailed c h e m i c a l kinetics.

Problem Formulation S t r a i n e d p r e m i x e d l a m i n a r flames have played a n i m p o r t a n t role i n recent theories of t u r b u l e n t p r e m i x e d c o m b u s t i o n . T h e reacting surface i n such models c a n be v i e w e d as b e i n g composed of a n u m b e r o f l a m i n a r flamelets (see e.g., L i b b y a n d B r a y , (4), L i b b y a n d W i l l i a m s , (5) a n d B r a y , L i b b y a n d M o s s , (6). T h e flamelets are often m o d e l e d by considering counterflowing streams of reactants or counterflowing streams of reactants a n d products. F r o m a n e x p e r i m e n t a l v i e w p o i n t , these flames c a n be p r o d u c e d w h e n a single reactant s t r e a m i m pinges o n an a d i a b a t i c w a l l or w h e n two counterflowing reactant streams (or a reactant s t r e a m a n d a p r o d u c t stream) emerge from two counterflowing c o a x i a l jets. I n the n e i g h b o r h o o d of the stagnation p o i n t p r o d u c e d by these flows, a chemically reacting b o u n d a r y layer is established. T h e governing conservation Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

SUPERCOMPUTER RESEARCH

406

equations c a n be reduced to a set of nonlinear two-point b o u n d a r y value p r o b lems along the s t a g n a t i o n p o i n t streamline. W h e n there is a single reactant jet, only one r e a c t i o n zone is p r o d u c e d . However, i f there are two reactant streams a n d each has the same exit velocity a n d equivalence r a t i o , t h e n a double flame is p r o d u c e d w i t h a plane of s y m m e t r y t h r o u g h the stagnation p o i n t a n d p a r a l l e l to the two jets (see F i g u r e 1). In this paper we consider d o u b l y p r e m i x e d l a m i n a r flames p r o d u c e d by two counterflowing c o a x i a l jets. B y a p p l y i n g a p p r o p r i a t e b o u n d a r y conditions at the plane of s y m m e t r y , the m o d e l we consider is, i n p r i n c i p l e , equivalent to t h a t of

Downloaded by RUTGERS UNIV on November 30, 2016 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch021

a single reactant s t r e a m i m p i n g i n g o n an adiabatic w a l l w i t h s l i p . O u r m o d e l for counterflowing p r e m i x e d l a m i n a r flames assumes a l a m i n a r , s t a g n a t i o n p o i n t flow.

We m o d e l the s y s t e m by e m p l o y i n g a b o u n d a r y layer a p p r o x i m a t i o n .

T h e governing equations for mass, m o m e n t u m , c h e m i c a l species a n d energy i n rectangular coordinates c a n be w r i t t e n i n the f o r m (2) «

«

+

dx pu

du

du +

/

ά η

Ά

9 +

dT

ρ

ν

dT

Ά

dp

0 d

m

dx ρ

= dy

dy

dx

d

ky

k

du \ \

= 0,

k

A

f,dT\

μ

dy \

+ ^-(pY V )-w W k

( Ι

(2.1) Λ λ

=

(2.2)

0

ay j k = l,2,...,K

(2.3)

A . _, .

dT

n

/ n

..

T h e s y s t e m is closed w i t h the ideal gas law,

P

_ Ρ™ ~ RT

(2.5)

In these equations χ a n d y denote independent s p a t i a l coordinates; T , the t e m perature;

1

Y * , the mass fraction of the k^

species; p, the pressure; u a n d ν

the t a n g e n t i a l a n d the transverse components of the velocity, respectively; p, 1

the mass density; W * , the m o l e c u l a r weight of the k^

species; W, the m e a n

m o l e c u l a r weight of the m i x t u r e ; R, the u n i v e r s a l gas constant; λ , the ther m a l c o n d u c t i v i t y of the m i x t u r e ; c , the constant pressure heat capacity of the p

1

m i x t u r e ; c , the constant pressure heat capacity of the k^ P A

1

m o l a r rate of p r o d u c t i o n of the k^ enthalpy of the

species; ti;*, the

species per u n i t volume; /ι*, the specific

species; μ the viscosity of the m i x t u r e a n d Vjb , the diffusion y

1

velocity of the k^

species i n the y d i r e c t i o n . T h e free s t r e a m tangential a n d

transverse velocities at the edge of the b o u n d a r y layer are given by u = ax a n d t

v = —ay, respectively, where a is the s t r a i n rate. T h e s t r a i n rate is a measure e

of the "stretch" i n the flame due to the i m p o s e d flow. T h e f o r m of the c h e m i c a l p r o d u c t i o n rates a n d the diffusion velocities can be found i n (7-8).



Counterflow

Premixed

Laminar

Flames


y

Flame

^Symmetry Plane

0

χ

F i g u r e 1. Schematic of the stagnation point flow configuration.


407


408 U p o n i n t r o d u c i n g the n o t a t i o n r

(

2

6

)


= h · αχ V = ρν (2.7) where / ' is related t o the derivative of a modified s t r e a m function (see e.g., D i x o n - L e w i s et a l . , (9)), the b o u n d a r y layer equations c a n b e transformed i n t o a s y s t e m of o r d i n a r y differential equations v a l i d along the stagnation-point stream line χ = 0. F o r a s y s t e m i n rectangular coordinates, we have

+ apf' = 0

^

d

= 0,

-^{pYkV ,)-V^+w W k

k

k

( >

k = 1,2,...,Κ

v

(2.10)

Wthk=0

i ( Ι) - > % λ

M

+«->

Y

Arclength

Continuation

Procedures e n a b l i n g the c a l c u l a t i o n of bifurcation a n d l i m i t points for systems of n o n l i n e a r equations have been discussed, for example, by K e l l e r (13) H e i n e m a n n et a l . (14-15) a n d C h a n (16). I n p a r t i c u l a r , i n the w o r k of H e i n e m a n et a i . , a version of K e l l e r ' s pseudo-arclength c o n t i n u a t i o n m e t h o d was used to calculate the m u l t i p l e steady-states of a m o d e l one-step, n o n a d i a b a t i c , p r e m i x e d l a m i n a r flame ( H e i n e m a n n et a l . , (14)) a n d a p r e m i x e d , n o n a d i a b a t i c , hydrogen-air s y s t e m ( H e i n e m a n n et a l . , (15)). In o u r c o m p u t a t i o n a l m o d e l the s t r a i n rate a n d the equivalence r a t i o are the n a t u r a l b i f u r c a t i o n parameters. If we denote either of these parameters by a , t h e n the s y s t e m of equations i n (3.2) c a n be w r i t t e n i n the f o r m F{U,a)=0

(4.1)

where specific reference t o the p a r a m e t r i c dependence of F o n α has been made. A s the flame nears e x t i n c t i o n (either from a n increase i n the s t r a i n rate or f r o m a change i n the equivalence r a t i o ) , the m a x i m u m value of the t e m p e r a t u r e de creases.

A t the e x t i n c t i o n p o i n t the J a c o b i a n of the system is singular.

To

alleviate the c o m p u t a t i o n a l difficulties, a modified f o r m of the governing equa tions is solved (13). W e introduce the parameter α as a new dependent variable. T

T h e vector of dependent variables (C7, a)

c a n now be considered functions of a

new independent parameter θ. If we define T

Z(s) = (U(s), a(s))

(4.2)

t h e n the new p r o b l e m we want to solve is given by G{Z,s)=0

(4.3)

where G(Z,s)

=

F(U(s),a(s)) N{U(s),a(s),s)

• F(Z(s)) N(Z(s),s)

•


(4.4)

21.


Counterflow

Premixed

Laminar

Flames 4 1 1

a n d where Ν is a n a r b i t r a r y n o r m a l i z a t i o n . T h e n o r m a l i z a t i o n is chosen such t h a t s approximates the arclength of the s o l u t i o n b r a n c h i n the space (C7, a ) . T h e J a c o b i a n o f the new s y s t e m c a n be w r i t t e n i n the f o r m J(Z,s)

= " Fu(U{8) a{*)) N {U{8)M*U)

= G (Z,s)

9

z

V

where Fu,F ,Nu Q

and N

a

F (U{s) a(*)) N (U(s),a(s),s) a

(4.5)

9

a

denote the appropriate p a r t i a l derivatives. It c a n be

s h o w n t h a t at a s i m p l e t u r n i n g p o i n t , even t h o u g h Fu is s i n g u l a r J is n o t . W e p o i n t o u t t h a t the J a c o b i a n of the system i n (3.2) is block t r i d i a g o n a l . However,


after i n t r o d u c t i o n o f α as a dependent variable along w i t h the e x t r a n o r m a l i z a t i o n c o n d i t i o n , the b l o c k t r i d i a g o n a l structure o f the J a c o b i a n is destroyed. F o r the n o r m a l i z a t i o n s considered b y H e i n e m a n n et a l . (14-15) (see also K e l l e r , (13)) this is the case. A l t h o u g h s o l u t i o n of the s y s t e m o f linear equations cor r e s p o n d i n g t o (4.3) c a n proceed b y methods discussed i n (16), we w o u l d like to m a i n t a i n the basic block t r i d i a g o n a l structure of the J a c o b i a n . I n this w a y we c a n u t i l i z e the s o l u t i o n m e t h o d used i n s o l v i n g a d i a b a t i c , p r e m i x e d , l a m i n a r flames (17), b u r n e r - s t a b i l i z e d , p r e m i x e d , l a m i n a r flames (10), counterflow, l a m i n a r , diffusion flames (18), a n d the e x t i n c t i o n problems we consider here (2). F o r each value o f the parameter θ, we want to o b t a i n the corresponding value o f α a n d the r e m a i n i n g dependent s o l u t i o n components.

W e p o i n t out

t h a t , for each value o f the pseudo- arc length, a is constant. It satisfies the t r i v i a l differential equation £ ( « ) = 0

(4.6)

Hence, we c a n m a i n t a i n the block t r i d i a g o n a l structure o f the J a c o b i a n i n (4.5) if we introduce the parameter α as a dependent variable at m o f the m + 1 g r i d points a n d i f we specify a n o r m a l i z a t i o n c o n d i t i o n at t h e r e m a i n i n g g r i d p o i n t t h a t does n o t introduce nonzero J a c o b i a n entries outside o f the three block diagonals. T h e success o f this procedure depends u p o n the choice o f the normalization condition. In flame e x t i n c t i o n studies the m a x i m u m t e m p e r a t u r e is used often as the o r d i n a t e i n b i f u r c a t i o n curves. I n the counterflowing p r e m i x e d flames we c o n sider here, the m a x i m u m t e m p e r a t u r e is a t t a i n e d at the s y m m e t r y plane y = 0. Hence, i t is n a t u r a l t o introduce the t e m p e r a t u r e at the first g r i d p o i n t along w i t h t h e r e c i p r o c a l o f the s t r a i n rate o r the equivalence r a t i o as t h e depen dent variables i n the n o r m a l i z a t i o n c o n d i t i o n . I n this way the block t r i d i a g o n a l structure of the J a c o b i a n c a n b e m a i n t a i n e d . T h e final f o r m of the governing equations we solve is given b y (2.8)-(2.18), (4.6) a n d the n o r m a l i z a t i o n c o n d i t i o n

Ν =

rfr(

βθ)

°'

( Γ ( 0 , s)- Γ ( 0 , so)) + j - ( α ( 0 , s )) ( α ( 0 , β ) - a(0,s ))-(s-s ) 0

0

0

= 0 (4.7)

where a represents either 1/a o r φ. T h i s n o r m a l i z a t i o n is s u c h t h a t s a p p r o x i mates the arclength o f the s o l u t i o n b r a n c h i n the space ( T ( 0 ) , a).



412

Numerical Results In this section we a p p l y the adaptive b o u n d a r y value s o l u t i o n procedure a n d the pseudo-arclength c o n t i n u a t i o n m e t h o d to a set of strained p r e m i x e d hydrogenair flames. O u r goal is to predict accurately a n d efficiently the e x t i n c t i o n be h a v i o r of these flames as a function of the s t r a i n rate a n d the equivalence r a t i o . D e t a i l e d t r a n s p o r t a n d c o m p l e x c h e m i c a l kinetics are i n c l u d e d i n a l l of the cal culations. T h e reaction m e c h a n i s m for the hydrogen-air system is listed i n T a b l e


I

02).

S t r a i n R a t e E x t i n c t i o n . W e performed a sequence of s t r a i n rate calculations for a n 8.4% a n d a 9.3% (mole fraction) hydrogen-air flame. T h e equivalence ratios of these flames are φ = 0.219 a n d φ = 0.245, respectively. In b o t h cases the L e w i s n u m b e r of the deficient reactant (hydrogen) was significantly less t h a n one. I n p a r t i c u l a r , at the i n p u t jet, the L e w i s numbers were equal to 0.29 for b o t h the 8.4% flame a n d the 9.3% flame. W e also found t h a t these values d i d not change by more t h a n 15% t h r o u g h the flame. A n u m b e r of theoretical (5), (19-23), e x p e r i m e n t a l (24-28) a n d c o m p u t a t i o n a l (2), (23), (29-32), studies of p r e m i x e d flames i n a stagnation p o i n t flow have appeared recently i n the literature. I n m a n y of these papers it was found t h a t the L e w i s n u m b e r of the deficient reactant played a n i m p o r t a n t role i n the behavior of the flames near e x t i n c t i o n . In p a r t i c u l a r , i n the absence of d o w n s t r e a m heat loss, it was s h o w n t h a t e x t i n c t i o n of s t r a i n e d p r e m i x e d l a m i n a r flames can be accomplished v i a one of the following two mechanisms. If the L e w i s n u m b e r (the r a t i o of the t h e r m a l diffusivity to the mass diffusivity) of the deficient reactant is greater t h a n a c r i t i c a l value, L e > 1, t h e n e x t i n c t i o n c a n be achieved by flame stretch alone. In such flames (e.g., r i c h methane-air a n d lean propane-air flames) e x t i n c t i o n occurs at a finite distance f r o m the plane of s y m m e t r y . However, i f the L e w i s number of the deficient reactant is less t h a n this value (e.g., lean hydrogen-air a n d lean methane-air flames), then e x t i n c t i o n occurs f r o m a c o m b i n a t i o n of flame stretch a n d incomplete c h e m i c a l reaction. B a s e d u p o n these results we anticipate t h a t the Lewis n u m b e r of hydrogen w i l l play a n i m p o r t a n t role i n the e x t i n c t i o n process. c

F o r a n i n i t i a l value of the s t r a i n rate, the adaptive b o u n d a r y value m e t h o d was used to o b t a i n a s o l u t i o n o n a c o m p u t a t i o n a l d o m a i n of one c m . In F i g ures 2 a n d 3 we p l o t the t e m p e r a t u r e a n d the n o r m a l velocity component as a function of the distance f r o m the s y m m e t r y plane for the 9.3% (mole fraction) flame w i t h a s t r a i n rate of a = 200 s e c " . In F i g u r e s 4 a n d 5 we illustrate the major a n d m i n o r species profiles. A s we found i n the s o l u t i o n of freely p r o p a gating p r e m i x e d hydrogen-air flames, the H0 and H 0 profiles peak before the flame zone. T h e r e m a i n i n g r a d i c a l species peak i n the region of m a x i m u m temperature. O n c e the i n i t i a l s t r a i n rate c a l c u l a t i o n was completed, the m o d i fied arclength c o n t i n u a t i o n procedure was i m p l e m e n t e d to o b t a i n profiles of the m a x i m u m t e m p e r a t u r e versus the inverse of the s t r a i n rate. F o r each hydrogenair m i x t u r e , the e q u i d i s t r i b u t i o n procedure discussed i n (3.4)-(3.5) was used to determine the mesh for the c o n t i n u a t i o n steps. E a c h s o l u t i o n contained between 1

2

2

2


21.


0"

.

Counterflow Premixed Laminar Flames 413

1

•

•

1

0.10

0.20

I

0.30


Y (CM)

F i g u r e 2. T e m p e r a t u r e profile i n Κ for the 9.3% (mole fraction) hydrogen-air flame w i t h a s t r a i n rate of a = 200 s e c . - 1

1E-041

•

1 0.10

• γ

1 0.20

X_i

1 0.30

(CM)

F i g u r e 4. M a j o r species profiles for the 9.3% (mole fraction) hydrogen-air flame w i t h a s t r a i n rate of a = 200 sec"" . 1



414 T A B L E

I

Hydrogen-Air Reaction Mechanism R a t e Coefficients I n T h e F o r m k = AT exp(—E /RT). U n i t s are moles, cubic centimeters, seconds, K e l v i n s a n d c a l o r i e s / m o l e fi

f

A

β

Ε

1.70E+13 1.17E+09 5.13E+16 1.80E+10 2.10E+18 6.70E+19 6.70E+19 5.00E+13 2.50E+14 4.80E+13 6.00E+08 2.23E+12 1.85E+11 7.50E+23 2.50E+13 2.00E+12 1.30E+17 1.60E+12 1.00E+13

0.000 1.300 -0.816 1.000 -1.000 -1.420 -1.420 0.000 0.000 0.000 1.300 0.500 0.500 -2.600 0.000 0.000 0.000 0.000 0.000

47780. 3626. 16507. 8826. 0. 0. 0. 1000. 1900. 1000. 0. 92600. 95560. 0. 700. 0. 45500. 3800. 1800.


REACTION 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. a

H + 0 ^ 20H ΟΗ + Η τ± H 0 + Η Η + 0 τ±ΟΗ + 0 0 + Η τ±ΟΗ + Η Η + 0 + Μ τ=± H0 + M Η + 0 + 0 τ± H0 + 0 H + 0 + N ^ H0 + i V OH + H0 ^ # 0 -h 0 H + #0 ^ 20# Ο + # 0 τ=±0 + ΟΗ 20H ^ 0 + # 0 2

2

2

2

2

2

a

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

iJ 0 # if

6

+ M ^ / f + # + M + M ^ O + 0 + M + 0 # + M ^ H20 + M + H02 ^H2 + 02

2 2

H0 H0 H0 H0 2

+# 0 ^ H0 + 0 + Μτ±ΟΗ + ΟΗ + Μ + Ητ=± H0 + H + ΟΗτ± H 0 + H0 2

2

2

2

2

2

2

T h i r d b o d y efficiencies: h{o2) = o.

b

c

T h i r d b o d y efficiencies: 3&i (Ar).

2

2

2

2

2

2

2

0

k (H 0) = 21k (Ar), k (H ) = 3 . 3 * ( A r ) , fc (iV ) = 5

2

5

5

2

5

k (H 0) = 6k (Ar), k (H) = 2k (Ar), 12

2

l2

12

l2

5

2

k (H ) = 12

2

2

c

T h i r d b o d y efficiency: k (H 0) = 2 0 f c ( A r ) . u

2

14

50-70 adaptively chosen g r i d points. I n a d d i t i o n , the E u l e r - N e w t o n c o n t i n u a t i o n procedure was used t o help o b t a i n solutions (physical a n d unphysical) for b o t h increasing a n d decreasing values o f the s t r a i n rate. In F i g u r e 6 we plot the m a x i m u m temperature versus t h e inverse o f the s t r a i n rate for the 8.4% a n d 9.3% (mole fraction) hydrogen-air flames. E a c h C - s h a p e d e x t i n c t i o n curve was the result o f 84-93 c o n t i n u a t i o n steps. If we consider the 9.3% curve, we see t h a t as we move along the upper b r a n c h i n the d i r e c t i o n o f increasing s t r a i n rate, the peak temperature first increases a n d then decreases. U l t i m a t e l y , as the value o f dT /d(l/a) —• oo, t h e flames e x t i n g u i s h (a = 1443 sec "" for the 9.3 % flame a n d a = 760 sec for the 8.3 % flame). W e c a n , however, continue past the e x t i n c t i o n p o i n t w i t h the arclength procedure. W e find that, as the s t r a i n rate begins t o decrease, the peak t e m p e r a t u r e continues t o fall. I n the absence of H o p f bifurcation, the upper b r a n c h represents p h y s i c a l solutions a n d the lower b r a n c h u n p h y s i c a l solutions. T h e increase i n the temperature can be a t t r i b u t e d t o the a d d i t i o n a l heat release maz

1

ext

_

exi


1


21.


Counterflow

Y

Premixed

Laminar

Flames

415

(CM)

F i g u r e 5. M i n o r species profiles for the 9.3% (mole fraction) hydrogen-air flame w i t h a s t r a i n rate of a = 200 s e c . - 1

Î.600,

5

1E-03

5

1E-02

l / a (sec)

F i g u r e 6. C - s h a p e d e x t i n c t i o n curves i l l u s t r a t i n g the m a x i m u m t e m p e r a t u r e versus the inverse of the s t r a i n rate for the 8.4% a n d 9.3% (mole fraction) hydrogen-air flames.


416



t h a t results f r o m the diffusion of hydrogen (the deficient reactant) f r o m the u n b u r n t m i x t u r e into the reaction zone. T h i s excess heat release is larger t h a n the conductive heat flow f r o m the reaction zone to the u n b u r n t m i x t u r e . A s the s t r a i n rate continues to increase, however, the peak temperature finally decreases. E v e n t u a l l y , as the reaction zone is pressed against the stagnation surface, incomplete c o m b u s t i o n occurs due to the decreased residence t i m e a n d the flame extinguishes. It is interesting to note t h a t the adiabatic flame t e m p e r a t u r e of a freely p r o p agating (unstrained) 9.3% (mole fraction) hydrogen-air flame is 1046 Κ w h i l e the peak t e m p e r a t u r e i n F i g u r e 6 is a p p r o x i m a t e l y 1324 K . S i m i l a r l y , the adiabatic flame t e m p e r a t u r e of a freely p r o p a g a t i n g 8.4% (mole fraction) hydrogen-air flame is 977 Κ w h i l e the peak temperature i n F i g u r e 6 is 1264 K . Hence, we have L e w i s n u m b e r t e m p e r a t u r e increases of 278 Κ a n d 287 K , respectively. These peak temperatures as well as the e x t i n c t i o n behavior i l l u s t r a t e d i n F i g u r e 6 are i n excellent qualitative agreement w i t h the a s y m p t o t i c results of L i b b y a n d W i l l i a m s (22). G o o d quantitative agreement, however, is, i n general, not feasible (Giovangigli a n d C a n d e l , (23)). If we continue the arclength process such t h a t the s t r a i n rate decreases even further, we w i l l continue to move along the u n p h y s i c a l b r a n c h a n d eventually reach another t u r n i n g p o i n t - a n i g n i t i o n p o i n t . After the i g n i t i o n p o i n t is passed, we w i l l move o n to the extinguished s o l u t i o n b r a n c h . W e have not c a l c u l a t e d the i g n i t i o n points since they o c c u r at s t r a i n rates so low t h a t the c o r r e s p o n d i n g flames are p h y s i c a l l y uninteresting a n d it is impossible to stabilize such a flame i n the l a b o r a t o r y (see also S m i t h et a l . , (29)). L e a n a n d R i c h E x t i n c t i o n . E x t i n c t i o n of s t r a i n e d p r e m i x e d l a m i n a r flames can also be o b t a i n e d (for a given s t r a i n rate) by adjusting the equivalence r a t i o . In p a r t i c u l a r , e x p e r i m e n t a l results show t h a t as the equivalence r a t i o is either increased or decreased past a c r i t i c a l point, the flame extinguishes. S t a r t i n g w i t h a hydrogen-air flame h a v i n g a n equivalence r a t i o φ = 0.245 (9.3 % mole fraction) a n d a s t r a i n rate of a = 1000 s e c " w h i c h was o b t a i n e d d u r i n g the s t r a i n rate e x t i n c t i o n calculations, the E u l e r - N e w t o n adaptive c o n t i n u a t i o n procedure was used to o b t a i n solutions (physical a n d unphysical) for a l l values of the equivalence ratio. E a c h s o l u t i o n contained between 80-120 adaptively chosen g r i d points. 1

S t a r t i n g f r o m stoichiometric conditions (φ = 1) a n d then proceeding i n the lean d i r e c t i o n (φ < 1), we anticipate t h a t the peak flame t e m p e r a t u r e w i l l be reduced gradually. In a d d i t i o n , as the m a x i m u m t e m p e r a t u r e is lowered a n d the corresponding adiabatic flame speed of a n u n s t r a i n e d (a = 0) flame is re duced, we anticipate t h a t the flame w i l l move closer to the plane of s y m m e t r y . U l t i m a t e l y , as the fuel to air r a t i o is lowered below a c r i t i c a l value, r a d i c a l pro d u c t i o n i n the flame w i l l be severely restricted a n d the flame w i l l e x t i n g u i s h (lean e x t i n c t i o n ) . T h e arclength c o n t i n u a t i o n procedure w i l l then generate u n p h y s i c a l solutions for a d d i t i o n a l c o n t i n u a t i o n steps u n t i l a m a x i m u m value of the Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

21.


Counterflow

Premixed

Laminar

Flames

417

equivalence r a t i o is reached. A t this point (rich extinction) p h y s i c a l solutions w i l l again be generated a n d the m a x i m u m flame temperature w i l l b e g i n to rise a n d the distance of the flame f r o m the plane of s y m m e t r y w i l l increase as w e l l . T h e adaptive c o n t i n u a t i o n m e t h o d was able to generate closed response curves i l l u s t r a t i n g these phenomena. In p a r t i c u l a r , i n F i g u r e 7 we illustrate the m a x i m u m t e m p e r a t u r e versus the equivalence r a t i o for a hydrogen-air system - 1

w i t h a = 1000 s e c .

We observe a 1500 Κ temperature v a r i a t i o n a m o n g a l l

the flames c o m p u t e d . T h e lean e x t i n c t i o n o c c u r r e d at φ = 0.227 w i t h Τ = 1177 Κ a n d the r i c h e x t i n c t i o n at φ = 4.6 w i t h Τ =

1256 Κ . V a r i a t i o n s i n the


peak t e m p e r a t u r e occur over a m u c h smaller range of equivalence ratios for lean systems t h a n under r i c h c o n d i t i o n s . In F i g u r e 8 we illustrate the distance of the flame f r o m the plane of s y m m e t r y . T h e flame c a n shift its p o s i t i o n by almost 0.9 c m f r o m lean to slightly r i c h conditions a n d then back again to very r i c h conditions. I n p a r t i c u l a r , the p o s i t i o n of the flame f r o m the s y m m e t r y plane was equal to 0.03556 c m at the lean e x t i n c t i o n l i m i t , 0.8639 c m at its furthest p o i n t a n d 0.1073 c m at the r i c h e x t i n c t i o n l i m i t . A t o t a l of 2500 c o n t i n u a t i o n steps were used i n generating the closed curves i n F i g u r e s 7 a n d 8. In a l l cases the adaptive m e s h was able to m a i n t a i n resolution i n the flame zone as the

flame's

p o s i t i o n adjusted accordingly.

2,400 2,200 2,000 g

1,800

I *~

1,600 1,400 1,200 1,000 800 0.00

0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Equivalence ratio F i g u r e 7. E x t i n c t i o n curve i l l u s t r a t i n g the m a x i m u m t e m p e r a t u r e versus the equivalence r a t i o for hydrogen-air flames w i t h a s t r a i n rate of a = 1000 sec . - 1



418


Equivalence ratio Figure 8. Distance of the flame from the symmetry plane versus the equiva lence ratio for hydrogen-airflameswith a strain rate of a = 1000 sec"". 1

Acknowledgment The research was supported in part by Direction des Recherches et Etudes Techniques and the Office of Naval Research. We would like to thank the Scien tific Council of the Centre de Calcul Vectoriel pour la Recherche for providing access to a CRAY-1 computer and the Centre Inter Regional de Calcul Electron ique for providing computational resources. In addition, we would like to thank Professors S. Candel and G. Duvaut for their support and helpful discussions. References 1. Giovangigli, V. and Smooke, M. D. (1987). J. Comp. Phys. 68, p. 327. 2. Giovangigli, V. and Smooke, M. D. (1987). Comb. Sci. and Tech. 53, p. 23. 3. Smooke, M. D., Reuven, Y., Rabitz, H. and Dryer, F. L., to be published in Comb. Sci. and Tech., (1987). 4. Libby, P. A. and Bray, K. N. C (1980). Comb. and Flame 39, p. 33. 5. Libby, P. A. and Williams, F. A. (1982). Comb.andFlame 44, p. 287. 6. Bray, K. N. C., Libby, P. A. and Moss, J. B. (1984). Comb. Sci. and Tech. 41, p. 143. 7. Kee, R. J., Miller, J. A. and Jefferson, T. H. (1980). Sandia National Laboratories Report, SAND80-8003. 8. Kee, R. J., Warnatz, J., and Miller, J. A. (1983). Sandia National Labora tories Report, SAND83-8209. 9. Dixon-Lewis, G., David, T., Haskell, P. H., Fukutani, S., Jinno, H., Miller, J. Α., Kee, R. J., Smooke, M. D., Peters, N., Effelsberg, E., Warnatz, J. and Behrendt, F. (1982). Twentieth Symposium (International) on Com bustion, Reinhold, New York, p. 1893. 10. Smooke, M. D. (1982). J. Comp. Phys. 48, p. 72. 11. Deuflhard, P. (1974). Numer. Math. 22, p. 289. Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

21. SMOOKE & GIOVANGIGLI

Counterflow Premixed Laminar Flames 419

12. Smooke, M. D. (1983). J. Opt. Theory and Appl. 39, p. 489. 13. Keller, H . B . (1977). in Applications of Bifurcation Theory, P. Rabinowitz,


E d . , Academic Press, New York.

14. Heinemann, R. F., Overholser, K. A. and Reddien, G. W. (1979). Chem Eng. Sci. 34, p. 833. 15. Heinemann, R. F., Overholser, K. A. and Reddien, G. W. (1980). AIChE Journal 26, p. 725. 16. Chan, T. (1984). in Numerical Methods for Bifurcation Problems, T. Kupper, H. Mittelmann and H. Weber, Eds., Birkhauser Verlag, Basel. 17. Smooke, M. D., Miller, J. A. and Kee, R. J. (1983). Comb. Sci. and Tech. 34, p. 79. 18. Smooke, M. D., Miller, J. A. and Kee, R. J. (1985). in Numerical Boundary Value ODEs, U. M. Ascher and R. D. Russell, Eds., Birkhauser, Basel. 19. Sivashinsky, G. I. (1976). ACTA Astronautica 3, p. 889. 20. Buckmaster, J. (1982). Seventeenth Symposium (International) on Com bustion, Reinhold, New York, p. 835. 21. Libby, P. Α., Liñàn, A. and Williams, F. A. (1983). Comb. Sci. and Tech. 34, p. 257. 22. Libby, P. A. and Williams, F. A. (1984). Comb. Sci. and Tech. 37, p. 221. 23. Giovangigli, V. and Candel, S. (1986). Comb. Sci. and Tech., in press. 24. Fang, M., Schmitz, R. A. and Ladd, R. G. (1971). Comb. Sci. and Tech. 4, p. 143. 25. Law, C. K., Lshizuka , S. and Mizomoto, M. (1981). Eighteenth Symposium (International) on Combustion, Reinhold, New York, p. 1791. 26. Ishizuka, S. and Law C. K. (1982). Nineteenth Symposium (International) on Combustion, Reinhold, New York, p. 327. 27. Sato, J. (1982). Nineteenth Symposium (International) on Combustio Reinhold, New York, p. 1541. 28. Tsuji, H. and Yamaoka, I. (1982). Nineteenth Symposium (International) on Combustion, Reinhold, New York, p. 1533. 29. Smith, H. W., Schmitz, R. A. and Ladd, R. G. (1971). Comb. Sci. and Tech. 4, p. 131. 30. Sato, J. and Tsuji, H. (1978). in Proceedings of the Sixteenth Japanese Symposium on Combustion, p. 13. 31. Sato, J. and Tsuji, H. (1983). Comb. Sci. and Tech. 33, p. 193 32. Darabiha, N., Candel, S. and Marble, F. E. (1985). Comb, and Flame 64, p. 203. RECEIVED

June

15, 1987


Supercomputer Research in Chemistry and Chemical Engineering

Recommend Documents