Adaptive continuation algorithms with application to combustion problems

Applied Numerical Mathematics 5 (1989) 305-331 North-Holland

305

ADAPTIVE CONWWJATION ALGOIWIl-IMS WITH APPLICATION TO COMBUSTION PROBLEMS V. GIOVANGIGLI Laboratoire d’Acoustique et Mhcanique, Universit6 Paris VI, Paris, France; and Laboratoire EM2C du CSRS et de C’ECP Ecole Centrale des Arts et Manufactures, 92290 ChBtenay-Malabry, France

M.D. SMOOKE Department of Mechanical Engineering, Yale University, New Haven, CT 06520-2159, U.S.A.

We investigate the solution of parameterized nonlinear two-point boundary value problems by phase-space, pseudo-arclength, continuation methods that employ Newton-like iterations and adaptive gridding techniques. In particular, attention is focused on a modification of the basic continuation algorithm so that the Jacobian structure can be maintained both before and after the continuation reparameterization and on the determination of an adaptive mesh at each continuation step using equidistribution techniques. The numerical method is applied in the solution of counterflow premixed laminar flames with complex chemistry. We predict both strain rate and equivalence ratio extinction for a hydrogen-air system.

1. Introduction 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 govcming 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 parameters. For example, the parameters may include the strain rate or the equivalence ratio in a counterflow premixed laminar flame [17,18]. 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 a first-order sensitivity analysis to the physical system [%I. In other cases, the researcher may want to know how the system accuah’y 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 rate; this solution could then be used as a starting estimate for a predictor-corrector continuation process in the same parameter. In many applications this procedure works well. However, when there are multiple solutions and the solution curve, when plotted versus the parameter that is varied, turns back on itself by passing through a singular point, this procedure encounters difficulty. It is at these singular points, however, that the system often experiences special behavior. For example, in the solution of counterflosv premixed laminar flames, as the strain rate is increased past a critical value, the flame will extinguish. Similarly, for a given strain rate, as the equivalence ratio is either increased or decreased past a critical value, the flame wifl also extinguish. Hence, accurate knowledge of these points provides information about the flame’s flammability limits0168-9274/89/$3.50

0 1989, Elsevier Science Publishers 8-V. (North-Holland)

VI Govangigli, M.D. Smooke / Adaptive continuation algorithms

306

TO solve problems of physical interest numericaP one often discretizes the spatial operators of the governing ordinary or partial differential eal .G:LU- with the result that the solution of the original continuous problem is reduced to finding d * rp ,proximation to this solution at each point of a discrete mesh. If the solution contains 6gions of high spatial activity, it is also essential for reasons of both efficiency and accuracy that adaptive spatial grids be employed. The corresponding discrete nonlinear equations can then be solved by Newton’s method, or a variant thereof. However, in the neighborhood of singular points, i.e., turning points, where the derivative of the solution curve with respect to the parameter possesses a vertical tangent, Newton’s method encounters difficulties. To be able to follow the solution in the neighborhood of turning points, a modification of the basic algorithnr must be introduced. Two approaches have been discussed extensively in the literature. One method involves solving iteratively an inflated system with the turning point a unique and isolated solution and the other approach involves a path tracing continuation method in which points on the arc of solutions are computed [8,13,15,25,29,33,40.44]. In this paper we iocus our attention on the latter class of methods. We illustrate how a powerful combination of continuation techniques, Newton-like iterations and adaptive gridding can be used to solve laminar flame extinction problems. Interaction between the individual techniques is examined. In particular, we focus our development on the structure of the Jacobian matrices arising from the continuation process and on the influence of continuation on the adaptive gridding procedure. The resulting method can also be used to solve a larger class of parameterized nonlinear two-point boundary value problems with high gradients in the solution. In the next section the basic continuation algorithms are described. In Section 3 we consider two-point boundary value problems with a block tridiagonal structure and we show how to modify the formulation of the continuation process to maintain the same Jacobian structure both before and after the continuation reparameterization. In Section 4 we investigate the use of adaptive gridding continuation procedures that equidistribute a positive weight function with either a fixed number of points or a fixed equidistribution constant. Computational considerations are discussed in Section 5 and, finally, in Section 6, we illustrate the techniques with a set of detailed transport, complex chemistry, laminar flame extinction calculations.

2. Continuation algorithms Continuation usually consists in following a solution path by using an Euler method as a predictor and a Newton-like method as a corrector [8,13,29&l] with an appropriate reparameterization of the solution path. This reparameterization is in terms of a new independent parameter which locally ..pproximates the arclength of the solution in some phase space with the original parameter bp.oming an eigenvalue. The reparam\?terization is introduced in the form of an extra scalar equation approximating the corresponding arclength condition. The advantage of this procedure is that it suppresses Jacobian singularities at turning points where the solution path bends back on itself. More specifically, we consider a system of nonlinear equations in the general form

F(X,A)=O, where ,??R”xIR

-+R”

(2.1) is an appropriately

smooth

function,

X E IR!’ represents

the discrete

K Giovangigli, M.D. Smooke / Adaphve continuation algorirkms

307

solution to a physical problem and h E R is related to a physical parameter. We are interested in the dependence of the solution X on the parameter X, i.e., in tracing the solution branches ( X(h), A). Starting from a Freviously obtained solution ( X( A,), A,) one can extend the solution branch by using a one step Euler method as a predictor, i.e.,

x0@)= x@,) + (h -

h,)g(A,),

(2 .2)

to supply the initial iterate for the Newton corrector at X

F,(X'(X), X)(X'+'(x)-X'(X))= -F(X'(X),A), v=OJ,...~ where Fx is the Jacobian matrix. However, this type of continuation procedure fails points where F, is not invertible- in particular, near turning points. To overcome difficulty one has to reparameterize the solution arc (X, X) in terms of a new parameter s which usually approximates the arclength of the solution curve in some More precisely, we seek the solution

y(s) = ( X(S)T, Us))‘, of the augmented

(2.3) near singular this type of independent phase space.

(2.4)

system

W(49 Us)) N(X(s),

where the normalization write

X(s), s) 1

(2.5)

=”

I!’ : A” x 88 x R -+ BB determines

N=jF(s,))T(PY(s)-PY(so))-(s-so),

the dependence

of s on (X, A). We

(2.6)

in the neighborhood of a previously obtained solution Y(so) where P :R *+'---B R P is an appropriately smooth operator from IR‘+I into a phase IRP and where the projected tangent vector is normalized, i.e.,

(2.7) where all norms are taken to be Euclidean. Equation (2.7) shows that s is a local approximation of the arclength of the solution curve PY in the phase space R pL The Jacobian matrix of the augmented system in (2.5) can then be written in the form

(2.8) and it can be shown [29,44] tnat G, is always invertible if all of the singular points of Fx are turning points characterized by dim Null( Fx) = 1,

F,E Range( F,).

(2.9)

As a result, the predictor-corrector method in equations (2.2)-(2.3), or any related method, can be used in solving the augmented system in (2.5) thus allowing the calculation of solution branches with turning points.

K Giovangigii, M.D. Smooke / Adaptive continuation algorithms

308

3. Modified continuation algorithms for two-point boundary value problems The corrector step of the continuation procedure described in the preceding section is designed to solve the coupled system in (2.5). Most contmuation procedures use a Newton-like iterative method [8,12,13,29,44,45] (see [19] for least squares). A linear system invoIving an approximation of the Jacobian matrix G, must then be solved at each iteration. Our applications require direct solution methods, though we point out that in some applications iterative methods have been used (see, e.g. [5,9]). For systems of moderate size, any general solver can be used, but when F, is large and sparse one would like to exploit this structure in the solution algorithm. Although F, may have some initial structure (block-tridiagonal, banded, etc.) it is destroyed when the normalization is added to the original system (2.8). As a result, it is natural to seek modifications of the continuation process to maintain the Jacobian structure of Fx both before and after the reparameterization. In this section we derive such a modified continuation procedure for two-point boundary value probler7ls whose Jacobian matrix is block-tridiagonal. In Section 5 we show that the corresponding computational costs are comparable to the block elimination algorithms introduced by Keller [29-311 and Chan [5,7], providing the number of solution components is not to small. The method could also be apphed to different Jacobian structures, multidimensional problems, singular points of higher nullity, etc.. Finally, since we are interested in adaptive gridding, we only focus on grid independent normalizations. We consider the solution of a two-point boundary value problem depending on a parameter X on the interval [ a$]. The problem is written symbolically 9(5%-, X) = 0,

(3 . 1)

where %:[a,b]-,@ is the unknown solution and L is the number of components Z’, I= 1, 2 , . . . , L, of 5. In the modified procedure, we only consider normalizations which depend on %(t*) and X, where 5* is a given point of [a, b], which can be written symbolically as X(5%-(<*),

A, s) =o.

(3 .2)

3r(Lr, s) = 1

solving the coupled system in (3.1) and (3.2), we replace X by a function 2 of E and we let the unknown .5= (ZT, -Ep)T be the solution of a three-point limit value problem Rather

than

d.2’

= 0.

d5

i ,/tr(~(s*).

=H5*),

(3 .3)

s)

