Calculation of the structure and extinction limit of a methane-air counterflow diffusion flame in the forward stagnation region of a porous cylinder

Twentieth Symposium (International) on Combustion/The Combustion Institute, 1984/pp. 1893-1904

CALCULATION OF T H E STRUCTURE A N D EXTINCTION LIMIT OF A METHANE-AIR C O U N T E R F L O W D I F F U S I O N F L A M E IN THE FORWARD STAGNATION REGION OF A POROUS CYLINDER G. DIXON-LEWIS, T. DAVID, P. H. GASKELL

Department of Fuel and Energy, University of Leeds, Leeds LS2 9JT, U.K. S. FUKUTANI, H. JINNO

Department of Industrial Chemistry, Faculty of Engineering, University of Kyoto, Japan J. A. MILLER, R. J. KEE, M. D. SMOOKE

Sandia National Laboratories, Livermore, CA 94550 N. PETERS, E. EFFELSBERG

Institut fiir AUgemeine Mechanik, R.W.T.H., Aachen, Germany J. WARNATZ, F. BEHRENDT

Institut fiir Angewandte Physikalische Chemie, Universitiit Heidelberg, Germany The paper addresses itself to the numerical representation of the structure of counterflow methane-air diffusion flames, with the use of complex chemistry and detailed formulation of the transport fluxes. A similarity solution is briefly outlined which allows the problem to be treated as one-dimensional in space. Following this, discretization and solution of the resulting algebraic equations has been achieved by five different methods. Although there are some differences between the results of the separate computations, these results agree within the probable experimental error with the observations of Tsuji and Yamaoka as regards the temperature and species mole fraction profiles in a particular flame. The more detailed chemical models are able to predict the profiles of C2 hydrocarbons with reasonable precision; but even somewhat less detailed models are able satisfactorily to predict the major structural features. Despite this agreement, the system overall does not behave as a straightforward boundary layer flow, and in order fully to match the solution to the experimental results it was necessary to modify the measured velocity gradient of 100 s -1 for the cold flow to near 130 s-1 in the flame region. The discrepancy is probably caused by modification of the pressure field due to apparatus effects or to the flame itself, and a full solution requires a twodimensional treatment to include this effect. As the velocity gradient which characterises these flames is increased, the increased strain in the flame reaction zone (or flame stretch) may lead to extinction. The effect of increasing strain on the flame structure and position near the porous cylinder is examined. There are differences between the two available predictions of the quenching velocity gradient, but both predictions are consistent with the measured, nominal gradient. It is shown that thermal quenching effects at the cylinder do not contribute appreciably to the observed limits. The significance of the extinction limits in the context of non-premixed turbulent combustion is briefly discussed.

1. Introduction During the past thirty years or so counterflow diffusion flames have been widely used experimentally both for the determination of the chemically controlled extinction limits of diffusion flames, and for the measurement of diffusion flame structure. The historical development has recently been excellently reviewed by Tsuji./~/The present paper is

