- Email: [email protected]
Strain rate effects on the nonlinear development of hydrodynamically unstable flames
Available online at www.sciencedirect.com
Proceedings of the
Proceedings of the Combustion Institute 33 (2011) 1087–1094
Combustion Institute www.elsevier.com/locate/proci
Strain rate effects on the nonlinear development of hydrodynamically unstable flames F. Creta *, M. Matalon Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Available online 14 August 2010
Abstract In this study we numerically implement the hydrodynamic model for a premixed flame as a nonlinear free boundary problem where the flame is tracked via a level set equation and the flow is described by a solution of the variable density Navier–Stokes equations. Unlike an earlier similar study, the present model is enriched by fully accounting for hydrodynamic strain in the flame stretch relation which, in turn, affects the local flame speed. The objective is to comprehensively analyze the effect of strain on the onset of the hydrodynamic instability and on the nonlinear development that takes place beyond its inception. The initial evolution is corroborated with the results of a linear stability analysis for which strain rate effects are fully included. We show that while strain provides an additional stabilizing effect on the short wavelength disturbances, thereby delaying the onset of the hydrodynamic instability, it acts to sharpen the cusps near the troughs of the corrugated flame that develops beyond the stability threshold resulting in a larger flame surface area and a higher propagation speed. Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Premixed flames; Thermal expansion; Hydrodynamic instability; Strain rate; Flame stretch
1. Introduction The study of flame dynamics in the past decades has been conducted utilizing a wide range of numerical and perturbative techniques. Early linear stability analyses have led Darrieus [1] and Landau [2] to establish the unconditional instability of a planar flame front propagating at constant speed. Although stable flames were observed in the laboratory owing to stabilizing influences of
* Corresponding author. Address: Dept. of Mechanical Science and Engineering, University of Illinois at Urbana Champaign, 1206 W. Green St., MC244, Urbana, IL 61801, USA. Fax: +1 217 333 1942. E-mail addresses: [email protected] (F. Creta), matalon @illinois.edu (M. Matalon).
diffusion, hydrodynamic instabilities, caused by the substantial thermal expansion, have been observed in experiments [3–5] as large scale cusp-like structures essentially controlled by the system’s domain size. Further insight on the occurrence of instabilities and on the stability limits was obtained by Markstein [6] who introduced the dependence on curvature of the flame speed through a phenomenological model. The more rigorous asymptotic treatment of multi-dimensional premixed flames exploits the multi-scale nature of the problem, which is characterized by two disparate length scales: the diffusion length ‘f that characterizes the flame thickness and the hydrodynamic length L associated with the geometrical dimensions of the vessel within which combustion takes place. Such an analysis was developed by Clavin and Williams
1540-7489/$ - see front matter Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2010.06.029
1088
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
[7] for a weakly wrinkled flame in low intensity turbulence and Oð1Þ expansion ratios r of the unburnt to burnt densities. Matalon and Matkowsky [8] eliminated the assumptions of small flame displacements and weak flow nonuniformity, and formulated a general nonlinear hydrodynamic model in which the flame is a gasdynamic discontinuity separating burnt and unburnt gases. The derived flame speed was found to depend on both curvature and hydrodynamic strain, effects that constitute the total stretch rate experienced by the flame. Stability analyses within this context provided more comprehensive results, incorporating the earlier ideas and replacing phenomenology by explicit dependence on physico-chemical parameters [8–10]. The nonlinear development of a hydrodynamically unstable flame can be studied via the complete reactive Navier–Stokes equations, often a prohibitive task, or more systematically through much simpler models valid in the weakly nonlinear regime, such as the Michelson–Sivashinsky (MS) equation [11,12], valid for r 1 1. A first attempt to account for realistic values of thermal expansion ratios, r 1 ¼ Oð1Þ, in the fully nonlinear regime, is the recent work of Rastigejev and Matalon [13,14] in which the hydrodynamic model is implemented in a hybrid Navier– Stokes/front capturing numerical scheme. For simplicity reasons they adopted a sub-model of the hydrodynamic theory that neglects correction terms on the order of the flame thickness in the jump relations and utilizes a Markstein-type assumption for the flame speed, so that only curvature effects are retained. The main objective of the present work is to enrich this model with a flame speed relation that includes both, curvature and hydrodynamic strain, so that the complete influence of stretch on the nonlinear evolution of the flame can be investigated. In assuming a weakly disturbed flow field the MS equation effectively excludes strain rate effects and thus cannot be used to assess the full influence of stretch. At the other end of the spectrum, complex reactive DNS or even experimental studies suffer from an intrinsic difficulty in uniquely identifying the flame surface, in determining the relative hydrodynamic quantities and ultimately in establishing an interpretable a posteriori correlation between stretch and relevant parameters such as flame speed. In this respect the present study is a first attempt to discuss strain rate effects on flame propagation in a systematic and comprehensive manner.
the much larger hydrodynamic scale, the flame is an interface that separates the fresh cold mixture from the hot combustion products. The flow on either side of the flame is described by solving the incompressible Navier–Stokes equations rv¼0 q @ t v þ qðv rÞv ¼ rp þ lr2 v
ð1Þ
where v is the velocity vector, p the pressure, l the gas viscosity and q the density which takes the values qu ; qb for the unburnt/burnt gas. Consistency requires that the viscous term constitutes an OðdÞ correction. Asymptotic matching provides expressions that relate the pressure and velocities across the entire flame, and an equation for the propagation speed relative to the incoming flow, commonly referred to as the flame speed. The jump conditions take the form of Rankine–Hugoniot relations with OðdÞ corrective terms which, in the present study, will be neglected as in [13,14]. If the flame sheet is described by a function wðx; tÞ ¼ 0, such that w < 0 identifies the fresh mixture and w > 0 the burnt gas (see Fig. 1), the relations across w ¼ 0 are qðv n V f Þ ¼ 0 ð2Þ ½½n ðv nÞ ¼ 0 ð3Þ p þ qðv nÞðv n V f Þ ¼ 0 ð4Þ where ½½ is the jump operator, n ¼ rw=jrwj is the unit normal to the flame surface directed towards the burnt gas, and V f ¼ wt =jrwj is the propagation speed (in the laboratory frame) back along the normal. The flame speed S f v n V f , where v is the gas velocity just ahead of the flame ðw ¼ 0 Þ, is given by S f ¼ S L LK; K ¼ SLj þ K S ð5Þ where S L is the laminar flame speed, K is the flame stretch that includes the effects of curvature j ¼ r n and strain K S ¼ n E n with E the rate of strain tensor, and L is a coefficient on the order of the flame thickness, commonly known as the Markstein length. An explicit relation is available [15] for L in terms of the effective Lewis number of the mixture (an average of the individual Lewis numbers with a heavier weight on the deficient component) and the density ratio r.
2. Governing equations The hydrodynamic model [8,15] is based on a multi-scale analysis in which the flame, where diffusion and chemical reactions occur, is a boundary layer of thickness d ‘f =L 1. Viewed on
Fig. 1. Schematic of a flame as a gasdynamic discontinuity.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
The hydrodynamic model is a nonlinear free boundary problem, which will be addressed numerically following the methodology proposed by Rastigejev and Matalon [13]. The flame is tracked via a level set equation obtained by recasting the flame speed relation in terms of the function wðx; tÞ, namely wt þ v rw ¼ S f jrwj:
ð6Þ
The determination of v on the Lagrangian mesh representing the surface wðx; tÞ ¼ 0 requires a careful extrapolation of the velocity field that was effectively accomplished using an immersed boundary method. The ensuing w-field, suitably redistanced [16,17] to yield a signed distance function, is used to approximate the piecewise constant density field across the flame surface by means of a tanh-like function. This effectively smears the density jump between the unburnt and burnt regions over a thin numerical flame thickness, which albeit remaining hydrodynamically unresolved, retains the diffusive processes acting within it through the Markstein-like relation (5) for the flame speed. Approximating the density as a continuous function enables us to solve for the velocity field over the entire domain by discretizing the variable density Navier–Stokes equations on a uniform grid. To ensure that the jump relations (2)–(4) are satisfied, a source term has been introduced for the divergence of the velocity as discussed in [13]. A modified version of the IAMR variable density incompressible Navier–Stokes fractional step solver developed at Lawrence Berkeley National Laboratory [18] was used to solve the governing equations in a two-dimensional configuration, with appropriate conditions imposed on the boundaries. A Reynolds number Re ¼ 105 was found sufficiently accurate to simulate the required near-inviscid flows.
1089
which suggests that a propagating flame can be stretched for two reasons: its normal propagating motion ðaÞ, which can expand or contract the flame depending on the sign of the curvature, and nonuniformities in the flow giving rise to hydrodynamic strain ðb þ cÞ. There are two distinct contributions to hydrodynamic strain: a normal straining term ðbÞ, caused by the normal component of the fluid velocity in the unburnt region contracting or expanding the flame depending, again, on the sign of the local curvature and a tangential straining term ðcÞ represented by the divergence of the tangential velocity vector, which, depending on its sign, can have a compressing or expanding effect. The action of such terms is illustrated in Fig. 2. An expanding or contracting spherical flame will only be subject to the normal straining term ðbÞ, in addition to the propagation term ðaÞ, the tangential contribution being zero because the induced flow is radial. The planar flame in a stagnation-point flow will only be subject to the tangential straining term ðcÞ, the normal contribution as well as the propagation term ðaÞ being zero due to the absence of curvature. The latter two paradigmatic examples are representative of conditions that will generally coexist in a generic flame. The stretch rate experienced by a spherically expanding/converging flame is easily determined from the specific rate of flame area increase as _ K ¼ 2R=R, where RðtÞ is the flame radius and the “dot” represents differentiation with respect to time. For an outwardly propagating flame (Fig. 3, left panel), R_ ¼ rS L and j ¼ 2=R. Then, according to Eq. (7), K ¼ S L j þ ðr 1ÞS L j, where the propagation term (first term) is canceled by part of normal straining (second term) so that flame stretch K ¼ rS L j originates from residual
3. Flame stretch It is constructive to examine the constituents of the expression for flame stretch and their physical interpretation from an intrinsic reference. A curvilinear coordinate system attached to the flame surface is introduced, as in [15], with n the distance from the reaction sheet, and the local velocity vector is decomposed into components normal and tangential to the surface using the notation v ¼ vs þ vn n. The strain rate K S then takes the form vn j þ rs vs , where the gradient operator in the intrinsic coordinates is given by rs þ n @ n . The expression for flame stretch becomes
ð7Þ Fig. 2. Effects of normal and tangential straining.
1090
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
thereafter referred to as (C + S), yields a closed form dispersion relation ½r þ 1 þ kLðr 1Þx2 þ 2ð1 þ kLrÞrkS L x ½r 1 kLð3r 1Þrk 2 S 2L ¼ 0
ð10Þ
that can be contrasted to the equivalent expression based on the Markstein model [6] which contains the effect of curvature only ðr þ 1Þx2 þ 2ð1 þ kLÞrkS L x Fig. 3. (Left) Outwardly propagating spherical flame. (Right) Inwardly propagating spherical flame.
normal straining. For an inwardly propagating flame (Fig. 3, right panel), R_ ¼ S L , j ¼ 2=R and the straining term is absent. Here K ¼ S L j and stretch originates solely from the propagation term. The constituents of flame stretch can be also viewed from the Cartesian reference of Fig. 1 in which the flame sheet is described by x ¼ f ðy; tÞ. If the local velocity field is decomposed as v ¼ ðu; vÞ, then ðv nÞj ¼
vfy fyy ufyy
ð8Þ
ð1 þ fy2 Þ2
rs vs ¼ þ
vfy fyy ufyy ð1 þ
fy2 Þ2
fy2 ux ð1 þ fy2 Þ1=2
þ
þ
vy fy uy þ 1 þ fy2 1 þ fy2 fy vx
ð1 þ fy2 Þ1=2
:
ð9Þ
These expressions apply for any unfolded single valued flame profile and show that the normal straining is canceled by part of the tangential straining. The net straining effect is therefore due to the surviving tangential straining which, as we will see further on, is a pattern typical of a cusplike unstable flame front. For a planar front, such as a flame in a twodimensional stagnation-point flow parallel to the x-axis, the incoming flow is given by v ¼ ðx; yÞ, and the strain rate K S ¼ vy ¼ and is uniform along the flame surface. Similarly, for small deflections and velocity perturbations, linearization yields K S vy , so that the transverse velocity gradient is the main source of strain experienced by the flame. For a Bunsen burner flame, if the incoming flow is assumed parallel ðv ¼ 0Þ, the ensuing strain rate is K S ¼ fy uy =ð1 þ fy2 Þ. The source of straining originates here from the nonuniformity ðuy – 0Þ of the incoming flow.
ðr 1 2LrkÞrk 2 S 2L ¼ 0
ð11Þ
and is referred to as (C). Here x is the growth rate of a disturbance and k is the wavenumber. Both relations reduce to the classical Darrieus–Landau result for L ¼ 0. Note that the straining effect is a direct consequence of thermal expansion, effectively disappearing when r 1 ! 0, thus reducing Eq. (10) to Eq. (11). The solutions to the two dispersion relations, together with the Darrieus–Landau relation, are displayed in Fig. 4 for a representative L ¼ 0:01 and r ¼ 6. While the Darrieus–Landau result implies instability ðx > 0Þ to disturbances of all wavelength, the influence of diffusion through stretch provides a short wavelength stabilization, so that a limited range of the long wavelength (small wavenumber k) are unstable. It is evident from the figure that the inclusion of straining effects reduces the interval of unstable modes, and thus has an additional stabilizing effect on the short wavelength disturbances. The dispersion relation (10) provides an expression for the critical wavenumber k c or, equivalently, the critical wavelength kc 2p=k c , obtained by setting x ¼ 0; namely kc ¼ 2pð3r 1ÞL=ðr 1Þ. Only disturbances of wavelength larger than kc are amplified, so that in domains of lateral size L < kc the planar flame is unconditionally stable. Expressed in terms of the reciprocal scaled Markstein number c ðr 1ÞL=L, the condition for absolute stability is c < cc where cc ¼
4. The planar flame and its stability The sub-model under consideration, inclusive of both curvature and straining effects, here and
Fig. 4. Growth rate as a function of the wavenumber.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
2pð3r 1Þ. The equivalent expression in the absence of strain stemming from Eq. (11) is cc ¼ 4pr, which reduces to the former for r ¼ 1. When varying c, by increasing the size of the domain for example, the initially stable planar flame becomes hydrodynamically unstable at cc and a new solution evolves for c > cc . The solution beyond the bifurcation point c ¼ cc will be discussed in the next section. It can be easily verified that when expanded in powers of k, the results of the Markstein-type models Eqs. (10) and (11) take a form identical to the more rigorous asymptotic expression derived from the complete hydrodynamical model [8–10], namely
cal scheme described in Section 2. The flame is considered to propagate in a quiescent gas in a domain ðLx ; Ly Þ where Lx is large enough for the flame to attain steady conditions. Inflow/outflow conditions are specified in the streamwise direction and periodic boundary conditions in the transverse direction. Figure 5 displays the evolution to a steady state of an initially perturbed flame for increasing values of c or, equivalently, increasing domain sizes. Case (a) corresponds to c < cc and, as expected, the flame after an initial transient flattens out and propagates at a speed S L . Hydrodynamic instabilities arise when c > cc or, equivalently, when the transverse domain size Ly > kc as in cases (b) and (c). Here the flame gradually acquires a cusp-like conformation as it freely propagates in the quiescent mixture at a speed U larger than the laminar flame speed. For larger domains (larger c), the crest is increasingly sharper. The flow pattern for steady propagating conditions is shown in Fig. 6 for representative values
x ¼ x0 S L k AS L ðL Lc Þk 2 þ . . . where x0 ¼ ðr þ ðr3 þ r2 rÞ1=2 Þ=ðr þ 1Þ is the Darrieus–Landau growth rate and A ¼ AðrÞ. While the Markstein models yield Lc ¼ 0 as the critical value separating stabilizing ðL > Lc Þ from destabilizing ðL < Lc Þ diffusive-thermal effects, the exact dispersion relation implies that Lc is a function of r, which is found to be always negative in accordance to the classical linear diffusivethermal (constant density) stability analysis [19]. The expression for A obtained from the Markstein model (C + S) with straining effects included, A¼
1091
rðr þ x0 Þ ð1 þ x0 Þ ðr þ ðr þ 1Þx0 Þ
is identical to that derived from the full hydrodynamic model, whereas the last factor ð1 þ x0 Þ is missing in the original Markstein model (C) which only accounted for curvature effects. This suggests that the strain rate effects on the flame speed play a role in the onset of instability, but the corrective OðdÞ transverse convection terms appearing in the jump relations of the full hydrodynamic model – neglected in the present model – may only affect the determination of the critical Lc .
Fig. 6. Streamline pattern across a steadily propagating (downwards) cusp-like flame structure; calculated for r ¼ 4 and c ¼ 160. Bold dashed line: flame profile (level set w ¼ 0); solid lines: streamlines; in color: velocity magnitude scaled with respect to S L . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
5. Nonlinear development The nonlinear evolution of an unstable flame is analyzed based on the methodology and numeri-
a
b
c
Fig. 5. Initially perturbed flame profiles shown at consecutive instants developing in domains of size (a) Ly ¼ 1, (b) Ly ¼ 2 and (c) Ly ¼ 4; calculated for r ¼ 4 and L ¼ 0:075.
1092
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
of the parameters. A vortical flow is induced in the unburned gas by the propagating flame which gives rise to normal and tangential velocity fields similarly to the scheme illustrated in Fig. 2. As the flame propagates downwards, the expanding gas moves upwards at a significantly large speed, on the order of r. The streamlines are deflected upon crossing the flame surface and become parallel at large distances downstream. The ensuing straining field along the surface of such an unstable flame is nontrivial and depends on the given parameters, such as the expansion ratio r and the reciprocal of the scaled Markstein number c. Patterns resulting from two sets of computations, one retaining only curvature (C) and the other curvature and strain (C + S), for two values of c are displayed in Fig. 7 showing the flame shape (top panel), the velocity components in intrinsic coordinates (middle panel) and the various contributions to the strain rate K S (bottom panel). At the highly curved crest of the flame profile (where j < 0) the normal, vn j, and tangential, rs vs , straining terms have visibly opposing stretching effects with a residual compressing effect, K S < 0. At the troughs, away from the inflection points (where j > 0), the overall effect is expanding, K S > 0. The expanding effect is predominant at higher values of c, effectively “pushing” the crest upwards and resulting in a sharper “cusp” of
greater surface area compared to the profile (C) obtained when only curvature effects are considered. At lower values of c the relative compressing effect at the crest increases, resulting in a smoother, less pronounced profile of smaller surface area. This effect keeps increasing as c is lowered further, resulting in a planar profile for c < cc . The stabilization mechanism of short wavelength disturbances due to strain that was identified earlier based on linear theory, as also seen in Fig. 4, may be interpreted in a way similar to the stabilization mechanism due to curvature. The enhanced local flame speed at the crests resulting from flame compression ðK S < 0Þ and the reduced speed at the troughs resulting from flame stretching ðK S > 0Þ tend to reduce the amplitude of the perturbation, which would otherwise grow due to Darrieus–Landau instability. It was noted earlier that the flame in the nonlinear regime gradually acquires a cusp-like shape that propagates at a constant speed U without further changes in shape. At steady state the flame profile may thus be expressed as x ¼ Ut þ f ðyÞ which, upon substitution into Eq. (6) and using Eq. (5), yields U ¼ S L ð1 þ fy2 Þ1=2 LKð1 þ fy2 Þ1=2 ðu fy v Þ:
0.3
0.3
γ = 80
γ = 160 (C+S)
(C)
0.2
0.2
0
(C)
(C+S)
0.1
0.1
0
1
1
0
0
1
γ = 80
0.5
1 γ = 160
0.5
vn
0
vn
0
vτ
-0.5
-0.5
vτ -1
0
20
1
0
20
γ = 80
10
1 γ = 160
10
-vnκ
Although the individual terms on the right-hand side of this equation vary from point to point along the flame surface, their sum remains constant at steady state. The contributions to the constant propagation speed U are due to (i) area variations (first term on right-hand side), (ii) strain rate and curvature effects, through their influence on the flame speed S f (second term), and (iii) induced flow as a result of thermal expansion (last two terms). These three contributions are displayed in Fig. 8 for the representative choice c ¼ 80 with r ¼ 4 and, as seen, sum up to U 1:09. The propagation speed U was computed for various values of the expansion ratio r as a function of the bifurcation parameter c and the results are summarized in Fig. 9. The two sets of curves refer to the Markstein model (C), similarly to
-vnκ
0
0
∇τ. vτ κ
-10
Ks
Ks
-10
∇τ . vτ
-20
-20 -30
-1
ð12Þ
0
y
1
-30
κ 0
y
1
Fig. 7. (Top) Steady state flame profiles. (Middle) Normal and tangential velocity components along the flame surface. (Bottom) Contributions to the total strain rate K S ðr ¼ 4Þ.
Fig. 8. Constitutive terms of propagation speed U.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
Fig. 9. Bifurcation diagram showing the (rescaled) propagation speed as a function of parameter c. Solid red lines and filled symbols (C + S), dashed lines and empty symbols (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the findings in [14], and to the model that incorporates all the effects of stretch (C + S), which is the novel aspect of the present work. We note that for modest deformations, i.e., when the increase in surface area is insignificant but stretch effects are already active, the propagation speed U for L > 0 will be smaller than the laminar flame speed. This effect is due to the sum of the last two terms in Eq. (12) and persists until the area contribution (first term) prevails at greater flame deformations, which occurs for higher values of c. To remove this effect, which is irrelevant to the present discussion and adds unnecessary complications to the bifurcation diagram of Fig. 9, we have rescaled the propagation speed in such a way e is always positive. that the normalized value U The graph shows that the propagation speed increases abruptly as the flame spontaneously acquires a cusp-like structure and for realistic values of r reaches values exceeding above the laminar flame speed by over 20%, in accord with experimental observations [4]. Figure 9 also identifies the bifurcation point c ¼ cc , which has been determined directly by extrapolation of least square fits. The values cc obtained numerically using an extrapolation of the form ða þ bc þ cc2 Þ=ðd þ ecÞ, are displayed in Fig. 10 as a function of r together with their linear fits and the analytical relations obtained from the linear theory for the two models (C + S) inclusive of strain and (C) deprived of strain. The excellent correspondence serves as a validation of the numerical implementation of the hydrodynamic model. The results confirm that strain has an additional stabilizing effect beyond the stabilizing effects of curvature, and causes an increase in the critical value cc at a given r.
1093
Fig. 10. Critical value of the bifurcation parameter cc as a function of the thermal expansion parameter r. Dashed lines: analytical expressions for Markstein model (C) cc ¼ 4pr and present model (C + S) cc ¼ 2pð3r 1Þ.
We note that similar curves to those shown in Fig. 9 for model (C) are also obtainable [14] for small values of the expansion ratio ðr 1 1Þ using the weakly nonlinear MS equation [12,14]. For such equation, one defines U as based on contributions of area variations alone, and as a result the propagation speed curves tend to systematically plateau at high values of c. This effect occurs in correspondence to the cusp-like profiles reaching a limiting conformation, which is a behavior confirmed by the analytical pole solutions of the MS equation. Such plateauing effect, as can be seen in Fig. 9, is absent when both curvature and strain are accounted for through Eq. (12), all the more so when the expansion ratio becomes substantial ðr P 4Þ. From Fig. 9, we also note that at high postcritical values of c > cc (and r P 4), strain rate effects, when added to curvature effects, tend to yield higher propagation speeds. This is due to the prevailing expansion effect of strain and ensuing higher flame surface area as was explained in the analysis of Fig. 7. 6. Conclusions Darrieus and Landau have both concluded in their seminal papers that the instability of a flame front modeled as a surface of density discontinuity leads in itself to turbulence. Recent evidence, also supported by this work, is that the hydrodynamic instability leads to large corrugated structures, essentially controlled by the system’s domain size, that propagate steadily at a speed much larger than the laminar flame speed. These large cells may split as a result of the hydrodynamic instability, or background noise amplified by the instability, resulting in small-scale wrinkles superimposed
1094
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
on the large cells. This self-wrinkling phenomenon may eventually lead to turbulence. The focus in this work is on the nonlinear development of a hydrodynamically unstable flame, particularly the large structure that evolves beyond the instability threshold. Novel aspects of this study include accounting for the large density variations that are observed in real flames, and modeling the flame propagation with an equation that has been derived from physical first principles and accounts for both, flame front curvature and hydrodynamic strain. The effect of strain on the propagation speed, which has been neglected in previous studies, is (i) in providing an additional stabilizing effect beyond that of curvature that delays the onset of the hydrodynamic instability and (ii) acting to sharpen the amplitude of a corrugated flame near the crest, an effect that may be balanced in part by the overall influence of curvature. The evolving flame beyond the instability threshold has wide troughs and a relatively sharp, cusp-like, crest. The amplitude of the flame and the propagation speed are found to increase with increasing unburnt-to-burnt density ratio (or thermal expansion) and size of the domain (or cell size) and to decrease with increasing Markstein length. The latter depends on the mixture composition; it takes on lower values when a hydrocarbon–air mixture is made richer, or a hydrogen–air mixture is made leaner. The increase in speed amounts to over 20% above the laminar flame speed, which is a substantial contribution to the overall increase on speed observed in highly corrugated flames. The methodology of simulating a flame in the context of a hydrodynamic theory is an ideal tool for the investigation of flame–turbulence interaction. It is well-suited for the systematic investigation of turbulent flames in the flamelet regime, i.e., where no interaction exists between the smallest turbulent scales and the internal flame structure. Of great importance is the pursuit of a functional relationship adequately determining the turbulent flame speed in terms of turbulence and flame characteristics, including turbulence intensity, integral scale, heat release (or thermal expansion) and Markstein length (which incorporates mixture composition and stoichiometry). In this respect the flame sheet model permits the independent selection of all the main functional parameters, thereby allowing for extensive parametric investigation of their individual role on
the flame propagation. Such correlations may be more difficult to obtain from fully resolved DNS because this would require identifying ad-hoc iso-contours representing the flame surface from which to extract the necessary information.
Acknowledgment This work has been partially supported by the National Science Foundation under Grant CBET0943094.
References [1] G. Darrieus, Unpublished work, presented at La Technique Moderne (Paris), and in 1945 at Congre`s de Mec´anique Appliquee´ (Paris), 1938. [2] L. Landau, Acta Phys. Chim. USSR 19 (1944) 77– 85. [3] A. Palm-Leis, R.A. Strehlow, Combust. Flame 69 (1969) 111–129. [4] C. Lind, J. Whitson, Phase III Report No. CG-D85-77, Department of Transportation, US Coast Guard Final Report, 1977. [5] D. Bradley, T. Cresswell, J. Puttock, Combust. Flame 124 (2001) 551–559. [6] G. Markstein, Nonsteady Flame Propagation, MacMillan, 1964. [7] P. Clavin, F.A. Williams, J. Fluid Mech. 116 (1982) 251–282. [8] M. Matalon, B. Matkowsky, J. Fluid Mech. 124 (1984) 239–259. [9] P. Pelce, P. Clavin, J. Fluid Mech. 124 (1982) 219– 237. [10] M. Frankel, G. Sivashinsky, Combust. Sci. Technol. 29 (1982) 207–224. [11] G. Sivashinsky, Acta Astronaut. 4 (1977) 1177– 1206. [12] D. Michelson, G. Sivashinsky, Acta Astronaut. 4 (1977) 1207–1221. [13] Y. Rastigejev, M. Matalon, Combust. Theory Model. 10 (2006) 459–481. [14] Y. Rastigejev, M. Matalon, J. Fluid Mech. 554 (2006) 371–392. [15] M. Matalon, C. Cui, J. Bechtold, J. Fluid Mech. 487 (2003) 179–210. [16] J. Sethian, Level Set Methods and Fast Marching Methods, Cambridge University Press, 1996. [17] M. Sussman, E. Fatemi, SIAM J. Sci. Comput. 20 (1999) 1165–1191. [18] A. Almgren, J. Bell, P. Colella, L. Howell, M. Welcome, J. Comput. Phys. 142 (1998) 1–46. [19] G.I. Sivashinsky, Combust. Sci. Technol. 15 (1977) 137–145.
Proceedings of the
Proceedings of the Combustion Institute 33 (2011) 1087–1094
Combustion Institute www.elsevier.com/locate/proci
Strain rate effects on the nonlinear development of hydrodynamically unstable flames F. Creta *, M. Matalon Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Available online 14 August 2010
Abstract In this study we numerically implement the hydrodynamic model for a premixed flame as a nonlinear free boundary problem where the flame is tracked via a level set equation and the flow is described by a solution of the variable density Navier–Stokes equations. Unlike an earlier similar study, the present model is enriched by fully accounting for hydrodynamic strain in the flame stretch relation which, in turn, affects the local flame speed. The objective is to comprehensively analyze the effect of strain on the onset of the hydrodynamic instability and on the nonlinear development that takes place beyond its inception. The initial evolution is corroborated with the results of a linear stability analysis for which strain rate effects are fully included. We show that while strain provides an additional stabilizing effect on the short wavelength disturbances, thereby delaying the onset of the hydrodynamic instability, it acts to sharpen the cusps near the troughs of the corrugated flame that develops beyond the stability threshold resulting in a larger flame surface area and a higher propagation speed. Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Premixed flames; Thermal expansion; Hydrodynamic instability; Strain rate; Flame stretch
1. Introduction The study of flame dynamics in the past decades has been conducted utilizing a wide range of numerical and perturbative techniques. Early linear stability analyses have led Darrieus [1] and Landau [2] to establish the unconditional instability of a planar flame front propagating at constant speed. Although stable flames were observed in the laboratory owing to stabilizing influences of
* Corresponding author. Address: Dept. of Mechanical Science and Engineering, University of Illinois at Urbana Champaign, 1206 W. Green St., MC244, Urbana, IL 61801, USA. Fax: +1 217 333 1942. E-mail addresses: [email protected] (F. Creta), matalon @illinois.edu (M. Matalon).
diffusion, hydrodynamic instabilities, caused by the substantial thermal expansion, have been observed in experiments [3–5] as large scale cusp-like structures essentially controlled by the system’s domain size. Further insight on the occurrence of instabilities and on the stability limits was obtained by Markstein [6] who introduced the dependence on curvature of the flame speed through a phenomenological model. The more rigorous asymptotic treatment of multi-dimensional premixed flames exploits the multi-scale nature of the problem, which is characterized by two disparate length scales: the diffusion length ‘f that characterizes the flame thickness and the hydrodynamic length L associated with the geometrical dimensions of the vessel within which combustion takes place. Such an analysis was developed by Clavin and Williams
1540-7489/$ - see front matter Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2010.06.029
1088
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
[7] for a weakly wrinkled flame in low intensity turbulence and Oð1Þ expansion ratios r of the unburnt to burnt densities. Matalon and Matkowsky [8] eliminated the assumptions of small flame displacements and weak flow nonuniformity, and formulated a general nonlinear hydrodynamic model in which the flame is a gasdynamic discontinuity separating burnt and unburnt gases. The derived flame speed was found to depend on both curvature and hydrodynamic strain, effects that constitute the total stretch rate experienced by the flame. Stability analyses within this context provided more comprehensive results, incorporating the earlier ideas and replacing phenomenology by explicit dependence on physico-chemical parameters [8–10]. The nonlinear development of a hydrodynamically unstable flame can be studied via the complete reactive Navier–Stokes equations, often a prohibitive task, or more systematically through much simpler models valid in the weakly nonlinear regime, such as the Michelson–Sivashinsky (MS) equation [11,12], valid for r 1 1. A first attempt to account for realistic values of thermal expansion ratios, r 1 ¼ Oð1Þ, in the fully nonlinear regime, is the recent work of Rastigejev and Matalon [13,14] in which the hydrodynamic model is implemented in a hybrid Navier– Stokes/front capturing numerical scheme. For simplicity reasons they adopted a sub-model of the hydrodynamic theory that neglects correction terms on the order of the flame thickness in the jump relations and utilizes a Markstein-type assumption for the flame speed, so that only curvature effects are retained. The main objective of the present work is to enrich this model with a flame speed relation that includes both, curvature and hydrodynamic strain, so that the complete influence of stretch on the nonlinear evolution of the flame can be investigated. In assuming a weakly disturbed flow field the MS equation effectively excludes strain rate effects and thus cannot be used to assess the full influence of stretch. At the other end of the spectrum, complex reactive DNS or even experimental studies suffer from an intrinsic difficulty in uniquely identifying the flame surface, in determining the relative hydrodynamic quantities and ultimately in establishing an interpretable a posteriori correlation between stretch and relevant parameters such as flame speed. In this respect the present study is a first attempt to discuss strain rate effects on flame propagation in a systematic and comprehensive manner.
the much larger hydrodynamic scale, the flame is an interface that separates the fresh cold mixture from the hot combustion products. The flow on either side of the flame is described by solving the incompressible Navier–Stokes equations rv¼0 q @ t v þ qðv rÞv ¼ rp þ lr2 v
ð1Þ
where v is the velocity vector, p the pressure, l the gas viscosity and q the density which takes the values qu ; qb for the unburnt/burnt gas. Consistency requires that the viscous term constitutes an OðdÞ correction. Asymptotic matching provides expressions that relate the pressure and velocities across the entire flame, and an equation for the propagation speed relative to the incoming flow, commonly referred to as the flame speed. The jump conditions take the form of Rankine–Hugoniot relations with OðdÞ corrective terms which, in the present study, will be neglected as in [13,14]. If the flame sheet is described by a function wðx; tÞ ¼ 0, such that w < 0 identifies the fresh mixture and w > 0 the burnt gas (see Fig. 1), the relations across w ¼ 0 are qðv n V f Þ ¼ 0 ð2Þ ½½n ðv nÞ ¼ 0 ð3Þ p þ qðv nÞðv n V f Þ ¼ 0 ð4Þ where ½½ is the jump operator, n ¼ rw=jrwj is the unit normal to the flame surface directed towards the burnt gas, and V f ¼ wt =jrwj is the propagation speed (in the laboratory frame) back along the normal. The flame speed S f v n V f , where v is the gas velocity just ahead of the flame ðw ¼ 0 Þ, is given by S f ¼ S L LK; K ¼ SLj þ K S ð5Þ where S L is the laminar flame speed, K is the flame stretch that includes the effects of curvature j ¼ r n and strain K S ¼ n E n with E the rate of strain tensor, and L is a coefficient on the order of the flame thickness, commonly known as the Markstein length. An explicit relation is available [15] for L in terms of the effective Lewis number of the mixture (an average of the individual Lewis numbers with a heavier weight on the deficient component) and the density ratio r.
2. Governing equations The hydrodynamic model [8,15] is based on a multi-scale analysis in which the flame, where diffusion and chemical reactions occur, is a boundary layer of thickness d ‘f =L 1. Viewed on
Fig. 1. Schematic of a flame as a gasdynamic discontinuity.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
The hydrodynamic model is a nonlinear free boundary problem, which will be addressed numerically following the methodology proposed by Rastigejev and Matalon [13]. The flame is tracked via a level set equation obtained by recasting the flame speed relation in terms of the function wðx; tÞ, namely wt þ v rw ¼ S f jrwj:
ð6Þ
The determination of v on the Lagrangian mesh representing the surface wðx; tÞ ¼ 0 requires a careful extrapolation of the velocity field that was effectively accomplished using an immersed boundary method. The ensuing w-field, suitably redistanced [16,17] to yield a signed distance function, is used to approximate the piecewise constant density field across the flame surface by means of a tanh-like function. This effectively smears the density jump between the unburnt and burnt regions over a thin numerical flame thickness, which albeit remaining hydrodynamically unresolved, retains the diffusive processes acting within it through the Markstein-like relation (5) for the flame speed. Approximating the density as a continuous function enables us to solve for the velocity field over the entire domain by discretizing the variable density Navier–Stokes equations on a uniform grid. To ensure that the jump relations (2)–(4) are satisfied, a source term has been introduced for the divergence of the velocity as discussed in [13]. A modified version of the IAMR variable density incompressible Navier–Stokes fractional step solver developed at Lawrence Berkeley National Laboratory [18] was used to solve the governing equations in a two-dimensional configuration, with appropriate conditions imposed on the boundaries. A Reynolds number Re ¼ 105 was found sufficiently accurate to simulate the required near-inviscid flows.
1089
which suggests that a propagating flame can be stretched for two reasons: its normal propagating motion ðaÞ, which can expand or contract the flame depending on the sign of the curvature, and nonuniformities in the flow giving rise to hydrodynamic strain ðb þ cÞ. There are two distinct contributions to hydrodynamic strain: a normal straining term ðbÞ, caused by the normal component of the fluid velocity in the unburnt region contracting or expanding the flame depending, again, on the sign of the local curvature and a tangential straining term ðcÞ represented by the divergence of the tangential velocity vector, which, depending on its sign, can have a compressing or expanding effect. The action of such terms is illustrated in Fig. 2. An expanding or contracting spherical flame will only be subject to the normal straining term ðbÞ, in addition to the propagation term ðaÞ, the tangential contribution being zero because the induced flow is radial. The planar flame in a stagnation-point flow will only be subject to the tangential straining term ðcÞ, the normal contribution as well as the propagation term ðaÞ being zero due to the absence of curvature. The latter two paradigmatic examples are representative of conditions that will generally coexist in a generic flame. The stretch rate experienced by a spherically expanding/converging flame is easily determined from the specific rate of flame area increase as _ K ¼ 2R=R, where RðtÞ is the flame radius and the “dot” represents differentiation with respect to time. For an outwardly propagating flame (Fig. 3, left panel), R_ ¼ rS L and j ¼ 2=R. Then, according to Eq. (7), K ¼ S L j þ ðr 1ÞS L j, where the propagation term (first term) is canceled by part of normal straining (second term) so that flame stretch K ¼ rS L j originates from residual
3. Flame stretch It is constructive to examine the constituents of the expression for flame stretch and their physical interpretation from an intrinsic reference. A curvilinear coordinate system attached to the flame surface is introduced, as in [15], with n the distance from the reaction sheet, and the local velocity vector is decomposed into components normal and tangential to the surface using the notation v ¼ vs þ vn n. The strain rate K S then takes the form vn j þ rs vs , where the gradient operator in the intrinsic coordinates is given by rs þ n @ n . The expression for flame stretch becomes
ð7Þ Fig. 2. Effects of normal and tangential straining.
1090
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
thereafter referred to as (C + S), yields a closed form dispersion relation ½r þ 1 þ kLðr 1Þx2 þ 2ð1 þ kLrÞrkS L x ½r 1 kLð3r 1Þrk 2 S 2L ¼ 0
ð10Þ
that can be contrasted to the equivalent expression based on the Markstein model [6] which contains the effect of curvature only ðr þ 1Þx2 þ 2ð1 þ kLÞrkS L x Fig. 3. (Left) Outwardly propagating spherical flame. (Right) Inwardly propagating spherical flame.
normal straining. For an inwardly propagating flame (Fig. 3, right panel), R_ ¼ S L , j ¼ 2=R and the straining term is absent. Here K ¼ S L j and stretch originates solely from the propagation term. The constituents of flame stretch can be also viewed from the Cartesian reference of Fig. 1 in which the flame sheet is described by x ¼ f ðy; tÞ. If the local velocity field is decomposed as v ¼ ðu; vÞ, then ðv nÞj ¼
vfy fyy ufyy
ð8Þ
ð1 þ fy2 Þ2
rs vs ¼ þ
vfy fyy ufyy ð1 þ
fy2 Þ2
fy2 ux ð1 þ fy2 Þ1=2
þ
þ
vy fy uy þ 1 þ fy2 1 þ fy2 fy vx
ð1 þ fy2 Þ1=2
:
ð9Þ
These expressions apply for any unfolded single valued flame profile and show that the normal straining is canceled by part of the tangential straining. The net straining effect is therefore due to the surviving tangential straining which, as we will see further on, is a pattern typical of a cusplike unstable flame front. For a planar front, such as a flame in a twodimensional stagnation-point flow parallel to the x-axis, the incoming flow is given by v ¼ ðx; yÞ, and the strain rate K S ¼ vy ¼ and is uniform along the flame surface. Similarly, for small deflections and velocity perturbations, linearization yields K S vy , so that the transverse velocity gradient is the main source of strain experienced by the flame. For a Bunsen burner flame, if the incoming flow is assumed parallel ðv ¼ 0Þ, the ensuing strain rate is K S ¼ fy uy =ð1 þ fy2 Þ. The source of straining originates here from the nonuniformity ðuy – 0Þ of the incoming flow.
ðr 1 2LrkÞrk 2 S 2L ¼ 0
ð11Þ
and is referred to as (C). Here x is the growth rate of a disturbance and k is the wavenumber. Both relations reduce to the classical Darrieus–Landau result for L ¼ 0. Note that the straining effect is a direct consequence of thermal expansion, effectively disappearing when r 1 ! 0, thus reducing Eq. (10) to Eq. (11). The solutions to the two dispersion relations, together with the Darrieus–Landau relation, are displayed in Fig. 4 for a representative L ¼ 0:01 and r ¼ 6. While the Darrieus–Landau result implies instability ðx > 0Þ to disturbances of all wavelength, the influence of diffusion through stretch provides a short wavelength stabilization, so that a limited range of the long wavelength (small wavenumber k) are unstable. It is evident from the figure that the inclusion of straining effects reduces the interval of unstable modes, and thus has an additional stabilizing effect on the short wavelength disturbances. The dispersion relation (10) provides an expression for the critical wavenumber k c or, equivalently, the critical wavelength kc 2p=k c , obtained by setting x ¼ 0; namely kc ¼ 2pð3r 1ÞL=ðr 1Þ. Only disturbances of wavelength larger than kc are amplified, so that in domains of lateral size L < kc the planar flame is unconditionally stable. Expressed in terms of the reciprocal scaled Markstein number c ðr 1ÞL=L, the condition for absolute stability is c < cc where cc ¼
4. The planar flame and its stability The sub-model under consideration, inclusive of both curvature and straining effects, here and
Fig. 4. Growth rate as a function of the wavenumber.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
2pð3r 1Þ. The equivalent expression in the absence of strain stemming from Eq. (11) is cc ¼ 4pr, which reduces to the former for r ¼ 1. When varying c, by increasing the size of the domain for example, the initially stable planar flame becomes hydrodynamically unstable at cc and a new solution evolves for c > cc . The solution beyond the bifurcation point c ¼ cc will be discussed in the next section. It can be easily verified that when expanded in powers of k, the results of the Markstein-type models Eqs. (10) and (11) take a form identical to the more rigorous asymptotic expression derived from the complete hydrodynamical model [8–10], namely
cal scheme described in Section 2. The flame is considered to propagate in a quiescent gas in a domain ðLx ; Ly Þ where Lx is large enough for the flame to attain steady conditions. Inflow/outflow conditions are specified in the streamwise direction and periodic boundary conditions in the transverse direction. Figure 5 displays the evolution to a steady state of an initially perturbed flame for increasing values of c or, equivalently, increasing domain sizes. Case (a) corresponds to c < cc and, as expected, the flame after an initial transient flattens out and propagates at a speed S L . Hydrodynamic instabilities arise when c > cc or, equivalently, when the transverse domain size Ly > kc as in cases (b) and (c). Here the flame gradually acquires a cusp-like conformation as it freely propagates in the quiescent mixture at a speed U larger than the laminar flame speed. For larger domains (larger c), the crest is increasingly sharper. The flow pattern for steady propagating conditions is shown in Fig. 6 for representative values
x ¼ x0 S L k AS L ðL Lc Þk 2 þ . . . where x0 ¼ ðr þ ðr3 þ r2 rÞ1=2 Þ=ðr þ 1Þ is the Darrieus–Landau growth rate and A ¼ AðrÞ. While the Markstein models yield Lc ¼ 0 as the critical value separating stabilizing ðL > Lc Þ from destabilizing ðL < Lc Þ diffusive-thermal effects, the exact dispersion relation implies that Lc is a function of r, which is found to be always negative in accordance to the classical linear diffusivethermal (constant density) stability analysis [19]. The expression for A obtained from the Markstein model (C + S) with straining effects included, A¼
1091
rðr þ x0 Þ ð1 þ x0 Þ ðr þ ðr þ 1Þx0 Þ
is identical to that derived from the full hydrodynamic model, whereas the last factor ð1 þ x0 Þ is missing in the original Markstein model (C) which only accounted for curvature effects. This suggests that the strain rate effects on the flame speed play a role in the onset of instability, but the corrective OðdÞ transverse convection terms appearing in the jump relations of the full hydrodynamic model – neglected in the present model – may only affect the determination of the critical Lc .
Fig. 6. Streamline pattern across a steadily propagating (downwards) cusp-like flame structure; calculated for r ¼ 4 and c ¼ 160. Bold dashed line: flame profile (level set w ¼ 0); solid lines: streamlines; in color: velocity magnitude scaled with respect to S L . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
5. Nonlinear development The nonlinear evolution of an unstable flame is analyzed based on the methodology and numeri-
a
b
c
Fig. 5. Initially perturbed flame profiles shown at consecutive instants developing in domains of size (a) Ly ¼ 1, (b) Ly ¼ 2 and (c) Ly ¼ 4; calculated for r ¼ 4 and L ¼ 0:075.
1092
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
of the parameters. A vortical flow is induced in the unburned gas by the propagating flame which gives rise to normal and tangential velocity fields similarly to the scheme illustrated in Fig. 2. As the flame propagates downwards, the expanding gas moves upwards at a significantly large speed, on the order of r. The streamlines are deflected upon crossing the flame surface and become parallel at large distances downstream. The ensuing straining field along the surface of such an unstable flame is nontrivial and depends on the given parameters, such as the expansion ratio r and the reciprocal of the scaled Markstein number c. Patterns resulting from two sets of computations, one retaining only curvature (C) and the other curvature and strain (C + S), for two values of c are displayed in Fig. 7 showing the flame shape (top panel), the velocity components in intrinsic coordinates (middle panel) and the various contributions to the strain rate K S (bottom panel). At the highly curved crest of the flame profile (where j < 0) the normal, vn j, and tangential, rs vs , straining terms have visibly opposing stretching effects with a residual compressing effect, K S < 0. At the troughs, away from the inflection points (where j > 0), the overall effect is expanding, K S > 0. The expanding effect is predominant at higher values of c, effectively “pushing” the crest upwards and resulting in a sharper “cusp” of
greater surface area compared to the profile (C) obtained when only curvature effects are considered. At lower values of c the relative compressing effect at the crest increases, resulting in a smoother, less pronounced profile of smaller surface area. This effect keeps increasing as c is lowered further, resulting in a planar profile for c < cc . The stabilization mechanism of short wavelength disturbances due to strain that was identified earlier based on linear theory, as also seen in Fig. 4, may be interpreted in a way similar to the stabilization mechanism due to curvature. The enhanced local flame speed at the crests resulting from flame compression ðK S < 0Þ and the reduced speed at the troughs resulting from flame stretching ðK S > 0Þ tend to reduce the amplitude of the perturbation, which would otherwise grow due to Darrieus–Landau instability. It was noted earlier that the flame in the nonlinear regime gradually acquires a cusp-like shape that propagates at a constant speed U without further changes in shape. At steady state the flame profile may thus be expressed as x ¼ Ut þ f ðyÞ which, upon substitution into Eq. (6) and using Eq. (5), yields U ¼ S L ð1 þ fy2 Þ1=2 LKð1 þ fy2 Þ1=2 ðu fy v Þ:
0.3
0.3
γ = 80
γ = 160 (C+S)
(C)
0.2
0.2
0
(C)
(C+S)
0.1
0.1
0
1
1
0
0
1
γ = 80
0.5
1 γ = 160
0.5
vn
0
vn
0
vτ
-0.5
-0.5
vτ -1
0
20
1
0
20
γ = 80
10
1 γ = 160
10
-vnκ
Although the individual terms on the right-hand side of this equation vary from point to point along the flame surface, their sum remains constant at steady state. The contributions to the constant propagation speed U are due to (i) area variations (first term on right-hand side), (ii) strain rate and curvature effects, through their influence on the flame speed S f (second term), and (iii) induced flow as a result of thermal expansion (last two terms). These three contributions are displayed in Fig. 8 for the representative choice c ¼ 80 with r ¼ 4 and, as seen, sum up to U 1:09. The propagation speed U was computed for various values of the expansion ratio r as a function of the bifurcation parameter c and the results are summarized in Fig. 9. The two sets of curves refer to the Markstein model (C), similarly to
-vnκ
0
0
∇τ. vτ κ
-10
Ks
Ks
-10
∇τ . vτ
-20
-20 -30
-1
ð12Þ
0
y
1
-30
κ 0
y
1
Fig. 7. (Top) Steady state flame profiles. (Middle) Normal and tangential velocity components along the flame surface. (Bottom) Contributions to the total strain rate K S ðr ¼ 4Þ.
Fig. 8. Constitutive terms of propagation speed U.
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
Fig. 9. Bifurcation diagram showing the (rescaled) propagation speed as a function of parameter c. Solid red lines and filled symbols (C + S), dashed lines and empty symbols (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the findings in [14], and to the model that incorporates all the effects of stretch (C + S), which is the novel aspect of the present work. We note that for modest deformations, i.e., when the increase in surface area is insignificant but stretch effects are already active, the propagation speed U for L > 0 will be smaller than the laminar flame speed. This effect is due to the sum of the last two terms in Eq. (12) and persists until the area contribution (first term) prevails at greater flame deformations, which occurs for higher values of c. To remove this effect, which is irrelevant to the present discussion and adds unnecessary complications to the bifurcation diagram of Fig. 9, we have rescaled the propagation speed in such a way e is always positive. that the normalized value U The graph shows that the propagation speed increases abruptly as the flame spontaneously acquires a cusp-like structure and for realistic values of r reaches values exceeding above the laminar flame speed by over 20%, in accord with experimental observations [4]. Figure 9 also identifies the bifurcation point c ¼ cc , which has been determined directly by extrapolation of least square fits. The values cc obtained numerically using an extrapolation of the form ða þ bc þ cc2 Þ=ðd þ ecÞ, are displayed in Fig. 10 as a function of r together with their linear fits and the analytical relations obtained from the linear theory for the two models (C + S) inclusive of strain and (C) deprived of strain. The excellent correspondence serves as a validation of the numerical implementation of the hydrodynamic model. The results confirm that strain has an additional stabilizing effect beyond the stabilizing effects of curvature, and causes an increase in the critical value cc at a given r.
1093
Fig. 10. Critical value of the bifurcation parameter cc as a function of the thermal expansion parameter r. Dashed lines: analytical expressions for Markstein model (C) cc ¼ 4pr and present model (C + S) cc ¼ 2pð3r 1Þ.
We note that similar curves to those shown in Fig. 9 for model (C) are also obtainable [14] for small values of the expansion ratio ðr 1 1Þ using the weakly nonlinear MS equation [12,14]. For such equation, one defines U as based on contributions of area variations alone, and as a result the propagation speed curves tend to systematically plateau at high values of c. This effect occurs in correspondence to the cusp-like profiles reaching a limiting conformation, which is a behavior confirmed by the analytical pole solutions of the MS equation. Such plateauing effect, as can be seen in Fig. 9, is absent when both curvature and strain are accounted for through Eq. (12), all the more so when the expansion ratio becomes substantial ðr P 4Þ. From Fig. 9, we also note that at high postcritical values of c > cc (and r P 4), strain rate effects, when added to curvature effects, tend to yield higher propagation speeds. This is due to the prevailing expansion effect of strain and ensuing higher flame surface area as was explained in the analysis of Fig. 7. 6. Conclusions Darrieus and Landau have both concluded in their seminal papers that the instability of a flame front modeled as a surface of density discontinuity leads in itself to turbulence. Recent evidence, also supported by this work, is that the hydrodynamic instability leads to large corrugated structures, essentially controlled by the system’s domain size, that propagate steadily at a speed much larger than the laminar flame speed. These large cells may split as a result of the hydrodynamic instability, or background noise amplified by the instability, resulting in small-scale wrinkles superimposed
1094
F. Creta, M. Matalon / Proceedings of the Combustion Institute 33 (2011) 1087–1094
on the large cells. This self-wrinkling phenomenon may eventually lead to turbulence. The focus in this work is on the nonlinear development of a hydrodynamically unstable flame, particularly the large structure that evolves beyond the instability threshold. Novel aspects of this study include accounting for the large density variations that are observed in real flames, and modeling the flame propagation with an equation that has been derived from physical first principles and accounts for both, flame front curvature and hydrodynamic strain. The effect of strain on the propagation speed, which has been neglected in previous studies, is (i) in providing an additional stabilizing effect beyond that of curvature that delays the onset of the hydrodynamic instability and (ii) acting to sharpen the amplitude of a corrugated flame near the crest, an effect that may be balanced in part by the overall influence of curvature. The evolving flame beyond the instability threshold has wide troughs and a relatively sharp, cusp-like, crest. The amplitude of the flame and the propagation speed are found to increase with increasing unburnt-to-burnt density ratio (or thermal expansion) and size of the domain (or cell size) and to decrease with increasing Markstein length. The latter depends on the mixture composition; it takes on lower values when a hydrocarbon–air mixture is made richer, or a hydrogen–air mixture is made leaner. The increase in speed amounts to over 20% above the laminar flame speed, which is a substantial contribution to the overall increase on speed observed in highly corrugated flames. The methodology of simulating a flame in the context of a hydrodynamic theory is an ideal tool for the investigation of flame–turbulence interaction. It is well-suited for the systematic investigation of turbulent flames in the flamelet regime, i.e., where no interaction exists between the smallest turbulent scales and the internal flame structure. Of great importance is the pursuit of a functional relationship adequately determining the turbulent flame speed in terms of turbulence and flame characteristics, including turbulence intensity, integral scale, heat release (or thermal expansion) and Markstein length (which incorporates mixture composition and stoichiometry). In this respect the flame sheet model permits the independent selection of all the main functional parameters, thereby allowing for extensive parametric investigation of their individual role on
the flame propagation. Such correlations may be more difficult to obtain from fully resolved DNS because this would require identifying ad-hoc iso-contours representing the flame surface from which to extract the necessary information.
Acknowledgment This work has been partially supported by the National Science Foundation under Grant CBET0943094.
References [1] G. Darrieus, Unpublished work, presented at La Technique Moderne (Paris), and in 1945 at Congre`s de Mec´anique Appliquee´ (Paris), 1938. [2] L. Landau, Acta Phys. Chim. USSR 19 (1944) 77– 85. [3] A. Palm-Leis, R.A. Strehlow, Combust. Flame 69 (1969) 111–129. [4] C. Lind, J. Whitson, Phase III Report No. CG-D85-77, Department of Transportation, US Coast Guard Final Report, 1977. [5] D. Bradley, T. Cresswell, J. Puttock, Combust. Flame 124 (2001) 551–559. [6] G. Markstein, Nonsteady Flame Propagation, MacMillan, 1964. [7] P. Clavin, F.A. Williams, J. Fluid Mech. 116 (1982) 251–282. [8] M. Matalon, B. Matkowsky, J. Fluid Mech. 124 (1984) 239–259. [9] P. Pelce, P. Clavin, J. Fluid Mech. 124 (1982) 219– 237. [10] M. Frankel, G. Sivashinsky, Combust. Sci. Technol. 29 (1982) 207–224. [11] G. Sivashinsky, Acta Astronaut. 4 (1977) 1177– 1206. [12] D. Michelson, G. Sivashinsky, Acta Astronaut. 4 (1977) 1207–1221. [13] Y. Rastigejev, M. Matalon, Combust. Theory Model. 10 (2006) 459–481. [14] Y. Rastigejev, M. Matalon, J. Fluid Mech. 554 (2006) 371–392. [15] M. Matalon, C. Cui, J. Bechtold, J. Fluid Mech. 487 (2003) 179–210. [16] J. Sethian, Level Set Methods and Fast Marching Methods, Cambridge University Press, 1996. [17] M. Sussman, E. Fatemi, SIAM J. Sci. Comput. 20 (1999) 1165–1191. [18] A. Almgren, J. Bell, P. Colella, L. Howell, M. Welcome, J. Comput. Phys. 142 (1998) 1–46. [19] G.I. Sivashinsky, Combust. Sci. Technol. 15 (1977) 137–145.