Statistical model of turbulent premixed combustion with interacting flamelets

Statistical Model of Turbulent Premixed Combustion with Interacting Flamelets A. A. BURLUKA, M. A. GOROKHOVSKI, AND R. BORGHI* ESM-2, IMT-Technop~le de Chateau-Gombert, 13451 Marseille, C#dex 20, France A statistical model is proposed for the turbulent premixed combustion regime when chemical reactions occur in interacting thin layers (flamelets). The statistics of the reacting turbulent scalar field are described in terms of a reference scalar field (rsf) with an extra probability dimension. A transport equation for rsf is derived from that for the probability density function of the reacting scalar, and the presence of interacting flamelets is reckoned through a closure hypothesis on the conditionally averaged scalar dissipation. It appears that molecular diffusivity, chemistry, and flamelet interaction terms enter explicitly in the resulting rsf equation, and in this framework microscale mixing is described within a well-posed Cauchy problem. The model proposed has been applied to simulations of homogeneous mixing with and without chemical reactions. Predicted trends for scalar variance and the flatness factor quantatively follow the Direct Numerical Simulations (DNS) data. For the reacting case the rsf model results in a nearly bimodal probability density function shape. Simulations of V-shaped lean methane- and hydrogen-air flames have also been performed with the rsf equation coupled with the standard k - e model. The results fit the experimental data fairly well in spite of simple flow dynamics modelling within a boundary layer approximation. © 1997 by The Combustion Institute

INTRODUCTION A well-known problem of averaging of chemical reaction rates arises in turbulent reacting flow modelling. As early as at the beginning of the 1960s methods based on the balance equation for the probability density function (pdf) of temperature a n d / o r reactive species concentrations were proposed [1, 2]. These methods, in principle, treat the influence of turbulence on chemistry with exact nonclosed terms, and numerous turbulent combustion studies have been made within the pdf approach [3]. In the pdf transport equation there is microscale mixing (MSM) term caused by molecular diffusion in turbulent media. This term is responsible for pdf evolution in concentration space toward the perfect mixing state and it is nonclosed [2]. Many models have been proposed to describe MSM within the single-point pdf approach [4, 5], but it still remains a difficult problem [4]. The difficulty in closure of MSM terms is twofold. First, the conditionally averaged scalar dissipation (or diffusion), which is determined with two-point statistics of the turbulent scalar field [2], must be known. Second, the particular exact form of the MSM

* Corresponding author.

COMBUSTION AND FLAME 109:173 187 (1997) © 1997 by The Combustion Institute Published by Elsevier Science Inc.

term, which is written in terms of conditional dissipation, yields an ill-posed Cauchy problem for the pdf transport equation for any initial pdf distribution, but Gaussian [6], even where the conditional dissipation is known. This is why the exact expression for the MSM term is replaced with some closure formulae in existing models [4] and external information for turbulent scales is applied [3-6]. For the most part these models are physically sound for fully developed turbulence without consideration of any internal structures (e.g., flamelets). In [5] it was proposed to describe the MSM processes with some reference scalar field (rsD obeying the molecular diffusion equation, coupled with the block-inversion mechanism to represent the turbulence effects. In [7] a similar idea was proposed to simulate MSM with either a molecular diffusion equation or with an equation of "telegraph" type for rsf. In the present paper an equation for the reference scalar field with an extra probabilistic dimension is derived from the pdf transport equation itself. A simple hypothesis on conditional scalar dissipation is made for the case of combustion in interacting thin layers (flamelets), the structure of which is determined with molecular diffusion and chemistry. This approach avoids negative diffusivity in the Cauchy problem for MSM description and it is able to incorporate 0010-2180/97/$17.00 PI1 S0010-2180(96)00147-2

174

A . A . BURLUKA ET AL.

the physical characteristics of small structures (flamelets) into the modelling. Simulations of both the homogeneous case and realistic V-shaped premixed flames were performed to illustrate the capabilities of the model. For the latter case one can demonstrate reasonable agreement with the experimental data in spite of the rather approximate numerical model of turbulent flow. PHYSICAL REASONING BEHIND THE MODEL

Consider the turbulent premixed flame "brush" which consists of thin reacting layers, the thickness of which is much smaller than turbulent integral length scale. These layers form the interface between hot products and cold fresh mixture. For flow with sufficiently large Reynolds number the distance between the adjacent parts of the instantaneous flame front in the turbulent flame brush may be so small that those layers cannot be considered as independent of each other. Their proximity enhances heat and mass transfer processes and may change the structure of flamelets in comparison with an isolated laminar flame. For example, when a narrow tongue of fresh mixture is surrounded by the hot products, the preheat zones of neighbouring flamelets interact with each other, which changes the local conditions of flame propagation. When the narrow tongue of flame front is surrounded by

~

0

x

the cold fresh mixture, the heat transfer from beyond the reaction zone is locally augmented and the temperature of the products drops. Below, such processes will be referred to as the flamelets interaction. This physical picture is illustrated in Fig. 1, which presents the instantaneous scalar signal c(x, t o) along some line crossing the turbulent flame brush. The solid line shows the signal c(x, t o) which would exist if there were no flamelets interaction and the dashed line sketches the possible effects of flamelet interaction on both the cold and hot sides of the flame. It is clear that this interaction is the only phenomenon which is able to decrease the amplitude of concentration fluctuations in comparison with the case of an infinitely thin flame front. So it should play a principal role in the physics of microscale mixing. The occurrence of this interaction has not so far been considered in the so-called laminar flamelet models. On the contrary, in classical pdf models, which employ the notions of "fluid particles" and their interaction, the presence of flamelet-like structures has never been taken into account. In fact, in reacting flow with sufficiently fast chemistry both phenomena exist and it will be seen from the subsequent analysis that both should be taken into account. To describe the statistical properties of the scalar field within the turbulent flame the probability density function (pdf) description is

to)

