COMPUTATIONS
OF LAMINAR FLAME EXPLICIT NUMERICAL
PROPAGATION METHOD
USING
AN
R O L F D. R E I T Z Research Scientist, Courant Institute of Mathematical Sciences, New York University, New York, N.Y. 10012
Numerical methods are applied to a coupled system of onedimensional unsteady reactiondiffusion equations to seek propagating wave solutions. These equations model flame propagation in certain combustion systems w h e n constant pressure combustion is assumed, with Lewis number unity, and when a Lagrangian coordinate transformation is introduced. In the numerical integration schemes the reaction terms are computed noniteratively, using two different secondorder accurate methods. In the first, the numerical timestep is limited by the Lipschitz timestep constraint as is the case with standard explicit schemes. The second scheme is constructed in such a way that this timestcp limitation is not present and computations are made with timesteps which are orders of magnitude larger than those possible with standard methods. In both schemes the diffusion terms are differenced using the unconditionally stable first order accurate explicit method introduced by Saul'yev. The leading truncation error of this method is a numerical dispersion error which has been found to reduce accuracy less than the numerical diffusion error which arises in more standard explicit methods. These numerical schemes are applied to an ozone decomposition flame model problem. An adaptive grid procedure is also implemented for further increases in computational efficiency and the results show good agreement with the fourth order accurate results of Margolis. This indicates that efficient lower order numerical methods can be sufficiently accurate for practical computations of combustion.
Introduction Detailed mathematical m o d e l s are now b e i n g applied to practical problems in combustion. Attention is being given to models w h i c h include comprehensive elementary chemical kinetic a n d turbulence mechanisms. This trend has increased the n u m b e r of equations to be solved and has stimulated research in numerical methods for the conservation equations of reacting gases. In this paper we discuss numerical results for onedimensional unconfined deflagration wave propagation. The problem is formulated in such a way that the momentum and overall mass conservation equations need not be considered. This simplification focuses attention on the treatment of the important reaction rate and diffusion terms in the governing equations. The presence of the reaction terms can cause timestep restrictions in computations of combustion. Each chemical species evolves with its own characteristic time scale, and the numerical timestep may
be required to be comparable to the smallest of these time scales for reasons of stability a n d / o r accuracy. This problem has led many researchers to use implicit solution techniques in calculations of stiffly coupled combustion equations. Reactiondiffusion equations have also been solved following this approach. This includes methods in which a system of coupled ordinary differential equations (O.D.E.) are solved using stiff O.D.E. routines (methods of lines, t'* operator splitting3); and methods in w h i c h block tridiagonal matrices which result from an implicit treatment of the diffusion terms, are computed. 4 In this work we formulate and apply two explicit methods to a onedimensional propagating flame calculation. The motivation for the use of these methods is the potential for increases in computational efficiency a n d accuracy. These benefits are important considerations in problems which involve a large n u m b e r of participating chemical species, due to the n u m b e r of equations to be solved. The numerical method used for the diffusion terms 433
434
PREMIXED FLAME STUDIES
does not have a diffusional stability restriction and can be extended for use in multidimensional computations. For the reaction terms, we present a new stifflystable noniterative method which allows large timesteps to be used. The construction of this method exploits the physically meaningful requirement that the mass fraction of the chemical species and the gas temperature be nonnegative. Further increases in computational efficiency are obtained through the introduction of an adaptive grid procedure which concentrates numerical grid points within the flame. Theoretical results are obtained using these numerical methods for the propagating ozone decomposition flame. This model problem has been studied in detail by other authors,~'~'~ and thus results are available for comparison with the present results.
Model Formulation The governing equations for onedimensional flame propagation can be written as follows: e
a(pu)
ap + #t
a(pu)
a(pu ~)
#t
ax
a(ph.)
a(puh.)
at
o~x
= o
ax
+
ax
+
p
(1)

(2)

ax
 ~ h (k)
T,p, P, u, ~,"'
T.~,
X=  ~
at
ax
+to
u.~. k~ XI
oo
FIG. 1. Schematic diagram of onedimensional flame propagation problem.
equations for the N chemical species having mass fractions Y(*). h, is the total mixture enthalpy and in Eq. (5), the thermal equation of state, it is assumed that the mixture is a perfect gas. Fiek's law has been used to describe species diffusion with species diffusion coefficients .~k>. Thermal radiation, body forces, Soret, Dufour, and pressure gradient diffusion, and bulk viscosity effects have been assumed to be negligible. In unconfined flame propagation the combustion wave separates the unburned gases at x = oo and the burned gases at x = o0 as depicted in Fig. 1. The unburned mixture has constant density Po, temperature T Oand species mass fractions Y~) and velocity Uo = 0. In the burned gases p = p~, T = T~ and Y(*) = Y~). However, the quantities at x = 0o and the wave speed, S, are unknown in general and are to be computed. The governing equations are simplified by assuming that p = const., which is a good approximation for low Math number flames, e This allows the momentum equation (2) to be dropped, and a further consistent assumption is to neglect the kinetic energy (  u 2 / 2 ) and its derivatives in Eqs. (3) and (6). The convective terms in Eqs. (3) and (4) are eliminated by defining i" = t and
(3)
a~
p
  =
. ~ (k)
UNBURNED gASES (PREMIXEO)
P $
Opy(k) _
(FLAME)
BUR NEO EASES
ax _ _
REAOTION ZONE
j,~/
ax
(*),
a~ ,
po
.
.
p .
at
.
po
u
(7)
ax k = 1..... N (4)
p = p~ T E
Y~)/MW(~)
This transformation satisfies Eq. (1) identically and Eq. (4) becomes
(5)
= 
aE
kffil
.~(*)
as
+  
a~
p
(8)
where The energy equation (3) is written in terms of temperature by using Eqs. (8) and (6) and, when it is assumed that for all k, c (k~ = %, ~(k) = .~ where % and.~are constants wit~ p.~ = h/% (Lewis number = 1), Eq. (3) becomes
N
h, = E Y(k)h(k)+ u2/2 kI
E y~k) k  1
(*)dT' + h(,k) + u'/2
Cp
(6)
N
Tr
aT = __a[ ( L ) 2 ~ aT l The first four equations are the overall mass, momentum, energy and the N species conservation
EtO~k)k, h~k) (9)
a{
a~
a~ J
p%
435
COMPUTATIONS OF LAMINAR FLAME PROPAGATION Equations (8) and (9) are supplemented by the equation of state (5) and the boundary conditions which specify zero diffusion of heat or mass at +w, i.e. aY ~k~
OT =0;
0s
=0
at
i=_+oo
(10)
0s
The reaction mechanism for the ozone decomposition flame model problem involves atomic oxygen (O), molecular oxygen (O2) and ozone (03) through the reactions
alized with To, Po, to and poTo%, respectively, t o is a characteristic time scale chosen as 5.878 • 10 5 see. and a masslength scale xo = 5.048 • 10e g / cm 2 is also introduced. The computational domain is 0 < x < 50 and the initial conditions were prescribed using the data of Margolisfl The problem which was studied models flame propagation in a mixture of 2 / 3 0 2 , 1/3 O3 (wt.) with T O = 300*K and Po = 0.821 atm. In summary then, Eqs. (8) and (9) are written in the generalized form
~U(k)
k[ 03+M~O2+O+M k~
k,' 0+0~0~+0~ k: k: 02+M~O+O+M
k: M stands for any of the three chemical species (hence there are 7 reactions) and the superscripts, f and b, refer to the forward and backward reactions. The reaction rates are given by 7
r176 = E
MW'k' (vk, . .  v,.,)" {Rim R~} (11)
where
R'. = k'.
~=, \ M W ~") ]
and R~  0 is given correspondingly with f replaced by b and the stoichiometric coefficients v',.~ replaced by v','.~. The rate constants are given by the Arrhenius expression k~ = B~ T *~ e es with a corresponding expression for k~. The constants B,,, Sm and E , and other relevant constants are given in Margolis. 2 The dimensionless governing model reactiondiffusion equations are as given in Eqs. (8) and (9) but with p 2 ~ assumed equal to a constant = d = x~o/to, and x = po~/Xo, t = ?/to and N = 3. yr y~2~ and y~a~ are the mass fractions of O, 02 and 03, and T, p, r and h~k~ are nondimension
~2Utk)
= d Ot
+ F ..... u~X~);
(12)
ax 2
k = 1, 2 . . . . . K. Here the F
Numerical Method Two finite difference schemes were used for the numerical integration of Eq. (12). In both schemes the diffusion terms are differenced using the unconditionally stable first order accurate explicit method first introduced by Saul'yev. 7 Variations of this method which allow the order of accuracy to be increased are given in Reitz. s In the present case the leading truncation error of the method can be ~k~ w here r = d A t / A x 2, shown to be of the form rAtu,~, At is the numerical timestep, Ax is the mesh spacing and the subscripts denote differentiation. This error corresponds to a numerical dispersion error in steady wave propagation. In scheme I, the reaction terms are computed using a second order accurate explicit extrapolation formula. In scheme II the reaction terms use a variation of this method which maintains the same order of accuracy, and is stifflystable. The numerical schemes are summarized in Eqs. (13) below. The basic method consists of two predictor sweeps of the mesh (in opposite directions) in Eqs. (13a,b), followed by a correetor step in Eq. (13e). The finite difference approximation to the continuum solution U~.~ = u a~ (x, nat) is computed at grid points x~ from the set of equations V,  U~.,
2d(1 + ~/)
At
hlh2h
{h 1 (U~t+,  U~,.,)
 h2(V,  V,_I)} + a G ~  13G~' (13a) w,

