- Email: [email protected]
Transient laminar flame propagation in confined premixed gases: Numerical predictions
T R A N S I E N T LAMINAR F L A M E P R O P A G A T I O N IN C O N F I N E D P R E M I X E D GASES: N U M E R I C A L P R E D I C T I O N S DOUGLAS E. KOOKER US Army Ballistic Research Laboratory, Aberdeen Proving Ground, MD 21005
A model is developed for ignition and laminar flame propagation in premixed gases confined in a closed chamber. Predictions follow from a numerical solution to the one-dimensional compressible Navier-Stokes Equations. The solutior~ procedure properly accounts for multispecies diffusion, finite-rate chemical reactions, and interaction of the combustion wave with the confining geometry. A prediction of the flammability limit in ozone/oxygen is consistent with experimental data. Computations for a mixture of 18 mole % 03 in 0 2 indicate that the confined flame (1) remains far from chemical equilibrium, (2) pulsates with significant amplitude at a slowly increasing frequency predicted by Jost's acoustic analysis, and (3) accelerates to speeds which are an order of magnitude greater than a steady state flame at the pressure levels involved.
Introduction Combustion chamber flow field oscillations are known to occur in systems as diverse as ramjets, turbo-jet afterburners, industrial furnaces, rocket engines and guns. Although small amplitude organized disturbances are nearly impossible to eliminate, they are tolerable below a certain level without performance degradation. On occasion, however, these oscillations can be driven to extreme amplitudes which lead to catastrophic system failure. Despite the wide range of flow speeds and dependences on mixing and chemical kinetics, each of the above systems shares the common feature of partial confinement. The flow field interaction between the confining geometry and the combustion process exerts a critical influence on combustion behavior. Flame speed provides a dramatic illustration. Streng and Grosse i report that pure ozone burns at 475 c m / s e c as an unconfined flame sustained above a porous plate burner, but achieves detonation velocity when confined in a horizontal glass tube. Although this difference is known to be the result of confinement and influenced by a multitude of related factors such as wall roughness, tube diameter, etc., no analysis is yet capable of predicting the details of such a transition. Dynamic stability of the combustion process becomes the central issue. The numerous observations z-~ of cellular flame structure and oscillations in burner flames and partially confined premixed
gases represent a mild interaction with acoustic disturbances. When the limiting amplitudes are small, linearized analyses 7 can yield valuable insight. If the process in question cannot be represented as small departures from a stable flow, then a nonlinear analysis is required to predict limiting amplitudes. As mentioned in Ref. 1, any point in the region between the two extreme burning velocities is attainable during unstable or transient combustion. The present study is motivated by a desire to model solid propellant combustion in a gun-like environment, and to predict possible instabilities resulting from confinement. As a first step, the investigation is restricted to ignition and transient combustion of premixed gases confined in a closed one-dimensional chamber (See Fig. 1). One chamber boundary is assumed adiabatic and the other is subjected to a prescribed time-dependent heat flux which serves as the ignition stimulus. Ignition and transient laminar flame propagation are described with a rmmerical solution to the one-dimensional, compressible, Navier-Stokes Equations for a reacting flow field. Similar approaches have been reported in Refs. 8 and 9. A mixture of ozone and oxygen was selected as the reactive gas. Although the decomposition of ozone does not involve a branching chain as in many explosive systems, the simple three-species mechanism is well established and experimental data are available for the reaction rates. Furthermore, model
1329
1330
COMBUSTION OSCILLATIONS
03/02 z--0
I
(2)
pt+m~=0 Momentum
z--I
(3)
Fro. 1. Problem configuration.
z
z
Energy predictions can be compared to measurements of laminar flame speed and flammability limit.
E,+
[
(E+p)7
The reactive gas confined in the chamber is assumed to be a mixture of perfect gases, each of which may have a different molecular weight, specific heat capacity, heat of formation, etc. Mixture viscosity, thermal conductivity, and species diffusion coefficients are described at the molecular level, i.e., all transport is laminar. Mixture viscosity is found from Wilke's semi-empirical formula, used here as
~L= E
N
}
(P~),+
Pi
= t o , i = 1 , . . . , N - 1 (5)
+Vi z
State N
p = p~T, (1)
,=1 ~ Y k %
(4)
Species Continuity
y,p., N"
.
= (k~T. + Ixauu* + P E h~Y~V,
Analysis
N
m]
~ = Ro E
Y,/..s j
(6A)
t=1 N
h=EY,
h ,, h , = h ' ~ ' + c ; , ( T - T " ) Tn
where
+ IT~ % ' d T
(6B)
where m=. pu, E =- ph - p + pu2 /2, Y,=-p,/p.
with the analogous expression for thermal conductivity, kg. No constraints are imposed on the magnitude or variability of the mixture Lewis number or Prandtl number. The influence of radiation, gravity, bulk viscosity, the Dufour effect, the Sorer effect, and diffusion due to pressure gradients is neglected. Predictions of laminar flame propagation follow from the numerical solution of the one-dimensional, compressible, Navier-Stokes Equations for a reacting flow. In contrast to the common combustion assumption of uniform pressure, the present analysis retains the full momentum equation to account for the influence of the finite-size chamber. Furthermore, no artificial distinction is made between ignition and unsteady combustion; both are governed by identical equations, in strict conservation form, they are: Global Continuity
Yet unspecified is the species diffusion velocity, V,, which for a two-species system is often described by Fick's Law, i.e., W, =- Y, V, = -D,~Y, , i = 1,2
(7)
Attempts to generalize this concept to a multispecies mixture are based on similar coefficients (e.g., D~,~) which govern diffusion of the i th species into the "rest of the mixture." The fact that these coefficients, in general, are not equal creates the following dilemma. The trivial identities, X,~ Y, = 1 and p X,~ W, = 0, lead immediately to N
Y,~ = o
(8A)
D,~, Y,, = 0
(8B)
N
E
both of which cannot be satisfied. Some previous
1331
NUMERICAL PREDICTIONS analyses ~~ apparently violate Eq. (8B), while others ~'~3 adjust the coefficients D~2in a somewhat arbitrary manner to satisfy Eq. (8B). Obviously, internal consistency requires that the sum of the diffusion mass fluxes be zero. The diffusion velocities, V,, can be found from the binary diffusion coefficients in the manner given by Hirschfelder, ~4 - v,), i = 1 ..... N
--
(9)
where the contributions due to pressure gradients, body forces, and thermal gradients are neglected here. This same concept was employed by Edelman. L5Equation (9) is used in the present analysis in terms of species mass fractions, i.e., ~r
tr 0
I
-%
I $
.Of
~' X, Xk x,,
r
N
trn q
*
*
$
tmq+L ~a r
t*(ms) F1G. 2. Time-dependent heat flux imposed at the z = 0 boundary. tion enforces a statement of mechanical equilibrium, i.e.,
N
Y, W k - YkW, = Y', -"gr Y' E Y**/~s162 k k~ D,kcKCk k=l
i= 1..... N - 1
p~=[Ixa(~)
] atz=0,1 z
(11)
z
(10A) to determine boundary values of pressure ~6 consistent with the internal flow field. This concept is also presented in Ref. 17.
along with N
E Wk=0
(lOB)
Chemical Reaction Sequence
k~l
As noted by Hirschfelder, ~4 the N components of Eq. (9) are not linearly independent, and hence, Eq. (10B) is the required constraint. Solution of Eq. (10) at each flow field point yields N values of W, which are functions of the instantaneous composition of the mixture, the local mass fraction gradients, and the binary diffusion coefficients between all pairs of species.
Boundary Conditions It should be noted that the present problem involves two boundaries at finite locations; this rules out any transformation of the z coordinate with a density integral (e.g., Howarth transformation) for the purpose of eliminating direct dependence on the global continuity equation. In the inertial z coordinate, boundary values of the dependent variables must be specified at z = 0 and z = 1 (See Fig. 1). The assumption of impenetrable walls requires the mass average velocity and the species diffusion velocities to vanish at both boundaries. The adiabatic boundary at z = 1 requires T~ = 0 while the heated boundary at z = 0 is described by T, --- -qw/k~, where q *w(t) is shown in Fig. 2. Although no physical boundary condition is missing, the numerical solution procedure also requires boundary values of pressure (or density). Avoiding "reflection" conditions or other numerical artifices, the present solu-
Ozone decomposition is controlled by three reversible reactions: 0.~ + M kr~--5---02 + 0 + M
(12A)
kh I
k[ 2
03 + 0 ~ 202
(12B)
kb 2
02+M~iO+M (12C) kb3 where M represents any third body. In the present study, by assumption, the chamber initially contains a mixture of 03 and 02. Hence, only three species (03, 02, and O) ever participate in the reaction mechanism. Since the concentration level of O is low, only O~ and 03 are considered as M. Thus, Eq. (12) leads to five reversible reactions, whose rates were specified by the recommended values in the review by Johnston, ~s and essentially confirmed by the Leeds review. ~9 Temperature dependent specific heats for the three species were determined from the dual range curve fits given by Ref. 20. Transport properties were provided by Ref. 21 and are listed in Table I.
Numerical Method The conservative integration technique is constructed from the balance equations written for a
1332
C O M B U S T I O N OSCILLATIONS TABLE I Transport property data ~' Curve-fit expression
Ref. data
Binary Diffusion Coefficients (arm cma /sec) p D~, - 1 . 3 2 • 10-~T*"vv4 (T * ~ p D,3 - 1.6636 • 1 0 - ~ T *1~65 p* D~3 = 1.1779 • 10 -5 T *l"~
31 32 33
o
o
o
Q
--
Thermal Conductivity (cal / cm-sec- ~ k~o = 1.6037 • 10 -~ T ,o.7t k~o 2 = 6.5559 • 10 -7 T*~176 k~o ~ = 4.5904 X 10 -7 T ~ 0str~ Viscosity (g/ cm-sec) 2.6693 • 10 ~a N/T *~r o
34 35 33
12(~'21 = 1.6507 - 0.6688r + 0.10725r z r = I n ( T * / ~ ) , ~ = 106.7~ a = 3.05A Ix~z = 2.4929 • 10~" T ,0.7733, T * ~ 550 ~ K = 4.9219 • 10 -n T *~ T* ~ 550~ P~3 = 5.85 • 10 -~ T *~
36 32 32 33
control volume of length Az centered about each grid point. With the exception of transport properties and diffusion velocities, values of the flow field d e p e n d e n t variables are computed only at integer grid point locations; values required at half-integer locations are found from an arithmetic average of the adiacent neighbors, Adopting the "delta linearization" idea of Beam and Warming, ~ Eq. (2) for example becomes fPi + [~(fm~+, - fmj_~) = - a ( m ~ + , - m~_,)"
(13)
where fp = p,+t
_
9"
a =-=-A t / 2 A z , ~ = aO Equations (3) and (4) are written in an analogous manner after expanding
~p = 3P ~p + OP f m 3p dm N Op
with a similar expression for ST. The indicated partial derivatives are straightforward to compute from Eq. (6). The method is a two time level scheme, where the centering parameter 0 = 0.55 for all
computations shown. Accuracy is O(&zZ), but not quite O(At2). Full O(At ~) accuracy could be achieved by operating at the stability boundary, 0 = 1/2, or by using a three time level scheme. ~7 At this point, the mechanical flow field variables fp, fro, and fiE are coupled together in a ] • block (3 • 3) tridiagonal matrix, which yields the advanced time solution with a single inversion. This is easily accomplished with H i n d m a r s h ' s 22 routine. Thus the question of convergence never arises, as it would in an iterative or successive substitution method. However, still missing is a description of the chemical changes, i.e., the fY,'s. If the species continuity equations, containing "stiff" production rate terms, are simply expanded with delta linearization and then integrated along with the other equations at the Courant time step (or larger), significant errors can result. Elimination of this error may require a prohibitively small time step dictated by the fastest chemical reaction in the flow field. An alternative is operator-splitting, 2a where the chemical production rate terms are decoupled from the rest of the equation system and integrated separately as a system of ordinary differential equations. Dwyer ~4 has pointed out that, during the integration of the chemical system, each grid cell is momentarily a constant volume reactor without knowledge of local convection or diffusion. Thus, 2~ constant-volume heat release will serve to increase the pressure, even though it must decrease in the deflagration. The potential exists for numerical oscillation. A partial remedy is to form the driving terms for the system of ordinary differential equations by combining the chemical production rate
NUMERICAL PREDICTIONS terms with an explicit evaluation of convection and diffusion. These equations, coupled with the Energy equation, predict the chemical changes at the advanced time level "n + 1." The computational sequence is then: 1. With all coefficients evaluated at "n," solve Eq. (10) at each flow field point for W,"., i = 1,..., N 2. At each flow field point, integrate 15,=r
i = 1..... N (15A)
[mY,+pW,]",
k
+ IxRttu~ +
p
2
h ~W
~=l
(15B) z
from n to n + 1. This is done with Hindmarsh's E P I S O D E code, 25 where all Jacobian derivatives are evaluated analytically. 3. Compute y.+li
= P~n+l
~kn+l
and BY,:YT+'-Y7
i=
1..... N
4. With the bY~ known, invert the block tridiagohal matrix to obtain bp, Bin, BE. Although inversion of the block tridiagonal matrix is theoretically stable for an arbitrary time step, Eq. (15) contains explicit terms which are only conditionally stable. However, the controlling physical process in the present problem is wave motion occurring at the local sound speed. Hence, all predictions were generated by integration at the Courant time step; no instabilities were evident. All solutions were obtained without explicit numerical smoothing terms, which pose the threat of obliterating the combustion oscillations of interest here. However, a cell Reynolds number limit of two severely restricts the maximum size of the chamber.
I ~ w~-
As formulated here, the present model contains no adjustable constants which could be used to match theoretical predictions with experimental data. Assessment of predictive capability began with the flammability limit of an O:~/Oz mixture at atmospheric pressure. This limit is reported in the early experimental work of Streng and Grosse' to be 16/17 mole % O:~ in Oz, i.e., 16 would not sustain a flame. The theoretical prediction for 16 mole %, illustrated in Fig. 3, clearly shows a sustained flame.
'
'
I~n--T~f
;"
9001-
,ooF\\\\\\ 0
'
~
-
a 27.52
_~
I 32.32
"
a
700
~,~
'
x 9.70
+ 14.77 " v 18.12 "
.2
.4 DISTANCE,Z
34.73
~7 37.14 9 39.55
' '
.6 (N-D)
. . . . . . . .8
"
-
1.0
Fie. 3. Chamber temperature distribution timehistory for initial mixture containing 16 mole % 0 3 in O z. Flame is sustained. Initial conditions for all the examples are p ~ ~- 1 atm., u ~ ~- 0, T ~ = 298~ L ~ = 0.50 mm, J = 100, with q,~ given by Fig. 2. Sequential numbers listed on Figs. 3-6 are non-dimensional times from the onset of external heating. The apparent contradiction prompted further study which uncovered a later paper by Streng and Grosse 2a on the subject of quenching diameter and its influence on flammability limit. Without explanation, the revised data 2~ obtained in the 0.95 cm diameter horizontal tube of the earlier work' indicate a failure to propagate at 10 mole % O:~/0 2. All theoretical attempts to generate a sustained flame in this mixture were unsuccessful, including raising the boundary gas temperature 300 ~ above the constant volume explosion temperature of 966~ with an external heat addition more than an order of magnitude greater than that required to ignite a 25 mole % mixture. Typical results are shown in Fig. 4. Failure occurs simply because heat is conducted away from the reaction zone faster than it is generated by chemical
1100I~r-q~'r~1 i r~r~r-~r--Fn--~ ~+x114.629.i~-561l ~" 9 0 0 t ' - ~ ~,~ ~
Results
I
11001"~-~
300
}o
E = / - u ( E + p) + k~T~
~ --w-r~
1333
v 16.99
A 19.56 o 22.29 z~ i~m3 31.08 26.59 i 5.71 9 42.06
700 ~
oof\\ 3001 j 0
j'--~
..... .2
I ..... t , . L , , , .4 .6 .B 1.0 DISTANCE,Z (N-D)
F I G . 4. C h a m b e r temperature d i s t r i b u t i o n timehistory for initial mixture containing 10 mole % 0 3 in 0 z, Flame fails.
1334
COMBUSTION OSCILLATIONS
reaction. Although the present one-dimensional ,200 model simulates an infinite diameter tube, the ex1100 ~ 31.39 ~ . . ~ ~ 9 33.62 perimental flammability limit could be influenced looo \\\\\\\\\ ~, 35.85 by two dimensional effects such as lateral wall heat "\\\\\\\\\ ~ 38.08 loss. Neglecting these losses, however, makes the w 900 \\\\//I// o4o~o one-dimensional prediction conservative. The im8oo \\IIIIIII .4231 plication is that the experimental flammability limit 7ooj \\\///t/l ,,4472 ///////// ~46.~1 in the 0.95 cm diameter tube contains only a minimum influence from wall effects. Flame propagation experiments ~in confined mixtures containing more than 70 mole % O,s in O 2 have observed destruction of the tube, suggesting 3001 detonation. Although the present model is illequipped to predict a detonation wave, it should X 31.39 be capable of showing shock wave formation. The '* 3 3 . 6 2 pressure-time history resulting from the successful '9 3 5 . 8 5 ignition of a 70 mole % mixture is shown in Fig. A 38.08 5. The trend is clear, but the second-order finite-difO 40.30 ference solution has insufficient resolution to predict 9 42.51 /3 4 4 . 7 2 the actual shock wave. This combined with cell ~7 4 6 . 9 1 Reynolds numbers much greater than two eventually z 49.10 halts the computation. However, the above results demonstrate consistency with experimental observation. They are necessary but not sufficient proof of accuracy. ~ J J L i ~ h 0 I , ~ , I L The primary objective here is a study of transient -, , i i i flame propagation (without shock waves) in con.0025 x 31.39 fined Oa/O2 mixtures. Shown in Fig. 6 are simulta+ 33.62 neous time-histories of temperature, atomic oxygen V 35.85 mass fraction, and density distributions for a flame ~'E .0020 38.08 a 40.30 / propagating in 18 mole % O, in O 2. The continuous 9 42.51 temperature increase behind the flame front (see /3 4 4 . 7 2 Fig. 6A) is attributable to the behavior of atomic ~ .0015 ~7 4 6 . 9 1 . ~ oxygen. Similar to the steady-state ozone flame predictions of Hirschfelder 27 and Margolisy the present transient calculation shows a pronounced maximum in Yo just downstream of the leading temperature front. Since the flame speed is fast in .0005-.2 .4 .6 .8 1.0 comparison to the recombination rate of atomic DISTANCE,Z {N-D) oxygen, a reservoir of O is created behind the flame. Fic. 6. (A) Chamber temperature distribution, (B) atomic oxygen mass fraction distribution, and (C) , I , I , , I I I I I ] I I I density distribution time-histories during flame 1.40 x 12.28 propagation in 18 mole % O.s in O~. + 12.30
e6 o o 1 ~
o. o o 3 o ~
\\
I
J
O .OLO
'~ .~ 1.30 ::) 09
1.20
~
A 12,32 12.33 o T2.35
v
9 12.37 " 12.39 ~2 1 2 . 4 0
9 12.42 ~ 12.43
1.10 0
I
I
I
J
,2
L
,
,
|
.4
,
~ ,
I
~ ,
.6
t
I
.8
,
,
i
1.0
DISTANCE ,Z (N-D)
FIG. 5. Chamber pressure distribution time-history just after ignition of 70 mole % O~ in O 2mixture. Build-up to shock formation is evident.
Values of Yo near z = 0 (in Fig. 6B) are six orders of magnitude greater than those computed from an equilibrium thermodynamic calculation, i.e., the flame chemistry is far from equilibrium. However, the constant-volume explosion temperature of 1334~ will eventually be attained by subsequent recombination. The increasing density and temperature fields ahead of the flame (see Figs. 6A, C) are a direct consequence of confinement. Compression of the unburned mixture, trapped between the advancing flame and the fixed boundary at z = 1, is isentropic to within a fraction of one percent. Not unexpectedly, this compression of the reactants
1335
NUMERICAL PREDICTIONS will influence the rate of flame propagation. Flame speed predictions must be preceded by a criterion for flame location, which appears to have no unique definition. The present investigation chose an arbitrary value of temperature, 700~ in the high gradient region (see Fig. 6A) and then tracked the movement of this value in time and space. The resulting flame-front locus is plotted in Fig. 7A. Numerical time differentiation of this path produces the plot of flame velocity, as seen by an inertial observer, shown in Fig. 7B. The prediction is a pulsating laminar flame, whose velocity varies from approximately 300-600 cm/sec. Although the oscillation is a complicated function of several modes, the dominant mode is characterized by a slowly increasing frequency which can be predicted
.60
~ ~
:~ ~" ~S .50 5.40 ..
t
,
i
I
j
I
I
I
h
i
I
,
,
,
with an acoustic analysis originally developed by ]ostY Assuming stationary gases contained in a unit length closed chamber subdivided by a flame surface at z, where burned gas occupies the region 0 ---> z, Guenoche 3~ gives pbab tan
a,P
---pua,, tan La--'~] (16)
where P is the period of oscillation. For the flame movement shown in Fig. 6, the root of Eq. (16) for the period of the fundamental mode is within a few percent of the value predicted by the numerical flame tracking routine. Visualization of this mode is easier with a plot of the chamber velocity distribution (see Fig. 8A) during one-half cycle of the flame oscillation, i.e., the time period indicated by arrows in Fig. 7B. Superposed on Fig. 8A is the temperature distribution during the same time interval. A coupling between the velocities of the burned and unburned gas is involved here. When the flame velocity is minimum, the unburned gas is nearly stationary while the burned gas achieves its maximum velocity away from the flame (negative velocities are directed toward z = 0). When the flame velocity is maximum, the burned gas is at minimum velocity while the
,
y y!yV!
~410
390
~
~ 2.1 ~
~
t.t.i u-
380
350
'i
,
,
i 35
,
h
~
,
i 40
,
*
*
,
{ 45
. . . .
9
2.1,
5O
TIME I N - D )
F{o. 7. (A) Flame position (based on 700~ vs. time, (B) flame velocity (as seen by inertial observer) vs. time, and (C) flame velocity minus local gas velocity at flame location vs. time. Corresponds to results shown in Fig. 6.
,2
.4
DISTANCE,Z
.6
(N-O)
,8
1.0
Fie. 8. (A)Chamber velocity distribution, and (B) chamber pressure distribution during one-half cycle oscillation (indicated by arrows in Fig. 7B).
1336
COMBUSTION OSCILLATIONS
unburned gas achieves its maximum velocity away from the flame. The chamber pressure distributions for the identical time interval are shown in Fig. 8B. The behavior of the pressure profile on the burned gas side is a result of a sound speed which is more than twice the value in the unburned gas. The origin of the flame oscillation, in the present example, can be traced to the heating phase before ignition of the mixture. Rapid heat transfer at one boundary excites acoustic wave motion in the chamber as pressure increases are equilibrated, i.e., thermally-driven convection. Thus, small amplitude standing waves exist in the confined gas mixture before ignition. After ignition, these waves are amplified and driven by the flame as it propagates away from the wall boundary. In other problems where the time to ignition is on the order of the acoustic travel time in the chamber, the onset of oscillation is delayed until the flame propagates some distance into the chamber. Although not required in the present solution, there is considerable interest in the question of burning velocity, i,e., the flame velocity relative to the unburned gas velocity. A trivial concept in steady flame propagation, it becomes ambiguous in a transient problem where the unburned gases are in non-uniform motion. Without resolving this issue, it was decided to reference the flame velocity to the local gas velocity at the flame location point. The numerical solution predicts this local gas velocity is nearly inphase with the flame velocity, again confirming the acoustic origin of the oscillation. The difference between the two velocities is plotted in Fig. 7C for the flame path in Fig, 7A. Note that this velocity difference is positive, i.e., the flame always moves into the unburned gas. Since the flame-tracking routine has an estimated error bound of +__5cm/sec, it is not possible to make definitive statements about the fine structure. Qualitatively, however, the curve indicates an overall flame acceleration which is not evident in the flame velocity seen by an inertial observer (Fig. 7B). As the flame propagates through the chamber, the average pressure level continuously increases. During the time interval of Figs. 6 and 7, the mean chamber pressure varies from 1.5 to 2.2 atm. The pressurization rate when the flame is near midchamber is approximately 3 • 105 psi/sec, and increasing. The steady state laminar burning velocity in 18 mole % O~ in O~ at 1 atm is 12-24 cm/sec ~, depending upon the flame shape assumed, If the velocity difference discussed above is assigned the role of a burning velocity, then clearly the transient flame is propagating into the unburned mixture more than an order of magnitude faster than a steady state flame at the pressure levels involved. The explanation is based on the dynamics of the process, A direct consequence of confinement is compression of the unburned mixture ahead of the flame. Thus,
chemical conversion in the flame releases essentially a fixed amount of energy into a continuously increasing temperature and density field. This increases the rate of heat release, which in turn steepens the temperature profile in the flame (this is evident from a close inspection of Fig. 6A). The cumulative effect leads to an acceleration of the flame.
Concluding R e m a r k s A one-dimensional model, based on a numerical solution of the compressible Navier-Stokes Equations, was developed to study transient laminar flame propagation in confined premixed gases. The solution procedure properly accounts for multispecies diffusion, finite-rate chemical reactions, and interaction of the combustion wave with the confining geometry. In this initial study, using independentlydetermined values for the chemical reaction rates and other gas properties, the model demonstrated encouraging predictive capability, Obviously, much remains to be done. Although the present investigation examined only ozone/oxygen mixtures, other exothermic premixed systems should exhibit similar oscillation and acceleration behavior. Theoretically, because of the nonlinear coupling between the equations governing conservation of momentum and energy, confined combustion is inherently oscillatory. In practice, the limiting amplitudes may be negligible as the result of system losses, etc, However, any analysis of confined or partially confined combustion based on assumptions which preclude oscillatory behavior will lead to predictions with unknown error, Nomenclature
a*
reference speed of sound (cm/sec). , " i c heat at constant pressure. cvjT 9, r, z ; speclf I*,, binary diffusion coefficient beD~k D*~ /a*r ~ tween i th and k'" species. E ph - p + pu 2 /2; total energy per unit volume. h h*/a*2; enthalpy, ~=~ Y,h,. h~' heat of formation of the i 'h species at T ~ J total number of computational cells. k~ k~ T*/O~* a*~L*; thermal conductivity. L* reference length (cm). m pu, mass flux. ,,~', molecular weight of the i 'h species (g/g-mole). 1/(~,Y,/~r mixture molecular weight (g/g-mole). N total number of species. ,1 ,a, 2 P P /Or 9 ; pressure. P a * / f ~ L *; nondimensional period (f * = cyeles/seq). q,~ qw/p~a r ; external heat flux, q * = c a l /
cp,
em
-sec.
NUMERICAL PREDICTIONS Re
p* a* L * / I x * ; reference value of Reynolds number. Ro R*o T r* f a r ,2 ; universal gas constant. Ro E ~_, Y,/.~r gas constant. t t * a * / L * ; t i m e ( h e r e , L * / a * = 1.638• 10 -6 sec). tin,, cut-off time of external heat flux (See Fig. 2). T T * / T *; temperature. T* reference value of temperature (~ u u * / a *~; velocity. V~ V ~ / a * ; diffusion velocity of i 'h species. W, Y,V,. X, mole fraction of i th species. Y~ mass fraction of i th species. z z * / L*; distance coordinate. a At/2Az, OR. Ix Ix */Ix *; coefficient of viscosity. Ix* reference value of coefficient of viscosity (g/cm-sec). I~, (4/3) ~ / R e . p p * / p * ; density. p* reference value of density (g/cm:~), 0 weight factor, O = 0.55 here. to 0~"L /p~ a~ ; chemical production rate for species i. Superscripts * n
dimensional quantity. time level,
Subscripts i i ( )~ (), (), ( )~
pertaining to the i th species, grid point location. partial derivative with respect to z. partial derivative with respect to t. u n b u r n e d gas. b u r n e d gas. Acknowledgment
The author gratefully acknowledges the substantial contribution of Mr. R. D, Anderson who executed the computer program and constructed the CALC O M P plots.
REFERENCES I. STRENG, A. G. ANt~GROSSE, A. V.: Sixth Symposium (International) on Combustion, p, 264, Reinhold, 1957. 2. MARXSTEIN, G. H.: Fourth S y m p o s i u m (International) on Combustion, p. 44, Williams and Wilkins, 1953. 3. MAR~;STEIN, G. H. AND SOMERS, L. H,: ibid, p. 527.
1337
4. JosT, W., KauG, J., AND S1EG, L.: ibid, p. 535. 5. KASrAN, W. E.: ibid, p. 575. 6. MARKSTEIN,G. H. (Eo.): Non-Steady Flame Propagation, Chapters D, G, MacMillan, 1964. 7. JONES, H.: Proc. R. Soe. Load. A, 353, p. 459 (1977). 8, BUTLER, T. D. AND O'RoURKE, P. J.: Sixteenth S y m p o s i u m (International) on Combustion, p. 1503, The Combustion Institute, 1977. 9. WESTBROOK,C. K. AND CHANG, J. S.: "A Theoretical Study of Flame Propagation Through Stratified Media," Lawrence Livermore Laboratory Report UCRL-78925, 1977. 10. KOTrlARI, A. P., ANDERSON, J. D., JR., AND JONES, E.: AIAA J,, 15, p. 92 (1977). 11. CRAMAROSSA,F. ANDDIXON-LEwIS, G.: Comb. and Flame, 16, p. 243 (1971). 12. RWARO, W. C., FARMER, O. A., AND BUTLER, T. D.: "'RICE: A Computer Program for Multicomponent Chemically Reactive Flows at All Speeds," Los Alamos Scientific Laboratory Report LA-5812, 1975. 13. KEE, R. J. AND MILLER, J. A.: "A Split-Operator, Finite Difference Solution for Axisymmetrie Laminar Jet Diffusion Flames," Sandia Livermore Laboratories, SAND77-8502, 1977. 14. HIRSCHFELDER, J. O., CURTISS, C. F., AND BIRD, R. B.: Molecular Theory o f Gases and Liquids, Chapter 7, John Wiley, 1954. 15. EDELMAN, R. B., FORTUNE, O. F., WEILERSTEIN, G., COCHRAN, T. H., AND HAGGARD, J. B., JR.: Fourteenth S y m p o s i u m (International) on Combustion, p. 399, The Combustion Institute, 1973. 16. KoogE~, D. E.: "Numerical Predictions for Ignition and Flamespreading in a Confined O z o n e / O x y g e n Mixture," 14th J A N N A F Combustion Meeting, CPIA Publication 292, p. 335, Vol. II, 1977, 17. BEA~I, R. M. aND WARMING, R. F.: AIAA J., 16, p. 393 (1978), See also Warming, R. F. and Beam, R. M.: SIAM-AMS Proceedings, 11, 1977. 18. JOHNSTON, H. S.: "'Gas Phase Reaction Kinetics of Neutral Oxygen Species," National Bureau of Standards, NSRDS-NBS-20, 1968. 19. BAULCH, D. L., DRYSDALE, n . D., HORNE, D. G., AND LLOYD, A. C.: Evaluated Kinetic Data for High Temperature Reactions, Butterworths, 1973. 20. GORDON, S. AND McBRIDE, B. J.: "Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks, and Chapman-Jouget Detonations," NASA/ Lewis, NASA-SP-273, 1971 (1976 program version). 21. HEIMERL, J. M. AND COVVEE, T. P.: Private Communication 1978 (to be published). 22. H1NOMARSH, A. C.: "'Solution of Block-Tridiagonal Systems of Linear Algebraic Equa-
1338
23. 24.
25.
26. 27. 28. 29. 30.
COMBUSTION OSCILLATIONS
tions," Lawrence Livermore Laboratory, UCID30150, 1977. YANENKO,N. N.: The MethodofFractional Steps, Springer-Verlag, 1971. DwYER, H.: Private communication Dec 1977, See also Dwyer, H. and Otey, G.: "A Numerical Study of the Interaction of Fast Chemistry and Diffusion," AIAA Preprint 78-946, 1978. HINDMARSH,A. C. ANOBYRNE,G. D.: "'EPISODE: An Experimental Package for the Integration of Systems of Ordinary Differential Equations," Lawrence Livermore Laboratory, UCRL-75868, 1975. STRENG, A. G. AND GROSS, A. V.: Comb. and Flame, 5, p. 81 (1961) (See Fig. 2). HIRSCHFELDER,]. O., CURTISS, C. F., AND CAMeBELL, D. E.: J. Phys. Chem., 57, p. 403 (1953). MARGOLIS, S. B.: J. Comp. Phys., 27, p. 410 (1978). ]OST, W.: Explosion and Combustion Processes in Gases, p. 103, McGraw-Hill, 1946. GUENOCHE, H.: "Flame Propagation in Tubes
31. 32. 33.
34. 35.
36.
and Closed Vessels," p. 149 in Nonsteady Flame Propagation, G. H. Markstein (ed.), MacMillan, 1964. MARaEaO,T. R. ANDMASON,E. A.: J. Phys., Chem. Ref. Data, 1, No. I (1972). HANLEr, H. ]. M. AND ELY, ]. F.: J. Phys. Chem. Ref. Data, 2, No. 4 (1973). TRETTEa, V. J., JR.: "'An Experimental Determination of the Viscosity of Gaseous Ozone," M.S. Thesis, Chem. Engineering, University of Idaho, 1966. DALGARNO,A. AND SMITH, F. J.: Planetary Space Sci., 9, p. 1 (1962). TOULOUKIAN,Y. S., L1LEY, P. E., AND SAXENA,S. C.: Thermophysical Properties of Matter, Vol. 3; IFF-Plenum, 1970. EDELMAN,R. B., FORTUNE, O., AND WEILERSTEIN, G.: "Analytical Study of Gravity Effects on Laminar Diffusion Flames," General Applied Science Laboratory, GASL TR-771, 1972 (also NASA CR-120921).
COMMENTS B. T. Zinn, Georgia Institute of Technology, USA. Your results indicate the formation of a shock wave, and I wonder if you have identified one and have you followed its evolution with time? Also, do your results indicate the presence of the fundamental and higher harmonics of the acoustic modes of the chamber? Also, in previous combustion oscillation studies the onset of oscillations has been related to the presence of some physical/kinetic mechanism whose characteristic time is of the same order of magnitude as the period of oscillation of one of the combustor's modes, Can you identify such a process for or a physical-acoustic interaction from your results? Finally, will you comment upon the length of your computations and the major contributors to this time.
Authors'Reply. O f t h e c o m p u t a t i o n s m a d e t o date, only mixtures containing 70 mole percent 0 3 / 0 2 and higher indicate the presence of a shock wave, and this occurs at ignition. The gas velocity ahead of the forming shock wave quickly exceeds the maximum value dictated by a cell Reynolds number equal to two. Furthermore, without a prohibitively fine grid mesh system, the numerical solution is incapable of resolving the interior structure of the shock wave. Momentarily ignoring these difficulties leads to a steep-fronted wave containing oscillations
on the high-pressure side. Although tolerable in a non-reacting inviscid "shock-capturing" solution, these numerical oscillations quickly destroy the present calculation w h e n coupled into the finite-rate chemistry package. Several frequencies are evident in the plot of flame velocity versus time (see Fig. 7B). As shown in the paper, the dominant mode is the fundamental acoustic mode of the chamber, assuming the chamber is subdivided by a flame surface separating burned and unburned gases. The first harmonic of this mode is also present, until the flame passes midehamber. These modes are not stationary in the usual acoustic sense, since their frequencies change with the advancing flame. Most likely, the combustion oscillation studies referred to in your question are the instability analyses of solid and liquid rocket engines. Note that both of these analyses include a basic combustion mechanism which is independent of the flow field in which oscillations occur. The solid propellant regression rate is customarily assumed to be an exponential function of surface temperature which is controlled by the time-dependent thermal wave in the u n b u r n e d solid; hence, there is a characteristic relaxation time which is only a function of the propellant. In liquid combustion instability analyses, a host of unknowns involving droplet vaporization, mixing, and chemical kinetics are often lumped into two parameters called " n " and '"r";
NUMERICAL PREDICTIONS both are specified independently of the rocket engine flow field. Then the question is posed whether these driving sources can be "in-phase" or "out-of-phase" with the natural acoustic modes of the confining geometry. The present problem is different in the following sense. The driving source for the premixed flame is the heat release distribution which follows directly from the chemical production rate terms; these are sensitive functions of density and temperature in the local flow field. What can and does happen is that the flow field oscillations alter the heat release
1339
distribution. This coupling mechanism allows considerable latitude in which modes can be driven. O f the cases run to date, it appears that Rayleighs' criterion is satisfied by an easily accessible mode and hence the flame oscillation is sustained even though the frequency may change. The present computations require over 3 minu t e s / 1000 steps on a CDC 7600. Depending on flame speed, between 8000-12,000 steps are necessary to traverse the chamber. More than two-thirds of this time is consumed by the stiff integrator chemistry package.
A model is developed for ignition and laminar flame propagation in premixed gases confined in a closed chamber. Predictions follow from a numerical solution to the one-dimensional compressible Navier-Stokes Equations. The solutior~ procedure properly accounts for multispecies diffusion, finite-rate chemical reactions, and interaction of the combustion wave with the confining geometry. A prediction of the flammability limit in ozone/oxygen is consistent with experimental data. Computations for a mixture of 18 mole % 03 in 0 2 indicate that the confined flame (1) remains far from chemical equilibrium, (2) pulsates with significant amplitude at a slowly increasing frequency predicted by Jost's acoustic analysis, and (3) accelerates to speeds which are an order of magnitude greater than a steady state flame at the pressure levels involved.
Introduction Combustion chamber flow field oscillations are known to occur in systems as diverse as ramjets, turbo-jet afterburners, industrial furnaces, rocket engines and guns. Although small amplitude organized disturbances are nearly impossible to eliminate, they are tolerable below a certain level without performance degradation. On occasion, however, these oscillations can be driven to extreme amplitudes which lead to catastrophic system failure. Despite the wide range of flow speeds and dependences on mixing and chemical kinetics, each of the above systems shares the common feature of partial confinement. The flow field interaction between the confining geometry and the combustion process exerts a critical influence on combustion behavior. Flame speed provides a dramatic illustration. Streng and Grosse i report that pure ozone burns at 475 c m / s e c as an unconfined flame sustained above a porous plate burner, but achieves detonation velocity when confined in a horizontal glass tube. Although this difference is known to be the result of confinement and influenced by a multitude of related factors such as wall roughness, tube diameter, etc., no analysis is yet capable of predicting the details of such a transition. Dynamic stability of the combustion process becomes the central issue. The numerous observations z-~ of cellular flame structure and oscillations in burner flames and partially confined premixed
gases represent a mild interaction with acoustic disturbances. When the limiting amplitudes are small, linearized analyses 7 can yield valuable insight. If the process in question cannot be represented as small departures from a stable flow, then a nonlinear analysis is required to predict limiting amplitudes. As mentioned in Ref. 1, any point in the region between the two extreme burning velocities is attainable during unstable or transient combustion. The present study is motivated by a desire to model solid propellant combustion in a gun-like environment, and to predict possible instabilities resulting from confinement. As a first step, the investigation is restricted to ignition and transient combustion of premixed gases confined in a closed one-dimensional chamber (See Fig. 1). One chamber boundary is assumed adiabatic and the other is subjected to a prescribed time-dependent heat flux which serves as the ignition stimulus. Ignition and transient laminar flame propagation are described with a rmmerical solution to the one-dimensional, compressible, Navier-Stokes Equations for a reacting flow field. Similar approaches have been reported in Refs. 8 and 9. A mixture of ozone and oxygen was selected as the reactive gas. Although the decomposition of ozone does not involve a branching chain as in many explosive systems, the simple three-species mechanism is well established and experimental data are available for the reaction rates. Furthermore, model
1329
1330
COMBUSTION OSCILLATIONS
03/02 z--0
I
(2)
pt+m~=0 Momentum
z--I
(3)
Fro. 1. Problem configuration.
z
z
Energy predictions can be compared to measurements of laminar flame speed and flammability limit.
E,+
[
(E+p)7
The reactive gas confined in the chamber is assumed to be a mixture of perfect gases, each of which may have a different molecular weight, specific heat capacity, heat of formation, etc. Mixture viscosity, thermal conductivity, and species diffusion coefficients are described at the molecular level, i.e., all transport is laminar. Mixture viscosity is found from Wilke's semi-empirical formula, used here as
~L= E
N
}
(P~),+
Pi
= t o , i = 1 , . . . , N - 1 (5)
+Vi z
State N
p = p~T, (1)
,=1 ~ Y k %
(4)
Species Continuity
y,p., N"
.
= (k~T. + Ixauu* + P E h~Y~V,
Analysis
N
m]
~ = Ro E
Y,/..s j
(6A)
t=1 N
h=EY,
h ,, h , = h ' ~ ' + c ; , ( T - T " ) Tn
where
+ IT~ % ' d T
(6B)
where m=. pu, E =- ph - p + pu2 /2, Y,=-p,/p.
with the analogous expression for thermal conductivity, kg. No constraints are imposed on the magnitude or variability of the mixture Lewis number or Prandtl number. The influence of radiation, gravity, bulk viscosity, the Dufour effect, the Sorer effect, and diffusion due to pressure gradients is neglected. Predictions of laminar flame propagation follow from the numerical solution of the one-dimensional, compressible, Navier-Stokes Equations for a reacting flow. In contrast to the common combustion assumption of uniform pressure, the present analysis retains the full momentum equation to account for the influence of the finite-size chamber. Furthermore, no artificial distinction is made between ignition and unsteady combustion; both are governed by identical equations, in strict conservation form, they are: Global Continuity
Yet unspecified is the species diffusion velocity, V,, which for a two-species system is often described by Fick's Law, i.e., W, =- Y, V, = -D,~Y, , i = 1,2
(7)
Attempts to generalize this concept to a multispecies mixture are based on similar coefficients (e.g., D~,~) which govern diffusion of the i th species into the "rest of the mixture." The fact that these coefficients, in general, are not equal creates the following dilemma. The trivial identities, X,~ Y, = 1 and p X,~ W, = 0, lead immediately to N
Y,~ = o
(8A)
D,~, Y,, = 0
(8B)
N
E
both of which cannot be satisfied. Some previous
1331
NUMERICAL PREDICTIONS analyses ~~ apparently violate Eq. (8B), while others ~'~3 adjust the coefficients D~2in a somewhat arbitrary manner to satisfy Eq. (8B). Obviously, internal consistency requires that the sum of the diffusion mass fluxes be zero. The diffusion velocities, V,, can be found from the binary diffusion coefficients in the manner given by Hirschfelder, ~4 - v,), i = 1 ..... N
--
(9)
where the contributions due to pressure gradients, body forces, and thermal gradients are neglected here. This same concept was employed by Edelman. L5Equation (9) is used in the present analysis in terms of species mass fractions, i.e., ~r
tr 0
I
-%
I $
.Of
~' X, Xk x,,
r
N
trn q
*
*
$
tmq+L ~a r
t*(ms) F1G. 2. Time-dependent heat flux imposed at the z = 0 boundary. tion enforces a statement of mechanical equilibrium, i.e.,
N
Y, W k - YkW, = Y', -"gr Y' E Y**/~s162 k k~ D,kcKCk k=l
i= 1..... N - 1
p~=[Ixa(~)
] atz=0,1 z
(11)
z
(10A) to determine boundary values of pressure ~6 consistent with the internal flow field. This concept is also presented in Ref. 17.
along with N
E Wk=0
(lOB)
Chemical Reaction Sequence
k~l
As noted by Hirschfelder, ~4 the N components of Eq. (9) are not linearly independent, and hence, Eq. (10B) is the required constraint. Solution of Eq. (10) at each flow field point yields N values of W, which are functions of the instantaneous composition of the mixture, the local mass fraction gradients, and the binary diffusion coefficients between all pairs of species.
Boundary Conditions It should be noted that the present problem involves two boundaries at finite locations; this rules out any transformation of the z coordinate with a density integral (e.g., Howarth transformation) for the purpose of eliminating direct dependence on the global continuity equation. In the inertial z coordinate, boundary values of the dependent variables must be specified at z = 0 and z = 1 (See Fig. 1). The assumption of impenetrable walls requires the mass average velocity and the species diffusion velocities to vanish at both boundaries. The adiabatic boundary at z = 1 requires T~ = 0 while the heated boundary at z = 0 is described by T, --- -qw/k~, where q *w(t) is shown in Fig. 2. Although no physical boundary condition is missing, the numerical solution procedure also requires boundary values of pressure (or density). Avoiding "reflection" conditions or other numerical artifices, the present solu-
Ozone decomposition is controlled by three reversible reactions: 0.~ + M kr~--5---02 + 0 + M
(12A)
kh I
k[ 2
03 + 0 ~ 202
(12B)
kb 2
02+M~iO+M (12C) kb3 where M represents any third body. In the present study, by assumption, the chamber initially contains a mixture of 03 and 02. Hence, only three species (03, 02, and O) ever participate in the reaction mechanism. Since the concentration level of O is low, only O~ and 03 are considered as M. Thus, Eq. (12) leads to five reversible reactions, whose rates were specified by the recommended values in the review by Johnston, ~s and essentially confirmed by the Leeds review. ~9 Temperature dependent specific heats for the three species were determined from the dual range curve fits given by Ref. 20. Transport properties were provided by Ref. 21 and are listed in Table I.
Numerical Method The conservative integration technique is constructed from the balance equations written for a
1332
C O M B U S T I O N OSCILLATIONS TABLE I Transport property data ~' Curve-fit expression
Ref. data
Binary Diffusion Coefficients (arm cma /sec) p D~, - 1 . 3 2 • 10-~T*"vv4 (T * ~ p D,3 - 1.6636 • 1 0 - ~ T *1~65 p* D~3 = 1.1779 • 10 -5 T *l"~
31 32 33
o
o
o
Q
--
Thermal Conductivity (cal / cm-sec- ~ k~o = 1.6037 • 10 -~ T ,o.7t k~o 2 = 6.5559 • 10 -7 T*~176 k~o ~ = 4.5904 X 10 -7 T ~ 0str~ Viscosity (g/ cm-sec) 2.6693 • 10 ~a N/T *~r o
34 35 33
12(~'21 = 1.6507 - 0.6688r + 0.10725r z r = I n ( T * / ~ ) , ~ = 106.7~ a = 3.05A Ix~z = 2.4929 • 10~" T ,0.7733, T * ~ 550 ~ K = 4.9219 • 10 -n T *~ T* ~ 550~ P~3 = 5.85 • 10 -~ T *~
36 32 32 33
control volume of length Az centered about each grid point. With the exception of transport properties and diffusion velocities, values of the flow field d e p e n d e n t variables are computed only at integer grid point locations; values required at half-integer locations are found from an arithmetic average of the adiacent neighbors, Adopting the "delta linearization" idea of Beam and Warming, ~ Eq. (2) for example becomes fPi + [~(fm~+, - fmj_~) = - a ( m ~ + , - m~_,)"
(13)
where fp = p,+t
_
9"
a =-=-A t / 2 A z , ~ = aO Equations (3) and (4) are written in an analogous manner after expanding
~p = 3P ~p + OP f m 3p dm N Op
with a similar expression for ST. The indicated partial derivatives are straightforward to compute from Eq. (6). The method is a two time level scheme, where the centering parameter 0 = 0.55 for all
computations shown. Accuracy is O(&zZ), but not quite O(At2). Full O(At ~) accuracy could be achieved by operating at the stability boundary, 0 = 1/2, or by using a three time level scheme. ~7 At this point, the mechanical flow field variables fp, fro, and fiE are coupled together in a ] • block (3 • 3) tridiagonal matrix, which yields the advanced time solution with a single inversion. This is easily accomplished with H i n d m a r s h ' s 22 routine. Thus the question of convergence never arises, as it would in an iterative or successive substitution method. However, still missing is a description of the chemical changes, i.e., the fY,'s. If the species continuity equations, containing "stiff" production rate terms, are simply expanded with delta linearization and then integrated along with the other equations at the Courant time step (or larger), significant errors can result. Elimination of this error may require a prohibitively small time step dictated by the fastest chemical reaction in the flow field. An alternative is operator-splitting, 2a where the chemical production rate terms are decoupled from the rest of the equation system and integrated separately as a system of ordinary differential equations. Dwyer ~4 has pointed out that, during the integration of the chemical system, each grid cell is momentarily a constant volume reactor without knowledge of local convection or diffusion. Thus, 2~ constant-volume heat release will serve to increase the pressure, even though it must decrease in the deflagration. The potential exists for numerical oscillation. A partial remedy is to form the driving terms for the system of ordinary differential equations by combining the chemical production rate
NUMERICAL PREDICTIONS terms with an explicit evaluation of convection and diffusion. These equations, coupled with the Energy equation, predict the chemical changes at the advanced time level "n + 1." The computational sequence is then: 1. With all coefficients evaluated at "n," solve Eq. (10) at each flow field point for W,"., i = 1,..., N 2. At each flow field point, integrate 15,=r
i = 1..... N (15A)
[mY,+pW,]",
k
+ IxRttu~ +
p
2
h ~W
~=l
(15B) z
from n to n + 1. This is done with Hindmarsh's E P I S O D E code, 25 where all Jacobian derivatives are evaluated analytically. 3. Compute y.+li
= P~n+l
~kn+l
and BY,:YT+'-Y7
i=
1..... N
4. With the bY~ known, invert the block tridiagohal matrix to obtain bp, Bin, BE. Although inversion of the block tridiagonal matrix is theoretically stable for an arbitrary time step, Eq. (15) contains explicit terms which are only conditionally stable. However, the controlling physical process in the present problem is wave motion occurring at the local sound speed. Hence, all predictions were generated by integration at the Courant time step; no instabilities were evident. All solutions were obtained without explicit numerical smoothing terms, which pose the threat of obliterating the combustion oscillations of interest here. However, a cell Reynolds number limit of two severely restricts the maximum size of the chamber.
I ~ w~-
As formulated here, the present model contains no adjustable constants which could be used to match theoretical predictions with experimental data. Assessment of predictive capability began with the flammability limit of an O:~/Oz mixture at atmospheric pressure. This limit is reported in the early experimental work of Streng and Grosse' to be 16/17 mole % O:~ in Oz, i.e., 16 would not sustain a flame. The theoretical prediction for 16 mole %, illustrated in Fig. 3, clearly shows a sustained flame.
'
'
I~n--T~f
;"
9001-
,ooF\\\\\\ 0
'
~
-
a 27.52
_~
I 32.32
"
a
700
~,~
'
x 9.70
+ 14.77 " v 18.12 "
.2
.4 DISTANCE,Z
34.73
~7 37.14 9 39.55
' '
.6 (N-D)
. . . . . . . .8
"
-
1.0
Fie. 3. Chamber temperature distribution timehistory for initial mixture containing 16 mole % 0 3 in O z. Flame is sustained. Initial conditions for all the examples are p ~ ~- 1 atm., u ~ ~- 0, T ~ = 298~ L ~ = 0.50 mm, J = 100, with q,~ given by Fig. 2. Sequential numbers listed on Figs. 3-6 are non-dimensional times from the onset of external heating. The apparent contradiction prompted further study which uncovered a later paper by Streng and Grosse 2a on the subject of quenching diameter and its influence on flammability limit. Without explanation, the revised data 2~ obtained in the 0.95 cm diameter horizontal tube of the earlier work' indicate a failure to propagate at 10 mole % O:~/0 2. All theoretical attempts to generate a sustained flame in this mixture were unsuccessful, including raising the boundary gas temperature 300 ~ above the constant volume explosion temperature of 966~ with an external heat addition more than an order of magnitude greater than that required to ignite a 25 mole % mixture. Typical results are shown in Fig. 4. Failure occurs simply because heat is conducted away from the reaction zone faster than it is generated by chemical
1100I~r-q~'r~1 i r~r~r-~r--Fn--~ ~+x114.629.i~-561l ~" 9 0 0 t ' - ~ ~,~ ~
Results
I
11001"~-~
300
}o
E = / - u ( E + p) + k~T~
~ --w-r~
1333
v 16.99
A 19.56 o 22.29 z~ i~m3 31.08 26.59 i 5.71 9 42.06
700 ~
oof\\ 3001 j 0
j'--~
..... .2
I ..... t , . L , , , .4 .6 .B 1.0 DISTANCE,Z (N-D)
F I G . 4. C h a m b e r temperature d i s t r i b u t i o n timehistory for initial mixture containing 10 mole % 0 3 in 0 z, Flame fails.
1334
COMBUSTION OSCILLATIONS
reaction. Although the present one-dimensional ,200 model simulates an infinite diameter tube, the ex1100 ~ 31.39 ~ . . ~ ~ 9 33.62 perimental flammability limit could be influenced looo \\\\\\\\\ ~, 35.85 by two dimensional effects such as lateral wall heat "\\\\\\\\\ ~ 38.08 loss. Neglecting these losses, however, makes the w 900 \\\\//I// o4o~o one-dimensional prediction conservative. The im8oo \\IIIIIII .4231 plication is that the experimental flammability limit 7ooj \\\///t/l ,,4472 ///////// ~46.~1 in the 0.95 cm diameter tube contains only a minimum influence from wall effects. Flame propagation experiments ~in confined mixtures containing more than 70 mole % O,s in O 2 have observed destruction of the tube, suggesting 3001 detonation. Although the present model is illequipped to predict a detonation wave, it should X 31.39 be capable of showing shock wave formation. The '* 3 3 . 6 2 pressure-time history resulting from the successful '9 3 5 . 8 5 ignition of a 70 mole % mixture is shown in Fig. A 38.08 5. The trend is clear, but the second-order finite-difO 40.30 ference solution has insufficient resolution to predict 9 42.51 /3 4 4 . 7 2 the actual shock wave. This combined with cell ~7 4 6 . 9 1 Reynolds numbers much greater than two eventually z 49.10 halts the computation. However, the above results demonstrate consistency with experimental observation. They are necessary but not sufficient proof of accuracy. ~ J J L i ~ h 0 I , ~ , I L The primary objective here is a study of transient -, , i i i flame propagation (without shock waves) in con.0025 x 31.39 fined Oa/O2 mixtures. Shown in Fig. 6 are simulta+ 33.62 neous time-histories of temperature, atomic oxygen V 35.85 mass fraction, and density distributions for a flame ~'E .0020 38.08 a 40.30 / propagating in 18 mole % O, in O 2. The continuous 9 42.51 temperature increase behind the flame front (see /3 4 4 . 7 2 Fig. 6A) is attributable to the behavior of atomic ~ .0015 ~7 4 6 . 9 1 . ~ oxygen. Similar to the steady-state ozone flame predictions of Hirschfelder 27 and Margolisy the present transient calculation shows a pronounced maximum in Yo just downstream of the leading temperature front. Since the flame speed is fast in .0005-.2 .4 .6 .8 1.0 comparison to the recombination rate of atomic DISTANCE,Z {N-D) oxygen, a reservoir of O is created behind the flame. Fic. 6. (A) Chamber temperature distribution, (B) atomic oxygen mass fraction distribution, and (C) , I , I , , I I I I I ] I I I density distribution time-histories during flame 1.40 x 12.28 propagation in 18 mole % O.s in O~. + 12.30
e6 o o 1 ~
o. o o 3 o ~
\\
I
J
O .OLO
'~ .~ 1.30 ::) 09
1.20
~
A 12,32 12.33 o T2.35
v
9 12.37 " 12.39 ~2 1 2 . 4 0
9 12.42 ~ 12.43
1.10 0
I
I
I
J
,2
L
,
,
|
.4
,
~ ,
I
~ ,
.6
t
I
.8
,
,
i
1.0
DISTANCE ,Z (N-D)
FIG. 5. Chamber pressure distribution time-history just after ignition of 70 mole % O~ in O 2mixture. Build-up to shock formation is evident.
Values of Yo near z = 0 (in Fig. 6B) are six orders of magnitude greater than those computed from an equilibrium thermodynamic calculation, i.e., the flame chemistry is far from equilibrium. However, the constant-volume explosion temperature of 1334~ will eventually be attained by subsequent recombination. The increasing density and temperature fields ahead of the flame (see Figs. 6A, C) are a direct consequence of confinement. Compression of the unburned mixture, trapped between the advancing flame and the fixed boundary at z = 1, is isentropic to within a fraction of one percent. Not unexpectedly, this compression of the reactants
1335
NUMERICAL PREDICTIONS will influence the rate of flame propagation. Flame speed predictions must be preceded by a criterion for flame location, which appears to have no unique definition. The present investigation chose an arbitrary value of temperature, 700~ in the high gradient region (see Fig. 6A) and then tracked the movement of this value in time and space. The resulting flame-front locus is plotted in Fig. 7A. Numerical time differentiation of this path produces the plot of flame velocity, as seen by an inertial observer, shown in Fig. 7B. The prediction is a pulsating laminar flame, whose velocity varies from approximately 300-600 cm/sec. Although the oscillation is a complicated function of several modes, the dominant mode is characterized by a slowly increasing frequency which can be predicted
.60
~ ~
:~ ~" ~S .50 5.40 ..
t
,
i
I
j
I
I
I
h
i
I
,
,
,
with an acoustic analysis originally developed by ]ostY Assuming stationary gases contained in a unit length closed chamber subdivided by a flame surface at z, where burned gas occupies the region 0 ---> z, Guenoche 3~ gives pbab tan
a,P
---pua,, tan La--'~] (16)
where P is the period of oscillation. For the flame movement shown in Fig. 6, the root of Eq. (16) for the period of the fundamental mode is within a few percent of the value predicted by the numerical flame tracking routine. Visualization of this mode is easier with a plot of the chamber velocity distribution (see Fig. 8A) during one-half cycle of the flame oscillation, i.e., the time period indicated by arrows in Fig. 7B. Superposed on Fig. 8A is the temperature distribution during the same time interval. A coupling between the velocities of the burned and unburned gas is involved here. When the flame velocity is minimum, the unburned gas is nearly stationary while the burned gas achieves its maximum velocity away from the flame (negative velocities are directed toward z = 0). When the flame velocity is maximum, the burned gas is at minimum velocity while the
,
y y!yV!
~410
390
~
~ 2.1 ~
~
t.t.i u-
380
350
'i
,
,
i 35
,
h
~
,
i 40
,
*
*
,
{ 45
. . . .
9
2.1,
5O
TIME I N - D )
F{o. 7. (A) Flame position (based on 700~ vs. time, (B) flame velocity (as seen by inertial observer) vs. time, and (C) flame velocity minus local gas velocity at flame location vs. time. Corresponds to results shown in Fig. 6.
,2
.4
DISTANCE,Z
.6
(N-O)
,8
1.0
Fie. 8. (A)Chamber velocity distribution, and (B) chamber pressure distribution during one-half cycle oscillation (indicated by arrows in Fig. 7B).
1336
COMBUSTION OSCILLATIONS
unburned gas achieves its maximum velocity away from the flame. The chamber pressure distributions for the identical time interval are shown in Fig. 8B. The behavior of the pressure profile on the burned gas side is a result of a sound speed which is more than twice the value in the unburned gas. The origin of the flame oscillation, in the present example, can be traced to the heating phase before ignition of the mixture. Rapid heat transfer at one boundary excites acoustic wave motion in the chamber as pressure increases are equilibrated, i.e., thermally-driven convection. Thus, small amplitude standing waves exist in the confined gas mixture before ignition. After ignition, these waves are amplified and driven by the flame as it propagates away from the wall boundary. In other problems where the time to ignition is on the order of the acoustic travel time in the chamber, the onset of oscillation is delayed until the flame propagates some distance into the chamber. Although not required in the present solution, there is considerable interest in the question of burning velocity, i,e., the flame velocity relative to the unburned gas velocity. A trivial concept in steady flame propagation, it becomes ambiguous in a transient problem where the unburned gases are in non-uniform motion. Without resolving this issue, it was decided to reference the flame velocity to the local gas velocity at the flame location point. The numerical solution predicts this local gas velocity is nearly inphase with the flame velocity, again confirming the acoustic origin of the oscillation. The difference between the two velocities is plotted in Fig. 7C for the flame path in Fig, 7A. Note that this velocity difference is positive, i.e., the flame always moves into the unburned gas. Since the flame-tracking routine has an estimated error bound of +__5cm/sec, it is not possible to make definitive statements about the fine structure. Qualitatively, however, the curve indicates an overall flame acceleration which is not evident in the flame velocity seen by an inertial observer (Fig. 7B). As the flame propagates through the chamber, the average pressure level continuously increases. During the time interval of Figs. 6 and 7, the mean chamber pressure varies from 1.5 to 2.2 atm. The pressurization rate when the flame is near midchamber is approximately 3 • 105 psi/sec, and increasing. The steady state laminar burning velocity in 18 mole % O~ in O~ at 1 atm is 12-24 cm/sec ~, depending upon the flame shape assumed, If the velocity difference discussed above is assigned the role of a burning velocity, then clearly the transient flame is propagating into the unburned mixture more than an order of magnitude faster than a steady state flame at the pressure levels involved. The explanation is based on the dynamics of the process, A direct consequence of confinement is compression of the unburned mixture ahead of the flame. Thus,
chemical conversion in the flame releases essentially a fixed amount of energy into a continuously increasing temperature and density field. This increases the rate of heat release, which in turn steepens the temperature profile in the flame (this is evident from a close inspection of Fig. 6A). The cumulative effect leads to an acceleration of the flame.
Concluding R e m a r k s A one-dimensional model, based on a numerical solution of the compressible Navier-Stokes Equations, was developed to study transient laminar flame propagation in confined premixed gases. The solution procedure properly accounts for multispecies diffusion, finite-rate chemical reactions, and interaction of the combustion wave with the confining geometry. In this initial study, using independentlydetermined values for the chemical reaction rates and other gas properties, the model demonstrated encouraging predictive capability, Obviously, much remains to be done. Although the present investigation examined only ozone/oxygen mixtures, other exothermic premixed systems should exhibit similar oscillation and acceleration behavior. Theoretically, because of the nonlinear coupling between the equations governing conservation of momentum and energy, confined combustion is inherently oscillatory. In practice, the limiting amplitudes may be negligible as the result of system losses, etc, However, any analysis of confined or partially confined combustion based on assumptions which preclude oscillatory behavior will lead to predictions with unknown error, Nomenclature
a*
reference speed of sound (cm/sec). , " i c heat at constant pressure. cvjT 9, r, z ; speclf I*,, binary diffusion coefficient beD~k D*~ /a*r ~ tween i th and k'" species. E ph - p + pu 2 /2; total energy per unit volume. h h*/a*2; enthalpy, ~=~ Y,h,. h~' heat of formation of the i 'h species at T ~ J total number of computational cells. k~ k~ T*/O~* a*~L*; thermal conductivity. L* reference length (cm). m pu, mass flux. ,,~', molecular weight of the i 'h species (g/g-mole). 1/(~,Y,/~r mixture molecular weight (g/g-mole). N total number of species. ,1 ,a, 2 P P /Or 9 ; pressure. P a * / f ~ L *; nondimensional period (f * = cyeles/seq). q,~ qw/p~a r ; external heat flux, q * = c a l /
cp,
em
-sec.
NUMERICAL PREDICTIONS Re
p* a* L * / I x * ; reference value of Reynolds number. Ro R*o T r* f a r ,2 ; universal gas constant. Ro E ~_, Y,/.~r gas constant. t t * a * / L * ; t i m e ( h e r e , L * / a * = 1.638• 10 -6 sec). tin,, cut-off time of external heat flux (See Fig. 2). T T * / T *; temperature. T* reference value of temperature (~ u u * / a *~; velocity. V~ V ~ / a * ; diffusion velocity of i 'h species. W, Y,V,. X, mole fraction of i th species. Y~ mass fraction of i th species. z z * / L*; distance coordinate. a At/2Az, OR. Ix Ix */Ix *; coefficient of viscosity. Ix* reference value of coefficient of viscosity (g/cm-sec). I~, (4/3) ~ / R e . p p * / p * ; density. p* reference value of density (g/cm:~), 0 weight factor, O = 0.55 here. to 0~"L /p~ a~ ; chemical production rate for species i. Superscripts * n
dimensional quantity. time level,
Subscripts i i ( )~ (), (), ( )~
pertaining to the i th species, grid point location. partial derivative with respect to z. partial derivative with respect to t. u n b u r n e d gas. b u r n e d gas. Acknowledgment
The author gratefully acknowledges the substantial contribution of Mr. R. D, Anderson who executed the computer program and constructed the CALC O M P plots.
REFERENCES I. STRENG, A. G. ANt~GROSSE, A. V.: Sixth Symposium (International) on Combustion, p, 264, Reinhold, 1957. 2. MARXSTEIN, G. H.: Fourth S y m p o s i u m (International) on Combustion, p. 44, Williams and Wilkins, 1953. 3. MAR~;STEIN, G. H. AND SOMERS, L. H,: ibid, p. 527.
1337
4. JosT, W., KauG, J., AND S1EG, L.: ibid, p. 535. 5. KASrAN, W. E.: ibid, p. 575. 6. MARKSTEIN,G. H. (Eo.): Non-Steady Flame Propagation, Chapters D, G, MacMillan, 1964. 7. JONES, H.: Proc. R. Soe. Load. A, 353, p. 459 (1977). 8, BUTLER, T. D. AND O'RoURKE, P. J.: Sixteenth S y m p o s i u m (International) on Combustion, p. 1503, The Combustion Institute, 1977. 9. WESTBROOK,C. K. AND CHANG, J. S.: "A Theoretical Study of Flame Propagation Through Stratified Media," Lawrence Livermore Laboratory Report UCRL-78925, 1977. 10. KOTrlARI, A. P., ANDERSON, J. D., JR., AND JONES, E.: AIAA J,, 15, p. 92 (1977). 11. CRAMAROSSA,F. ANDDIXON-LEwIS, G.: Comb. and Flame, 16, p. 243 (1971). 12. RWARO, W. C., FARMER, O. A., AND BUTLER, T. D.: "'RICE: A Computer Program for Multicomponent Chemically Reactive Flows at All Speeds," Los Alamos Scientific Laboratory Report LA-5812, 1975. 13. KEE, R. J. AND MILLER, J. A.: "A Split-Operator, Finite Difference Solution for Axisymmetrie Laminar Jet Diffusion Flames," Sandia Livermore Laboratories, SAND77-8502, 1977. 14. HIRSCHFELDER, J. O., CURTISS, C. F., AND BIRD, R. B.: Molecular Theory o f Gases and Liquids, Chapter 7, John Wiley, 1954. 15. EDELMAN, R. B., FORTUNE, O. F., WEILERSTEIN, G., COCHRAN, T. H., AND HAGGARD, J. B., JR.: Fourteenth S y m p o s i u m (International) on Combustion, p. 399, The Combustion Institute, 1973. 16. KoogE~, D. E.: "Numerical Predictions for Ignition and Flamespreading in a Confined O z o n e / O x y g e n Mixture," 14th J A N N A F Combustion Meeting, CPIA Publication 292, p. 335, Vol. II, 1977, 17. BEA~I, R. M. aND WARMING, R. F.: AIAA J., 16, p. 393 (1978), See also Warming, R. F. and Beam, R. M.: SIAM-AMS Proceedings, 11, 1977. 18. JOHNSTON, H. S.: "'Gas Phase Reaction Kinetics of Neutral Oxygen Species," National Bureau of Standards, NSRDS-NBS-20, 1968. 19. BAULCH, D. L., DRYSDALE, n . D., HORNE, D. G., AND LLOYD, A. C.: Evaluated Kinetic Data for High Temperature Reactions, Butterworths, 1973. 20. GORDON, S. AND McBRIDE, B. J.: "Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks, and Chapman-Jouget Detonations," NASA/ Lewis, NASA-SP-273, 1971 (1976 program version). 21. HEIMERL, J. M. AND COVVEE, T. P.: Private Communication 1978 (to be published). 22. H1NOMARSH, A. C.: "'Solution of Block-Tridiagonal Systems of Linear Algebraic Equa-
1338
23. 24.
25.
26. 27. 28. 29. 30.
COMBUSTION OSCILLATIONS
tions," Lawrence Livermore Laboratory, UCID30150, 1977. YANENKO,N. N.: The MethodofFractional Steps, Springer-Verlag, 1971. DwYER, H.: Private communication Dec 1977, See also Dwyer, H. and Otey, G.: "A Numerical Study of the Interaction of Fast Chemistry and Diffusion," AIAA Preprint 78-946, 1978. HINDMARSH,A. C. ANOBYRNE,G. D.: "'EPISODE: An Experimental Package for the Integration of Systems of Ordinary Differential Equations," Lawrence Livermore Laboratory, UCRL-75868, 1975. STRENG, A. G. AND GROSS, A. V.: Comb. and Flame, 5, p. 81 (1961) (See Fig. 2). HIRSCHFELDER,]. O., CURTISS, C. F., AND CAMeBELL, D. E.: J. Phys. Chem., 57, p. 403 (1953). MARGOLIS, S. B.: J. Comp. Phys., 27, p. 410 (1978). ]OST, W.: Explosion and Combustion Processes in Gases, p. 103, McGraw-Hill, 1946. GUENOCHE, H.: "Flame Propagation in Tubes
31. 32. 33.
34. 35.
36.
and Closed Vessels," p. 149 in Nonsteady Flame Propagation, G. H. Markstein (ed.), MacMillan, 1964. MARaEaO,T. R. ANDMASON,E. A.: J. Phys., Chem. Ref. Data, 1, No. I (1972). HANLEr, H. ]. M. AND ELY, ]. F.: J. Phys. Chem. Ref. Data, 2, No. 4 (1973). TRETTEa, V. J., JR.: "'An Experimental Determination of the Viscosity of Gaseous Ozone," M.S. Thesis, Chem. Engineering, University of Idaho, 1966. DALGARNO,A. AND SMITH, F. J.: Planetary Space Sci., 9, p. 1 (1962). TOULOUKIAN,Y. S., L1LEY, P. E., AND SAXENA,S. C.: Thermophysical Properties of Matter, Vol. 3; IFF-Plenum, 1970. EDELMAN,R. B., FORTUNE, O., AND WEILERSTEIN, G.: "Analytical Study of Gravity Effects on Laminar Diffusion Flames," General Applied Science Laboratory, GASL TR-771, 1972 (also NASA CR-120921).
COMMENTS B. T. Zinn, Georgia Institute of Technology, USA. Your results indicate the formation of a shock wave, and I wonder if you have identified one and have you followed its evolution with time? Also, do your results indicate the presence of the fundamental and higher harmonics of the acoustic modes of the chamber? Also, in previous combustion oscillation studies the onset of oscillations has been related to the presence of some physical/kinetic mechanism whose characteristic time is of the same order of magnitude as the period of oscillation of one of the combustor's modes, Can you identify such a process for or a physical-acoustic interaction from your results? Finally, will you comment upon the length of your computations and the major contributors to this time.
Authors'Reply. O f t h e c o m p u t a t i o n s m a d e t o date, only mixtures containing 70 mole percent 0 3 / 0 2 and higher indicate the presence of a shock wave, and this occurs at ignition. The gas velocity ahead of the forming shock wave quickly exceeds the maximum value dictated by a cell Reynolds number equal to two. Furthermore, without a prohibitively fine grid mesh system, the numerical solution is incapable of resolving the interior structure of the shock wave. Momentarily ignoring these difficulties leads to a steep-fronted wave containing oscillations
on the high-pressure side. Although tolerable in a non-reacting inviscid "shock-capturing" solution, these numerical oscillations quickly destroy the present calculation w h e n coupled into the finite-rate chemistry package. Several frequencies are evident in the plot of flame velocity versus time (see Fig. 7B). As shown in the paper, the dominant mode is the fundamental acoustic mode of the chamber, assuming the chamber is subdivided by a flame surface separating burned and unburned gases. The first harmonic of this mode is also present, until the flame passes midehamber. These modes are not stationary in the usual acoustic sense, since their frequencies change with the advancing flame. Most likely, the combustion oscillation studies referred to in your question are the instability analyses of solid and liquid rocket engines. Note that both of these analyses include a basic combustion mechanism which is independent of the flow field in which oscillations occur. The solid propellant regression rate is customarily assumed to be an exponential function of surface temperature which is controlled by the time-dependent thermal wave in the u n b u r n e d solid; hence, there is a characteristic relaxation time which is only a function of the propellant. In liquid combustion instability analyses, a host of unknowns involving droplet vaporization, mixing, and chemical kinetics are often lumped into two parameters called " n " and '"r";
NUMERICAL PREDICTIONS both are specified independently of the rocket engine flow field. Then the question is posed whether these driving sources can be "in-phase" or "out-of-phase" with the natural acoustic modes of the confining geometry. The present problem is different in the following sense. The driving source for the premixed flame is the heat release distribution which follows directly from the chemical production rate terms; these are sensitive functions of density and temperature in the local flow field. What can and does happen is that the flow field oscillations alter the heat release
1339
distribution. This coupling mechanism allows considerable latitude in which modes can be driven. O f the cases run to date, it appears that Rayleighs' criterion is satisfied by an easily accessible mode and hence the flame oscillation is sustained even though the frequency may change. The present computations require over 3 minu t e s / 1000 steps on a CDC 7600. Depending on flame speed, between 8000-12,000 steps are necessary to traverse the chamber. More than two-thirds of this time is consumed by the stiff integrator chemistry package.