Fig. 1. Instantaneous scalar signal c(x, t o) in the presence of reacting flamelets.

T U R B U L E N T P R E M I X E D COMBUSTION

175

usually employed. The "geometrical" interpretation of pdf [6, 8] can be written, after having chosen some x direction along the flame brush,

2 for the simple case corresponding to a single pass of the instantaneous flame front through the turbulent brush. From Eq. 2 one can approximately estimate that

P(c;x't°)

= ~fs~

dc - t d n c(x,t,,)=?

O c ( x , t ) ,,(x, fl ---Ox t)=? l

1]ec I

5-" l d~ c(x,to)=~~'~cOs

--

-

l

--

(1)

,~

where S~ is the area of isoscalar surface c(x, t 0) = ~ in some large volume V and l -(SjV) ~ is some turbulent length scale which is of the order of the integral length scale, ~ is the angle between the normal n to the instantaneous flame surface and chosen x direction, and /3 is an a priori unknown function related to the distribution of the angles a as well as the curvature and stretch effects. If local extinctions caused by stretch are neglected, /3 necessarily takes a value of the order of unity. This value is adopted below. Equation 1 means that the probability density P(~, t 0) is nothing but the mean part of the flame brush volume which is "occupied" by the scalar values ~ < c(x, t o) <__~ + dE. Furthermore, the integral probability distribution X ( ? , t o) can also be related to the volume share x(~) in which ~ < c ( x , t o) < 1:

Oc(X,t) OX

c(x.

t)=e

.

(3)

By using Eqs. 2 and 3, and assuming constant density and molecular diffusivity DM, the diffusion-reaction equation for the instantaneous field c( x, t ), Oc(x,t) ,)2c(x,t) - - DM• + W(c(x,t)), Ot cgx?

(4)

may be transformed into an equation for rsf c( X , t) like Oc(X,t)

3t

/32 O % ( X , t ) -- DM" l 2 t~X 2

(5)

+ W ( c ( X , t)),

where W ( c ) is an instantaneous reaction rate. The solution of an equation similar to Eq. 5 would give all the single-point statistics of the scalar turbulent field, e.g., ( c ) = f01c(X) d X , (c'") = [ 1 ( c ( X ) go

- ( c ) ) 2 dX.

