Energy & Fuels 2009, 23, 2903–2907
2903
Study of Methane Propagating Flame Characteristics Using PDF-Monte Carlo Model and Reduced Chemical Kinetic Scheme Ridha Ennetta,* Marzouk Lâjili, and Rachid Said Ionized and ReactiVe Media Studies (EMIR), Preparatory Institute of Engineers Studies of Monastir (IPEIM), AVenue of Ibn El Jazzar, 5019 Monastir, Tunisia ReceiVed January 17, 2009. ReVised Manuscript ReceiVed April 25, 2009
The main purpose of this work is to simulate turbulent premixed and expanding flame in an adiabatic combustion chamber without taking into account flame-walls interactions. The turbulence is supposed to be isotropic and homogeneous with no decay or decreasing, but spatially correlated. The chemistry is represented by a four-step scheme of methane-air combustion. A one-dimensional simulation is considered because of the spherically symmetric of the problem. Ignition occurs in the middle of the domain thanks to an energy source such a spark. The simulation, based on Monte Carlo scalar probability density functions (PDF) transport method, is used under different equivalence ratios (ER) and different turbulence intensities (u′). We have placed our emphasis on some flame characteristics such as the flame mean radius, the turbulent flame radius, the flame propagation velocity, and the flame brush thickness. The results of our simulations, carried out in similar conditions to some available experiments, are in good agreements. Indeed, we notice that the flame radius is enhanced by ER at constant u′, and by u′ at constant ER. The flame brush thickness shows a quick growth at the first stage and a moderate growth in the intermediate regime, but a bending effect is observed in the final stage of the so-called “fully developed flame”.
Due to their fundamental importance for premixed combustion theory, turbulent flame characteristics were a subject of a large number of numerical investigations for many decades.1 Reliable and flexible computational modeling is the key in achieving the objectives of such investigations. Direct numerical simulation (DNS), Large eddy simulation (LES), and Reynoldsaveraged Navier-Stokes (RANS) are widely considered as three principal approaches in computational turbulent combustion.2 DNS provides a very accurate, model-free representation of the unsteady evolution of turbulent flows. However, applications are limited by the computational power.3-5 LES has been the subject of much modern research and is increasingly becoming more popular.6 RANS is the most popular approach for engineering applications.7 Various models have been proposed to describe turbulent premixed flames.8 Probabilistic approaches have proven to be effective in this regard.3,6 This approach follows from the
definition of the fine-grained density function.9 The ensemble average of the fine-grained density function is considered and is termed the probability density function (PDF).10 The primary advantages of probability methods are: (i) they provide closed-form representation of chemical source terms, and (ii) they are applicable to both premixed and non-premixed flames. The probability density functions (PDF) methods have been popular since 1970s.11,12 Development of the Lagrangian Monte Carlo particle methods13 has enabled PDF calculations to be conducted of a variety of flame configurations.14 The main purpose of this work is to simulate a premixed turbulent flame in a constant volume vessel. The adopted conditions are similar to those realized in many experiments. The Monte Carlo scalar PDF transport method that is the basis of our numerical simulation has been used for few decades by many authors to simulate ignition time and the flame growth in turbulent combustion case.15-18 Nowadays, this method continues to be efficient for calculating flame characteristics taking
* To whom correspondence should be addressed. Phone: (+216)98487430; fax: (0216)73500277; E-mail:
[email protected], Ridha.Ennetta@ isetjb.rnu.tn. (1) Lipatnikov, N.; Chomiak, J. Progr. Energ. Combust. Sci. 2002, 28, 1–74. (2) Givi, P. Spectral and Random Vortex Methods in Turbulent Reacting Flows. In Turbulent Reacting Flows; Libby, P. A., Williams, F. A. Eds.; Academic Press: London, UK, 1994; Ch 8, 475-572. (3) Pope, S. B. Proc. Combust. Inst. 1990, 23, 591–612. (4) Poinsot, T.; et al. Prog. Energy Combust. Sci. 1996, 21, 531–576. (5) Vervisch, L.; Poinsot, T. Annu. ReV. Fluid Mech. 1998, 30, 655– 691. (6) Givi, P. AIAA J. 2006, 44 (1), 16–23. (7) Pope, S. B. Advances in PDF Methods for Turbulent Reactive Flows. In AdVances in Turbulence X; Andersson, H. I., Krogstad, P. A. Eds.; CIMNE: 2004; pp 529-536. (8) Veynante, D.; Vervisch, L. Prog. Energy Combust. Sci. 2002, 28 (3), 193–266.
(9) Lundgren, T. S. Phys. Fluids 1967, 10 (5), 969–975. (10) O’Brien, E. E. The Probability Density Function (PDF) Approach to Reacting Turbulent Flows. In Turbulent Reacting Flows, Topics in Applied Physics; Libby, P. A., Williams, F. A., Eds.; Springer-Verlag: Heidelberg, 1980; Ch 5, 44, 185-218. (11) Dopazo, C.; O’Brien, E. E. Combust. Sci. Technol. 1976, 13, 99– 112. (12) Pope, S. B. Combust. Flame 1976, 27, 299–312. (13) Pope, S. B. Annu. ReV. Fluid Mech. 1994, 26, 23–63. (14) Bilger, R. W.; et al. Proc. Combust. Inst. 2005, 30, 21–42. (15) Haworth, D. C.; Pope, S. B. J. Comput. Phys. 1987, 72, 311–346. (16) Galzin, F. Contribution à la modélisation de la combustion dans les moteurs à allumage commandé. Ph.D. Thesis, University of Rouen: France, 1996. (17) Fruchard, N. L’allumage dans les moteurs à essence: une modélisation et des applications. Ph.D. Thesis, University of Rouen: France, 1995. (18) Borghi, R.; Champion, M. Mode´lisation et the´orie des flammes; TECHNIP: Paris, France, 2000.
1. Introduction
10.1021/ef9000456 CCC: $40.75 2009 American Chemical Society Published on Web 05/13/2009
2904
Energy & Fuels, Vol. 23, 2009
Ennetta et al.
into account turbulence and chemistry interaction using reduced chemical kinetic mechanisms.19,20 For our case, turbulence is supposed to be homogeneous and isotropic, and its characteristics17 are those of the frozen k-ε model. Moreover, the turbulence does not present any decay or dissipation rate decrease because we suppose that the calculation domain is an open space without flame-walls interactions, nut turbulent velocities are spatially correlated, and this fact makes our study different from classical PDF methods. In the present study, chemistry is described by the four-step chemical kinetic mechanism of Jones and Lindstedt.21 The mixture is composed of methane as fuel and air as oxidizer. The results of our simulations are in good agreements with some available experimental data obtained in fan-stirred spherical bombs (experimental apparatus and measurement methods are well described in literature30-32,43,44). We have notice an enhancement, by the equivalence ration (ER) and the turbulent intensity (u′), of the flame radius and the turbulent velocity. The flame brush thickness shows a quick growth at the first stage, when the turbulent velocity grows too, and a moderate growth in the intermediate regime after the turbulent speed reaches its asymptote, but a bending effect is observed in the final stage of the so-called “fully developed flame”. 2. Theoretical Basis 2.1. Formulation. The equations that determine our problem are, respectively: • the mass conservation equation written in spherical coordinates, whose solution gives us the expansion velocities due to temperature gradient between hot products and fresh gas: 2 1 ∂(r Fu˜r) ∂F + 2 )0 ∂t ∂r r
(1)
• the Lagrangian joint PDF transport equation: ∂ (F(ψ)f (V_,ψ)) + ∂t _ u_,Φ_ _
3
∑ ∂∂ [F(ψ_ )Vif
_,ψ _ )] u_,Φ _ (v
)
xi
i)1
3
-
∑ ∂∂ [F(ψ_ )〈A |V_,ψ_ 〉f i)1 N
Vi
R)1
ψR
i
_,ψ _ )] u_,Φ _ (V
(2)
∑ ∂∂ [F(ψ_ )〈θ |V_,ψ_ 〉f
_,ψ _ )] u_,Φ _ (V
R
where
(
)
(
)
Ai )
∂P 1 ∂τij + FFi ; F ∂xj ∂xj
θR )
R 1 ∂Jk + Fω ˙ R + SR ; R ) 1, N F xk
i, j ) 1, 3
and
• the perfect gas equation: ˜ j ) F RT P M
(3)
where F is the density, u˜r is the Favre averaged radial velocity corresponding to the expansion velocity due to the temperature (19) Ren, Z.; Pope, S. B. Combust. Flame 2004, 136, 208–216. (20) Garrick, S. C.; and Interante, V. 9th Int. Symp. Flow Visualisat., Heriot-Watt University: Edinburgh, 2000; 1-10. (21) Jones, W. P.; Lindstedt, R. P Combust. Flame 1988, 73, 233–249.
gradient between hot products and fresh gas. In eq 2, Φ and u are, respectively, a scalar vector and a velocity vector in the physical space and to which correspond respectively the vectors of random values, Ψ and V in the conditional space. Fu,φ is the PDF of velocities and scalars. The change terms that characterized the stochastic process Ai and θR present the following terms: τij that is the strain tensor, Fi is the stirred force per volume ˙ R is unity, P is the pressure, JRk represents the diffusive fluxes, ω the reaction rate, and, finally, SR is the source term. In eq 3, T˜ is Favre averaged temperature, and M is the molar mass of the mixture. 2.2. PDF-Monte Carlo Method. The use of PDF constitutes a potential solution to describe the evolution of turbulent reactive flows in which fluctuation terms need statistical treatment. The type of PDF that we use in this work is the evolution PDF (transported PDF) called Pope’s method.15 This method uses a Monte Carlo particle solver, and the form of the PDF may freely evolve. It can be highlighted that the high dimensionality of underlying PDF scalar transport equation requires Monte Carlo stochastic solution methods. Monte Carlo method evokes the representation of the PDF with a whole of elements distributed throughout the flow field, and from which the moments of interest may be calculated. In the Lagrangian case, the elements (particles) are free to roam the physical domain as dictated by the hydrodynamic field, and the composition of the elements changes only due to mixing and reaction.22 The calculation domain is divided into a given number of cells Nc. Initially, each one contains Ni particles. These particles move in the domain thanks to the following velocities: • gas expansion velocity due to temperature gradient between burned gases (hot products) and “fresh gases”. • turbulent diffusion where a correlation velocity deduced from turbulence spectrum is respected. 2.3. Simulations Parameters. The physical environment of the current computation is a one-dimensional (1D) spherically symmetric domain with changes only in the radial direction being considered. The premixed methane-air mixture is initially at rest, and flame initiation is triggered by a localized deposit of thermal energy generated by a hot source (spark, laser ignition, pyrotechnic devices, etc.), making the temperature grow highly in this region,17 as is the case in internal combustion engines or other such devices. The computational environment in the radial direction is composed of Nc ) 500 cells each with a specified size ∆r. Initially, each cell contains a certain number of particles (Ni ) 1000 for example). Each particle has a velocity and some themochemical variables, such as concentration and temperature. The pressure is held constant throughout the calculation at 1 atm. The numerical procedure tracks each fluid particle during each time step to account for convective transport by the mean flow and by turbulence in space and in time. Chemical and thermal composition of particles changes due to chemical reaction and molecular mixing. To detect more precisely the flame front evolution, we have chosen a cell size (∆r ) 0.25 mm) and a calculation time step (∆t ) 0.1 ms) which remain constant along the time and overall the domain. These lasts parameters are chosen from literature. Indeed, in a recent work, Raman et al.23 have performed some numerical tests to determine the optimal global time step, the optimal space step, and to estimate the increase in computation (22) Lajili, M.; et al. Prog. Comput. Fluid Dynamics 2008, 8 (6), 331– 341. (23) Raman, V.; et al. Combust. Flame 2004, 136, 327–350.
Methane Flame Characteristics with PDF-MC Method
Energy & Fuels, Vol. 23, 2009 2905
load with the number of particles per cell and consequently the number of particles in the whole domain. It can be concluded that the tracking time scales linearly with particles number, size grid, and global time step. Moreover, it is found for the methane-air studied that a time step of 1.0 × 10-4 s yields consistent results for both reacting and nonreacting cases.23 Turbulent length scale and turbulent time scale values are equal to those used in experimental cases and given by the following formulas: 3
lt ) 0.164
k /2 ε
k τt ) 0.3 ε
flame turbulent velocity, and the flame-brush thickness, which represents a great challenge for the turbulent combustion community. The predicted values are compared to the experimental results realized by many authors. The geometry studied was a parallelepiped vessel with a constant volume. Ignition occurs in the middle of the domain then a flame kernel grows spherically. 3.1. Flame Mean Radius. The flame mean radius (Rf) is defined by Lecordier,30 in 2D configuration, as the radius of a circle that contains the same surface of burned gas (Sb).
(4) Rf ) (5)
The following input values are required for the simulation: time step ∆t, cell size ∆r, turbulent kinetic energy k, dissipation of turbulent kinetic energy ε, normalized mean fuel concentration, Cfini, ignition energy x, pressure, ignition time tignit, ignition gap rignit, initial number of particles in each cell Ni, adiabatic flame temperature Tadiab, and initial temperature Tini. 2.4. Chemical Kinetics. Detailed chemical kinetic descriptions of hydrocarbon combustion may require the tracking of hundreds of chemical species and thousands of reaction steps. For the foreseeable future, CPU time and computer memory limitations will prohibit implementation of fully detailed descriptions of combustion chemistry into CFD simulations of combustion hardware.24-26 As a result, it is important to minimize this number while retaining essential features of the detailed chemistry. Reduced chemical kinetic mechanisms, which can represent important aspects of the behavior of these detailed mechanisms using few enough scalars that they can be implemented into CFD simulations, offer large potential improvement in the modeling of practical combustion devices.27 The feasibility of these mechanisms in the simulation of internal combustion engines (ICE) was demonstrated in a previous study.28 In the present study, the four-step reaction mechanism of Jones and Lindstedt21 was chosen thanks to its good results in 1D propagating flame21,29 and in ICE simulations:28
π1 S
b
(6)
However, in our case (1D configuration), the flame mean radius is defined as the distance between the ignition center and the flame front position. The latter is determined when finding a “pitchfork” of temperature in the vicinity of 600 K, which corresponds to silicon oil vaporisation temperature in the experimental case.
Figure 1. Flame mean radius evolution vs time.
CH4 + 1/2O2 f CO + 2H2 CH4 + H2O f CO + 3H2 H2 + 1/2O2 T H2O CO + H2O T CO2 + H2 3. Results and Discussion In this section we try to study more precisely the turbulent flame front characteristics such as the flame mean radius, the (24) Montgomery, C. J.; et al. 46th AIAA Aerospace Sciences Meeting and Exhibit; AIAA: Reno, Nevada, 2008; p 1014. (25) Glassman, I.; Yetter, R. A. Combustion, 4th ed.; Elsevier: 2007. (26) Skevis, G.; Goussis, D. A.; Mastorakos, E. Int. J. Altern. Propul. 2007, 1(2/3), 216-227. (27) Montgomery, C. J.; et al. J. Propul. Power 2002, 18 (1), 192–198. (28) Ennetta, R.; Hamdi, M.; Said, R. Therm. Sci. 2008, 12 (1), 43–51. (29) Truffin, K. Modelling of Methane/Air Flames Using Detailed and Reduced Kinetic Schemes; CERFACS, INP Toulouse: France, 2001. (30) Lecordier, B. Etude de l’interaction d’une flamme prémélangée avec le champ aérodynamique, par association de la tomographie Laser et de la P.I.V. Ph.D. Thesis, University of Rouen, France, 1997. (31) Bradley, D.; et al. Combust. Flame 1994, 99, 562–572. (32) Hainsworth, D. Study of Free Turbulent Premixed Flames. M.Sc. Thesis, Massachusetts Institute of Technology, USA, 1985.
Figure 2. Flame mean radius evolution vs time.
Figures 1 and 2 show the flame mean radius evolution function of time for different turbulence intensities and different equivalence ratios. We can see that the predicted values found by our simulations are in good agreements with the experimental ones. In order to study the equivalence ratio (ER) effect on flame mean radius, one maintains the turbulence intensity u′ constant and varies the ER. The effect of this last parameter on flame
2906
Energy & Fuels, Vol. 23, 2009
Ennetta et al.
mean radius is shown in Figures 1 and 2. It is worth noting that for lean flames the flame mean radius is quite parabolic. However, as we increase ER to stoichiometry; the flame mean radius becomes strongly parabolic. Also, Figure 1 shows that the effect of turbulence intensity is more significant than the ER one, on the flame mean radius. These results are in perfect agreement with those announced by Bradley31 and Hainsworth.32 3.2. Flame Propagation Velocity. The flame propagation velocity is defined as the differential of the flame mean radius versus time. Figures 3-5 compare the flame propagation velocities found by our simulations and experimental results presented by Abdel-Gayed et al.,33 Kobayashi et al.,34 and Aldredge et al.35 There is a very good agreement between experimental and simulation results. Figure 5. Normalized turbulent velocity vs normalized rms turbulent velocity.
Figure 3. Normalized turbulent velocity vs normalized rms turbulent velocity in spherical flame.
Figure 4. Dependences of turbulent flame speed on rms turbulent velocity.
Our comparisons were made with different experiments. Each one describes a certain stage of flame propagation. Description of the different stages of the turbulent flame will be discussed in details in Section 3.3. It is obviously clear that under Kobayashi experimental34 conditions, the flame propagation is in a transient regime (33) Abdel-Gayed, R. G.; Bradley, D.; Lawes, M. Proc. R. Soc. Lond., Ser. A 1987, 414, 389–413. (34) Kobayashi, H.; et al. 26th Symp. (Int.) Combust.; The Combustion Institute: Pittsburgh; PA, 1996; pp 389-396. (35) Aldredge, R. C.; et al. Combust. Flame 1998, 115, 395–405.
characterized by a constant propagation velocity (but at the same time a growing flame width). This phenomenon seems to be in contradiction with Abdel-Gayed et al. experiments33 in the spherical bomb, which only describe the initial stage of the turbulent flame characterized by a growth in both the flame speed and width. However, in our case, we’ve prolonged simulations to enough time until obtaining a constant turbulent flame speed. So, St values are taken when the propagation velocity reaches its asymptote. Thus, our comparison with the experimental results of Kobayashi et al.34 is legitimate. In this following section, we place our emphasis on equivalence ratio and turbulence intensity effects on flame propagation velocity. Concerning equivalence ratio effect on flame propagation velocity, it has been known for many decades, in laminar combustion regimes, that the flame burning velocity is maximal when the equivalence ratio is around stoichiometry value.36,37 In recent work, Bradley38 affirmed that for turbulent case the flame propagation velocity, which strongly depends on burning velocity, increases with ER. These declarations are justified in Figure 4. Also, we can remark that in the case of lean flames (ER ) 0.9), the flame propagation velocity increases slowly. To discover the effect of the turbulence on turbulent premixed flame characteristics, we have represented the flame velocity ratio St/SL versus u′/SL (see Figures 3 and 5). It is obviously clear that the flame propagation velocity increases with turbulence intensity. These results are confirmed by theory. In fact, a high level of turbulence intensity makes microscale mixing more efficient and turbulence diffusivity stronger, so the flame will be much faster (see Figures 3 and 5.) It is worth noting that this big scatter of experimental data (see Figure 3) is a widely recognized fact.1 It is attributed to both the experimental error and inadequate or incomplete theoretical assumptions.32,39 We can affirm that, at weak and moderate turbulence, the flame-burning velocity ratio varies function of the dimensionless turbulence intensity, responding to Damko¨hlers’ model in the case of n ) 1. Thus, the expression takes the following form: (36) Metgalchi, M.; Keck, J. C. Combust. Flame 1980, 38, 111–129. (37) Konnov, A. A.; Dyakov, I. V. Combust. Flame 2004, 136, 371– 376. (38) Bradley, D.; et al. Combust. Flame 2003, 133, 415–430. (39) Lapatnikov, A. N.; Chomiak, J. Prog. Energy Combust. Sci. 2005, 31, 1–73.
Methane Flame Characteristics with PDF-MC Method St
/SL ) 1 + C1u′/SL
Energy & Fuels, Vol. 23, 2009 2907
(7)
This formula is similar to that given by Gu¨lder40 for weak Reynolds number. 3.3. Flame Brush Thickness. Figure 6 shows the flamebrush thickness δt (called also flame width) evolution vs time. This flame width is defined by: δt )
1 max |∇c|
(8)
where c ) (T - Tu)/(Tb - Tu) is the progress variable defined as a normalized temperature (u for unburned and b for burned gases).
equilibrium of both small-scale and large-scale wrinkles, in accordance with the Taylor theory.47 Thus, we notice an increasing of the brush width δt like u′t and relatively fast growth of St. Concerning the intermediate stage, it is marked by the quasiequilibrium of the small-scale wrinkles while the large scales still practically nonequilibrium. In fact, we have in this case an intermediate asymptotic state. So, these transient flames are characterized by nearly constant flame speed and at the same time by increasing brush width, and that is why they are called intermediate steady propagation (ISP) flames. As for the final steady stage, both large- and small-scale wrinkles are in equilibrium; it is the final asymptotic state. It is steady state (SS) flames analyzed in many publications (usually in the form of the classical 1D stationary turbulent flame). Indeed, this stage cannot be reached by experiments because it needs a long period time of flame study, a thing that is not possible experimentally. However, by simulation and especially by Monte Carlo simulations, this study is possible with current CPU capacities. 4. Conclusion
Figure 6. Dimensionless mean flame brush thickness (δt/Lt) vs dimensionless time (t/τt).
The agreement between simulation and experimental results is satisfactory. The asymptotic tendency of δt, affirmed by Galzin16 and Liptanikov and Chomiak,1 and observed by Moreau,41 is predicted by our simulation. Although, Karlovitz42 and Scurlock43 have pointed out that the increase of δt is mainly controlled by the turbulent diffusion law, whereas flame propagation reduces δt and it can reach approximately constant values after the development phase. However, in some experiments, such as engines, for example, this constancy of δt can be caused by other effects (e.g., turbulence decay or wall influence) rather than by reaching the regime of turbulent flame propagation characterized by a fully developed δt. Indeed, Zimont, in his last book,46 affirms that the first initial stage of flame kernel formation is characterized by the nonˆ . L. Combust. Flame 2000, 120, 407–416. (40) Gu¨lder, O (41) Moreau, P. Aerospace Science Meeting, AIAA Paper 1977, 77, 49. (42) Karlovitz, B.; et al. J. Chem. Phys. 1951, 19, 541–547. (43) Scurlock, A. L.; Grover, J. H. Proc. Combust. Inst. 1953, 4, 645– 658. (44) Renou, B.; et al. 18th Int. Colloq. on the Dynamic of Explosions and ReactiVe Systems: Seattle, Washington, USA, 2001. (45) Atashkari, K. Towards a general correlation of turbulent premixed flame wrinkling. In Engineering Turbulence Modelling and Measurements; Rodi, W., Laurence, D. Eds.; Elsevier: Amsterdam, 1999; Vol. 4, pp 805814. (46) Zimont, V. Kolmogorov’s legacy and turbulent premixed combustion modelling. In New DeVelopments in Combustion Research; William J. Carey Inc, Ed.; Nova Science Publishers: New York, 2006; Ch 1, pp 1-93.
Turbulent premixed combustion was numerically simulated in an adiabatic constant volume vessel. Simulation was based on the PDF-Monte Carlo method. Turbulence was supposed to be homogeneous and isotropic. The chemistry was described by a four-step reaction mechanism of methane combustion. Confrontation with many experimental results gave satisfactory agreements, especially on mean flame radii, turbulent propagation velocity, and flame-brush thickness. The current model correctly predicts trends such as the effect of increasing turbulence intensity and equivalence ratio on the premixed methane-air flame propagation characteristics, such as the mean radius and the turbulent speed. Indeed, we have notice an increase of the mean flame radius and the turbulent speed when ER or u′ increases, but the effect of turbulence intensity was more significant. Also, simulations show that turbulence intensity is the dominant factor in the determination of the turbulent flame speed. After ignition and during the first stage of the flame propagation, both the turbulent velocity and the brush width are growing. After that, a bending effect was observed in the turbulent flame speed evolution (characterizing the second stage of the flame propagation called also the ISP regime). Then, PDFMonte Carlo simulations demonstrate that the flame brush thickness has a tendency to be constant at increasing time (this stage is called the SS regime). Experiments were unable to reach this last regime. Thus, the PDF-Monte Carlo method could be a very useful tool to explain theoretical prediction of this regime. Finally, it will be very interesting, in a next work, to study the propagation flame characteristics by more considering turbulence-chemistry interactions. Also, it would be useful to take into account more detailed chemistry allowing us to evaluate pollutant emissions. EF9000456 (47) Taylor, G. I. Proc. R. Soc. Lond., Ser. A 1935, 151, 465–478.