uL
At
2a(1 +~) =
h',h'2h'
{h',(w,+,
 h'~ ( U ~ , ,  U ~ , , , _ , ) }
 w,)
+ etG~  1 3 G ~  '
v ; . ; ~ = (v, + w , ) / 2
(13b)
(13c)
436
PREMIXED FLAME STUDIES
)
where
h,
= x j  x~_,;
h~ = x i + l  x j ;
h'1 = x ,  x,_l;
h = h~+h 2
h l = x , + ~  x , ; h' = h ' l + h ' 2
j = 1,2 ..... M;
i = Mi+I
k = 1,2 ..... K;
l = 1,2,...,K
where uo is u(t = 0), while the method in scheme II leads to the difference formula n + l

un
u"(1  u "§ and M is the number of computational points in the domain. In scheme I the parameter /in the diffusion term (the first term on the right hand side of Eqs. (13a,b)) is set equal to zero and the reaction variables a, 13, G'~ and G~ ~ in Eq. (13a) are given by
At which can be shown by induction to have the analytical solution u" = i
/(1. 1+
G V ' = F,(UT.;') ct = 3 / 2 ; 13 = 1 / 2
(14a)
(The subscript j in Eq. (14a) is replaced by i for the case of Eq. (13b).) For scheme II, the variables in Eq. (13a) are
F2
.~_
)
.
Uo
G~ = F~(U~.,)
,
(1 + A t )  "
+
n
n
This is readily seen to be a first order accurate approximation to the solution of the differential equation and the method permits arbitrarily large timesteps At to be used. Notice that the construction of the method requires that u > 0 and F(u)/u < oo as u ~ 0. These conditions are met for first or higher order reactions. In the case of scheme I, the Lipschitz timestep constraint const
(U2,)
At __
(16) 0F~ (U, ,)
max I ~ 1 j )+ G~_ ~ = . F k _( U t ,,_~
k,vt.j ,
u~,~'
ct = ~ + 3 / 2 ; 13 = 1 / 2
Vj (14b)
For the case of Eq. (13b) the subscript j in Eqs. (14b) is replaced by i and V is replaced by W. The superscripts + and  in Eqs. (14b) refer to the strictly positive and negative contributions to the overall source term F k, i.e. Fk(V~,,) = F ~ + (u,,,)   F/(U L) +
where F k and F~ are each > 0 and are each made up of terms which are only positive (compare Eq. (11)). The use of scheme II is illustrated by applying it to a single nonlinear O.D.E.
tin(t)
= F(u) = F§
 F(u)
(15)
dt We consider, for example, the case F(u) = u  u 2 with 0 < u and take d = 0, ct = 1 and 13 = 0 in Eqs. (13a,b). This corresponds to a first order accurate version of the method. In this case Eq. (15) has the exact solution
limits the timestep of the method. 9 In practice, the computation of the dominant eigenvalue in the Jacobian matrix of Eq. (16) was avoided by requiring that the species mass fractions remained positive. If negative mass fractions were found, the timestep was halved and the computational cycle was repeated; after 20 such steps the timestep was doubled. The adaptive grid proeedure is similar to that used by Lund, 4 but is simpler to implement since the nonuniform mesh is prescribed analytically. A new mesh was generated when the numerical solution changed beyond a prescribed tolerance, and the solution variables were interpolated onto the updated mesh. The details are as follows. If z(x) is a continuous variable which represents the mesh spacing at point x, a variational problem is formulated which minimizes the integral
i(z) I xM Jx,
z \dx/
(17)
subject to the constraint
I xM dx M' =
xI
;~
(18)
COMPUTATIONS OF LAMINAR FLAME PROPAGATION Minimizing the integral in Eq. (17) minimizes the change in mesh spacing across the domain where x~ and xM are fixed boundary points, while Eq. (18) conserves the number of zones M' = M  1. Equations (17) and (18) lead to an Euler equation which may be solved to give A z(x) =  
(x 
x.,,~ ~ +
Zm,~
which is exact for exponential functions. Equations (21) and (22) were chosen in anticipation of solutions which have approximately exponential decay outside of a narrow flame zone. Finally, mass conservation was enforced on each newly defined mesh by adjusting the species mass fractions at each grid point x~ such that
~I k>
(19)
4Zmi.
yl, k,
Here the constants of integration were specified by requiring that the minimum value of the mesh function, Zm, . , coincides with the point of maximum temperature gradient, x ~ , . , where
] T2

Tt ]
(20)
Zmi n = C
c is a constant which specifies the number of zones in the flame and T~ and T~ are the maximum and minimum temperatures in the domain respectively. The eigenvalue was determined by integrating Eq. (18) with z ( x ) given in (19). The special case A = 0 corresponds to an equally spaced mesh, while A = (2~r/M') 2 corresponds to the case xM * 0% x~
437
N
<23)
kI
where ~is obtained from the interpolation Eqs. (21) and (22). In practice, this correction was negligible in the computations with scheme I. The correction became more significant as the timestep was increased in the computations with scheme II.
Discussion of Results Numerical results for the propagating ozone decomposition flame model problem are shown in Figs. 2 to 6. The results in Figs. 2, 3, 4 and 5 were I
2 ~ 900K o
00. 0.0
The mesh was redefined if the new x .... and z . , . violated either of the inequalities
0.4~~ L
o Xmi n   Zmi n ~ Xmi n ~~ Xmi n "~ Zmi n
0
1,2
zOain
<2
~
,
I
,
I
,
I
I
i
I
~
{3 0
o
z ~ , ~ < 2 z ~
where X~176 ZOm,n defined the mesh at the previous remap. The solution variables (e.g. the temperature T, at x~) were interpolated onto the new mesh xj with the nonlinear formula Tj
I
% 9
I D
D
S >"
T~ + f~T 2
0
1.2
o L
08
'
0
'
(21)
l+f, where T, f j = g,
g,+~; g, 
_

T1
_
T 2  T,
and 8 = (xj  x , ) / ( x , § 1  x,) for x, < x, < x , § and T, # TI, T v This interpolation is exact for a class of nonlinear model reactiondiffusion equations which are presented in Reitz." For T i = T~ (1 = 1,2) we used instead
o
,
0
I0
20
30
,
o
40
50
X
FIG. 2. Profiles of temperature [] and atomic oxygen mass fraction O in the propagating ozone decomposition flame for t = 0, 10, 20 and 30. 15 grid points, c = 4 in Eq. (20) using scheme I.
P R E M I X E D F L A M E STUDIES
438
T
y(l}
i. O0
yt2)
o 0.75
x
Z
2
o
6
4
I
I
8
Io
x FIG. 3. P r o f i l e s o f t e m p e r a t u r e [ ] and a t o m i c oxygen 9 molecular oxygen ~ , ozone A mass fractions
for t = 0.1. 15 grid points, c = 4 in Eq. (20) using scheme L   r e s u l t s of Margolis. 2 1.25
y(I)
1.00
y(2)
0 rO
Z []
'o X
obtained using scheme I. Figure 2 shows the computed temperature and atomic oxygen mass fraction profiles through the flame in four separate graphs, each of which spans the computational domain. The graphs show results which correspond to t = 0, 10, 20 and 30 where the nondimensional time increases from top to bottom in the figure. The adaptive grid used only the 15 computational points w h i c h are shown in each frame with c = 4 in Eq. (20). The initial conditions are such that the ignition process was essentially confined to 0.8 < x ~ 1.2. During the computation the minimum grid spacing z.,,, changed from .07 at t = 0 to about .5 for t > 4. The slope of the oblique line w h i c h connects the reference point in the temperature profile (T = 900~ of each graph is proportional to the flame speed since the steady speed is invariant under the transformation (7). 5 The fact that the line is straight for t > 10 indicates that steady wave propagation has been reached. However, the mass fraction of atomic oxygen at x = 0, which is diminishing in the successive frames, is still not in equilibrium.
Y
O. 7 5
(3)
I
.>. 0"
0,50 o x
O 0.25
A
0
i 20
22
24
26
28
30
oJ 32
X FIG. 4. Profiles of temperature [] and atomic oxygen O, molecular oxygen • ozone /~ mass fractions for t = 35. 15 grid points, c = 4 in Eq. (20) using scheme L   r e s u l t s of Margolis. 2
COMPUTATIONS OF LAMINAR FLAME PROPAGATION 50


. . . . . . . . . . . . . .
~
0
49
LB
47
NO S T E A D STATE ~
46
1,5
I 5
I 10 I
I
I 15 M I
2
3
4
I 20
1 25
I 30
I
I
6
8
C
F1c. 5. Variation of computed steady flame speed with number of grid points M using scheme I. cconstant in Eq. (20), At = 0.005.  wave speed of Margolis. 2 Details of the flame's transient development and its steady state structure are shown in Figs. 3 and 4 respectively. The figures show temperature and atomic oxygen, molecular oxygen and ozone mass fraction profiles in portions of the computational domain. The solid curves show the numerical results obtained by Margolis, 2 who used an equally spaced mesh (Ax = 0.2) and a fourth order accurate method of lines. A comparison of the results shows excellent agreement of the profiles in Fig. 3 at t = 0.1. In Fig. 4 (t = 35) the computed profiles also show good agreement, but the flame location differs due to a difference in the steady wave speed in the two calculations. This difference in flame speed can be seen in Fig. 5 which shows the variation in flame speed as a function of the number of grid points used in five separate computations. Here the zone constant c was varied in proportion with the number of grid points as is indicated in the figure. The flame speed was computed using the method which was indicated in the discussion of Fig. 2 and has the value S = 49.8 + 0.1 c m / s for 30 grid points. This agrees well with S = 49.7 c m / s which was found by Margolis. ~ The first order accurate method of lines of Bledjian' gave 54.3 c m / s (with &x = 0.4). The rate of convergence of the results with an increase in the number of grid points in Fig. 5 indicates that scheme I behaves as if it were second order accurate. This is explained by the fact that the factor r in the first order truncation error of the method is small (less than 0.1) due to the size of the timestep, which is limited by the reaction rate terms in scheme I. The computer execution time
439
varied nearly linearly with the number of grid points; the computation with M = 15 took 7.3 min., CDC 6600, which compares favorably with the 20 min. quoted by Margolis. ~ The time step ranged from At = 0.0025 in the transient to At = 0.005 (i.e. At 0.3 its) for t ~> 4. Time step values were not given by Margolis; ~ Bledjian' gave At < 0.005 and Lund, 4 in a similar ozone model problem, used the value At = 0.007, or At = 0.4 its. In Lund's implicit method the timestep was controlled by specifying a maximum allowable number of iterations for convergence of the solution in each computational cycle. The agreement between these timesteps indicates that explicit methods can be competitive with implicit methods in certain combustion problems. The results which are shown in Fig. 6 summarize flame propagation data obtained using scheme II. The plot shows the error in the computed steady state flame speed for 21 different runs as a function of the numerical timestep which was kept fixed for each run. Here percentage error is defined as (S  49.8) • 100/49.8, where S is the computed flame speed in c m / s for each run. The three curves correspond to 11, 15 and 30 grid points, The initial data was given with the profiles computed using scheme I at t = 10, since the resolution of the transient is not required for steady flame speed results, and the computations were continued up to t = 30. The results in Fig. 6 show that scheme II permits timesteps two orders of magnitude larger than those allowable with scheme I. The computer time with At = 0.5 was correspondingly two orders of magni
50

40
t~ 3O O t~ w  20
10
0
I
I
I
I
I
0.I
0.2
0.3
0.4
0.5
OOO5
At
FIG. 6. Percentage error in computed steady flame speed relative to S = 49.8 c m / s as a function of the timestep At using scheme II. [] O A correspond to 30, 15 and 11 grid points respectively.
440
PREMIXED FLAME STUDIES
tude shorter. The computed flame speed with scheme II is seen to decrease with increasing At. A t the stability limit of scheme I, At = 0.005, the results with scheme II were virtually identical to those shown in Fig. 5 for scheme I. However, the error was about 50% at At = 0.5 and at these large values of At the flame speed was not truly steady but changed by up to 5% during the course of a run. The factor r in the truncation error of scheme II with M = 30 and At = 0.5 is  8. This would help explain the increased error in this case. Notice that this is roughly an order of magnitude larger than the standard diffusional stability limit, r < 1/2. In spite of the relatively large error in flame speed at At = 0.5, the flame temperature and the species mass fractions at x = 0 differed from those with At = 0.0')5 by less than 1%, with the exception of the O atom mass fraction which differed by about 9%. The mass conservation correction in Eq. (23) had a maximum value of 4%. These results show that substantial savings in computer time can be gained with scheme II but with an accompanying decrease in the accuracy of the computations as the timestep is increased.
Conclusions
The results demonstrate that increases in computational efficiency have been gained with the use of the present explicit methods. With the adaptive grid, calculations were made with an order of magnitude fewer grid points than used by Margolis. z The stifflystable scheme II could use timesteps two orders of magnitude larger than that of standard methods, with the diffusion parameter, r, an order of magnitude larger than that usually permitted for numerical stability. The good agreement between the present results and the fourth order accurate results of Margolis, 2 shows that lower order numerical methods can be sufficiently accurate for practical computations. The methods presented here show improved accuracy over the first order results of Bledjian.~ This indicates that a first order numerical dispersion truncation error reduces accuracy less than the diffusive error in his method. This finding confirms the results of other studies in Reitz. s The fact that the timestep of scheme I, which is limited by the dominant eigenvalue of the reaction terms, was comparable to the timesteps used by Lund 4 and Bledjian 1 indicates that explicit methods can be competitive with implicit methods in certain combustion problems. The results show that substantial savings in computer time can be gained with scheme II, but at the expense of a compromise in the accuracy of the computations. In practice the magnitude of the error could be controlled through the choice of the
timestep, and an error control procedure based on the mass conservation correction, Eq. (23), is currently being developed. Scheme II could be expected to perform well in problems which involve trace species which do not influence the energy release and/or flame speed significantly, but which do have short characteristic evolution times, and thus would limit the timestep of more standard methods. Work is in progress on higher order accurate versions of scheme II and the method is also being applied to systems of O.D.E.'s.
Nomenclature
B c ~,d
preexponential constant zone specifying constant in Eq. (20) specific heat at constant pressure species diffusion coefficients E activation energy F source term in Eqs. (12) and (13) h,h~,h~ grid point spacing in Eqs. (13) h specific enthalpy k rate constant K number of reactiondiffusion equations in (12) M number of computational points MW species molecular weight n running index; n a t = t in Eqs. (13) N number of chemical species p pressure universal gas constant r
dAt/Ax 2
s S, S t, t T u
constant in Eq. (11) wave speed in (x, t) and (s t) planes time temperature generalized dependent variable, mass averaged gas velocity numerical approximation to u variables in Eqs. (13) space coordinate in (x, t) and if, i) planes species mass fraction nonuniform mesh spacing = hx for equally spaced mesh.
U V,W x,~ Y z
Greek symbols 7 At hx h A bt V
P (O
parameters in source terms of Eqs. (13) parameter in diffusion terms of Eqs. (13) numerical timestep mesh spacing for uniform mesh thermal conductivity of mixture eigenvalue in Eq. (19) mixture viscosity stoichiometric index in elementary reaction step fluid density mass rate of species production/depletion
COMPUTATIONS O F LAMINAR FLAME PROPAGATION
Subscripts 0 1 i,j k, l, m m r
REFERENCES
quantities in the unburned mixture quantities in the burned mixture running indices 1..... M running indices 1..... K running index for elementary reaction steps reference value
Superscripts b f k n ' "
441
backward elementary reaction step forward elementary reaction step running index 1..... N numerical time level nat = t forward reaction backward reaction
1. BLEDJIAN,L.: Combust. Flame 20, 5 (1973). 2. MARGOLlS,S. B.: J. Comp. Phys. 27, 410 (1978). 3. OTEY, G. R. AND DWYER,H. A.: A.I.A.A.J. 17, 606 (1979). 4. LUND, C. M.: A General Computer Program for
Calculating timedependent Phenomena involving Onedimensional Hydrodynamics, Transport and Detailed Chemical Kinetics, Lawrence Livermore Laboratory Preprint UCRL52504 (1978). 5. HIRSCHFELDER,J. O., CURTlSS,C. F. ANDCAMPBELL, D. E.: J. Phys. Chem. 57, 403 (1953). 6. WILUAMS,F. A.: Combustion Theory, p. 95, AddisonWesley, 1965. 7. SAUL'YEV, V. K.: Integration of Equations of Parabolic Type by the Method of Nets, p. 52, Pergamon Press, 1964. 8. REITZ, R. D.: The Application of an Explicit
Acknowledgments The author would like to thank Professors S. Burstein and E. Isaacson for helpful discussions and comments. This work was supported under D.O.E. contract EY76C023077.
Numerical Method to a ReactionDiffusion System in Combustion, Courant Institute Report COO3077162, 1979. 9. HOFF, D.: Siam J. Numer. Anal. 15, 1161 (1978).
COMMENTS R. E. Mitchell, Sandia National Laboratories, USA. Have you applied this explicit procedure to a realflame problem in which you have variable thermodynamic and transport properties and detailed finiterate chemistry? Would you care to comment on its possible applicability?
schemes do not work in these systems due to extreme stiffness problems [2, 3]. Did you try at any time to simulate more complex flames than the ozone flame considered in your presentation?
REFERENCES
Author's Reply. The method is presently being extended to include variable thermodynamic and transport properties. Encouraging results have been obtained with a model scalar equation in which the diffusion term is of the form (du~) x. The method differs from that given in Eqs. (13) in that the diffusion coefficients dk.j are now evaluated at three neighboring points and two time levels (n and n  1) in the difference scheme.
[1] J. WABNATZ,Berichte der Bunsengescllscha~fur PhysikalischeChemie82, 193 (1978). [2] Ibid. 82, 643 (1978). [3] I b i d 8 2 , 8 3 4 (1978).
Author's Reply. It is planned to include the modeling of more complex flames in future studies as a further test of the method. It is expected that problems due to stiffness, which are commonly associated with explicit methods, will be lessened with the use of the stiffstable scheme II.
1. Warnatz, Institut fuer Physikalische Chemie, West Germany. In connection with implicit calculations of ozone decomposition flames [1] the structure of these flames could be shown to be calculated by explicit methods likewise. Therefore it would be interesting to extend your explicit method to more complex flames including chain branching (e.g. H 2  O~ flames). To my opinion explicit
H. Krier, University of Illinois, USA. It would appear that your assumption of a flame with constant specific heat, Cp, and a constant mass diffusivity, D, and then the assumption of a unity Lewis number,
442
PREMIXED FLAME STUDIES
that y o u are treating a c o n s t a n t d e n s i t y flame. W o u l d y o u e x p a n d on this, if i n d e e d y o u are a l l o w i n g the d e n s i t y to vary in t h e flame, or do y o u let the gas c o n d u c t i v i t y adjust w i t h the d e n s i t y ?
Author's Reply. T h e d e n s i t y is a l l o w e d to vary
in the flame. T h e t e m p e r a t u r e a n d s p e c i e s m a s s fractions are o b t a i n e d from t h e s o l u t i o n of Eqs. (8) a n d (9), a n d t h e d e n s i t y is t h e n c o m p u t e d from the e q u a t i o n of state (5). T h i s is p o s s i b l e s i n c e the p r e s s u r e is a s s u m e d c o n s t a n t t h r o u g h o u t the flame in the model.