- Email: [email protected]
Premixed turbulent combustion modeling using tabulated detailed chemistry and PDF
Proceedings of the
Combustion Institute
Proceedings of the Combustion Institute 30 (2005) 867–874
www.elsevier.com/locate/proci
Premixed turbulent combustion modeling using tabulated detailed chemistry and PDF B. Fiorinaa,b,*, O. Gicquela, L. Vervischc, S. Carpentierb, N. Darabihaa a b
Laboratoire EM2C, UPR 288 CNRS, Ecole Centrale Paris, F92295 Chatenay Malabry, France Poˆle CCMF, Direction de la Recherche, Gaz de France, F93211 La Plaine Saint Denis, France c CORIA, UMR 6614 CNRS, INSA de Rouen, F76800, Saint Etienne du Rouvray, France
Abstract Tabulated chemistry and presumed probability density function (PDF) approaches are combined to perform RANS modeling of premixed turbulent combustion. The chemistry is tabulated from premixed flamelets with three independent parameters: the equivalence ratio of the mixture, the progress of reaction, and the specific enthalpy, to account for heat losses at walls. Mean quantities are estimated from presumed PDFs. This approach is used to numerically predict a turbulent premixed flame diluted by hot burnt products at an equivalence ratio that differs from the main stream of reactants. The investigated flame, subjected to high velocity fluctuations, has a thickened-wrinkled structure. A recently proposed closure for scalar dissipation rate that includes an estimation of the coupling between flame wrinkling and micromixing is retained. Comparisons of simulations with experimental measurements of mean velocity, temperature, and reactants are performed. 2004 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Turbulent combustion; Premixed flame; Modeling; PDF; Tabulated chemistry
1. Introduction The increase of computational resources allows the intensive use of numerical tools to predict chemical structure of flames, understand flame stabilization processes, and help in burner design. Several strategies have been proposed to model and to compute turbulent flames [1–5]. In the context of premixed turbulent combustion, Bradley et al. [6,7] have tabulated detailed chemistry with a collection of one-dimensional premixed laminar flames, then estimating the mean flame structure
*
Corresponding author. Fax: +33 0 1 49 22 53 73. E-mail address: benoit.fi[email protected] (B. Fiorina).
after averaging with presumed probability density function (PDF). Bray et al. [8] have formulated the hypothesis of a fully bimodal PDF of progress variable to compute flames. A combined strategy is adopted in this paper, to perform RANS modeling of a premixed turbulent flame piloted by hot burnt products [9]. Detailed chemistry effects are included in the turbulent combustion model by the means of premixed laminar flamelets, following a recently developed procedure for chemistry tabulation [10–12], called FPI. The chemical flame structure is stored inside a look-up table in terms of three coordinates, which represent the progress of reaction, the local dilution of reactants with products at a different equivalence ratio, and the enthalpy level. PDFs are used to presume
1540-7489/$ - see front matter 2004 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2004.08.062
868
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
the mean flame structure. Two closures are tested for the scalar dissipation rate of the progress variable, which is one of the fundamental ingredients of the presumed PDFs. A classical linear relaxation formulation and a novel expression are introduced, the latter derives from the Bray–Moss–Libby (BML) analysis to combine information on the flame topology with turbulent micromixing [13]. An ONERA combustion rig that consists of a premixed turbulent flame piloted by hot gases [9] is considered as a test case. This flame is subjected to high velocity fluctuations, leading to strong interactions between chemistry and turbulence.
2. Modeling premixed turbulent combustion 2.1. Chemistry tabulation Predicting the flame shape, its stabilization process, and pollutant emissions in practical burners requires the use of detailed chemistry. Various techniques have been proposed to account for the complexity of chemistry in turbulent flames. ISAT (In situ adaptive tabulation) [14] is based on in situ generation of look-up tables, which are constructed from the direct solving of the time evolution of species concentrations. This method allows for highly accurate simulations of jet [2] and lifted flames [3]. The ILDM methodology [15] is another alternative to chemistry tabulation, which is grounded on a direct mathematical analysis of the dynamic behavior of the non-linear response of the chemical system. Relevant attractive subspaces are constructed from the direct solving of the time evolution of the species concentration. This chemistry tabulation method has been applied to a large variety of problems, including nonpremixed turbulent jet flames [4]. ILDM is an exact description of the chemical behavior of the system. Nevertheless, to provide accurate results with basic hydrocarbon fuels (methane or propane), the dimension of the manifold should at least be around four or five. To combine the advantages of mathematically reduced databases with a good estimation of the whole chemical flame structure, a tabulation strategy called FPI has recently been proposed [11]. FPI is based on physical considerations and is similar to FGM [16]. It determines the reduced chemical subspace effectively accessed in a complex geometry configuration using a set of reference laminar flame computations. In the case of premixed laminar flames, premixed flamelets were used to build the database [10], while auto-ignition computation was retained when the estimation of the ignition delay is crucial [17,18]. The FPI method was also successfully applied to non-adiabatic partially premixed laminar flames [12].
The objective of this paper was to model a premixed turbulent flame diluted by a pilot flame at a different equivalence ratio, this pilot stream containing only hot burnt products. Therefore, the equivalence ratio only weakly varies over the reaction zone surface. There are no reactants in the pilot stream, and the premixed character of the reactive layer is thus mainly preserved. It is then assumed that the chemical flame structure can be reproduced by a collection of one-dimensional laminar premixed flames. A set of laminar premixed flames is computed at various equivalence ratios and for different enthalpy levels, using detailed chemistry. A FPI look-up table is constructed to store all chemical information, with three coordinates: Yc, a progress of reaction; Yz, a mixture fraction characterizing the equivalence ratio; and h, the enthalpy. Yc and Yz are both suitable combinations of species mass fraction, which, respectively, reproduce the progress of reaction and the local mixture of fuel and oxidizer. The third coordinate h, the mixture specific enthalpy, enables to take account for heat losses. Details about the database generation are provided in [12]. All chemical species involved in the detailed kinetics scheme, the temperature, and the thermochemical data are then tabulated in terms of the three coordinates (Yc, Yz, h). If u denotes a species mass fraction or any thermochemical data, then ð1Þ
u ¼ uðY c ; Y z ; hÞ:
This tabulation appears as an interesting tool to perform very fast calculations including complex chemistry effects, specifically when very complex geometries are involved. 2.2. Presumed PDF modeling The statistical properties of the reactive and diffusive layers are easily described from the joint PDF, Pe ðY c ; Y z ; h Þ of the variables chosen to describe the internal flame structure: Z had Z Y c Z Y z f b ~¼ uðY c ; Y z ; h Þ u hmin
Y cf
Y zo
Pe ðY c ; Y z ; h Þ dY c dY z dh ;
ð2Þ
where hmin and had are, respectively, the minimum and the adiabatic values of enthalpy, Y cf and Y cb are, respectively, the values of Yc in fresh and burned gases and Y zo and Y zf correspond to values of Yz in leanest and richest inlets of the physical domain. In fact, the lowest enthalpy limit is the one obtained by computing the chemical equilibrium at T = Tmin, where Tmin is the minimum value of temperature encountered in the turbulent flow. The three coordinates of the FPI database are normalized according to the considered problem, leading to z, hn, and c:
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
z¼
Y z Y zo ; Y zf Y zo
hn ¼
c¼
h hmin ðY z Þ ; had ðY z Þ hmin ðY z Þ
Y c Y cf ðY z ; hÞ : Y cb ðY z ; hÞ Y cf ðY z ; hÞ
ð3Þ
one can readily 002 ; and h : c002 ; ez ; zf n
ð4Þ
c002
ð5Þ
Due to its careful choice and normalization, the progress of reaction c and the mixture fraction z can be considered as two statistically independent variables as mentioned in [20] and further tested from measurements in [13]. It is then possible to introduce specific turbulent model based on this assumption [13]. In the present paper, this model is extended using the normalized enthalpy hn that is also assumed to be statistically independent from other variables. The introduction of the enthalpy as an additional coordinate is motivated by the need to account for flame wall interaction, driving heat transfer to the wall and consequently, the mean turbulent flame shape. Notice that the assumed statistical independence would not be true for the primitive variables Yc, Yz, and h. One may then write for any quantity u (mass fractions or source terms) Z 1Z 1Z 1 ~¼ uðc ; z ; hn Þ u 0
0
0
P ðc Þ Pe ðz ÞP ðhn Þ dc dz dhn ;
ð6Þ
where P ðc Þ and Pe ðz Þ take the form of b-functions that are defined by their first and sec002 Þ, respectively. ond moments ðc; c002 Þ and ðez ; zf Because enthalpy fluctuations are not that large in the considered test case, P ðhn Þ dðhn hn Þ is assumed, and the PDF of the enthalpy becomes a delta function centered on hn . ~ using Estimation of any averaged quantity u relation (6) requires the evaluation of first and second moments of c and z, and of the averaged normalized enthalpy. Starting from relation (5), one can write the conditional average of c at z = z* and hn ¼ hn
ðc j z ¼ z ;hn ¼ hn Þ ¼
ðY c j z ¼ z ;hn ¼ hn Þ Y cf ðz ;hn Þ : Y cb ðz ;hn Þ Y cf ðz ;hn Þ ð7Þ
derive
869
relations
g e 2 f2 ð Ye c Ye cf Þ2 Y 002 c þ Y c Y cf ¼ g ð Ye cb Ye cf Þ2 ðY cb Y cf Þ2 ð Ye c Ye cf Þ qðY cb Y cf ÞY cf 2 ; g ðY cb Y cf Þ2 qY cb qY cf
ez ¼
Ye z Ye zo ; e Y zf Ye zo
002 ¼ zf
hn ¼
g Y 002 z ; ð Ye zf Ye zo Þ2
e he h min : e h ad e h min
for
ð10Þ
ð11Þ
ð12Þ
ð13Þ
2.3. Balance equations In addition to the classical mass, momentum, mixture fraction (mean, Ye z , and fluctuations, g Y 002 z , energy, and RANS k balance equations, Y 002 transport equations for Ye c and g c need to be solved. They may be written as: l o q Ye c _ Y ; þ r ð qe u Ye c Þ ¼ r t r Ye c þ x c ot r l o qg Y 002 t c þ r ð qe ug Y 002 rg Y 002 c Þ ¼ r c ot r l þ 2 t j r Ye c j2 2qD j rY 00c j2 þ 2Y 00c x_ Y c ; r
ð14Þ
ð15Þ
where usual gradient transport is assumed for modeling convection by velocity fluctuations, lt is the eddy viscosity, and r = 0.6 is the turbulent Schmidt _ Y and Y 00 x_ Y number. The chemical source terms x c c c are obtained after averaging the FPI table with the PDFs via Eq. (6). The scalar dissipation rate of the progress of reaction may be write ten q v Y c ¼ qD j rY c j2 ¼ qD j r Ye c j2 þ qDj rY 00c j2 , whose fluctuating part appears in Eq. (15) and requires additional modeling.
As c, z and hn are supposed to be independent ðc j z ¼ z ; hn ¼ hn Þ ¼ c:
ð8Þ
Therefore, integrating relation (7) with Pe ðz Þ leads to the following approximation: c’
Ye c Ye cf : e Y cb Ye cf
ð9Þ
Averaged quantities Ye cf and Ye cb are estimated with relation (6). Following a similar analysis,
2.3.1. Linear relaxation for e vY c Assuming equilibrium between production and dissipation in the balance equation for the scalar dissipation rate, the widely used linear relaxation type closure appears [19] 002 e q g ð16Þ Y ; q v Y c ¼ C k c where C is a modeling constant. Because the premixed turbulent flame contains thin reaction
870
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
zones that strongly impact on the small scale gradients, the constant C should depend on the local flame structure to properly express the scalar dissipation rate. However, its exact determination is far from being straightforward. 2.3.2. Flame surface closure for e vY c A new formulation for the scalar dissipation rate coupling thin reaction zones with turbulent micromixing has recently been proposed [13]. It follows a BML type derivation [21]. The balance equation for c (1 c) is first written. Then, one introduces the fact that the PDF of c is close to a fully bimodal shape, an assumption motivated by the strong Yc gradients imposed by the thin reaction zones. After averaging, the following relation is obtained [13]: 2 _ c qe vc ð17Þ qx 2cm 1 v Y c =ð Ye cb Ye cf Þ2 . The coefficient cm dewith e vc ’ e pends on both the PDF of the progress variable and the FPI database, it may be written R1 cx_ c Pe ðc Þ dc : ð18Þ cm ¼ R0 1 x_ c Pe ðc Þ dc 0
To account for thin reaction zones, the burning rate in the expression of the scalar dissipation rate is cast in [21] _ c ¼ qf e qx S L j rc j;
ð19Þ
where qf is the density of the fresh gases and e S L is the mean characteristic laminar flame speed of the mixture, which is estimated from the burning velocity distribution SL (z) for a given enthalpy Z 1 e S L ðz Þ Pe ðz Þ dz : ð20Þ SL ¼ 0
The generalized flame surface density j rc j ¼ N j rc j [19] is conveniently expressed as a function of the flame wrinkling factor N ¼ j rc j= j rc j. Combining Eq. (17) with Eq. (19) leads to 1 qe v Y c ¼ ð Ye cb Ye cf Þ ð2cm 1Þ qf e S L N j rY c j 2 g Y 002 c : ð21Þ g 002 ð Y c Þmax
its BML value cm 3/4. Eq. (21) has been validated using Direct Numerical Simulation of a premixed turbulent V-flame [13]. 3. RANS computation of a premixed turbulent flame 3.1. Introduction The FPI chemistry tabulation along with the presumed PDF formulation, including both closures for the scalar dissipation rate, have been implemented into the commercial CFD software CFX [22]. The measurements performed by Moreau and Boutier [9] are retained as a test case. The premixed rig consists of a ducted turbulent premixed flame at atmospheric pressure. The duct-shaped combustor has a rectangular cross section. Figure 1 is a sketch of the experimental setup. The main stream is a preheated mixture of air and methane (600 K) at an equivalence ratio of 0.83. It is injected into the combustion chamber at a bulk velocity of the order of 65 m s1. The flame is ignited and stabilized by a jet of hot combustion products, which is injected parallel to the fresh mixture flow under higher velocity conditions, of the order of 130 m s1. The composition of the burned gases is equivalent to the equilibrium state of a stoichiometric methane/air mixture chemically frozen at 2000 K. One of the main characteristics of the device is the high level of turbulence of the flow. Velocity fluctuations result from the rolling up of the shear layer, located between fresh and burned gas inlets, and they are coupled with the strong thermal expansion. An analysis based on measured turbulent kinetic energy and on the estimation of turbulent scales [23] has provided some information on the combustion regime. The characteristic ratio u0 = e S L is about 10. The Kolmogorov microscale is of the order of the lamiThe nar flame thickness, 1 · 104 m. Kolmogorov time is of the order of 1 · 105 s and is smaller than the flame time of the order of 1 · 103 s. Strong interaction between the premixed reaction zone and turbulence is thus expected.
g 002 that vanishes when The term g Y 002 c =ð Y c Þmax g 002 Y c ! 0 makes Eq. (21) a realizable closure. ðg Y 002 Þ , the maximum value reached by g Y 002 , is c
max
c
computed using Eq. (10) with c002 ¼ ðc002 Þmax ¼ cð1 cÞ. Using Eq. (21), the modification of the mixing speed, induced by the premixed flamelets that introduce stiff gradients at the smallest scales, is directly accounted for. In the RANS computations reported below, cm stays close to
Fig. 1. Combustion chamber geometry and boundary conditions.
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
871
3.2. Numerical procedure A set of one-dimensional premixed methane/ air flames has been computed to generate the FPI database. Unity Lewis number is chosen, and the equivalence ratio varies from 0.83 to 1, with the enthalpy level also changed. Flamelets are computed using the PREMIX code from CHEMKIN [24]. The chemical mechanism of Lindstedt [25], including 29 species and 300 reactions, is adopted. This mechanism neglects NOx formation and, under unity Lewis assumption, the mass fraction of nitrogen remains constant across premixed flames. Then, Y z ¼ Y N2 is directly related to the equivalence ratio of the local fuel/air mixture. Preliminary studies have shown that an efficient definition of Yc can be obtained from a linear combination of CO2 and CO mass fractions [12]: Y c ¼ Y CO þ Y CO2 . To save computational time, PDF integrations (Eq. (6)) are fully pre-tabulated inside a look-up table with five entries: 002 ; c; c002 ; and h . ez ; zf n The computational domain has been divided into 12,700 elements. A constant wall temperature of 600 K is used in the simulation. Boundary conditions of temperature and velocity are indicated in Fig. 1, and no-slip conditions with wall-functions [22] were applied for capturing the velocity near the wall. The inlet velocity distribution allows for reproducing the velocity profile located at 39 mm (Figs. 2A and B). The inlet values of the turbulent kinetic energy (k) and its dissipation rate (), are estimated assuming a turbulent intensity of 10% for the main stream and 15% for the hot burnt product inlet, and a characteristic length scale corresponding to the channel height. The standard k model is chosen [22] . 3.3. Results The classical linear relaxation type closure for the scalar dissipation rate is first tested, i.e., Eq. (16) with C = 1. Figures 2A, 3A, 4A, and 5A, respectively, provide transverse profiles of mean axial velocity, temperature, CH4, and O2 mass fraction taken at various streamwise positions, x. RANS computations (lines) and measurements [9,23,26,27] (symbols) are displayed. Transverse profiles located at x < 400 mm are in good agreement with experimental data but discrepancies are observed in the second part of the flow. These errors are observed in the profile of O2 at x = 522 mm, Fig. 5A and that of the temperature at x = 922 mm, Fig. 3A. As the aerodynamic field is strongly related to the thermal expansion, discrepancies are also observed at x = 500 and 650 mm in the streamwise component of the velocity field displayed in Fig. 2A. It is found that the level of fluctuations of progress variable is too high, indicating a very low level of scalar dissipa-
Fig. 2. Transverse profiles of mean axial velocity for different streamwise positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, Simulation. Full lines and squares, x = 39 mm; dashed lines and triangles up, x = 251 mm; dashed dot lines and triangles down, x = 500 mm; and dotted lines and diamonds, x = 650 mm.
tion rate. This is illustrated by isocontours of g Y 002 c _ Y in Fig. 6. Even and of the production term x c when varying the constant C in Eq. (16), it was not possible to reproduce the flame structure. Either the PDF of Yc is fully bimodal, and the burning rate vanishes, or the PDF features a quasi-gaussian shape, and the flame is not properly described. The challenge consists of capturing, through the appropriate level of scalar dissipation rate, the correct quasi-bimodal nature of the PDF. This approximation must include the correct level of the PDF peaks and of the PDF plateau, which appears between fresh and burned gases, and where the chemical burning rates are nonzero. The linear relaxation type closure, which incorporates a turbulent cascade time only, can hardly capture this behavior. This difficulty has been overcome with the flame surface density expression for the scalar dissipation rate (Eq. (21)), in which an additional flame length scale is imbedded [13]. To complete a computation with this type of closure for e vY c , it is necessary to estimate the flame wrinkling S L , where ST is a turbuN ¼ j rc j= j rc j S T = e lent flame speed. Under strong turbulence conditions, i.e., when the ratio u0 =S L is very high, the turbulent flame speed saturates [28] before the local quenching (which is accompanied by the socalled bending effect) can develop [29]. Prediction
872
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
Fig. 3. Transverse profiles of mean temperature for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 42 mm; dashed lines and triangles up, x = 122 mm; dashed dot lines and triangles down, x = 322 mm; and dotted lines and diamonds, x = 922 mm.
Fig. 4. Transverse profiles of mean CH4 mass fraction for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 130 mm; dashed lines and triangles up, x = 330 mm; and dashed dot lines and triangles down, x = 420 mm.
Fig. 5. Transverse profiles of mean O2 mass fraction for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 122 mm; dashed lines and triangles up, x = 522 mm.
Fig. 6. Isocontours of g Y 002 c (A) and of the mean source v Y c is modeled term of Yc expressed in kg m3 s1 (B). e using the linear relaxation closure (Eq. (16) with C = 1).
of the saturation level of the ST curve is very difficult, and both theoretical and empirical models show considerable scattering [29]. For the level of u0 /SL found in the experiment, it can be assumed that the turbulent flame speed has
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
873
4. Conclusions
Fig. 7. Isocontours of g Y 002 c (A) and of the mean source v Y c is modeled term of Yc expressed in kg m3 s1 (B). e using the flame surface expression (Eq. (21) with N = 8).
reached its higher level and then that the wrinkling factor N is almost constant. Following an empirical approach, we have performed various computations varying the wrinkling factor, and it was found that the calculations are weakly sensitive to this parameter when 6 < N < 10. The results presented in this paper have been obtained with N = 8. Figure 7 shows isocontours of g Y 002 c e and of the Y c production term x_ Y c , using relation (21) to estimate the Yc variance dissipation rate. The level of g Y 002 c remains high, with a quasi-bimodal PDF allowing for the computation of the chemical source. Due to constant wall temperature, conductive heat losses are observed near the wall and are coupled with the FPI table through the enthalpy coordinate. This does not affect the main part of the flame but plays an important role on the behavior of the flame near the wall, preventing it from propagating in the boundary layer towards fresh gases. Figure 2B shows the streamwise component of the velocity field that is well estimated. Temperature profiles, Fig. 3B, are also in good agreement with experimental data. Small discrepancies observed at x = 322 mm and at x = 922 mm at the top of the burner may be due to errors in the mean flame thickness. Indeed, at the location x = 922 mm combustion is not completely achieved in the calculation whereas the measurements report a uniform temperature profile. Also note that the equilibrium value of temperature is slightly overestimated by the simulation, probably because of the lack of radiative effects in the computations. Mass fractions of CH4 and O2, respectively, displayed in Figs. 4B and 5B are in good agreement with experimental data.
(1) The FPI method for tabulating chemistry has been coupled with presumed PDFs to perform fast RANS calculations of a premixed turbulent flame. The equivalence ratio of the mixture, the progress of reaction, and the enthalpy are the three control parameters of the chemistry tabulation. (2) A recently proposed closure for the scalar dissipation rate of the progress variable is used. This closure is based on flame surface density and explicitly introduces thin reaction zones in the modeling of turbulent micro-mixing. (3) RANS computations of a ducted premixed turbulent flame developing in a stream of premixed reactant mixing with a stream of hot burnt products are performed. The modeling of the scalar dissipation rate from flame surface density is then evaluated. (4) The heat transfer at the wall is included in the modeling procedure, and detailed comparisons with measurements are reported. Major flame behaviors are recovered, and mean species concentrations are fully reproduced.
References [1] N. Peters, Proc. Combust. Inst. 21 (1986) 1231–1250. [2] Q. Tang, J. Xu, S.B. Pope, Proc. Combust. Inst. 28 (2000) 133–139. [3] A. R. Masri, R. Cao, S.B. Pope, G.M. Goldin, in: Third Mediteranean Combustion Symposium. The Combustion Institute, Marrakech, Morocco, 2003, pp. 1027–1038. [4] T. Landefeld, A. Sadiki, J. Janicka, Flow, Turbulence Combust. 68 (2002) 111–135. [5] S. Carpentier, P. Meunier, F. Aguile´, A new combustion model with detailed chemistry for the prediction of pollutants in gas fired industrial furnaces, Proc. Clean Air, Lisbon, Portugal, 2003. [6] D. Bradley, H. Gaskell, A.K.C. Lau, Proc. Combust. Inst. 23 (1990) 685–692. [7] D. Bradley, H. Gaskell, X.J. Gu, Proc. Combust. Inst. 27 (1998) 1199–1206. [8] K. Bray, M. Champion, P.A. Libby, Combust. Flame 120 (1) (2000) 1–18. [9] P. Moreau, A. Boutier, Proc. Combust. Inst. 16 (1977) 1747–1756. [10] O. Gicquel, De´veloppement dÕune nouvelle me´thode de re´duction des sche´mas cine´tiques: Application au me´thane. Thesis, Ecole Centrale de Paris, Laboratoire EM2C, 1999. [11] O. Gicquel, N. Darabiha, D. The´venin, Proc. Combust. Inst. 28 (2000) 1901–1908. [12] B. Fiorina, R. Baron, O. Gicquel, D. The´venin, S. Carpentier, N. Darabiha, Combust. Theory Modell. 7 (2003) 449–470. [13] L. Vervisch, R. Hauguel, P. Domingo, M. Rullaud, J. Turbulence 5 (4) (2004) 1–36. [14] S.B. Pope, Combust. Theory Modell. 1 (1997) 41–63.
874
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
[15] U. Maas, S.B. Pope, Proc. Combust. Inst. 24 (1992) 103–112. [16] J.A. van Oijen, L.P.H. de Goey, Combust. Sci. Technol. 161 (2000) 113–137. [17] M. Embouazza, O. Gicquel, N. Darabiha, Third Mediteranean Combustion Symposium. The Combustion Institute, Marrakech, Morocco, 2003, pp. 444–454. [18] T. Tura`nyu, Proc. Combust. Inst. 25 (1994) 949–955. [19] D. Veynante, L. Vervisch, Prog. Energy Combust. Sci. 28 (2002) 193–266. [20] L. Vervisch, Study and modeling of finite rate chemistry effects in turbulent non-premixed flames, Annual Research Briefs, CTR, Stanford U., 1992, pp. 411–429. [21] K.N.C. Bray, Proc. Comb. Inst. 26 (1996) 1–26. [22] ANSYS. CFX5.6 user manual.
[23] P. Magre, P. Moreau, G. Collin, R. Borghi, M. Pe´alat, Combust. Flame 71 (2) (1988) 147–168. [24] R.J. Kee, J.F. Grcar, M.D. Smooke, J.A. Miller, A Fortran Program for Modeling a Steady Laminar One-Dimensional Premixed Flame, Report No SAND85-8240, Sandia National Laboratories, 1985. [25] P. Lindstedt, 12 Month Progress Report, Report No TR-96 009, TOPDEC, Brite Euram Program Project, Contract CT 95 0056, 1997. [26] R. Borghi, P. Moreau, C. Bonniot, La recherche ae´rospatiale (3) (1978) 105–116. [27] A.C. Benim, K.J. Syed, Appl. Math. Modell. 22 (1998) 113–136. [28] R. Borghi, M. Champion, Mode´lisation et The´orie des flammes Ed TECHNIP, Paris, France, 2000. [29] T. Poinsot, D. Veynante, in: R.T. Edwards (Ed.), Theoretical and Numerical Combustion, Philadelphia, USA, 2001.
Comments A. Lipatnikov, Chalmers University of Technology, Sweden. What is the difference between your approach and the models developed by Zimont (1977) and by Weller (1993)? These models provide the following equation
@qc @c þ qu V t jrcj þ @=@xi ðqui cÞ ¼ @=@xi qDt @t @xi Reply. In the Zimont modeled equation, complex chemistry effects are included via the ad-hoc expression for Vt the displacement velocity [1]. In our approach, detailed chemistry effects are accounted for using a fully detailed chemical scheme, after using the FPI tabulation. Moreover, this approach brings the possibility of predicting intermediate chemical species and it accounts for partial premixing. The term quVtŒ$cŒ enters the expression for the scalar dissipation rate only, after using a BML type description of micro-mixing in the presence of premixed flamelets; it is not introduced in an equation for the progress of reaction. Therefore the links with the above equation are quite weak.
Reference [1] V.L. Zimont, Exp. Thermal and Fluid Science 21 (2000) 179–186. d
bending of the curve linking the turbulent flame speed ST to velocity fluctuations uÕ. The prediction of the bending point is a critical subject, leading to high discrepancies between different relations proposed in the literature [28 in paper] (both empirical and theoretical). However experiments show that in this limit regime, the wrinkling factor N saturates and is therefore constant in terms of the ratio u0 /SL. Moreover, we have performed various computations by varying the wrinkling factor and it was found that the calculations are weakly sensitive to this parameter when 6 < N < 10. The ultimate state of development of the proposed closure would be to solve a transport equation for the flame surface density. Then, the modeling procedure would be able to handle regimes below the saturation of the turbulent burning velocity versus velocity fluctuations. d
J. H. Chen, Sandia National Laboratories, USA. Did you check the statistical independence of the progress variable, (CO + CO2), and the enthalpy? Reply. The statistical behavior of the progress variable c and the mixture fraction z has been investigated previously [13 in paper] by post-processing measurements of flame D and F performed by Barlow et al. [1]. It appears that the statistical independence is guaranteed by the use of CO and CO2 mass fractions to normalize the progress variable c. Unfortunately, we have neither DNS nor experimental results of flames submitted to sufficient heat losses to check the statistical behavior of all three FPI coordinates: c, z and the enthalpy h. Nevertheless, in a first stage, one can assume that, as it has been observed for the progress variable, the normalization of the enthalpy will minimize its dependence to c and z.
F. Dinkelacker, University of Erlangen, Germany. 1. You closed the scalar dissipation term obviously with a flame surface density term, being proportional to ST/ SL. How did you determine ST? When you assumed this ratio to be constant then this would probably lead to only very special solutions, not being universal. Or otherwise the scalar dissipation term is maybe of low sensitivity?
Reference
Reply. In the considered experiment, the ratio u0 /SL is very high and operates in the regime where occurs the
[1] R.S. Barlow, J.H. Frank, Proc. Combust. Inst. 27 (1998) 1087–1095.
Combustion Institute
Proceedings of the Combustion Institute 30 (2005) 867–874
www.elsevier.com/locate/proci
Premixed turbulent combustion modeling using tabulated detailed chemistry and PDF B. Fiorinaa,b,*, O. Gicquela, L. Vervischc, S. Carpentierb, N. Darabihaa a b
Laboratoire EM2C, UPR 288 CNRS, Ecole Centrale Paris, F92295 Chatenay Malabry, France Poˆle CCMF, Direction de la Recherche, Gaz de France, F93211 La Plaine Saint Denis, France c CORIA, UMR 6614 CNRS, INSA de Rouen, F76800, Saint Etienne du Rouvray, France
Abstract Tabulated chemistry and presumed probability density function (PDF) approaches are combined to perform RANS modeling of premixed turbulent combustion. The chemistry is tabulated from premixed flamelets with three independent parameters: the equivalence ratio of the mixture, the progress of reaction, and the specific enthalpy, to account for heat losses at walls. Mean quantities are estimated from presumed PDFs. This approach is used to numerically predict a turbulent premixed flame diluted by hot burnt products at an equivalence ratio that differs from the main stream of reactants. The investigated flame, subjected to high velocity fluctuations, has a thickened-wrinkled structure. A recently proposed closure for scalar dissipation rate that includes an estimation of the coupling between flame wrinkling and micromixing is retained. Comparisons of simulations with experimental measurements of mean velocity, temperature, and reactants are performed. 2004 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Turbulent combustion; Premixed flame; Modeling; PDF; Tabulated chemistry
1. Introduction The increase of computational resources allows the intensive use of numerical tools to predict chemical structure of flames, understand flame stabilization processes, and help in burner design. Several strategies have been proposed to model and to compute turbulent flames [1–5]. In the context of premixed turbulent combustion, Bradley et al. [6,7] have tabulated detailed chemistry with a collection of one-dimensional premixed laminar flames, then estimating the mean flame structure
*
Corresponding author. Fax: +33 0 1 49 22 53 73. E-mail address: benoit.fi[email protected] (B. Fiorina).
after averaging with presumed probability density function (PDF). Bray et al. [8] have formulated the hypothesis of a fully bimodal PDF of progress variable to compute flames. A combined strategy is adopted in this paper, to perform RANS modeling of a premixed turbulent flame piloted by hot burnt products [9]. Detailed chemistry effects are included in the turbulent combustion model by the means of premixed laminar flamelets, following a recently developed procedure for chemistry tabulation [10–12], called FPI. The chemical flame structure is stored inside a look-up table in terms of three coordinates, which represent the progress of reaction, the local dilution of reactants with products at a different equivalence ratio, and the enthalpy level. PDFs are used to presume
1540-7489/$ - see front matter 2004 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.proci.2004.08.062
868
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
the mean flame structure. Two closures are tested for the scalar dissipation rate of the progress variable, which is one of the fundamental ingredients of the presumed PDFs. A classical linear relaxation formulation and a novel expression are introduced, the latter derives from the Bray–Moss–Libby (BML) analysis to combine information on the flame topology with turbulent micromixing [13]. An ONERA combustion rig that consists of a premixed turbulent flame piloted by hot gases [9] is considered as a test case. This flame is subjected to high velocity fluctuations, leading to strong interactions between chemistry and turbulence.
2. Modeling premixed turbulent combustion 2.1. Chemistry tabulation Predicting the flame shape, its stabilization process, and pollutant emissions in practical burners requires the use of detailed chemistry. Various techniques have been proposed to account for the complexity of chemistry in turbulent flames. ISAT (In situ adaptive tabulation) [14] is based on in situ generation of look-up tables, which are constructed from the direct solving of the time evolution of species concentrations. This method allows for highly accurate simulations of jet [2] and lifted flames [3]. The ILDM methodology [15] is another alternative to chemistry tabulation, which is grounded on a direct mathematical analysis of the dynamic behavior of the non-linear response of the chemical system. Relevant attractive subspaces are constructed from the direct solving of the time evolution of the species concentration. This chemistry tabulation method has been applied to a large variety of problems, including nonpremixed turbulent jet flames [4]. ILDM is an exact description of the chemical behavior of the system. Nevertheless, to provide accurate results with basic hydrocarbon fuels (methane or propane), the dimension of the manifold should at least be around four or five. To combine the advantages of mathematically reduced databases with a good estimation of the whole chemical flame structure, a tabulation strategy called FPI has recently been proposed [11]. FPI is based on physical considerations and is similar to FGM [16]. It determines the reduced chemical subspace effectively accessed in a complex geometry configuration using a set of reference laminar flame computations. In the case of premixed laminar flames, premixed flamelets were used to build the database [10], while auto-ignition computation was retained when the estimation of the ignition delay is crucial [17,18]. The FPI method was also successfully applied to non-adiabatic partially premixed laminar flames [12].
The objective of this paper was to model a premixed turbulent flame diluted by a pilot flame at a different equivalence ratio, this pilot stream containing only hot burnt products. Therefore, the equivalence ratio only weakly varies over the reaction zone surface. There are no reactants in the pilot stream, and the premixed character of the reactive layer is thus mainly preserved. It is then assumed that the chemical flame structure can be reproduced by a collection of one-dimensional laminar premixed flames. A set of laminar premixed flames is computed at various equivalence ratios and for different enthalpy levels, using detailed chemistry. A FPI look-up table is constructed to store all chemical information, with three coordinates: Yc, a progress of reaction; Yz, a mixture fraction characterizing the equivalence ratio; and h, the enthalpy. Yc and Yz are both suitable combinations of species mass fraction, which, respectively, reproduce the progress of reaction and the local mixture of fuel and oxidizer. The third coordinate h, the mixture specific enthalpy, enables to take account for heat losses. Details about the database generation are provided in [12]. All chemical species involved in the detailed kinetics scheme, the temperature, and the thermochemical data are then tabulated in terms of the three coordinates (Yc, Yz, h). If u denotes a species mass fraction or any thermochemical data, then ð1Þ
u ¼ uðY c ; Y z ; hÞ:
This tabulation appears as an interesting tool to perform very fast calculations including complex chemistry effects, specifically when very complex geometries are involved. 2.2. Presumed PDF modeling The statistical properties of the reactive and diffusive layers are easily described from the joint PDF, Pe ðY c ; Y z ; h Þ of the variables chosen to describe the internal flame structure: Z had Z Y c Z Y z f b ~¼ uðY c ; Y z ; h Þ u hmin
Y cf
Y zo
Pe ðY c ; Y z ; h Þ dY c dY z dh ;
ð2Þ
where hmin and had are, respectively, the minimum and the adiabatic values of enthalpy, Y cf and Y cb are, respectively, the values of Yc in fresh and burned gases and Y zo and Y zf correspond to values of Yz in leanest and richest inlets of the physical domain. In fact, the lowest enthalpy limit is the one obtained by computing the chemical equilibrium at T = Tmin, where Tmin is the minimum value of temperature encountered in the turbulent flow. The three coordinates of the FPI database are normalized according to the considered problem, leading to z, hn, and c:
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
z¼
Y z Y zo ; Y zf Y zo
hn ¼
c¼
h hmin ðY z Þ ; had ðY z Þ hmin ðY z Þ
Y c Y cf ðY z ; hÞ : Y cb ðY z ; hÞ Y cf ðY z ; hÞ
ð3Þ
one can readily 002 ; and h : c002 ; ez ; zf n
ð4Þ
c002
ð5Þ
Due to its careful choice and normalization, the progress of reaction c and the mixture fraction z can be considered as two statistically independent variables as mentioned in [20] and further tested from measurements in [13]. It is then possible to introduce specific turbulent model based on this assumption [13]. In the present paper, this model is extended using the normalized enthalpy hn that is also assumed to be statistically independent from other variables. The introduction of the enthalpy as an additional coordinate is motivated by the need to account for flame wall interaction, driving heat transfer to the wall and consequently, the mean turbulent flame shape. Notice that the assumed statistical independence would not be true for the primitive variables Yc, Yz, and h. One may then write for any quantity u (mass fractions or source terms) Z 1Z 1Z 1 ~¼ uðc ; z ; hn Þ u 0
0
0
P ðc Þ Pe ðz ÞP ðhn Þ dc dz dhn ;
ð6Þ
where P ðc Þ and Pe ðz Þ take the form of b-functions that are defined by their first and sec002 Þ, respectively. ond moments ðc; c002 Þ and ðez ; zf Because enthalpy fluctuations are not that large in the considered test case, P ðhn Þ dðhn hn Þ is assumed, and the PDF of the enthalpy becomes a delta function centered on hn . ~ using Estimation of any averaged quantity u relation (6) requires the evaluation of first and second moments of c and z, and of the averaged normalized enthalpy. Starting from relation (5), one can write the conditional average of c at z = z* and hn ¼ hn
ðc j z ¼ z ;hn ¼ hn Þ ¼
ðY c j z ¼ z ;hn ¼ hn Þ Y cf ðz ;hn Þ : Y cb ðz ;hn Þ Y cf ðz ;hn Þ ð7Þ
derive
869
relations
g e 2 f2 ð Ye c Ye cf Þ2 Y 002 c þ Y c Y cf ¼ g ð Ye cb Ye cf Þ2 ðY cb Y cf Þ2 ð Ye c Ye cf Þ qðY cb Y cf ÞY cf 2 ; g ðY cb Y cf Þ2 qY cb qY cf
ez ¼
Ye z Ye zo ; e Y zf Ye zo
002 ¼ zf
hn ¼
g Y 002 z ; ð Ye zf Ye zo Þ2
e he h min : e h ad e h min
for
ð10Þ
ð11Þ
ð12Þ
ð13Þ
2.3. Balance equations In addition to the classical mass, momentum, mixture fraction (mean, Ye z , and fluctuations, g Y 002 z , energy, and RANS k balance equations, Y 002 transport equations for Ye c and g c need to be solved. They may be written as: l o q Ye c _ Y ; þ r ð qe u Ye c Þ ¼ r t r Ye c þ x c ot r l o qg Y 002 t c þ r ð qe ug Y 002 rg Y 002 c Þ ¼ r c ot r l þ 2 t j r Ye c j2 2qD j rY 00c j2 þ 2Y 00c x_ Y c ; r
ð14Þ
ð15Þ
where usual gradient transport is assumed for modeling convection by velocity fluctuations, lt is the eddy viscosity, and r = 0.6 is the turbulent Schmidt _ Y and Y 00 x_ Y number. The chemical source terms x c c c are obtained after averaging the FPI table with the PDFs via Eq. (6). The scalar dissipation rate of the progress of reaction may be write ten q v Y c ¼ qD j rY c j2 ¼ qD j r Ye c j2 þ qDj rY 00c j2 , whose fluctuating part appears in Eq. (15) and requires additional modeling.
As c, z and hn are supposed to be independent ðc j z ¼ z ; hn ¼ hn Þ ¼ c:
ð8Þ
Therefore, integrating relation (7) with Pe ðz Þ leads to the following approximation: c’
Ye c Ye cf : e Y cb Ye cf
ð9Þ
Averaged quantities Ye cf and Ye cb are estimated with relation (6). Following a similar analysis,
2.3.1. Linear relaxation for e vY c Assuming equilibrium between production and dissipation in the balance equation for the scalar dissipation rate, the widely used linear relaxation type closure appears [19] 002 e q g ð16Þ Y ; q v Y c ¼ C k c where C is a modeling constant. Because the premixed turbulent flame contains thin reaction
870
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
zones that strongly impact on the small scale gradients, the constant C should depend on the local flame structure to properly express the scalar dissipation rate. However, its exact determination is far from being straightforward. 2.3.2. Flame surface closure for e vY c A new formulation for the scalar dissipation rate coupling thin reaction zones with turbulent micromixing has recently been proposed [13]. It follows a BML type derivation [21]. The balance equation for c (1 c) is first written. Then, one introduces the fact that the PDF of c is close to a fully bimodal shape, an assumption motivated by the strong Yc gradients imposed by the thin reaction zones. After averaging, the following relation is obtained [13]: 2 _ c qe vc ð17Þ qx 2cm 1 v Y c =ð Ye cb Ye cf Þ2 . The coefficient cm dewith e vc ’ e pends on both the PDF of the progress variable and the FPI database, it may be written R1 cx_ c Pe ðc Þ dc : ð18Þ cm ¼ R0 1 x_ c Pe ðc Þ dc 0
To account for thin reaction zones, the burning rate in the expression of the scalar dissipation rate is cast in [21] _ c ¼ qf e qx S L j rc j;
ð19Þ
where qf is the density of the fresh gases and e S L is the mean characteristic laminar flame speed of the mixture, which is estimated from the burning velocity distribution SL (z) for a given enthalpy Z 1 e S L ðz Þ Pe ðz Þ dz : ð20Þ SL ¼ 0
The generalized flame surface density j rc j ¼ N j rc j [19] is conveniently expressed as a function of the flame wrinkling factor N ¼ j rc j= j rc j. Combining Eq. (17) with Eq. (19) leads to 1 qe v Y c ¼ ð Ye cb Ye cf Þ ð2cm 1Þ qf e S L N j rY c j 2 g Y 002 c : ð21Þ g 002 ð Y c Þmax
its BML value cm 3/4. Eq. (21) has been validated using Direct Numerical Simulation of a premixed turbulent V-flame [13]. 3. RANS computation of a premixed turbulent flame 3.1. Introduction The FPI chemistry tabulation along with the presumed PDF formulation, including both closures for the scalar dissipation rate, have been implemented into the commercial CFD software CFX [22]. The measurements performed by Moreau and Boutier [9] are retained as a test case. The premixed rig consists of a ducted turbulent premixed flame at atmospheric pressure. The duct-shaped combustor has a rectangular cross section. Figure 1 is a sketch of the experimental setup. The main stream is a preheated mixture of air and methane (600 K) at an equivalence ratio of 0.83. It is injected into the combustion chamber at a bulk velocity of the order of 65 m s1. The flame is ignited and stabilized by a jet of hot combustion products, which is injected parallel to the fresh mixture flow under higher velocity conditions, of the order of 130 m s1. The composition of the burned gases is equivalent to the equilibrium state of a stoichiometric methane/air mixture chemically frozen at 2000 K. One of the main characteristics of the device is the high level of turbulence of the flow. Velocity fluctuations result from the rolling up of the shear layer, located between fresh and burned gas inlets, and they are coupled with the strong thermal expansion. An analysis based on measured turbulent kinetic energy and on the estimation of turbulent scales [23] has provided some information on the combustion regime. The characteristic ratio u0 = e S L is about 10. The Kolmogorov microscale is of the order of the lamiThe nar flame thickness, 1 · 104 m. Kolmogorov time is of the order of 1 · 105 s and is smaller than the flame time of the order of 1 · 103 s. Strong interaction between the premixed reaction zone and turbulence is thus expected.
g 002 that vanishes when The term g Y 002 c =ð Y c Þmax g 002 Y c ! 0 makes Eq. (21) a realizable closure. ðg Y 002 Þ , the maximum value reached by g Y 002 , is c
max
c
computed using Eq. (10) with c002 ¼ ðc002 Þmax ¼ cð1 cÞ. Using Eq. (21), the modification of the mixing speed, induced by the premixed flamelets that introduce stiff gradients at the smallest scales, is directly accounted for. In the RANS computations reported below, cm stays close to
Fig. 1. Combustion chamber geometry and boundary conditions.
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
871
3.2. Numerical procedure A set of one-dimensional premixed methane/ air flames has been computed to generate the FPI database. Unity Lewis number is chosen, and the equivalence ratio varies from 0.83 to 1, with the enthalpy level also changed. Flamelets are computed using the PREMIX code from CHEMKIN [24]. The chemical mechanism of Lindstedt [25], including 29 species and 300 reactions, is adopted. This mechanism neglects NOx formation and, under unity Lewis assumption, the mass fraction of nitrogen remains constant across premixed flames. Then, Y z ¼ Y N2 is directly related to the equivalence ratio of the local fuel/air mixture. Preliminary studies have shown that an efficient definition of Yc can be obtained from a linear combination of CO2 and CO mass fractions [12]: Y c ¼ Y CO þ Y CO2 . To save computational time, PDF integrations (Eq. (6)) are fully pre-tabulated inside a look-up table with five entries: 002 ; c; c002 ; and h . ez ; zf n The computational domain has been divided into 12,700 elements. A constant wall temperature of 600 K is used in the simulation. Boundary conditions of temperature and velocity are indicated in Fig. 1, and no-slip conditions with wall-functions [22] were applied for capturing the velocity near the wall. The inlet velocity distribution allows for reproducing the velocity profile located at 39 mm (Figs. 2A and B). The inlet values of the turbulent kinetic energy (k) and its dissipation rate (), are estimated assuming a turbulent intensity of 10% for the main stream and 15% for the hot burnt product inlet, and a characteristic length scale corresponding to the channel height. The standard k model is chosen [22] . 3.3. Results The classical linear relaxation type closure for the scalar dissipation rate is first tested, i.e., Eq. (16) with C = 1. Figures 2A, 3A, 4A, and 5A, respectively, provide transverse profiles of mean axial velocity, temperature, CH4, and O2 mass fraction taken at various streamwise positions, x. RANS computations (lines) and measurements [9,23,26,27] (symbols) are displayed. Transverse profiles located at x < 400 mm are in good agreement with experimental data but discrepancies are observed in the second part of the flow. These errors are observed in the profile of O2 at x = 522 mm, Fig. 5A and that of the temperature at x = 922 mm, Fig. 3A. As the aerodynamic field is strongly related to the thermal expansion, discrepancies are also observed at x = 500 and 650 mm in the streamwise component of the velocity field displayed in Fig. 2A. It is found that the level of fluctuations of progress variable is too high, indicating a very low level of scalar dissipa-
Fig. 2. Transverse profiles of mean axial velocity for different streamwise positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, Simulation. Full lines and squares, x = 39 mm; dashed lines and triangles up, x = 251 mm; dashed dot lines and triangles down, x = 500 mm; and dotted lines and diamonds, x = 650 mm.
tion rate. This is illustrated by isocontours of g Y 002 c _ Y in Fig. 6. Even and of the production term x c when varying the constant C in Eq. (16), it was not possible to reproduce the flame structure. Either the PDF of Yc is fully bimodal, and the burning rate vanishes, or the PDF features a quasi-gaussian shape, and the flame is not properly described. The challenge consists of capturing, through the appropriate level of scalar dissipation rate, the correct quasi-bimodal nature of the PDF. This approximation must include the correct level of the PDF peaks and of the PDF plateau, which appears between fresh and burned gases, and where the chemical burning rates are nonzero. The linear relaxation type closure, which incorporates a turbulent cascade time only, can hardly capture this behavior. This difficulty has been overcome with the flame surface density expression for the scalar dissipation rate (Eq. (21)), in which an additional flame length scale is imbedded [13]. To complete a computation with this type of closure for e vY c , it is necessary to estimate the flame wrinkling S L , where ST is a turbuN ¼ j rc j= j rc j S T = e lent flame speed. Under strong turbulence conditions, i.e., when the ratio u0 =S L is very high, the turbulent flame speed saturates [28] before the local quenching (which is accompanied by the socalled bending effect) can develop [29]. Prediction
872
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
Fig. 3. Transverse profiles of mean temperature for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 42 mm; dashed lines and triangles up, x = 122 mm; dashed dot lines and triangles down, x = 322 mm; and dotted lines and diamonds, x = 922 mm.
Fig. 4. Transverse profiles of mean CH4 mass fraction for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 130 mm; dashed lines and triangles up, x = 330 mm; and dashed dot lines and triangles down, x = 420 mm.
Fig. 5. Transverse profiles of mean O2 mass fraction for different axial positions. (A) e v Y c is modeled using the linear relaxation closure (Eq. (16) with C = 1); (B) e v Y c is modeled using the flame surface expression (Eq. (21) with N = 8). Symbols, experimental data; lines, simulation. Full lines and squares, x = 122 mm; dashed lines and triangles up, x = 522 mm.
Fig. 6. Isocontours of g Y 002 c (A) and of the mean source v Y c is modeled term of Yc expressed in kg m3 s1 (B). e using the linear relaxation closure (Eq. (16) with C = 1).
of the saturation level of the ST curve is very difficult, and both theoretical and empirical models show considerable scattering [29]. For the level of u0 /SL found in the experiment, it can be assumed that the turbulent flame speed has
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
873
4. Conclusions
Fig. 7. Isocontours of g Y 002 c (A) and of the mean source v Y c is modeled term of Yc expressed in kg m3 s1 (B). e using the flame surface expression (Eq. (21) with N = 8).
reached its higher level and then that the wrinkling factor N is almost constant. Following an empirical approach, we have performed various computations varying the wrinkling factor, and it was found that the calculations are weakly sensitive to this parameter when 6 < N < 10. The results presented in this paper have been obtained with N = 8. Figure 7 shows isocontours of g Y 002 c e and of the Y c production term x_ Y c , using relation (21) to estimate the Yc variance dissipation rate. The level of g Y 002 c remains high, with a quasi-bimodal PDF allowing for the computation of the chemical source. Due to constant wall temperature, conductive heat losses are observed near the wall and are coupled with the FPI table through the enthalpy coordinate. This does not affect the main part of the flame but plays an important role on the behavior of the flame near the wall, preventing it from propagating in the boundary layer towards fresh gases. Figure 2B shows the streamwise component of the velocity field that is well estimated. Temperature profiles, Fig. 3B, are also in good agreement with experimental data. Small discrepancies observed at x = 322 mm and at x = 922 mm at the top of the burner may be due to errors in the mean flame thickness. Indeed, at the location x = 922 mm combustion is not completely achieved in the calculation whereas the measurements report a uniform temperature profile. Also note that the equilibrium value of temperature is slightly overestimated by the simulation, probably because of the lack of radiative effects in the computations. Mass fractions of CH4 and O2, respectively, displayed in Figs. 4B and 5B are in good agreement with experimental data.
(1) The FPI method for tabulating chemistry has been coupled with presumed PDFs to perform fast RANS calculations of a premixed turbulent flame. The equivalence ratio of the mixture, the progress of reaction, and the enthalpy are the three control parameters of the chemistry tabulation. (2) A recently proposed closure for the scalar dissipation rate of the progress variable is used. This closure is based on flame surface density and explicitly introduces thin reaction zones in the modeling of turbulent micro-mixing. (3) RANS computations of a ducted premixed turbulent flame developing in a stream of premixed reactant mixing with a stream of hot burnt products are performed. The modeling of the scalar dissipation rate from flame surface density is then evaluated. (4) The heat transfer at the wall is included in the modeling procedure, and detailed comparisons with measurements are reported. Major flame behaviors are recovered, and mean species concentrations are fully reproduced.
References [1] N. Peters, Proc. Combust. Inst. 21 (1986) 1231–1250. [2] Q. Tang, J. Xu, S.B. Pope, Proc. Combust. Inst. 28 (2000) 133–139. [3] A. R. Masri, R. Cao, S.B. Pope, G.M. Goldin, in: Third Mediteranean Combustion Symposium. The Combustion Institute, Marrakech, Morocco, 2003, pp. 1027–1038. [4] T. Landefeld, A. Sadiki, J. Janicka, Flow, Turbulence Combust. 68 (2002) 111–135. [5] S. Carpentier, P. Meunier, F. Aguile´, A new combustion model with detailed chemistry for the prediction of pollutants in gas fired industrial furnaces, Proc. Clean Air, Lisbon, Portugal, 2003. [6] D. Bradley, H. Gaskell, A.K.C. Lau, Proc. Combust. Inst. 23 (1990) 685–692. [7] D. Bradley, H. Gaskell, X.J. Gu, Proc. Combust. Inst. 27 (1998) 1199–1206. [8] K. Bray, M. Champion, P.A. Libby, Combust. Flame 120 (1) (2000) 1–18. [9] P. Moreau, A. Boutier, Proc. Combust. Inst. 16 (1977) 1747–1756. [10] O. Gicquel, De´veloppement dÕune nouvelle me´thode de re´duction des sche´mas cine´tiques: Application au me´thane. Thesis, Ecole Centrale de Paris, Laboratoire EM2C, 1999. [11] O. Gicquel, N. Darabiha, D. The´venin, Proc. Combust. Inst. 28 (2000) 1901–1908. [12] B. Fiorina, R. Baron, O. Gicquel, D. The´venin, S. Carpentier, N. Darabiha, Combust. Theory Modell. 7 (2003) 449–470. [13] L. Vervisch, R. Hauguel, P. Domingo, M. Rullaud, J. Turbulence 5 (4) (2004) 1–36. [14] S.B. Pope, Combust. Theory Modell. 1 (1997) 41–63.
874
B. Fiorina et al. / Proceedings of the Combustion Institute 30 (2005) 867–874
[15] U. Maas, S.B. Pope, Proc. Combust. Inst. 24 (1992) 103–112. [16] J.A. van Oijen, L.P.H. de Goey, Combust. Sci. Technol. 161 (2000) 113–137. [17] M. Embouazza, O. Gicquel, N. Darabiha, Third Mediteranean Combustion Symposium. The Combustion Institute, Marrakech, Morocco, 2003, pp. 444–454. [18] T. Tura`nyu, Proc. Combust. Inst. 25 (1994) 949–955. [19] D. Veynante, L. Vervisch, Prog. Energy Combust. Sci. 28 (2002) 193–266. [20] L. Vervisch, Study and modeling of finite rate chemistry effects in turbulent non-premixed flames, Annual Research Briefs, CTR, Stanford U., 1992, pp. 411–429. [21] K.N.C. Bray, Proc. Comb. Inst. 26 (1996) 1–26. [22] ANSYS. CFX5.6 user manual.
[23] P. Magre, P. Moreau, G. Collin, R. Borghi, M. Pe´alat, Combust. Flame 71 (2) (1988) 147–168. [24] R.J. Kee, J.F. Grcar, M.D. Smooke, J.A. Miller, A Fortran Program for Modeling a Steady Laminar One-Dimensional Premixed Flame, Report No SAND85-8240, Sandia National Laboratories, 1985. [25] P. Lindstedt, 12 Month Progress Report, Report No TR-96 009, TOPDEC, Brite Euram Program Project, Contract CT 95 0056, 1997. [26] R. Borghi, P. Moreau, C. Bonniot, La recherche ae´rospatiale (3) (1978) 105–116. [27] A.C. Benim, K.J. Syed, Appl. Math. Modell. 22 (1998) 113–136. [28] R. Borghi, M. Champion, Mode´lisation et The´orie des flammes Ed TECHNIP, Paris, France, 2000. [29] T. Poinsot, D. Veynante, in: R.T. Edwards (Ed.), Theoretical and Numerical Combustion, Philadelphia, USA, 2001.
Comments A. Lipatnikov, Chalmers University of Technology, Sweden. What is the difference between your approach and the models developed by Zimont (1977) and by Weller (1993)? These models provide the following equation
@qc @c þ qu V t jrcj þ @=@xi ðqui cÞ ¼ @=@xi qDt @t @xi Reply. In the Zimont modeled equation, complex chemistry effects are included via the ad-hoc expression for Vt the displacement velocity [1]. In our approach, detailed chemistry effects are accounted for using a fully detailed chemical scheme, after using the FPI tabulation. Moreover, this approach brings the possibility of predicting intermediate chemical species and it accounts for partial premixing. The term quVtŒ$cŒ enters the expression for the scalar dissipation rate only, after using a BML type description of micro-mixing in the presence of premixed flamelets; it is not introduced in an equation for the progress of reaction. Therefore the links with the above equation are quite weak.
Reference [1] V.L. Zimont, Exp. Thermal and Fluid Science 21 (2000) 179–186. d
bending of the curve linking the turbulent flame speed ST to velocity fluctuations uÕ. The prediction of the bending point is a critical subject, leading to high discrepancies between different relations proposed in the literature [28 in paper] (both empirical and theoretical). However experiments show that in this limit regime, the wrinkling factor N saturates and is therefore constant in terms of the ratio u0 /SL. Moreover, we have performed various computations by varying the wrinkling factor and it was found that the calculations are weakly sensitive to this parameter when 6 < N < 10. The ultimate state of development of the proposed closure would be to solve a transport equation for the flame surface density. Then, the modeling procedure would be able to handle regimes below the saturation of the turbulent burning velocity versus velocity fluctuations. d
J. H. Chen, Sandia National Laboratories, USA. Did you check the statistical independence of the progress variable, (CO + CO2), and the enthalpy? Reply. The statistical behavior of the progress variable c and the mixture fraction z has been investigated previously [13 in paper] by post-processing measurements of flame D and F performed by Barlow et al. [1]. It appears that the statistical independence is guaranteed by the use of CO and CO2 mass fractions to normalize the progress variable c. Unfortunately, we have neither DNS nor experimental results of flames submitted to sufficient heat losses to check the statistical behavior of all three FPI coordinates: c, z and the enthalpy h. Nevertheless, in a first stage, one can assume that, as it has been observed for the progress variable, the normalization of the enthalpy will minimize its dependence to c and z.
F. Dinkelacker, University of Erlangen, Germany. 1. You closed the scalar dissipation term obviously with a flame surface density term, being proportional to ST/ SL. How did you determine ST? When you assumed this ratio to be constant then this would probably lead to only very special solutions, not being universal. Or otherwise the scalar dissipation term is maybe of low sensitivity?
Reference
Reply. In the considered experiment, the ratio u0 /SL is very high and operates in the regime where occurs the
[1] R.S. Barlow, J.H. Frank, Proc. Combust. Inst. 27 (1998) 1087–1095.