Chemical closure model for fractal flamelets

COMBUSTION

AND

FLAME

77:241-259

(1989)

241

Chemical Closure Model for Fractal Flamelets F. C. G O U L D I N Cornell University, Ithaca, N Y 14850

K. N. C. B R A Y Cambridge University, Cambridge, England

J.-Y. CHEN Sandia National Laboratories, Livermore, CA 94550

A chemical closure model for premixed turbulent flames is proposed and tested by analysis and numerical computation for flames with vanishingly small density change. The model is based on the assumption that reaction zones can be modeled as thin sheets--flamelets--and that the geometry of these sheets can be represented by fractal surfaces. The model expression for mean fuel consumption rate is ( ,~) = CRpo( a YfUL)f( ly/ n )O- 2lF - ' ( C)(1 -- ( c) ),

with f given by f = [ l - ( I - A t - I / 4 R 1 - 3 / 4 ) exp(--Atl/4Rl-l/4U'/(UL)f)], and I/~=Atlt4Ri3/4 ,

where D is the fractal dimension of the flamelet surface and is the new parameter introduced by the fractal geometry assumption. This model is tested in simplified analyses of normal and oblique flames with good results. The oblique flame analysis provides new insight into the definition of the turbulent burning velocity. Numerical computations are performed with a conditioned second-order closure scheme, and the chemical closure model performance is found to be good. Computed results with a gradient transport model for species diffusion show that turbulent fluxes are significantly under predicted in comparison with the second-order closure results.

NOMENCLATURE

C~

A At

CR

A~, Ao Ca

surface area e m p i r i c a l constant in the expression f o r / / ~ ; see T a b l e 1 surface areas m e a s u r e d at the inner and o u t e r c u t o f f scales constant in t u r b u l e n c e transport m o d e l ; see T a b l e 1

Copyright © 1989 by The Combustion Institute Published by Elsevier Science Publishing Co., Inc. 655 Avenue of the Americas, New York, NY 10010

c, c ' D Dt

n o r m a l i z a t i o n constant in e x p r e s s i o n for Pc constant in c h e m i c a l closure m o d e l ; see T a b l e 1 the reaction p r o g r e s s v a r i a b l e and p r o g r e s s v a r i a b l e fluctuation the fractal d i m e n s i o n turbulent diffusivity

0010-2180/89/$03.50

F. C. GOULDIN ET AL.

242

f &,L,& h, hx, by, g, r k L 1

IF t~ M Min, Mout

P

P gl

U,u U' MLj U0 Ut

Ayf

A0

function accounting for flamelet stretch, see Ref. 1 integrals used in the control volume analysis functions used in similarity analysis turbulence kinetic energy length scale used in discussion of fractal surfaces turbulence integral scale the turbulent flame brush thickness length scale introduced in the discussion of the closure model turbulence generating grid scale mean product mass flow into and out of control volume mean rate of product formation in control volume analysis the probability of finding a flamelet in a small volume, lc3 the ratio of product mass formed per unit mass of fuel reacted the turbulent Reynolds number velocity and velocity fluctuation the absolute turbulent intensity strained and unstrained laminar burning velocity the turbulent burning velocity change in fuel mass fraction across flamelet flame brush half angle

Greek Symbols

Ei, ~o Ed

O, OR, 0p p, p ' Po 00f

(), (/o~

control volume width measurement length scale inner and outer c u t o f length scales mean dissipation rate of turbulence energy flame brush oriented coordinates the Kolmogorov microscale flame brush angles; see Fig. 2 density and density fluctuation reactant gas density fuel formation rate in mass units mean and conditioned (in products) mean quantity

INTRODUCTION In Ref. 1 a model for turbulent burning velocity is proposed based on the assumption that chemical reaction and heat release rates can be modeled as occurring in thin sheets--flamelets--and that the geometry of these sheets can be represented as fractal. It is argued that flamelet surfaces are extremely rough, with multiple scales of wrinkling, and that mathematical relationships true for fractal surfaces can be used to advantage in characterizing flamelet surface area. Fractal surfaces are rough and a significant feature of this roughness is the regular variation of measured surface area with measurement scale [13]. For an isotropic fractal surface in a volume L 3, A / L 3 ~ e 2 - O / L 3,

(1)

where A is the measured area and ~ is the measurement scale. D is the fractal dimension of the surface and is bounded such that 2 < D < 3. Fractal surfaces are extremely wrinkled, and according to Mandelbrot [3] they appear to be space filling. The larger D the more space filling a fractal surface will appear. For flamelets, as for constant-property surfaces in turbulent flow, there is a limit to the range of length scales over which fractal behavior is observed. The variation of A with ~ will not follow Eq. 1 for ~ outside the fractal range, and for surfaces with a limited range of fractal behavior, Gouldin [1] hypothesizes a simple behavior for A versus ~, which is depicted in Fig. 1. At this point two new parameters enter the fractal description, the inner and outer cutoffs--el and eo [3]. Flamelet surface area is finite and thus A must approach a limit, Ai, as e --~ 0. For large e and anchored flames such as a bunsen flame, we expect that on large measurement scales the surface will appear smooth and A ~ Ao as e ~ 0o. This limiting behavior at both large and small e is represented by horizontal line segments, with ei and eo defined such that A o / L a _ % 2 - D / L 3 and A i / L 3 ei 2 - D / L 3. Mean total flamelet surface area is estimated in Ref. 1 using the results displayed in Fig. 1. Here we use similar reasoning to obtain a model for the - -

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS log

A / L3

~

2

=

°

D

log

(a)

e/L

,og A/L 3

Ae/d Ao/1~

F InnerCutoff ~ , ~ l o p e =2- D ~ , , ~ Outer Cutoff I eI

/L

~/L

(b)

log e/L

Fig. 1. Variation with measurement scale of measured area in a cubic volume of side L for a surface exhibiting fractal character. (a) No inner or outer cutoffs. (b) Inner and outer cutoffs to fractal behavior.

local mean flamelet area per unit volume, and hence the local mean rate of fuel consumption. Then normal and oblique flames are analyzed with the proposed rate model under the assumption that density change is negligible, and a burning velocity expression is obtained that is identical to that of Ref. 1. A numerical analysis of anchored flames in grid turbulence is performed with the new closure model and with a conditioned, second-order transport model developed by Chen, Lumley, and Gouldin [4]. The calculation results are realistic and consistent with those of Ref. 1. Finally, flame calculation results with a gradient transport model for turbulent species flux are compared to results from the second order model. The gradient model predictions are not in agreement with those of the second-order model.

243

where CR is a model "constant," UL is the local laminar burning velocity of the flamelet, and A Yf is the change in fuel mass fraction across the flamelet. The notation ()f represents an average over the flamelet surface, 1~ is the local, turbulent flame brush thickness, while lc will be defined below. Finally, Po is the cold gas density and (c) is the mean degree of reaction, where c = (Yf Yf°)/A Yf and yf0 is the initial fuel mass fraction. This expression is justified using the following arguments. Consider a cubic volume centered in space on a point x with side lc that is small enough that the ensemble mean conditions inside the volume can be assumed to be uniform with a mean degree of reaction (c). Assume that the flamelet geometry can be represented by a fractal surface of fractal dimension D, inner cutoff ei, and outer cutoff co. Here, eo is representative of the larger scales of wrinkling of the flamelet surface, and for the flamelet surface to be distributed on average in a uniform manner over l¢3, we require that 1~ _> ~0. Let Pc be the probability that a flamelet is present in the volume l~3. According to Mandelbrot [1, 3] a constant property surface (e.g., a constant temperature surface) in homogenous, isotropic turbulence is fractal. In a cube of volume L 3, the mean surface area per unit volume of a fractal surface measured at a length scale e is [1] A / L 3 = C ( e / L ) - ~ e2 / L ~,

where C is a constant termed by Mandelbrot the lacunarity [2]. As noted, Gouldin [1] has argued that if there is a characteristic minimum length scale of wrinkling, which is defined as the inner cutoff, the actual mean surface area per unit volume is finite and can be written as A / L 3 = C ( e . i / L ) -D f:i2 / Z 3.

C H E M I C A L CLOSURE M O D E L The proposed rate expression for fuel consumption is of the form (6Of) = CRP0 ( A YfUL)f(ei/eo)2-Dlc -1

x (c)(1 - (C))lc/IF,

(2)

Mandelbrot considered homogeneous, isotropic turbulence. Then L 3 is a volume embedded in turbulent fluid and the surface is dispersed throughout this volume in a statistically uniform manner. Flamelets are not expected to fdl space in the same manner as Mandelbrot's constant prop-

244

F. C. GOULDIN ET AL.

erty surfaces. To adapt Mandelbrot's fractal surface concepts to the present problem it is assumed that, on average, when there is a flamelet present in lc3, it fills lc3 in the same manner as a constant property surface in isotropic, homogeneous turbulence fills L 3. Thus the average flamelet surface area in lc3 can be expressed as Ac = C(Ei/lc) - O ¢ i 2 P c ,

(3)

where Pc, as noted, is the probability of finding a flamelet in the volume Ic3. For modeling we propose that the mean rate of fuel consumption can be written as a product of A~ and the mean rate of fuel consumption per unit flamelet area--po(A YfUL)f. Flamelet surface wrinkling and other effects of turbulence will perturb the flamelet structure and therefore alter fuel consumption rates. PoAYfuL represents the local instantaneous consumption rate in terms of the reactant density, the change in fuel mass fraction across the flamelet and the laminar burning velocity. This expression may be interpreted as a definition for the laminar burning velocity. It is implicit in our modeling approach that good estimates of po(A YfUL)f can be obtained and our focus is on the modeling ofA¢. A fuller discussion of possible modeling for po(AYfUL)f is given in Ref. 1. With our expression for mean fuel consumption per unit area, the mean consumption rate in lc3 is just po( A YfuL)fAc and (~f) --p0( A YfML)f~.i2-OlcO-3Pc, with Pc, el, and lc yet to be determined. There is experimental evidence for v-flames [5] that the probability density function describing flamelet location on a line that is normal to the flame brush is Gaussian. If flamelet motion is the result of many random impulses the central limit theorem also suggests that this probability will be Ganssian. For thin flamelets a reasonable approximation [5] to the value of the Gaussian function at a point x is given by Cn(c(x))(1 - (c(x))), where Cn is a normalization constant. For clarity the functional dependency of (c) on position, x, is shown explicitly. The probability of finding the flamelet along a segment of the normal which is of length lc is just Cn(c)(1 - (c))lc. The integrated

probability over the flame brush thickness, IF, should be one, and therefore Cn - IF- l, where the constant of proportionality depends on the definition for IF. For Pc it follows that Pc - (c) (1 - (c)) lc ~IF.

(4)

The inner cutoff is a measure of the smallest scales of wrinkling. For unit Schmidt number we expect that the smallest possible scale of wrinkling due to turbulent velocity fluctuations is related to the Kolmogorov microscale. Gouldin [1] has hypothesized that for flamelets the inner cutoff is equal to or greater than the Kolmogorov scale and that the ratio of the two scales depends on u ' / (/dE) f. He suggests an expression for ei [1]-(:i ~- q / f =

)7/[1 - (1 - A t - l/4R1-3/4)

x exp( - At l/4Rl - 1/4u'/(UL)

)],

(5)

where )7 is the Kolmogorov microscale, Rm the turbulent Reynolds number, and u ' the root mean square velocity fluctuations. )7 represents the minimum possible scale of surface wrinkling a n d f accounts for smoothing at the smallest scales of wrinkling due to fiamelet propagation relative to the reactant flow. Although this relationship for f is speculative and future research may well lead to an alternate expression for ei, the above form is adopted here. The arguments in Ref. 1 suggest that lc should be a multiple of l, the turbulence integral scale, if not equal to i. As noted, the arguments leading to Eq. 3 are based on the assumption that, on average, when present the wrinkled flamelet surface fills lc3 as Mandelbrot's fractal, constant property surface fills L 3. In order for this to be true we expect lc ~ l, where I is taken as the outer cutoff (l = ~o). Our argument here is based on the assumption that, although the flamelet can extend over large distances in two directions, its extent in the third direction is limited. At any instant, flamelet sheets can be bounded by two smooth surfaces with a local separation distance that is proportional to the outer cutoff, which in turn is assumed to be I. If flamelet wrinkling is such that i Peters [17] has suggestedthat the Gibsonscalebe used as the inner cutoff;see Ref. 1 for a furtherdiscussionof the Gibson scale.

245

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS

l o

Fig. 2. Schematic representation of fractal flamelet curve generated by the intersection of a flamelet surface with a smooth planer surface.

flamelets double back on themselves many times, the separation distance of these two bounding surfaces could be larger than 1, and then lc would be larger than the outer cutoff. We assume that such doubling back on large scales is not important. Doubling back on scales up to and including the order of the outer cutoff is expected and accounted for by the fractal description of flamelets. These structures are illustrated qualitatively in the sketch of Fig. 2. It is apparent from Fig. 2 that there are scales of flamelet curvature that are much larger than 1. These scales are assumed to be nonfractal and their effect on (oJf) is accounted for by the Pc term in the closure model. At this point a digression on the length scales introduced in this work may be helpful. ~/and I are traditional turbulence length scales and are widely used for characterizing turbulence. Ei and eo are scales introduced with the fractal description of rough surfaces, and they are related by hypothesis (which is supported in part by experiment [9, 18]) to the turbulence length scales, ~/ and 1. lc is introduced for purposes of exposition, while IF is the mean thickness of the turbulent flame brush (Fig. 2).

A major difficulty in adapting Mandelbrot's results [3] for surfaces in turbulent flow to the present problem arises because Mandelbrot assumes isotropic surface distributions while our flamelet distributions are highly anisotropic. To overcome this problem we focus attention on a small volume, lc3. When the flamelet is in this volume its distribution in the volume is nearly isotropic, and an ensemble of such occurrences should have an even more isotropic distribution. Clearly, for the arguments leading to Eq. 4 to hold lc must be larger than the scales of surface wrinkling or the surface contained in lc3 will not be representative of all scales of wrinkling. At the same time lc cannot be too large or the isotropic character will be lost. Thus/c = 1 = eo appears to be the obvious choice. This choice for lc implies that IF > 1, which may surprise some people.2 Instantaneously the flamelet is not wrinkled sufficiently to occupy a volume comparable to the entire turbulent flame brush. Instead it occupies at 2 Note that, for a v-flame stabilized on a transverse rod of diameter small in comparison with !, the condition iv > I will not be met for the initial region of flame development.

246

F. C. G O U L D I N ET AL.

any instant a smaller volume, and it occupies the entire flame brush only in an average sense. In summary, the final expression for (o:f) is

f(lf/,7) ° -

2IF - '

(1 -
Assume that

f,~, f and /c are uniform in space and dependent on Conditions in the reactants. Also note that the fractal modeling in Ref. 1 gives UtI
f ~- ([f/17) D - 2 . It f o l l o w s t h a t

(2 ')

=CR

tF-'

f is given by f = [1 - (1 - A t -'I4R1-3/4) x exp(-Atl/4Ri

×

ANALYSES OF NORMAL AND OBLIQUE FLAMES Normal and oblique flames are analyzed using the proposed rate model as presented in Eq. 2 ' . Results for the burning velocity are compared to the burning velocity model o f Ref. 1. For a normal flame assuming high Reynolds number, one-dimensional flow in the mean and A yf constant (for convenience), one can develop an equation for
aoutd
(6)

where ut is the mean reactant flow velocity approaching the flame brush or time mean reaction zone and is the turbulent burning velocity. Primes denote fluctuations in c and p, and U and u are the instantaneous and fluctuating velocities, respectively. If this equation is integrated from - oo to +oo with
I"


Yf dx,

and substituting for <~f) P0Ut = p0

i

(7)

-I/4U'I
Also note that for homogeneous, isotropic turbulence l/~l = h t l l 4 R i 314 [1].

poUt=

(l-)ax.

For this relationship to be consistent with the results of Ref. 1 it is required that lc = l, as suggested in the previous section, and

CR=IF/ ITo. (C>(1--(c>) dx.

(8)

Clearly,

I)I~
'UfT is a constant for a given geometry and initial conditions, and the application of the present model to the normal flame geometry is selfconsistent with previous work [1]. 3 3 A referee has noted that the Koimogorov-Piskunov-Petrovski theory [21J for gradient transport with a constant turbulent diffusivity scaled as u 'l gives Ut///o -- (ill'q) (/3-2)/2(U ' IMo)lt2

with the proposed closure model provided that IF -- I. If this relationship between I and IF is not assumed, the expression for ut/Uois Ut//,/° ~ (f//~)(D- 2)/2(/./#//3o) I/2(W/F) 1/2.

Furthermore, if one assumes that #,/IF- (fll)i)D- 2uol u ", then

oa

CRf(IcfI~)D-2
x ( I - (c>)IF-i

dx.

It is implicit in the formulation presented in the previous section that most of the spatial variation in ), i.e., it is in the probability of fmding a flamelet in l~3 at x.

utluo -

07/,)) D- 2,

as obtained here. A physical interpretation of the above expression for I/IF can be given. The arguments used to develop the closure model imply that l/(fl/~)D-2Uois an appropriate time scale for the mean fuel consumption rate, while l/u' is a typical turbulence time scale. Thus it follows that I/IFis proportional to a ratio of turbulent time scale to mean chemical reaction rate time scale.

247

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS

approximated locally by plane surfaces and, to first order, that the flame brush thickness grows linearly with distance along the flame brush. Finally, we assume slow growth in flame thickness and parabolic flow. The geometry to be studied is shown in Fig. 3. Note that by definition the ~"axis lies along (c) = 0.5, and ~ is perpendicular to ~', while the x and y axes lie along and perpendicular to the approach flow, respectively. A control volume analysis is pursued and the net mean flux of products out of the control volume is set equal to the mean rate of production in the volume. For a rectangular volume of unit depth, the rate of product formation can be written as

In view of the preceding discussion it makes sense to define IF by ~

(c)(1-{c))dx.

With this definition for IF, Cn is just IF-1, and if the probability density function of flamelet position is Gaussian, one finds that IF = (w/8) l/2cr, where o is the standard deviation of the flamelet position probability density function. With the above definition for flame brush thickness Eq. 2' may be written as

,o(a

(I - (c)).

The above analysis for normal flames can be extended to oblique flames without great difficulty and it is interesting to see if the notion of a turbulent burning velocity may be applied to oblique flames. As above, assume constant density and, hence, straight streamlines. We are interested in cases where the flame brush thickness increases in the streamwise direction as is observed experimentally, and we assume a priori that the mean flame brush structure follows similarity at least locally. We assume that (c) constant surfaces can be

where ~a and ~p define the edges of the flame brush and p is the ratio of the amount of product mass formed per unit mass of fuel reacted. By referring to the above model for (c0f/ one notes that, if (c / satisfies similarity, that is, (c) is a function only of 7//(~" + ~'0), where ~'0 is a suitably chosen constant, and eo/ei is constant, the integral of (~0f) over 7/ is a constant and production can be written as P = ~pI~.

= 0

Reactants

~-/2 pu

~UCts

,,. ~

0p

l

~,¢.¢>

-I

Fig. 3. Schematicdiagram of flame brush showingcoordinatesystemsused in the similarity analysis.

248

F. C. GOULDIN ET AL.

where I,~

=

I 'IR (wf)

can obtain such a result in two ways: assuming either Ic = 1/2 or A0 = 0. In either case the same expression for ut is obtained.

d~,

t/p

Out = sin OOu = p I , o. which is constant, i.e., independent of g'. For the product fluxes into and out of the control volume only convection is considered for faces 1 and 2, and turbulent fluxes across the other two faces are zero since conditions are uniform in the vicinity of both faces. Thus only convective fluxes need be estimated and no turbulence modeling is needed at this point. The neglect of turbulent fluxes of product across faces 1 and 2 implies a boundary layer approximation, namely, parabolic flow and small dlF/d~. Product flux into the control volume is M i , = o ( U ) cos O[IF(1)I~+tan(O+OF)6], where similarity has been used and

(10)

Thus a meaningful burning velocity, defined as the component of the mean approach flow velocity perpendicular to ~', is obtained and is the integral of the mean product production rate across the flame brush. The control volume width, 6, does not appear in the final expression for ut, as would be expected. However, implicit in the analysis is the assumption that similarity is valid for a control volume with finite & Using the closure model presented above one can evaluate I,~ and find by substitution Ut =

CR(UL)f(fl/tl)D-2plR.

(11)

where

to = e],, (c) d,lltF. tip

The angles 0 and 0p are defined in Fig. 3 and IF(l) denotes the flame brush thickness at 1. For the flux out one obtains

Mo~t = o ( U) cos O[ IF(2) I~ + tan 0 6]. The second term in this expression is the flux across the control volume face, which is totally in products. Equating the difference of Mo,t and Mi, to production one obtains

6plo, = cos 0O (U) [( IF(2) - IF(1)) I~ + (tan 0 - tan(0 + 0v))~]. For unit Schmidt number and (oJf) an even function of ~t, by symmetry OR -- 0 = 0 + 0p = A0. Then

plo,= sin 00(U)[1 + ( 2 L - 1) tan A0/tan 0],

(9)

where we have made use of the relationship (/F(2) - /F(1)) = 26 tan A0. Without chemical reaction, symmetry requires that 0 = 0 and Ic = 1/2. Here we look for a solution that gives a meaningful definition for the turbulent burning velocity. One

IR= f ~YIp (c)(1 - (c)) d,1/t~. From the analysis of normal flames CR = 1/Is, and one sees that u t as defined here for oblique flames is equal to that for normal flames and that ut is independent of ~" provided similarity holds, which implies nondecaying turbulence. The major assumption used to obtain this result is that of similarity. This assumption plays a major role in our evaluation of Min and Mout. It should be noted, on the other hand, that the constancy of I,~ is inherent in the modeling of (oaf), provided Co/el is constant across the flame brush, and does not require an additional assumption. The quantity (c)(1 - (C))/IF is proportional to the probability of finding the flamelet at a point where the mean degree of reaction is (c), and thus the integral I s is constant because the probability of finding the flamelet somewhere across the flame brush is unity. The condition A0 = 0 corresponds to a constant flame brush thickness, which is contrary to observations, while with Ic = 1/2 the thickness continues to increase with increasing ~'. If (o~f) is an even function of ~/, by symmetry, we expect Ic,

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS which is the mean value of (c) when averaged across the flame brush, to be approximately one half. (N.B.: our numerical calculations give Ic --1/2.) A continually growing flame brush might be thought to imply an increasing combustion rate because of an increase in flamelet area. However, as noted, our rate model scales (of) by IF, and 1,0 is independent of ~'. This scaling arises via Eq. 4 according to which the probability Pc of finding a flamelet in a volume lc3 decreases as ~"increases. A physical interpretation of this behavior is that the flame brush thickness grows with increasing ~" due to the cumulative effects of the many, random upstream displacements of the flamelet surface(s) in a process that may have some likeness to a random walk. These displacements cause the flame brush to grow but they do not increase flamelet area and thereby combustion rate. Flamelet area is increased by the fractal wrinkling of flamelet surfaces due to the reactant flow turbulence. This area generation is dependent on lf/~l and not the flame brush thickness. NUMERICAL ANALYSIS Results of calculations using the proposed rate model of Eq. 2 ' and the second-order transport model with conditional statistics of Chen, Lumley, and Gouldin [4] are now presented and discussed. The transport model is based on model equations for conditioned moments, where conditioning is for being either in reactants, c = 0, or in products, c = 1. The instantaneous progress variable acts as the indicator function. An equation for the mean progress variable, for negligible density change and high Reynolds number, can be written [4] as o(c)/ot+ (u) •

v(c)

249

~'J



{((U)-(U)(I))(c)}+(6of)/AYf,

Fig. 4. Sketch of flow configuration assumed for the numerical calculations. consequently, gradient transport modeling is avoided. These model equations (consult Ref. 4 for more detail of the model) are combined with the new rate model, and calculations are performed for an oblique flame in decaying grid turbulence. The flow and flame configuration are shown schematically in Fig. 4. This is the same configuration as considered in Ref. 4, and comparison of the present results with those in Ref. 4 can be made by the reader to see the influence of changing rate models. For calculations the flow is assumed to be parabolic. Upstream boundary conditions and model constants for the calculations are given in Table 1; note that D = 2.4.. The upstream value of IF is taken to be the integral scale. For these calculations IF is defined as the local TABLE 1

Standard inlet conditions: (u) = 7 nes t21(V) = 0.07 12 = 49 cm/s l, = 7.27 mma u, = 1 cm/s (laminar flamespeed)°

(12)

where (U/(1) is the conditioned velocity in the products, and the first term on the right-hand side of the equation is the turbulent flux of products. This term would be given by - V - D t V ( c / if a gradient transport model were used, Dt being the turbulent diffusivity. Instead we solve a model equation for (U/o) as well as finding (U/, and

Flamelet

<(::>=I

v =

=V"

~

Rl =

1.85 x

I0 -s m21s

229

M = 0.025 m Model constants:

CR = 4.0 D = 2.4 A, = 0.37 [1] Ca = 0.09 [7] For other constants see Ref. 4

a r denotes standard reference value.

F. C. GOULDIN ET AL.

250 -2

10

Model Predictions

163

10 ~

I

I

I

I

I

10

I

I

I

1

10

2

x/M

Fig. 5. Test of conditional, second-order turbulence model against data for decaying grid turbulence obtained by Warhaft [6]. (a) Decay of normal stresses in flow direction.

separation distance between the (c) = 0.1 and (c) = 0.9 surfaces. It can be shown, for a Gaussian distribution of flamelet position, that this distance is 4.15 times the thickness defined, as proposed in the previous section, by the integral of/c)(1 (c)) across the flame brush. It follows that, if the former definition of iF is used in the rate expression of Eq. 2 ', Ca will be 4.15 instead of 1.0, as suggested in Eq. 8. As noted below, a value of 4.0 for Ca is found to give good results in the numerical calculations when IF is defined by the distance between the 0.1 and 0.9 (c) constant surfaces, a value that is within 5 % of the expected value. The difference in these two values may be due to the fact that (c)(1 - (c))/IF is an approximate expression for the flamelet position probability density function measured in Ref. 5. As a test of the turbulence model for nonreacting flow, the decay of grid turbulence was calculated, and the results for (U2), 1, and Ri compared to the experimental data of Warhaft [6]. Comparisons are made in Fig. 5 where it can be seen that model predictions for (u 2) and Rl are very good. If the assumptions of the analysis of the previous section are valid one would expect that the turbulent burning velocity defined as the component of the mean reactant flow velocity normal to

the (c) = 0.5 surface would equal that predicted by the fractal model of Ref. 1; see Eq. 11. In Fig. 6 calculation results for ut/uo are compared to those predicted by the burning velocity model. For this comparison Ca is chosen to give the best overall agreement; CR = 4. The comparisons shown in Fig. 6 for a broad range of conditions are quite good both with regard to magnitude and variation of ut/uo with x/M. We attribute the small observed differences, which are present, between burning velocity model and calculation results to a departure from flow similarity that is assumed in the analysis of the last section. A lack of similarity can be seen in the variation of lateral (c) profiles with x/M, and several typical profiles are shown in Fig. 7. By similarity one would expect to see iF grow linearly with x / M, but such behavior is not found in the numerical results. Instead iF is found to grow very rapidly at small x / M and then to grow much more slowly at larger x/M; see Fig. 8, which presents IF/X versus

x/M. In addition a similarity analysis of Eq. 12 can be performed and a result equivalent of Eq. 11 is obtained. For stationary flow and (U) = constant, Eq. 12 becomes v • [{u)"J{c)]---

re.

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS

251

200.0

180.0

ModelP r e d l c t l o n s 160.0 re 140.0

120.0

100.0

I

,

20.0

0.0

I

,

40.0

I

,

60.0

'

'

80.0

100.0

x/M

Fig. 5.b. Variation of turbulence Reynolds number. M is the scale of the turbulence generating grid; see Table 1.

10.0 Calculation

~\

results

BumingVelocitym o d e l

8.0

~" ~

~

RI = 2290

6.0

4.o

Ii

~

g 2.O

I

I I I I '

0.0 0.0

RRI - 229

,_

~

~

'

'

~

'

~ 50.0

Rl= 22.9

. . . .

I 100.0

(a)

. . . .

i 150.0

. . . . 200.0

X~

Fig. 6. Comparison of calculated burning velocity with values predicted by the fractal burning velocity model for conditions given in Table 1 and for (a) various Reynolds number, (b) various laminar burning velocities, and (c) various initial integral length scales. In (a) the kinematic viscosity only is varied, whereas in (c) the kinematic viscosity is varied to maintain a constant Reynolds number. Ur is arbitrarily set to 1 cm/s.

252

F. C. G O U L D I N ET AL. 6.0

5.0

~ 2 u

Calculation results

~.-£ ~.--~-~.:~.~_'.:

Burning Velocity Model

r 4.0

3.0

2.0 i

1.0 i

o

0.0

i

0.0

I

50.0

'

'

*

'

t

'

100.0

'

'

'

I

150.0

200.0

X/M (b) 6.o

5.0

4.0

~

2offi10 I r

2.O x"x. "

20 ffi 2 r 110 ""'"

1.0

0.0

i 0.0

. . . .

50.0

I

. . . .

100.0

i

'

150.0

200.0

(c) Fig. 6. (Continued).

For similarity (U)o) = h(y/If) with components hx and hy; (c) = g(y/lf) and (of) = r(y/lf). Substitution into the above gives a

a

~x ( hxg) + ~y ( hyg)= r/A YF It follows that

h" dtf+ (hyg+hyg ') -(xg+hxg')(y/lf) dx = l f r / A YF,

where the primes denote differentiation with respect to y/lF. From the later equation it can be seen that for similarity dlF/dx = constant and r must scale as IF- l, as is the case. Also note that for dlE/dx = constant the line (c) = 0.5, i.e., the ~" axis, is straight and 0 is constant. Now let Y/IF = )7 and integrate across the flame brush at some point, x = constant. One obtains

dlf/dx I (hxg)'.~ d]+ I =(lf/A Yf) ~ r dfi.

(hyg)' d.)7 (13)

253

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS ....... ..°..°.. --

1.0

x/M = ""'""

0.8

•.... %50

A

0.6

".%

0A ~t "°°.,

%%'-. ~t

0.2

°'%

\ \ "....

\% 0.0

°%,, "",, %'.'N, """

,

....... L

.0.05

-0.03

-0.01

0.01

0.03

0.05

y/x-(y/x) .(C> = 0.5 Fig. 7. Lateral profiles of mean progress variable at three streamwise locations for the reference conditions given in Table 1. (y/x)/~/=05 denotes the value of y/ x at (c) = 0.5.

For parabolic flow (0 small)

constant for given initial conditions since for similarity the left-hand side of Eq. 13 is constant. Combining this result with results of the previous control volume analysis one obtains out = 1,0. Thus, a similarity analysis of the model equa-

I (wf) dr=I~=lf l r d.v ( y and r/ nearly coincide) which one sees is a 0.06

i/s'

0.03

~.

i7

~'~

R

2290

0.04

.J" 0.03

0.00

'

0.0

'

'

'

I

50.0

'

'

'

'

i

100.0

'

'

'

'

I

150.0

'

'

'

'

200.0

XJM Fig. 8. Streamwise variation in calculated flame brush thickness for different Reynolds number. Conditions are the same as in Fig. 5a.

F. C. GOULDIN ET AL.

254 6.0

5.0

4.0 I

~

I

3.0

Gradient Transport

°

...........

. .....

. ..........

2.0 I

I I

1.0

I

i

0.0

0.0

,

i

,

I

.

,

,

$0.0

,

I

100.0

,

,

•

•

I

150.0

.

.

.

.

200.0

X/M

Fig. 9. Comparisonof burning velocitiescomputedusing the full, conditional,second-order transport model, a gradienttransportmodeland the fractalburningvelocitymodel. Conditions are the referenceconditions given in Table 1. Note that all three models give very similar results. tion for (c) shows that I~ and both dlf/dx and Oare independent of x and that the control volume analysis of the previous section is applicable to the model equation. With similarity the model equation gives the expected expression for burning velocity. In performing the calculations for this work some results were also obtained with a gradient transport model for the (c) equation while a full second order model was retained for the turbulent Reynold stresses. For these calculations Dt = Cdk2/ed, where k is the turbulence kinetic energy and ed is the turbulence energy dissipation rate; the constant Co is set equal to 0.09 [7]. We find that the results for burning velocity are similar for the two transport models (Fig. 9), but that there is a marked difference in IF (Fig. 10). From Fig. 11, it can be seen that gradient transport greatly under predicts turbulent fluxes in comparison with the second order model. In addition, countergradient diffusion is not observed. Libby et al. [8] have shown that the interaction between the mean pressure gradient and density fluctuations plays a crucial role in causing countergradient diffusion.

Pressure gradient fluctuations may also contribute to this phenomenon. Thus it is not surprising that when combustion with negligible density change is considered countergradient diffusion is not observed. That both transport models give essentially the same burning velocity may at first be surprising. The reason for this behavior is as follows. From Eqs, 9 and 10 we see that with similarity put is given by pIo, for Ic = 1/2. As discussed above, with our rate model I~ is dependent only on (po A YfUL)(fl/~I) 0-2 and hence independent of the species transport model. This result follows from our expression for Pc and its scaling by IF- 1. Then since both transport models give Ic = 1/2, they give similar ut values. Conditional and unconditional moments of velocity are presented in Figs. 12 and 13. It can be seen that unconditional mean velocities in the transverse ( y ) direction remain close to zero. They should equal zero in the absence of density change. Conditional mean velocities (Fig. 12) for the burned and unburned zones are of equal magnitude but opposite sign, and their difference

255

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS 0.06

0.05

/ I..

0.04

~.... Condlllonal Model

I

:,< 0.03

I

0.02

"°'-.~

....

. ...........................

0.01

Gradient T r a n s p o r t

0.00

, 0.0

,

I

.

.

.

.

50.0

I 100.0

,

,

,

,

I 150.0

,

,

,

, 200.0

xJM

Fig. I0. Comparisonof flamebrush thicknessas calculatedusingthe fullconditionaltransport model and the gradient transport model. decreases in the flow direction as the turbulence decays. The relative magnitude and sign of the conditional transverse velocities must undergo significant change in the presence of combustion with heat release and density change if countergradient diffusion is to be observed. Transverse turbulent normal stress profdes are presented in Fig. 13. The unconditioned stresses are nearly uniform, as is required in the absence of density change, while conditional stresses vary significantly across the flame brush. From the above presentation and the figures it is clear that the new rate model gives consistent and plausible results. Predicted burning velocities are in the expected range. However, exact correspondence to the burning velocity model of the last section is not achieved. This lack of agreement is apparently the result of the departure from flow similarity in the numerical results, which in turn we attribute to the decay of turbulence. In the model calculations a constant density has been assumed and no countergradient diffusion is observed. At present the conditional modeling approach has not been adapted to account for density changes across the flamelet. Work to allow for density change is needed.

SUMMARY AND CLOSING COMMENTS A closure model for the mean fuel consumption rate in premixed turbulent flames has been proposed and applied to several test problems. The model expresses the mean rate in terms of the product of the mean fuel consumption rate per unit flamelet surface area and of an expression for the mean flamelet area. The latter expression is obtained under the assumption that flamelet surface geometry can be represented by fractal surfaces. In addition the model uses an empirical expression from Ref. 5 for the probability of finding the flamelet in a volume 43. The new parameter appearing with the model is the fractal dimension. Based on the results reported in Ref. 1 and experimental [9, 10, 19] and theoretical results [11-13] all for high turbulence Reynolds number we would recommend a value for D between approximately 2.33 and 2.4. On the other hand results of a preliminary nature for low Reynolds number premixed flames [16-18, 20] indicate a D between 2.1 and 2.2. Thus it appears that D may be a function of Reynolds number at least for low turbulence Reynolds number. Test problem results are mutually consistent and

256

F. C. GOULDIN ET AL. $.0

o°..°" °.. oo".....°

5.0 %

,eA

/

4.0

",.

° - x / M = 200

/ /

v A "o > V

,,.. ,,

3.0

/

lOO

-

.

,

.......

2.0 1.0 0.0 -0.05

•0.03

-0.01

0.01

ylx - (y/X)

< ©> =

0.03

0.05

0.5

(a) 0.0

5.0

&

4,0

v A

3.0 ...,. x/M = 200

V

,."

°',

.. 2.0

.:

1.0 0.0

-0.05

~ - ' " ' ~

,.

i . . . . . . l .........

-0.03

"'"'"- -

-0.01

0.01

y l x - (y/x) < c > :

I

0.03

t

0.05

0.5

(b)

Fig. 11. Comparison of lateral turbulent fluxes of progress variable as calculated for the reference conditions by (a) the full conditional model and (b) the gradient transport model.

plausible. The rate model is easy to use and gives burning velocity values that are equal to the predictions of a new burning velocity model [1] that also uses fractal concepts. The oblique flame analysis suggests that a

meaningful burning velocity can be defined for oblique flames, and that ut values can be related to the traditional definition of the burning velocity of normal flames. However, the oblique flame analysis rests on the assumption of similarity, and the

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS

257

0.3

0.2

0.1

> v

~

0.0

~

~

~

O)

'*'"'-°... ..... . ......... "-... .... ...........................

41.1

< V >(°)

..... ~°*---....,.......°.°....°°°°

°'%. -0.2

•"0.3

I 0.0

I 0.2

I

I 0.4

,

I

I

0.6

I

I 1.0

0.8



Fig. 12. Conditional and unconditional lateral mean velocities at x / M = 50 calculated for the reference conditions in Table 1. Conditional mean velocities in the reactants, (V) (0), are nearly equal in magnitude to those in products, ( V ) % as expected.

36.0

< w >1< U >2 32.0

~"

e~ ^ v

A

%".°°

28.0

/

°%

< w>(1)/< U >:

"%%°°°° %°'°° °°'"°'*,.° °

24.O

V

/

/J ~" ~

J"

20.0

16.0

- "

/

...............

< w >(mI< U >2

/

12.0

I

0.0

I 0.2

i

I

i

OA

I 0.6

i

I 0.8

1.0


Fig. ]3. Conditional and unconditional lateral normal su'esses calculated at x / M = 50 for reference conditions. The unconditioned normal stresses are essentially uniform across the flow as expected. Note that (0) and (1) denote conditioning in reactants and products respectively.

258 expression for burning velocity requires a supplemental condition--either A0 = 0 or lc = 1/2. The numerical calculations of the last section show that departures from similarity can affect the burning velocity variation. In the absence of similarity it may be possible to develop an alternative, useful definition for the burning velocity. Until such a definition has been established, turbulent burning velocity data should be interpreted with caution, especially if flow similarity is not observed. The empirical basis for the relationship Pc (c)(1 - (c)) can and should be made more substantial. For example the passage time model of Ref. 14 suggests P - (c)(1 - (c)) but not the normalization b y / F - l . Furthermore, in Ref. 5 the turbulence Reynolds number is low and multiple flamelet crossings of lines parallel to the ~ axis are rare events. It is reasonable to expect that future work will lead to an improved expression of Pc. Developments of the passage time analyses of Bray, Libby, and Moss [14, 15] may help in this respect. The possible occurrence of large-scale doubling back, as described in the second section needs further study. If such events are important it may be necessary to change our scaling of lc and our expression for Pc. To avoid confusion it should be noted that the neglect of large-scale doubling back does not preclude doubling back on scales up to the order of the outer cutoff. Thus along lines parallel to the ~ axis multiple flamelet crossings at any instant will occur over a length of order of the outer cutoff. The region of these crossings will fluctuate in time and thus form a flame brush that is thicker than the outer cutoff. The ratio ~ / f has been used to represent the inner cutoff. There is some question as to the appropriateness of this representation [1, 16]. Our expression for ~ / f h a s been taken from Ref. 1 and the present work provides no test of its validity. Should a better representation of ei be developed it can easily be incorporated in the model. The fractal character of flamelets is open to direct experimental verification by methods such as those described by Sreenivasan and Maneveau [9]. Information on the fractal dimension and the cutoffs can also be obtained. Premixed turbulent, burner flames have been studied using imaging techniques. The results available to date are

F. C. G O U L D I N ET AL. preliminary [16-18, 20], and they are for low turbulence Reynolds numbers. While they show fractal behavior, as noted above inferred D values are low--between 2.1 and 2.2. In addition, data from a spark ignition engine are also available [19] that show that at high Reynolds a value for D of approximately 2.33 is approached. More measurements are required to substantiate the existence of fractal character for flamelets and to determine the Reynolds number dependency of D.

This research was initiated while the first author (FCG) was a Science and Engineering Research Council Visiting Fellow at the University Engineering Department, Cambridge, England. Support from the SERC and from the US A rm y Research Office (Contract number DAAG28-82-K-O187) is gratefully acknowledged. The work o f the second author (KNCB) was supported by a grant from the SERC, and this support is gratefully acknowledged. The work o f the third author (JVC) was supported by the United States Department o f Energy, Office o f Basic Energy Sciences, Division o f Chemical Sciences. REFERENCES 1. Gouldin, F. C., Combust. Flame 68:249-266 (1987). 2. Mandelbrot, B. B., The Fractal Geometry of Nature, Freeman, San Francisco, 1983. 3. Mandelbrot, B. B., J. FluidMech. 72:401-416 (1975). 4. Chen, J.-Y., Lumley, J. L., and Gouidin, F. C.,

Twenty-first Symposium (InternationaO on Combustion, The Combustion Institute, Pittsburgh, 1986, pp. 1483-1491. 5. Miles, P., and Gouldin, F. C., Simultaneous measurements of flamelet position and gas velocity in premixed turbulent flames, College of Engineering, Energy Report E-86-03, Cornell University, Ithaca, NY, 1986. Presented at 2nd ASME-JSME Thermal Engineering Conference, Honolulu, Hawaii, March 1987. 6. Warhaft, Z., J. FluidMech. 144:363-387 (1984). 7. Jones, W. P., and Launder, B. E., Int. J. Heat Mass Transfer 15:301 (1972). 8. Libby,P. A., and Bray, K. N. C., AIAA J. 19:205-213 (1981). 9. Sreenivasan, K. R., and Meneveau, C., J. FluidMech. 173:356-386 (1986). 10. Lovejoy,S., Science216:1185 (1982). 11. Hentsehel, H. G. E., and Procaccia, I., Phys. Rev. A 29:1461-1470 (1984).

CHEMICAL CLOSURE MODEL FOR FRACTAL FLAMELETS 12.

13. 14. 15. 16.

17.

18.

Kingdon, R. D., and Ball, R. C., The fractal dimension of clouds, Physics Department, University of Cambridge, submitted. Gouldin, F. C., A I A A J. 26:1405-1407 (1988). Bray, K. N. C., Libby, P. A., and Moss, J. B., Combust. Sci. Tech. 41:143 (1984). Bray, K. N. C., and Libby, P. A., Combust. Sci. Tech. 47:253-274 (1986). North, G. L., and Santavicca, D. A., Fractal analysis of premixed turbulent flame structure, presented at the Fall Technical Meeting of the Eastern Section of the Combustion Institute, San Juan, Puerto Rico, Dec. 15-17, 1986. Peters, N., Twenty-first Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1986, pp. 1231-1250. Gouldin, F. C., Hilton, S. M., and Lamb, T., Experi-

19. 20.

21.

259

mental evaluation of the fractal geometry of flamelets, presented at the 22rid Symposium (International) on Combustion, Seattle, Washington, August, 1988. Bracco, F. V., Princeton University, private communication, March 1988. Murayama, M., and Takano, T., Fractal-like character of flamelets in turbulent premixed combustion, presented at the 22rid Symposium (International) on Combustion, Seattle, Washington, August, 1988. Zel'dovich, Ya. B., Barenblatt, G. I., Librovich, V. B., and Makhviladze, G. M., The Mathematical Theory o f Combustion and Explosions, Consultants Bureau, New York, 1985.

Received 21 January 1988; revised 14 September 1988