Influence of laminar flame speed on turbulent premixed combustion

COMBUSTION

AND FLAME

6 4 : 3 5 3 - 3 6 7 (1986)

353

Influence of Laminar Flame Speed on Turbulent Premixed Combustion GORO MASUYA National Aerospace Laboratory, Kakuda Branch, Ohgawara, Miyagi, Japan

A laminar flamelet model utilized with the Bray-Moss-Libby model to analyze premixed turbulent combustion is revised to include effects of pressure change across flamelets as well as those of chemical production and molecular transports. The revised model introduces the laminar flame speed as an important parameter and exposes a new mechanism of a countergradient diffusion flux of products. In contrast to the original model, contributions from the flamelets to the budget of the velocity fluctuation intensity are not necessarily dissipative but possibly productive. In planar turbulent flames normal to oncoming flows, transverse components of the velocity fluctuation are amplified only by this effect.

1. I N T R O D U C T I O N In order to analyze the structure o f turbulent flames in a premixed reactant flow with high turbulent Reynolds and Damkoehler numbers, Bray, Libby, and coworkers have developed a second-order closure method [1-5] utilizing the B r a y - M o s s (BM) model o f thermochemistry [61, the B r a y - M o s s - L i b b y (BML) model o f aerothermochemistry [1, 6], and a laminar flamelet (LF) model [7]. The BM model along with an assumption o f high Damkoehler number leads to such a mean flame structure as the reactants and the products are separated by flamelets or thin burning sheets. The instantaneous thermochemical state at a certain location is expressed by a scalar variable whose probability density function (pdf) is effectively bimodal. The velocity statistics conditioned by this variable play an essential role in the B M L model. The LF model deals with the contribution from the flamelets to dissipation o f fluctuations and with chemical effects. Their results predict that an interaction between density fluctuation and mean pressure gradient causes countergradient Copyright © 1986 by The Combustion Institute Published by Elsevier Science Publishing Co., Inc. 52 Vanderbdt Avenue, New York, NY 10017

diffusion [1, 2] in flames normal to oncoming reactant flows. In addition, turbulence is intensified in the flames by the same mechanism [2]. Although the predictions are in " e n c o u r a g i n g " agreement [2] with the appropriately interpreted experimental data [8], some o f the closure assumptions must be revised or improved to set the analysis on a firmer basis and to make a better agreement with experimental data. For example, several alternatives of an empirical relation for one of conditional velocity fluctuation intensities, i.e., those within the reactants or the products, are tested in Refs. [1-5]. In order to eliminate such an empiricism, Masuya and Libby [9] add a conservation equation for a third-order correlation, which is a product of two fluctuating velocity components and a scalar fluctuation. Their formulation is adopted in the present study. Another example is the effects of fluctuating pressure gradients. The previous closure methods [1-3, 5, 9] neglect these effects, reasoning that the thickness o f the turbulent premixed flame, which is on the order o f the integral scale o f the o n c o m i n g flow, is insufficient for them to

0010-2180/86/$03.50

354

GORO MASUYA

be significant. This neglect is rather controversial [10]. Strahle and Chandran [11, 12] take a different approach for the bulk pressure gradient terms including both mean and fluctuation. These terms are modeled considering only unconditioned quantities. Although they adopt the same closure assumption as in Ref. [2] for the other terms, their results are quite different from those in Ref. [2]; turbulence is produced in a flame by the bulk pressure gradient terms but the countergradient diffusion is not predicted [12]. In the present work, we revise the LF model [71 following the way taken in Ref. [13] to deal with the effects of fluctuating pressure gradient across the flamelets, which is an essential difference between turbulent premixed flames and constant density flows.

2. DESCRIBING EQUATIONS 2.1. The Thermochemistry Since the propagation Mach number of the turbulent premixed flame is usually much lower than unity, the thermochemical effect due to the pressure change is negligible. Adopting the BM model, we can express all the thermochemical states in term of a scalar reaction progress variable, c(x, t) [6], which may be considered a normalized concentration of the products. The equation of state leads to the following expression:

T/T~= p~/p = 1 + rc,

(1)

where o and T are the density and the temperature respectively, where suffix r denotes the condition within the reactants, and where r is a heat release parameter such that the temperature of the products T o = (I + r) Tr. Taking the average of Eq. (1), we obtain a mean equation of state:

T/T~=p,/~= 1 + r 6 ,

(2)

where the tilde means Favre-average and the overbar indicates the conventional Reynoldsaverage.

The pdf of c can be written, in general, as P(c; x) = ~(x)8(c) + ~(x)6(1 - c)

+y(x)H(c)H(1-c)f(c;

x),

(3)

where c~, /3, and Y are relative intensity of the reactant (c = 0), product (c = 1), and burning (0 < c < 1) modes, respectively, f is a normalized pdf of c in the burning mode, and 6 and H are the Dirac delta and Heaviside functions, respectively. Introducing the fast chemistry or the thin flame assumption, we assume that the value of y is much smaller than unity throughout the flame. Then the pdf P(c; x) becomes bimodal and the intensities of the other modes, ct and 13, as well as all the nondimensional correlations involving only state variables such as T, p, and c become simple algebraic functions of 6 and r [1, 4]. The fast chemistry assumption is valid for a large value of Damkoehler number such as Da = (lo/U~)/(IL/UL) where 1o and u o are the integral length scale and the root mean square of a streamwise velocity fluctuation intensity in the oncoming flow, and where IL and UL are the thickness and the propagation speed of an unstrained laminar flame in the same reactants, respectively. Bimodal pdfs of state variables are observed in various type of flames [8, 14-18]. Thus the fast chemistry assumption is useful to analyze turbulent flames, at least as a first approximation.

2.2. The Aerothermoehemistry The BML model utilizes joint pdf of u, and c to deal with the velocity statistics [1]. As the fast chemistry assumption is introduced, two conditional pdfs within the reactants and the products are of essential importance. The unconditioned statistics containing velocity components are readily rewritten in terms of 6, r, and the conditional velocity statistics. Most such relations used in the present paper are exhibited in Ref. [4], so that we write only a new relation which appears in the additional conservation

355

T U R B U L E N T PREMIXED FLAME equation for the third order correlation: pu,"uTu~'c"/pr = h-J[(1 - 3,5+ 3,se)(u,p -

Utr )

x ( u j p - ujr)(u~p- ukr)

+ {(1 - ,5) u,'ul' o + ,su; uL}(Ukv- U,r) -4- {(1 -- ,5)//f g/~p q'- C U j t U/~r} (//t p -- Utr )

+ {(1 --~)U;~U,'p+,SU~U,'~}(Ujp--Ulr) + u,'ujU/~p- u , ' u j u ~ l + 0(30,

(4)

where h = (1 + r`5)/{,5(1 - `5)} and ( ) ' denotes fluctuation from the conditional mean value.

2.3. The Conservation Equations The following Favre-averaged conservation equations are used to describe a turbulent flame in a premixed reactant flow with low Mach number and high turbulent Reynolds number:

O(~a~)/Oxk =O,

(5)

O(tStik,5 + p U/; c" )/Oxk = ~i,,

(6)

a(,,5~ka, + pu/," u,")lOxk = - Ol~lax,,

(7)

c)(~kpu," c" +pu/[u," c")/Oxk + pu/,' u," O,5/Oxk + pu/," c" 06,/cDxk = u," ( w + (9O~k/OXk) + C" O( -- 6,kp + (rlk)/'OXk,

(8)

C3(tJkpUt" Uj" + p tile' U," Uj" )/OX k + pU/,' U7 0~,/OXk + pU/," Ui" at2~/Oxk = u f o( - 6,kp + O,k)/OXk + U/' O( -- 6~p + %k)/CDXk,

(9)

O(9~pU,"UTC" +pU//U,"Uf' C")/OXk +(Ukpt, - l ~" Uj" + pu/; u," u T) ae/OXk A¢_ ( ~kpl4j"C" ..b pldl~ " LIj " c " )

O£t,/Oxx

+ ( ~ p u , " c" + pu;/ u," c " ) O~j/Ox~

where x, is a Cartesian coordinate, u, is the x,component of velocity, p is the pressure, w is the chemical production rate of c, a U is the viscous stress tensor, ac, is the molecular diffusion flux of c, 6,~ is the Kronecker's delta, and ( )" is a fluctuation from the Favre-average. Terms relating to mean molecular transports are omitted because of high turbulent Reynolds number. Equation (10) is added to the describing equations used in Refs. [1]-[5] to eliminate the need of an empirical relation for conditioned intensities [9]. Although one of such relations used in Ref. [2] is qualitatively supported by a recent experiment [19], it is essentially an ad hoc assumption and difficult to apply to multidimensional flows. Employing Eq. (I0), we are free from the empiricism on the second-order conditioned velocity correlations. In this sense, the lebel of the present closure is the same as that of Spalding's two-fluid model [20], and corresponds to the Reynolds stress closure for a constant density flow.

3. CLOSURE ASSUMPTIONS 3.1. Outline of Closure Among terms in Eqs. (5)-(10), five kind of terms are required to be modeled to close the equations; (1) correlations including third-order velocity fluctuations such as p u " uj"u~c" ; (2) chemistry terms such as u / ' w ; (3) molecular diffusion terms such as u " 3ack/Ox~; (4) pressure gradient terms such as c" Op/Sx,; and (5) viscous stress terms such as c" Oo,~/OXk. The BML model leads the closure problem of the unconditioned correlations including thirdorder velocity fluctuations to that of the conditional third-order correlations. A conventional gradient transport model might be applicable to the latter correlations. However, in the present study, calculation is made only for the normal flames, so that a simpler closure assumption of zero conditional third-order correlations [1, 2, 9] is adopted here, namely,

= Uf'C" O( -- 6,kp + O,k)/OXk t t l U t ldj U R r = l i t

+ U," C" O(--bjkp + ajk)/aXk + U," Uf (W + O(rJOxD,

(10)

O.

(ll)

In the special case o f i

= j = k, Eq. (11) is

p

t

t

lt t U k p -~-

356

GORO MASUYA

satisfied by symmetrical conditional velocity pdfs on their mean values. Such pdfs are observed in experiments [8, 22]. Equation (11) generally means that there appears no crosscorrelation in the normal flame. Since unconditioned cross-correlations are zero in the normal flame, this assumption is plausible. The experimental results of Cheng [23] seems to indirectly support Eqs. (11). The chemistry and molecular diffusion terms have no contributions from the reactant and product modes. Those from the burning mode are estimated by the revised LF model in the next section. The pressure gradient terms are represented only by their mean values in our previous works [1-3, 5, 9]. Effects of the fluctuating pressure gradient are neglected because the thickness of the turbulent flame is insufficient for them to be significant as in incompressible flows [1]. Such an argument may be reasonable in the reactant and the product modes where the density is constant, and we neglect the contributions from these modes. In the burning mode, on the other hand, the instantaneous pressure gradient due to acceleration caused by density change should be strongly correlated to the changes in c and velocity components, and their correlations may not be negligible. It is evaluated by the revised LF model, too. The viscous stress terms have nonzero contributions from the reactant and the product modes but these contributions in the case of the normal flame are presumed to be negligible compared with that from the burning model [1, 7]. Here we adopt this hypothesis and evaluate the contribution from the burning mode using the revised LF model.

(13)

puo d u , / d ~ = d ( - p

(14)

+ o~O/d~l,

where ~7is a coordinate perpendicularly fixed on the flamelet. Only Eqs. (12) and (13) are used in Refs. [7] and [13]. Equation (14) is added to deal with dynamic effects such as viscous stresses and pressure. Comparing terms in Eqs. (8)-(10) and the right hand sides (RHS) of Eqs. (13) and (14), we find that the modeling is greatly facilitated by use of the following coupled variables: qc = w + OOck/OXk, q, = O( - •,kP + (r,k ) / OXk .

(15)

A similar coupled variable is used by Kuznetsov [21] in estimating the interaction effect of fluctuating pressure and velocity divergence on the turbulent kinetic energy. Equations (1) and (12) lead to a simple relation between the gradient of c and these variables: q~= p~UL dc/d~?,

q, = rPrUL2q~, d c / d ~ ,

(16)

where the @i are direction cosines specifying instantaneous direction of the flamelet. Since assumed large values of Da mean that the time scale of turbulence is much larger than that of the laminar flame, the direction and the upstream velocity of each flamelet are almost unchanged during its passage of a certain point. Then the burning mode pdf f ( c , x) becomes independent of x and is determined only by the laminar flame structure [6]: f ( c ) = IL/( dc/drl),

(17)

where IL is defined as

3.2. The L a m i n a r F l a m e l e t M o d e l

The LF model [7] considers the turbulent flame as an ensemble of unstrained planar laminar flames which are governed by the following equations: dCo u O / d , = 0,

p u,~ dc/d~l = w + dac~/d~l,

(12)

IL=

f

l

(d~l/dc) dc.

(18)

0

Namazian pdf of the Eq. (17). written in

et al. [22] observed the burning mode density which is well represented by The mean reaction rate ~ can be terms of the burning mode intensity

T U R B U L E N T PREMIXED FLAME

357

and the laminar flame parameters: (I 9)

I,P = "~(Pr UL//L).

Equations (1), (12), (16), (18), and (19) lead to the following expressions of the contributions from the burning mode, which is indicated by suffix b, to the various terms in Eqs. (8)-(10): 3/c"q,b = I'PT"HLq~,(1 -- 26)/2, "Yu['qcb = ~'{(~a--t/,) + ZULqS,/2}, vtli"qJb = 1'~'7"UL{(tdta --/dt)6y + UtadPf + rULC~,d~j/2 } ,

3'U," Ug"qcb =

I'v[(Uta

-

-'l- 7"/dL{(Uta

ldi)(Uda -- ldd) d- Utt uf a --

Idi)6d -}- (~ja

--

aj)~t

+u,'Oj~+ujO,a}/2 + (rUL)Z6i%/3}, q/U," c " q:, = ~rUe[(1/2 -- C-){(U,~ -- t~,)~j + U,'~j~} + ZUL(a,4~j(1/3 -- ~/2)],

(20)

where suffix a denotes a condition in the burning mode with c --, 0, that is, one at the cold edge of the flamelet. The terms modeled in Eqs. (20) are called flamelet terms. Since the left hand side (LHS) of Eqs. (20) also contributes to the mean pressure gradient terms, we must subtract such contributions as

"y(Op/Ox,)b =

-

~'ruL6,

(21)

from the mean pressure gradient. Here we used a relation for a gradient of c appearing later. Although Eqs. (20) are not yet closed because of unknown correlations among the cold edge velocity and the direction cosines, it seems better to review here their features. The most important point of the revised model compared with the original one [7] is the appearance of the laminar flame speed UL, which is justified by the following argument. I f an activation energy of the flamelet is sufficiently high and both the Prandtl and Lewis numbers are near unity, the turbulent Reynolds number Re is on the order of (UJo)/(ULIL). The assump-

tion that the Reynolds and Damkoehler numbers are very large is then expressed by the following inequalities: IL/lo "~. Uo/UL ~. lo/IL.

(22)

We find that the large value of Io/IL is an essential requirement while the value of Uo/UL is not necessarily large. Therefore it is not surprising to see that only the laminar flame speed appears in the model as a parameter and that the laminar flame thickness does not. In Eqs. (20), UL is always multiplied by the heat release parameter 7" to represent the acceleration in the flamelet. The difference of two conditional mean velocities within the reactant and the product modes, which is used instead of 7UL in the original LF model [7], depends upon the history of each mode and is not suitable to represent the local acceleration by the flamelets. We can see that the third of Eqs. (20), which is for the Reynolds stress transport equation, is proportional to rUL. If the value of rUL is much smaller than the velocity scale of the flow such as Uo, then the contribution from the burning mode to the Reynolds stress is negligible. However, the others in Eqs. (20) have terms not proportional to rUL, and they do not vanish even if T//L//./° is much smaller than unity. The second feature of Eqs. (20) is the existence of the flamelet direction effects. They were suppressed in the dissipations of scalar quantities, e.g., the fluctuation of c [13], and, with some additional assumptions, the turbulent kinetic energy [21] because of such an identical relation for direction cosines as

~kk

=

1.

(23)

The flamelet direction affects vector or tensor quantities. Among the statistics of the ~b,, the mean value 4]; is exactly related to the mean spatial gradient of c by the LF model. Since the value of c varies only within the burning mode, we have Oe/Ox, = ~(OC/3 X,)b = ~ , ( d c / d ~ ) b :" "Y~,/IL = I¢~,/ PrUL.

(24)

358

GORO MASUYA

We need experimental data on the direction and the movement of flamelets to model higher order correlations, but there are only few available experiments [24, 25]. The last feature of Eqs. (20) is the inclusion of the flamelet cold edge velocity U,a in the model. The statistics of U,a are not necessarily equal to those of u,r except at the hot boundary of the mean reaction zone. We need experimental data on U,a, too, but they would be more difficult to measure than the flamelet directions are. In order to demonstrate the influence of the burning mode on the budget of the turbulence intensities, we inspect a special case of a normal planar flame. The transverse component is of interest. The transverse coordinate and the velocity component in this direction are denoted by y and v, respectively. Any mean value of v is zero in this case. If we assume that the fluctuations of ~by and Va are not correlated, then v-V~yy= ~ {v--~+ (rUL)2~---7}.

In our previous works on the normal flames [2, 9], it is not necessary to solve numerically the equations for the transverse velocity fluctuation, because there is no net source or sink term in those equations under the previous closure assumptions. As mentioned in the preceding section, however, the revised LF model brings in a net source term. Thus those equations should be numerically integrated, too. We take the x-coordinate normal to the mean reaction zone and the y- and z-coordinates in transverse directions. The structure of the mean reaction zone is uniform and isotropic in the transverse plane, so that either the y- or the Zcomponent can represent the transverse components. The x- and y-components of the velocity are denoted by u and v, respectively. The constraints on the flamelet direction cosines are as follows:

~v~

0,

(25)

This relation clearly shows that there is a net change in the fluctuation intensity of the transverse velocity component. Since the RHS of above equation is nonnegative, the flamelet terms do not dissipate turbulence but produce it through acceleration across flamelets in different directions. The same mechanism of turbulence production is predicted by Clavin and Williams [26] for wrinkled laminar flames. In our previous works [1-3, 5, 9], only the viscous stress is taken into account at the flamelets, and they dissipate velocity fluctuations. Thus the pressure gradient at the flamelets included in the present model causes this difference. 3 . 3 . T h e Final C l o s u r e for the N o r m a l Flames

Since our knowledge on the directions of the flamelets and the cold edge velocity in it is quite limited, our work on the final closure is restricted to the simplest case of the normal planar flames for which several relations for unknown correlations are constrained by its geometrical simplicity. We adopt Eq. (11) for conditioned third-order correlations.

0y 2 = (1 - 0x2)/2.

(26)

We see that only ~x2 should be modeled. There are two extreme cases. The first case is that ~x2 takes its minimum value ~x2 and the pdf of ~x is a single delta function at ~x. Such a pdf, if possible, may be observed for very weak upstream turbulence. The other extreme is that ~x2 takes its maximum value of unity and the pdf of 0x is represented by two delta functions at 1 and - 1 , which means that the flamelets are always normal to the mean flow and propagate upstream or downstream. Such a flame structure is not likely to be realized even if the initial turbulence is very high. A measurement of instantaneous flame angle was conducted by Suzuki and Hirano [25] using an electrostatic probe with three sensors. Their results are essentially two-dimensional because of the number of sensors, and have inherent distortion from their three-dimensional counterparts. For example, even if the three-dimensional pdf of ~x is a single delta function at a certain value, the two-dimensional pdf of ~x becomes a continuous function between that value and unity. Never-

T U R B U L E N T PREMIXED FLAME

359

theless their result shows important qualitative features of the pdf of the flamelet direction and suggests that it is a continuous function with a single peek near its mean value. I f we assume that the pdf of ~bxbecomes uniform for ~¢er__yhigh upstream turbulence, then ~x = 0 and ~bx2 = 1/ 2. Considering these results, we tentatively adopt the following model for ~bx2:

~x"--2--~x2=(1--~x 2)

exp(--kjUo/UL)/2.

(27)

This model should be revised in the light of future experimental data. Since there are no available data on the cold edge velocity components Ua and va, it is better to assume a simple model for them. The simplest form of linear change between the boundary values, which is specified later, does not work well, so we assume a slightly complicated form: /~a --/~ao : (/~r -Oa - - Vao = ( O r -

uLZ)g(c;

k3), (28)

where g is a function of c which monotonically changes from zero to unity, and where k2 and k3 are constants whose values are determined to match the solutions with the boundary conditions by iteration. However, since the choice of the value of the constant k3 is found to have only little influence on the solutions, it is equated to k2 to reduce the complexity of iterations. In Eq. (28), a local value of the conditional mean velocity within the reactants is used instead of its hot boundary to avoid further iteration. Three functions are tested as g: g= 1-(l-Ok2, g = 1 - exp{ -

k26/(1

-

We assume that the turbulence of the oncoming flow is isotropic for the sake of further simplicity. The solutions of Eqs. (5)-(10) must satisfy the following conditions at the both boundaries of the mean reaction zone. At the cold boundary o f x ~ --0o, 6=0, pti" C" = p U " 20" : D O " 2C'; = 0 , p U " 2 = p v " 2 = PrUo 2,

(30a)

and at the hot boundary of x --, oo,

p u " c " = p u " 2c" = p v " 2c" = O.

k2),

u~' 2 _ U,o2 = (U/2 _ U/o2)g(c; k3), Ua 2 __ U~o2 _~_(O r 2 __

3.4. The Boundary Conditions

6=1,

tJao)g(c; k2),

17ao)g(c;

no correlation between the cold edge velocity and the direction cosine.

6)},

g = m i n i m u m of k2c~ and 1 - ( 1 - ~ ' ) / k 2 .

(29)

Since we use the unstrained planar laminar flame to approximate the flamelets, the direction of the flamelet does not affect the properties within the flamelet. Thus we assume that there is

(30b)

We see that there are too many conditions. They are satisfied by specifying suitable values for parameters, namely eigenvalues. In addition to the conditions above, we can impose those for the direction cosines and the cold edge velocity of the flamelets. Many flame visualization experiments have shown that a cusp structure of flame sheets does not appear near the cold boundary while it is often observed near the hot boundary. Therefore the mean direction angle of the flamelets should be normal to the x-coordinate at the cold boundary. No simple boundary condition for the mean flamelet direction exists at the hot boundary because of the cusp structure observed. At the cold boundary, the mean propagation speed of the flamelets in the x-direction must be equal to the streamwise component of the mean velocity within the reactants entering the flamelets which is the cold edge velocity, because the flamelets do not go further upstream. As the mean direction of the flamelets is normal to the x-direction there and the flamelet is approximated by an unstrained laminar flame, tZa = UL. At the hot boundary, all the reactants enter the flamelets and therefore aa = t~r as mentioned

360

GORO MASUYA

previously. The mean value and the fluctuation intensity of the transverse component of the cold edge velocity at the both boundary are equated with those in the reactants, that is, 0 a = O r and v~'z = v/2 throughout the flame. Since there appears the product of {6(1 c')} -~ and p u " c " , p u " l c ", or Ov"2c " in the equations, the dependent variables at both boundaries must be expanded. Expansion is made to the second order o f d o r (1 - c'), but the noninteger order used in Refs. [1]-[3] and [9] is not introduced. Coefficients of expansion are determined by substituting the expanded variables in the equations and collecting the same order terms. We find that u~'2 and u p,2 at the cold boundary and u~'2 at the hot boundary are zero. We also find that the values of p u " c " / ( 1 - c') and p v " 2 c " / ( 1 - c') at the hot boundary can be arbitrarily specified. This may be the reason why we need only one eigenvalue, k2, to satisfy all the boundary conditions.

Q.~ = - 2"y{(u " qx)b -- U" (Op/OX)b} / {prl2o3(dc/dx)}, Qvy = - 2v(o" qy)b/{pdio3(d(/dx)}, Qxxc = - -y{(2u"c"qx+ u" 2qc)b - 2u "c" (Op/ax)b} /{prao3(d(/dx)},

Qyyc = -'Y(v" 2qc)b/{o, aJ(de/dx)}.

Using Eqs. (5)-(7) to eliminate t/, ~, and dlM dE, we have the final form of the equations: {(1 + r d ) F ~ } ' + K L + r F x + I ~ (34)

= r(r + I / x ) / h - Qxc,

{(1 + r d ) I x x } ' +Lxxx+2rIxx (35)

= 2 r F j r + l/x) - Q ~ ,

{(1 + rd)K~x}' + N ~ +

2r{(1 + r 6 ) F x + gx~}

+ (1 + r ( ) I ~ + Lxxx

3.5. The Numerical Calculations

= 2(1 + r - r 6 ) F x ( r + I]x) - Q~,., In order to facilitate the numerical treatment, the independent variable x is replaced by c. We introduce the following nondimensional dependent variables:

Fx= Pu" c" /prl2o,

/p~ t~o2,

Kyy --~po" 2C"/Prl~o2,

(31) and auxiliary variables: L x ~ = p u " 3/,Or/.~o3,

Lxyy=pUnU " 2/pr/.~o3 '

N~rx = p u " 3c"/prt2o 3,

Nxyy = p u "v " 2c tt/Or ~o 3, and flamelet terms (FL): Qxc = - 3' { ( c" qx + u "q~)u - c" (Op/0X)b } / { p , ao2(d?/dx)},

{(1 + r6)Iey } ' + L~ey= - Qyy,

(32)

(36) (37)

{(1 + rd)Kyy}'+N/~yy + (1 + rd)Iyv + Lxyy = - Qyyc,

lxx= pU" 2/pr1~o2,

Iyy = pV " 2/pr/~o 2, K~ = ~

(33)

(38)

where ( ) ' denotes the derivative of e. Each of Eqs. (34)-(38) may be identified by the variable involved in the first term of the LHS; for example, Eq. (34) is the Fx-equation. In the LHS of Eqs. (34)-(38), the first and the second term are called as a convection term (CV) and a diffusion term (DF), respectively, and the term with a factor r is a dilatation term (DL), while the others are production terms (PR). In the RHS of those equations, the term including a factor (r + I~x) is a negative of a mean pressure gradient term (MPG). As in our previous works [1-3, 5], the present formulation does not predict the turbulent flame speed Uo but can specify Io which is an upstream value of I~x and is related to that speed. When the solution is compared with a specific experiment, measured values of Io and Uo/UL in that

T U R B U L E N T PREMIXED FLAME

361

experiment are used in the calculation. In calculations to show the effects of laminar flame speed and heat release, we use the turbulent flame speed data compiled by Abdel-Gayed and Bradley [27]. Among the data, those with the largest Reynolds number are adopted to satisfy the condition of large Reynolds and Damkoehler numbers. They are represented by the following linear relation: /to : UL "1- 1 . 9 6 7

u'.

(39)

Among Eqs. (34)-(38), the Fx-, l~x-, and Kxxequations do not contain Iyy or Kyy. Therefore we integrate only those three equations during the iteration to determine the value of the constant k2. The iteration starts by assuming a value for the constant kz. In every iteration, the values of the dependent variables at a small value of ,5 (typically ,5 = 0.0001) are calculated utilizing the expanded form of those variables. The equations are numerically integrated with the Runge-Kutta-Gill method up to the hot boundary (typically 6 = 0.9999). The expanded form of the dependent variables is used to calculate their values there for comparison with those obtained by numerical integration. Until both values agree with a certain error, the value of the constant k2 is corrected and the iteration continues. Since the values of variables compared at the hot boundary are very small, the numerical procedure is carried out with double precision. The solutions are usually obtained after several iterations except for the case of a large value of the ratio UJUL and a near zero value of r, for which we cannot have a solution.

for r = 5 and U'/UL = 1.5. The transverse component Ivy is most significantly affected because main part of the net source term FL in the Iyy-equation is proportional to 1 - ~bx2. The model that ~bx2 = 1 may not be realistic, as mentioned previously, and thus the ranges of realistic solutions are somewhat narrower than that shown in Fig. 1. The streamwise component Ixx is not so much affected as lyy is because of the existence of other source terms in the Ixxequation. The turbulent diffusion flux Fx, which is not shown here, is least affected, because its conservation equation does not explicitly include ~bx2. Note that the turbulent kinetic energy is also affected by the choice of the model for ~b,:2, because there are other source terms in the Ix~-equation. Equation (26) with kl = 1 is used to obtain the predictions shown below. The model for the cold edge velocity in the flamelet Ua affects the solution far less significantly than that for Sx2 does. Since the constant k2 appearing in the model is not arbitrarily specified but determined as a part of the solution, the distributions of/~a and Ua 2 with various models are nearly the same. The results shown hereafter is obtained utilizing the second model in Eqs. (28).

4.2. Comparison with the Experiments The results of the present closure are shown and compared with the experimental data of Moss

05

r

/////

4. R E S U L T S A N D D I S C U S S I O N 03

4.1. The Model Sensitivity

\,

0.2

In order to check the sensitivity of the solutions to the choice of the models for the flamelet direction cosines, the calculations are made with three models; one is expressed by Eq. (26) with kt = 1 and the others are two extreme cases of Sx2 = I and ~bx2 --- q~x2. The latter bound the ranges of possible solutions. The Favre-averaged velocity fluctuations are shown in Fig. 1

0 --

0'2

0[4

0'6

018

1.0

F i g . 1. S e n s i t i v i t y o f the s o l u t i o n s to the m o d e l f o r f l a m e l e t d i r e c t i o n : z = 5, U J U L = 1 , . _ _ _ _ 5 : - - , Ox 2 = c~,.2 + (1 ~.2) e x p ( - u o / u L ) / 2 ; ----, ~)x 2 = 1; - ' - - Ox -'-5 = 6 2 .

362

GORO MASUYA

[8] and our previous predictions [2, 9] in Figs. 2a-2c. The values of r and Uo/UL are estimated to be 6.5 and 1.2, respectively [8]. In Refs. [2] and [9], data of Moss are interpreted for the velocity statistics normal to the mean reaction zone under the assumption that the transverse 0,20

w

0 IE

i

,

I

,~PRESENT

0-0~

o.o,

t o • MOSSN O'2

t

o4

L

06

J

0.8

1.0

Fig. 2a. Comparison of predictions and e×perlment: r = 6.5, /o = 0.]6, u'/uL = 1.2. The nondimensional mean flux of product: Q, calculated from the conditioned mean velocities; O , calculated from the measured correlation

pU" c". 16

o:2 o:~

0:6 o2a

~.o

Fig. 2b. The normalized unconditional intensity in direction measured by Moss [8].

3

05

0'5

Fig. 2c. The normalized conditional intensity m direction measured by Moss 18].

velocity components has zero mean and frozen fluctuation intensity. The present closure adopts the former assumption but not the latter. Therefore the velocity fluctuation intensity in the direction measured by Moss [8], which is denoted by U, is utilized for comparison. Figure 2a shows the nondimensional Favreaveraged flux of c normal to the mean reaction zone, i.e., Fx. We see that there is little difference among the predictions and that they are in good agreement with the experimental data. The only one qualitative difference between the results of the present model and the previous ones is the sign of the flux near both boundaries. The former predicts positive countergradient diffusion flux while the latter's prediction is negative. The previous closures [2, 9] require negative flux regions at both boundaries, but the present one with specified values of the parameters does not. This difference leads to a new physical explanation of the mechanism of the countergradient diffusion, which is discussed later. The nondimensional velocity fluctuation intensity pU"Z/Ora2° is shown in Fig. 2b. Though the predictions qualitatively agree with the experiment, the result of the present model is higher than the experiment, while those of the previous models are lower. Neglect of the contributions from the reactant and the product modes to the dissipation term might be one of the reasons for overestimation by the present model. Instead of pU"2c"/pra2o, we show the conditional intensities U~'2 and Up 2 in Fig. 2c. The flux oU"2c" can be easily calculated from the latter and the solutions in Figs. 2a and 2b using a relation in Ref. [4]. The prediction of the present model for Up 2 is in better agreement with the experiment than those of the previous models are, but that for U~'2 is less favorable. The present model predicts a slight decrease in U '2 through the mean reaction zone. This is in contrast to the data of Moss and the other predictions. However, there is an experiment [23] which supports the present result for U~'z. We need further experimental data to examine the validity of the predictions.

T U R B U L E N T PREMIXED FLAME

363

4.3. The Effects of Heat Release The heat release parameter ~ has played the most important role in our previous models [1-5, 9]. In the present closure, it is also involved in FL. Its influence on the solutions is exhibited in Figs. 3a and 3b for / , / o / / / L = 1.5. The former figure shows the streamwise flux F~ for various value of 7. Similar results are obtained in Refs. [1] and [2]. When the heat release is very low, the flux is entirely negative and thus in accord with a gradient diffusion hypothesis. The countergradient diffusion of c is observed near the hot boundary for a r of unity and throughout the mean reaction zone for r more than three. Figure 3b exhibits the effect of • on the amplification factors of the velocity fluctuation intensities across the flame u'/Uo and v'/Uo, where the subscript oo denotes the condition downstream of the mean reaction zone. In the limiting case of zero heat release, which is not of practical importance but is mathematically

0

0.2

0-4

015

O.B

Fig 3a. The influence of heat release: n o n d i m e n s i o n a l mean flux o f product.

8

,

,

,

6

1.0

U,~/UL =

1.5. The

,

uMu'

2 i

i

T

Fig. 3b. The ratlo of upstream and downstream intensitles.

interesting, the velocity fluctuation intensities upstream and downstream of the mean reaction zone are exactly the same. This indicates that the flame affects turbulence through the density change caused by heat release. As ~- increases the streamwise intensity first decreases and then increases. Such behavior is closely related to the sign of the diffusion flux of c. The transverse intensity monotonically increases because the source term is proportional to ~.2.

4.4. The Effect of the Laminar Flame Speed A parameter representing the thermochemical effects of the reactants in the present model is the nondimensional acceleration across the laminar flame rUL/ao. An alternative parameter is the ratio of the upstream turbulence and the laminar flame speed Uo/ULwhich may be more useful to compare with experiments. Since such parameters as those above have not been introduced in the previous versions of the theory [ 15, 9], we examine their effects in detail. Figures 4a-4e show the solutions for various value of U'/UL while the value of r is kept equal to five. The nondimensional flux Fx shown in Fig. 4a is entirely positive and in a countergradient direction for very low upstream turbulence. As U'/UL increases, Fx decreases. There appears a region of negative flux near the cold boundary for Uo/UL = 10. It is possible to give a simple physical explanation of the mechanism of the countergradient diffusion for the low initial turbulence. In this case, the mean velocity within the reactant at the cold boundary tJro, which is identical with t~o, is on the order of UL. The mean velocity at the cold edge of the flamelets t3a is always assumed equal to UL at the cold boundary. Thus the conditional mean velocity within the products t~po becomes (1 + r)uL and is greater than t~ro, that is, Fx is positive. For higher upstream turbulence, the value of Uro becomes higher than (1 + 7")ULbut apo is unchanged, so that the flux becomes negative. Apart from the cold boundary, ur has higher value than uro and the mean direction of the flamelets is not along with x-axis. The acceleration across the flamelets no longer

364

GORO MASUYA

!

'~I u;luL=O "1

0~

0-2

0.4

0.6

0,8

1.0

Fig. 4a. The influence of the ratio of the laminar flame speed and the upstream fluctuation velocity: z = 5. The nondimensional mean flux of product.

Fig. 4d. The nondimensional streamwise mean flux of product flux.

012 01C

12

_=

u~,/ut=02

1.0

0.0~

0.8

.,~ 0-0(

0.6

O.Ot 0-02

0

0-2

0./,

0.6

0.8

1.0

Fig. 4b. The nondimensional streamwise intensity.

05. 04

03 ~-, 02

"~ 05

:,

o

0

02

04

0-6

0.8

tO

Fig. 4c. The nondimensional transverse intensity.

0

02

04

i

0-6

08

1.0

Fig. 4e. The nondimensional transverse mean flux of product flux.

makes an important direct contribution to the countergradient diffusion there, but the transverse components of fluctuating velocity at the hot edge of the flamelets are converted into the streamwise component within the products and accelerate ap. Note that the both mechanisms are related to the mean pressure gradient. Such mechanisms of the countergradient diffusion as those above are different from that explained by Libby and Bray [1], namely, packets of the lighter products in the matrix of the heavier reactants are preferentially accelerated and packets of the heavier reactants in the matrix of the lighter products are preferentially left behind by the mean pressure gradient. The present

T U R B U L E N T PREMIXED FLAME mechanism may be applicable to the lower upstream turbulence cases for which the flamelet is a continuous surface, and that of Libby and Bray [1] to the higher turbulence for which burning occurs at the surface of the discrete packets. The contribution from each term, which is denoted by the abbreviation defined previously, in the Fx-equation is shown in Fig. 5a for typical values of UJUL. In the figure, the negative portion of each term except CV acts as a source. The balance of terms in the case that uo/u L = 10 is essentially similar to that of our previous result [2] where the flamelet term (FL) is divided into molecular dissipation and chemistry terms. The result for UJUL = 0.1 is quite different from those for higher values and thus from that of previous model [2]. We see that the large negative value of FL in the lower upstream turbulence case modifies the distributions of the other terms and produces a significant positive flux of c near the cold boundary. In the latter cases, the distributions of all the terms are almost the same excep: near the cold boundary. The main source of the'. countergradient diffusion flux of c is the mea:, pressure gradient term (MPG), as explained in Ref. [1]. However, when Uo/UL is unity there still remains a region of negative FL. It disappears and then a region in which the flux of c is negative appears near

365 the cold boundary as the value of UJUL increases to ten. The effects of Uo/UL on the nondimensional streamwise velocity fluctuation lxx is shown in Fig. 4b. Note that its cold boundary value Io varies following Eq. (39) as the value of UJUL changes. Although the apparent effects are similar to those on Fx, the contribution of each term in the equation for Ixx is not similar to that in the F~-equation, as can be seen in Fig. 5b. Again, for the highest value of Uo/UL, the balance of terms except FL is similar to the one shown in Ref. [2] where FL is a sink throughout the reaction zone and is not negligible. In the present results, the magnitude of MPG plays an important role. In the case that Uo/UL = 0.1, though FL has a very large value, MPG also becomes very large because of the large positive value of F~. In the cases that Uo/UL is unity and ten, the direct contribution of FL is negligible. The effect of Uo/UL indirectly appears through

2(}

U:/UL= 01

]

lUL=I

1 0

05

U;/uI. = 10

o

-10 -20

U;IUL =0-1

U'/U L =1

[~u;/ut = 10

0'5

0-5

Fig. 5b. The l~-equation. 0

4

-2 MPG

• /Cv,IO0

I

-4 :k 0

0'5

I 1 0

C~5

L o

0'.5

Fig. 5a. The contributions to the conservation equations for various value of UTUL: r = 5; CV, the convection terms; DF, the turbulent diffusion term; DL, the ddatation terms; PR, the production terms; MPG, the mean pressure gradient terms; FL, the flamelet terms The F~-equation.

0

U~/UL=I

l

U~IUL=IO

05

0

05

Fig. 5c. The Ivv-equation.

366

GORO MASUYA

the increase of Fx which results in larger negative value of MPG for a Ho/UL of unity. In contrast to/xx, the transverse intensity/~y is directly affected by FL as shown in Figs. 4c and 5c. For low upstream turbulence, Iyy increases across the mean reaction zone. As the upstream turbulence intensifies,/yy approaches its asymptotic curve of Io/(1 + rc-), which is the frozen transverse velocity fluctuation intensity assumed in our previus models [1,2, 9]. Figure 5c clearly shows that all the terms in the /yyequation become negligibly small as Uo/UL increases. The distributions of two other dependent variables Kxx and Kyy are exhibited in Figs. 4d and 4e, respectively. Since they are the turbulent diffusion terms in the F~- and Fy-equations, respectively, we see that the gradient diffusion hypothesis is not applicable to them too. Although the shapes of their distributions are not similar to those of I~x and Iyy, the efferts of Uo/UL appear in qualitatively similar ways to those for the turbulence intensities of the same component. Finally, in order to relate the original LF model and the revised one, we discuss the limiting case of very high upstream turbulence. The value of k2 increases as the upstream turbulence intensifies, so that the statistics of u~ is essentially same as those of ur in this limit. The terms including rUL in Eqs. (20) become negligible. Then we find that the original LF model for FL [7, 9] recovered in this case. However, the values of empirical constants appearing in the original model are no longer arbitrary but specified by the revised model. Those values chosen in Refs. [1]-[5] and [9] are not far from the specified values.

from the flamelet to the velocity fluctuation budget is not necessarily dissipative. This is in contrast with the dissipative nature of the original model [7]. Numerical results are obtained only for planar flames normal to the uniform oncoming reactant flow. The predictions are in reasonable agreement with the experiment data of Moss [8]. As in the predictions utilizing the original model [2, 9], the countergradient diffusion flux and the streamwise turbulence production are predicted for practical values of the heat release parameter. In addition, transverse turbulence, which is predicted to be frozen by the previous model [2, 9], is amplified solely by the contribution from the flamelets. When the turbulence upstream of the flame is small compared with the laminar flame speed, the countergradient diffusion and the turbulence production are significantly enhanced by the pressure gradient due to the acceleration across the flamelet. The revised laminar flamelet model introduces several empirical relations on the direction of the flamelet and the velc,'ity within it. We need experimental results f~zr them to improve these relations. In many practical cases, turbulence is so intensive that the adiabatic unstrained planar laminar flame no longer sufficiently approximates the flamelet. Even in such cases, however, the present model would be useful as a first approximation and a basis for installing the various effects [28-30].

5. C O N C L U D I N G R E M A R K S

REFERENCES

The laminar flamelet model [7] for turbulent premixed flames is revised to include the effects of pressure gradient at the flamelets as well as those of molecular transports and chemistry there. The laminar flame speed is introduced as the second parameter of the theory, while the first one is the heat release parameter. Inspection of the revised model exposes that the contribution

The author gratefully acknowledges helpful and stimulating comments from Professors P. A. Libby and K. N. C. Bray.

1. 2. 3. 4. 5.

Libby, P. A., and Bray, K. N. C., A I A A J. 19:205213 (1981). Bray, K. N. C., Libby, P. A., Masuya, G., and Moss, J. B., Combust. Sei. Tech. 25:127-140 (1981). Masuya, G., and Llbby, P. A., A I A A J. 19:15901599 (1981). Bray, K. N. C., Libby, P. A., and Moss, J. B., Combust. Flame 61:87-102 (1985). Libby, P. A., Prog. Combust. Energy Sci. 11'83-96 (1985).

TURBULENT 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18.

19.

PREMIXED

FLAME

Bray, K. N. C., and Moss, J. B., Acta Astronautica 4:291-319 (1977). Libby, P. A., and Bray, K. N. C., Combust. Flame 39:33-41 (1981). Moss, J. B , Combust. Sci. Tech. 22:119-129 (1980). Masuya, G., and Libby, P. A., National Aerospace Laboratory Technical Report NAL TR-802T, 1984. Jones, W. P., and Whitelaw, J. H., Combust. Flame 48:1-26 (1982). Strahle, W. C., and Chandran, S. B. S., A I A A J. 20:129-135 ('1982). Chandran, S. B. S., and Strahle, W. C., Combust. Flame 51:313-324 (1983). Masuya, G., Combust. Flame 56:123-124 (1984). Yoshlda, A., and Guenther, R., Combust. Flame 38:249-258 (1980). Bill, R. G., Jr., Namer, I., Talbot, L., and Robben, F., Combust. Flame 44:277-285 (1982). Tanaka, H., and Yanagi, T., Cornbust. Flame 51:183-191 (1983). Moss, J. B., AIAA Paper 84-1198, 1984. Shephard, I. G., Moss, J. B., and Bray, K. N. C., Nineteenth Syrup. (lnt.) Comb., Combustion Institute, Pittsburgh, 1983, pp. 423--431. Cheng, R. K., Talbot, L., and Robben, F., to appear in Twentieth Symp. (Int.) Comb., Combustion Institute, Pittsburgh, 1985.

367 20. 21. 22.

23. 24

25.

26. 27. 28. 29 30.

Spalding, D. B., AIAA Paper 84-0476, 1984. Kuznetsov, V R., Fluid Dynamics 14:328-334 (1979). Namazlan, M., Talbot, L., and Robben, F., to appear in Twentieth Syrup. (Int.) Comb., Combustion Institute, Pittsburgh, 1985. Cheng, R. K., Combust. Sci. Tech. 41 109-142 (1984). Suzuki, T., Hirano. T., and Tsujl, H., Seventeenth Syrup. (Int.) Comb., Combustion Institute, Pittsburgh, 1979, pp. 289-297. Suzuki, T., and Hirano. T., to appear m Twentteth Syrup. (Int.) Comb., Combustion Institute, Pittsburgh, 1985. Clavm, P., and Wilhams, F. A., J. Flutd Mech. 116:251-282 (1982). Abdel-Gayed, R. G., and Bradley, D., Proc. Roy. Soc. London A301.1-25 (1981) Libby, P. A., and Wilhams, F. A., Combust. Flame 44:287-303 (1982). Libby, P. A., and Williams, F. A., Combust. Sct. Tech. 31:1-42 (1983). Libby, P. A., Linan, A., and Wilhams, F. A., Combust. Sci. Tech. 34257-293 (1983).

Received 13 May 1985; revised 3 February 1986