X(~, t o) = 1 - f(~P(~, t o) d~ _ [~fl --

1

--

Tx(e).

J0

dc -1

Tl~l,,(~,,,,)=~

dc

/3

(2) i

i

¢

This relationship illustrates that the statistics of the real scalar field c(x, t) in physical space may be connected with some reference scalar field (rsf) c ( X , t) evolving in probabilistic space X. This rsf c ( X , t) is obtained from integral probability distribution X(~, t) on changing the dependent function and independent variable: c ( X , t) = X -l ( X ( ~ , t)).

All these transformations are illustrated in Fig.

Ax i(

~i

I

x

21,

Fig. 2. Illustration of geometrical treatment of probabilities.

176

A. A. B U R L U K A ET AL.

However one should remark that there are no terms in Eq. 5 which account for flamelet interactions, (e.g., for interaction of preheat zones of adjacent flame front filaments). Those terms should be inversely proportional to the integral time scale ~tc of the scalar field, as it is large-scale turbulence that curves and bends the flame front. Later, when the balance equation for rsf c(X, t) is derived rigorously from the pdf equation, those terms will appear as an "out-of-flamelet" contribution [9] in a rather simple form. Here it would be convenient to make only some qualitative estimations of the relative importance of this "out-of-flamelet" term which has an order of magnitude 1/ztc. The diffusion term in Eq. 5 may be estimated with Eq. 2 as

DM'92c(X,t) 12 ,gX 2

DM ('gZc(x,t)) -- l-T-- • l 2 ,gx 2 DM --

DM ~

m

Ac2 u'

12

l

,gP(?)

,9 2

,gt

,9

,gdz xeP( 6 ) - --~ W( ~) P( 6 ). (7)

Here Xe = (DM(VC(x, t))21c(x,t) = ~) is a scalar dissipation conditionally averaged at fixed concentration ~, W(~) is the instantaneous reaction rate, and ~ is an independent variable in composition space. Equation 7 is a diffusion-reaction equation of inverse-parabolic type for probability P(~), where the diffusivity - h e is always negative. The boundary conditions at the edges ~ = 0, ~ = 1 are unknown. The integral probability distribution

X(~) = 1 - foeP(~) de

Re

1 -

written as [2]

(6)

rtc '

where AC is the Taylor microscale for the scalar field, Re is the Reynolds number, and u' is turbulent velocity fluctuation level. Therefore, it can be seen that flamelet interaction would have the same order of magnitude as the flamelet diffusion-like term. In other words, relatively small changes (of temperatures, concentrations, etc.) which occur in thicker preheat zones occupying a rather greater volume, would have the same effects as the stronger changes in thinner reaction zones occupying a relatively small volume. MODEL FORMULATION Consider the case when the thermochemical state of flow can be described with the only variable c having the distribution P(6) (for a turbulent premixed flame c may be the progress variable). In particular, the Lewis number is supposed to be unity. A balance equation for P(c) for the homogeneous case is

is monotonically decreasing, bounded within 0 < X(~) < 1, and has no delta-peak singularities. X(~) is defined as the probability to meet at a given point (x-', t) the scalar value equal to or greater than c(X; 3, t). By integration of Eq. 7 from 0 to ~ the balance equation for X(~) becomes

`gx(~,) ,92x(,~) ,gx(~) ,9-----S--= x~ ,9~-~-- w(~) ,9~ Oxe ox(~) -

--

,9~

00~ '

(8)

where the last term is due to the nonuniform distribution of Xe over ~. A neglect of this term corresponds to the assumption of Gaussian pdf distributions [4, 6] only. However the Xe dependence ~ should be taken into account for reactive flow because of strongly non-Gaussian pdf shape [11]. Equation 8 should be supplied with boundary conditions X ( 0 ) = 1, X ( 1 ) = 0 resulting from the definition of X(d). In addition, to provide satisfying conservation of the mean scalar in the nonreacting case, the conditions of zero probabilistic flux Xe(aX(~))/O~ at the edges ~ = 0, ~ = 1 must be satisfied as well [12]. For an arbitrary X(6) distribution the

T U R B U L E N T P R E M I X E D COMBUSTION only way to satisfy both boundary conditions simultaneously is to require Xe = 0 at 6 = 0 and ~ = 1. Both the equations for pdf P(6) and for integral probability distribution X(3) are exact, but the microscale mixing (MSM) terms containing Xe require closure assumptions [4]. Among the proposed MSM models (see [4]) only the mapping closure [13] uses the notion of conditional scalar dissipation, while the others, such as the coalescence-dispersion models [4, 10, 14] and relaxation "return-to-mean" models [1, 3] avoid it by reformulating Eq. 7 to a form in which Xe is no longer applied. However, the mapping closure [13] has an applicability constrained to the Gaussian ensemble of scalar field realisations, so it seems ill-suited for the premixed turbulent combustion where the pdfs reveal a strongly bimodal structure [11,151. For this particular combustion regime an equation for the reference scalar field (rsf) c(X, t) with an extra probability dimension X is derived below from Eqs. 7 and 8. The pdf P(~) and rsf c(X, t) are related in a straightforward manner and contain the same statistics on the turbulent scalar field. Consider one-to-one reverse mapping c(X) = X - I ( x ( ~ ) ) , where X becomes an independent variable in some probabilistic space. In this space all the mean values may be obtained as ( f ( c ) ) = jdf(c(X))dX, e.g., (W(c)> = /1W(c(X)) dX. An example of such one-to-one reverse mapping c ( X ) = X - I ( x ( ~ ) ) for a bimodal pdf is depicted in Fig. 3. The values of at which c(X) = 6 within some X interval of nonzero length correspond to the 6-peak singularities in pdf (P(~). The lengths of such X intervals where c(X) is constant are the mag-

P(~)

!! X(~)

177

nitudes of the corresponding &peaks in the pdf P(6) distribution. However, where X(~) has the usual (that is, not in a sense of generalised functions) derivatives, so does c(X):

ax(e) _ (oc(X) 7i

)-*

t---T£-

O2X(c) D

'

[ OX(c)

I 3

/

/

02c(X)

c(X)=~" (9)

Thus the initial pdf P(t?) and new reference field c(X) are related by P ( & t)

OX(~,t) _ ( .Oc(X,t) _ _ O~ c~X c(x,t)=e ac(X,t)

if

0X

~6(~ - c(X,t)),

1

4=O, if

Oc(X,t) #X

O, (10)

where a is a length of the X range over which

c(X) = ~. By noting that

ac(X)

ac(X)

ax(~)

3t

OX

Ot

one can obtain a balance equation for c(X, t)

!:~ ' ~

Fig. 3. Sample of the transformation of the nearly bimodal pdf P(~) into the reference scalar field c( X ).

178

A . A . BURLUKA ET AL. pothesis on Xe. Following [9] the conditional scalar dissipation may be written as a sum of two terms: due to flamelets (from within the reacting layers) and "out of flamelet" (due to the large-scale effects) contributions gl. The latter may be suggested as the first approximation:

from Eq. 8 using Eq. 9:

c~c(X,t)ot

Xc (Oc(X't)) -2c9X "

O2c(X,t) OX2 + W(c(X, t)) c~Xc (Oc(X,t))-' + a--X"

aX

,

1

(11)

where the last term is caused by the nonuniform distributions of conditional scalar dissipation Xe. Then by using the relationship

0

c~c(X) I

-

2]

toc 1-2 t---22-- j

oc(x) 1-3 o2c(X)

-2"Xc. ------~-~]

ogX 2

one can rewrite Eq. 11 as 1

Oc(X,t)at = -21OXc(°c(X't) ) OX "" OX 1 Oc(X,t) + W(c(X, t)) + - . 2 8X "OX Xc" ~

•

= - - .A(c(X), (c),

where A(c(X), (c), (c'2)) is a function of order unity. In fact, Xl may depend on many other parameters, for example, on the scales of the velocity field or higher moments of the concentration fields, etc. Some of the estimations below are valid for an arbitrary shape of A(c(X),(C),(C'2)) which, however, for the sake of modelling completeness, is taken in the simplest "relaxation-like" form. Here the physical ideas behind the model, which were presented previously, are reformulated in a statistical sense in order to analyse more rigorously the contribution to Xc from within the flamelets. Consider the ensemble of realisations of the stochastic scalar field c(x, t) in physical space over the integral scale 1. Then according to the definition of pdf P(~) and its "geometrical" interpretation [6, 8], one has

=

_

(oc(X)

ccx : ) l

(12) =

At this point a closure assumption for needed for Eq. 12.

(c'2)),

Z,c

Xc

lira

Ac--,0

is

CLOSURE ASSUMPTION FOR SCALAR DISSIPATION As mentioned above, we are interested in the case when the chemical reactions occur in interacting thin layers, the internal structure of which is determined with chemistry and molecular diffusivity. However, the interaction of those layers is due to the large-scale turbulence and, hence, it should have the time scale almost the same as the integral scale of turbulence rtc. These two phenomena should be taken into account by means of a closure hy-

=

dc(x)

AC

]-1),

(1(_.__~ c,,,=e]

(13)

where the angle brackets stand for ensemble averaging and (Ax(& Ac))/l is the relative part of the x axis (see Fig. 2) occupied by the scalar values in a (~, 6 + Ac) interval; see Fig. 2. In Eq. 13 the existence of a nonzero limit limac _~0[(Ax(O, Ac))/Ac] is implicitly assumed (see details of the procedure in [6]). From the formal properties of the solution of Eq. 4 (Green's function is zero nowhere), this limit is always nonzero. In reality, such a procedure corresponds to consideration of a turbulent flame brush in which the flamelets are not

TURBULENT PREMIXED COMBUSTION infinitely thin, though their thickness is much smaller than turbulence integral length scale l, and there is always a nonzero probability to meet, at a given point, the flamelet containing intermediate concentrations, even if it is very small• By this geometrical analogy the contribution from within the flamelets to the conditional scalar dissipation )6, can be determined as

179 By substituting these relationships into Eq. 12 one obtains

3c(X,t)at - [~-~ - 2 • A(c(X))r,,, (3c(X,t))

2]

3X

~2 c 3X 2 !

+ W(c(X, t)) + --Ttc

D M ( 3 c ( X ) ]2 --- 1-5-.

~-X---] .

Thus, by summarising all the considerations above, a formula for the conditional scalar dissipation becomes

DM ( O c ( X , t ) ) 2 Xc = 12 " OX 1

+ - - . A ( c ( X ) , (c), (c'2)).

(15)

r,c

A possible normalising factor, which would be analogous to /3 in Eqs. 1-5, is omitted here and therefore stretch and curvature effects are neglected. Now with Eq. 15 we can obtain the auxiliary relationships 3X,, - -

3X

DM =

2 • - -

l2

Oc(X,t) 3X

OA(c(X)) 3c

(14)

(17)

If the thickness of reacting zones 6. and their velocity u. are close to those of a laminar flame, one can write D i to be of the order of 6. • u.. Then, by using the relationship

1 3c(X,t) l 3X

~ [Oc(x,t)} \ Ox

and by estimating the mean scalar gradient as ~, - 1 (that also results from geometrical pdf treatment), the term in the square brackets in Eq. 17 may be rewritten as

DM A(c(X))(3c(X,t)) l--5- - 2 " rtc OX

2

,21,2)

O2c(X,t) 3X 2

l---2-" 1 - 2D---~ " r t ~ - ' - - 2I 1

OA(c(X))

+ --

rtc

3c

3c(X,

t)

-- i f - 1-2"u--~" 7

and O -0X -

u)

(16)

3X

[

Xc"

----T-" 1 - 2 "

3X

1 rtc

3A(c(X)) 3c 2

•

- R e -1

(18)

U t7

A(c(X)) --

u2) 2

rtc

a2c(X,t) OX z

(3c(X,t))-' 3X (Üc(X,t))

3

3X (16')

where u' is the r.m.s, turbulent velocity, rtc is estimated as I. u'-l, and Re = (u'. l ) / O M is the turbulent Reynolds number• For the typical values u ' = 1 m / s , l = 1 cm, u~-~ 25 cm/s, and D M ~-0.2 cm2/s, the ratio 2. (u'Z/u,2) • Re 1 = 0.06 is small compared with unity and can be neglected• It is worth noting, that in [16] a similar parameter H = (u'/u,).

180

A . A . BURLUKA ET AL.

Re -1/2, which is the ratio of chemical time to the mean eddy lifetime, appeared from the experimental data analysis of flame quenching. In terms of that H parameter the premixed turbulent flame extinction occurs in experiments when

Finally the rsf equation for the homogeneous case may be written Oc(X,t)

DM

02c(X,t)

Ot

12

OX 2

+ W ( c ( X , t)) -

c ( X , t ) - (c)

Ut2

(19)

2H2 = 2. ---5- • Re 1 = 1 . Note that results of Eq. 19 are:

Un

The consequence of this for the rsf Eq. 17 is that it has a positive diffusivity coefficient in MSM description for any regime of premixed combustion that has been experimentally observed. Though Eq. 18 holds for arbitrary function A ( c ( X ) , (c), (c'2)) of the order of unity, the rsf Eq. 17 requires it to be applied for modelling. Now, to find the possible form of the A ( c ( X ) ) function standing for the flamelets interaction, we apply one of the results in kinetic theory [17]. Namely, in [17] it was shown that the contribution to the dissipation function from any &correlated disturbances is proportional to the square of the deviation from mean divided by the relaxation time. (This may be obtained by developing of Xc into Taylor series and then by the usual estimation of turbulent gradient fluctuations [17].) For the present case it gives A ( c ( X ) , (c)) -- O - ½(c(X,t) - ( c ) ) 2, where ® is some constant level (an estimation may be obtained X < O < 2X) and 1 Ttc

3A(c(X)) 0C

rtc

c ( X , t) - ( c ) Ttc

The last formula appears to be similar to the well-known IEM-LMSE model [1, 3]. In other words, the IEM-LMSE model may be treated as a description of turbulent combustion in the case of vanishing molecular diffusivity, where turbulence induces the interaction between the infinitely thin flamelets.

(a) the computed pdf is always normalised; (b) equations for ( c ) and (c '2) are consistent, namely,

d(c) -

-

dt

=

(w),

d ( c '2 )

(x> + (W'c'>;

(20)

dt (c) the chemical terms are always closed; (d) the Cauchy problem is well-posed for microscale effects description; (e) easy implementation in modelling. Below some results obtained with Eq. 19 are presented. RESULTS FOR THE HOMOGENEOUS CASE

As the first test problem, the MSM processes with and without chemical reaction in a homogeneous gaseous medium are considered with initially bimodal pdf. To obtain the discrete analogy of Eq. 19 a central difference numerical scheme is used and the source terms are linearised. The resulting algebraic system is solved with the three-diagonals matrix algorithm. In Fig. 4 the results are presented for the pure mixing case with ( c ) = 0.5. One can see the symmetric movement of the 6-peaks in the pdf temporal evolution with the progressive decrease of their heights and production of the concentrations in between. In a sense this picture corresponds to the LMSE-IEM model modified to allow for molecular diffusivity, which diminishes the 6-peaks heights and produces concentrations between them. This nonreactive case is, however, beyond the domain for which this model is established and these

TURBULENT PREMIXED COMBUSTION

181

,,' i-] 1.11• ,",/ .~;''"~" 0.1t.

|

i 1.It

"¢'%"~;" .,,"

O.TO-

~

,..-

u

~ . . , ,"° "° O.OS.-.

i

-" ~ ~

,/

"1 1.1

/

P(c) ~,~

/

~.,

GO~,-

Fig. 4. Probability density function (pdf) temporal evolution and corresponding behaviour of second- and fourth-order moments of the scalar field for the pure mixing case with ( c ) = 0.5. Dimensionless parameter D M • " r z c / l 2 = 0.02.

taneous reaction rate taken as in [19]: results are presented only for the sake of completeness. Predicted scalar variance decay has a form exp(-at/rtc), where a monotonically decreases from the initial value a -- 2 to rather small values a ,~ 1, while the scalar flatness factor F = ( c ' 4 ) / ( c ' 2 ) 2 grows steadily from the initial unity level. Such a behaviour of scalar variance and flatness factor qualitatively follow the DNS results [19, 20]. In Fig. 5 the results for reacting case with the initially bimodal pdf with (c)0 = 0.05 are presented. They are obtained with the instan-

c 5 • (1 - c) W(c)

rr

where z r is a characteristic chemical time. One can see that the evolution of the pdf preserves a nearly bimodal shape with low intermediate concentrations probability (compare with Fig. 3). That is, for a remarkable time (near 4~-t¢ in the case at hand, w h e r e DMT"t,./I 2 = 0.02 and r t c / r r = 40.0) there is no motion of &peaks at = 0 and ~ = 1 toward each other. The chemical reaction rate is large enough at ~ ~ 1 to

182

A . A . B U R L U K A ET AL.

7~ P(c)

3o

% ;0

Fig. 5. Probability density function (pdf) temporal evolution and corresponding behaviour of second- and fourth-order moments of the scalar field reacting with (c)0 = 0.05. Dimensionless parameters: DM • r t c / l 2 = 0.02 and r t c / r r = 40.0.

prevent effectively &peak at 6 = 1 from moving and becoming smeared. At the same time, model 19 shows that during the time t = 4rtc the motion of the 6(~) peak is negligible until the heights of peaks 6(~) and c5(6- 1) become nearly the same, and (c '2) attains the maximum (C'2)max = 0.2. It is worth noting that strictly bimodal pdf P(6) = (1 - (c))6(~) + (c)8(~ - 1) results in (c'2)max = 0.25; hence, Eq. 19 gives rise to a noticeable probability of intermediate concentrations. Scalar variance (c 'z) decreases during some small time until the mixing produces sufficient intermediate concentrations to initiate the

chemical reaction. Then (C '2) rises due to the chemistry and then decays. For the reacting case the flatness factor F = (c'4)/(c'2) 2 reveals a complex behaviour. It drops to unity until (c 'z) achieves the maximum; then it rests at unity for one turbulence time, after which it rises rapidly attaining big values, and at the last stage, it drops again. The trends predicted are quite consistent with Eq. 20 and DNS results [19], though the latter are available only for low Damk~Shler numbers. Thus in Fig. 6 temporal evolution of scalar variance (c '2) is presented in comparison with the DNS data [19] and that given by the IEM-LMSE model

TURBULENT PREMIXED COMBUSTION

183


012"

o = 0.048



0~ ~

Da

-. '-.

=

20.

Da=2.

0•010

010

LMSE

005

RSF

RSF

0

°%0

015

o ,~

1'o

1'$

t1%o ° = 0.018 Da = 4.0

o.1o-

000

oo6 oo.

/'

,,,/

'"LMSE

,, 0

t/,l~. Fig. 6. Damk6hler number effects on the temporal evolution of the scalar variance. Dotted lines correspond to the LMSE model results, solid lines correspond to the rsf model Eq. 19, and symbols are DNS results taken from [19].

for different Damk6hler numbers. The initial pdf distribution is the clipped Gaussian, and corresponding values are shown in the legend. One can see the substantial improvement given by taking account of molecular transfer by Eq. 19.

Fig. 7. Typical inlet parameters are: uniform mean flow velocity U0 = 5.0 m / s ; initial level of turbulence u ( = 0.25 m / s ; integral length scale of turbulence l -~ 4 ram; the cold flow temperature T~ = 300 K.

RESULTS F O R V-FLAME To illustrate the abilities of the proposed model a stationary V-flame with low turbulence level is simulated in a configuration corresponding to the experiments [21, 22]. Lean h y d r o g e n - a i r (7% of hydrogen) and m e t h a n e - a i r (5.5% of methane) mixtures were used. The flame stabilised on a catalytic wire propagates upward in freely expanding turbulent flow. The overall combustion flow configuration is sketched in

With these parameters one might expect a considerable influence of chemical kinetics on the global flame configuration [3]; in particular, it can be seen from the estimations of the ratios u ' / u , -- 2.5 and 71/6, = ( u ~ / u ' ) R e I/4 -- 1, where -q is Kolmogorov's microscale of turbulence. Besides, the significance of molecular transfer of mass may be much grater than that of molecular friction, as in a lean hydrog e n - a i r flame the molecular diffusivity D M exceeds more than twice the viscosity u.

184

A. A. B U R L U K A ET AL.

i~ ~ if ~ r~ °

applying a boundary layer approximation for stationary V-shaped flame the rsf c(X; x,y) equation may finally be written as

instantaneous n

~i .~ "&; 1~ ....

t 2m . _e-. ; nes with important ~ flamelet interaction

V,y

(0,0

+

c9C

0 D 7- 3c --o--'-+ W(c) 3y Sc r 3y

fresh mixture

U,x ~

c~c

+ °(c)D'-F " a% OX 2

c - (c) p(c) . - -rtc, (22)

catalyticwire

Fig. 7. Overall configurationof the V-shaped flame. The rsf c(X; ~, t) transport equation for the inhomogeneous case is

Op(c)c(X; ~, t)

+

where p denotes the mean density and Sc r is the turbulent Schmidt number. The instantaneous chemical reaction rate W(c) is taken to be

ap(c)(ui)cc(X; ~, t)

at

W(c) = B p ( c ) . c ( 1 - c ) . exp

OX i

1 + ~bc '

DM 0 2 c ( X ; x , t ) = p(c) 12

p(C)

o3X2

c(X; Z, t) - (c)

+ W ( c ( x ; Z, t)),

zt¢ (21) where (Ui)c is the ith component of flow velocity, conditionally averaged at the scalar value c(X; ~, t). In Eq. 21 the molecular diffusion is neglected in physical 2' space, though it is emphasized to be important in probability X space for MSM modelling. Here, as usual, (u~)c is expressed using the eddy diffusivity concept [8] as

p(c)(Ui)c . f = ( p ) ( u i ) "f + ( p ) Dr Scf

Of Oxi '

where Scf is the turbulent Schmidt number for the f variable. This is adopted here as a first attempt, though there is some evidence that the expression above for (u i)c should be modified for reacting turbulent flow [3]. Then by

where B and /3~ = Ta/T o are the pre-exponential and reduced activation temperatures in the Arrhenius law and ~b = Po/Pb -- 1 is the volume expansion ratio. The values of parameters for hydrogen-air mixture are D M = 0.6 cmZ/s, /3A = 40.0, and ~b = 2.8; for methane-air mixture are D M = 0.2 cm2/s, /3A = 50.0, and ~b = 4.2. The pre-exponential factor is taken to be B = 8.1.109 s -1 for both mixtures. These constants approximate reasonably well the experimental conditions. In particular, they give a laminar burning velocity u, = 12.0 c m / s for methane- and u, = 10.0 c m / s for hydrogenair mixture at these low given equivalence ratios. The standard variant of the "k - e " model [23] of turbulence is used. The turbulent time scale for scalar field and that given by the k - e model were taken to be equal: z,c = k/6. The length scale l in Eq. 22 is taken to be equal to that of the standard k - e model [23]. On the free boundary the parameters of the cold fresh mixture are used. The zero gradient conditions on the symmetry axis were imposed for all variables.

TURBULENT PREMIXED COMBUSTION

The governing set of equations is solved with the G E N M I X code [24] enhanced for rsf calculations. Note that the rsf equation with an extra probabilistic dimension is solved with the usual finite-difference technique. In Fig. 8 radial distribution of mean temperature and methane mass fraction given by the rsf model are compared with the measurements [22] for two distances downstream at x = 40.0 mm and x = 80.0 mm. They are calculated from the computed rsf as

c(X; x,y) c(X; x, y)

c(X; x, y)

(T(x,y))= To.(l + ChfolC(X;x,y)dX), (CcH,) = C °cn4

"(1-folc(x;x,y)dX ) •

They demonstrate reasonable agreement for both profiles. The results for the H z - a i r flame are compared with experiments in Fig. 9. There are two general remarks which are applicable equally well to the modelling of both flames. First, as was emphasized in [21], the standard k - c model requires that additional terms be introduced to take into account

CH,

185

variable density and pressure field interactions. (For example, even if there is no mean pressure gradient in the cold flow, the recirculation zone behind the flame holder may induce it.) The importance of "pressure-velocity dilatation" effects is stressed in [6] as well. Further discussion of these questions as well as more detailed comparisons of predictions obtained with Eq. 22 and the results from other models may be found in [22], including mean velocity, turbulent kinetic energy, and its dissipation in several cross sections downstream. Second, as was found in [3], the turbulent Schmidt number may depend strongly on the local turbulence characteristics as well as on combustion. For the flame considered it may vary by a factor of 3 [21]. However, it can be seen from Fig. 9, that the agreement of the model predictions with experiments is not as good for the H 2 - a i r flame as it is for CHa-air. Both of these flames are rather similar, i.e., they have nearly the same values of the dimensionless numbers, except Lewis (Le) number. While for the C H a - a i r flame Le is very close to unity, this lean H 2 - a i r mixture has a Lewis number L e - - - 3 , which results in the temperature and concentration

TOK J 1500

%] mm

,•

•

X : 0.04 rn 34

'

A

X = 80 m m 23

I-.°',-i

points points

X= .4000E-01 ........

~

X= .8000E-01

mm

poonls

•

X = 0 . 0 4 m 35

l~

X = 0.00 m 20 Oo~ts X= .4000E-01

.......

X= ,8000E-01

;

A

1000

1

0.000

=w

4',

50(1

I

0.010

.

.

.

.

I

0.020

'

Y,m

0.1 O0

)

A

\ ' 0.~10

r

0.020

Y, m

Fig. 8. C o m p u t e d r a d i a l profiles (lines) of the m e a n m e t h a n e v o l u m e fraction a n d t e m p e r a t u r e versus e x p e r i m e n t a l d a t a (symbols) for two d i s t a n c e s x d o w n s t r e a m : • x = 4.0 cm; zx x = 8.0 cm.

186

A . A . BURLUKA ET AL.

/ ~ . - . . ....

.... A •=. "

........

H,%

~/,"

A

~

o°°

6

A

T'K

AA'

1000

•

i 'i i

A A

A:

/ o/ o -: _ - . , x o •

;

Zk

:

A

iI

m

X = 0 . 0 8 m 17 points

~

x= .4oooE-ol

........

X= i 81~0E-01

800

o' ~A _

•

•

,',,

i

•

-

•

X = 0.04 m 29 points

A

X = 0,08 m 26 points

-

X= .4000E-01

x..SO0OE-Ol

........ A

', A i t

600

i, !

\.

A

•

•

400 i

)

•

A

:

A

'..

%

~mmm ",,

~kAA . . . . . • AA#~X~x

u~,...~.. 0

,

,¢

0.000

'

I

0.005

'

'

'

'

I

0.010

.

.

.

.

'

I

0.015 Y, m

0.~0

I

0.005

'

'

'

~

I

0.010

.

.

.

. .

I

0.015 Y, m

Fig. 9. Computed radial profiles (lines) of the mean hydrogen volume fraction and temperature versus experimental data (symbols) for two distances x downstream: • x = 4.0 cm; z~ x = 8.0 cm.

diffusivity and flamelet interaction enter the model explicitly. The computations performed for the homogeneous case show the rsf model to retain the bimodal shape of the pdf for the reacting case. Even when it is coupled with an approximate turbulence model, the approach proposed results in a reasonable agreement with experiments in simulations of V-shaped flames. Though at present turbulent diffusion is treated with a single diffusivity coefficient, a more rigorous turbulence description may be included in the rsf model.

fields being dissimilar, that, in turn, casts doubt even on the possibility of introducing a single progress variable c. In this case the details of flame stabilisation and chemical kinetics may be important. Thus, all the visible differences in behaviour of these flames should be attributed to Lewis number effects [18]. This hypothesis requires further experimental verification which should clear up how the molecular transfer (microscale effects) may change the overall macroscale combustion flow structure. However, taking into account the high complexity of the phenomena, one may estimate the prediction of the present model as reasonable, in spite of a quite approximate parabolic model of the flow dynamics.

The authors are grateful to Prof. C. Dopazo, Prof. V. A. Sabel'nikov, and Prof. V. A. Frost for valuable discussions.

CONCLUSION

REFERENCES

A statistical description of microscale phenomena in turbulent premixed flames is proposed in terms of a reference scalar field. Its transport equation is derived from that for a pdf. A closure assumption on the conditional scalar dissipation is made for the flamelet combustion regime. From this assumption, molecular

1.

Frost, V. A., Transactions of the Third All-Union Conference on Combustion, Izd-vo AN SSSR, Moscow, 1960, pp. 121-125 (in Russian). 2. O'Brien, E. E., in Turbulent Reacting Flows (P. A. Libby and F. A. Williams, Eds.), Springer-Verlag, Heidelberg, 1980, vol. 44, pp. 185-218. 3. Borghi, R., Prog. Energy Combust. Sci. 14:245-292 (1988). 4. Dopazo, C., in Turbulent Reacting Flows (P. A. Libby

TURBULENT PREMIXED COMBUSTION

5. 6. 7. 8. 9. 10.

11. 12.

13. 14. 15.

and F. A. Williams, Eds.). Academic Press, London, 1994, pp. 375-474. Kerstein, A., Combust. Sci. Technol. 60:391-421 (1988). Kuznetsov, V. R., and Sabel'nikov, V. A., Turbulence and Combustion. Hemisphere, New York, 1990. Zimont, V. L., private communication. Tennekes, H., and Lumley, J. L., A First Course in Turbulence. MIT Press, London, 1972. Borghi, R., Combust. Flame 80:304-312 (1990). Emelianov, V. M., Frost, V. A., and Nedoroob, S. A., Twenty Fourth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1992, pp. 435 -442. Moss, J. B., Combust. Sci. Technol. 22:119-129 (1980). Klimenko, A. Yu., and Bilger, R. W., Relationship between Consented Scalar PDFS and Scalar Dissipation in Turbulent Flows, Charles Colling Research Laboratory, Report No. TN F-100, 1993. Kraichnan, R. H., Bull. Amer. Phys. Soc. 34:2298 (1989). Curl, R. L., AIChEJ. 9:175-181 (1963). Yoshida, A., and Tsuji, H., Seventeenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1978, p. 945.

187 16. Abdel-Gayed, R. G., and Bradley, D., Combust. Flame 62:61-68 (1985). 17. Landau, L. M., and Lifshitz, E. M., in Course of Theoretical Physics. Statistical Physics, vol. 9 (E. M. Lifshitz and L. P. Pitaevskii, Eds.), Pergamon Press, Oxford, 1980. 18. Karpov, K. P., and Severin, E. S., Comb~t. Explos. Shock Waves 16:45-51 (1980). 19. Picart, A., Borghi, R., and Chollet, J. P., Sixth Symposium on Turbulent Shear Flows, Toulouse, Sept. 1987. 20. Pope, S. B., and Eswaran, V., Phys. Fluids 31:506-520 (1988). 21. Borghi, R., and Escudie, D., Combust. Flame 56:149-164 (1988). 22. Burluka, A. A., and Borghi, R., in Proceedings of First ETCM Workshop, Ecole Centrale de Lyon, Aussois, 5-10. February 1995. 23. Launder, B. E., and Spalding, D. B., Mathematical Models of Turbulence. Academic Press, London, 1972. 24. Spalding, D. B., A General Computer Program for Two-Dimensional Parabolic Phenomena. Pergamon Press, London, 1978, p. 378. Received 18 May 1995