- Email: [email protected]
Asymptotic flame shapes and speeds of hydrodynamically unstable laminar flames
Twenty-Seventh Symposium (International) on Combustion/The Combustion Institute, 1998/pp. 545–553
ASYMPTOTIC FLAME SHAPES AND SPEEDS OF HYDRODYNAMICALLY UNSTABLE LAMINAR FLAMES LI-ZHENG MA and JERZY CHOMIAK Department of Thermo and Fluid Dynamics Chalmers University of Technology Go¨teborg, Sweden
The self-induced baroclinic instability of flames, the Landau-Darrieus instability, is studied numerically in the nonlinear range. A level set (G equation) based approach accounting for heat-release effects is used to follow the flame response to initial perturbation and shape evolution. It is shown that the instability leads to the development of product bubbles moving into the unburned mixture and cold mixture spikes penetrating into the burned gases similar to the bubble-spike configuration produced by the nonlinear Rayleigh–Taylor instability of interfaces in a gravitational field. Independently of the initial perturbation, the flame bubble approaches asymptotically a shape close to a paraboloid. The substantial growth of the flame-surface area due to the instability increases the flame propagation speed to the asymptotic value: (1 ` 0.29!a 1 1)Sl, where a is the density ratio and Sl the laminar burning velocity. The asymptotic amplitude of the flame is approximately A 4 0.37d!a 1 1 , where d is the flame width. The burning velocity has a minor effect on the asymptotic shape of the flames. When the turbulence scale is much smaller than the size of the combustion apparatus, the results can be directly applied to turbulent flames by replacing Sl by the turbulent burning velocity.
Introduction Laminar flame fronts are known to be intrinsically unstable due to thermodiffusional and baroclinic effects [1–4]. Although they can be stabilized by gravity, stretch, and curvature effects, there almost always exist ranges of perturbations for which the instability is observed causing severe difficulties in predicting the response of the flames to forcing, for instance, by turbulence. Our goal is to focus on the nonlinear range of the hydrodynamic instability where some asymptotic shapes and speeds of general importance to practical flames appear. The number of studies of hydrodynamic instability in the nonlinear range is very limited. Forced flames were studied recently by Cambray and Joulin [5,6] and nonlinear growth rates by Denet and Haldenwang [7]. Method of Solution—The Level Set Approach To study nonlinear 3-D flames, a simple and efficient method of front tracking is needed. Mashayek and Ashgriz [8] provide a recent review of the numerical techniques used for that purpose. Here, we use a level set formulation first introduced by Markstein [1] and then developed by Williams [9], Osher and Sethian [10] and Sethian [11] to cover fronts propagating with curvature-dependent speeds. The formulation was applied to studies of zero heat-release turbulent flames by Kerstein et al. [12] and zero heat-release flames subject to periodic flows by
Ashurst and Sivashinsky [13], Aldredge [14,15], and Yu et al. [16]. In short, the method is based on solving an equation for a scalar G in the entire domain of study so that the level set G 4 0 always corresponds to the position of the flame front. The equation reads ]G ` (u¹)G ` Sl|¹G| 4 0 ]t
(1)
where u is the velocity vector and Sl the laminar flame speed. Equation 1 is a hyperbolic, partial differential equation known as a Hamilton-Jacobi equation (for Sl 4 constant) to be solved together with fluid-flow equations since u is an unknown a priori. The solution then allows us to evaluate the interface position. Obviously, all the numerical methodologies for equations of motion can be used to solve this equation that, on the one hand, allows us to evaluate the interface position as long as the velocity field has been obtained from the equation of motion. On the other hand, such an embedding of the level function does not affect the original equation system, except for changing the isentropic exponent, the speed of sound, viscosity, and density when G changes sign. An important advantage of the method is that it is applicable to an arbitrary number of space dimensions, and because it does not rely on the discrete parameterization of the front itself, it can be solved numerically on an Eulerian grid exactly the same as for fluid-flow equations. The fluid flow can be
545
546
LAMINAR PREMIXED FLAMES
described by the continuity, momentum, energy, and equation of state of an ideal gas. We also assume that all the heat release occurs at the flame interface for which we then have
1
2
1 dq qu dAf 4 Sl 11 q dt qb dv
(2)
where dAf is the area element of the flame front, dv the volume element in which the flame is present, qu the unburned mixture density, and qb the burned gas density. As we can see, the flame in our approach is a site of volumetric expansion and, thus, becomes dynamically active introducing its own pressure field. This flame feature makes our approach different from the many known applications of the G-equation–based front tracking technique. The finite-volume method [17] is applied to transform the partial differential equations to algebraic relations that link the values of the dependent variables at the nodes of the computational grid. The TDMA (Tri-Diagonal Matrix Algorithm) is then adapted to solve the resulting algebraic relations. In this method, the Navier–Stokes equations are cast into the general transport equation in Cartesian coordinates
#
kv
](qU) dv ` ]t
#
dX
IdA 4
#
dv
SUdv
(3)
where the transport variable U can be equal to U, V, and W, and I, the total flux vector, containing convection and diffusion terms, is defined as ]U Ii 4 qUiU 1 CU ]xi
(4)
CU is the exchange coefficient that is equal to viscosity in momentum equations and is equal to zero in the continuity equation. SU denotes the source per unit volume of the variable U. By using different discretization schemes in both spatial dimensions and time, such as hybrid, MUSCL, and Van Leer for convection, the central differencing for the diffusive and backward Euler implicit scheme for the time difference terms, the discretized equation can be written aPUp 4
o anbUnb ` S c U
(5)
where ap 4 ( anb 1 SUp, the coefficients anb contain contributions from convection and diffusion, and source terms SUp and SUc contain the remaining terms. The continuity equation is transformed into a pressure correction equation as follows:
#
dv
](q) dv ` ]t
# 1a
dv
dX
p
2
¹P8 dA ` Dm ˙p 4 0
(6)
where Dm ˙ p is the continuity error and ap is the diagonal coefficient in momentum equations. Equation 6 can also be cast into the form of equation 5.
The numerical strategy used for solving the flame tracking problem in our approach was described in some detail by Chomiak and Zhou [18]. Here, we would like to add a short description of the area element calculation dAf at G 4 0 necessary for equation 2 and of the reinitialization procedure of the level set function. For the area calculation, we first redistribute G from the cell centers to the vertices using a linear approximation. Then we obtain the location of G 4 0 on the cell edges by the linear interpolation of two adjacent values of G that have opposite signs. The points give us the intersection of the flame surface with the cell edges and allow the simple calculation of dAf, assuming that the flame is locally flat. The reinitialization of the level set function is necessary for the G equation to remain a distance function and to reduce the values and irregularities of G away from the interfaces that cause difficulties in long time computations. To eliminate the difficulties, we applied the Sussman et al. [19] approach and extended it to three-dimensional computations. At each time step, the following equation is solved by an iteration method to reinitialize G as a distance function G 4 Se(G0)(1 1 |¹G|)
(7)
where Se(G0) is the following function: Se(G0) 4
G0
!
G02
` e2
(8)
and where e is a small constant. Equation 7 has the property of G remaining unchanged at the interface. Therefore, the zero level of G0 and G are the same. Away from the interface, G will converge to |¹G| 4 1. In practical calculations, one iteration sweep was carried out at each time step, a global smoothing was also applied to regions not close to the flame interface, because oscillations were introduced due to the reinitialization in long computations. Reduction of the Problem The replacement of a real flame by a density interface provides very substantial savings in computational time. However, the approach requires a prescription of the laminar flame speed in terms of mixture composition, flame stretch, and curvature. The effects of flame stretch and curvature make the problem scale dependent and, thus, difficult to generalize. However, we shall note that the physical scale of importance in instability problems is the flame thickness that is typically very small in comparison with the dimensions of the combustion apparatus. This allows us to restrict the analysis to constant flame speeds, that is, scales that are much larger than flame thickness. An additional simplification is provided by the assumption that initial perturbations are isotropic, that is, the same in the two
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
547
Fig. 1. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 2 and burning speed of 7 cm/s.
directions tangent to the flame. The problem is then reduced to a solution of flame behavior in a square space segment normal to the flame separating one perturbation from the other. The use of such periodic boundaries is allowed because the interactions of different segments of the flame are weak. Thus, we have reduced the problem to a square tube case where the tube width is much larger than the flame thickness and where slip is allowed at the walls. The calculations were performed for a 5 2 5 2 30 cm tube. One end of the tube is closed, and the other, open with flame propagation toward the closed end. At the open end, a Neumann condition has been used. The computational domain is divided into 20 2 20 2 80 grid nodes. The spatial grid size in the z direction is 3.75 2 1013 m, and the time step is 2.5 2 1014 s. Results Figure 1 shows the calculation results for a hypothetical flame with the laminar burning velocity Sl 4 7 cm/s and a density ratio of 2. The instability is triggered by sinusoidal perturbations of the same length in the x and y directions with a maximum velocity at the center of the tube pointing upward. Figure 2 shows the development of the flame interface with a burning speed Sl 4 7 cm/s and a density ratio of 5 corresponding, approximately, to a lean limit methane-air flame. Figure 3 shows the development of the flame interface with a burning speed Sl 4 7 cm/s and a density ratio of 10. Figures 4 and 5 illustrate the evolution of flames having laminar
burning velocities of 20 and 40 cm/s, respectively, and a density ratio of 10. The development of the flame is faster, but the final flame surface is almost identical, thus showing that the relative growth of the flame area is not a function of the flame speed. Similarly, we found that for larger tubes, the flame area normalized by the tube cross section does not change with tube size. The flame geometry is illustrated in more detail in Fig. 6, where two views of the flame are presented, one normal to the tube and the other normal to the surface passing through the corners of the tube for the density ratio of 10 and flame speeds equal to 7, 20, and 40 cm/s. The flame shapes have very much in common. In particular, the asymptotic flame surface and, consequently, the final normalized flame speed are almost constant. A particular feature of high-speed flames is the development of the bubble–spike system familiar from the studies of upward-propagating flames in a gravity field [18] and nonlinear development of Rayleigh– Taylor instabilities [20]. It may, thus, be surprising that the flame surface does not change with the burning velocity, but we can notice that the “bubble” geometry is very similar for all the cases and the spikes, even if their length grows with the burning velocity, become, at the same time, narrower without increasing the flame-surface area. Figure 7 summarizes the calculations showing the asymptotic flame-speed dependence on density ratio. The dependence can be approximated by a very simple equation Ua 4 1 ` 0.29!a 1 1 Sl
(9)
548
LAMINAR PREMIXED FLAMES
Fig. 2. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 5 and burning speed of 7 cm/s.
Fig. 3. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 7 cm/s.
where a is the density ratio across the flame. In Fig. 7, the asymptotic amplitude of the flame perturbation is also shown. Once again, a very simple equation can describe the dependence of the flame amplitude A on the density ratio A 4 0.37d !a 1 1
(10)
where d is the tube width. An interesting problem
is to study the flame behavior for multiple superimposed perturbations. A study of this problem was performed using 2-D simulations to reduce computational time. Figure 8 shows the results for a density ratio of 5 and two wavelengths: one basic with the initial amplitude Ao 4 0.024d and one superimposed with amplitude A1 4 0.003d and length d1 4 d/6. We can see that the stretch generated by the
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
549
Fig. 4. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 20 cm/s.
Fig. 5. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 40 cm/s.
long wave suppresses short perturbations that shows that nonlinear growth is dominated by instabilities having the longest wavelength and confirms the generality of our single perturbation mode analysis. The gravity effects on the flame shapes are illustrated in Fig. 9 for a density ratio of 10 and a flame with a burning speed of 40 cm/s. Comparison of Fig.
9 with Fig. 5 shows that for the conditions considered, gravity effects are clearly visible but are not dominant. Stronger gravity effects can be expected for larger flames. To check the accuracy of the method, we compared our calculations with measurement data for downward-propagating flames. The comparison is
550
LAMINAR PREMIXED FLAMES
Fig. 6. Asymptotic flame shapes seen (a) perpendicular to the tube wall and (b) perpendicular to the surface passing through the tube corners for the density ratio 10 and burning speeds 7(1), 20(2), and 40(3) cm/s. Final flame surface area (1) 48.50 cm2, (2) 46.71 cm2, and (3) 51.73 cm2.
Fig. 7. Asymptotic flame speeds and amplitudes for zero gravity and downward-propagating flames in 1 g environment as a function of density ratio. A: asymptotic flame speed; o: numerical data, B: asymptotic downward flame speed; -.-.: calculations; *: experimental data; C: asymptotic flame amplitude; o: numerical data; D: calculated asymptotic amplitude of downward-propagating flame.
also shown in Fig. 7. Further comparison is shown in Fig. 10, where a schlieren picture of a downwardpropagating stoichiometric methane–air flame is compared with calculations. The accuracy is quite good, although determining where the flame ends is not easy from the schlieren pictures. In general, the oversimplified representation of the flame close to
Fig. 8. 2-D nonlinear Landau–Darrieus flame instability development for two initial perturbations superposed ko 4 d, Ao 4 0.024d and k1 4 d/6, A1 4 0.003d; density ratio, 5; burning speed, 7 cm/s.
the boundaries of the computational domain (walls) where the flame speed is variable is the main source of errors in our approach. Discussion Zeldovich [2,21,22] was the first to provide an estimate of the asymptotic perturbation amplitude and flame speed resulting from the hydrodynamic instability of the flame. The asymptotic amplitude is determined from the equation dA 4 XA 1 cA2 dt
(11)
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
551
Fig. 9. Asymptotic flame shape of a downward-propagating flame for density ratio of 10 and burning speed 40 cm/s.
with dA/dt → 0 as t → `, where X is the perturbation growth parameter taken from the Landau–Darrieus linear theory [1–4] and c 4 2k2Sl /p2 results from the Huygens wave interaction effect reducing the amplitude. The particular value of c in which k is the wave number (k 4 2p/k) is obtained assuming a parabolic shape of the perturbation. After substitution, the asymptotic amplitude becomes 3 2 p !a ` a 1 a 1 a k 4 a`1 p 4 f (a)k 4
Aa 4
(12)
using the above value, we can calculate the asymptotic flame propagation speed by comparing the length of the parabola arc with the amplitude Aa and wavelength k with k u8a 2 4 Sl k
3!A 1A a
a
2
`
k 4pf (a)
`
k Arsh!p2f (a)2 4pf (a)
4
(13)
For a k 1, we have Aa '
values as they should be, because the growth rate of the perturbation diminishes with the growth of the amplitude, whereas in the Zeldovich estimate, it is assumed to be constant and equal to the initial growth rate. The second difference is that a 2-D case was considered using the Zeldovich approach. It appears that our estimate of the asymptotic shapes and speeds of the unstable laminar flames has a general importance. To show this, let us consider a flame perturbed by turbulent fluctuations of scale, l, much smaller than the size of the combustion apparatus, L(l K L). The flame is now characterized by the turbulent burning velocity St, but because l K L, it can still be treated as a thin interface. Then, the interface will behave exactly as shown in our analysis, with the only difference that Sl will now be replaced by St. Thus, a self-induced deformation of the interface will appear, as before, leading to an increase in the flame propagation speed by a (1 ` 0.29!a 1 1) factor, which, consequently, characterizes the instability effects on flame propagation in a universal way.
Conclusions p !ak 4
(14)
and u8a p ' !a Sl 2
(15)
As we can see, our results are much below these
1. Asymptotic flame shapes and speeds of hydrodynamically unstable flames were studied numerically using a level set method (G equation) accounting for heat-release effects, to track the flame. 2. The accuracy of the method was tested by a comparison with experimental data and found to be good.
552
LAMINAR PREMIXED FLAMES
a
Acknowledgment This work was supported by the Commission of the European Union within the BRITE EURAM Program.
REFERENCES
b Fig. 10. Comparison of calculated flame shape of a downward-propagating stoichiometric methane–air flame with a schlieren photograph of the flame. (a) Experimental flame shape and (b) numerical flame shape.
3. The instability leads to the development of product bubbles moving into the unburned mixture and cold mixture spikes moving into the burned gases, similar to the bubble–spike configuration produced by the Rayleigh–Taylor instability of interfaces. 4. The asymptotic amplitude of the flame can be approximated by A 4 0.37d!a 1 1, where d is the tube width and a is the density ratio. 5. The asymptotic flame speed due to the instability is uF 4 (1 ` 0.29!a 1 1)Sl, where Sl is the burning velocity. 6. The burning velocity has a minor effect on the asymptotic shape of the flames. 7. The results can be directly applied to turbulent flames when the turbulence scale is much smaller than the size of the combustion apparatus.
1. G. H. Markstein, ed., Non Steady Flame Propagation, Pergamon Press, Oxford, 1964. 2. Zeldovich, Ya. B., Barenblatt, G. I., Librovich, W. B., and Makhviladze, G. M., in Mathematical Theory of Combustion and Explosions, Plenum Press, New York, 1985, chap. 6.3. 3. Clavin, P., Prog. Energy Combust. Sci. 11:1–58 (1985). 4. Williams, F. A., in Combustion Theory, Benjamin/ Cummings, Menlo Park, CA, 1985, chap. 9.5. 5. Cambray, P. and Joulin, G., in Twenty-Fourth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1992, pp. 61–67. 6. Cambray, P. and Joulin, G., Combust. Sci. Technol. 97:405–428 (1994). 7. Denet, B. and Haldenwang, P., Combust. Sci. Technol. 104:143–167 (1995). 8. Mashayek, F. and Ashgriz, N., Int. J. Num. Meth. Fluids 20:1337–1361 (1995). 9. Williams, F. A., in Mathematics of Combustion, (J. D. Buckmaster, ed.), SIAM, Philadelphia, 1985, p. 97. 10. Osher, S. and Sethian, J. A., J. Comput. Phys. 79:12– 49 (1988). 11. Sethian, J. A., J. Diff. Geo. 31:131–161 (1990). 12. Kerstein, A., Ashurst, W. T., and Williams, F. A., Phys. Rev. A 37:2728–2731 (1988). 13. Ashurst, W. T. and Sivashinsky, G. I., Combust. Sci. Technol. 80:159–164 (1991). 14. Aldredge, R. C., Combust. Flame 90:121–133 (1992). 15. Aldredge, R. C., Combust. Flame 106:23–40 (1996). 16. Yu, K. M., Sung, C. J., and Law, C. K., Combust. Flame 97:375–383 (1994). 17. Davidson, L. and Farhanieh, B., “A Finite-Volume Code Employing Collocated Variable Arrangement and Cartesian Velocity Components for Computation of Fluid Flow and Heat Transfer in Complex ThreeDimensional Geometries,” report 92/4, Thermo and Fluid Dynamics, Chalmers University of Technology, Gothenburg, 1992. 18. Chomiak, J. and Zhou, G., in Twenty-Sixth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1996, pp. 883–889. 19. Sussman, M., Smereka, P., and Osher, S., J. Comput. Phys. 114:146–159 (1994). 20. Sharp, D. H., Physica D 12:3–18 (1984). 21. Zeldovich, Ya. B., PMTF 1:102 (1960). 22. Zeldovich, Ya. B., Istratov, A. G., Kidin, N. I., and Librovich, V. B., Combust. Sci. Technol. 24:1–13 (1980).
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
553
COMMENTS Ralph C. Aldredge III, University of California, Davis, USA. 1. Does your study consider flame propagation from the open or closed end of the tube? 2. How do your predicted flame-surface-area growth rates, at early times, compare with those predicted by early investigations by Sivashinsky, [1,2] and Clavin [3]?
REFERENCES 1. Sivashinsky, G. I., Acta Astronautica 4:1177–1206 (1997).
2. Sivashinsky, G. I., Annual Review of Fluid Mech. 115:179–199 (1983). 3. Clavin, P., and Garcia, P., J. de Me´canique The´orique et Applique´e 1-1:245–236 (1983). Author’s Reply. The study considers, both in the upward and downward propagating flame cases, propagation from the open end of the tube toward a closed end. To reduce the computational time and approach the asymptotic propagation faster, we started our computations with a substantial initial flame perturbation amplitude, typically of the order of 10% of the tube width. We did not study the early development assuming that it will follow the linear theory.
ASYMPTOTIC FLAME SHAPES AND SPEEDS OF HYDRODYNAMICALLY UNSTABLE LAMINAR FLAMES LI-ZHENG MA and JERZY CHOMIAK Department of Thermo and Fluid Dynamics Chalmers University of Technology Go¨teborg, Sweden
The self-induced baroclinic instability of flames, the Landau-Darrieus instability, is studied numerically in the nonlinear range. A level set (G equation) based approach accounting for heat-release effects is used to follow the flame response to initial perturbation and shape evolution. It is shown that the instability leads to the development of product bubbles moving into the unburned mixture and cold mixture spikes penetrating into the burned gases similar to the bubble-spike configuration produced by the nonlinear Rayleigh–Taylor instability of interfaces in a gravitational field. Independently of the initial perturbation, the flame bubble approaches asymptotically a shape close to a paraboloid. The substantial growth of the flame-surface area due to the instability increases the flame propagation speed to the asymptotic value: (1 ` 0.29!a 1 1)Sl, where a is the density ratio and Sl the laminar burning velocity. The asymptotic amplitude of the flame is approximately A 4 0.37d!a 1 1 , where d is the flame width. The burning velocity has a minor effect on the asymptotic shape of the flames. When the turbulence scale is much smaller than the size of the combustion apparatus, the results can be directly applied to turbulent flames by replacing Sl by the turbulent burning velocity.
Introduction Laminar flame fronts are known to be intrinsically unstable due to thermodiffusional and baroclinic effects [1–4]. Although they can be stabilized by gravity, stretch, and curvature effects, there almost always exist ranges of perturbations for which the instability is observed causing severe difficulties in predicting the response of the flames to forcing, for instance, by turbulence. Our goal is to focus on the nonlinear range of the hydrodynamic instability where some asymptotic shapes and speeds of general importance to practical flames appear. The number of studies of hydrodynamic instability in the nonlinear range is very limited. Forced flames were studied recently by Cambray and Joulin [5,6] and nonlinear growth rates by Denet and Haldenwang [7]. Method of Solution—The Level Set Approach To study nonlinear 3-D flames, a simple and efficient method of front tracking is needed. Mashayek and Ashgriz [8] provide a recent review of the numerical techniques used for that purpose. Here, we use a level set formulation first introduced by Markstein [1] and then developed by Williams [9], Osher and Sethian [10] and Sethian [11] to cover fronts propagating with curvature-dependent speeds. The formulation was applied to studies of zero heat-release turbulent flames by Kerstein et al. [12] and zero heat-release flames subject to periodic flows by
Ashurst and Sivashinsky [13], Aldredge [14,15], and Yu et al. [16]. In short, the method is based on solving an equation for a scalar G in the entire domain of study so that the level set G 4 0 always corresponds to the position of the flame front. The equation reads ]G ` (u¹)G ` Sl|¹G| 4 0 ]t
(1)
where u is the velocity vector and Sl the laminar flame speed. Equation 1 is a hyperbolic, partial differential equation known as a Hamilton-Jacobi equation (for Sl 4 constant) to be solved together with fluid-flow equations since u is an unknown a priori. The solution then allows us to evaluate the interface position. Obviously, all the numerical methodologies for equations of motion can be used to solve this equation that, on the one hand, allows us to evaluate the interface position as long as the velocity field has been obtained from the equation of motion. On the other hand, such an embedding of the level function does not affect the original equation system, except for changing the isentropic exponent, the speed of sound, viscosity, and density when G changes sign. An important advantage of the method is that it is applicable to an arbitrary number of space dimensions, and because it does not rely on the discrete parameterization of the front itself, it can be solved numerically on an Eulerian grid exactly the same as for fluid-flow equations. The fluid flow can be
545
546
LAMINAR PREMIXED FLAMES
described by the continuity, momentum, energy, and equation of state of an ideal gas. We also assume that all the heat release occurs at the flame interface for which we then have
1
2
1 dq qu dAf 4 Sl 11 q dt qb dv
(2)
where dAf is the area element of the flame front, dv the volume element in which the flame is present, qu the unburned mixture density, and qb the burned gas density. As we can see, the flame in our approach is a site of volumetric expansion and, thus, becomes dynamically active introducing its own pressure field. This flame feature makes our approach different from the many known applications of the G-equation–based front tracking technique. The finite-volume method [17] is applied to transform the partial differential equations to algebraic relations that link the values of the dependent variables at the nodes of the computational grid. The TDMA (Tri-Diagonal Matrix Algorithm) is then adapted to solve the resulting algebraic relations. In this method, the Navier–Stokes equations are cast into the general transport equation in Cartesian coordinates
#
kv
](qU) dv ` ]t
#
dX
IdA 4
#
dv
SUdv
(3)
where the transport variable U can be equal to U, V, and W, and I, the total flux vector, containing convection and diffusion terms, is defined as ]U Ii 4 qUiU 1 CU ]xi
(4)
CU is the exchange coefficient that is equal to viscosity in momentum equations and is equal to zero in the continuity equation. SU denotes the source per unit volume of the variable U. By using different discretization schemes in both spatial dimensions and time, such as hybrid, MUSCL, and Van Leer for convection, the central differencing for the diffusive and backward Euler implicit scheme for the time difference terms, the discretized equation can be written aPUp 4
o anbUnb ` S c U
(5)
where ap 4 ( anb 1 SUp, the coefficients anb contain contributions from convection and diffusion, and source terms SUp and SUc contain the remaining terms. The continuity equation is transformed into a pressure correction equation as follows:
#
dv
](q) dv ` ]t
# 1a
dv
dX
p
2
¹P8 dA ` Dm ˙p 4 0
(6)
where Dm ˙ p is the continuity error and ap is the diagonal coefficient in momentum equations. Equation 6 can also be cast into the form of equation 5.
The numerical strategy used for solving the flame tracking problem in our approach was described in some detail by Chomiak and Zhou [18]. Here, we would like to add a short description of the area element calculation dAf at G 4 0 necessary for equation 2 and of the reinitialization procedure of the level set function. For the area calculation, we first redistribute G from the cell centers to the vertices using a linear approximation. Then we obtain the location of G 4 0 on the cell edges by the linear interpolation of two adjacent values of G that have opposite signs. The points give us the intersection of the flame surface with the cell edges and allow the simple calculation of dAf, assuming that the flame is locally flat. The reinitialization of the level set function is necessary for the G equation to remain a distance function and to reduce the values and irregularities of G away from the interfaces that cause difficulties in long time computations. To eliminate the difficulties, we applied the Sussman et al. [19] approach and extended it to three-dimensional computations. At each time step, the following equation is solved by an iteration method to reinitialize G as a distance function G 4 Se(G0)(1 1 |¹G|)
(7)
where Se(G0) is the following function: Se(G0) 4
G0
!
G02
` e2
(8)
and where e is a small constant. Equation 7 has the property of G remaining unchanged at the interface. Therefore, the zero level of G0 and G are the same. Away from the interface, G will converge to |¹G| 4 1. In practical calculations, one iteration sweep was carried out at each time step, a global smoothing was also applied to regions not close to the flame interface, because oscillations were introduced due to the reinitialization in long computations. Reduction of the Problem The replacement of a real flame by a density interface provides very substantial savings in computational time. However, the approach requires a prescription of the laminar flame speed in terms of mixture composition, flame stretch, and curvature. The effects of flame stretch and curvature make the problem scale dependent and, thus, difficult to generalize. However, we shall note that the physical scale of importance in instability problems is the flame thickness that is typically very small in comparison with the dimensions of the combustion apparatus. This allows us to restrict the analysis to constant flame speeds, that is, scales that are much larger than flame thickness. An additional simplification is provided by the assumption that initial perturbations are isotropic, that is, the same in the two
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
547
Fig. 1. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 2 and burning speed of 7 cm/s.
directions tangent to the flame. The problem is then reduced to a solution of flame behavior in a square space segment normal to the flame separating one perturbation from the other. The use of such periodic boundaries is allowed because the interactions of different segments of the flame are weak. Thus, we have reduced the problem to a square tube case where the tube width is much larger than the flame thickness and where slip is allowed at the walls. The calculations were performed for a 5 2 5 2 30 cm tube. One end of the tube is closed, and the other, open with flame propagation toward the closed end. At the open end, a Neumann condition has been used. The computational domain is divided into 20 2 20 2 80 grid nodes. The spatial grid size in the z direction is 3.75 2 1013 m, and the time step is 2.5 2 1014 s. Results Figure 1 shows the calculation results for a hypothetical flame with the laminar burning velocity Sl 4 7 cm/s and a density ratio of 2. The instability is triggered by sinusoidal perturbations of the same length in the x and y directions with a maximum velocity at the center of the tube pointing upward. Figure 2 shows the development of the flame interface with a burning speed Sl 4 7 cm/s and a density ratio of 5 corresponding, approximately, to a lean limit methane-air flame. Figure 3 shows the development of the flame interface with a burning speed Sl 4 7 cm/s and a density ratio of 10. Figures 4 and 5 illustrate the evolution of flames having laminar
burning velocities of 20 and 40 cm/s, respectively, and a density ratio of 10. The development of the flame is faster, but the final flame surface is almost identical, thus showing that the relative growth of the flame area is not a function of the flame speed. Similarly, we found that for larger tubes, the flame area normalized by the tube cross section does not change with tube size. The flame geometry is illustrated in more detail in Fig. 6, where two views of the flame are presented, one normal to the tube and the other normal to the surface passing through the corners of the tube for the density ratio of 10 and flame speeds equal to 7, 20, and 40 cm/s. The flame shapes have very much in common. In particular, the asymptotic flame surface and, consequently, the final normalized flame speed are almost constant. A particular feature of high-speed flames is the development of the bubble–spike system familiar from the studies of upward-propagating flames in a gravity field [18] and nonlinear development of Rayleigh– Taylor instabilities [20]. It may, thus, be surprising that the flame surface does not change with the burning velocity, but we can notice that the “bubble” geometry is very similar for all the cases and the spikes, even if their length grows with the burning velocity, become, at the same time, narrower without increasing the flame-surface area. Figure 7 summarizes the calculations showing the asymptotic flame-speed dependence on density ratio. The dependence can be approximated by a very simple equation Ua 4 1 ` 0.29!a 1 1 Sl
(9)
548
LAMINAR PREMIXED FLAMES
Fig. 2. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 5 and burning speed of 7 cm/s.
Fig. 3. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 7 cm/s.
where a is the density ratio across the flame. In Fig. 7, the asymptotic amplitude of the flame perturbation is also shown. Once again, a very simple equation can describe the dependence of the flame amplitude A on the density ratio A 4 0.37d !a 1 1
(10)
where d is the tube width. An interesting problem
is to study the flame behavior for multiple superimposed perturbations. A study of this problem was performed using 2-D simulations to reduce computational time. Figure 8 shows the results for a density ratio of 5 and two wavelengths: one basic with the initial amplitude Ao 4 0.024d and one superimposed with amplitude A1 4 0.003d and length d1 4 d/6. We can see that the stretch generated by the
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
549
Fig. 4. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 20 cm/s.
Fig. 5. 3-D nonlinear Landau–Darrieus flame instability development for a density ratio of 10 and burning speed of 40 cm/s.
long wave suppresses short perturbations that shows that nonlinear growth is dominated by instabilities having the longest wavelength and confirms the generality of our single perturbation mode analysis. The gravity effects on the flame shapes are illustrated in Fig. 9 for a density ratio of 10 and a flame with a burning speed of 40 cm/s. Comparison of Fig.
9 with Fig. 5 shows that for the conditions considered, gravity effects are clearly visible but are not dominant. Stronger gravity effects can be expected for larger flames. To check the accuracy of the method, we compared our calculations with measurement data for downward-propagating flames. The comparison is
550
LAMINAR PREMIXED FLAMES
Fig. 6. Asymptotic flame shapes seen (a) perpendicular to the tube wall and (b) perpendicular to the surface passing through the tube corners for the density ratio 10 and burning speeds 7(1), 20(2), and 40(3) cm/s. Final flame surface area (1) 48.50 cm2, (2) 46.71 cm2, and (3) 51.73 cm2.
Fig. 7. Asymptotic flame speeds and amplitudes for zero gravity and downward-propagating flames in 1 g environment as a function of density ratio. A: asymptotic flame speed; o: numerical data, B: asymptotic downward flame speed; -.-.: calculations; *: experimental data; C: asymptotic flame amplitude; o: numerical data; D: calculated asymptotic amplitude of downward-propagating flame.
also shown in Fig. 7. Further comparison is shown in Fig. 10, where a schlieren picture of a downwardpropagating stoichiometric methane–air flame is compared with calculations. The accuracy is quite good, although determining where the flame ends is not easy from the schlieren pictures. In general, the oversimplified representation of the flame close to
Fig. 8. 2-D nonlinear Landau–Darrieus flame instability development for two initial perturbations superposed ko 4 d, Ao 4 0.024d and k1 4 d/6, A1 4 0.003d; density ratio, 5; burning speed, 7 cm/s.
the boundaries of the computational domain (walls) where the flame speed is variable is the main source of errors in our approach. Discussion Zeldovich [2,21,22] was the first to provide an estimate of the asymptotic perturbation amplitude and flame speed resulting from the hydrodynamic instability of the flame. The asymptotic amplitude is determined from the equation dA 4 XA 1 cA2 dt
(11)
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
551
Fig. 9. Asymptotic flame shape of a downward-propagating flame for density ratio of 10 and burning speed 40 cm/s.
with dA/dt → 0 as t → `, where X is the perturbation growth parameter taken from the Landau–Darrieus linear theory [1–4] and c 4 2k2Sl /p2 results from the Huygens wave interaction effect reducing the amplitude. The particular value of c in which k is the wave number (k 4 2p/k) is obtained assuming a parabolic shape of the perturbation. After substitution, the asymptotic amplitude becomes 3 2 p !a ` a 1 a 1 a k 4 a`1 p 4 f (a)k 4
Aa 4
(12)
using the above value, we can calculate the asymptotic flame propagation speed by comparing the length of the parabola arc with the amplitude Aa and wavelength k with k u8a 2 4 Sl k
3!A 1A a
a
2
`
k 4pf (a)
`
k Arsh!p2f (a)2 4pf (a)
4
(13)
For a k 1, we have Aa '
values as they should be, because the growth rate of the perturbation diminishes with the growth of the amplitude, whereas in the Zeldovich estimate, it is assumed to be constant and equal to the initial growth rate. The second difference is that a 2-D case was considered using the Zeldovich approach. It appears that our estimate of the asymptotic shapes and speeds of the unstable laminar flames has a general importance. To show this, let us consider a flame perturbed by turbulent fluctuations of scale, l, much smaller than the size of the combustion apparatus, L(l K L). The flame is now characterized by the turbulent burning velocity St, but because l K L, it can still be treated as a thin interface. Then, the interface will behave exactly as shown in our analysis, with the only difference that Sl will now be replaced by St. Thus, a self-induced deformation of the interface will appear, as before, leading to an increase in the flame propagation speed by a (1 ` 0.29!a 1 1) factor, which, consequently, characterizes the instability effects on flame propagation in a universal way.
Conclusions p !ak 4
(14)
and u8a p ' !a Sl 2
(15)
As we can see, our results are much below these
1. Asymptotic flame shapes and speeds of hydrodynamically unstable flames were studied numerically using a level set method (G equation) accounting for heat-release effects, to track the flame. 2. The accuracy of the method was tested by a comparison with experimental data and found to be good.
552
LAMINAR PREMIXED FLAMES
a
Acknowledgment This work was supported by the Commission of the European Union within the BRITE EURAM Program.
REFERENCES
b Fig. 10. Comparison of calculated flame shape of a downward-propagating stoichiometric methane–air flame with a schlieren photograph of the flame. (a) Experimental flame shape and (b) numerical flame shape.
3. The instability leads to the development of product bubbles moving into the unburned mixture and cold mixture spikes moving into the burned gases, similar to the bubble–spike configuration produced by the Rayleigh–Taylor instability of interfaces. 4. The asymptotic amplitude of the flame can be approximated by A 4 0.37d!a 1 1, where d is the tube width and a is the density ratio. 5. The asymptotic flame speed due to the instability is uF 4 (1 ` 0.29!a 1 1)Sl, where Sl is the burning velocity. 6. The burning velocity has a minor effect on the asymptotic shape of the flames. 7. The results can be directly applied to turbulent flames when the turbulence scale is much smaller than the size of the combustion apparatus.
1. G. H. Markstein, ed., Non Steady Flame Propagation, Pergamon Press, Oxford, 1964. 2. Zeldovich, Ya. B., Barenblatt, G. I., Librovich, W. B., and Makhviladze, G. M., in Mathematical Theory of Combustion and Explosions, Plenum Press, New York, 1985, chap. 6.3. 3. Clavin, P., Prog. Energy Combust. Sci. 11:1–58 (1985). 4. Williams, F. A., in Combustion Theory, Benjamin/ Cummings, Menlo Park, CA, 1985, chap. 9.5. 5. Cambray, P. and Joulin, G., in Twenty-Fourth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1992, pp. 61–67. 6. Cambray, P. and Joulin, G., Combust. Sci. Technol. 97:405–428 (1994). 7. Denet, B. and Haldenwang, P., Combust. Sci. Technol. 104:143–167 (1995). 8. Mashayek, F. and Ashgriz, N., Int. J. Num. Meth. Fluids 20:1337–1361 (1995). 9. Williams, F. A., in Mathematics of Combustion, (J. D. Buckmaster, ed.), SIAM, Philadelphia, 1985, p. 97. 10. Osher, S. and Sethian, J. A., J. Comput. Phys. 79:12– 49 (1988). 11. Sethian, J. A., J. Diff. Geo. 31:131–161 (1990). 12. Kerstein, A., Ashurst, W. T., and Williams, F. A., Phys. Rev. A 37:2728–2731 (1988). 13. Ashurst, W. T. and Sivashinsky, G. I., Combust. Sci. Technol. 80:159–164 (1991). 14. Aldredge, R. C., Combust. Flame 90:121–133 (1992). 15. Aldredge, R. C., Combust. Flame 106:23–40 (1996). 16. Yu, K. M., Sung, C. J., and Law, C. K., Combust. Flame 97:375–383 (1994). 17. Davidson, L. and Farhanieh, B., “A Finite-Volume Code Employing Collocated Variable Arrangement and Cartesian Velocity Components for Computation of Fluid Flow and Heat Transfer in Complex ThreeDimensional Geometries,” report 92/4, Thermo and Fluid Dynamics, Chalmers University of Technology, Gothenburg, 1992. 18. Chomiak, J. and Zhou, G., in Twenty-Sixth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1996, pp. 883–889. 19. Sussman, M., Smereka, P., and Osher, S., J. Comput. Phys. 114:146–159 (1994). 20. Sharp, D. H., Physica D 12:3–18 (1984). 21. Zeldovich, Ya. B., PMTF 1:102 (1960). 22. Zeldovich, Ya. B., Istratov, A. G., Kidin, N. I., and Librovich, V. B., Combust. Sci. Technol. 24:1–13 (1980).
FLAME SPEEDS OF HYDRODYNAMICALLY UNSTABLE FLAMES
553
COMMENTS Ralph C. Aldredge III, University of California, Davis, USA. 1. Does your study consider flame propagation from the open or closed end of the tube? 2. How do your predicted flame-surface-area growth rates, at early times, compare with those predicted by early investigations by Sivashinsky, [1,2] and Clavin [3]?
REFERENCES 1. Sivashinsky, G. I., Acta Astronautica 4:1177–1206 (1997).
2. Sivashinsky, G. I., Annual Review of Fluid Mech. 115:179–199 (1983). 3. Clavin, P., and Garcia, P., J. de Me´canique The´orique et Applique´e 1-1:245–236 (1983). Author’s Reply. The study considers, both in the upward and downward propagating flame cases, propagation from the open end of the tube toward a closed end. To reduce the computational time and approach the asymptotic propagation faster, we started our computations with a substantial initial flame perturbation amplitude, typically of the order of 10% of the tube width. We did not study the early development assuming that it will follow the linear theory.