concerned with the detailed numerical representation of such flames with complex chemistry and detailed transport properties. For this purpose the methane-air counterflow diffusion flames studied experimentally by Tsuji and Tamaoka(2-~1 have been chosen for investigation, since both extinction limits and detailed profiles of temperature and species concentrations have been measured. Physically, methane is blown from a porous cylinder into an 1893

1894

LAMINAR FLAMES

oncoming air stream, so that a free stagnation line parallel with the cylinder axis builds up in front of the porous surface. This line is in the symmetry plane of a stagnation point boundary layer, and the combustion takes place within a thin flame zone around the location of the stoichiometric mixture. The overall flow configuration is shown in Fig. 1. The location of the flame zone is within the boundary layer but outside the stagnation point, and there is a convective backflow of the combustion products into the fuel-rich region of the flame. (~) The paper summarizes results obtained by several groups working independently. The major part of the results was initially discussed at a GAMM Workshop entitled "Numerical Calculation of Nonequilibrium Diffusion Flamelet Structure," organized by J. Warnatz and N. Peters, and held at the University of Heidelberg on 1-3 December, 1983. 2. Equations to be Solved

tum, energy, and chemical species in stationary, laminar, plane stagnation point flow are:

a/ax(pu) + O/ay(pv) = 0

p(u(Ou/ax) + v(Ou/O#

v

p{u(Oh/Ox) + v(Oh/ay)} = u(dp/dx) - OQ/dy

'

(1)

(2)

Ue = ax

in the free stream; but the second equation requires modification everywhere. The boundary, layer continuity equations for mass, x-direction momen-

,',, \ "~,, 'V boundary layer edge ~1'~

I,

J

IX

02, N2

[.]~,,,.11

|

p{u(O(ri/Ox ) + v(O(rJOy)} = -O/Oy(ji/mi) + Ri (i = 1,2 . . . . . N)

(6)

where p is the density, IX is the viscosity of the mixture, h is its specific enthalpy, cri = Yi/mi where Yi and mi are mass fraction and molecular weight respectively of the i-th component, and Ri is the molar chemical production rate of the i-th species within the control volume. The specific enthalpy h, the pressure p, the energy transport flux Q and the diffusive fluxes Ji are given respectively by: N

h = ~ (r

(7)

i=1

p = pRT~((r,)

(8)

/

Q = -

'yh(oh/c)y) +

flame

zone

jl = -

13~h(oh/OY) +

(9)

[3umj(Offj/OY)

(10)

j=l

where the coefflcients' ,~h, ,~j,(6)131banandH 13~ are derived from gas kinetic theory, d i is the molar enthalpy of the i-th species. The H~ are expressed as polynomial functions of the temperature T. In the following, the contributions of the pressure gradient and viscous terms in the enthalpy equation (5) will be neglected. It can also be shown (7'8) that the symmetry of the system is such that all the dependent variables except the velocity u have zero x-direction derivatives at the stagnation point. It is noted at this stage that Eq. (3) is identically satisfied if the following relations hold:

y

pl, t

peXS'(y)

pv = -Pef(Y).

(11)

(12)

Next, by introducing the similarity variable xI such that: n =

FIG. 1. F l o w configuration.

"yjm~(&rj/Oy j=l

'

/J"/!,; .'~ f ree stagnation /'/ / point ,/'

(5)

}

-a(U-Uo)

where a is the velocity gradient. For viscous, variable density flow around a boundary layer, the first of these equations still remains valid in the free stream region of the flow (subscript e), that is:

CH4

-clpldx + alay(ix(aulay)) (4)

=

If u and v represent the tangential (x) and radial (y) direction velocities in the neighbourhood of the stagnation point at x = 0, y = Yo, then for frictionless, steady state, plane potential flow of an incompressible fluid the velocity distribution in this neighbourhood is given by:

u = ax

(3)

L

(a/Ve) 1/2" ( p / p e ) d y

(13)

METHANE-AIR COUNTERFLOW DIFFUSION FLAMES where v = tx/p is the kinematic viscosity, and putting: f(y) = (av~) 1/2 . ~b(-q),

(~,)2 _ ~b~" = Pe/P q- d/d~(Cd~"),

3. Methods of Solution and Input Data (15)

where C = (plx)/(p}x)~. Equations (11) and (12) then lead respectively to: Cb' = UlUe -

(16)

pv l{(plx)ea} 1/~,

(17)

and finally introducing V = -~b, Eq. (2) to (5) transform into the o.d.e's (18) to (21): d V / d ~ + d~' --- 0

(18)

v(dC'/dn) = dlld'l]{C(ddp'/dlq)} -}- Pe/P

appropriate Eq. (20) or (21), and do represents the source term Ro/p. These forms will be encountered below.

(14)

the momentum equation (4) transforms into:

d~ =

1895

- -

(~b') 2

V(dh/drl) = -{(plx)~a}-I/2 dQ/d~l

(19)

(20)

V(dcri/d~) = -{(plx)ea} -I/2 d/d'q(ji/mi) + Ri/ap (~ = 1,2 . . . . . N)

(21)

with boundary conditions (22) at the wall (subscript w), and (23) at the outer edge of the free stream: aq = 0, d~' = 0, V = Vw, T = Tw, Vw(Gi/mi - ~i) = {(p~)ea} -1/2" (ji/mi)w

(22)

71 = 'l']e, ~b' = 1, T = Te, cri = cri,~,

The system of differential equations (20) and (21) is usually solved by some form of finite difference procedure; and since the high temperature oxidation of methane in flames involves upwards of forty or so elementary chemical steps with characteristic times differing by several orders of magnitude, it becomes necessary to use an implicit approach. For a system with N chemical components this leads at each spatial node to (N + 1) non-linear algebraic equations in (N + 1) unknowns. The equations are linearized about guessed initial values of all the dependent variables (or a previous solution), and the resulting set of linear equations is solved to give the corrections ~Ol~. However, because of the form of the transport flux terms Q and ji (Eq. (9) and (10)), the set of equations at each node is coupled by the gradient terms with the sets at neighbouring nodes to give, at its most formidable, a complete system having a tridiagonal block matrix of left-hand side coefficients, and covering the whole flame. Uncoupling from neighbouring nodes is accomplished by a Ganssian elimination procedure,/9/but an inversion of an (N + 1) x (N + 1) matrix is still required at each node. The whole process of calculation of ~O1~ is repeated until convergence is achieved, though such convergence is not, of course, guaranteed independently of the starting profiles. Two common approaches to the iterative solution are (a) Newton's method, and (b) the time-dependent approach. In the former,/m'lt,l~/an initial solution estimate O m to the set of equations:

zero gradients of all dependent variables except V.

(23)

Here Gi denotes the mass fraction of the i-th component in the input stream. From the experimental data:/5) Gca, = 1.0, Gi,,ca4 = 0, cro~.e = 7.278 x 10-3, CrN2.e = 2.738 X 10-z, and all other boundary cri,e are zero. The temperatures at the boundaries are Tw = 332.16 K and Te = 283.i6 K. To conclude this section, it should be noted that, subject to the constraints which determine the steady state flow in the region of the stagnation point streamline, the time-dependent forms of the energy and species equations become, for the general dependent variable 0:

(lla)OOlOt + V(O0/Oxl) = -{(pl.t)ea} -1/2 acolO~ + dola

(24)

where co represents the transport flux term of the

f,(01,02,03, 9 .) = 0 leads to the equations:

corrections

(i = 1,2,-.,N) (~O) m+l

(25)

by way of the

j ( o m ) ( ~ o ) m + 1 = _ h m f ( o m)

(26)

where J(O m) = (OfilOOk) m for all i and k is the Jacobian matrix, and h m (0 < h -~< 1) is the m-th damping parameter, h = 1 corresponds with a full Newton step, and gives rapid convergence provided that O " is within an appropriate domain around the true solution, h < 1 increases the domain of convergence. For the present problem this method caused difficulty due to ill-conditioning of the Jacobian matrix. In the time dependent method the O/J+l)stln ap" propriate to Eq. (24) can be written as (OJ( 8 0 + 8OJ/St/)l~. With linearization of the chemical source term around its value at j(~t), the equations to be

1896

LAMINAR FLAMES

solved for the 5 0 become exactly the same as Eq. (26), except that h is always 1, and damping is introduced instead by increasing the diagonal elements only of the ]acobian matrix by (a.~t) -l. This increases the likelihood of convergence compared with the Newton method, but decreases the rate of convergence. The quantity (a.St) -1 is gradually decreased as the computation proceeds, and a minimum over-ride value may be imposed so as to preserve computational stability. The full, undamped Newton method is equivalent to the time-dependent approach with infinite time step. Both methods have the capability of adaptive gridding so as to concentrate the nodes at the positions of maximum activity. The specific approaches employed by the various groups on the present problem were as follow: A. F. Behrendt, J. Warnatz (Heidelberg): Solution of time-dependent equations for ~b', species and temperature, with V by numerical integration from V = Vw. Solution by finite difference method, (13) second order in space (assuming parabolic behaviour between three neighbouring grid points), and first order in time; off-diagonal elements in Jacobian matrix neglected (thus uncoupling the (N + 1) equations at each node). Thirty five grid points, nonuniformly distributed, with concentration in the hot part of the flame. First solution of the corresponding cold diffusion problem, then introduction of bellshaped profiles of temperature (Tm~ = 2200 K) and H atom mole fraction (Xmax = 0.02) to initiate reaction. Reaction mechanism as in ref. (14), with C3, Ca hydrocarbons removed. Transport model as described by Kee, Warnatz and Miller./15) B. T. David, G. Dixon-Lewis, P. H. Gaskell (Leeds): Solution by finite domain method, (6/ with micro-integration over computational cells on the assumption of linear variation of the dependent variables with xl between nodes. Enthalpy and species by block implicit solution of time-dependent equations, d~' by Newton method, V by numerical integration from V = Vw. Sixty grid points eventually in a = 100 s -1 flame, after start with 33 and adaptive placing of remainder in two stages in hotter part of flame. Initial profiles of similar shape to final result, but (owing to error in setting up) grossly too large values of lq. Bell-shaped initial H, O and OH profiles, with maximum equivalent to one or two percent H atoms by volume. Eighteen species, forty five reaction mechanism (t~'17) with C~ hydrocarbons truncated at ethylene stage. Multicomponent diffusion and thermal diffusion coefficients, and thermal conductivity, by detailed calculation. (6'1s) Viscosity of mixtures by Buddenberg and Wilke formulafl 9) Lennard-Jones potential parameters (mostly from Svehla(2~ given in ref (6), except that here ~rH = 0.30 nm. Studies on flames with higher velocity gradients were initiated from developed profiles of previously

converged solutions. Reproducible determination of extinction limits required further adaptive placing of nodes in the hotter region, so as to give adequate spatial resolution in and near the very narrow region of net radical production (see w C. S. Fukutani, H. Jinno (Kyoto): Finite difference method. Seventy nine grid points, with distribution proportional to the temperature gradient. Initial S-shaped distributions for CH4, 02 and Na, with Gaussian distribution for temperature. Sixty one reaction mechanism, with most of rate parameters from Westbrook. ~211Multicomponent diffusion velocities by Curtiss-Hirschfelder approximation. Thermal diffusion not included. D. J. A. Miller, R. J. Kee, M. D. Smooke (Sandia, Livermore): Finite differences to discretize the conservation equations; central differences for diffusive terms, upwind differences for convective terms. Solution by combination of (i) time-dependent method of lines approach using the computer code DASSL, a variable order, variable step size, implicit differential-algebraic equation solver written by Petzold,/a2) (ii) a first order, fixed step size, implicit time integration procedure, and (iii) a steady state, modified, damped Newton's method algorithm as used previously in premixed flame problems/w'11A2/ (see comments above). The time-de~ pendent approaches were used to bring the solution into the domain of convergence of Newton's method. Sixty six grid points, with adaptive placing, mostly concentrated in the regions of high temperature. Two reaction mechanisms employed: (a) a "full" mechanism, including Ca species; and (b) a "short" mechanism, excluding all reactions involving CH2, CH, and all C2 and higher hydrocarbons. The "short" mechanism is regarded as probably sufficient to predict most of the gross flame properties, though not of course the formation of higher hydrocarbons. Multicomponent diffusion velocities by CurtissHirschfelder approximation. Thermal diffusion not included. Mixture viscosities by Wilke's formula. Stockmayer potentials from Kee et al. 051 E. N. Peters, E. Effelsberg (Aachen): Hermitian finite difference method./z3/This is a fourth order implicit scheme that leads to a tridiagonal block matrix. Solution by Newton iteration. Thirty six equally spaced grid points, with ~1 = 0.125. Initial profiles were inserted on the assumption of constant enthalpy, constant pressure equilibrium across the mixing zone. (24) Forty eight reaction mechanism, including C2 hydrocarbons. Multicomponent diffusion and thermal diffusion computed as in ref.

(25). 4. Results

(a) Flame with Velocity Gradient a = 100 s -1 This flame was studied initially since the extensive experimental observations of Tsuji and Ya-

,oo

METHANE-AIR C O U N T E R F L O W DIFFUSION FLAMES

1897

N2

t- 10-I o

18

02

"S

v O

T'-

o E

-i

0-~

I

~10 10-40

E

2

4

6

y/(mrn) FIG. 4. Mole fraction profiles of major species in flame with a = 100 s-l, comparing computed lines with points representing observations of Tsuji and Yamaoka. 15) Computations by group B. Points: 0, CH4; O, N2; &, Oz; V, H20; ~7, CO2; A, CO; II,

2 0

3 y/(mm)

H~.

6

FIC. 2. T e m p e r a t u r e profiles. Points represent observations of Tsuji and Yamaoka 15/ on flame with nominal a = 100 s -t. Lines computed by group B: line 1, a = 100 s -l, Vw = 1.5; line 2, a = 100 s -~, V,, = 2.45; line 3, a = 130 s -l, V,. = 2.15; line 4, a = 320 s -l, Vw = 1.5.

-60

I

,~i, ~ 3 ~,~_.,~,

i

maoka(5) are nominally available for comparison. The initial computations were performed for a dimensionless methane blowing velocity Vw = 1.5. The temperature profile, the velocity profile and the major species mole fraction profiles across the flame, as computed by group B, are shown by the indicated curves in Fig. 2 to 4. Figure 5 shows the molecular hydrogen and C2 hydrocarbon profiles as

-

10~I

~

J

10-~-

-20

0

~/. /

-

-

~ 10-2

1"

2 0 ~ 1 0

2

I 4 y / (mm)

' 6

FIG. 3, Velocity profiles, Points represent observations of Tsuji and Yamaoka (s) on flame with nominal a = 100 s -t. Lines computed by group B: line 1, a = 100 s -I, Vw = 1.5; line 2, a = 100 s -l, V~, = 2.45; line 3, a = 130 s -l, Vw = 2.15.

y/(mm) FIC. 5. Mole fraction profiles of molecular hydrogen and Cz hydrocarbons in flames with a = 100 s -l, comparing computed lines with points representing observations of Tsuji and YamaokaJ 5) Computations by group D. Points: II, H2; 0 , C2H2; 0 , CzH4; O, C2H6.

1898

LAMINAR FLAMES

I

1 0 -2

t- 1 0 -3

.9

the radical mole fraction profiles from group B. There are no measurements with which to compare these. The profiles obtained by all five groups are close to those shown in Figs. 2 to 5. A comparison of the major features of these computational results is available from the first five columns of Table I, in which the maximum temperatures and the maximum mole fractions of several individual species are compared. The following features emerge:

-

_o10 -4 o

E10-5 0

/ 2

y/(mm)

4

FIG. 6. Computed mole fraction profiles of CH20 and radicals in flame with a = 100 s-~. Computations by group B.

computed by group D. In each case the computed profiles are compared with the experimental observations. Figure 6 shows the computed CH20 and

(i) The computed peak temperature is in all cases higher (by between 100 and 250 K) than that presented in Tsuji and Yamaoka's paper./~ Much of this discrepancy is due to lack of a radiation correction to the reported thermocouple measurements. For a 0.2 mm bead, calculation suggests a correction AT equal to 1440 ~, where ~ is the emissivity of the thermocouple material. More disturbing is the lack of agreement between the computed peak temperatures themselves, particularly since Miller et al (group D) found that changing from a "full" to a "short" mechanism only led to a 20 K difference in the peak value. This aspect needs further investigation. (ii) The physical distance from the wall at which the maximum temperature occurs is in good agree-

TABLE I Computed maximum temperatures and species mole fractions in counterflow methane-air diffusion flames Group

a/s i TmJK At y/mm

A

B

100 1879 3.42

100 1994 3.52

Max. mole fractions: 10 H20 102 CO2 102 CO 102 H2 103 H 103 O 103 OH 103 CH3 103 C2H2 103 Call4 104 C2H6 104 CH20 103 HO2 106 CHO

1.83 7.58 4.56 3.65 4.75 3.13 5.84 2.34 9.92 3.63 5.58 3.26 5.30 1.27

1.72 6.79 3.94 2.14 3.71 3.11 6.48 7.55 -0.66 42.7 2.97 6.31 16.1

C 100 2047 -1.79 6.05 4.09 2.54 2.74 2.59 6.24 5.35 2.75 5.72 11.6 1.66 -17.6

D* 100 1991 3.45 1.71 6.56 4.71 3.23 4.29 3.14 6.29 3.25 3.71 1.64 9.98 5.99 3.37 8.09

E 100 1904 -1.83 6.8 3.36 2.29 3.25 2.7 5.97 1.56 8.85 2.35 16.5 1.77 5.7 0.38

B

B

320 1864 1.88

410 1769 1.65

1.65 5.35 4.53 2.03 3.96 3.46 5.26 7.69 -1.32 41.8 8.63 10.4 39.5

1.57 4.73 4.68 1.73 3.21 2.96 4.02 6.10 -1.37 36.2 11.8 11.6 53.7

*This peak temperature and the mole fractions refer to the "full" reaction mechanism. For the "short" mechanism (see Sec. 3) Tm~ = 2011 K.

METHANE-AIR COUNTERFLOW DIFFUSION FLAMES ment with experiment in those cases where it is recorded (Table I and Fig. 2). (iii) The appropriate computed velocity profile in Fig. 3 is in fair agreement with observation, but shows a lower velocity gradient at large distances from the wall. (iv) For those mechanisms (in particular groups A, C and D) which attempt realistically to model the C2 chemistry, all the species profiles are in fair agreement with experiment. With the less precise chemical model of group B, unrealistically high concentrations of C2H6 (and probably CH3) are obtained, combined with low C~H4 and zero CzH2. The computed maximum molecular hydrogen may also be somewhat on the low side. However, the overall level of agreement even in this case confirms the finding of group D that the Cz chemistry is not essential for the prediction of the gross flame properties. (v) An additional feature of the species profiles, which confirms an observation made some years ago by Gordon, Smith and NcNesby, (2~) is the breakthrough of molecular oxygen into the fuel-rich side of the flame (Fig. 4). Such breakthrough is also predicted by the computations, but not to the extent observed experimentally. Curiously, the predictions when C2 chemistry is absent appear here to be better than when it is included.

(b) Effects of Changes in Methane Blowing

Velocity and in the Velocity Gradient The basis for the non-dimensionalization used by Tsuji and Yamaoka(5) in the presentation of their experimental results is not the same as that employed here. As a consequence, direct comparability of the computational results with their experiments demands the use, in combination with a = 100 s -I of a dimensionless velocity V,~. = 2.48 at the wali boundarv. The computed temperature and velocity profiles for a = 100 s- 1 and V,~ = 2.45 are shown as curves 2 in Figs. 2 and 3. Comparison with cmwes 1 (Vw = 1.5) shows that essentially the effect of the increased wall velocity is to move the stagnation point, the flame and the far flow field bodily away from the cylinder, with the structure of the combined system being virtually undisturbed. On the other hand, the computed and measured profiles are now no longer aligned on the physical distance scale. The straightforward similarity solution breaks down. In order to restore alignment whilst remaining within the constraints of the similarity solution, it becomes necessary to increase the velocity gradient in the flame region, and at the same time to reduce Vw somewhat so as to preserve the correct methane blowing velocity. The changes are related by Eq. (17). The measured velocity profile in the region to

1899

the right of Fig. 3 suggests a modified velocity gradient near 130 s -I. The curves 3 in Figs. 2 and 3 show the temperature and velocity profiles computed with a = 130 s -1 and the corresponding Vw = 2.15. The alignment of the predicted and measured profiles on the physical distance scale is restored, and the magnitudes of the velocities in the near flow field simultaneously agree better than before with the measurements. It must be concluded that, in relation to the similarity solution in the stagnation region, the flame and the flow field near it behave approximately as if the velocity gradient there were near 130 s -1 instead of the nominal 100 s -1. The nominal, or measured, value is that for potential flow of an inviscid, incompressible fluid with free stream velocity v past a transversely oriented cylinder of radius R, and is given by 2 v/R. The change from the nominal value is probably largely associated with pressure gradient effects inherent in the specific combustion chamber geometry. For the particular experiment under discussion the chamber consisted of a rectangular air channel of 3 cm • 12 cm cross-section, into which was inserted symmetrically a porous cylinder 6 cm in diameter and 3 cm in length. The additional pressure gradient effects may .thus arise from the partial blockage of the air flow channel by the cylinder. (2/ Further, there may be pressure gradient effects due to the flame itself. In either case the effect invalidates the boundary layer approach as a means of obtaining a complete solution, and necessitates a full treatment of the two-dimensional pressure field. The computed maximum temperature in the flame with a = 130 s -1 is 1970 K. The predicted species profiles become slightly narrower than when a = 100 s-I, but are not otherwise significantly changed.

(c) Extinction Limit of Flame The velocity gradient a may be regarded as the inverse of a characteristic residence time in the reaction zone of the stagnation point diffusion flame. As a is increased, the maximum temperature in the flame drops until eventually the flame is extinguished. This extinction limit occurs at the point of vertical tangency on the upper branch of the wellknown S-shaped curve of peak temperature versus DamkShler number due to Fendell (e7~ and Lifian.(z8~ Experimentally the extinction limit for the methane-air flames occurs at a nominal velocity gradient (see w of 330 s -1 for a cylinder diameter of 6 cm, (41 and at 375 s -1 for the smaller diameter of 3 cm. (1) In view of the comments of the preceding section, the higher result is likely to be the more meaningful value with which to compare computational predictions from the similarity solution. The computations of Miller et al (group D) gave

1900

LAMINAR FLAMES

extinction at a -> 350 s -1. With their "short" mechanism, and with a grid spacing ~'q = 0.025 in the reaction zone, they compute a stable flame at a = 330 s -1, and for a = 100, 200, 300, 310 and 330 s-1 the peak temperatures were 2011, 1945, 1870, 1860 and 1825 K respectively. Again with the "short" mechanism, but with ~q = 0.00625, the a = 330 s -1 case extinguished. With this higher resolution the maximum temperatures were 2009, 1939, 1836 and 1796 K for the remaining four flames. Initial exploratory computations by David et al (group B) used only 33 grid points, with ~-q = 0.144 near Tmax. These also led to extinction at a -> 350 s-1, and gave a stable flame at a = 300 s -1. However, with ~aq = 0.072 near Tmax (48 nodes), stable solutions were obtained not only for a = 320 and 330 s -1, but also for 350 and 380 s-i; and with 8"q = 0.018 (72 nodes in total) further stable solutions were obtained for a = 400 and 410 s -1, but not for a = 415 s -1. Clearly very high spatial resolution near Tmax is necessary for the reliable determination of the extinction limit. At a = 100, 300, 320, 380, 400 and 410 s -1 the computed peak temperatures were 1994, 1890, 1864, 1817, 1793 and 1769 K. All the computations in this section had Vw = 1.5. To give some further indication of the properties of the nearly quenched flames, the temperature profile computed by group B for the flame with a = 320 s -1 and Vw = 1.5 is shown as curve 4 in Fig. 2, and the maximum mole fractions for the flames with a = 320 and 410 s -1 are given in the final two columns of Table I. Figure 2 shows that the position of maximum temperature moves appreciably towards the porous surface at the higher velocity gradient. The maximum and minimum ydirection velocities in the flames are also found to increase compared with Fig. 3, to 85 and 28.6 cm s -1 respectively in the a = 320 s -1 flame. The levels of molecular oxygen breakthrough in the a = 320 and 410 s -1 flames are respectively about five and ten times the levels in the a = 100 s -1 flame. The difference between the computed extinction limits of groups B and D is clearly of considerable interest, and no firm explanation can be offered at this time. To retain consistency with a reaction mechanism previously calibrated on premixed methane-air flames, David et al used a somewhat low diffusion coefficient for atomic hydrogen, by setting the Lennard-Jones force constant CrH = 0.30 nm. A sensitivity test showed that if crH was reduced to 0.225 nm, thereby raising the H atom diffusion velocity by some 20 percent, the predicted extinction velocity gradient was reduced to b e tween 370 and 380 s -I. For crr~ = 0.225 nm, a--370 s-1 and Vw = 1.5 the predicted maximum temperature was 1781 K.

5. Discussion The structure of strained laminar flames, particularly in conditions which may lead to extinction, is of great interest not only in its own right, but also in relation to fluctuating laminar flamelet theories of turbulent burning./a9'3~ Theoretical treatments of non-premixed turbulent combustion yield the local scalar dissipation rate • as the most useful parameter relating flame stretch effects with turbulent mixing. • is defined by: x = 2D(OZ/Oxk) ~

(27)

where D is a diffusion coefficient and Z is the mixture fraction (a conserved scalar) to which the element mass fractions and the enthalpy are linearly related. X can be regarded as the inverse of a characteristic diffusion time. Within the counterflow flame geometry it is related to the velocity gradient a by: X = 2a(C/Pr)(dZ/d~l) 2

(28)

where Pr = Izc~/k. Peters/z9'3~ proposed the use of the measured quenching velocity gradients aq in the counterflow diffusion flame, in conjunction with the gradient (dZ/drl)st at the stoichiometric mixture fraction Zst, to define a quenching limit • q in general situations which may arise in turbulent combustion. The early part of the present paper suggests that there may be some practical difficulty in that the measured aq is only a nominal value. However, apart from this, Liew, Bray and Moss(31) have questioned the proposal since (i) by an alternative Lagrangian approach to the computational consideration of a transient diffusion flame they find much higher values of x~than are predicted by Peters' proposal, and (ii) they then argue that the value of • should be determined at the peak temperature (giving Xmax,q), and not at the stoichiometric composition. Their computation suggested that the two values Zst and Zm~x (at Tma~) become appreciably different at high rates of strain, with Tm~x moving towards richer compositions where Z and dZ/d~ are higher. They also argue that because of the closeness of approach of the nearly quenched flames to the porous surface, thermal quenching as well as aerodynamic quenching must become prominent in such flames. To quantify the role of thermal quenching, a further computational determination of the extinction velocity gradient in the methane-air flame was carried out by group B, but with the methane blowing velocity increased to give Vw = 7.0. For ~xI = 0.036 near Tmax, extinction occurred at a = 400 s-l; but stable flames were obtained at a = 375 s -1 with both ~1 = 0.036 and 0.072. For 8~1 = 0.018, a

METHANE-AIR COUNTERFLOW DIFFUSION FLAMES

1901

TABLE II Scalar dissipation rates in counterflow methane-air diffusion flames Group A

a/s -1 X,r/s- 1 Xm~/S-~ Zm~

100 6.1 7.4 --

B

C

D

E

100 100 100 4.3 4.43 3.25 5.5 --+ 0.3* 10.0 5.21 0.067 -+ 0.002 -0.071 -+ 0.002*

Mole fraction at Tm~: 10a O2 2.68 6.3 --- 1.0

5.76

4.56 -+ 1.0

B

100 320 4.75 12.3 15.9 19 --+ 1" 0.114 0.068 - 0.002 5.1

17.5 -+ 1.0

B 410 16.7 22.8 -- 0.5* 0.066 -----0.002 29~6 -+ 0.5

*These limits allow for uncertainty in the location of the peak temperature. Node separation in this region is 0.1 mm for 100 s-1 flame, 0.05 mm for 320 s -1 flame, and 0.023 mm for 410 s -l flame (group

B).

stable flame with Tm~x at a distance of 4 mm from terflow diffusion flame by means of a one-dimenthe surface was obtained for a = 400 s-l; but exsional stagnation point similarity solution, provided tinction occurred at a = 405 s -1. The extinction that the velocity gradient in the flame region itself occurs at a velocity gradient very slightly less than can be estimated. However, (ii) The flow field around such a flame is comwhen Vw = 1.5, and the contribution of thermal plex, so that this effective velocity gradient is not quenching in the latter case must therefore be relthe nominal gradient for a transversely oriented atively small. cylinder of radius R, which is given by 2 v/R, where On the assumption that the stoichiometric mixv is the free stream approach velocity. Because of ture is that which burns completely to CO2 and the pressure field induced by the combustion water, the stoichiometric methane-air mixture fracchamber geometry, or by the flame itself, a full sotion takes the value Zst = 0.055. On this assumplution cannot be obtained by use of the boundary tion, Table II gives values of • and • for the layer assumptions, and a two-dimensional represencounterflow diffusion flames with a = 100, 320 and tation of the flow field is required. 410 s -1, and Vw = 1.5. Although the values of • (iii) The extinction limits of the flames, as meaare of the same order as those estimated by Pesured by Tsuji and Tamaoka, are aerodynamic pheters, (a~ in no case is the difference between Xst and nomena, free of thermal quenching effects at the Xm~x, or even Zst and Zm~, very appreciable. It turns cylinder wall. They can thus be used for reasonable out for all the flames that Zm~x is very close to the estimates of the scalar dissipation rates correspondvalue Zst' = 0.071 which results from the assumping with aerodynamic quenching under adiabatic tion of stoichiometric burning to CO and water, and conditions. it is possible on this basis to discuss the structure of the flames in some detail (3z) in" terms of the gen-(28) eral diffusion flame model put forward by Lifian. When the final profiles of the a = 410 s -1, Vw = 1.5 flame are used as the starting point for a comREFERENCES putation of the extinction at a = 415 s-I, it is also found that Zm= does not exceed 0.07 until after the 1. TSUII, H.: Progr. in Energy and Combustion net radical production has ceased everywhere, that Science, 8, 93 (1982). is, until after the flame is already virtually extin2. Tsull, H. AND YAMAOKA,I.: Eleventh Sympoguished. sium (International) on Combustion, p. 979. The Combustion Institute, 1967. 3. TsuJI, H. AND YAMAOKA, I.: A gas dynamic 6. Conclusions analysis of the counterflow diffusion flame in The major conclusions from the present study are

thus: (1) It is possible satisfactorily to model the measured features of the structure of the present coun-

the forward stagnation region of a porous cylinder, Inst. Space and Aeronaut. Sci., Univ. Tokyo Rep. 404 (1966). 4. TSUJI, H. AND YAMAOKA,1. : Twelfth Symposium (International) on Combustion, p. 997. The

1902

LAMINAR FLAMES

Combustion Institute, 1969. 5. TsuII, H. AND YAMAOKA, I.: Thirteenth Symposium (International) on Combustion, p. 723. The Combustion Institute, 1971. 6. DIXON-LEwiS, G. : Computer modelling of com-

18. DlXON-LEWlS, G.: Proc. Roy. Soc. Lond. A 307, 111 (1968). 19. BUDDENBERG, J. W. AND WILKE, C. R.: Ind. Eng. Chem. 41, 1345 (1949). 20. SVErlLA, R. A. : Estimated viscosities and ther-

bustion reactions in flowing systems with transport, In "Combustion Chemistry" (W. C.

mal conductivities of gases at high temperatures. NASA Tech. Rep. R-132, 1962.

Gardiner Jr., Ed.), Springer-Verlag, in press. 7. LEES, L. : Jet Propulsion, 26, 259 (1956). 8. FAr, J. A. AND RIDDELL, F. R.: J. Aeronaut. Sci., 25, 73 (1958). 9. RICHTMVER, R. D. AND MORTON, K. W.: Difference Methods for Initial Value Problems, 2nd. ed. Interscience Publishers, 1967. 10. SMOOKE, M. D., MILLER, J. A. AND KEE, R. J.:

21. WESTBROOK, C. K.: Combust. Flame, 46, 191 (1982). 22. PETZOLD, L. R.: Proc. IMACS World Conference, Montreal, and Sandia National Laboratories Rept. SAND 82-8637. 23. PETERS, N.: Proc. Fourth International Conference on Numerical Methods in Fluid Dynamics. In Lecture Notes in Physics, 35, 313 (1975). 24. PETERS, N.: Berechnung einer Methan-Luft-

Numerical solution of burner-stabilized laminar premixed flames by an efficient boundary value method. In "Numerical Methods in Lam-

11. 12. 13. 14. 15.

inar Flame Propagation" (N. Peters and J. Warnatz, Eds.), Friedr. Vieweg & Sohn, Wiesbaden, 1982. SMOOKE, M. D.: J. Comput. Phys. 48, 72 (1982). SMOOKE, M. D., MILLER, J. A. AND KEE, R. J.: Comb. Science and Tech. 34, 79 (1983). WARNATZ,J.: Ber. Bunsenges. Phys. Chem. 82, 193 (1978). WPmNATZ,J.: Comb. Science and Tech. 34, 177, (1983). KEE, R. J., MILLER, J. A. AND WARNATZ, J.: A

Fortran computer code package for the evaluation of gas phase viscosities, conductivities and diffusion coefficients. Sandia National Laboratories Report SAND83-8209, 1983. 16. DIXON-LEwlS, G.: First Specialists' Meeting (International) of the Combustion Institute, p. 284. Section Fran~aise du Combustion Institute, 1981. 17. DIXON-LEwiS, G. AND ISLAM, S. M.: Nineteenth Symposium (International) on Combustion, p. 283. The Combustion Institute, 1982.

Diffusionsflamme im 6rtlichen Gleichgewicht und im Nichtgleichgewicht. VDI Berichte Nr. 246, 1975. 25. PETERS, N.: Int. J. Heat Mass Transfer, 19, 385 (1976). 26. GORDON, A. S., SMITH, S. a. AND NcNEsBY, J. R.: S e v e n t h Symposium (International) on Combustion, p. 317. Butterworths, 1959. 27. FENDELL, F. E.: J. Fluid Mech. 21, 291 (1965). 28. LIIqAN, A.: Acta Astronautica, 1, 1007 (1974). 29. PETERS, N.: Combs. Science and Tech. 30, 1 (1983). 30. PETERS, N.: Laminar diffusion flamelet models in non-premixed turbulent combustion. To appear in Progr. in Energy and Combustion Science, 1984. 31. LIEW, S. K., BraY, K. N. C. AND MOSS, J. B.:

The predicted structure of stretched and unstretched methane-air diffusion flames. Paper presented at the 9th International Colloquium on Dynamics of Explosions and Reactive Systems, Poitiers, France, July, 1983. 32. DAVID, T., DIXON-LEWXS, G. AND GASKELL, P. H. : In course of publication.

COMMENTS R. Bilger, University of Sidney, Australia. My congratulations to the various teams, the organizers of the GAMM workshop on a most significant effort. I hope these efforts will be continued as an understanding of this flame is essential to formulation of simple models for use in turbulent combustion studies. My t h r e e questions are: 1) Does argon explain the discrepancy in the predicted O~ leakage to the rich side? 2) Is it not the H + O2---, OH + O reaction which breaks down at extinction? 3) Reaction of fuel on the rich side, ~b > 1.5 ap-

pears to be prevented by lack of radicals. Even as • ---> 0 this will be so as the fuel sequesters itself from these radicals by chain terminating reactions at d) = 1.4. There will thus be no approach to equilibrium on rich side as X '-> 0. Do you agree?

Authors" Reply. Replying to the comments in order: (1) Yes, the discrepancy may be very largely explained in this way (see comment by Pitz, and reply). (2) Because of space limitations it has not been

METHANE-AIR C O U N T E R F L O W DIFFUSION FLAMES possible to discuss the chemical structure of the flames, and this will be done elsewhere. Briefly, the effect of c h a n g i n g velocity gradient on the structure of the flame shows that in Lifian's terms (Ref. 28 of paper), the 'equilibrium' branch of the flame is on the air side (where no primary fuel penetrates) and the frozen branch on the fuel side. Departure from equilibrium starts to occur at high temperatures on the air side, where net radical production occurs by the reactions of the H e / C O / O2 system, the H2 and CO having diffused across from the richer region where they are produced. The important step controlling the radical production is the reaction H + 05 --) OH + O. It is this reaction that must break down at extinction. (3) As X tends to zero, or the velocity gradient tends to zero, there will still be frozen flow on the fuel side of the flame. However without detailed figures at low velocity gradients being available, I would hesitate to quote a stoichometric ratio above which the freezing occurs. For the velocity gradient a = 100 s -~ the computations of group B show that the rate of consumption of primary fuel by radical attack is still, when XcH4 = 4 x 10 -z and the mole fraction ratio Xcn4/Xo2 = 10, some 8% of its maximum value.

W. Pitz, Lawrence Livermore National Lab, USA. The calculated value of 05 concentration on the fuel side of the diffusion flame was lower than the measured value. Methane pyrolysis on the fuel side may be affected by the amount of O2 present. Can you speculate concerning the reasons for the difference in 05 concentrations?

Authors" Reply. It now appears that the measured molecular oxygen breakthrough concentrations reported by Tsuji and Yamaoka also include inert gases (argon) to the extent of approximately one percent of the local nitrogen content. If correction is made for this, the agreement between the measured and calculated concentrations on the fuelrich side of the flame is greatly improved.

T. Takeno, University of Tokyo, Japan. It has been known that major structures of the counterflow diffusion flame can be predicted well by the simplified flame surface model. L The details which could not be correlated satisfactorily with the experimental observation were the boundary layer thickness, the maximum flame temperature, and a dent which appears in the fuel side temperature distribution and was attributed to fuel pyrolysis. The introduction of detailed kinetics and transport properties were expected to explain these discrepancies. However, in

1903

your results the improvement is not particularly significant. How do you explain this?

REFERENCE 1. TAKENO, T., Combus. Sci. Techn., 5, 99 (1972).

Authors" Reply. The flame surface or Burke-Schumann model indeed explains the major features of the diffusion flame structure to fair approximation, and the first two of the departures indicated by Professor Takeno have nothing to do with detailed kinetic and transport properties. The discrepancy in the boundary layer thickness is due to the use by Takeno of an incorrect velocity gradient in his calculation. As shown in Fig. 3 of the paper, ff 'allowance is made for some particle lag in following the higher accelerations in the flame, the measured velocity profile is reproduced very well with a velocity gradient of 130 s -l. This velocity gradient is an eigenvalue of the similarity solution of the experimental system; and as noted in the paper, the departure from the nominal gradient of 100 s -~ is probably associated with pressure gradient effects inherent in the specific combustion chamber geometry. The discrepancy between the predicted and measured maximum temperatures is due to the omission of a radiation correction to the reported thermocouple measurements. This explanation is further supported by the agreement of the maximum computed velocity in curve 3 of Fig. 3 with the experimental observation. The third discrepancy, that concerning the dent in the fuel side temperature distribution, would indicate a remaining shortcoming in the detailed reaction mechanism, in particular in relation to pyrolysis reactions. However, although the 'full' reaction mechanism of Miller et al. (group D) predicts the observed concentrations of C2H6, C2H4 and C2H2 reasonably well, it does not predict such a dent. The nature of the possible shortcoming thus remains uncertain. It may also be observed that such steep temperature profiles are extremely difficult to measure precisely by probing flame systems with fine thermoeouples (cf. Dixon-Lewis and Isles, Proc. Roy. Soc. Lond. A 308, 517 (1969)), and the possibility of the dent being due to some flame disturbance by the thermocouple should not therefore be completely discounted.

B. Raghunandan, Indian Institute of Science, India. We have conducted some experiments on fiatplate boundary layer flames in a combustion chamber. As you mentioned in your presentation, the confined nature of the flow is found to affect the flame characteristics. In addition we have observed

1904

LAMINAR FLAMES

that buoyancy is also important, depending on how close the diffusion flame is to the injecting surface. I expect that these two effects would be significant in these types of counterflow diffusion flame experiments you discussed. In particular the velocity profiles across the boundary layer velocity gradient induced by buoyancy and combustion chamber geometry. Please comment?

Authors" Reply. The paper does indeed, in Sec. 4(b), draw a t t e n t i o n to an effect of c o m b u s t i o n chamber geometry in determining the effective velocity gradient of the system, and in Sec. 4(c) attention is drawn to an experimental observation by Tsuji and Yamaoka of a change in extinction limit with cylinder radius. These are both effects associated with combustion chamber geometry.

H. Mukunda, Indian Institute of Science, India. I thought it would be appropriate to make a comparison of inert (Nz) profiles including wall mass fractions for many conditions of test before embarking on comparisons of reactive species so that the effect of kinetic data and quality of experimental data can be distinguished. Could you indicate such a comparison of inert profiles at various values of free steam speeds?

Authors" Reply. I am not quite clear about the purpose of this question. For a constant Vw, the computed nitrogen mole fraction profile when plotted as a function of the similarity variable "q, is changed but little by alteration of the velocity gradient.