- Email: [email protected]
Modification of the turbulent burning velocity by gas expansion
Proceedings of the Combustion Institute, Volume 28, 2000/pp. 235–243
MODIFICATION OF THE TURBULENT BURNING VELOCITY BY GAS EXPANSION N. PETERS,1 H. WENZEL1 and F. A. WILLIAMS2 1
Institut fu¨r Technische Mechanik RWTH Aachen, D-52056 Aachen, Germany 2 Center for Energy and Combustion Research Department of Mechanical and Aerospace Engineering University of California at San Diego La Jolla, California 92093-041, USA
Premixed turbulent combustion in the corrugated flamelets regime has been addressed by modeling of a G-equation containing a gas-expansion term, the form of which was derived by generalizing results of Sivashinsky along lines initiated by Frankel. Previous analyses by the first author of the mean, variance, and two-point correlation spectrum of G, as well as their relation to the turbulent burning velocity, were extended to include the influences of gas expansion. It is shown that this effect has the same form as kinematic restoration but acts in the opposite direction. Because of this similarity in form, no new modeling approximations are required. On average, kinematic restoration dominates gas expansion. In physical terms, this means that the energy-containing eddies at the integral scale are strong enough to suppress gasdynamic instabilities at smaller scales. Gas expansion lessens the smoothing effect of kinematic restoration, thereby increasing the variance of G and the turbulent flame brush thickness. It also increases the turbulent burning velocity by a factor that increases with increasing ratio of unburned to burned gas density and with decreasing ratio of turbulence intensity to laminar burning velocity. The effect on the burning velocity, however, is not large, remaining less than 40% even as the turbulence intensity begins to decrease below the laminar burning velocity.
Introduction
sT ⳱ b1v⬘
In recent years, several attempts have been made to derive models for premixed turbulent combustion based on a level-set approach using the so-called Gequation [1,2],
where b1 is a constant of proportionality which has been determined from experimental data as b1 ⳱ 2.0. More recently [4], the level-set approach was used to describe a new regime in turbulent premixed combustion called the thin reaction zones regime, where the smallest turbulent structures may enter the preheat zone of the flame but not the reaction zone. In the thin reaction zones regime, the diffusive effects become dominant. In equation 1, gas-expansion effects enter through changes of the velocity field v induced by heat release at the flame surface. The impact of gas expansion can be analyzed from jump conditions for the flow variables across the flame surface. Numerical methods using local reconstruction within a numerical cell have been derived [5,6] to do this accurately. The objective of the present paper is to derive closure relations as in Ref. [3] but with the effect of gas expansion included. For that purpose, we start from a modified version of the G-equation. The foundation for this modification was laid by Sivashinsky [7]. The original Sivashinsky equation has been subsequently analyzed and extended in a number of papers (e.g., [8–13]).
G Ⳮ v • G ⳱ sL|G| t
(1)
The G-field is convected by the flow velocity v and propagates into the unburned mixture with the laminar burning velocity sL. An isosurface G(x, t) ⳱ G0, where G0 is a fixed constant, describes the location of the flame front in an arbitrary flow field. The levelset approach was used in Ref. [3] to derive a twopoint closure relation in a regime of premixed turbulent combustion that we call here the corrugated flamelets regime. It has also been called the multiple-sheet or fractal regime in the literature. That regime is characterized by a local kinematic balance between flow velocity and burning velocity. It was also shown in Ref. [3] that the turbulent burning velocity sT becomes independent of the laminar burning velocity sL and is proportional to the turbulence intensity v⬘ in the limit of large values of v⬘/ sL,
235
(2)
236
TURBULENT COMBUSTION—Premixed Flame Modeling
To facilitate analysis, we adopted a model that does not account for the all physical effects that occur in a premixed turbulent flame. Gas-expansion effects are treated as volume sources located on the flame front. The volume sources generate a potential flow field which then is simply superposed on a turbulent flow field that is described by the NavierStokes equations for constant density flows. The resulting velocity field does not satisfy the non-constant density Navier-Stokes equations and should be considered as a model. This model does not account, for instance, for the vorticity generated by the flame. While this exerts influences, it may be expected that it will not qualitatively affect the results. Recently, Bychkov [14] has shown that inclusion of vorticity generation in his model equation for curved stationary laminar flames reduces the burning velocity by 40% for an unburned-to-burned density ratio qu/qb of 5 but does not introduce important qualitative modifications. The modeling of gas expansion by a source term in the G-equation has the advantage that its effect can be easily analyzed and quantified, a task that is not trivial if one analyses data of direct numerical simulations with simplified or full chemistry and heat release [15]. The modeling of the additional source term in the G-equation is addressed first. A numerical analysis is then presented in which a volume source at each point of the flame surface, previously used in a numerical study by Ashurst [16], is employed to circumvent the restriction of small departures of the flame sheet from a plane. Frankel [11] has shown for the two-dimensional case that in the limit of small flame front perturbations, the use of volume sources is equivalent to the non-local expansion term appearing in the Sivashinsky equation [7]. Finally, the main results are summarized and interpreted in terms of the turbulent burning velocity.
mixture immediately ahead of the flame structure, sL ⳱ sLu and v ⳱ vu in equation 4. If, on the other side, the flame front position is defined to lie in the burned gas immediately behind the flame structure, sL ⳱ sLb and v ⳱ vb in equation 4. Introducing equation 4 into equation 3 gives, for the velocity jump normal to the flame front, the expression c s ⳱ csLb 1 ⳮ c Lu
n • (vu ⳮ vb) ⳱
where c ⳱ (qu ⳮ qb)/qu is the gas-expansion parameter, which is zero for constant density but close to unity for large values of qu/qb. Equation 5 shows that gas-expansion affects scale with the laminar burning velocity in the sense that the velocity jump is proportional to sL and suggests that they will be important only if sL is of the same order as the turbulent velocity fluctuations v⬘ or larger. These effects are important in the corrugated flamelets regime but are likely to be insignificant in the thin reaction zones regime, which therefore will not be considered. In the work of Sivashinsky [7], an equation for small perturbations of a planar two-dimensional flame front in the limit of small c was derived, accounting for stabilizing or destabilizing diffusionalthermal instabilities associated with negative Markstein length and the Darrieus-Landau instability arising from gas expansion. For simplicity, we do not address the former here, setting the Markstein length Lb small but positive. The latter appears in the equation for the perturbation F(y, z, t) of the flame front location as an integral over transverse coordinates. We introduce G ⳱ x ⳮ F(y, z, t) and write Sivashinsky’s integral as
1 I(G,x,t) ⳱ (2p)3
冮冮冮冮冮冮 ⳮ
ik•(xⳮx˜)
|k| e Formulation and Modeling In the thin flame approximation, the the mass balance across the flame is
冢
qun • vu ⳮ
冣
冢
冣
dxf dx ⳱ qbn • vb ⳮ f dt dt
(3)
where qu and qb, and vu and vb are the densities and velocities of the unburned and burned gases, respectively, n is the local normal vector of the flame front pointing toward the unburned mixture, and dxf ⳱ v Ⳮ sLn dt
(4)
is the propagation velocity of the front into the unburned mixture with respect to the unburned gas, with sL being the laminar burning velocity. If the position of the front is defined to lie in the unburned
(5)
G(x˜,t) dx˜dk
(6)
As a generalization of Sivashinsky’s equation, the resulting G-equation on which we base our modeling is G Ⳮ v • G ⳱ s0Lb|G| ⳮ DLbj|G| t Ⳮ
c 0 s I(G,x,t) 2 Lb
(7)
where s0Lb is the burning velocity of an unstretched laminar flame with respect to the burned gas, 0 DLb ⳱ sLb Lb is the corresponding Markstein diffusivity, and j is the curvature. In a numerical study, Dandekar and Collins [17] used a similar formulation. When the assumption of small perturbations about a plane is removed, the integral becomes a volume integral [11]. Further details of the derivation may be found in Ref. [18].
TURBULENT BURNING VELOCITY
¯ , the To develop equations for the mean value of G variance G⬘2, and the two-point correlation function g2(r, t) ⳱ G⬘(x, t)G⬘(x Ⳮ r, t), we regard equation 7 as a three-dimensional field equation in a constantdensity flow field and split the variables G, r ⳱ |G| ¯ and v into mean and fluctuating parts [3], G ⳱ G Ⳮ G⬘, r ⳱ r¯ Ⳮ r⬘, v ⳱ v¯ Ⳮ v⬘. The equation for ¯ is readily obtained as G ¯ G ¯ Ⳮ • v⬘G⬘ ⳱ Ⳮ v¯ • G t c ¯ I(G,x,t) (8) 2 ¯ in the argument of I(G ¯ , x, t) is where the use of G due to the linearity of the expansion integral I(G, x, t) with respect to G. This linearity is an important feature and is exploited in analyzing the variance equation and the equation for the two-point correlation. A surprising result of equation 8 is that for a freely propagating planar turbulent flame, the expansion ¯ , x, t), being an odd function of the coorterm I(G dinate with respect to the flame in the direction of propagation, vanishes and has no direct influence on the propagation speed of that flame. The expansion ¯ only by modifying effects enter the equation for G the flame surface area ratio r. ¯ In this formulation, the variance G⬘2 of G is related to the turbulent flame brush thickness lF,t (cf. equation 30 below). The equation for its dynamic evolution is obtained by subtracting the equation for ¯ (equation 8) from the instantaneous G-equation G (equation 7) to get the equation for the fluctuation of G. Multiplying this by G⬘ and then averaging yields the equation for the variance of G, namely,
p¯ ⳱
s0Lb (2p)3
冮冮冮冮冮冮|k| e
ik•(xⳮx˜)
G⬘(x˜,t)G⬘(x,t)dx˜dk
(9)
(10)
Numerical simulations [19,20] have shown that these terms always act as sinks for the variance. The physical interpretation of kinematic restoration is that wrinkles introduced by turbulent velocity fluctuations are removed by the smoothing effect of local laminar flame propagation, tending to restore a flatter surface. Scalar dissipation is a higher-order
(11)
ⳮ
Like kinematic restoration, this quantity is proportional to the laminar burning velocity, but differently from x, ¯ it acts as a source term for the variance and increases irregularities of the flame front similar to the way that the Darieus-Landau instability increases initial small disturbances of a planar flame in a quiescent medium. At this stage, there is no information on the relative orders of magnitude of these two counteracting terms and consequently no forecast can be made as to which will dominate in a turbulent premixed flame. Some information relevant to this can be gained by proceeding with the formulation of the conservation equation for the two-point correlation function along the lines described previously [3]. Following standard methods [21] for homogeneous isotropic turbulence, all spatial dependencies are reduced to the correlation distance r ⳱ |r|, and the correlation function obeys g2(r,t) t Ⳮ 2 • v⬘(x Ⳮ r,t)G⬘(x,t)G⬘(x Ⳮ r,t) 0 Ⳮ 2s0LbS1 Ⳮ 2DLbS2 ⳮ csLb S3 ⳱ 0
The terms on the left-hand side are local accumulation of variance, transport by the mean velocity, and turbulent transport. The first term on the righthand side is the production of variance by the action of velocity fluctuations on the mean gradient of G. The next terms on the right-hand side are the kinematic restoration and the scalar dissipation [3], defined here as x ¯ ⳱ ⳮ2s0LbG⬘r⬘, v¯ ⳱ 2DLb(jr)⬘G⬘
term in the corrugated flamelets regime as shown in Ref. [4]. It becomes more important, however, in cases where gas-expansion effects tend to cancel kinematic restoration, as is seen below. The last term on the right-hand side of equation 9, the covariance p¯ of fluctuations of the G-field and the expansion integral I, is defined as
0 s0Lbr¯ ⳮ DLbj|G| Ⳮ sLb
G⬘2 Ⳮ v¯ • G⬘2 Ⳮ • v⬘G⬘2 ⳱ t ¯ ⳮx ⳮ2v⬘G⬘ • G ¯ ⳮ v¯ Ⳮ cp¯
237
(12)
Here the two-point correlations between r⬘ and G⬘ and between (jr)⬘ and G⬘ have been denoted, respectively, by [3] S1(r,t) ⳱ ⳮr⬘(x,t)G⬘(x Ⳮ r,t), S2(r,t) ⳱ ⳮ(j(x,t)r(x,t))⬘G⬘(x Ⳮ r,t)
(13)
The last term of equation 12 is defined as S3(r,t) ⳱ I(G,x,t)G⬘(x Ⳮ r,t)
1 ⳱ (2p)3
冮冮冮冮冮冮 |l| e
il•(xⳮx˜)
ⳮ
G⬘(x˜,t)G⬘(x Ⳮ r,t) dx˜dl
(14)
and describes the influence of gas expansion on the evolution g2(r, t). To close the equation for the two-point correlation function, we follow Ref. [3] and transform it into
238
TURBULENT COMBUSTION—Premixed Flame Modeling
wave-number space k by applying Fourier transforms to obtain the scalar spectrum function C(k, t). With use made of the transform pair
冮冮冮 f(r,t)e
ikr
dr,
ⳮ
f(r,t)⳱
冮冮冮
ⳮikr ˆ f(k,t)e dk
(15)
ⳮ
for any function f(r, t), the spectrum function C(k, t) of fluctuations of G is related to the Fourier transform of g2(r, t) by
冮
C(k,t) ⳱ k2 䡬 gˆ2 (k,t)dX ⳱ 4pk2gˆ2(k,t)
(16)
where integration is carried over the solid angle dX, and k ⳱ |k| is the absolute value of the wave number. Taking the Fourier transform of the divergence of the triple correlation v⬘(x Ⳮ r, t) G⬘ (x, t) G⬘ (x Ⳮ r, t) leads to a spectral transport term T(k, t), which was approximated in Ref. [3] by a gradient transport assumption, thereby making it local in wave-number space. Dimensional arguments and the assumption of locality led in Ref. [3] to linear approximations for the Fourier transforms Sˆ1 and Sˆ2, ⳮ1
2ˆ
ⳮ1 2
2ˆ
8pk S1 ⳱ c1Cs kC, 8pk S2 ⳱ c2Cs k C
(17)
where Cs is the universal constant of the scalar spectrum, and c1 is a modeling constant that is assumed to be independent of the wave number to agree with the locality hypothesis. To derive the Fourier transform of S3, we introduce vector s, defined by x Ⳮ r ⳱ x˜ Ⳮ s, into equation 14 and change the integration variable x˜ to s and l to ⳮk. Since x and r are fixed, dx˜ ⳱ ⳮds and dl ⳱ ⳮdk, and the definition for S3 becomes
S3 ⳱
冮冮冮 ⳮ
C(k,t) 0 ⳮ T(k,t) Ⳮ c1Csⳮ1sLb kC(k,t) t 1 2 0 Ⳮ c2Cⳮ S DLbk C(k,t) ⳮ csLbkC(k,t) ⳱ 0
G⬘(x˜,t)G⬘(x˜ Ⳮ s,t) dsdk
iks
The noteworthy point about this equation is that the term that accounts for gas-expansion effects has the same wave-number dependence as the kinematic restoration. These two terms may therefore be combined to give C(k,t) ⳮ1 0 ⳮ T(k,t) Ⳮ c*C 1 s sLbkC(k,t) t 1 2 Ⳮ c2Cⳮ S DLbk C(k,t) ⳱ 0
G⬘2 ⳱ ⳮx* ⳮ v¯ t
⳱
冮冮冮 |k|gˆ (k,t)e
dk
(22)
is recovered, where the variance, the modified kinematic restoration, and the scalar dissipation have been defined as G⬘2 ⳱
冮
k0
C dk,
ⳮ1 0 x* ⳱ x ¯ ⳮ cp¯ ⳱ c*C 1 s sLb
1 v¯ ⳱ c2Cⳮ s DLb
冮
k0
冮
k0
kC dk,
k2C dk
x* ⳱ c*x
ikr
(21)
with c*1 ⳱ c1 ⳮ cCs. Upon integrating equation 21 over wave-number space, the spectral transfer term vanishes, and the variance equation for homogeneous isotropic turbulence,
ⳮ
2
(20)
(23)
respectively. The lower limit of integration in these integrals is taken to be the inverse of the integral turbulent length scale k0 ⳱ lⳮ1. Provided that c*1 0 for k lⳮ1, the modified kinematic restoration reduces fluctuations of the flame front. It then can be modeled in the same way as x ¯ in Ref. [3], namely, as
|k|eikr (2p)3
冮冮冮 e
(19)
We are now able to write the dynamical equation for the spectrum function C(k, t) by taking the Fourier transform of equation 12 and applying the definition of C(k, t) in equation 16, obtaining
1 ˆ f(k,t) ⳱ (2p)3
C(k,t) Sˆ 3(k,t) ⳱ |k|gˆ2(k,t) ⳱ 4pk
(18)
e G⬘2 k
(24)
ⳮ
Comparing this with the definition of the inverse Fourier transform, equation 15, gives, due to linearity of S3 in g2, with the use of equation 16 the closed form
Numerical Simulation Numerical simulations can check the validity of some of the very strong assumptions that were made in the previous section. Results from two different
TURBULENT BURNING VELOCITY
239
the cusps without adding numerical diffusion to the problem, because no flux limiting is used, and the arising Riemann problem is solved exactly. This procedure has been described elsewhere [20,27]. Fig. 1 shows a sketch of the model used to incorporate gas expansion in the numerical simulation. At every point on the flame surface G ⳱ G0, there is a volume source that induces a velocity into the whole flow field. We are interested in the induced normal velocity vex(x) ⳱ n(x) • vn(x) only at the flame front itself, which can be calculated as vex(x) ⳱ c Fig. 1. A schematic illustration of the velocity induced at the flame front position x by the effect of gas expansion originating from x¯.
numerical simulations are presented, one using the limit c ⳱ 0 to test the model that has led to the first of equation in equation 17 and to calculate the numerical value of c1, the other including the effects of gas expansion for values of c between 0.5 and 0.8 corresponding to density ratios qu/qb between 2.0 and 5.0. It is found that that these effects become less important for larger values of the turbulence intensity v⬘/sL. The model introduced above allows for the turbulent flow field to be computed independently of the G-field. This flow field is calculated in a cubic box using a pseudospectral code developed by Ruetsch [22] which employs forcing of the large scales to allow for a statistically stationary homogeneous isotropic velocity field following a method developed by Pope and Eswaran [23]. Thereby, velocity fluctuations induce fluctuations of G. The forcing of the scalar field is accomplished by imposing a mean gradient in the x direction [22,24], G(x,y,z,t) ⳱ x Ⳮ g(x,y,z,t)
冮
A
(x˜ ⳮ x) • n(x) dA |x˜ ⳮ x|3
(26)
in the three-dimensional case. These integrals of the Biot-Savart type have to be taken over the whole flame front. Although it is not immediately clear that the integrand stays bounded in the limit x˜ → x, it can be shown to converge to a finite value, which turns out to be that of the local curvature [20]. Similarly to Ashurst [16], we include the expansion effects into the G-equation by modifying the laminar burning velocity with the induced velocity,
冢
sL ⳱ sLb 1 ⳮ
冣
c Ⳮ vex(x) 2
(27)
and insert it into the G-equation (equation 1). The factor 1 ⳮ c/2 originates according to equation 4 from the definition of the flame surface G(x, t) ⳱ G0 at the midpoint of the flame structure. A direct evaluation of the integral in equation 26 by triangulation proved to be numerically too expensive. Therefore, a transformation introduced by Maz’ja [28] was used to express the surface integral as a volume integral according to
冮冮 f(r,x) dS ⳱ 冮冮冮 f(r,x)d(G(r) ⳮ G ) S
v
(25)
thereby simulating a planar turbulent flame. Note that this decomposition allows for a multivalued solution in the direction of the mean flame propagation and is therefore more general than the decomposition employing F(y,z,t). Since the solution for g(x, y, z, t) is periodic in all three space dimensions, the same pseudospectral method that is used to calculate the velocity field can be used to compute the scalar field. Equations of the G-equation type (equation 1) describe in the limit of constant sL ⳱ s0L a front propagation of the Huygens type [25]. The front develops cusps which lead to unstable solutions if a standard pseudospectral method is used to calculate the absolute value of the gradient of G on the right-hand side of the G-equation. Following Ref. [26], an upwind method was developed to account for this behavior. This upwind scheme is capable of capturing
sLb 4p
|G(r)|drxdrydrz
0
(28)
where d(G(r) ⳮ G0) is the Dirac delta function, and f(r) is the integrand in equation 26. An equivalent procedure to evaluate surface integrals was used by Vervisch et al. [29] to compute the surface area of premixed turbulent flames.
Results 3
On a 64 grid on a cubic domain of edge length 2p, the evolution of the G-equation was computed first for c ⳱ 0 for three different ratios of the turbulence intensity v⬘ to the laminar burning velocity s0L. The root-mean-square value of the normalized velocity fluctuations was v⬘ ⳱ 17.4, the normalized integral length scale was l ⳱ 1.08, and the normalized kinematic viscosity was m ⳱ 0.24. This results in a Reynolds number of 42 based on the Taylor
240
TURBULENT COMBUSTION—Premixed Flame Modeling
Fig. 2. Evaluation from DNS of the coefficient c1 in the model for the Fourier transform of the two-point correlation representing the kinematic restoration. Dotted line, v⬘/sL ⳱ 1.0; solid line, v⬘/sL ⳱ 1.6; dashed line, v⬘/sL ⳱ 2.0.
microscale. After a statistical stationary stage was reached, the spectrum function and the Fourier transform Sˆ1 of the kinematic restoration was calculated, and the coefficient c1 was then evaluated from equation 17 using CS ⳱ 1.2 for the universal constant of the scalar spectrum. The results, given in Fig. 2, show that this coefficient is larger than unity for k k0. Therefore, c*1 ⳱ c1 ⳮ cCs is positive in the inertial range of the scalar spectrum c ⱕ 0.8, corresponding to density ratios qu/qb ⱕ 5. Even though these results were obtained at c ⳱ 0, they provide a first confirmation of the hypothesis that kinematic restoration dominates gas-expansion effects in most cases of practical interest. Since the numerical evaluation of the non-local expansion integral in equation 28 is expensive, computations for c ⬆ 0 were confined to a 323 grid. Attention was focused first on the behavior of the mean gradient r, ¯ which can be interpreted as the turbulent flame surface area ratio, which in turn is proportional to the ratio of the turbulent to the laminar burning velocity [2,3,19]. Fig. 3 shows the ratio of the turbulent burning velocity thus calculated with and without gas-expansion effects as a function of c and v⬘/sLm with
冢
sLm ⳱ sLb 1 ⳮ
冣
c 2
(29)
The increase of the turbulent burning velocity with increasing c is substantial only for the smallest value of the turbulence intensity shown, where it is still less than 40%; it is only 12.5% for v⬘/sL ⳱ 3.0. To further test the hypothesis that the kinematical restoration dominates the destabilizing effect of the expansion term, the time series of the source terms of the variance equation (equation 9) is plotted in
Fig. 3. DNS results for the ratio of the turbulent burning velocity sT to the value sT1c. Both burning velocities were calculated as the long time averages for flame surface area ratio r, ¯ sT with volume sources in the G-equation included, sT,c ⳱ 0 without volume sources.
Fig. 4. Time series of the source terms in the variance equation. The turbulence intensity v⬘/sLm is 2/3 and c ⳱ 0.8. The time is normalized with the integral time scale and the other terms with sLm and lF,t.
Fig. 4. For numerical reasons, it is necessary to include the curvature term with a positive Markstein diffusivity DLb of 0.48 in the G-equation (equation 7). As can be seen from that figure, the scalar dissipation v, ¯ which accounts for curvature effects, is significantly smaller than the other two terms, which is consistent with the analysis in Refs. [4] and [3], where it was shown that scalar dissipation is a higherorder contribution in the corrugated flamelets regime. Fig. 4 also shows that the kinematical restoration term x ¯ closely tracks the expansion term p, ¯ dominating it most of the time. It is perhaps surprising how often p¯ exceeds x, ¯ leading to a relatively small value c*x ⳱ 0.5 in equation 24 for this example. The Turbulent Flame Brush Thickness and the Turbulent Burning Velocity It can be shown [4,18] that in the limit of a stationary planar premixed flame, all gradients of G⬘2
TURBULENT BURNING VELOCITY
vanish, and the production term in equation 9 may ¯ |2, where Dt ⳱ a4v⬘l is the be modeled as 2Dt|G turbulent diffusivity. Neglecting the scalar dissipation term and using the model in equation 24 for x* with k/e ⳱ a3l/v⬘, one obtains for the turbulent flame brush thickness ᐉF,t ⳱
冢
1/2
冣
(G⬘2)1/2 2a3a4 ⳱ ¯ c*x |G|
ᐉ
(30)
The coefficients a3 ⳱ 4.05 and a4 ⳱ 0.78 have been evaluated in Ref. [18] using standard models of turbulence. For constant density, the coefficient c*x becomes equal to the value cx ⳱ 1.62 obtained in Ref. [3]. Since gas-expansion effects lead to smaller values of c*x, it follows from equation 30 that the turbulent flame brush thickness increases accordingly. The model for the kinematic restoration x ¯ has been used together with dimensional arguments in Refs. [3] and [4] to relate the turbulent burning velocity sT and the flame surface area ratio r¯ for a stationary planar flame in the corrugated flamelets regime to x ¯ as 1/2
冢 ke冣
sT ⳱ sLr¯ x ¯
(31)
¯ is replaced by x* for the case where expansion If x effects are included, it follows with equations 24 and 30 that sT
2 1/2
冢c*a G⬘ᐉ 冣 x
2 3
2
1/2
v⬘
冢2aa 冣 4
v⬘
(32)
3
¯ | ⳱ 1 has been used without loss of genwhere |G erality (cf. [18]). The last expression in this relation shows that in the limit of large values of v⬘/sL, the turbulent burning velocity is independent of c*x and therefore of the gas-expansion effects modeled by it. The modeling constant b1 in equation 2 is therefore expected to be independent of the gas-expansion parameter c. Conclusions The preceding simulation results support the general validity of the model in Ref. [3] that led to equation 2 in the corrugated flamelets regime. Kinematic restoration counteracts and dominates gas-expansion effects on the time average, but their balance may instantaneously be reversed. Kinematic restoration is the basis of the cascade hypothesis for scalar fluctuations; its reduction by gas-expansion effects can become substantial, leading to a turbulent flame thickness appreciably larger than obtained without gas expansion. This is consistent with an earlier simulation [16]. The associated increase in the wrinkled flame area provides, however, only a slightly higher turbulent burning velocity. The simulation results shown in Fig. 3 provide this enhancement in
241
the turbulent burning velocity as the gas-expansion parameter c increases. For the realistic value c ⳱ 0.8, gas expansion increases sT by about 30% for v⬘ ⳱ sLm and only about 12% for v⬘ ⳱ 3sLm. The influence of gas expansion on sT therefore is significantly reduced by increasing the turbulence intensity. This is also shown to be consistent with previous modeling. In high-intensity turbulence, gas expansion therefore can be neglected in models of the turbulent burning velocity. The effect of gas expansion does, however, increase as the turbulence intensity decreases, and in the weak turbulence limit the effect can be dominant, as earlier studies have suggested [30]. A consequence of these findings is that the turbulent burning velocity is a well-defined turbulent quantity that depends on other local mean quantities such as the local turbulence intensity. Non-local effects associated with gas expansion have no influence on scales smaller than the integral length scale. On larger scales, however, where kinematic restoration is unable to control gas-expansion effects, large-scale instabilities may still develop. This may be deduced from Fig. 2 where c1 1 for k lⳮ1, leading to c*1 0 for realistic values of c. Large-scale instabilities of a turbulent front are therefore equivalent to those of the Darieus-Landau type of a laminar front, with the laminar burning velocity being replaced by the turbulent burning velocity, and with a stabilizing effect occurring at the integral length scale. Since only positive values of the Markstein length were used for numerical reasons, it is not clear, whether the destabilizing effect of a negative Markstein length, which would lead to a negative scalar dissipation rate, could unsettle the delicate balance between kinematic restoration and gas-expansion effects.
Nomenclature Symbols a3, a4, b1, c1, c2, cx CS DL G g2 I k k, l l lF,t n r s S1 • • • 3 T
modeling constants universal constant of the scalar spectrum Markstein diffusivity scalar to define the flame position two point correlation of fluctuations of G expansion integral turbulent kinetic energy wave number vector integral length scale turbulent flame brush thickness normal vector correlation distance burning velocity source terms in the equation for the correlation function transfer function
242
TURBULENT COMBUSTION—Premixed Flame Modeling
t v v⬘ vex x e c C j p¯ q r¯ v¯ x ¯
time flow velocity turbulence intensity induced velocity by expansion effects space vector dissipation expansion parameter scalar spectrum function curvature expansion term in the variance equation density flame surface area ratio scalar dissipation kinematic restoration
Subscripts L T b m u
laminar turbulent burned mid point of the turbulent flame brush unburned
Superscripts 0 reference value * modified value Acknowledgments The work of the first and second author has been financially supported by the Deutsche Forschungsgemeinschaft under grant Pe 241/13-2.
REFERENCES 1. Williams, F. A., and Buckmaster, J. (eds.), The Mathematics of Combustion, SIAM, Philadelphia, 1985, pp. 99–131. 2. Kerstein, A., Ashurst, W., and Williams, F. A., Phys. Rev. A 37(7):2728–2731 (1988). 3. Peters, N., J. Fluid Mech. 242:611–629 (1992). 4. Peters, N., J. Fluid Mech. 384:107–132 (1999). 5. Smiljanovski, V., Moser, V., and Klein, R., Combust. Theory Modelling 1:183–215 (1997). 6. Schneider, T., Terhoeven, P., and Klein, R., “Numerical Investigations of Acoustically Driven Instablilities of Premixed Flames,” Seventeenth Int. Colloq. Dyn. Expl. React. Syst., Heidelberg, Germany, July 25–30, 1999. 7. Sivashinsky, G. I., Acta Astronaut. 4:1177–1206 (1977). 8. Buckmaster, J. D., and Ludford, G. S. S., Theory of Laminar Flames, Cambridge University Press, Cambridge, 1982.
9. Michelson, D. M., and Sivashinsky, G. I., Combust. Flame 48:211–217 (1982). 10. Sivashinsky, G. I., and Clavin, P., J. Phys. 48:193–198 (1987). 11. Frankel, M. L., Phys. Fluids 2:1879–1883 (1990). 12. Matalon, M., and Metzener, P., J. Fluid Mech. 336:331–350 (1997). 13. Bychkov, V. V., Kovalev, K. A., and Liberman, M. A., Phys. Rev. E 60:2897ff (1999). 14. Bychkov, V. V., Phys. Fluids 10(8):2091–2098 (1998). 15. Veynante, D., Trouve´, A., Bray, K. N. C., and Mantel, T., J. Fluid Mech. 332:263–293 (1997). 16. Ashurst, W., Combust. Theory Modelling 1:405–428 (1997). 17. Dandekar, A., and Collins, L. R., Combust. Flame 101:428–440 (1995). 18. Peters, N., Turbulent Combustion, Cambridge University Press, Cambridge, 2000. 19. Wenzel, H., and Peters, N., “Direct Numerical Simulation of the Propagation of a Flame Front in a Homogenous Turbulent Flow Field,” Seventeenth Int. Colloq. Dyn. Expl. React. Syst, Heidelberg, Germany, 1999. 20. Wenzel, H., “Direkte numerische Simulation der Ausbreitung einer Flammenfront in einem homogenen Turbulenzfeld,” Ph.D. thesis, RWTH Aachen, Aachen, Germany, 2000. 21. Batchelor, G. K., The Theory of Homogeneous Turbulence, Cambridge University Press, Cambridge, 1953. 22. Ruetsch, G. R., and Maxey, M. R., Phys. Fluids 3, 6:1587–1597 (1991). 23. Eswaran, V., and Pope, S. B., Comp. Fluids 16:257– 278 (1988). 24. Ulitsky, M., and Collins, L. R., Combust. Flame 111:257–275 (1997). 25. Sethian, J., Level Set Methods (Cambridge Monographs on Applied and Computational Mathematics), Cambridge University Press, Cambridge, 1996. 26. Fornberg, B., A Practical Guide to Pseudo-Spectral Methods (Cambridge Monographs on Applied and Computational Mathematics), Cambridge University Press, Cambridge, 1996. 27. Wenzel, H., Annual Research Briefs, NASA Ames/ Center for Turbulence Research, 1997, Stanford, CA, pp. 237–252. 28. Maz’ja, V. G., Sobolev Spaces, Springer-Verlag, Berlin, 1985, chapters 1, 2, and 4. 29. Vervisch, L., Bidaux, E., Bray, K. N. C., and Kollmann, W., Phys. Fluids 7(10):2496–2503 (1995). 30. Cambray, P., and Joulin G., Proc. Combust. Inst. 24:61–67 (1992).
TURBULENT BURNING VELOCITY
243
COMMENTS John Dold, UMIST, UK. The assumption that strain-rate corrections to laminar flame speed are small is very useful because it means that the solenoidal component of velocity does not need to be considered directly. You said that you have good reasons for this assumption; could you please explain them? Author’s Reply. Modifications of the laminar burning velocity by strain and curvature effects represent higher order terms that cannot reasonably be combined with traditional turbulence modeling which can only take leading order terms into account. As far as curvature is concerned, it becomes a leading order term in the thin reactions zones regime and therefore is worth being considered in a model that covers this regime as well. An additional argument against the inclusion of strain effects in the corrugated flamelets regime is based on scaling. The Markstein length multiplying this term is of the order of the flame thickness which is smaller than the Gibson scale in this regime. Since the Gibson scale represents the spectral cutoff, the strain term is active beyond the cutoff only and has little influence on kinematic restoration [1].
the G-equation model developed by you, in terms of v⬘/sL; L/QL ⳱ L SL/m; Le; t/s1 ⳱ tv⬘/L. Author’s Reply. The G-equation is valid as long as the flame front is continuous, that is, as long as no local quenching occurs. A conservative estimate is that the turbulent Karlovitz number must be smaller than 100. For L/ lF ⳱ 100, this corresponds, for instance, to v⬘/SL 100. In the model the Lewis number is assumed unity. ● J. Chomiak, Chalmers University of Technology, Sweden. You have analyzed the pressure fluctuation effects generated by the laminar flamelet on your turbulent combustion model. However, you have neglected two other effects, the dilatation effects on turbulence and the baroclinic turbulence generation effect, which appear to be much stronger than the pressure fluctuations. In general, the pressure fluctuation generated by the laminar flame fronts are much smaller than the turbulent pressure fluctuations especially for v⬘/SL k 1.
REFERENCE 1. Peters, N., Turbulent Combustion, Cambridge University Press, New York, 2000, p. 108. ● Andrei Lipatnikov, Chalmers University of Technology, Sweden. Could you please specify the domain of validity of
Author’s Reply. We have not analyzed pressure fluctuation effects but the influence of gas expansion modeled by volume sources on the flow field. This accounts for dilatation but not for baroclinic turbulence generation. We believe that since the baroclinic effect scales with SL the main conclusion of the paper, namely that heat release effects do not strongly influence the turbulent burning velocity for v⬘/SL k 1, remains valid.
MODIFICATION OF THE TURBULENT BURNING VELOCITY BY GAS EXPANSION N. PETERS,1 H. WENZEL1 and F. A. WILLIAMS2 1
Institut fu¨r Technische Mechanik RWTH Aachen, D-52056 Aachen, Germany 2 Center for Energy and Combustion Research Department of Mechanical and Aerospace Engineering University of California at San Diego La Jolla, California 92093-041, USA
Premixed turbulent combustion in the corrugated flamelets regime has been addressed by modeling of a G-equation containing a gas-expansion term, the form of which was derived by generalizing results of Sivashinsky along lines initiated by Frankel. Previous analyses by the first author of the mean, variance, and two-point correlation spectrum of G, as well as their relation to the turbulent burning velocity, were extended to include the influences of gas expansion. It is shown that this effect has the same form as kinematic restoration but acts in the opposite direction. Because of this similarity in form, no new modeling approximations are required. On average, kinematic restoration dominates gas expansion. In physical terms, this means that the energy-containing eddies at the integral scale are strong enough to suppress gasdynamic instabilities at smaller scales. Gas expansion lessens the smoothing effect of kinematic restoration, thereby increasing the variance of G and the turbulent flame brush thickness. It also increases the turbulent burning velocity by a factor that increases with increasing ratio of unburned to burned gas density and with decreasing ratio of turbulence intensity to laminar burning velocity. The effect on the burning velocity, however, is not large, remaining less than 40% even as the turbulence intensity begins to decrease below the laminar burning velocity.
Introduction
sT ⳱ b1v⬘
In recent years, several attempts have been made to derive models for premixed turbulent combustion based on a level-set approach using the so-called Gequation [1,2],
where b1 is a constant of proportionality which has been determined from experimental data as b1 ⳱ 2.0. More recently [4], the level-set approach was used to describe a new regime in turbulent premixed combustion called the thin reaction zones regime, where the smallest turbulent structures may enter the preheat zone of the flame but not the reaction zone. In the thin reaction zones regime, the diffusive effects become dominant. In equation 1, gas-expansion effects enter through changes of the velocity field v induced by heat release at the flame surface. The impact of gas expansion can be analyzed from jump conditions for the flow variables across the flame surface. Numerical methods using local reconstruction within a numerical cell have been derived [5,6] to do this accurately. The objective of the present paper is to derive closure relations as in Ref. [3] but with the effect of gas expansion included. For that purpose, we start from a modified version of the G-equation. The foundation for this modification was laid by Sivashinsky [7]. The original Sivashinsky equation has been subsequently analyzed and extended in a number of papers (e.g., [8–13]).
G Ⳮ v • G ⳱ sL|G| t
(1)
The G-field is convected by the flow velocity v and propagates into the unburned mixture with the laminar burning velocity sL. An isosurface G(x, t) ⳱ G0, where G0 is a fixed constant, describes the location of the flame front in an arbitrary flow field. The levelset approach was used in Ref. [3] to derive a twopoint closure relation in a regime of premixed turbulent combustion that we call here the corrugated flamelets regime. It has also been called the multiple-sheet or fractal regime in the literature. That regime is characterized by a local kinematic balance between flow velocity and burning velocity. It was also shown in Ref. [3] that the turbulent burning velocity sT becomes independent of the laminar burning velocity sL and is proportional to the turbulence intensity v⬘ in the limit of large values of v⬘/ sL,
235
(2)
236
TURBULENT COMBUSTION—Premixed Flame Modeling
To facilitate analysis, we adopted a model that does not account for the all physical effects that occur in a premixed turbulent flame. Gas-expansion effects are treated as volume sources located on the flame front. The volume sources generate a potential flow field which then is simply superposed on a turbulent flow field that is described by the NavierStokes equations for constant density flows. The resulting velocity field does not satisfy the non-constant density Navier-Stokes equations and should be considered as a model. This model does not account, for instance, for the vorticity generated by the flame. While this exerts influences, it may be expected that it will not qualitatively affect the results. Recently, Bychkov [14] has shown that inclusion of vorticity generation in his model equation for curved stationary laminar flames reduces the burning velocity by 40% for an unburned-to-burned density ratio qu/qb of 5 but does not introduce important qualitative modifications. The modeling of gas expansion by a source term in the G-equation has the advantage that its effect can be easily analyzed and quantified, a task that is not trivial if one analyses data of direct numerical simulations with simplified or full chemistry and heat release [15]. The modeling of the additional source term in the G-equation is addressed first. A numerical analysis is then presented in which a volume source at each point of the flame surface, previously used in a numerical study by Ashurst [16], is employed to circumvent the restriction of small departures of the flame sheet from a plane. Frankel [11] has shown for the two-dimensional case that in the limit of small flame front perturbations, the use of volume sources is equivalent to the non-local expansion term appearing in the Sivashinsky equation [7]. Finally, the main results are summarized and interpreted in terms of the turbulent burning velocity.
mixture immediately ahead of the flame structure, sL ⳱ sLu and v ⳱ vu in equation 4. If, on the other side, the flame front position is defined to lie in the burned gas immediately behind the flame structure, sL ⳱ sLb and v ⳱ vb in equation 4. Introducing equation 4 into equation 3 gives, for the velocity jump normal to the flame front, the expression c s ⳱ csLb 1 ⳮ c Lu
n • (vu ⳮ vb) ⳱
where c ⳱ (qu ⳮ qb)/qu is the gas-expansion parameter, which is zero for constant density but close to unity for large values of qu/qb. Equation 5 shows that gas-expansion affects scale with the laminar burning velocity in the sense that the velocity jump is proportional to sL and suggests that they will be important only if sL is of the same order as the turbulent velocity fluctuations v⬘ or larger. These effects are important in the corrugated flamelets regime but are likely to be insignificant in the thin reaction zones regime, which therefore will not be considered. In the work of Sivashinsky [7], an equation for small perturbations of a planar two-dimensional flame front in the limit of small c was derived, accounting for stabilizing or destabilizing diffusionalthermal instabilities associated with negative Markstein length and the Darrieus-Landau instability arising from gas expansion. For simplicity, we do not address the former here, setting the Markstein length Lb small but positive. The latter appears in the equation for the perturbation F(y, z, t) of the flame front location as an integral over transverse coordinates. We introduce G ⳱ x ⳮ F(y, z, t) and write Sivashinsky’s integral as
1 I(G,x,t) ⳱ (2p)3
冮冮冮冮冮冮 ⳮ
ik•(xⳮx˜)
|k| e Formulation and Modeling In the thin flame approximation, the the mass balance across the flame is
冢
qun • vu ⳮ
冣
冢
冣
dxf dx ⳱ qbn • vb ⳮ f dt dt
(3)
where qu and qb, and vu and vb are the densities and velocities of the unburned and burned gases, respectively, n is the local normal vector of the flame front pointing toward the unburned mixture, and dxf ⳱ v Ⳮ sLn dt
(4)
is the propagation velocity of the front into the unburned mixture with respect to the unburned gas, with sL being the laminar burning velocity. If the position of the front is defined to lie in the unburned
(5)
G(x˜,t) dx˜dk
(6)
As a generalization of Sivashinsky’s equation, the resulting G-equation on which we base our modeling is G Ⳮ v • G ⳱ s0Lb|G| ⳮ DLbj|G| t Ⳮ
c 0 s I(G,x,t) 2 Lb
(7)
where s0Lb is the burning velocity of an unstretched laminar flame with respect to the burned gas, 0 DLb ⳱ sLb Lb is the corresponding Markstein diffusivity, and j is the curvature. In a numerical study, Dandekar and Collins [17] used a similar formulation. When the assumption of small perturbations about a plane is removed, the integral becomes a volume integral [11]. Further details of the derivation may be found in Ref. [18].
TURBULENT BURNING VELOCITY
¯ , the To develop equations for the mean value of G variance G⬘2, and the two-point correlation function g2(r, t) ⳱ G⬘(x, t)G⬘(x Ⳮ r, t), we regard equation 7 as a three-dimensional field equation in a constantdensity flow field and split the variables G, r ⳱ |G| ¯ and v into mean and fluctuating parts [3], G ⳱ G Ⳮ G⬘, r ⳱ r¯ Ⳮ r⬘, v ⳱ v¯ Ⳮ v⬘. The equation for ¯ is readily obtained as G ¯ G ¯ Ⳮ • v⬘G⬘ ⳱ Ⳮ v¯ • G t c ¯ I(G,x,t) (8) 2 ¯ in the argument of I(G ¯ , x, t) is where the use of G due to the linearity of the expansion integral I(G, x, t) with respect to G. This linearity is an important feature and is exploited in analyzing the variance equation and the equation for the two-point correlation. A surprising result of equation 8 is that for a freely propagating planar turbulent flame, the expansion ¯ , x, t), being an odd function of the coorterm I(G dinate with respect to the flame in the direction of propagation, vanishes and has no direct influence on the propagation speed of that flame. The expansion ¯ only by modifying effects enter the equation for G the flame surface area ratio r. ¯ In this formulation, the variance G⬘2 of G is related to the turbulent flame brush thickness lF,t (cf. equation 30 below). The equation for its dynamic evolution is obtained by subtracting the equation for ¯ (equation 8) from the instantaneous G-equation G (equation 7) to get the equation for the fluctuation of G. Multiplying this by G⬘ and then averaging yields the equation for the variance of G, namely,
p¯ ⳱
s0Lb (2p)3
冮冮冮冮冮冮|k| e
ik•(xⳮx˜)
G⬘(x˜,t)G⬘(x,t)dx˜dk
(9)
(10)
Numerical simulations [19,20] have shown that these terms always act as sinks for the variance. The physical interpretation of kinematic restoration is that wrinkles introduced by turbulent velocity fluctuations are removed by the smoothing effect of local laminar flame propagation, tending to restore a flatter surface. Scalar dissipation is a higher-order
(11)
ⳮ
Like kinematic restoration, this quantity is proportional to the laminar burning velocity, but differently from x, ¯ it acts as a source term for the variance and increases irregularities of the flame front similar to the way that the Darieus-Landau instability increases initial small disturbances of a planar flame in a quiescent medium. At this stage, there is no information on the relative orders of magnitude of these two counteracting terms and consequently no forecast can be made as to which will dominate in a turbulent premixed flame. Some information relevant to this can be gained by proceeding with the formulation of the conservation equation for the two-point correlation function along the lines described previously [3]. Following standard methods [21] for homogeneous isotropic turbulence, all spatial dependencies are reduced to the correlation distance r ⳱ |r|, and the correlation function obeys g2(r,t) t Ⳮ 2 • v⬘(x Ⳮ r,t)G⬘(x,t)G⬘(x Ⳮ r,t) 0 Ⳮ 2s0LbS1 Ⳮ 2DLbS2 ⳮ csLb S3 ⳱ 0
The terms on the left-hand side are local accumulation of variance, transport by the mean velocity, and turbulent transport. The first term on the righthand side is the production of variance by the action of velocity fluctuations on the mean gradient of G. The next terms on the right-hand side are the kinematic restoration and the scalar dissipation [3], defined here as x ¯ ⳱ ⳮ2s0LbG⬘r⬘, v¯ ⳱ 2DLb(jr)⬘G⬘
term in the corrugated flamelets regime as shown in Ref. [4]. It becomes more important, however, in cases where gas-expansion effects tend to cancel kinematic restoration, as is seen below. The last term on the right-hand side of equation 9, the covariance p¯ of fluctuations of the G-field and the expansion integral I, is defined as
0 s0Lbr¯ ⳮ DLbj|G| Ⳮ sLb
G⬘2 Ⳮ v¯ • G⬘2 Ⳮ • v⬘G⬘2 ⳱ t ¯ ⳮx ⳮ2v⬘G⬘ • G ¯ ⳮ v¯ Ⳮ cp¯
237
(12)
Here the two-point correlations between r⬘ and G⬘ and between (jr)⬘ and G⬘ have been denoted, respectively, by [3] S1(r,t) ⳱ ⳮr⬘(x,t)G⬘(x Ⳮ r,t), S2(r,t) ⳱ ⳮ(j(x,t)r(x,t))⬘G⬘(x Ⳮ r,t)
(13)
The last term of equation 12 is defined as S3(r,t) ⳱ I(G,x,t)G⬘(x Ⳮ r,t)
1 ⳱ (2p)3
冮冮冮冮冮冮 |l| e
il•(xⳮx˜)
ⳮ
G⬘(x˜,t)G⬘(x Ⳮ r,t) dx˜dl
(14)
and describes the influence of gas expansion on the evolution g2(r, t). To close the equation for the two-point correlation function, we follow Ref. [3] and transform it into
238
TURBULENT COMBUSTION—Premixed Flame Modeling
wave-number space k by applying Fourier transforms to obtain the scalar spectrum function C(k, t). With use made of the transform pair
冮冮冮 f(r,t)e
ikr
dr,
ⳮ
f(r,t)⳱
冮冮冮
ⳮikr ˆ f(k,t)e dk
(15)
ⳮ
for any function f(r, t), the spectrum function C(k, t) of fluctuations of G is related to the Fourier transform of g2(r, t) by
冮
C(k,t) ⳱ k2 䡬 gˆ2 (k,t)dX ⳱ 4pk2gˆ2(k,t)
(16)
where integration is carried over the solid angle dX, and k ⳱ |k| is the absolute value of the wave number. Taking the Fourier transform of the divergence of the triple correlation v⬘(x Ⳮ r, t) G⬘ (x, t) G⬘ (x Ⳮ r, t) leads to a spectral transport term T(k, t), which was approximated in Ref. [3] by a gradient transport assumption, thereby making it local in wave-number space. Dimensional arguments and the assumption of locality led in Ref. [3] to linear approximations for the Fourier transforms Sˆ1 and Sˆ2, ⳮ1
2ˆ
ⳮ1 2
2ˆ
8pk S1 ⳱ c1Cs kC, 8pk S2 ⳱ c2Cs k C
(17)
where Cs is the universal constant of the scalar spectrum, and c1 is a modeling constant that is assumed to be independent of the wave number to agree with the locality hypothesis. To derive the Fourier transform of S3, we introduce vector s, defined by x Ⳮ r ⳱ x˜ Ⳮ s, into equation 14 and change the integration variable x˜ to s and l to ⳮk. Since x and r are fixed, dx˜ ⳱ ⳮds and dl ⳱ ⳮdk, and the definition for S3 becomes
S3 ⳱
冮冮冮 ⳮ
C(k,t) 0 ⳮ T(k,t) Ⳮ c1Csⳮ1sLb kC(k,t) t 1 2 0 Ⳮ c2Cⳮ S DLbk C(k,t) ⳮ csLbkC(k,t) ⳱ 0
G⬘(x˜,t)G⬘(x˜ Ⳮ s,t) dsdk
iks
The noteworthy point about this equation is that the term that accounts for gas-expansion effects has the same wave-number dependence as the kinematic restoration. These two terms may therefore be combined to give C(k,t) ⳮ1 0 ⳮ T(k,t) Ⳮ c*C 1 s sLbkC(k,t) t 1 2 Ⳮ c2Cⳮ S DLbk C(k,t) ⳱ 0
G⬘2 ⳱ ⳮx* ⳮ v¯ t
⳱
冮冮冮 |k|gˆ (k,t)e
dk
(22)
is recovered, where the variance, the modified kinematic restoration, and the scalar dissipation have been defined as G⬘2 ⳱
冮
k0
C dk,
ⳮ1 0 x* ⳱ x ¯ ⳮ cp¯ ⳱ c*C 1 s sLb
1 v¯ ⳱ c2Cⳮ s DLb
冮
k0
冮
k0
kC dk,
k2C dk
x* ⳱ c*x
ikr
(21)
with c*1 ⳱ c1 ⳮ cCs. Upon integrating equation 21 over wave-number space, the spectral transfer term vanishes, and the variance equation for homogeneous isotropic turbulence,
ⳮ
2
(20)
(23)
respectively. The lower limit of integration in these integrals is taken to be the inverse of the integral turbulent length scale k0 ⳱ lⳮ1. Provided that c*1 0 for k lⳮ1, the modified kinematic restoration reduces fluctuations of the flame front. It then can be modeled in the same way as x ¯ in Ref. [3], namely, as
|k|eikr (2p)3
冮冮冮 e
(19)
We are now able to write the dynamical equation for the spectrum function C(k, t) by taking the Fourier transform of equation 12 and applying the definition of C(k, t) in equation 16, obtaining
1 ˆ f(k,t) ⳱ (2p)3
C(k,t) Sˆ 3(k,t) ⳱ |k|gˆ2(k,t) ⳱ 4pk
(18)
e G⬘2 k
(24)
ⳮ
Comparing this with the definition of the inverse Fourier transform, equation 15, gives, due to linearity of S3 in g2, with the use of equation 16 the closed form
Numerical Simulation Numerical simulations can check the validity of some of the very strong assumptions that were made in the previous section. Results from two different
TURBULENT BURNING VELOCITY
239
the cusps without adding numerical diffusion to the problem, because no flux limiting is used, and the arising Riemann problem is solved exactly. This procedure has been described elsewhere [20,27]. Fig. 1 shows a sketch of the model used to incorporate gas expansion in the numerical simulation. At every point on the flame surface G ⳱ G0, there is a volume source that induces a velocity into the whole flow field. We are interested in the induced normal velocity vex(x) ⳱ n(x) • vn(x) only at the flame front itself, which can be calculated as vex(x) ⳱ c Fig. 1. A schematic illustration of the velocity induced at the flame front position x by the effect of gas expansion originating from x¯.
numerical simulations are presented, one using the limit c ⳱ 0 to test the model that has led to the first of equation in equation 17 and to calculate the numerical value of c1, the other including the effects of gas expansion for values of c between 0.5 and 0.8 corresponding to density ratios qu/qb between 2.0 and 5.0. It is found that that these effects become less important for larger values of the turbulence intensity v⬘/sL. The model introduced above allows for the turbulent flow field to be computed independently of the G-field. This flow field is calculated in a cubic box using a pseudospectral code developed by Ruetsch [22] which employs forcing of the large scales to allow for a statistically stationary homogeneous isotropic velocity field following a method developed by Pope and Eswaran [23]. Thereby, velocity fluctuations induce fluctuations of G. The forcing of the scalar field is accomplished by imposing a mean gradient in the x direction [22,24], G(x,y,z,t) ⳱ x Ⳮ g(x,y,z,t)
冮
A
(x˜ ⳮ x) • n(x) dA |x˜ ⳮ x|3
(26)
in the three-dimensional case. These integrals of the Biot-Savart type have to be taken over the whole flame front. Although it is not immediately clear that the integrand stays bounded in the limit x˜ → x, it can be shown to converge to a finite value, which turns out to be that of the local curvature [20]. Similarly to Ashurst [16], we include the expansion effects into the G-equation by modifying the laminar burning velocity with the induced velocity,
冢
sL ⳱ sLb 1 ⳮ
冣
c Ⳮ vex(x) 2
(27)
and insert it into the G-equation (equation 1). The factor 1 ⳮ c/2 originates according to equation 4 from the definition of the flame surface G(x, t) ⳱ G0 at the midpoint of the flame structure. A direct evaluation of the integral in equation 26 by triangulation proved to be numerically too expensive. Therefore, a transformation introduced by Maz’ja [28] was used to express the surface integral as a volume integral according to
冮冮 f(r,x) dS ⳱ 冮冮冮 f(r,x)d(G(r) ⳮ G ) S
v
(25)
thereby simulating a planar turbulent flame. Note that this decomposition allows for a multivalued solution in the direction of the mean flame propagation and is therefore more general than the decomposition employing F(y,z,t). Since the solution for g(x, y, z, t) is periodic in all three space dimensions, the same pseudospectral method that is used to calculate the velocity field can be used to compute the scalar field. Equations of the G-equation type (equation 1) describe in the limit of constant sL ⳱ s0L a front propagation of the Huygens type [25]. The front develops cusps which lead to unstable solutions if a standard pseudospectral method is used to calculate the absolute value of the gradient of G on the right-hand side of the G-equation. Following Ref. [26], an upwind method was developed to account for this behavior. This upwind scheme is capable of capturing
sLb 4p
|G(r)|drxdrydrz
0
(28)
where d(G(r) ⳮ G0) is the Dirac delta function, and f(r) is the integrand in equation 26. An equivalent procedure to evaluate surface integrals was used by Vervisch et al. [29] to compute the surface area of premixed turbulent flames.
Results 3
On a 64 grid on a cubic domain of edge length 2p, the evolution of the G-equation was computed first for c ⳱ 0 for three different ratios of the turbulence intensity v⬘ to the laminar burning velocity s0L. The root-mean-square value of the normalized velocity fluctuations was v⬘ ⳱ 17.4, the normalized integral length scale was l ⳱ 1.08, and the normalized kinematic viscosity was m ⳱ 0.24. This results in a Reynolds number of 42 based on the Taylor
240
TURBULENT COMBUSTION—Premixed Flame Modeling
Fig. 2. Evaluation from DNS of the coefficient c1 in the model for the Fourier transform of the two-point correlation representing the kinematic restoration. Dotted line, v⬘/sL ⳱ 1.0; solid line, v⬘/sL ⳱ 1.6; dashed line, v⬘/sL ⳱ 2.0.
microscale. After a statistical stationary stage was reached, the spectrum function and the Fourier transform Sˆ1 of the kinematic restoration was calculated, and the coefficient c1 was then evaluated from equation 17 using CS ⳱ 1.2 for the universal constant of the scalar spectrum. The results, given in Fig. 2, show that this coefficient is larger than unity for k k0. Therefore, c*1 ⳱ c1 ⳮ cCs is positive in the inertial range of the scalar spectrum c ⱕ 0.8, corresponding to density ratios qu/qb ⱕ 5. Even though these results were obtained at c ⳱ 0, they provide a first confirmation of the hypothesis that kinematic restoration dominates gas-expansion effects in most cases of practical interest. Since the numerical evaluation of the non-local expansion integral in equation 28 is expensive, computations for c ⬆ 0 were confined to a 323 grid. Attention was focused first on the behavior of the mean gradient r, ¯ which can be interpreted as the turbulent flame surface area ratio, which in turn is proportional to the ratio of the turbulent to the laminar burning velocity [2,3,19]. Fig. 3 shows the ratio of the turbulent burning velocity thus calculated with and without gas-expansion effects as a function of c and v⬘/sLm with
冢
sLm ⳱ sLb 1 ⳮ
冣
c 2
(29)
The increase of the turbulent burning velocity with increasing c is substantial only for the smallest value of the turbulence intensity shown, where it is still less than 40%; it is only 12.5% for v⬘/sL ⳱ 3.0. To further test the hypothesis that the kinematical restoration dominates the destabilizing effect of the expansion term, the time series of the source terms of the variance equation (equation 9) is plotted in
Fig. 3. DNS results for the ratio of the turbulent burning velocity sT to the value sT1c. Both burning velocities were calculated as the long time averages for flame surface area ratio r, ¯ sT with volume sources in the G-equation included, sT,c ⳱ 0 without volume sources.
Fig. 4. Time series of the source terms in the variance equation. The turbulence intensity v⬘/sLm is 2/3 and c ⳱ 0.8. The time is normalized with the integral time scale and the other terms with sLm and lF,t.
Fig. 4. For numerical reasons, it is necessary to include the curvature term with a positive Markstein diffusivity DLb of 0.48 in the G-equation (equation 7). As can be seen from that figure, the scalar dissipation v, ¯ which accounts for curvature effects, is significantly smaller than the other two terms, which is consistent with the analysis in Refs. [4] and [3], where it was shown that scalar dissipation is a higherorder contribution in the corrugated flamelets regime. Fig. 4 also shows that the kinematical restoration term x ¯ closely tracks the expansion term p, ¯ dominating it most of the time. It is perhaps surprising how often p¯ exceeds x, ¯ leading to a relatively small value c*x ⳱ 0.5 in equation 24 for this example. The Turbulent Flame Brush Thickness and the Turbulent Burning Velocity It can be shown [4,18] that in the limit of a stationary planar premixed flame, all gradients of G⬘2
TURBULENT BURNING VELOCITY
vanish, and the production term in equation 9 may ¯ |2, where Dt ⳱ a4v⬘l is the be modeled as 2Dt|G turbulent diffusivity. Neglecting the scalar dissipation term and using the model in equation 24 for x* with k/e ⳱ a3l/v⬘, one obtains for the turbulent flame brush thickness ᐉF,t ⳱
冢
1/2
冣
(G⬘2)1/2 2a3a4 ⳱ ¯ c*x |G|
ᐉ
(30)
The coefficients a3 ⳱ 4.05 and a4 ⳱ 0.78 have been evaluated in Ref. [18] using standard models of turbulence. For constant density, the coefficient c*x becomes equal to the value cx ⳱ 1.62 obtained in Ref. [3]. Since gas-expansion effects lead to smaller values of c*x, it follows from equation 30 that the turbulent flame brush thickness increases accordingly. The model for the kinematic restoration x ¯ has been used together with dimensional arguments in Refs. [3] and [4] to relate the turbulent burning velocity sT and the flame surface area ratio r¯ for a stationary planar flame in the corrugated flamelets regime to x ¯ as 1/2
冢 ke冣
sT ⳱ sLr¯ x ¯
(31)
¯ is replaced by x* for the case where expansion If x effects are included, it follows with equations 24 and 30 that sT
2 1/2
冢c*a G⬘ᐉ 冣 x
2 3
2
1/2
v⬘
冢2aa 冣 4
v⬘
(32)
3
¯ | ⳱ 1 has been used without loss of genwhere |G erality (cf. [18]). The last expression in this relation shows that in the limit of large values of v⬘/sL, the turbulent burning velocity is independent of c*x and therefore of the gas-expansion effects modeled by it. The modeling constant b1 in equation 2 is therefore expected to be independent of the gas-expansion parameter c. Conclusions The preceding simulation results support the general validity of the model in Ref. [3] that led to equation 2 in the corrugated flamelets regime. Kinematic restoration counteracts and dominates gas-expansion effects on the time average, but their balance may instantaneously be reversed. Kinematic restoration is the basis of the cascade hypothesis for scalar fluctuations; its reduction by gas-expansion effects can become substantial, leading to a turbulent flame thickness appreciably larger than obtained without gas expansion. This is consistent with an earlier simulation [16]. The associated increase in the wrinkled flame area provides, however, only a slightly higher turbulent burning velocity. The simulation results shown in Fig. 3 provide this enhancement in
241
the turbulent burning velocity as the gas-expansion parameter c increases. For the realistic value c ⳱ 0.8, gas expansion increases sT by about 30% for v⬘ ⳱ sLm and only about 12% for v⬘ ⳱ 3sLm. The influence of gas expansion on sT therefore is significantly reduced by increasing the turbulence intensity. This is also shown to be consistent with previous modeling. In high-intensity turbulence, gas expansion therefore can be neglected in models of the turbulent burning velocity. The effect of gas expansion does, however, increase as the turbulence intensity decreases, and in the weak turbulence limit the effect can be dominant, as earlier studies have suggested [30]. A consequence of these findings is that the turbulent burning velocity is a well-defined turbulent quantity that depends on other local mean quantities such as the local turbulence intensity. Non-local effects associated with gas expansion have no influence on scales smaller than the integral length scale. On larger scales, however, where kinematic restoration is unable to control gas-expansion effects, large-scale instabilities may still develop. This may be deduced from Fig. 2 where c1 1 for k lⳮ1, leading to c*1 0 for realistic values of c. Large-scale instabilities of a turbulent front are therefore equivalent to those of the Darieus-Landau type of a laminar front, with the laminar burning velocity being replaced by the turbulent burning velocity, and with a stabilizing effect occurring at the integral length scale. Since only positive values of the Markstein length were used for numerical reasons, it is not clear, whether the destabilizing effect of a negative Markstein length, which would lead to a negative scalar dissipation rate, could unsettle the delicate balance between kinematic restoration and gas-expansion effects.
Nomenclature Symbols a3, a4, b1, c1, c2, cx CS DL G g2 I k k, l l lF,t n r s S1 • • • 3 T
modeling constants universal constant of the scalar spectrum Markstein diffusivity scalar to define the flame position two point correlation of fluctuations of G expansion integral turbulent kinetic energy wave number vector integral length scale turbulent flame brush thickness normal vector correlation distance burning velocity source terms in the equation for the correlation function transfer function
242
TURBULENT COMBUSTION—Premixed Flame Modeling
t v v⬘ vex x e c C j p¯ q r¯ v¯ x ¯
time flow velocity turbulence intensity induced velocity by expansion effects space vector dissipation expansion parameter scalar spectrum function curvature expansion term in the variance equation density flame surface area ratio scalar dissipation kinematic restoration
Subscripts L T b m u
laminar turbulent burned mid point of the turbulent flame brush unburned
Superscripts 0 reference value * modified value Acknowledgments The work of the first and second author has been financially supported by the Deutsche Forschungsgemeinschaft under grant Pe 241/13-2.
REFERENCES 1. Williams, F. A., and Buckmaster, J. (eds.), The Mathematics of Combustion, SIAM, Philadelphia, 1985, pp. 99–131. 2. Kerstein, A., Ashurst, W., and Williams, F. A., Phys. Rev. A 37(7):2728–2731 (1988). 3. Peters, N., J. Fluid Mech. 242:611–629 (1992). 4. Peters, N., J. Fluid Mech. 384:107–132 (1999). 5. Smiljanovski, V., Moser, V., and Klein, R., Combust. Theory Modelling 1:183–215 (1997). 6. Schneider, T., Terhoeven, P., and Klein, R., “Numerical Investigations of Acoustically Driven Instablilities of Premixed Flames,” Seventeenth Int. Colloq. Dyn. Expl. React. Syst., Heidelberg, Germany, July 25–30, 1999. 7. Sivashinsky, G. I., Acta Astronaut. 4:1177–1206 (1977). 8. Buckmaster, J. D., and Ludford, G. S. S., Theory of Laminar Flames, Cambridge University Press, Cambridge, 1982.
9. Michelson, D. M., and Sivashinsky, G. I., Combust. Flame 48:211–217 (1982). 10. Sivashinsky, G. I., and Clavin, P., J. Phys. 48:193–198 (1987). 11. Frankel, M. L., Phys. Fluids 2:1879–1883 (1990). 12. Matalon, M., and Metzener, P., J. Fluid Mech. 336:331–350 (1997). 13. Bychkov, V. V., Kovalev, K. A., and Liberman, M. A., Phys. Rev. E 60:2897ff (1999). 14. Bychkov, V. V., Phys. Fluids 10(8):2091–2098 (1998). 15. Veynante, D., Trouve´, A., Bray, K. N. C., and Mantel, T., J. Fluid Mech. 332:263–293 (1997). 16. Ashurst, W., Combust. Theory Modelling 1:405–428 (1997). 17. Dandekar, A., and Collins, L. R., Combust. Flame 101:428–440 (1995). 18. Peters, N., Turbulent Combustion, Cambridge University Press, Cambridge, 2000. 19. Wenzel, H., and Peters, N., “Direct Numerical Simulation of the Propagation of a Flame Front in a Homogenous Turbulent Flow Field,” Seventeenth Int. Colloq. Dyn. Expl. React. Syst, Heidelberg, Germany, 1999. 20. Wenzel, H., “Direkte numerische Simulation der Ausbreitung einer Flammenfront in einem homogenen Turbulenzfeld,” Ph.D. thesis, RWTH Aachen, Aachen, Germany, 2000. 21. Batchelor, G. K., The Theory of Homogeneous Turbulence, Cambridge University Press, Cambridge, 1953. 22. Ruetsch, G. R., and Maxey, M. R., Phys. Fluids 3, 6:1587–1597 (1991). 23. Eswaran, V., and Pope, S. B., Comp. Fluids 16:257– 278 (1988). 24. Ulitsky, M., and Collins, L. R., Combust. Flame 111:257–275 (1997). 25. Sethian, J., Level Set Methods (Cambridge Monographs on Applied and Computational Mathematics), Cambridge University Press, Cambridge, 1996. 26. Fornberg, B., A Practical Guide to Pseudo-Spectral Methods (Cambridge Monographs on Applied and Computational Mathematics), Cambridge University Press, Cambridge, 1996. 27. Wenzel, H., Annual Research Briefs, NASA Ames/ Center for Turbulence Research, 1997, Stanford, CA, pp. 237–252. 28. Maz’ja, V. G., Sobolev Spaces, Springer-Verlag, Berlin, 1985, chapters 1, 2, and 4. 29. Vervisch, L., Bidaux, E., Bray, K. N. C., and Kollmann, W., Phys. Fluids 7(10):2496–2503 (1995). 30. Cambray, P., and Joulin G., Proc. Combust. Inst. 24:61–67 (1992).
TURBULENT BURNING VELOCITY
243
COMMENTS John Dold, UMIST, UK. The assumption that strain-rate corrections to laminar flame speed are small is very useful because it means that the solenoidal component of velocity does not need to be considered directly. You said that you have good reasons for this assumption; could you please explain them? Author’s Reply. Modifications of the laminar burning velocity by strain and curvature effects represent higher order terms that cannot reasonably be combined with traditional turbulence modeling which can only take leading order terms into account. As far as curvature is concerned, it becomes a leading order term in the thin reactions zones regime and therefore is worth being considered in a model that covers this regime as well. An additional argument against the inclusion of strain effects in the corrugated flamelets regime is based on scaling. The Markstein length multiplying this term is of the order of the flame thickness which is smaller than the Gibson scale in this regime. Since the Gibson scale represents the spectral cutoff, the strain term is active beyond the cutoff only and has little influence on kinematic restoration [1].
the G-equation model developed by you, in terms of v⬘/sL; L/QL ⳱ L SL/m; Le; t/s1 ⳱ tv⬘/L. Author’s Reply. The G-equation is valid as long as the flame front is continuous, that is, as long as no local quenching occurs. A conservative estimate is that the turbulent Karlovitz number must be smaller than 100. For L/ lF ⳱ 100, this corresponds, for instance, to v⬘/SL 100. In the model the Lewis number is assumed unity. ● J. Chomiak, Chalmers University of Technology, Sweden. You have analyzed the pressure fluctuation effects generated by the laminar flamelet on your turbulent combustion model. However, you have neglected two other effects, the dilatation effects on turbulence and the baroclinic turbulence generation effect, which appear to be much stronger than the pressure fluctuations. In general, the pressure fluctuation generated by the laminar flame fronts are much smaller than the turbulent pressure fluctuations especially for v⬘/SL k 1.
REFERENCE 1. Peters, N., Turbulent Combustion, Cambridge University Press, New York, 2000, p. 108. ● Andrei Lipatnikov, Chalmers University of Technology, Sweden. Could you please specify the domain of validity of
Author’s Reply. We have not analyzed pressure fluctuation effects but the influence of gas expansion modeled by volume sources on the flow field. This accounts for dilatation but not for baroclinic turbulence generation. We believe that since the baroclinic effect scales with SL the main conclusion of the paper, namely that heat release effects do not strongly influence the turbulent burning velocity for v⬘/SL k 1, remains valid.