- Email: [email protected]
A numerical study of large amplitude baroclinic instabilities of flames
Twenty-Sixth Symposium (International) on Combustion/The Combustion Institute, 1996/pp. 883–889
A NUMERICAL STUDY OF LARGE AMPLITUDE BAROCLINIC INSTABILITIES OF FLAMES JERZY CHOMIAK and GANG ZHOU Chalmers University of Technology Go¨teborg, Sweden
A steady propagation of a premixed flame subject to large amplitude perturbations generated by baroclinic instabilities was studied numerically. The square flammability limit test tubes were used to model the effects of confinement. The pure Landau-Darrieus and Rayleigh-Taylor instabilities as well as the complex case of superposed instabilities were considered with different initial perturbations leading to the formation of mushroom- and tulip-shaped flames. Three dimensional simulations of Navier-Stokes equations were performed to obtain a solution of the fluid flow problem, and a Eulerian interface tracking technique was used to follow the flame, based on the level set methodology (G equation) with expansion effects included. The effects of different factors, such as gravity, density ratio, burning velocity, initial perturbations, and tube end effects were investigated. When the results of the predictions are compared with experimental data, a good agreement is obtained in spite of the disregarded heat loss and stretch effects on flames. Due to a self-adaptive nature of the flame-flow interaction reflected in the generation of extremely complex 3-D vortices, the flow behind the flames does not contain separation regions, which maximizes the flame speed. Introducing restriction on the flow structure by imposing 2-D geometry leads to much lower flame velocities. Both calculations and experiments show that spontaneous conversion of mushroom to tulip flames occurs in open vertical tubes when the flame surface is flat. The conversion is triggered by the Landau-Darrieus instability and is very sensitive to flow asymmetries. Slight inclination of the tubes prevents the conversion.
Introduction Flame propagation is strongly influenced by hydrodynamic instabilities generated by the pressuredensity interaction, often called baroclinic effects. Maxworthy [1] first attempted to include these effects for the prediction of flames and provided a solution for flame shapes and speeds in tubes, accounting for Landau-Darrieus instability [2,3] as modified by Markstein [4] to include gravity and flame curvature effects on burning velocity. Simple physical consideration shows that the flame front speed must be very sensitive to the flame shape and, thus, we are not surprised that the agreement of the solution provided by Maxworthy [1] based on the single Fourier-Bessel component of the true shape was very poor, something he himself admits. The problem was reconsidered by Zeldovich et al. [5] for the case of pure Landau-Darrieus instability. By considering the shape effects, they note that it is not the shape of the leading edge but of the trailing edge of the flame front that is important for the front velocity, since it determines the size of the flow separation zones formed in the combustion products. Although maximum and minimum speed limits were evaluated, the characteristic speed was not determined. The approach of Zeldovich et al. [5] was generalized by Pelce´ [6] by including gravity. Noting that the flow structure is extremely complex and that it is
unlikely that an analytical solution for flames propagating in tubes can be found, he proceeded to develop some integral relations involving the upstream velocity field and flame shape to be used for the determination of the front speed and the size of the flow separation stagnation zones. Although the solution is qualitatively correct in a range of Froude numbers Fr (Fr 4 gR/S2l , where g is gravity, R the tube radius, and sl the laminar flame speed), it leads to some paradoxical results as well. For instance, for negative Froude numbers (upward propagation) below 140, the only solution is a flat flame, which is clearly unphysical. Also, the upward propagation speed, in general, is much lower than the experimental values [7]. The conclusion Pelce´ reaches is that the disagreement is a result of the upstream flow structure assumed to be composed of a single Fourier mode. Again, the problem of the flow structure and resulting flame shape becomes critical. Since the analytical solutions appear to be so difficult and ineffective, progress can be expected by performing numerical simulations. Several attempts to predict the flow in vertical tubes are known from the literature. Awn [8], Awn and Spalding [9], and Spalding and Wu [10] studied a two-dimensional case. The results of Spalding and Wu [10] deal only with a square box case of minor practical importance for which no experimental data exist, and the results of Awn and Spalding [10] for the open tube case do not
883
884
LAMINAR PREMIXED FLAMES
agree well with other experiments, particularly, concerning flame shapes and speeds. Several numerical studies were also devoted to the case of flames propagating in closed tubes and, particularly, to the case of tulip flame formation which occurs in these tubes when the flames approach the closed end [11–14]. Thus, numerical studies have not contributed much to the solution of the problem, which is still open. This is quite surprising since, in addition to pure academic interest, the solution of the problem is important for flammability limit studies in which tubes are used as standard equipment [15,16], understanding of fires in enclosures, and understanding of flame acceleration including deflagration to detonation transitions, flames in nonstationary flows in general, swirl flows, and so forth. The Physics of the Stability Loss of FluidFluid Interfaces Due to Pressure Gradients To understand the problem better, let us first give a physical explanation of the stability loss of fluidfluid interfaces caused by pressure gradients. Let us use, for that purpose, the vorticity equation. The equation is fully equivalent to the Navier-Stokes equation and can be written in the form: d X 4 (X • ¹)U dt ` m¹2X 1 X(¹ • U) 1 ¹
1 2 ¹P q
(1)
where X 4 ¹ 2 U is the vorticity and d/dt the Stokes material derivative. The physical meaning of the right-hand side components of equation 1 is as follows: • (X • ¹)U 4 X grad U represent the change in vorticity resulting from the stretching of vortices; • m¹2X is the viscous dissipation term; • X(¹ • U) 4 X div U corresponds to the effects of change in density for the vortices, which rotate more slowly if their size increases by conservation of angular momentum; and • ¹1/q 2 ¹P 4 grad 1/q 2 grad p represents the vorticity generation caused by the existence of pressure gradients in a medium of varying density, often called the baroclinic effect. Equation 1 is remarkable because of the clear physical meaning of all the components and due to the lack of interference between the terms. It also shows that there is no other term except the last one that can cause a flame (or density interface) to become unstable. Thus, the baroclinic torque is the only source of the instability. Let’s take the LandauDarrieus case first, starting with some small perturbations of the flame. The perturbed flame creates a pressure gradient determined by the mean flame po-
sition, whereas the orientation of the density gradients is determined by the surface perturbation. The misalignment of the two gradients generates the baroclinic torque and, consequently, growth of the perturbation. Exactly the same happens in the RayleighTaylor instability. In this case, the direction of the pressure gradient is determined by gravity or acceleration of the frame of reference, whereas the density gradient is determined by the geometry of the interface. Thus, both gravitational (Rayleigh-Taylor) and hydrodynamic (Landau-Darrieus) flame instabilities driven by heat release are generated by the same mechanism and, consequently, are close relatives belonging to the same family of baroclinic instabilities, a fact overlooked in the literature. Although the mechanism of the instabilities is the same, the instabilities differ quantitatively due to differences in the pressure fields. Method of Solution—The G Equation Approach A crucial problem in simulating the interface propagation is the front tracking method. Many techniques were devised for that purpose in the past: vortex methods [17], boundary integral methods [18], fluid volume methods [19,20], and direct tracking methods using interface grids [21]. In nonstationary combustion problems, the vortex methods dominate [11,12,22,23]. Here, we use a level set formulation first introduced by Williams [24] and then developed by Osher and Sethian [25] and Sethian [26] to cover fronts propagating with curvature dependent speeds. The formulation was applied for studies of zero heat release turbulent flames by Kerstein et al. [27] and zero heat release flames subject to periodic flows by Ashurst and Sivashinsky [28] and Aldridge [29]. In short, the method is based on solving an equation for a scalar G, defined in all the domains of study so that the level set G 4 0 always corresponds to the position of the flame front. The equation reads [30]: ]G ` (U • ¹)G ` Sl(1 1 ej) |¹G| 4 0 ]t
(2)
where U is the velocity vector, Sl the laminar flame speed, e the Markstein constant [4], and j the local flame front curvature. Equation 2 is a hyperbolic partial differential equation known as a “HamiltonJacobi” equation (for Sl 4 constant) to be solved together with fluid flow equations since U is an unknown a priori and can be modified by the front, especially when dynamically active interfaces are concerned, such as with real flames with heat release. The solution then allows us to evaluate the interface position. To close the system, we assume that the flow is incompressible and adiabatic. In this case, the fluid flow can be described dynamically by
BAROCLINIC INSTABILITIES OF FLAMES
the continuity and Navier-Stokes equations, which for an arbitrary stationary control volume dv with boundary dX, an outer normal unit vector n at the surface element ds can be written in the integral form: ]q dv ` dv ]t
#
#
dX
(F) • nds 4
#
dv
Hdv
(3)
where q4 F4
3qqu4,
m , and H 4 3 4 3qu quu ` pI 1 s4 qgI e
with s 4 l [grad u ` (grad
u)T]
2
(4)
Numerical Strategy In this study, we use an algorithm based on the SIMPLE method in staggered grids. In this method, the Navier-Stokes equations are cast into the general transport equation in Cartesian coordinates as follows:
dv
](qU) dv ` ]t
#
dX
J • nds 4
#
dv
SUdv
aP 4
o anbUnb ` SC U
(7)
(5)
where the transport variable U can be equal to 1
o anb 1 SP
U
The coefficients anb contain contributions from convection and diffusion, and the source terms SUP and SUC contain the remaining terms. The continuity equation is transformed into a pressure correction equation as follows:
#
dv
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 a consequence of the localized heat release, the density and viscosity outside the interface are constant. The equation system (2–4) is closed and can be used for flame simulation. 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 methodology.
#
CU is the exchange coefficient and is equal to viscosity in the 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 in the time dimension, such as hybrid, MUSCL [31,32], and CHARM [33] for convection, central differencing for diffusive and the backward Euler implicit scheme for time difference terms, the final discretized equation will have the form: where
where q is the vector of the conserved dependent variables, F the total flux tensors, H the forcing term, s the stress tensor, q the density, l the dynamic viscosity coefficient, p the pressure, and u the velocity vector. The gravity is assumed to be known and has the form g 4 (0, 0, g)T. Let us further assume that all the heat release occurs at the interface, then we have
1
(continuity) or u, v, or w, and J denotes the total flux vector, containing convection and diffusion terms: ]U Ji 4 qUiU 1 CU (6) ]xi
aPUP 4
2 1 l div uI 3
Sl(1 1 ej)dAf qu me 4 11 dv qb
885
]q dv ` ]t
# 1a
dv
dX
P
2
¹P8 • nds ` Dm ˙P40
(8)
where Dm ˙ P is the continuity error and aP is the diagonal coefficient in the momentum equations. Equation 8 can also be cast into the form of equation 7. The calculation was performed for a 5 2 5 2 30 cm cubic rectangular tube which is used in standard flammability limit studies. Since the computation domain corresponds to a physical device, the wall effect cannot be avoided. We use a no-slip boundary condition at the walls. At the open end, a Neumann condition has been used. The computational domain is divided into 20 2 20 2 180 grid nodes. The spatial grid size in the Z-direction in the computation is 0.0017 m, roughly corresponding to the physical thickness of a flame, and the time step is 0.002. All the times are in seconds. The no-slip condition at the wall is not important for the results of the calculation, as the boundary layer in the tube is very thin. For the G equation determining the flame shape, we used a second-order extrapolation value at the wall which further reduced the effects of the boundary on the calculations. The calculation is also very insensitive to the grid size due to smoothing of the flame interface performed after each time step and due to insensitivity of the baroclinic terms to the actual density distribution in the interface. It can be shown that the total baroclinic torque across the interface is determined by the pressure gradient and total density difference.
886
LAMINAR PREMIXED FLAMES
Fig. 1. Interface evolution for a pure Rayleigh-Taylor instability.
Fig. 2. Flame front evolution under the influence of a Landau-Darrieus instability.
Results Figure 1 shows the calculation results for a pure Rayleigh-Taylor instability under the influence of natural gravity of 9.81 m/sec2 for a density ratio of two fluids equal to 5. The instability is triggered by a single-wave initial perturbation with a maximum velocity at the center of the tube pointing upward. The well-known features of the Rayleigh-Taylor instabilities—the formation of a bubble of light fluid rising with a constant speed and spikes of heavy fluid approaching a free-fall condition—are well reproduced in the calculation. Figure 2 shows the calculation results for a pure Landau-Darrieus instability occurring when gravity is switched off under the influence of a lean limit methane flame having Sl 4 7 cm/s and producing a density ratio of 5. The instability is triggered by the same initial wave as is the case in Fig. 1. The insta-
Fig. 3. Flame front evolution under the influence of combined Rayleigh-Taylor and Landau-Darrieus instabilities.
Fig. 4. “Tulip” flame evolution due to combined Rayleigh-Taylor and Landau-Darrieus instabilities.
bility develops slower but the final flame shape is similar, with a constant speed of hot bubble front and cold spike. Figure 3 shows the calculation results for the instabilities superimposed. The speed of the bubble of 20 cm/s, composed of the speed of 13 cm/s for the Rayleigh-Taylor case and 7 cm/s due to flame propagation, compares well with the measured data of Jarosinski [16], Pelce´-Savornin et al. [6], and von Lavante and Strehlow [34]. It is interesting that the Landau-Darrieus instability does not contribute at all to the speed. The shape of the bubble was compared with experimental data and the achieved accuracy was very good. Figure 4 shows the flame shape for the same conditions as for Fig. 3 but with an initial perturbation wave having maximum upward speeds in the vicinity of the walls. A tulip flame develops in this case and propagates along the whole length of the tube. The speed of the tulip flame is slightly lower than for the previous “mushroom”
BAROCLINIC INSTABILITIES OF FLAMES
887
Fig. 5. Comparison of the streamlines of 3-D and 2-D flames. Y denotes the position of the longitudinal cross-section in which the streamlines are shown. The upper, gray part of the figure indicates the fresh mixture.
Fig. 6. Conversion of a flat stoichiometric methane “mushroom” flame starting from a closed end of the tube to a “tulip” flame, due to a Landau-Darrieus instability.
flames, due to slightly larger wall friction effects. Although the flame shapes look simple, the flame-induced flow is extremely complex. This complexity is illustrated by the streamlines in the different cross sections of the tube for the tulip case shown in Fig. 5. The figure also presents the 2-D calculations for the same conditions. It is interesting to note that the asymptotic 3-D bubble velocity, equal to 19.5 cm/s, is substantially higher than the 2-D bubble velocity of 16 cm/s. The flames having higher laminar burning velocities are much flatter and are prone to secondary instabilities with heavy fluid spikes penetrating through the top of the front. This instability can occur spontaneously, as shown in Fig. 6, or can be induced by the flame approaching the end of the tube. A mushroom flame having a laminar flame speed of
40 cm/s and density ratio of 7, starting from the closed end of the tube, is shown in the figure. As it approaches the open end, it converts to a tulip flame similar to those in the cases observed in the closed tubes [11–14]. The case was investigated experimentally and the results are shown in Fig. 7. As can be seen, the appearance of tulip flames in this case is confirmed by the experiment. The tulip flame formation is very sensitive to flow asymmetry—a slight inclination of the tube (158 off vertical) prevents it from appearing in both the experiment and the computations. Further studies of this effect are planned. In the computations, tulip shapes are not formed for zero heat release flames, which indicates that in flames the phenomenon is caused by the LandauDarrieus instability. A similar kind of instability for low burning velocity flames is also observed in wider tubes. In Fig. 8, simulations of multiwave initial perturbations are shown in a 10-cm wide tube with laminar flame speed of 7 cm/s and density ratio of 5. Whereas multiwave initial perturbation in a 5-cm tube leads always to single bubbles of identical shape, as in the single wave case, the 10-cm tube case produces unstable double bubbles in perfect qualitative agreement with experiments [34]. It is the multiple bubble formation that led to the selection of the 5-cm tube for flammability limit studies [15,16]. We may conclude that the simulations reproduce the physical behavior of flames with amazing correctness considering the crude nature of our assumptions, particularly concerning heat loss, flame thickness, and flame stretch effects. An important problem resolved by the calculation is the problem of stagnation zones in burned gases. The classical results [1,5,6] assumed that flow separation zones will occur in burned gases at the place where the flames meet the walls or behind flame segments concave to the mixture. These zones will substantially
888
LAMINAR PREMIXED FLAMES
reduce the global speed of the flame fronts. However, although stagnation zones were observed for 2-D flames stabilized in the tubes by external means [35], they were not observed in freely propagating flames [7,16,34]. The explanation is provided in Fig. 6. The formation of the separation zones behind the central concave fragment is suppressed by the strong three-dimensional side vortices. Thus, a freely moving flame creates a flow field that suppresses flow separation and, consequently, maximizes the flame propagation speed.
Conclusions
Fig. 7. Multiexposure photograph of the flame in a vertical tube with a time interval of 2 ms illustrating the tulip flame formation close to the open end of the tube. Conditions as in Fig. 6. In long vertical tubes, similar flame shape changes were observed at different locations, also very far from the open end.
1. A level set method based on the G equation was developed to simulate flame propagation in tubes for flames with heat release. 2. The flame geometry results from the nonlinear development of the Rayleigh-Taylor and LandauDarrieus instabilities of flames, both driven by the baroclinic effects. 3. The observed flame front speeds and shapes are close to experimental values within the range of the conditions tested. 4. 2-D and 3-D flame simulations lead to different speeds with substantially higher speeds for the 3D case. 5. Fast mushroom flames started at the closed end of vertical tube, and flames in general spontaneously formed “tulip” shapes. The effect is generated by the Landau-Darrieus instability and is very sensitive to flow asymmetry. A slight inclination of the tube (158 off vertical) prevents the conversion. 6. Whereas flames in small tubes (5 2 5 cm) form single bubbles even for multiwave initial perturbations, flames in larger tubes (10 2 10 cm) form multiple, nonstationary bubbles, in accordance with experiments. 7. Due to the self-adaptive nature of the flames reflected in the generation of complex 3-D vortices, the flow behind the flames does not contain separation regions that maximize the flame speed. Acknowledgments This work was supported by a grant from the Swedish Board for Technological Development (NUTEK).
REFERENCES
Fig. 8. Multiwave initial perturbation development in a 10 2 10 cm tube.
1. 2. 3. 4.
Maxworthy, T., Phys. Fluids 5:407–417 (1962). Landau, L. D., Acta Phys-Chem. URSS 19:77 (1944). Darrieus, G., Z. Phys. Chem. 88:641 (1944). Markstein, G., J. Aero. Sci. 18:199–209 (1951).
BAROCLINIC INSTABILITIES OF FLAMES 5. Zeldovich, Y. B. Istratov, A. G. Kidin, N. I., and Librovich, V. B., Comb. Sci. Technol. 24:1–13 (1980). 6. Pelce´, P., J. Phys. 46:503–510 (1985). 7. Pelce´-Savornin C., Quinard J., and Searby G., Combust. Sci. Technol. 58:337–346 (1988). 8. Awn, A. G., Ph.D. Thesis, London University, London 1979. 9. Awn, A. G., and Spalding, D. B., Second Conference of Mechanical Power Engineering, Air Shames University, Cairo, 1978. 10. Spalding, D. B., and Wu, J. Z. Y., PCH Phys. Chem. Hydrodynam. 7:353–384 (1986). 11. Dunn-Rankin, D., Barr, P. K., and Sawyer, R. F., Twenty First Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1986, p. 1291. 12. Rotman, D. A., and Oppenheim, A. K., Twenty First Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1986, pp. 1291–1303. 13. Gonzalez, M. Borghi, R., and Saouab, A. N., Combust. Flame 82:201–220 (1992). 14. Lee, S. T., and Tsai, Ch. H., Combust. Flame 99:484– 490 (1994). 15. Coward, H. F. and Jones, G. W., Bureau of Mines Bulletin 503, Washington, D.C., 1952. 16. Jarosinski, J., Prog. Energy Combust. Sci. 12:81–116 (1986). 17. Baker, G. R., Merion, D. I., and Orszag, S. A., J. Fluid Mech. 23:477–501 (1982). 18. Baker, G. R., and Moore, D. W., Phys. Fluids A 1:1451–1459 (1989). 19. Hirt, C. W., and Nichols, B. D., J. Comput. Phys. 39:201–225 (1981).
889
20. Fatemi, E., and Odeh, F., J. Comput. Phys. 108:209– 217 (1993). 21. Unverdi, S. O., and Tryggvason, G., J. Comput. Phys. 100:25–37 (1992). 22. Fines, A., Tsuruda, T., and Hirano, T., Combust. Flame 95:76–86 (1993). 23. Batley, G. A., McIntosh, A. C., Brindley, J., and Falle, S. A. E. G., J. Fluid Mech. 279:217–237 (1994). 24. Williams, F. A., in The Mathematics of Combustion, (J. D. Buckmaster, Ed.), SIAM, Philadelphia, 1985, p. 97. 25. Osher, S., and Sethian, J. A., J. Comput. Phys. 79:12– 49 (1988). 26. Sethian, J. A., J. Diff. Geometry 31:131–161 (1990). 27. Kerstein, A. R., Ashurst, W. T., and Williams, F. A., Phys. Rev. A 37:2728–2731 (1988). 28. Ashurst, W. T., and Sivashinsky, G. I., Combust. Sci. Technol. 80:159–164 (1991). 29. Aldridge, R. C.; Combust. Flame 90:121–133 (1992). 30. Zhu, J., and Sethian, J., J. Comput. Phys. 102:128–138 (1992). 31. van Leer, B., J. Comput. Phys. 14:361–370 (1974). 32. Lien, F. S. and Leschziner, M. A., Proceedings of the 5th International IAHR Symposium on Refind Flow Modeling and Turbulence Measurements, Paris, Sept. 1993. 33. Zhou, G., Davidson, L. and Olsson, E., Proceedings of the 14th ICNMFD, Lecture Notes in Physics, (S. M. Deshpande, S. S. Desai, and R. Narsimha, Eds.), Springer-Verlag, Berlin, 1995, No. 453, p. 372–377. 34. Von Lavante, E., and Strehlow, R. A., Combust. Flame 59:123–140 (1983). 35. Uberoi, M. S., Phys. Fluids 2:72–78 (1959).
COMMENTS Luc Bauwens, The University of Calgary, Canada. Did you have the opportunity to look at the effect of the Markstein length, especially when it becomes small? Can you say a few words on how the source of volume and vorticity at the front was implemented numerically in your code? (This strikes me as a very difficult problem.) Author’s Reply. We have not studied the effects of the Markstein length in the small value range. Our method does not allow us to resolve the flame structure and thus is unsuitable for that purpose. Furthermore, to avoid strong grid sensitivity, we had to smooth locally the flame interface after each time step using a sinusoidal smoothing function. To be able to evaluate the effects of volumetric expansion, we have to know the area element of flame front in each control volume. The computation of the area element is done using the following sequence. First we redistribute the G from the cell center to the cell vertexes. Then, we obtain G 4 0 by interpolation between two adjacent values of G which have opposite signs. After getting the data for G 4 0, we can reconstruct the local flame geometry. Then
the expansion source term is defined using Eq. (4) from the paper. The addition of the expansion source closes the hydrodynamic problem and there is no need to introduce vorticity at the interface since it is generated automatically by the pressure-density interaction terms in the N.S. equations. The baroclinic torque is insensitive to the actual density distribution in the flame and is determined by the pressure gradient and the total density difference across the flame. Practically, to avoid unphysical oscillations at the interface, we used the following equation for for q(w) in the flame zone:
5
qu qb pw q¯ ` Dqsin 2a
1 2
if w . a if w , a otherwise
(1)
where q¯ 4 (qu ` qb)/2 Dq 4 (qu 1 qb)/2 and a is a finite thickness of an order of 0(h), (where h is the mesh size).
A NUMERICAL STUDY OF LARGE AMPLITUDE BAROCLINIC INSTABILITIES OF FLAMES JERZY CHOMIAK and GANG ZHOU Chalmers University of Technology Go¨teborg, Sweden
A steady propagation of a premixed flame subject to large amplitude perturbations generated by baroclinic instabilities was studied numerically. The square flammability limit test tubes were used to model the effects of confinement. The pure Landau-Darrieus and Rayleigh-Taylor instabilities as well as the complex case of superposed instabilities were considered with different initial perturbations leading to the formation of mushroom- and tulip-shaped flames. Three dimensional simulations of Navier-Stokes equations were performed to obtain a solution of the fluid flow problem, and a Eulerian interface tracking technique was used to follow the flame, based on the level set methodology (G equation) with expansion effects included. The effects of different factors, such as gravity, density ratio, burning velocity, initial perturbations, and tube end effects were investigated. When the results of the predictions are compared with experimental data, a good agreement is obtained in spite of the disregarded heat loss and stretch effects on flames. Due to a self-adaptive nature of the flame-flow interaction reflected in the generation of extremely complex 3-D vortices, the flow behind the flames does not contain separation regions, which maximizes the flame speed. Introducing restriction on the flow structure by imposing 2-D geometry leads to much lower flame velocities. Both calculations and experiments show that spontaneous conversion of mushroom to tulip flames occurs in open vertical tubes when the flame surface is flat. The conversion is triggered by the Landau-Darrieus instability and is very sensitive to flow asymmetries. Slight inclination of the tubes prevents the conversion.
Introduction Flame propagation is strongly influenced by hydrodynamic instabilities generated by the pressuredensity interaction, often called baroclinic effects. Maxworthy [1] first attempted to include these effects for the prediction of flames and provided a solution for flame shapes and speeds in tubes, accounting for Landau-Darrieus instability [2,3] as modified by Markstein [4] to include gravity and flame curvature effects on burning velocity. Simple physical consideration shows that the flame front speed must be very sensitive to the flame shape and, thus, we are not surprised that the agreement of the solution provided by Maxworthy [1] based on the single Fourier-Bessel component of the true shape was very poor, something he himself admits. The problem was reconsidered by Zeldovich et al. [5] for the case of pure Landau-Darrieus instability. By considering the shape effects, they note that it is not the shape of the leading edge but of the trailing edge of the flame front that is important for the front velocity, since it determines the size of the flow separation zones formed in the combustion products. Although maximum and minimum speed limits were evaluated, the characteristic speed was not determined. The approach of Zeldovich et al. [5] was generalized by Pelce´ [6] by including gravity. Noting that the flow structure is extremely complex and that it is
unlikely that an analytical solution for flames propagating in tubes can be found, he proceeded to develop some integral relations involving the upstream velocity field and flame shape to be used for the determination of the front speed and the size of the flow separation stagnation zones. Although the solution is qualitatively correct in a range of Froude numbers Fr (Fr 4 gR/S2l , where g is gravity, R the tube radius, and sl the laminar flame speed), it leads to some paradoxical results as well. For instance, for negative Froude numbers (upward propagation) below 140, the only solution is a flat flame, which is clearly unphysical. Also, the upward propagation speed, in general, is much lower than the experimental values [7]. The conclusion Pelce´ reaches is that the disagreement is a result of the upstream flow structure assumed to be composed of a single Fourier mode. Again, the problem of the flow structure and resulting flame shape becomes critical. Since the analytical solutions appear to be so difficult and ineffective, progress can be expected by performing numerical simulations. Several attempts to predict the flow in vertical tubes are known from the literature. Awn [8], Awn and Spalding [9], and Spalding and Wu [10] studied a two-dimensional case. The results of Spalding and Wu [10] deal only with a square box case of minor practical importance for which no experimental data exist, and the results of Awn and Spalding [10] for the open tube case do not
883
884
LAMINAR PREMIXED FLAMES
agree well with other experiments, particularly, concerning flame shapes and speeds. Several numerical studies were also devoted to the case of flames propagating in closed tubes and, particularly, to the case of tulip flame formation which occurs in these tubes when the flames approach the closed end [11–14]. Thus, numerical studies have not contributed much to the solution of the problem, which is still open. This is quite surprising since, in addition to pure academic interest, the solution of the problem is important for flammability limit studies in which tubes are used as standard equipment [15,16], understanding of fires in enclosures, and understanding of flame acceleration including deflagration to detonation transitions, flames in nonstationary flows in general, swirl flows, and so forth. The Physics of the Stability Loss of FluidFluid Interfaces Due to Pressure Gradients To understand the problem better, let us first give a physical explanation of the stability loss of fluidfluid interfaces caused by pressure gradients. Let us use, for that purpose, the vorticity equation. The equation is fully equivalent to the Navier-Stokes equation and can be written in the form: d X 4 (X • ¹)U dt ` m¹2X 1 X(¹ • U) 1 ¹
1 2 ¹P q
(1)
where X 4 ¹ 2 U is the vorticity and d/dt the Stokes material derivative. The physical meaning of the right-hand side components of equation 1 is as follows: • (X • ¹)U 4 X grad U represent the change in vorticity resulting from the stretching of vortices; • m¹2X is the viscous dissipation term; • X(¹ • U) 4 X div U corresponds to the effects of change in density for the vortices, which rotate more slowly if their size increases by conservation of angular momentum; and • ¹1/q 2 ¹P 4 grad 1/q 2 grad p represents the vorticity generation caused by the existence of pressure gradients in a medium of varying density, often called the baroclinic effect. Equation 1 is remarkable because of the clear physical meaning of all the components and due to the lack of interference between the terms. It also shows that there is no other term except the last one that can cause a flame (or density interface) to become unstable. Thus, the baroclinic torque is the only source of the instability. Let’s take the LandauDarrieus case first, starting with some small perturbations of the flame. The perturbed flame creates a pressure gradient determined by the mean flame po-
sition, whereas the orientation of the density gradients is determined by the surface perturbation. The misalignment of the two gradients generates the baroclinic torque and, consequently, growth of the perturbation. Exactly the same happens in the RayleighTaylor instability. In this case, the direction of the pressure gradient is determined by gravity or acceleration of the frame of reference, whereas the density gradient is determined by the geometry of the interface. Thus, both gravitational (Rayleigh-Taylor) and hydrodynamic (Landau-Darrieus) flame instabilities driven by heat release are generated by the same mechanism and, consequently, are close relatives belonging to the same family of baroclinic instabilities, a fact overlooked in the literature. Although the mechanism of the instabilities is the same, the instabilities differ quantitatively due to differences in the pressure fields. Method of Solution—The G Equation Approach A crucial problem in simulating the interface propagation is the front tracking method. Many techniques were devised for that purpose in the past: vortex methods [17], boundary integral methods [18], fluid volume methods [19,20], and direct tracking methods using interface grids [21]. In nonstationary combustion problems, the vortex methods dominate [11,12,22,23]. Here, we use a level set formulation first introduced by Williams [24] and then developed by Osher and Sethian [25] and Sethian [26] to cover fronts propagating with curvature dependent speeds. The formulation was applied for studies of zero heat release turbulent flames by Kerstein et al. [27] and zero heat release flames subject to periodic flows by Ashurst and Sivashinsky [28] and Aldridge [29]. In short, the method is based on solving an equation for a scalar G, defined in all the domains of study so that the level set G 4 0 always corresponds to the position of the flame front. The equation reads [30]: ]G ` (U • ¹)G ` Sl(1 1 ej) |¹G| 4 0 ]t
(2)
where U is the velocity vector, Sl the laminar flame speed, e the Markstein constant [4], and j the local flame front curvature. Equation 2 is a hyperbolic partial differential equation known as a “HamiltonJacobi” equation (for Sl 4 constant) to be solved together with fluid flow equations since U is an unknown a priori and can be modified by the front, especially when dynamically active interfaces are concerned, such as with real flames with heat release. The solution then allows us to evaluate the interface position. To close the system, we assume that the flow is incompressible and adiabatic. In this case, the fluid flow can be described dynamically by
BAROCLINIC INSTABILITIES OF FLAMES
the continuity and Navier-Stokes equations, which for an arbitrary stationary control volume dv with boundary dX, an outer normal unit vector n at the surface element ds can be written in the integral form: ]q dv ` dv ]t
#
#
dX
(F) • nds 4
#
dv
Hdv
(3)
where q4 F4
3qqu4,
m , and H 4 3 4 3qu quu ` pI 1 s4 qgI e
with s 4 l [grad u ` (grad
u)T]
2
(4)
Numerical Strategy In this study, we use an algorithm based on the SIMPLE method in staggered grids. In this method, the Navier-Stokes equations are cast into the general transport equation in Cartesian coordinates as follows:
dv
](qU) dv ` ]t
#
dX
J • nds 4
#
dv
SUdv
aP 4
o anbUnb ` SC U
(7)
(5)
where the transport variable U can be equal to 1
o anb 1 SP
U
The coefficients anb contain contributions from convection and diffusion, and the source terms SUP and SUC contain the remaining terms. The continuity equation is transformed into a pressure correction equation as follows:
#
dv
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 a consequence of the localized heat release, the density and viscosity outside the interface are constant. The equation system (2–4) is closed and can be used for flame simulation. 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 methodology.
#
CU is the exchange coefficient and is equal to viscosity in the 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 in the time dimension, such as hybrid, MUSCL [31,32], and CHARM [33] for convection, central differencing for diffusive and the backward Euler implicit scheme for time difference terms, the final discretized equation will have the form: where
where q is the vector of the conserved dependent variables, F the total flux tensors, H the forcing term, s the stress tensor, q the density, l the dynamic viscosity coefficient, p the pressure, and u the velocity vector. The gravity is assumed to be known and has the form g 4 (0, 0, g)T. Let us further assume that all the heat release occurs at the interface, then we have
1
(continuity) or u, v, or w, and J denotes the total flux vector, containing convection and diffusion terms: ]U Ji 4 qUiU 1 CU (6) ]xi
aPUP 4
2 1 l div uI 3
Sl(1 1 ej)dAf qu me 4 11 dv qb
885
]q dv ` ]t
# 1a
dv
dX
P
2
¹P8 • nds ` Dm ˙P40
(8)
where Dm ˙ P is the continuity error and aP is the diagonal coefficient in the momentum equations. Equation 8 can also be cast into the form of equation 7. The calculation was performed for a 5 2 5 2 30 cm cubic rectangular tube which is used in standard flammability limit studies. Since the computation domain corresponds to a physical device, the wall effect cannot be avoided. We use a no-slip boundary condition at the walls. At the open end, a Neumann condition has been used. The computational domain is divided into 20 2 20 2 180 grid nodes. The spatial grid size in the Z-direction in the computation is 0.0017 m, roughly corresponding to the physical thickness of a flame, and the time step is 0.002. All the times are in seconds. The no-slip condition at the wall is not important for the results of the calculation, as the boundary layer in the tube is very thin. For the G equation determining the flame shape, we used a second-order extrapolation value at the wall which further reduced the effects of the boundary on the calculations. The calculation is also very insensitive to the grid size due to smoothing of the flame interface performed after each time step and due to insensitivity of the baroclinic terms to the actual density distribution in the interface. It can be shown that the total baroclinic torque across the interface is determined by the pressure gradient and total density difference.
886
LAMINAR PREMIXED FLAMES
Fig. 1. Interface evolution for a pure Rayleigh-Taylor instability.
Fig. 2. Flame front evolution under the influence of a Landau-Darrieus instability.
Results Figure 1 shows the calculation results for a pure Rayleigh-Taylor instability under the influence of natural gravity of 9.81 m/sec2 for a density ratio of two fluids equal to 5. The instability is triggered by a single-wave initial perturbation with a maximum velocity at the center of the tube pointing upward. The well-known features of the Rayleigh-Taylor instabilities—the formation of a bubble of light fluid rising with a constant speed and spikes of heavy fluid approaching a free-fall condition—are well reproduced in the calculation. Figure 2 shows the calculation results for a pure Landau-Darrieus instability occurring when gravity is switched off under the influence of a lean limit methane flame having Sl 4 7 cm/s and producing a density ratio of 5. The instability is triggered by the same initial wave as is the case in Fig. 1. The insta-
Fig. 3. Flame front evolution under the influence of combined Rayleigh-Taylor and Landau-Darrieus instabilities.
Fig. 4. “Tulip” flame evolution due to combined Rayleigh-Taylor and Landau-Darrieus instabilities.
bility develops slower but the final flame shape is similar, with a constant speed of hot bubble front and cold spike. Figure 3 shows the calculation results for the instabilities superimposed. The speed of the bubble of 20 cm/s, composed of the speed of 13 cm/s for the Rayleigh-Taylor case and 7 cm/s due to flame propagation, compares well with the measured data of Jarosinski [16], Pelce´-Savornin et al. [6], and von Lavante and Strehlow [34]. It is interesting that the Landau-Darrieus instability does not contribute at all to the speed. The shape of the bubble was compared with experimental data and the achieved accuracy was very good. Figure 4 shows the flame shape for the same conditions as for Fig. 3 but with an initial perturbation wave having maximum upward speeds in the vicinity of the walls. A tulip flame develops in this case and propagates along the whole length of the tube. The speed of the tulip flame is slightly lower than for the previous “mushroom”
BAROCLINIC INSTABILITIES OF FLAMES
887
Fig. 5. Comparison of the streamlines of 3-D and 2-D flames. Y denotes the position of the longitudinal cross-section in which the streamlines are shown. The upper, gray part of the figure indicates the fresh mixture.
Fig. 6. Conversion of a flat stoichiometric methane “mushroom” flame starting from a closed end of the tube to a “tulip” flame, due to a Landau-Darrieus instability.
flames, due to slightly larger wall friction effects. Although the flame shapes look simple, the flame-induced flow is extremely complex. This complexity is illustrated by the streamlines in the different cross sections of the tube for the tulip case shown in Fig. 5. The figure also presents the 2-D calculations for the same conditions. It is interesting to note that the asymptotic 3-D bubble velocity, equal to 19.5 cm/s, is substantially higher than the 2-D bubble velocity of 16 cm/s. The flames having higher laminar burning velocities are much flatter and are prone to secondary instabilities with heavy fluid spikes penetrating through the top of the front. This instability can occur spontaneously, as shown in Fig. 6, or can be induced by the flame approaching the end of the tube. A mushroom flame having a laminar flame speed of
40 cm/s and density ratio of 7, starting from the closed end of the tube, is shown in the figure. As it approaches the open end, it converts to a tulip flame similar to those in the cases observed in the closed tubes [11–14]. The case was investigated experimentally and the results are shown in Fig. 7. As can be seen, the appearance of tulip flames in this case is confirmed by the experiment. The tulip flame formation is very sensitive to flow asymmetry—a slight inclination of the tube (158 off vertical) prevents it from appearing in both the experiment and the computations. Further studies of this effect are planned. In the computations, tulip shapes are not formed for zero heat release flames, which indicates that in flames the phenomenon is caused by the LandauDarrieus instability. A similar kind of instability for low burning velocity flames is also observed in wider tubes. In Fig. 8, simulations of multiwave initial perturbations are shown in a 10-cm wide tube with laminar flame speed of 7 cm/s and density ratio of 5. Whereas multiwave initial perturbation in a 5-cm tube leads always to single bubbles of identical shape, as in the single wave case, the 10-cm tube case produces unstable double bubbles in perfect qualitative agreement with experiments [34]. It is the multiple bubble formation that led to the selection of the 5-cm tube for flammability limit studies [15,16]. We may conclude that the simulations reproduce the physical behavior of flames with amazing correctness considering the crude nature of our assumptions, particularly concerning heat loss, flame thickness, and flame stretch effects. An important problem resolved by the calculation is the problem of stagnation zones in burned gases. The classical results [1,5,6] assumed that flow separation zones will occur in burned gases at the place where the flames meet the walls or behind flame segments concave to the mixture. These zones will substantially
888
LAMINAR PREMIXED FLAMES
reduce the global speed of the flame fronts. However, although stagnation zones were observed for 2-D flames stabilized in the tubes by external means [35], they were not observed in freely propagating flames [7,16,34]. The explanation is provided in Fig. 6. The formation of the separation zones behind the central concave fragment is suppressed by the strong three-dimensional side vortices. Thus, a freely moving flame creates a flow field that suppresses flow separation and, consequently, maximizes the flame propagation speed.
Conclusions
Fig. 7. Multiexposure photograph of the flame in a vertical tube with a time interval of 2 ms illustrating the tulip flame formation close to the open end of the tube. Conditions as in Fig. 6. In long vertical tubes, similar flame shape changes were observed at different locations, also very far from the open end.
1. A level set method based on the G equation was developed to simulate flame propagation in tubes for flames with heat release. 2. The flame geometry results from the nonlinear development of the Rayleigh-Taylor and LandauDarrieus instabilities of flames, both driven by the baroclinic effects. 3. The observed flame front speeds and shapes are close to experimental values within the range of the conditions tested. 4. 2-D and 3-D flame simulations lead to different speeds with substantially higher speeds for the 3D case. 5. Fast mushroom flames started at the closed end of vertical tube, and flames in general spontaneously formed “tulip” shapes. The effect is generated by the Landau-Darrieus instability and is very sensitive to flow asymmetry. A slight inclination of the tube (158 off vertical) prevents the conversion. 6. Whereas flames in small tubes (5 2 5 cm) form single bubbles even for multiwave initial perturbations, flames in larger tubes (10 2 10 cm) form multiple, nonstationary bubbles, in accordance with experiments. 7. Due to the self-adaptive nature of the flames reflected in the generation of complex 3-D vortices, the flow behind the flames does not contain separation regions that maximize the flame speed. Acknowledgments This work was supported by a grant from the Swedish Board for Technological Development (NUTEK).
REFERENCES
Fig. 8. Multiwave initial perturbation development in a 10 2 10 cm tube.
1. 2. 3. 4.
Maxworthy, T., Phys. Fluids 5:407–417 (1962). Landau, L. D., Acta Phys-Chem. URSS 19:77 (1944). Darrieus, G., Z. Phys. Chem. 88:641 (1944). Markstein, G., J. Aero. Sci. 18:199–209 (1951).
BAROCLINIC INSTABILITIES OF FLAMES 5. Zeldovich, Y. B. Istratov, A. G. Kidin, N. I., and Librovich, V. B., Comb. Sci. Technol. 24:1–13 (1980). 6. Pelce´, P., J. Phys. 46:503–510 (1985). 7. Pelce´-Savornin C., Quinard J., and Searby G., Combust. Sci. Technol. 58:337–346 (1988). 8. Awn, A. G., Ph.D. Thesis, London University, London 1979. 9. Awn, A. G., and Spalding, D. B., Second Conference of Mechanical Power Engineering, Air Shames University, Cairo, 1978. 10. Spalding, D. B., and Wu, J. Z. Y., PCH Phys. Chem. Hydrodynam. 7:353–384 (1986). 11. Dunn-Rankin, D., Barr, P. K., and Sawyer, R. F., Twenty First Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1986, p. 1291. 12. Rotman, D. A., and Oppenheim, A. K., Twenty First Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1986, pp. 1291–1303. 13. Gonzalez, M. Borghi, R., and Saouab, A. N., Combust. Flame 82:201–220 (1992). 14. Lee, S. T., and Tsai, Ch. H., Combust. Flame 99:484– 490 (1994). 15. Coward, H. F. and Jones, G. W., Bureau of Mines Bulletin 503, Washington, D.C., 1952. 16. Jarosinski, J., Prog. Energy Combust. Sci. 12:81–116 (1986). 17. Baker, G. R., Merion, D. I., and Orszag, S. A., J. Fluid Mech. 23:477–501 (1982). 18. Baker, G. R., and Moore, D. W., Phys. Fluids A 1:1451–1459 (1989). 19. Hirt, C. W., and Nichols, B. D., J. Comput. Phys. 39:201–225 (1981).
889
20. Fatemi, E., and Odeh, F., J. Comput. Phys. 108:209– 217 (1993). 21. Unverdi, S. O., and Tryggvason, G., J. Comput. Phys. 100:25–37 (1992). 22. Fines, A., Tsuruda, T., and Hirano, T., Combust. Flame 95:76–86 (1993). 23. Batley, G. A., McIntosh, A. C., Brindley, J., and Falle, S. A. E. G., J. Fluid Mech. 279:217–237 (1994). 24. Williams, F. A., in The Mathematics of Combustion, (J. D. Buckmaster, Ed.), SIAM, Philadelphia, 1985, p. 97. 25. Osher, S., and Sethian, J. A., J. Comput. Phys. 79:12– 49 (1988). 26. Sethian, J. A., J. Diff. Geometry 31:131–161 (1990). 27. Kerstein, A. R., Ashurst, W. T., and Williams, F. A., Phys. Rev. A 37:2728–2731 (1988). 28. Ashurst, W. T., and Sivashinsky, G. I., Combust. Sci. Technol. 80:159–164 (1991). 29. Aldridge, R. C.; Combust. Flame 90:121–133 (1992). 30. Zhu, J., and Sethian, J., J. Comput. Phys. 102:128–138 (1992). 31. van Leer, B., J. Comput. Phys. 14:361–370 (1974). 32. Lien, F. S. and Leschziner, M. A., Proceedings of the 5th International IAHR Symposium on Refind Flow Modeling and Turbulence Measurements, Paris, Sept. 1993. 33. Zhou, G., Davidson, L. and Olsson, E., Proceedings of the 14th ICNMFD, Lecture Notes in Physics, (S. M. Deshpande, S. S. Desai, and R. Narsimha, Eds.), Springer-Verlag, Berlin, 1995, No. 453, p. 372–377. 34. Von Lavante, E., and Strehlow, R. A., Combust. Flame 59:123–140 (1983). 35. Uberoi, M. S., Phys. Fluids 2:72–78 (1959).
COMMENTS Luc Bauwens, The University of Calgary, Canada. Did you have the opportunity to look at the effect of the Markstein length, especially when it becomes small? Can you say a few words on how the source of volume and vorticity at the front was implemented numerically in your code? (This strikes me as a very difficult problem.) Author’s Reply. We have not studied the effects of the Markstein length in the small value range. Our method does not allow us to resolve the flame structure and thus is unsuitable for that purpose. Furthermore, to avoid strong grid sensitivity, we had to smooth locally the flame interface after each time step using a sinusoidal smoothing function. To be able to evaluate the effects of volumetric expansion, we have to know the area element of flame front in each control volume. The computation of the area element is done using the following sequence. First we redistribute the G from the cell center to the cell vertexes. Then, we obtain G 4 0 by interpolation between two adjacent values of G which have opposite signs. After getting the data for G 4 0, we can reconstruct the local flame geometry. Then
the expansion source term is defined using Eq. (4) from the paper. The addition of the expansion source closes the hydrodynamic problem and there is no need to introduce vorticity at the interface since it is generated automatically by the pressure-density interaction terms in the N.S. equations. The baroclinic torque is insensitive to the actual density distribution in the flame and is determined by the pressure gradient and the total density difference across the flame. Practically, to avoid unphysical oscillations at the interface, we used the following equation for for q(w) in the flame zone:
5
qu qb pw q¯ ` Dqsin 2a
1 2
if w . a if w , a otherwise
(1)
where q¯ 4 (qu ` qb)/2 Dq 4 (qu 1 qb)/2 and a is a finite thickness of an order of 0(h), (where h is the mesh size).