The main advantage of considering (3.3) is that, if a3.1) can be discretized in a block-tridiagonal form, then (3.3) has the same property. Indeed, on a given grid E Z= (a=(, let k


.-.

c[Jm-‘?),

(3 .4)

be the index such that

(3.5) and let us denote H = (HIT,. . . , k(?jT, (

x:.-* ’

,

X,‘)‘, A -(A,,...,

F

=

(FIT,.

. . , FJT)',

h,)T and N as discrete versions of X,

Z =

(z,7; . . . , z,‘)~, x =

9,

%“, %, L? and M,

V. Giovangigli, M.D. Smooke / Adap. we continuation alaorithms

respectively. The discrete problem corresponding be written H(Z,

to (3.3) then has ( L + l)J unknowns

s) =o

39

and can (3 .6)

with Zi = ( Xj’, X/)T and Jfj=(I;;T(xj-I,

xj9

xj+*9

xj)V

Mj)T9

(3

.7)

where Aj+l Mj=

-xj9

N( Xk_,, Xk, Xk+l, A,, 4, Xj- xj-],

if l
(3 .8)

if j= k if k
(3.9) (3.10)

The dummy indices 0 and J + 1 are introduced to avoid notational complexity. The matrix H, has a block-tridiagonal structure and one can prove easily that (3.1)-(3.2) and (3.3) are equivalent and that det H, = ( -l)k-*

det G,,

(3.11)

where G and Y are as in (2.4)-(2.5), thus allowing the use of the enlarged system in (3.3). As in (2.7), the discrete normalization N is taken to be (3.12) in the neighborhood of a previously converged solution Z( sO) where Q is an appropriately smooth operator from RtL+ I)’ into a phase space !# p which only depends on Xk_ *, Xk, X, + 1 and xk and where the projected tangent vector is normalized, i.e., ]IdQZ/ds(s,) 11= 1. We point out that the procedure of introducing the dummy equation dP/d< = 0 is similar to that used in 1591to simplify the Jacobian structure in nonlinear eigenvalue problems and to that used in [14] to reduce the sensitivity of multiple shooting problems with homotopy. In addition, Rheinboldt [44] has pointed out that normalizations involving only one component of the discrete solution have the numerical advantage of introducing more zeros in the Jacobian matrix. However, it is the combination of these techniques that simplifies the Jacobian structure after the reparameterization. Finally, we note that, if the normalization involves the solution at every point but can be put in the following integral form J(r,x(

jbaT(S)+(S)

(3.13)

de, A, s),

a

where ~:[a, b] +RL is an appropriately smooth function, then it is still possible to obtain a block&diagonal Jacobian structure by introducing In( 5) = /,Eo’( t )%(t) dt and discretizing the system

Sk(a) = 0, The corresponding

M&?(b),

2’(b),

s)=O.

discrete system has (L + 2) J equations and ( L + 2) J unknowns,.

(3.14)

310

V. Giovangigli. M.D. Smooke / Adaptive continuation ulgorithms

4. Adaptive continuation algorithms for two-point boundary value problems The adaptive selection of grid points is of critical importance to the ;r”ficiency of a continuation algorithm that is used to solve problems with regions of high spatial activity. In this section we summas-ize some of the essential features of adaptive grid selection procedures for two-point boundary value problems. We then combine the gridding and the continuation algorithm to obtain an adaptive continuation algorithm. We only focus on procedures that adaptively determine grid points by equidistributing a positive weight function (see, e.g., [18,23,26,35, 42,48,53-56)). We consider two-point boundary value problems depending on a parameter. After a grid-independent local parameterization the continuous problem can be written in the form X(9,

s)=o,

. (41)

where 9’: [a$] -B RL+’ is the unknown function and L + 1 is the number of solution components. We assume that a positive weight function w has been selected to determine the grid points using equidistribution. We say that a weight function w, discrete or continuous, is equidistributed by the grid Es on [ a,b] ?=

IQ=&*


with respect to the constant

r:+

1

/

w d[=C,

c

l

--

<&=b),

(4 .2)

C if j=1,2

(4.3)

,..., J*-1,

4:

where the integral sign denotes either a discrete or a continuous integration operator. The weight functions used in the adaptive selection of grid points for two-point boundary value problems are often related to the derivatives of the solution (e.g., gradient and/or curvature, local truncation errors, arclength of the solution graph, etc.) and to grid smoothness properties. Anticipating continuation parameter sensitivity problems, we let w be a function of 9 and its partial derivatives with respect to < and S, i.e.,

The discrete version of (4.1) on a grid E = ( a = 5, < & < H(Z,

l

-

l

< 55 = b) is denoted

S, -“) = 0,

where H: R(L+l)J x W x BP'--)WL+IJJ and 2 E W(L+l)J is a discrete corresponding discretization of YV is denoted by W w=

w(z,

E),

function

by (4 .5) on E. The

(4 .6) and its components are functions of Z and z. We assume that X and H are appropriately smooth and that they have locally smooth arcs of solutions which we denote by .9’(s) and Z#, respectively, for all the grids Z and all the local parameterizations s considered. We also assume that a Newton-like iterative method will be used to solve (4.5) for Z,(S) on a given grid z provided an initial guess Z:(s) is available in the domain of convergence of the method. Equidistribution can t=2 used to select grid points adaptively with either a fixed number of nodes or a fixed equidistribution constant. Indeed, for a given Z and z and for a given number

I/. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

311

of points J* there exists a unique grid Z* of J * points such that (4.3) holds with IV= W( 2, E) as in (4.6) and, alternatively, for a given equidistribution constant C, there exists a unique grid Z* such that (4.3) holds except for the last interval [&+ er* 1. In both cases Z* is obtained easily be inverting a discrete integration operator. Assuming that a strategy has been selected as to whether a fixed number of points or a fixed equidistribution constant is to be used, the following procedure can then be used to adapt the grid to the solution for a fixed value of s. We first select an initial grid Z0 and an initial guess Z&(s). Next we solve (4.5) for ZZO(s) starting from ZiO(s). We then compute Ei by equidistributing IV( Z(s), _‘) and we interpolate ZZO(s) onto Ei to provide an initial guess Z,(s) on Ei. Finally we solve (4.5) for Z,,(s) starting from Z:,(s). This process can be repeated a small number of times but during a continuation process one iteration is generally sufficient. We note that this adaptive grid procedure can be classified as static in the sense that the continuation parameter s is not changed during the procedure. A number of choices for the equidistribution weight function YY, have appeared in the literature. White has examined equidistributing the arclength of the solution graph [60,61], Pereyra and Sewell have equidistributed the local truncation errors [42] and Pearson [41) has subequidistributed the components 9’ of 9’ by employing YV= ]Id%“/d[ 11up to scaling factors. However, when L t- 1 is large and when all but a few of the variables peak in one region then, if the arclength or the Pearson procedure is employed, the main activity region will be well resolved but the secondary activity region will be underresolved. In particular, for combustion problems of physical interest, we often desire solution profiles to many different chemical species (e.g., methane-air combustion can require over 30 species with well over a hundred chemical reactions). Some of these species profiles may peak at spatially separated regions of the computational domain. As a result, to resolve each of the components, we employ a weight function that determines the grid point locations by combining gradient and curvature properties of each component of the solution. We define (4.7)

number smaller than one where 9’ is a family of weights w and E, is a corresponding (0 < e,,,< 1). This weight function 1Iy’ is used with a fixed equidistribution constant C = 1 and the family of weights 9’ is generally formed from the weight functions w = I, which prevents the subintervals from becoming too large, and w = d@/d[ and w = d*&/dt*, 1 = 1,. . . , L, which subequidistribute the solution’s components and their gradients, respectively 1181.It is important to note that the type of weight function in (4.7) can incorporate any new effect by enlarging the family of weights 9’ among which the maximum is taken. For instance, we have also used the weights w = 8*9“‘/a@s to account for continuation parameter sensitivity effects and to obtain more accurate discrete tangent vectors dZ,/ds. Note that the corresponding $V is independent of the locai paramekrization s because of the scaling factor in (4.7). In addition, we note that for multiple parameter problems, e.g., fold calculations [45] or critical point sensitivity analysis, one 2 / a@v for any control parameter v to obtain accurate can add families of weights like w = a2@01 tangent vectors in any direction. A potential problem with adaptive grid selection using equidistribution is the formation of a mesh in which the size of consecutive mesh intervals varies too rapidly. To alleviate this problem one would like the mesh to be locally bounded, i.e., the ratio of two consecutive intervals is

Y. Giovangigli, M.D. Smooke / Adaptive continuation aigorithms

312

bounded above and below by a constant [l&26.53]. We write

where R >, 1 is a given positive constant. When the fixed equidistribution used, one can check easily that (4.8) holds whenever the Lipschitz constant than log R/C. Since the weight functions w that we consider do not property, one has to replace them by the smallest fdnztion 3 greater than

s?-(E)= sup

{(~‘(l)+(logR/C)l~-tJ)-I),

constant strategy is of l/w is smaller usually satisfy this w which does, i.e.,

(4.9)

tE[U*h]

which was proposed by Kautsky and Nichols [26] using a somewhat different argument; equation (4.9) also holds for the corresponding discrete weight functions. Due to the complexity of the solution structure Z?‘, e.g., a complex chemistry laminar flame front with all the intermediate products peaking at different places, the final weight function @ may not be very smooth in 5. It is then useful to regularize # by using a convolution operator [34] with a Gaussian kernel. This procedure is employed for the combustion applications discussed in Section 5. We point out that there are other gridding strategies besides the global refinement of (4.3) with either a fixed number of points or a fixed equidistribution constant. For example, one can subequidistribute the weight function ‘#Y such that /

‘:“wdt
j=1,2

,...,

J*-1,

(4.10)

using successively refined grids Zm, m = 0, 1,. . . , where E,,, is obtained from Em_ 1 by only adding points, say in the middle of each subinterval of .E,,,_i when (4.10) does not hold [17,53,54,57]. These local refinement methods are very efficient if little is known about the solution. When calculating a first solution, the initial grids are relatively coarse since the positions of the steep fronts are not known a priori with high accuracy. The corresponding initial guesses Zim are, in general, not accurate and the greatest number of Newton-like iterations are then performed on these grids. As the mesh is refined, however, a more accurate starting estimate is available and so a smaller number of Newton iterations are crdinarily required. If this type of adaptive gridding procedure is employed in a continuation process, then the initial grid E0 for a new solution will be the previous solution’s grid [46]. One could also use a skeleton grid procedure [17,55,56] in which only a fraction of the points of the previous solution’s grid is contained in the starting mesh ZO. In the former case, this can quickly result in an inefficient algorithm with a constantly increasing number of nodes, e.g., a moving flame front and, in the latter case, one can obtain convergence difficulties with highly nonlinear problems [18]. In addition, in bcth of these cases neither the number of refinements nor the number of points are minimized. On the other hand, the global equidistribution procedure (4.3) can be applied directly during a continuation process without any convergence difficulties since information is known about the solution from previous continuation steps. As a result, the number of refinements and the number of points are optimized and the grid moves with the solution in a very natural way. We further remark that, in general, it is preferable to have an explicit adaptive gridding procedure, as discussed in the preceding sections, where the grid determination is uncoupled

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

313

from the system solution. This leads to low cost grids that are approximately equidistributed. It also allows flexibility in the choice of the weight functions and the use of a variable number of points [23,35,42,53-561. Finally, we note that the interesting idea of performing a continuation process on a coarse grid together with a refinement procedure only at a few selected solutions, as it is done in multigrid continuation techniques, for instance, cannot be used for our combustion applications. This procedure assumes that the qualitative behavior of the solution is the same on both the coarser and the finer grid [2]. This is not the case for laminar flame problems since, small temperature perturbations in a flame can induce large variations on the flow field and the front location” Hence, it is essential to use a continuation process on fine grids adapted at each continuation step. The ideas of the preceding sections can be combined to produce an adaptive continuation algorithm. Starting from a previously obtained discrete solution Z,,, = ZO,, with ZOoiaadapted to Zola, we want to obtain a new discrete solution Z,,,, on another discrete solution path ZZnM< l), with E_ adapted to Z,,. We have the following algorithm Adaptive Continuation Procedure 4.1. Step 1. Compute the oriented unit tangent vector 7Z0,dof the solution path Z_( 0) at Z,,,. step 2. Select a normalization N and set ZZO,,(sO) = Zold. Compute a steplength As = s - sO. Scale 7Z0,dto get d ZEO,,(s,,)/ds. Step 3. Compute the predictor ZP on a predicted grid EP. Step 4. Compute the corrector ZL,( s) on Znnewadapted to Z=“.,(s), starting from (E,,, 2,). If it fails to converge, adjust As and go to the predictor step. Step 5. If a turning point is located, compute it. A detailed study of each step is examined in the next section. 5. Computational considerations There are a number of computational details that are critical in the successful implementation the adapiive continuation algorithm described in the preceding two sections. We focus on five 1 specific areas, namely the calculation of the unit tangent vector, the formation of the predictor, the elimination algorithm used in the corrector the choice of the local parameterization, iterations and the calculation of turning points. I of

5.1. The unit tangent vector T, using any local parameterization

N, the COrreSpORdiRgtangent vector can be obtained from

(5.1) where Z denotes the mesh. If N is taken to be the local parametelization of the previous continuation step, then, since the corresponding Jacobian matrix H, has already been decomThe posed during the corrector iterations, only olle backsolve is needed to form dZ,/ds. oriented unit tangent vector 2, is obtained by normalizing dZ,/ds. Altough LZis adapted to Z,, it may not be adapted to (dZ,/ds) and the corresponding tangent vector will not be accurate.

314

V. Giovangigli, M.D. Smoke

/ Adaptive continuation algorithms

This may happen for instance with moving steep fronts. This phenomenon has already been observed in (471 for first-order sensitivity calculations. This problem can be solved by including the weights w = a2%o’/a@s in the family 9 in (4.7). However, if the problem is purely translational, i.e., the shape of the solution is shifted but unchanged during the continuation process, then the weights w = a2S”‘/at2 and w = a2.S”‘/a@s have the same effect. 5.2. Local parameterizations and steplengths The parameterizations found in the literature essentially correspond to different choices of the operator Q in (3.12), Awhich is often an orthogonal projection. For the simple choice Q = Q:R (t+l)J 3 Wtl+ 1, QZ=(XT, which corresponds to P = 1 in (2.6), equation (3.12) &IT, specializes to

(5.2) N={~~so)}T~X(s)-X(so~)+(~~so,](h which was introduced in the general framework of (2.5) by Keller [29]. In this situation s is a Iota! approximation of the arclength of the solution curve Y( ), i.e. !I(d Y/ds)(s,) I}= 1. We note, however, that (5.2) is not the discretization of normalizations like (3.2) or (3.13) but using L, norms for 5 leads to normalizations of the form l

I

(5 4 .

which can be used as in (3.14). If we denote (e:>, j= 1, 2 ,..., J, d= 1, 2 ,..., L + 1, the natural basis vector of RtL+l)j, choice Q : WL+ w - W, such that Q = el, QZ = z:, leads to N=(-~)“‘(z:(s)-z:(s,))

the

-(~--~)r

(5 .4) where ( - 1)” = (deLTZ/ds)(so). In this situation s approximates either a component of X or X up to d sign variation. This parameterization has been introduced in the general framework of (2.5) by Abbott [l], Kubicek [33j and Rheinboldt [44]. However, (5.4) is not grid-independent unless j = 1, but if a_,, a0 and a, are suit able elements of R , and, if we let Z’ = a _ ,eL_ I + a,eL + V:+P then the choice Q = CtT, QZ = W, leads to ~~=(-l)m(V(+cp(so))

-(s--01,

where 43= a_,~:_~ + a,zi + arz:,, and (- 1)” = (dq/ds)(s,). sign, an interpolant p, which approximates %I( 5 * ) with

In this situation

(5 .5) s is, up to a

[*=a

-A-* +ao5,+alrk+,* (5 .6) This type of grid-independent normalization is important for adaptive gridding problems and is used for the combustion applications discxssed in Section 6. If 1 < 1,~ L, the choice Q: R(t+l)J --) R2, Q = (e:, ek+*)T, QZ= (XL, X,)T, leads to

N=

(~(so)J(x:(s) --x:(so)l+ (~bo~](X,(s) -X,(s,))

-(s-so).

(5.7)

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

315

This normalization has been used by the authors [17,18] (Section 6), for combustion applications, with k= 1 in which case &=[* = u is also fixed, together with a nonlinear transformation for A. This choice was motivated by the fact that, as pointed out by Bank and Ghan [2], for many problems there are really only two parameters that affect the continuation process-one measure of X and one measure of X. Equation (5.7) can be ,made grid-independent by using Q = (g’, e,L+l)T. Aside from physical considerations, there are two important criteria in the selection of a normalization N. The first one ns related to the structure of the Jacobian matrix as shown in Section 3 and the second one is a condition number consideration as described below. Indeed, one can deduce from (3.11) and a result in [44] that at 2 = Z,,, we have det H,=

~-~)k~1~liQ~~~11/1t$2~~ll,

(5 .8)

where 9 is nonzero and is the determinant of G, when Q = & is as in (5.2), i.e., when P = I. Hence, not only ]]Q~T~ I] must be nonzero, but }IQzrs I] must not be too small. In particular, if QZ = eLTZ= z: as in (5.4), the indices k and I should be chosen such that ] T; 1 = 1eL“& 1 is maximized [l&t]. If QZ = zlTZ = cp as in (KS), 6 * and I should be chosen such that ] PiTrz ! is maximized and, similarly, if QZ = (XL, X,)T as in (5.7), then the indices k and I should be chosen such that ((T/)’ + (~,f+l)~)l/~ is not too small. We note also that among all linear orthogonal projectors Q which do not depend on Xj with j + k, ]det Hz 1 is maximized only for Q = Q as in (5.2). Steplength algorithms that automatically adjust As play an important role in continuation algorithms and generally are derived by bounding the distance between the solution and the predictor by an estimated safe radius of convergence of the corrector iterative process [ll-13,20,24] 5.3. The predictor The simplest choices for the predicted grid zr., and the predicted previously obtained grid and to use an Euler predictor

solution 2, are to take the

(5.9 the tangent vector being obtained from dZzO,d/ds = T&I] Q~T~*,,]I. Chan [8] has also introduced second-order Taylor expansion predictors and Bank and Chan [2] have used different steplengths As, and As, in the different directions X and X of 2. Furthermore, if zr., is different from zoold,then dynamic gridding procedures must be used (see, e.g., [55,56]). In our applications, however, we have only used the simple choices (5.9). 5.4. Corrector iterations If we denote the block structure of Hz by

(5.10)

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

316

then the following block elimination

.

.

can be used to decompose

A=

Hz in J%Y=LPP

.

%? k-l

1

41 1

% . .

. .

.

.

1 (5.12)

4Y= 1 .

.

.

where ~j, ~j and @, are given by

One can check easily that (5.13) is feasible or.Zy when det 2j+Oj

j=l,

2 ,...,

9.

(5.14 ,

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

317

Now if we introduce

d%Yj=

.

.

~j

9

=

9

iaJ-l

ksJ-l

(5.15) then we have

det 2j=

det Aj/det

~j_l,

det A/(det

Ak_l

if l
det O,,,),

if k
det @j/det oi+l, where we define det A,-, = det matrices of Fx by

A,

@J+

I

(5.16)

=

1. On the other hand, if we denote the corresponding

4

c2

F,=M= l

.. .. . .. .

BJ-I

cJ

L

(5.17)

,

AJ

and

4

Aj

4

c2

. . . . . . . .

Mj =

Cj

Bj

Cj+1 oj = Bj-l

Aj

.. . .. . .

9

l

B J-l cJ

‘J (5.18)

then one can show easily that det Mj,

det AV~=(-~)‘-’ det @j=det

Oj,

l
(5.20)

k
Hence we deduce that the decomposition

(5.n9j

in (5.14) can be done when (5.21)

det GY f 0, det MjZO,

l
(5.22)

det Oj # 0,

k
(5.23)

K Giovangigli, M.D. Smooke / Adaptive continuation algorithm

318

Equations (5.22) and (5.23) hold as pointed out by Keller [31], and (5.21) is a result of the reparameterization. Note also that the usual block&diagonal elimination algorithm corresponds to k =J. Several block elimination algorithms have been proposed to solve systems similar to (2.8) with the same solver used for both F’ and G,. These algorithms are the block eiimination algorithm (BE) of Keller [29], the singular block elimination algorithm (SBE) of Keller [30,31] and the deflated block elimination (DBE) of Chan [5,7]. Algorithm BE can always be used except when the singularity of F, is of the same order of magnitude as the square root of the relative precision. In this case algorithm SBE must be used instead. 0n the other hand, algorithm DBE is valid independently of the singularity of F,. In addition, LU decompositions with small last pivots and special pivoting strategies have to be used for SBE and DBE. The corresponding costs are one decomposition and two backsolves for BE and SBE and one decomposition and four backsolves for DBE. We refer to [5-9,29-321 for more details. Altough a corrector iteration for the system in (3.6) is simpler than a correction iteration for the one in (2.5) using algorithms BE and SBE or algorithm DBE, the linear system to solve is larger. Therefore it is necessary to compare the costs of the thret algorithms when applied to a block-tridiagonal matrix F, of size LI (L being the size of the blocks and J the number of block lines) and the cost of solving a block&diagonal system of size (L + 1)J. If we define one operation as an addition and a multiplication, then it can be shown that it takes, to leading order , 3L3J operations to factor a block-tridiagonal matrix of size LJ. A forward and backsolve takes, to leading order, 3L2J operations. Hence, the BE and SBE algorithms take a total of f L3J + 6L2J operations, the DBE algorithm takes $ L3J + 12 L2J operations and the algorithm proposed in (3.6) takes 5 L3J + 1OL’J operations. In particular, for L large, as is typically the case in combustion systems, we see that the operation costs for the four methods are comparable. 5.5. Turning points Simple quadratic turning points are characterized

by (2.9) and (5.24)

where cp and ri/ are right and left nullvectors check easily that dX ds- = 0

and

-d2X #O ds2 ’

of F, [45]. When (2.9) and (5.24) hold, one can

#O-25)

where s is any local parameteriz;.lon similar to .e ones in (2.6)-(2.7), or (3.12). As a result, these turning points are detected easily by looking at the sign of dX/ds and they can be computed by using Newton’s method to find s such that dh/ds = 0, the local parameterization being kept fixed during the process. This is essentially the method recommended by Chan [8]. Note, however, that for a finite difference approximation of d2X/ds2 the step in s must not be smaller than the square root of the relative precision of dX/ds, which may be affected by the gridding procedure. In this situation we can hold d2A/ds2 constant as in a fnodified Newton’s method. See [8] and [40] for excellent surveys on turning point calculations.

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

319

6. Application to combustion problems 6. I. Problem formulation

Strained premixed laminar flames have played an important role in recent theories of turbulent premixed combustion. The reacting surface in such models can be viewed as being composed of a number of laminar flamelets (see e.g., [3,36-381). The flamelets are often modeled by considering counterflowing streams of reactants or counterflowing streams of reactants and products. From an experiment viewpoint, these flames can be produced when a single reactant stream impinges on an adiabatic wall or when two counterflowing reactant streams (or a leactant stream and a product stream) emerge from two counterflowing opposed jets. In the neighborhood of the stagnation point produced by these flows, a chemically reacting boundary layer is established. The governing conservation equations can be reduced to a set of nonlinear two-point boundary value problems along the stagnation point streamline. If there are two reactant streams and each has the same exit velocity and equivalence ratio, then a double flame is produced with a plane of symmetry through the stagnation point and parallel to the surface of the two jets (see Fig. 1). Application of nonadaptive continuation techniques to combustion problems with simplified kinetics can be found in [10,21,22,43,50,51]. In this section we apply the adaptive continuation procedure to a set of doubly premixed laminar flames in a planar geometry. Detailed transport and complex chemical kinetics [16,27,28] are included in all of the calculations. The reaction mechanisms for the hydrogen-air [54] and methane-air [39] systems are listed in Tables 1 and 2, respectively. The governing equations for mass, momentum, chemical species and energy can be written in the form [18] dV +apU=O,

(6. 1)

dY

- Vg

(6.2)

+a(pJ+pU2)=0,

Plane

0

X

Fig 1. Schematic of the stagnation point flow configuration.

V. Giovangigli, MD. Smooke / Adaptive continuation algorithms

320

d -- dy(Py,v,V)-

dYk vdy +G,w,=O,

(6.3)

k=l,2,...,K

(6.4 p

(6.5)

=pW/RT,

where y denotes the spatial coordinate normal to the stagnation plane (see Fig. 1); T, the temperature; Yk, the mass fraction of the k th species; K the number of species; U a similarity function related to the tangential velocity; V the density times the normal velocity; p, the mass density; p, the pressure; wk, the molecular weight of the k th species; %, the mean molecular weight of the mixture; R, the universal gas constant; X, the thermal conductivity of the mixture; c,,, the constant pressure heat capacity of the mixture; cr,k, the constant pressure heat capacity of thic: k th species; 6k, the molar rate of production of the k th species per unit volume; p the viscosity of the mixture; h,, the specific enthalpy of the k th species; and vky, the diffusion velocity of the k th species in the y-direction. The velocity compor,ents u and u in the transverse (x) and normal (y) directions, respectively, can be obtained from the solution vector (K U, r,,..., Yk, T) by forming u = axU,

(6 .6)

v=

(6 .7)

v/p.

Complete specification end of the computational v=o,

of the problem requires that boundary domain. At y = 0 we have

dU dy=O,

dYk

dy

=0,

k=l,2

,...,

K,

conditions -= dT dy

0

be imposed at each

,

(6.8)

and as y --, 00

(6.9) where the temperature T= and mass fractions Yk,( (p), k = 1, 2,. . . , K, at the edge of the boundary layer are specified. We observe that the computational model depends on two parameters a and +. The quantity a denotes the strain rate and (P represents the equivalence ratio of the incoming reactant stream. The strain rate is a measure of the stretch in the flame due to the imposed flow and + is a measure of the relative proportion of fuel to air in the reactant stream. If the strain rate is increased past a critical value the flame will extinguish. Similarly, if the equivalence ratio is either increased or decreased past a critical point, the flame will extinguish. Hence, accurate knowledge of the behavior of counterflow premixed flames as these quantities are varied can provide information about the flame’s flammability limits. 0ur goal is to predict accurately and efficiently these extinction limits. This information can then be used in theories of turbulent premixed combustion. 6.2. Strain extinction 6.2. i. Hydrogen -air ,f&.m E.Y We performed a sequence of strain rate calculations for an 8.4% and a 9.3% (mole fraction) hydrogen-air flame. The equivalence ratios oii these flames are C#B = 0.219 and + = 0.245, respec-

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms Table 1 Hydrogen-air reaction mechanism rate coefficients centimeters, seconds, Kelvins and calories/mole

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

in the form kf = ATB exp( - E,/RT).

Reaction

A

Hz +O, + 20H OH+H,*H,O+H

1.70E+ 13 l.l7E+09

H+O,+OH+O O+H ,+OH+H H+O,+M=HO,+M” H+02+0,*H02+02 H+0,+N,*HO;!+N2 OH+HO, s H,O+O, H+HO,=20H O+HO,sO,+OH 20H+O+H,O H,+Mc=H+H+Mb O,+Mz=O+O+M H+OH+Mc=H20+Mc H+H02r--HHz+02 HO, + HO, * Hz02 +O, H,O,+M+OH+OH+M H,O,+H+HO,+H, H,O, +OH c= H,O+HO,

5.13E+ 16 1.80E + 10 2.10E+ 18 6.70E + 19 6.70E + 19 S.OOEc 13 2SOE + 14 4.80E + 13 6.OOE+ 08 2.23E + 12 1.85E+ 11 7.50E + 23 2.50E + 13 2.OOE+ 12 1.30E+ 17 1.60E+ 12 l.OOE+ 13

a Third body efficiencies: k,(H,O)

= 2lk,(Ar),

k,(H,)

b Third body efficiencies: k,,(H,O) = 6k12(Ar), k,,(H) ’ Third body efficiency: k,4(HZO) = Zlk,,(Ar).

= 3.3k,(AQ = 2k,,(Ar),

k,(N,) k,,(H,)

Units are moles, cubic E

P --

321

0.000 1.300

47780. 3626.

- 0.816 1.ooo - 1.ooo - 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

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

= k,(o,)

= 0.

= 3k,,(Ar).

tively. In both cases the Lewis number (the ratio of the thermal to mass diffusivity) of the deficient reactant (hydrogen) was approximately Le = 0.29. Based upon the experimental work of Ishizuka and Law [24], Sato [49], the theoretical work of Buckmaster [4], Libby and Williams, [37,38], Sivashinsky [52], and the computational work of Sato and Tsuji [50,51], we anticipate that the Lewis number of hydrogen will play an important role in the extinction process. In particular, in the absence of downstream heat loss, if the Lewis number of the deficient reactant is greater than a critical Lewis number Le, greater than one, then extinction can be achieved by flame stretch. In such flames extinction occurs at a finite distance from the plane of symmetry. However, if the Lewis number of the deficient reactant is less than this critical value Le,, as in the case of lean hydrogen-air flames, then extinction occurs from a combination of flame stretch and incomplete chemical reaction when the flame is in contact with the stagnation surface. For an initial value of the strain rate, an adaptive boundary value method 153,571 was used to obtain a solution on a computational domain of one cm. Once the initial strain rate calculation was completed, the Euler-Newton adaptive continuation procedure, with the equidistribution procedure discussed in (4.7)-(4.9), was implemented to obtain solutions (physical and unphysiCal) for both increasing and decreasing values of the strain rate. Each solution contained between 50-70 adaptively chosen grid points. A normalization of the type (5.7) was used with k = 1 and 2 = K + 3, namely, QZ = (T(O), l/a(o))T [18].

322

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

In Fig. 2 we plot the maximum temperature versus the inverse of the strain rate for the 8.4% and 9.3% (mole fraction) hydrogen-air flames. Each C-shaped extinction curve was the result of 84-93 continuation steps. Some inaccuracies are due to grid noise in the solutions. If we consider the 9.3% curve, we see that as we move along the upper brcr-h in the direction of increasing strain rate, the peak temperature first increases and then decreases. Ultimately, as the value of dT,,/d(l/u) + do, the flames extinguish ( aext = 1443 set-’ for the 9.3% flame and aext = 760 set -r for the 8.3% flame). We can, however, continue past the extinction point with the arclength procedure. We find that, as the strain rate begins to decrease, the peak temperature continues to fall. In the absence of Hopf bifurcation, the upper branch represents physical solutions and the lower branch unphysical solutions. The increase in the temperature can be attributed to the additional heat release results from the diffusion of hydrogen (the deficient reactant) from the unburnt mixture into the reaction zone. This excess heat release is larger than the conductive heat flow from the reaction zone to the unburnt mixture. As the strain rate continues to increase, however. the peak temperature finally decreases. Eventually, as the reaction zone is pressed against the stagnation surface, incomplete combustion occurs due to the decreased residence time and the flame extinguishes. It is interesting to note that the adiabatic flame temperature of a freely propagating (unstrained) 9.3% (mole fraction) hydrogen-air flame is 1046 K while the peak temperature in Fig. 2 is approximately 1324 K. Similarly, the adiabatic flame temperature of a freely propagating 8.4% (mole fraction) hydrogen-air flame is 977 K while the peak temperature in Fig. 2 is 1264 K. Hence, we have Lewis number temperature increases of 278 K and 287 K, respectively. These peak temperatures as well as the extinction behavior

Table 2 Methane-air reaction mechanism rate coefficients centimeters, seconds, Kelvins and calories/mole

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

in the form k, = ATp exp( - E,/RT).

Reaction

A

CH,+Mc’CH,+H+M CH, 90, G=CH, + HO, CH,+H*CHJ+H2 CH,-t-O*CH,+OH CH, +OH * CM, + H,O CH,O+OHF~?: HCO+H,O CH,O+Hz=HCO+H, CH,O+M=HCO+H+M CH,O+O= HC0-kOH HCO+OH=CO+H,O HCO+M+H+CO+M HCO+H=CO+H, HCO+O*OH+CO HCO + 0, =HO,+CO CO+O+M*COZ+M CO+OH=CO,+H cO+o,~co,+o CH, +O, * CH,O+O CH,O+M e CH,O+H+M

l.OOE + 17 7.9OE+ 13 2.20E + 04 1.60E + 06 1.60E + 06 7.53E + 12 3.31E+ 14 3.31E+ 16 1.81E+ 13 5.OOE+ 12 1.60E + 14 4.OOE+ 13 l.OOE+ 13 3.OOE+ 12 3.20E+ 13 1.51E+07 1.60E+ 13 7.OOE+ 12 2.4OE+ 13

Units are moles, cubic

E

P 0.000 0.000 3.000 2.360 2.100 0.000 0.000 0.000 0.000 0.000 0.000 O.OUO 0.000 0.000 0.000 1.300 0.000 0.00~ 0.000

86000. 56000. 8750. 7400. 2460. 167. 10500. 81000. 3082. 0. 14700. 0. 0. 0. - 4200. - 758. 41000. 25652. 28812.

K Giovangigli, M.D. Smooke / Adaptive continuation algorithms Table 2 (continued) Methane-air reaction mechanism rate coefficients centimeters, seconds, Kelvins and calories/mole Reaction 20. CH,O+H G=CH,O+H, 21. CH,O+OH * CH,@+H,O 22. CH,O+O+ CH,O+OH 23. CH,O+O, ;;i CH,O+HO, 24. CH, + 02 + CH,O+OH 25. CH,+O*CH,O+H 26. CHs +OH t CH,O+H, 27. HO2 + CO + CO2 + OH 28. Hz +02 it 20H 29. OH+H pH20+H 30. H+O,+OH+O 31. O+H ,+OH+H 32. H+O,+M+HO,-tMa 33. H+O,+Op=HOz+O, 34. H+02+N2+HOz+NZ 35. OH+HO, G=H,O+O, 36. H+HO,+20H 37. O+HO,+O,+OH 38. 20H+O+H,O 39. H,+M+H+H+Mb 40. O,+M*O+O+M 41. H+OH+M+H,O+MC 42. H+HO,+HH,+O, 43. HO2 + HOz * H,02 +O, 44. H,O,+M+OH+OH+M 45. H20,+H~H0,+H, 46. H,O, +OH+ H,O+HO, __-._ a Third oody efficiencies: kSz(H20) = 2lk,,(Ar), b Third body efficiencies: k,,(H,O) = 68&4r), ’ Third body efficiency: k,,(H,O) = 20k,,(Ar).

in the form k, = ATP exp(-

A 2.OOE+ 13 l.OOE+ 13 l.OOE+J3 6.30E + 10 5.20E+ 13 6.80E + 13 7.50E+ 12 5.80E + 13 1.70E+ 13 l.l7E+09 5.13E+ 16 1.80E + 10 2.10E + 18 6.70E + 19 6.70E + 19 5.OOE+ 13 2.50E + 14 4.80E+ 13 6.OOE+ 08 2.23E + 12 1.85E+ 11 7.50E + 23 2.50E + 13 2.OOE+ 12 1.30E + 17 1.60E+ 12 l.OOE+ 13

E,/RT).

323

Units are moles, cubic

E

P

-

-

O.OO!I 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.300 0.816 1.ooo 1.000 1.420 1.420 .300 0.000 0.000 1.300 0.500 0.500 2.600 0.000 0.000 0.000 0.000 0.000

k,,(H*) = 3.3k,,(Ar), kj2(N2) = k,,(O,) k&H) = 2k,,(Ar), k&HZ) = 3k,,(Ar).

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

= 0.

illustrated in Fig. 2 are in excellent qualitative agreement with the asymptotic and Williams [38].

results of Libby

6.2.2. Methane-air flames In addition to the hydrogen-air flames, we applied the adaptive continuation procedure to a stoichiometric and two lean methane-air flames. Strain rate calculations were performed for flames containing 6.5% 6.85% and 9.5% (mole fraction) methane. The equivalence ratios of these flames are + = 0.663, Q)= 0.703 and Q)= 1.0, respectively. The Lewis number of the deficient reactant (methane) in these flames is approximately 0.96. As a result, we anticipate the peak temperature to decrease monotonically for increasing values of the strain rate. For the initial value of the strain rate, the adaptive boundary value method was again used to obtain a solution on a computational domain of one cm. After the initial strain rate calculation was completed, the Euler-Newton adaptive continuation procedure was implemented to obtain

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

324

Fig. 2. C-shaped extinction cdrves illustrating the maximum tempetature versus the inverse of the strain rate for the 8.4% and 9.3% (mole fractior \ hydrogen-air flames. In the absence of Hopf bifurcation, the upper branch represents physical solutions and the lower branch unphysical solutions.

solutions (physical and unphysical) for both increasing and decreasing values of the strain rate. Between 60-90 adaptively chosen grid points were used in each calculation with some grid noise in the solutions. The same normalization was used in the metilane flames as in the hydrogen system. In Fig. 3 we plot the maximum temperature versus the inverse of the strain rate for each of the methane-air flames. Each C-shaped extinction curve was produced from 5: -162 continuation steps. The flames extinguished at uexc= 593 set-’ for the 6.5% flame, aext = 747 set-’ for the 6.85% flame, and a,,, = 1363sec-’ for the 9.5% flame.. Unlike the hydrogen-air extinction curves, we see that. as the strain rate is increased, the peak temperature gradually decreases. No increase

---

----

Fig. 3 C-shaped extinction curves illustrating the maximum temperature versus the inverse of the strain rate for the 6.5% 6.85% and 9.5% (mole fraction) methane-air flames. In the absence of Hopf bifurcation, the upper branch represents physical solutions and the lower branch unphysical solutions.

Fig. G. Comparison of the calculated (solid line) distance of the flames from the symmetry plane versus the jet velocity with the experimental measurements of Sato.

V. Giovangigli, M.D. Smooke / Adaptive continfration algorithms

325

in the peak temperature was observed for calculations with strain rates as low ;s a = 75 set-‘. The maximum temperatures in the three calculations were 1780 K for the 6.5% flame, 1846 K for the 6.85% flame and 2196 K for the 9.5% flame. These values compare favorably with the adiabatic flame temperatures of the corresponding freely propagating flames. In particular, we calculated adiabatic flame temperatures of 1780K, 1846 K and 2232 K for the 6.58, 6.85% and 9.5% flames, respectively. In all three cases, as the value of dZ’,,,/d( l/a) + 60, the flames extinguish. We can, however, continue past the extinction point with the arclength procedure. We find that, as the strain rate begins to decresease, the peak temperature continues to fall and again, in the absence of Hopf bifurcation, the upper branch represents physical solutions and the lower branch unphysical solutions. In Fig. 4 we plot the distance of the flames from the symmetry plane versus the velocity of the input jets. We observe that good agreement is obtained between our calculations and the experimental measurements of Sato [49]. Sato determined the separation distance by observing photographically the center of the luminous (CH emission) flame zone. In our calculations we define the flame separation distance as the distance from the plane of symmetry to the point at which the temperature is equal to 50% of its maximum value. In methane-air calculations with C, chemistry we have found that the luminous zone occurs at temperature values higher than 50% of the peak value. As a result, our separation distances are somewhat larger than those measured by Sate. Hence, our calculations and Sato’s measurements are actually in better agreement than those illustrated in Fig. 4. The small discrepancy in the extinction velocities is due in part to flow field perturbations that arise from the fact that a thermocouple probe was inserted into the flame when the distances were being determined. In addition, the outer flow model we have used at the edge of the boundary layer (u, = ax, uC= -ay) is not strictly the same as that in Sato’s experiments. 6.3. LRan and rich extinction Extinction of counterflow premixed ltinar flames can also be obtained (for a given strain rate) by adjusting the equivalence ratio. In particular, experimental results show that as the equivalence ratio is either increased or decreased past a critical point, the flame extinguishes. Starting with a hydrogen-air flame having an equivalence ratio Q)= 0.245 (9.3% mole fraction) and a strain rate of a - 10OO~ec-~ which was obtained during the strain rate extinction calculations (see Section 6.2), the Euler-Newton adaptive continuation procedure was used to obtain solutions (physical and unphysical) for all values of the equivalence ratio. Each solution contained between 80-120 adaptively chosen grid points. The normalization (5.5) was choosen with g* and I as described in (5.8). In Fig. 5 we illustrate the temperature profile for a hydrogen-air flame with + = 1.O (29.6% mole fraction). We observe the sharp rise of the temperature in the reaction zone with the peak value occurring at the plane of symmetry. The normal velocity is illustrated in Fig. 6. We observe immediately the sharp double hump profile that is characteristic of these flames. The dramatic velocity swing is a result of the strong density variations (see Fig. 7) along with the reduction of the normal component of the mass flow as the plane of symmetry is approached. In Figs. 8 and 9, we plot the mole fraction profiles of the chemical species in the hydrogen-air system. We observe that the radicals HO, and H,O, peak before the flame front (cf. equation (5.7))

V. Giovangigli. M.D. Smooke / Adaptive continuation algorithms

326

I.

000

‘*

0.10

‘*

070

‘.

0.30

’

*

040

’

05C

y

-

’

060

*

‘*‘*‘-

070

080

090

0 00 0 10

0 20 0 XI

(CN

0.40

OS0 y

Fig. 5. Temperature profile in K for a stoichiometric (29.6% mole fraction) hydrogen-air flame with a strain rate of a = lOOOsec-‘.

0.60

0.7c,

C.8C

PJo

Fig. 6. Normal velocity profile in cm/set for a stoichiometric (29.6% mole fraction) hydrogen-air flame with a strain rate of a = 1000sec-l.

Starting from stoichiometric conditions (+ = 1) and then proceeding in the lean direction 1). we anticipate that the peak flame temperature will be reduced gradually. In addition, as the maximum temperature is lowered and the corresponding adiabatic flame speed of an unstrained (Q = 0) flame is reduced, we anticipate that the flame will move closer to the piane of symmetry. Ultimately, as the fuel to air ratio is lowered below a critical value, radical production in the flame will be severely restricted and the flame will extinguish (lean extinctionj. The arclength continuation procedure will then generate unphysical solutions for additional continuation steps until a maximum value of the equivalence ratio is reached. At this point (rich extinction) physical solutions will again be generated and the maximum flame temperature will begin to rise and the distance of the flame from the plane of symmetry will increase as well. The adaptive continuation method was able to generate closed response curves illustrating these phenomena. In particular, in Fig. 10 we illustrate the maximum temperature versus the

(C/B<

0 30100 OOOOBOs

iz 55 2-

omo-

3g

oaoO40-

s 000020-

0.00

1. 0.10

I .I 0.20 0.30

.I

.,.I

0.40 y

Fig. 7. Density in g/cm’

for a stoichiometric

.I.,.,

0.50

060

0.70

.

0.80

0.90

lo

(CM)

(29.6% mole fraction) hydrogen-air a =lOOOsec-‘.

flame with a strain rate of

V. Giovangigli, M.D. Srnooke / Adaptive continuation algorithm

327

) I.

I

*

1

*

I

-

I

*

’

-

1%

0.00 0.10 0.20 0.30 0.40 0.50 0.60 y

C Xl 0.10 0 20 0.30 0.40 0 50 0 60

(CM)

y

Fig. 8. Major and minor species profiles for a stoichiometric (29.6% mole fraction) hydrogen-air flame with a strain rate of a = lOOOsec-‘.

080

10

@I

Fig. 9. Minor species profiles for a stoichiometric (29.6% mole fraction) hydrogen-air flame with a strain rate of a = lOOOsec-‘.

equivalence ratio for a hydrogen-air system with a = 1000 set- ‘. We observe a 1500 K temperature variation among all the flames computed. The lean extinction occurred at + = 0.227 with T = 117 7 K and the rich extinction at + = 4.6 with T = 1256 K. Variations in the peak temperature occur over a much smaller range of equivalence ratios for lean systems than under rich conditions. In Fig. 11 we illustrate the distance of the flame from the plane of symmetry. The flame can shift its position by almost 0.9 cm from lean to slightly rich conditions and then back again to very rich conditions. In particular, the position of the flame from the symmetry plane was equal to 0.03556 cm at the lean extinction limit, 0.8639 cm at its furthest point and 0.1073 cm at the rich extinction limit. A total of 2500 continuation steps were used in generating the closed curves in Figs. 10 and 11. In all cases the adaptive mesh was able to maintain resolution in the flame zone as the flame’s position adjusted accordingly.

z

E r+

1,600~ 1,400. 1,200 1,000

soot

I

I

I

I

I

I

I

I

I

I

0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Equivalence ratio Fig. 10. Extinction curve illustrating

the maximum temperature versus the equivalence with a strain rate of a = IOOOsec-‘.

ratio for hydrogen-air

flames

328

V. Giovangigii, M.D. Smooke / Adaptive continuation algorithms

g

2 +m 2 sE =

‘.OO( 0.60 0.80 t

E

s

2

0.40-

E _z

0.;3-

z 5 % iss

0.00I

f

f

f

I

I

I

f

f

0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5 IO

Equivalence ratio Fig. 11. Distance of the flame from the symmetry plane versus the equivalence ratro for hydrogen-air flames with a strain rzte of a = lOOOsec_‘.

7. Conclusion We have discussed the numerical solution of parameterized two-poiat boundary value problems by a combination of phase-space, pseudo-L rclength continuation methods, Newton-like iterations and global adaptive gridding. The coupling between these techniques has been examined. The numerical method was applied in rhe solution of counterflow premixed laminar flames. In particular, we predicted both strain rate and equivalence ratio extinction for a hydrogen-air system. The methodology is, however, quite general and could be applied to multiparameterized systems, including fold calculations, critical point sensitivity analysis, and multi-dimensional models.

Acknowledgment The research was supported in part by Direction des Recherches et Etudes Techniques and the Air Force Office of Scientific Research. We would like to thank the Scientific Council of the Centre de Calcul Vectoriel pour la Recherche for providing access to a CRAY-1 computer and the Centre Mer Regional de Calcul Electronique for providing computational resources. In addition, we would like to thank Professors S. Candel and G. Duvaut for their support and helpful discussions.

References [ l] J.P. Abbott, An efficient algorithm for the determination of certain bifurcation points, J. Comput. Appl. Math. 4 (1978) 19-27. [2) R.E. Bank and T. Chan, PLTMGC: A multigrid continuation program for parameterized nonlinear elliptic systems, SIAM J. Sci. Statist. Comput 7 (1986) 540.

V. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

329

131 K.N.C. Bray, PA. Libby and J.B. MOSS, Flamelet crossing frequencies and mean reaction rates in premix& turbulent combustion, Comb. Sci. Tech. 41 (1984) 143. [4] J. Buckmaster, The quenching of a deflagration wave held in front of a bluff body, in: proceed&s Seventeenrh Symposium (International) on Combustion (Reinhold, New York, 1979) 835. [5] T. Ghan, Techniques for large sparse systems arising from continuation methods, in: T. Kupper, H. Mittebmn and H. Weber, eds., Numerical Methods for Bifurcation Problems (Birkh;iuser, Base), 1984) 116. 161 T. Chan, On the existence and computation of LU factorization with small pivots, Math. Camp. 42 (1984) 535. [7] T. Chan, Deflation techniques and block elimination algorithms for solving bordered singular systems, SIAM J’. Sci. Statist. Comput. 5 (1984) 121. [8] T. Chm, Newton like pseudo arclength methods for computing simple turning points, SIAM J. Sci. Statist. Comput. 5 (1984) 135. [9] T. Chan and Y. Saad, iterative methods for solving bordered systems with application to continuation methods, SIAM J. Sci. Statist. Comput. 6 (1985) 438. lo] N. Darabiha, S. Candel and F.E. Marble, Numerical calculations of strained premixed laminar flames, Comb. Flame 64 (1985) 203. 11) P. Deuflhard, A stepsize control for continuation methods and its special application to multiple shooting, Numer. Math. 33 (1979) 115. [12] P. Deuflhard, B. Fiedler and P. Kunkel, Numerical path following beyond critical points in ODE models, in: P. Deuflhard and B. Engquist, eds., Large Scale Scientific Computing (BirkhPuser, Boston, MA, 1987) 97. [13) P. Deuflhard, B. Fiedler and P. Kunkel, Efficient numerical path following beyond critical points, SIAM 3. Numer. Anal. (to appear). [14] P. Deuflhard, H.J. Pesh and P. Rentrop, A modified continuation method for the numerical solution of nonlinear two-point boundary value problems by shooting techniques, Numer. Math. 26 (1976) 327. [15] E.J. Doedel and J.P. Kemevez, Software for continuation problems in ordinary differential equations, Tech. Rept., Applied Mathematics Caltech (1985). [16] V. Giovangigli and N. Darabiha, Vector computers and complex chemistry combustion, in: rroceedings Conference on Mathematical Modeling in Combustion, Lyon, France, NATO ASI Series (1987). [17] V. Giovangigli and M.D. Smooke, Calculation of extinction limits for premixed laminar flames in a stagnation point flow, J. Comput Phys. 68 j1987) 327. (181 V. Giovangigli and M.D. Smooke, Extinction of strained premixed laminar flames with complex chemistry, Comb. Sci. Tech. 53 (1987) 23. [19] R. Glowinski, H.B. Keller and L. Rheinhart, Continuation-conjugate gradient methods for the least squares solwtinn of nonlinear boundary value problems, SIAM J. Sci. Statist. Comput. 6 (1985) 793. [20] KX. Heijer and W.C. Rheinboldt, On steplength algofithms for a class of continuation methods, SIAM J. Numer. Anal. 18 (1981) 925. [2X] R.F. Heinemann, K.A. Overholser and G.W. Reddien, Multiplicity and stability of premixed laminar flames: An application of bifurcation theory, Chem. Engrg. Sci. 34 (1979) 833. [2i;] R.F. Heinemann, K.A. Overholser and G.W. Reddien, Multiplicity and stability of the hydrogen-oxygen-nitrogen flame: The influence of chemical pathways and kinetics on transitions between steady states, AIChE J. 25 (1980) 725. [23] J.M. Hyman and M.J. Naughton, Static rezone methods for tensor-product grids, LOS Alan~~ Rept. LA-UR-833245, Los Alamos, NM (1983). (241 S. Ishizuka and C.K. Law, An experimental study on extinction and stability of stretched premixed flames, in: Froceedings Nineteenth Symposium (International) on Combustion (Reinhold, New York, 1982) 327. 1251 A. Jepson and A. Spence, Singular points and their computations, in.- T. Kupper, H. Mittelmann and H. Weber, eds., Numerical Method for Bifurcation Problems (Birkhsuser, Basel, 1984). [26] J. Kautsky and N.K. Nichols Equidistributing meshps with constraints, SI,IM J. Sci. Statist. Compu:. 1 (1980) 499. A general-purpose, transpor table, Fortran chemical [27] R.J. Kee, J.A. Miller and T.H. Jefferson, CHEMKIN: kinetics code package, Sandia National Laboratories Rept. SAND80-8003, Livermore, CA (1980). [28] R.J. Kee, J. Wamatz and J.A. Miller, A Fortran computer code package for the evaluation of W-Phase viscosities, conductivities, and diffusion coefficients, Sandia National Laboratories Rept. SAND83-8209, Livermore, CA (1983).

330

If. Giovangigli, M.D. Smooke / Adaptive continuation algorithms

1291 H.B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, in: P. Rabinowitz, ed., App/ications of Bifurcation Theo7 (Academic Press, New York, 1977) 359. (301 H.B. Keller, Practical procedures in path following near limit points, in: R. Glowinski and J.L. Lions, eds., computing M&o& in Appiied Science and Engineering (North-Holland, Amsterdam, 1982) 177. [3I] H-B- Keller, me brrdering algorithm and path following near singular points of higher dlity. SIAM J. SCi. Statist. Comput. 4 (1983) 573. 1321 H.B. Keller and R.K.H. Szeto, Calculation of flows between rotating disks, in: R. Glowinski and J.L. Lions, eds., Computing Methods in Applied Science and Engineering (North-Holland, Amsterdam, 1980) 51. (33) M_ Kubicek, Dependence of solutions of nonlinear systems on a parameter, ACM Trans. Math. SofW- 2 (1976) 98. (341 B. Larrouturou, Utilisation de mailiages adaptatifs pour la simulation de flammes monodimensionnelles instationnaires, in: R. Glowinski, B. Larrouturou and R. ‘Temam, eds., Numerical Simulation of Combustion Phenomena (Springer, Heidelberg, 1985) 300. (351 M. Lentini and V. Pereyra, An adaptive finite difference solver for nonlinar two-point boundary value problems with mild boundary layers, SIAM J. Numer. Anal. 14 (1977) 91. [36] P.A. Libby and K.N.C. Bray, Implications of the laminar flamelet mode! in premixed turbulent combustion, Comb. Ffame 39 (1980) 33. (371 P.A. Libby and F.A. Williams, Structure of laminar flamelets in premixed turbulent flames, Comb. Flame 44 (1982) 287. 1381 P.A. Libby and F.A. Williams, Strained premixed laminar flames with two reaction zones, Comb. Sci. Tech. 37 (1984) 221. 139) J.A. Miller, R.J. Kee, M.D. Srrooke and ;.F. Grcar. The computation of the structure and extinction limit of a methane-air stagnation point diffusion flame, Paper #WSS/CI 84-10, presented at the 1984 Spring Meeting of the Western States Section of thz Combustion Institute, University of Colorado, Boulder, CO (1984). 1401 R.G. Melhem and W.C. Rheinboldt, A comparison of methods for determining turning points of nonlinear equations, Computing 29 (1982) 201. [41] C.E. Pearson, On a differential equation of boundary layer type, J. Math. Phys. 47 (1968) 134. [42] V. Pereyra and E.G. Sewell, Mesh selection for discrete solution of boundary value problems in ordinary differential equations, Numer. Math. 23 (1975) 261. [43] J.F. Rameau and Cl. Schmidt-Lain& Numerical bifurcation in chambered diffusion flames, in: R. Glowinski, B. Larrouturou and R. Ternam, eds., Numerical Simulation of Combustion Phenomena (Springer, Heidelberg, 1985) 335. 1441 W.C. Rheinboldt, Solution fields of nonlinear equations and continuation methods, SlAM J. Xumer. Anal. 17 (1980) 221. f45) W-C. Rheinboldt, Computation oi critical boundaries on equifibrium manifolds, SIAM J. Numer. Anal. 19 (1982) 653. [46] W-C. Rheinboldt, Numerical analysis of continuation methods for nonlinear structural problems, Cornput. Structures 13 (1981) 103. [47] Y- Reuven, M.D. Smooke and I-I. Rabitz, Sensitivity analysis of boundary value problems: Application to nonlinear reaction-diffusion systems, J. Comput. Phys. 64 (1986) 27. [48] RD. Russell and J. Christiansen, Adaptive mesh selection for solving boundary value problems, SIAM J. Numer. Anal- 15 (1978) 59. (491 J- Sate, Effects of Lewis number on extinction behavior of premixed flames in a stagnation flow, in: Proceedings N-+wteenth Symposium (%nternationa!) on Combustion (Reinhold, New York, 1982) 1541. [5O] J- Sato and H. Tsuji, Behavior and extinction of a premixed flame in a stagnation flow for general Lewis numbers, in: Pr0ceeding.s Sixteenth Japanese Symposium on Comoustion (1978) 13. [5I] J. Sato and I-I. Tsuji, Extinction of premixed flames in a stagnation flow considering general Lewis number, Comb. Sci. Tech. 33 (1983) 193. 1521 G-I. Sivashinsky, On a distorted flame front as a hydrodynamic discontinuity, Acta Astronautica 3 (1976) 889. (531 M-D. Smooke, Solution of burner-stabilized premixed laminar flames by boundary value methods, J. Comput. Phys. 48 (1982) 72. 1541 M-D. Smooke, On the use of adaptive grids in premixed combustion, AIChE J. 32 (1986) 1233. 1551 M-D- Smooke and M-L. Koszykowski, Fully adaptive solutions of one-dimensional mixed initial-boundary value problems with applications to unstable problems in combustion, SI.4 M J. Sci. Statist. Cornput. 7 (1986) 301.

K Giovangigli, M.D. Smooke / Adaptive continuation algorithms

331

1561 M 3. Smooke and M.L. Koszykowski, Two-dimensional fully adaptive solutions of solid-solid alloying reactions, 9. Comput. Phys. 62 (1986) 1. (571 M.D. Smooke and R.M.M. Mattheij, On the solution of nonlinear two-point boundary value problems on successively refined grids, Appl. Numer. Math. 1 (1985) 463. [58] M.D. Smooke, Y. Reuven, H. Rabitz and F.L. Dryer, Application of sensitivity analysis of premixed hydrogen-au flames, Comb. Sci. Tech. 59 (1988) 295. [59] M.D. Smooke, J.A. Miller and R.J. Kee, Determination of adiabatic flame speeds by boundary value methods, Comb. Sci. Tech. 34 (1983) 79. f6Q) A.B. White, On selection of equidistributing meshes for two-point boundary value problems, SIAM J. Numer. Anal. 16 (1979) 472. [al] A.B. White, On the numerical solution of initial/boundary value problems in one space dimension, SIAM J. Numer. Anal. 19 (1982) 683.