- Email: [email protected]
A Filtered Tabulated Chemistry model for LES of stratified flames
Combustion and Flame 159 (2012) 2704–2717
Contents lists available at SciVerse ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
A Filtered Tabulated Chemistry model for LES of stratified flames Pierre Auzillon, Olivier Gicquel, Nasser Darabiha, Denis Veynante, Benoit Fiorina ⇑ CNRS, UPR 288, Laboratoire d’Energétique Moléculaire et Macroscopique, Combustion (EM2C), Grande Voie des Vignes, 92290 Châtenay-Malabry, France Ecole Centrale Paris, Grande Voie des Vignes, 92290 Châtenay-Malabry, France
a r t i c l e
i n f o
Article history: Available online 23 April 2012 Keywords: Large Eddy Simulation Stratified combustion Tabulated chemistry Turbulent combustion modeling
a b s t r a c t The model called Filtered Tabulated Chemistry for LES (F-TACLES), recently developed for turbulent premixed combustion [1], is extended in the present article to stratified flames where the equivalence ratio is not spatially uniform. The principle is to generate a look-up table where thermo-chemical variables, computed from a set of filtered laminar unstrained premixed flamelets, are stored as a function of the ~ c , the filtered mixture fraction ~c, and the mixture fraction subgrid scale varifiltered progress variable Y 002 for a given value of the LES filter size D. From this database and using a classical subgrid scale ance zf flame surface wrinkling model, a closure model for the filtered progress variable balance equation is then 002 ) are added to the Navier–Stokes ~ c , ~z and zf proposed. As only three balance equations (respectively for Y system, this turbulent combustion model is straightforward to implement and is not expensive in terms of CPU time. After preliminary tests on 1-D planar flames, the proposed strategy is applied to a 3-D simulation of a swirled turbulent stratified flame. The agreements observed between numerical and measured data validate the modeling strategy. Finally, it is observed that the flame position cannot be captured without accounting for the subgrid scale mixture fraction heterogeneities that impact on the flame propagation. This result highlights the importance of the mixture fraction filtered density function on the filtered chemical reaction rate. Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Future aeronautical engines will introduce new combustion chambers based on lean-premixed, multi-injection, and/or prevaporized technologies. These systems promote combustion in lean premixed-like regimes to control the flame temperature and consequently the pollutant emissions such as nitrogen oxides [2]. For safety reasons, fuel and oxidizer are injected separately before entering the combustion chamber; therefore, in practice reactants are not perfectly mixed (non-spatially uniform equivalence ratio) and stratified combustion regimes are observed. The emergence of these complex advanced combustion devices introduces new challenges in terms of modeling [3]. A crucial issue is to capture the correct propagation speed of the flame front on grid coarser than the flame thickness. To overcome this difficulty, efforts have recently been made for fully premixed situations [1,4–8]. The popular thickened flame model (TFLES) initially developed by Colin et al. [8] is robust and can easily be applied to realistic combustion chambers [9,10]. However, a too large flame front thickening may affect the flame response to unsteady flow motions [11]. Another approach is to solve a smooth scalar field where a given iso-surface is identified to the instantaneous flame front posi⇑ Corresponding author at: Ecole Centrale Paris, Grande Voie des Vignes, 92290 Châtenay-Malabry, France. Fax: +33 1 47028035. E-mail address: benoit.fi[email protected] (B. Fiorina).
tion. In such technique, called G-equation model, the inner layer is tracked using a level-set technique. Initially developed in a RANS context [12], the G-equation has been reformulated for LES [4,7,13–15]. These methods are very efficient to capture the flame front propagation but, as pointed out in Refs. [16,17], the coupling between the flame front position and flow equations through heat expansion is very challenging. A recently-proposed alternative technique called Filtered Tabulated Chemistry for LES (F-TACLES) [1] consists in generating a chemical look-up database from filtered 1-D unstrained laminar premixed flames computed with detailed chemical schemes. This approach, straightforward to implement in a Navier–Stokes solver, has shown very promising results on turbulent premixed flames and provides a more accurate prediction of the DNS flame dynamics than the TFLES approach [11]. Unfortunately, in practical situations, the mixing between fuel and oxidizer is never perfect and the flame front propagates in an heterogeneous environment. In an LES context, the equivalence ratio fluctuations at the subgrid scale will affect the propagation of the filtered flame front. The G-equation formalism is well adapted to LES of premixed turbulent combustion because it is designed to accurately describe flame front displacements [7,15,18]. Interesting attempts have been made to perform LES of stratified systems with this approach [19,20], but modeling efforts are still required to apply this strategy in industrial combustors, especially when local exinctions occur.
0010-2180/$ - see front matter Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2012.03.006
2705
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Recent developments of the Linear Eddy Mixing (LEM) model combined with an Artificial Neural Network (ANN) enable LES of both premixed and non premixed systems with a skeletal chemical mechanism [21]. In the context of tabulated chemistry, the FPI-PCM [22] method accounts for mixture fraction subgrid scale fluctuations with the introduction of a presumed b filtered density function for the mixture fraction. However, this approach does not capture the correct propagation speed of the flame front. In particular, the laminar flame speed is not recovered when turbulence vanishes [1]. Finally, effects of subgrid scale stratification on the flame propagations are not predicted by the TFLES model that does not account for equivalence ratio fluctuations in chemical reaction rates closures. The objective of this paper is to draw a modeling route for LES of stratified flames, where the flame front propagates into nonperfectly mixed reactants. In the first section, the modeling issues for LES of stratified flames are presented. A tabulated chemistry formalism adapted to stratified combustion is introduced, and the propagation of a filtered flame front into a stratified environment is discussed. In a second section, the F-TACLES model initially developed for premixed combustion is extended to turbulent stratified flames. Numerical simulations of planar laminar flames that propagate in an heterogeneous equivalence ratio flow field are presented. Finally, 3-D LES of the experimental MOLECULES combustor where a swirled flame is stabilized behind a TURBOMECA injector are performed. Simulation results are compared with experimental data obtained by Janus et al. [23–25]. The mathematical derivation of the F-TACLES formalism is detailed in Appendix A while Appendix B is devoted to the practical implementation of tabulated chemistry. 2. Modeling issues for Large Eddy Simulation of stratified flames 2.1. Chemistry tabulation for stratified combustion A careful description of chemical reactions is required to accurately predict the flame structure and the pollutant formation. As hydrocarbon combustion generally involves hundreds of species, simulations of practical combustion systems including complex chemical schemes are not possible. Various studies have been performed during the last decades to account for detailed chemistry effects for a reasonable computational cost. One of the most popular techniques is the use of tabulated chemistry. Maas and Pope [26] have shown that in practical reactive systems chemical composition mainly evolves on a reduced subspace called manifold. Different strategies based on mathematical analysis such as Intrinsic Low Dimensional Manifold (ILDM) [26], or based on physical considerations like Flame Prolongation of ILDM (FPI) [27,28], or Flamelet Generated Manifold [29], or using artificial neural networks [21] have been suggested to identify these manifolds. For each of these approaches, thermochemical quantities are stored in a multi-dimensional look-up database. These techniques are able to predict regimes such as perfectly premixed or stratified combustion [30]. Recently, Nguyen et al. [31] have extended this formalism to diffusion flames. In this article, the FPI approach [27,28] is retained. The chemical subspaces accessed by a stratified flame in a complex geometry configuration is mapped by a collection of 1-D laminar premixed flames computed for various equivalence ratio within the flammability limits using detailed chemistry. The rich (lean) flammability limit is defined as the maximal (minimal) possible value of the equivalence ratio of a premixed flame. Above (bellow) that point, a freely propagating premixed planar laminar flame will extinguish. Note that this formalism excludes the possible combustion of reactant mixtures below or above flammability limits sustained by hot burnt gases issued from the combustion of flammable mix-
tures. In the present study, laminar thermochemical quantities are stored in a look-up table as a function of two coordinates: the progress variable Yc that evolves monotonically between fresh and burnt gases and a mixture fraction z equal to zero in pure oxidizer and unity in pure fuel. For methane–air combustion, preliminary studies have shown that an efficient definition of this progress variable is obtained from a linear combination of CO and CO2 mass fractions: Y c ¼ Y CO þ Y CO2 [28]. Such a definition of Yc ensures a one-to-one correspondence between species mass fractions and the progress variable. A normalized form of the progress variable c is introduced:
c¼
Y c Y fc
ð1Þ
Y bc Y fc
where Y fc and Y bc are the values of progress variable in fresh and burnt gases, respectively. The mixture fraction z is defined from Yz, an inert species or a linear combination of atomic elements representative of fuel and air mixing:
z¼
Y z Y ox z
ð2Þ
Y fz Y ox z
where Y fz and Y ox z denote values of Yz in fuel stream and oxidizer stream, respectively. 2.2. Propagation of a filtered flame front in a stratified mixture 2.2.1. Premixed flame front (uniform equivalence ratio) A premixed wrinkled reaction layer, propagating in perfectly mixed reactants (spatially uniform equivalence ratio), is first considered and illustrated in Fig. 1a. The solid black line represents the interface between fresh and burnt gases before applying the LES filter. The wide gray area shows the corresponding filtered flame front. As demonstrated in Appendix A, when the flame structure is assumed to be similar to the structure of a planar filtered flame, the progress variable reaction rate reads [1]:
e_ ¼ N q x Yc
Z
þ1
1
q x_ Y c x0n GD ðxn x0n Þdx0n
ð3Þ
where the subgrid scale flame surface wrinkling factor N ¼ jrcj=jrcj measures the ratio of the total flame surface to the resolved flame surface in the filter volume. The ⁄ superscript denotes quantities issued from computations of 1-D laminar unstrained premixed flames computed as function of the direction normal to the flame front, xn, for a constant equivalence ratio (or mixture fraction z). GD is a 1-D filter function of size D. The propagation speed SD of the filtered flame front is obtained by integrating the filtered reaction rate given in Eq. (3) along the normal direction to the reaction layer:
q0 SD ¼
Z
1 Y bc Y fc
þ1
Z
þ1
N
1
1
q x_ Y c x0n GD xn x0n dx0n dxn
ð4Þ
where q0 is the fresh gas density. This equation leads to, assuming that the flame wrinkling factor N is uniform over the flame brush:
SD ¼ NSl
ð5Þ
where Sl is the flame speed of a planar laminar premixed flame front:
q0 Sl ¼
Z
1 Y bc
Y fc
þ1
1
q x_ Y c ðxn Þdxn
ð6Þ
2.2.2. Stratified flame front (non-uniform equivalence ratio) The propagation of a flame front in a stratified environment where the equivalence ratio is not spatially uniform is illustrated in Fig. 1b where the dashed lines represent the iso-mixture fraction
2706
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Fig. 1. Propagation of a filtered flame front.
lines. It has been shown in Refs. [30,31] that the orientation of progress variable isosurfaces with mixture fraction isosurfaces, related to the cross-scalar dissipation rate vzc = qDrc rz, may affect the chemical flame structure. However, to take into account this effect in practical situations is not straightforward and considerably complicates the chemical look-up table generation. In order to obtain an easy-to-implement modeling strategy relevant for practical applications, the orientation of iso-c versus iso-z surfaces is assumed constant and perpendicular: vzc = 0. As indicated in Appendix B, this simplification supposes that the normalized form of the progress variable, c, and the mixture fraction, z, are independent variables at the subgrid scale level (unlike to Yc and z which are not independent, as the equilibrium value of Yc depends on the mixture fraction z [32]). Under these assumptions, Appendix A shows that the filtered reaction rate may be written as follows:
e_ ¼ N q x Yc
Z
1
0
_ Y c jz ¼ z0 iPðz0 Þdz hqx
ð7Þ
0
where P(z0 ) is the filtered density function of the mixture fraction _ c jz ¼ z0 i, the conditional filtered value of the reaction rate and hqx at z = z0 , is modeled by:
_ Y c jz ¼ z0 i ¼ hqx
Z
þ1
1
q x_ Y c
0 x0n ; z0 GD ðxn x0n Þdxn
ð8Þ
ð9Þ
where
1 e Sl ¼ q 0
Z
1 0
q0 ðz0 ÞSl ðz0 ÞPðz0 Þdz0
where q is the density, u the velocity vector, and D the diffusivity. x_ Y c is the progress variable reaction rate expressed in s1. Closure models have been proposed in [1] for the subgrid scale transport ~c q fc Þ, the filtered laminar diffusion term ~Y uY u terms r ðq e_ Y . This methodology is here qDrY c and the filtered source term x c extended to account for mixture fraction heterogeneities at the subgrid scale level. e_ Y 3.1. Filtered chemical reaction source term x c e_ Y reads as As shown in Appendix A the expression for x c follows:
1
e_ ¼ N x Yc q
Z
1
ð10Þ
where q0(z), the fresh gas density, is function of the mixture fracR 0 ¼ 01 q0 ðz0 ÞPðz0 Þdz0 . tion and q0 is then expressed by q 3. F-TACLES for stratified flames The balance equation for the filtered normalized progress variable ~c contains several cross-terms introduced by the dependency of Y bc and Y fc on the mixture fraction and is difficult to close in practice. Solving a transport equation for the filtered non-normalized e c is then preferred: progress variable Y
0
_ Y c jz ¼ z0 iPðz0 Þdz hqx
ð12Þ
0
_ Y c jz ¼ z0 i is estimated by filtering a 1-D laminar prewhere hqx mixed flamelet computed at a given mixture fraction z:
_ Y c jz ¼ z0 i ¼ hqx
Then, under these modeling assumptions, the propagation speed of a filtered flame front in a stratified mixture is given by:
SD ¼ N e Sl
fc Y @q e_ gc q fc Þ ¼ r ðqDrY c Þ r ðq fc Þ þ q ~Y ~Y u uY u x þ r ðq Yc @t ð11Þ
Z
þ1
1
q x_ Y c x0n ; z0 GD xn x0n dx0n
ð13Þ
In this article, P(z) is modeled by a b-function parametrized by the 002 . first- and second-order moments of the mixture fraction: ~z and zf e _ The filtered reaction rate x Y c is in practice pre-computed and e_ Y ¼ stored in a four-dimensional look-up database, x c e f f 002 x_ Y c ½ Y c ; ~z; Sz ; D where Sz ¼ z =ð~zð1 ~zÞÞ is the unmixedness factor, 002 ¼ 0) to unity going from zero in perfectly mixed situations ( zf when reactants are fully segregated (bimodal distribution z = 0 or z = 1). 3.2. Filtered laminar diffusion terms A model for the filtered molecular diffusion term has been discussed in [1] for premixed flames. The extension toward stratified combustion reads as follows:
fc Þ r ðqDrY c Þ ¼ r ðaY c q0 D0 r Y
ð14Þ
where q0 and D0 are reference values for the density and the molecular diffusion coefficient, respectively, for example, corresponding to fresh gas values. The correction factor aY c is estimated from 1D filtered premixed flames:
2707
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
aY c ½ Yfc ; ez ; Sz ; D ¼
c q D @Y @xn
ð15Þ
q0 D0 @@xYecn
e c =Y fc eq in Fig. 2 for Sz = 0 (perfect mixing) and a function of ~z and Y Sz = 0.02, respectively. 3.3. Unresolved transport term
The laminar diffusivity is computed assuming unity Lewis number, D ¼ k =q C p . In Eq. (15), q D @Y c =@xn is computed as follows:
q
@Y c
D
@xn
¼
Z
1
hq D
0
@Y c
@xn
0
0
0
jz ¼ z iPðz Þdz
ð16Þ
fc =@xn presents in the RHS denominator of The computation of @ Y Eq. (15), is decomposed as follows:
fc @Y @ q Y c ¼ @xn @xn q
!
The unresolved transport term is modeled as [1]:
1 @ q Y c f @ q Yc ¼ q @xn @xn
! ð17Þ
where
Z 1 @ q Y c @Y ¼ q c jz ¼ z0 Pðz0 Þdz0 @xn @xn 0 Z 1 @ q @q 0 ¼ jz ¼ z0 Pðz0 Þdz @xn @xn 0
ð18Þ ð19Þ
fc ; ~z; Sz ; D is As for the progress variable reaction rate, aY c ¼ aY c ½ Y stored in the look-up database. This correction factor is plotted as
gc q fc Þ fc ¼ NXY þ r ðaY ðN 1Þq D0 r Y ~Y uY u r q c c 0
ð20Þ
where the first RHS term, NXY c , corresponds to the thermal expansion and is estimated from 1-D detailed chemistry laminar premixed flame solutions as [1]:
XY c ¼ q0 ðzÞSl ðzÞ
fc @Y c @Y q0 ðzÞSl ðzÞ @xn @xn
ð21Þ
while the second RHS term in Eq. (20) is a gradient closure of the unresolved transport due to turbulence motions. Eqs. (12), (14), and (20) are designed to propagate a flame front at the turbulent flame speed SD ¼ Ne S l , according to Eq. (9). The first RHS term in Eq. (21) is estimated according to the procedure described in Appendix B.2:
q0 ðzÞSl ðzÞ
@Y c ¼ @xn
Z 0
1
q0 ðz0 ÞSl ðz0 Þ
@Y c 0 jz ¼ z0 Pðz0 Þdz @xn
ð22Þ
e c =Y fc eq and D ¼ 20dst for two values of mixture fraction unmixedness factor (dst is the flame thickness determined as the Fig. 2. Correction factor aY c as a function of ~z and Y r r full width at half maximum of the reaction rate in the stoichiometric conditions).
Fig. 3. 1-D simulations of stratified flame with F-TACLES for ~z ¼ zst . (a) Temperature profiles. (b) Flame speed; j: effective propagation speed, computed from Eq. (24); –: theoretical propagation speed given by Eq. (10).
2708
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
The second term is computed as follows:
q0 ðzÞSl ðzÞ
fc @ Y fc Z 1 @Y ¼ q ðz0 ÞSl ðz0 ÞPðz0 Þdz0 @xn @xn 0 0
5. LES of a 3-D turbulent stratified swirled flame
ð23Þ
fc =xn is estimated from Eq. (17). In practice, the subgrid with @ Y scale convection term of the progress variable is also stored in a fc ; ~z; Sz ; D. 4-D look-up database: XY c ¼ XY c ½ Y
4. Construction and validation of the filtered look-up table 4.1. Chemical database generation A chemical look-up table for methane/air combustion, representative of the conditions operated in the MOLECULES experimental setup [23–25], described in Section 5.1, is constructed. For that purpose, a collection of 1-D laminar premixed flames is computed using the GRI 3.0 mechanism [33] and the PREMIX [34] solver for equivalence ratio within the flammability limits. The unity Lewis assumption is retained for this study. The progress variable Yc is defined as Y c ¼ Y CO þ Y CO2 as in [28]. The mixture fraction z is defined from nitrogen assumed as inert as the NOx chemistry is not taken into account (Y z ¼ Y N2 ). This coordinate is normalized using Eq. (2). Y fz and Y ox denote values of z Y N2 in fuel stream (here pure methane at 388 K, Y fz ¼ 0) and oxidizer stream (here air at 623 K, Y ox z ¼ 0:77), respectively. From this flamelet library, the procedure proposed in Appendix B is followed to compute the filtered thermochemical quantities. The grid resolution of the filtered look-up table is 200 150 10 in e c , ~z, Sz) subspace. the ( Y
5.1. MOLECULES configuration The combustion chamber MOLECULES, designed to be representative of an actual aeronautical combustor, has been experimentally investigated by Janus et al. [23–25] at TU Darmstadt. The geometry shown in Fig. 4 consists of a plenum, a swirl-injector, and a combustion chamber. Preheated air is introduced in the swirler via the plenum. Fuel is directly injected through the center of the swirler. The combustion chamber pressure is 2 bar at the operating conditions. Air is injected at 623 K at a flow rate of 30 g/s. The global equivalence ratio is 0.8. The fuel used for the experiment is natural gas identified to pure methane in numerical simulations. The flow field has been characterized by Laser Doppler Anemometry for both non-reactive and reactive cases. The mixing has been
4.2. 1-D flame computations As explained in Section 2.2, the present modeling approach describes a flame front evolving in a stratified and turbulent environment at a velocity SD ¼ N Sel , where N is the subgrid flame wrinkling factor. As iso-c surfaces are assumed perpendicular to iso-z surfaces at the subgrid scale level, the propagation speed Sel of the filtered flame front is given by Eq. (10). To validate the model implementation, a series of planar flames evolving in a stratified mixture characterized by a constant subgrid scale variance are computed. The database described in Section 4.1 is used st with a filter size equal to D ¼ 20dst r (dr is the full width at half maximum of the progress variable reaction rate of the stoichiometric flame). A premixed flame under stoichiometric conditions without mixture fraction fluctuations at subgrid scale is computed as a reference. Then, 1-D stratified flames are computed with ~z ¼ zst and 0 < Sz < 1. Temperature profiles are plotted in Fig. 3a and the corresponding flame burning velocities are shown in Fig. 3b. Increasing the subgrid scale mixture fraction variance decreases the burnt gas temperature and increases the flame thickness as shown in Fig. 3a. The predicted flame propagation speed is computed after post-processing 1-D solution profiles data
1 Sel ¼
q0
Z
1
e_ ðxÞdx qx c
ð24Þ
0
and is compared with the flame speed given by Eq. (10). Fig. 3b shows that the effective flame propagation matches the one predicted by the model theory. The model implementation procedure is therefore validated. These tests also demonstrate that the filtered look-up table resolution is sufficient to capture the correct flame propagation speed.
Fig. 4. Schematic views of the MOLECULES burner.
Table 1 Computational parameters used for the simulations. MOLECULES combustion chamber Pressure Air temperature Air mass flow rate Global equivalence ratio Fuel temperature Fuel mass flow rate Reynolds number for air Laminar flame thicknesses Methane injector Swirler exit Chamber dimension
2 bar 623 K 30 g/s 0.8 388 K 1.4 g/s 46,000 5104 m < dc <1.5103 m 4 mm diameter (40 nodes in the radial direction) 28 mm diameter (112 nodes in the radial direction) 60 60 80 mm3
2709
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
measured by Planar Laser Induced Fluorescence (PLIF) technique using acetone as a fuel tracer. In addition, the flame front has been identified with PLIF of the OH radical. As fuel and oxidizer are not premixed before entering the combustion chamber, the modeling approach retained has to account for the mixture fraction stratification. The numerical setup is described in the following section. 5.2. Numerical setup The F-TACLES approach is implemented in the compressible LES code AVBP [35] using the Tabulated Thermochemistry for Compressible flow (TTC) formalism [36]. The filtered balance equation of the progress variable implemented in AVBP reads as follows: Fig. 5. Line positions where numerical simulations are compared with experimental data.
fc Y @q fc Þ fc ; ~z; Sz ; Dq D0 r Y fc Þ ¼ r ðNaY ½ Y ~Y u þ r ðq c 0 @t e_ ½ Y fc ; ~z; Sz ; D þ Nq fc ; ~z; Sz ; D ð25Þ x NXY c ½ Y Yc where the wrinkling factor N is estimated by the analytical model proposed by Charlette et al. [37]:
140
"
120
N¼
Axial velocity
100 80 60 40 20 0 -20 0
20
x [mm]
40
60
Fig. 6. Mean axial velocity component along the centerline for the non-reactive case. Symbols: measurements. Lines: simulation with F-TACLES.
D u0 1 þ min 0 ; C 0D dl Sl
#!0:5 ð26Þ
In this expression, d0l denotes the thermal flame thickness, u0D the subgrid-scale turbulent velocity, and C the spectral efficiency function (see [37]). In order to identify an adequate computational grid, LES of the cold flow has been initially performed on a coarse mesh. Then, the grid resolution has been increased until the convergence of the mean flow field. An unstructured 40 million cells mesh has then been identified and retained for the reactive flow simulation. The fuel nozzle is meshed with 40 cells along the diameter direction. st The filter size D ¼ 20dst r , where dr is the full width at half maximum of the reaction rate of the stoichiometric 1-D premixed flamelet, is retained to generate the filtered look-up table. It corresponds to a filter width 8 times larger that the laminar thermal layer. This choice ensures a computation of the filtered flame front over 5 nodes, which is the minimal resolution required to
Non reactive case 25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-20
-20
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
y [mm]
25
y [mm]
25
0
100
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
Fig. 7. Mean (left) and RMS (right) of axial velocity for the non- reactive case. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
2710
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Non reactive case x = 1mm
x = 5mm
x = 10mm
x = 15mm
x = 20mm
x = 5mm
x = 1mm
x = 10mm
x = 15mm
x = 20mm
25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-20
-20
-25
0 0.5
1
-25
0 0.5
1
-25
0
0.5
-25
0
0.5
-25
y [mm]
25
y [mm]
25
0
0.5
-25
0
0.3
-25
0
0.3
-25
0
0.3
-25
0
0.3
-25
0
0.3
Fig. 8. Mean (Left) and RMS (Right) of CH4 mass fraction for the non-reactive case. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
avoid numerical artifacts. More details and justifications may be found in Ref. [11]. The third-order finite element TTGC scheme [38] is used. Boundary conditions are prescribed using Navier–Stokes Characteristic Boundary Conditions [39]. The WALE turbulence model [40] closes the unresolved shear stress tensor. Boundary parameters are summarized in the Table 1 for both reactive and non-reactive simulations. For velocity field and fuel mass fraction, numerical simulations are compared with experimental data along the centerline and in the plane f = 0 at the locations indicated in Fig. 5.
5.3. Non-reactive case The non-reactive case is first computed. Mean and resolved Root Mean Square (RMS) quantities are computed by time averag140 120
5.4. Reactive case
100
Axial velocity
ing LES solutions over a physical time corresponding to 3 flow through times based on the air inlet velocity (14 ms). This observation time ensures sufficient convergence of the first-order statistical moments on the planes x = 1 mm to x = 20 mm, experimentally investigated by Janus et al. [23–25]. Mean axial velocity component is first compared with the experimental data along the centerline in Fig. 6. The fair agreement between measurements and simulation shows that the recirculation zone is well captured by the LES. Mean radial profiles of axial velocity and CH4 mass fraction are then plotted on Fig. 7 (left) and Fig. 8 (left), respectively. A correct agreement is observed between experimental and numerical profiles, showing that the mixing and the flow field are correctly predicted. RMS of axial velocity and CH4 mass fraction are plotted on Fig. 7 (right) and Fig. 8 (right), respectively. The general trends are well reproduced but an overestimation of the fluctuation amplitudes of CH4 by the simulations is observed at x = 5 mm and x = 10 mm (the relative differences between measured and computed data are about 30 %).
80 60 40 20 0 -20 0
20
x [mm]
40
60
Fig. 9. Mean axial velocity component along the centerline for the reactive case. Symbols: measurements. Lines: simulation with F-TACLES.
For the reactive case, mean and resolved root mean square (RMS) quantities are computed by time averaging LES solutions over a physical time corresponding to 8 flow through times (38 ms) based on the air inlet velocity. The performances of the F-TACLES model are assessed by comparing the predicted velocity and fuel mass fraction to the experimental data. Mean axial velocity component is compared with the experimental data along the centerline in Fig. 9. As for the cold case, the LES results match very well the velocity measurements. Mean axial and radial velocities are plotted in Fig. 10 (top) and Fig. 11 (top), respectively. A fair agreement is observed between experimental and numerical profiles, showing that the flow field is very well predicted in the reactive case. The corresponding RMS are plotted in Fig. 10 (bottom) and Fig. 11 (bottom). As these resolved RMS do not include the subgrid scale RMS’s conclusions regarding the model performance
2711
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
y [mm]
Reactive case 25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
y [mm]
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
0
100
50
Fig. 10. Mean (top) and RMS (bottom) of axial velocity in m s1. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
in terms of flame turbulence interactions are not straightforward. However, it is observed that LES RMS remains lower than measured RMS, as expected from theory, and that the evolution of the predicted RMS is very similar to the experimental results. Mean and RMS CH4 mass fraction are plotted in Fig. 12. The comparison shows that the mixture field is predicted in the reactive case excepted at x = 5 mm and x = 10 mm where deviations of about 30% are observed. A comparison of the predicted mean OH mass fraction versus the experimental PLIF signal is shown in Fig. 13. Both the lift-off height and the flame envelope are captured by the simulation.
The present modeling strategy relies on the assumption that the stratified flame is composed of a collection of flamelets that have a premixed-like structure. To investigate this hypothesis, scatter plots of filtered temperature versus the filtered mixture fraction, ~z, are shown in Fig. 14a where each computational node is represented by one dot. Statistics shown here are estimated from a single snapshot but the same conclusions are drawn when other time instants are analyzed. It shows that all the possible phase subspace is covered but no preferential flame structure is evidenced. Fig. 14b e ) colored by computational shows the reduced phase subspace (e z; T nodes population. This post-processing highlights two major
2712
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Reactive case 20
20
20
15
15
15
10
10
10
20
25 20
15
15
y [mm]
5
5
0
0
0
-5
-5
-5
-10
-10
-10
5
5
0
0 -5
-5
-10
-15
-20
-50 0
50
-15
-20
-50 0
50
20
20
15
15
15
10
10
10
-15
-20
-20
20
-50 0
-25
50
-50
0
50
-20
-50
0
50
20
25 20
15
15
y [mm]
5
5
0
0
0
-5
-5
-5
-10
-10
-10
5
5
0
0 -5
-5
-10
-20
-15
0
30
-20
-15
0
30
-20
20
20
15
15
15
10
10
10
5
5
5
0
0
0
-5
-5
-5
-10
-10
-10
-15
-15
-15
-20
-20
-20
-25
-50
0
50
-25
-50
0
50
-25
-50
25
25
25
20
20
20
15
15
15
10
10
10
5
5
5
0
0
0
-5
-5
-5
-10
-10
-10
-15
-15
-15
-20
-20
-20
0
50
-10 -15
-15
20
10 10
5
25
-10 -15
-15
25
10 10
5
25
-15
-20
0
-25
30
Fig. 11. Mean (top) and RMS (bottom) of radial velocity in m s locations indicated in Fig. 5.
0
30
-20
0
30
-25
0
30
-25
0
30
-25
0
30
1
. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at
trajectories: a line (1) and a narrow vertical strip (2). The line (1) shows that fuel and oxidizer are mixed together prior the chemical reactions. Indeed the transition between fresh and burnt gases occurs in the region (2) for mixture fractions that belong to the flammability limits. The same conclusion is observed in Fig. 15 where the joint filtered density function of temperature and mixture fraction is plotted. Statistical events are clustered in the vicinity of z = 0.05 plane. It indicates that the transition between fresh and burned gases follows premixed flames trajectories. Fig. 16 shows a scatter plot of the subgrid scale flame wrinkling factor N encountered in the flame front reaction zone. As N ranges
from 1 to 3, the impact of the subgrid turbulence on the flame propagation is important and therefore the flame front displacement speed is sensitive to the subgrid scale wrinkling factor model given in Eq. (26). A second reactive simulation has been performed without accounting for subgrid scale mixture fraction heterogeneities, that is, by modeling P(z) by a Dirac d function centered on ~z ðPðzÞ ¼ dðz ~zÞÞ. Figure 17 shows two 2-D snapshots of the instantaneous temperature iso-color field with P(z) modeled by a b-function (left) or by a Dirac function (right). The temperature, computed from the filtered chemical look-up table, is strongly
2713
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Reactive case
y [mm]
y [mm]
x = 1mm
x = 5mm
x = 10mm x = 15mm x = 20mm x = 30mm x = 40mm
20
20
20
20
20
20
20
10
10
10
10
10
10
10
0
0
0
0
0
0
0
-10
-10
-10
-10
-10
-10
-10
-20
-20
-20
-20
-20
-20
-20
0 0.5 1
0 0.5 1
x = 1mm
x = 5mm
0
0.6
0
0.4
0
0.3
0 0.03
0 0.1
x = 10mm x = 15mm x = 20mm x = 30mm x = 40mm
20
20
20
20
20
20
20
10
10
10
10
10
10
10
0
0
0
0
0
0
0
-10
-10
-10
-10
-10
-10
-10
-20
-20
-20
-20
-20
-20
-20
0 0.2
0 0.2
0
0.3
0
0.3
0
0.3
0
0.1
0 0.04
Fig. 12. Mean (top) and RMS (bottom) of CH4 mass fraction. Symbols: measurements. Lines: simulation with F-TACLES.
grid scale mixture fraction effects, the front displacement speed is too high. The flame is then stabilized at the burner lips and is not lifted anymore. Then, because of a misprediction of the flame anchoring mechanism, accurate statistics on temperature and velocity cannot be obtained.
5.5. Comparison with other modeling strategies applied in the MOLECULES combustor
Fig. 13. Iso-color of simulated (left) and measured (right) OH concentrations in the combustor. OH is made dimensionless by its maximal value.
dependent of the mixture fraction variance. As expected, peaks of temperature are smoothed when the impact of subgrid mixture fraction variance is considered. Also a strong influence of P(z) on the flame dynamics is observed. Indeed, without considering sub-
Numerical investigations of the MOLECULES combustion chamber have been performed in the past by other groups following different modeling strategies. U-RANS (Unsteady Reynolds Average Navier-Stokes) simulations have been conducted by Schneider et al. [41] with two turbulent combustion models: the G-equation [12] and the Bray-Moss-Libby approaches [42]. Results do not show specific advantage of one or other premixed combustion model in a RANS context and fair comparisons versus experiments has been achieved. However, the comparison between measured
2714
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Fig. 14. Left: Scatter plot of the filtered temperature versus the filtered mixture fraction in the MOLECULES case. The red line corresponds to the global equivalence ratio. Right: Combustion trajectories in the MOLECULES burner. Red color corresponds to a number of grid points larger than 10,000 and the blue lower than 50.
Fig. 15. Joint filtered density function of the mixture fraction and temperature.
e c =Y e eq Fig. 16. Scatter plot of the subgrid flame wrinkling for 0:1 < Y c < 0:9.
Fig. 17. 2-D snapshots of the temperature field. The Filtered Density Function P(z) is either modeled by a b function (left) or by a Dirac d function (right). The flame is then stabilized at the burner lips when the subgrid scale mixture fraction heterogeneities are not accounted for (right). The solid black line shows the stoichiometric line.
and simulated 2-D field of OH radical shows that the present LES results obtained with F-TACLES significantly improves the prediction of the flame front position. A possible explanation of the difference is the representation of the turbulence. Indeed, analyzing velocity data [41], the authors highlight the difficulty to accurately describe the turbulence with U-RANS. As shown both in the cold flow LES of the combustor performed by Sadiki et al. [43] and in the present non-reactive case simulation, a better representation of the flow dynamic is obtained with LES. LES of the reactive flow in the MOLECULES combustor has also been performed by Sadiki et al. [43] using two different turbulent combustion models: the RPV (Reaction Progress Variable) model and the steady diffusion flamelet model [12]. The diffusion flamelet model failed to reproduce the flame lift-off because the chemistry tabulation is not adapted to the stratified flame structure. In the RPV approach, the chemistry is tabulated from 1-D premixed flamelets as in FPI (Flame Prolongation of ILDM) [27] or FGM (Flamelet Generated Manifold) [29] approaches. However, unlike to F-TACLES, the filtered density function of the progress variable is simply
2715
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
modeled by a Dirac function. Subgrid scale flame front contributions are not accounted for and may cause errors in the flame front propagation. It may explain why predictions are closer to the experimental data with F-TACLES than with RPV.
6. Conclusions Issues relative to LES of stratified combustion are discussed in this article. An extension of the F-TACLES model, initially developed for premixed combustion, is proposed to tackle turbulent stratified flames. Detailed chemistry effects are accounted for in an LES context by storing in a look-up table a set of 1-D laminar filtered premixed flames. A mathematical formalism, including the impact of subgrid scale mixture fraction heterogeneities and subgrid scale flame wrinkling on the flame front propagation, is proposed to couple this filtered thermo-chemical table to the LES flow governing equations. The model is implemented in the compressible AVBP solver [35] to perform the 3-D simulations of a stratified swirled turbulent flame representative of an aeronautical combustor. Both reactive and non-reactive simulations are compared with experimental data. In general, a very good agreement between numerical simulations and measurements is observed. A comparison of the simulations performed with other techniques published in the literature demonstrates the performances of the F-TACLES approach.
e_ ¼ q x Yc
Appendix A. Modeling assumptions in F-TACLES formalism
jrcj
h/i ¼
Z
1
R
1
The filtered reaction rate of the progress variable is given by:
e_ ¼ q x Yc
Z
qðx0 Þx_ Y c ðx0 ÞGD ðx x0 Þdx0
ðA:1Þ
X
_ Y c ðx0 Þ is the progress variable reaction rate at location x0 where x and GD is the LES filter. This reaction rate may be rewritten as follows:
e_ ¼ q x Yc
Z
1
Z
qðx0 Þx_ Y c ðx0 Þ
0
jrcjdðcðx0 Þ c0 ÞGD ðx x0 Þdx0 dc jrcj 0 X Z 1 qx_ Y c ¼ Rðc0 Þdc0 ðA:2Þ jrcj c0 0
where d(c) is the Dirac delta-function and Rðc0 Þ ¼ jrcjdðc c0 Þ the density of the c = c0 -isosurface, that is, the available c = c0 -surface per unit filter volume. h/ic0 denotes the average of the quantity / over the c = c0 -surface in the filtering volume, defined as [44]:
R /jrcjdðcðx0 Þ c0 ÞGD ðx x0 Þdx0 RX / ¼ h ic0 0 0 0 0 X jrcjdðcðx Þ c ÞGD ðx x Þdx
ðA:3Þ
Following [44], a generalized surface density, R, is then introduced:
R¼
Z 0
1
Rðc0 Þdc0 ¼
Z
1
0
jrcjdðc c0 Þdc ¼ jrcj
0
and the filtered reaction rate expressed as follows:
ðA:4Þ
Rðc0 Þdc0 ¼
qx_ Y c
c0
jrcj
R¼
qx_ Y c jrcj
Njrcj ðA:5Þ
0
R /jrcjGD ðx x0 Þdx0 0 h/ic0 Rðc0 Þdc ¼ RX 0 0 X jrcjGD ðx x Þdx
ðA:6Þ
_ Y c =jrcji measures the mean reaction rate per unit of In Eq. (A.5), hqx flame area. The above relations are exact, according to the corresponding definitions. F-TACLES model assumptions are now discussed for two cases successively: (i) perfectly premixed combustion (i.e., the original F-TACLES model [1]); (ii) stratified combustion, an extension proposed in this paper. A.2. Perfectly premixed combustion By construction, the F-TACLES model is a flamelet model: the flame surface is identified to a collection of laminar flame elements. Surface densities R(c0 ) are assumed to be independent on c0 and the progress variable reaction rate per unit of flame surface, _ Y c =jrcji, is computed from one-dimensional laminar flame, hqx using Eq. (A.6):
q x0n x_ Y c x0n GD ðxn x0n Þdx0n
1
R þ1
jrcj
R þ1
0
j@c =@xn jGD ðxn x0n Þdxn Z þ1 1 q x0n x_ Y c x0n GD ðxn x0n Þdx0n j@c =@xn j 1 1
ðA:7Þ
where the ⁄ superscript denotes quantities issued from a 1-D laminar unstretched premixed flame computed in the xn direction. According to Eq. (A.5), the filtered progress variable reaction rate in turbulent flames is modeled as follows:
e_ Y q x c
Z
þ1
1
A.1. Definitions and general relations
where N ¼ jrcj=jrcj is the surface wrinkling factor (i.e., the ratio of the generalized surface divided by the resolved surface) and h/i the generalized surface average defined by:
qx_ Y c
We thank Prof. Andreas Dreizler and Dr. Simone Schoenfelder for providing us the burner geometry and the experimental data. This work was supported by the ANR-07-CIS7–008-04 Grant of the French Ministry of Research and accessed the high-performance computing resources of IDRIS and CINES under the allocation 2009-i2009020164 made by GENCI (Grand Equipement National de Calcul Intensif).
qx_ Y c
0
Acknowledgments
Z 1
q x0n x_ Y c x0n GD ðxn x0n Þdx0n N
ðA:8Þ
assuming that the resolved flame thickness are similar in 1-D laminar flames and in turbulent flames (i.e., j@c =@xn j jrcj). Eq. (A.8) is the F-TACLES model expression for the filtered progress variable reaction rate [1]. The quantity inside brackets is computed from 1-D laminar flames and tabulated as function of the filtered progress e c and the filter size D. variable Y A.3. Stratified combustion In stratified combustion regime, the reaction rate of the filtered progress variable is influenced by the mixture heterogeneities at the subgrid scale. The challenge is now to estimate the reaction _ Y c jrcji rate of the filtered progress variable per unit surface, hqx (Eq. (A.5)), in this situation. By definition:
qx_ Y c jrcj
R
¼
X
qðx0 Þx_ Y c ðx0 ÞGD ðx x0 Þdx0 R
R1 R ¼
0
R1
¼ R 10 0
X
jrcjGD ðx x0 Þdx0 0
0 _ 0 0 0 0 0 X qðx ÞxY c ðx Þdðzðx Þ z ÞGD ðx x Þdx dz R1 R 0 0 0 0 0 0 X jrcjdðzðx Þ z ÞGD ðx x Þdx dz
_ Y c jz ¼ z0 ÞPðz0 Þdz0 ðqx 0
ðjrcjjz ¼ z0 Þ Pðz0 Þ dz
ðA:9Þ
where z is the mixture fraction going from z = 0 for pure oxidizer to z = 1 for pure fuel. P(z0 ) is the subgrid scale filtered density function of the mixture fraction, defined by [45]:
2716
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Pðz0 Þ ¼
Z
dðzðx0 Þ z0 ÞGD ðx x0 Þ dx0
ðA:10Þ
X
The factor aY c correcting the filtered laminar diffusion term (Eq. (15)):
ð/jz ¼ z0 Þ is the filtered value of / conditioned on the iso-surface z = z0 :
R
ð/jz ¼
z0 Þ
¼
Þdðzðx0 Þ z0 ÞGD ðx x0 Þdx0 0 Þ z0 ÞG ðx x0 Þdx0 dðzðx D X
/ðx XR
0
_ Y c jz ¼ ðqx
Z
All these terms can be expressed in the form:
¼ Z ðjrcjjz ¼ z0 Þ ¼
þ1
1 þ1
1
ðx0n ; z0 Þ
0 0 Y c ðxn ; z ÞGD ðxn
q x_
@c 0 0 0 0
@x ðxn ; z Þ GD xn xn dxn n
x0n Þ
0 dxn
ðA:12Þ ðA:13Þ
Eq. (A.9) becomes R 1 R þ1 0 0 0 0 q ðx ; z Þx_ Y c ðxn ; z ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0 qx_ Y c
hR n
i 0 R1
1 þ1 @c 0 0 jrcj
ðx0 ; z0 Þ G ðx x0 Þ dx Pðz0 Þ dz 0
1 jrcj
1
Z
1
Z
0
@xn
þ1
1
D
n
n
n
q ðx0n ; z0 Þx_ Y c ðx0n ; z0 ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0
Assuming that jrcj extracted from one-dimensional laminar flame computations combined with presumed filtered density functions P(z0 ) is similar to jrcj in turbulent flame simulations (i.e., jrcj jrcj) allow to express the filtered reaction rate of the progress variable (Eq. (A.5)) as follows:
e_ N q x Yc
Z
0
þ1
1
Z
1
0
Z
q ðx0n ; z0 Þx_ Y c ðx0n ; z0 ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0
ðA:15Þ The quantity inside brackets is computed from 1-D laminar flames while the filtered density function P(z0 ) of the mixture fraction is modeled as a b-function and depends on the mean, ~z, and the 002 =~ unmixedness factor, Sz ¼ zf zð1 ~zÞ, of the mixture fraction z. The quantity inside parenthesis is then tabulated as function of e c , the filter size D, ~z and Sz, while the filtered progress variable Y the flame surface wrinkling factor, N, quantifies the flame/turbulence interactions.
B.1. Introduction This Appendix describes the practical construction of the filtered database used in the F-TACLES model. According to the model derivation, three terms are estimated from one-dimensional flame computations and then stored in the database: The contribution to the filtered reaction rate (Eq. (A.15)): 1
0
hqxY c jz ¼ z0 iPðz0 Þdz
0
ðB:1Þ
The contribution of the thermal expansion to the unresolved transport (Eqs. (21)–(23)):
XY c ¼
Z
A collection of 1-D premixed laminar flame is first computed with a detailed chemical scheme for equivalence ratios covering the mixture flammability limits. Thermo-chemical data u, such as mass fractions and reaction rates and solutions of the 1-D flam elets, are mapped in the 2-D x0n ; z0 subspace where x0n denotes the direction normal to the flamelet and z0 is the mixture fraction:
u ¼ u x0n ; z0
0
q
0
¼ z0 iPðz0 Þdz
@xn fc Z @Y @xn
ðB:6Þ
Outside the flammability range, chemical reaction rates are set to zero while species mass fractions are linearly interpolated between the point of the last tabulated flame on the lean (respectively rich) side and the air (respectively fuel) stream. For a given value of the filter size D, conditional filtered quantities are then computed for constant mixture fraction value z = z0 according to:
hujz ¼ z0 iðxn ; z0 ; DÞ ¼
Z
þ1
1
uðx0n ; z ¼ z0 ÞGD xn x0n dx0n
ðB:7Þ
where GD denotes the one-dimensional filter function, chosen here as Gaussian. The practical integration of the RHS in Eq. (B.7) is detailed in Appendix B.3. According to Eq. (1), conditional values of the progress variable hYcjz = z0 i and of the normalized progress variable hcjz = z0 i are related through the relation:
hY c jz ¼ z0 i Y fc ðz0 Þ Y bc ðz0 Þ Y fc ðz0 Þ
ðB:8Þ
which is used to map the conditional variables huj z = z0 i(xn, z0 , D) into the subspace (hcjz = z0 i,z0 , D) where hcjz = z0 i replaces the space coordinate xn. P(z0 ) is now presumed by a b-function P b ðz0 ; ~z; Sz Þ parametrized by the filtered mixture fraction ~z and the corresponding unmixed002 =~ ness factor Sz ¼ zf zð1 ~zÞ. Eq. (B.5) is then recast as follows: 002 ; DÞ ¼ ðc; ~z; zf u
Z
1
0 hujz ¼ z0 iðhcjz ¼ z0 i; z0 ; DÞPb ðz0 ; ~z; Sz Þ dz
ðB:9Þ
0
where
Z
1
0 hcjz ¼ z0 iPb ðz0 ; ~z; Sz Þ dz
ðB:10Þ
0
@Y c
0 0 0 ðz ÞSl ðz Þh
ðB:5Þ
and are easily extracted from one-dimensional laminar flames. When the mixture fraction is constant: P(z0 ) = d(z0 z0), where d is the Dirac delta-function. Then, only one flame is considered for a given mixture fraction value z0, as in the original F-TACLES model presented in Appendix A.2. For variable mixture fractions, several 1-D flames are combined through the filtered density function P(z0 ) as described in the following.
c ¼ 1
hujz¼z0 iðxn Þ
hcjz ¼ z0 i ¼
Z
u x0n ; z0 GD ðxn x0n Þ dx0n Pðz0 Þ dz0
1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
Appendix B. Construction of the filtered database
q x_ Y c ¼
þ1
B.2. General procedure
ðA:14Þ
1
¼ u
n
Z
ðB:3Þ
q0 D0 @@xYecn
ðA:11Þ
_ Y c jz ¼ z0 Þ, as well as The filtered progress variable reaction rate, ðqx the gradient of the progress variable, ðjrcjjz ¼ z0 Þ, both conditioned on the z = z0 isosurface, should now be determined (see Eq. (A.9)). F-TACLES assumes that these terms can be estimated from a one-dimensional laminar flame at constant mixture fraction z = z0 . Then,
z0 Þ
c q D @Y @xn
aY c ¼
jz
0
1
q0 ðz0 ÞSl ðz0 ÞPðz0 Þdz0
ðB:2Þ
The value of the conditional normalized progress variable hcjz = z0 i should now be specified for each mixture fraction value z0 to integrate Eqs. (B.9) and (B.10). Assuming that iso-c and iso-z surfaces are orthogonal at the subgrid scale level, the conditional filtered
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
GD ðxÞ ¼
6 p D2
1=2
6x2 exp 2 D
2717
ðB:18Þ
References
Fig. 18. Scheme for Gaussian convolution.
normalized progress variables hcjz = z0 i is kept constant during the integration. Note that this procedure corresponds to the following assumption: P(c, z) = P(c)P(z), where P(c, z) is the joint progress variable/mixture fraction FDF. Following the same methodology, Favre filtered thermochemical quantities are computed from the definiu ~ ¼ qu. In particular, the Favre filtered progress variable is tion q estimated by:
Z
fc ¼ 1 Y q
1
0 hqY c jz ¼ z0 iPb ðz0 ; ~z; Sz Þ dz
ðB:11Þ
0
e c , ~z, zf 002 , and u ~ are then mapped in the 4-D ( Y Filtered quantities u D) subspace to terminate the filtered chemical table generation. A constant value of D is chosen in the present article and balance 02 are solved to directly couple the filtered e c , ~z and zf equations for Y look-up table to the LES flow governing equations. B.3. Numerical integration of Eqs. (B.7) and (B.9) The F-TACLES model predictions are very sensitive to the quality of the numerical computations of (B.7) and (B.9). The integration of the b-FDF (filtered density function) P(z) in Eq. (B.9) is performed following the methodology proposed in [46]. The progress variable c(x) is a monotonic function that evolves between 0 (fresh gases) and 1 (burnt gases) as shown in Fig. 18. x0 and x1 are introduced such as c(x0) = 0 and c(x1) = 1. Thus, c is a bijection defined on [x0, x1] to [0, 1]. In order to compute the integral in Eq. (B.7), the flame is divided in three zones as shown in Fig. 18:
hujz ¼ z0 i ¼ I0 þ I þ I1
ðB:12Þ
where I0, I and I1 are as follows:
Z x0 I0 ¼ uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0 1 Z x1 uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0 I¼
ðB:13Þ ðB:14Þ
x0
I1 ¼
Z
þ1
uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0
ðB:15Þ
x1
I is numerically computed with a second-order integration scheme. The range ]1, x0] (resp. [x1, +1[) corresponds to fresh gases (resp. burnt gases); therefore, u(z = z0 , x0 ) is constant and equal to u(z = z0 , c = 0) (resp. u(z = z0 , c = 1)). I0 and I1 are then computed directly from the analytical solution:
!# pffiffiffi 6 1 erf I0 ¼ ðx x0 Þ 2 D " !# pffiffiffi 6 uðz ¼ z0 ; c ¼ 1Þ ðx x1 Þ I1 ¼ 1 erf 2 D
uðz ¼ z0 ; c ¼ 0Þ
"
the 1-D filter GD being defined by a Gaussian function:
ðB:16Þ ðB:17Þ
[1] B. Fiorina, R. Vicquelin, P. Auzillon, N. Darabiha, O. Gicquel, D. Veynante, Combust. Flame 157 (2010) 465–475. [2] A. Unal, Y. Hu, M. Chang, M.T. Odman, A. Russell, Atmos. Environ. 39 (3) (2005) 5787–5798. [3] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, R.T. Edwards, Inc., 2005. [4] H. Pitsch, Combust. Flame 143 (4) (2005) 587–598. [5] J. Galpin, A. Naudin, L. Vervisch, C. Angelberger, O. Colin, P. Domingo, Combust. Flame 155 (2008) 247–266. [6] C. Duwig, Combust. Theor. Modell. 13 (2) (2009) 251–268. [7] S. Menon, W. Jou, Combust. Sci. Technol. 75 (1991) 53–72. [8] O. Colin, F. Ducros, D. Veynante, T. Poinsot, Phys. Fluids 12 (7) (2000) 1843– 1863. [9] M. Boileau, G. Staffelbach, B. Cuenot, T. Poinsot, C. Berat, Combust. Flame 154 (1–2) (2008) 2–22. [10] J. Janicka, A. Sadiki, Proc. Combust. Inst. 30 (1) (2005) 537–547. [11] P. Auzillon, B. Fiorina, R. Vicquelin, N. Darabiha, O. Gicquel, D. Veynante, Proc. Combust. Inst. 33 (1) (2011) 1331–1338. [12] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [13] V.-K. Chakravarthy, S. Menon, Flow Turbul. Combust. 65 (2000) 133–161. [14] E. Knudsen, H. Pitsch, Combust. Flame 154 (4) (2008) 740–760. [15] V. Moureau, B. Fiorina, H. Pitsch, Combust. Flame 156 (4) (2009) 801–812. [16] J. Piana, D. Veynante, S. Candel, T. Poinsot, Direct numerical simulation analysis of the g-equation in premixed combustion, in: J. Chollet, P. Voke, L. Kleiser (Eds.), Direct and Large Eddy Simulation II, Kluwer Academic Publishers, 1997, pp. 321–330. [17] V. Moureau, P. Minot, H. Pitsch, C. Berat, J. Comput. Phys. 221 (2) (2007) 600– 614. [18] H. Pitsch, Annu. Rev. Fluid Mech. 38 (2006) 453–482. [19] B. Li, E. Baudoin, R. Yu, Z. Sun, Z. Li, X. Bai, M. Alden, M. Mansour, Proc. Combust. Inst. 32 (2) (2009) 1811–1818. [20] A. Roux, H. Pitsch, Large-eddy simulation of stratified combustion, in: Annual Research Briefs, Center For Turbulence Research, 2010, pp. 275–289. [21] B.A. Sen, S. Menon, Combust. Flame 157 (1) (2010) 62–74. [22] P. Domingo, L. Vervisch, D. Veynante, Combust. Flame 152 (3) (2008) 415–432. [23] B. Janus, A. Dreizler, J. Janicka, ASME Turbo Expo, Vienna (2004) 293–315. [24] B. Janus, A. Dreizler, J. Janicka, Flow Turbul. Combust. (75) (2004) 293–315. [25] B. Janus, A. Dreizler, J. Janicka, Proc. Combust. Inst. 31 (2) (2007) 3091–3098. [26] U. Maas, S. Pope, Combust. Flame 88 (1992) 239–264. [27] O. Gicquel, N. Darabiha, D. Thévenin, Proc. Combust. Inst. (2000) 1901–1908. [28] B. Fiorina, R. Baron, O. Gicquel, D. Thévenin, S. Carpentier, N. Darabiha, Combust. Theor. Modell. 7 (2003) 449–470. [29] J. van Oijen, L.P.H. de Goey, Combust. Theor. Modell. 8 (2004) 141–163. [30] B. Fiorina, O. Gicquel, L. Vervisch, S. Carpentier, N. Darabiha, Combust. Flame 140 (3) (2005) 147–160. [31] P. Nguyen, L. Vervisch, V. Subramanian, P. Domingo, Combust. Flame 157 (1) (2010) 43–61. [32] L. Vervisch, R. Haugel, P. Domingo, M. Rullaud, Flow Turbul. Combust. 5 (4) (2004) 1–36. [33] G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K. Hanson, S. Song, W.C. Gardiner, V.V. Lissianski, Z. Qin,
Contents lists available at SciVerse ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
A Filtered Tabulated Chemistry model for LES of stratified flames Pierre Auzillon, Olivier Gicquel, Nasser Darabiha, Denis Veynante, Benoit Fiorina ⇑ CNRS, UPR 288, Laboratoire d’Energétique Moléculaire et Macroscopique, Combustion (EM2C), Grande Voie des Vignes, 92290 Châtenay-Malabry, France Ecole Centrale Paris, Grande Voie des Vignes, 92290 Châtenay-Malabry, France
a r t i c l e
i n f o
Article history: Available online 23 April 2012 Keywords: Large Eddy Simulation Stratified combustion Tabulated chemistry Turbulent combustion modeling
a b s t r a c t The model called Filtered Tabulated Chemistry for LES (F-TACLES), recently developed for turbulent premixed combustion [1], is extended in the present article to stratified flames where the equivalence ratio is not spatially uniform. The principle is to generate a look-up table where thermo-chemical variables, computed from a set of filtered laminar unstrained premixed flamelets, are stored as a function of the ~ c , the filtered mixture fraction ~c, and the mixture fraction subgrid scale varifiltered progress variable Y 002 for a given value of the LES filter size D. From this database and using a classical subgrid scale ance zf flame surface wrinkling model, a closure model for the filtered progress variable balance equation is then 002 ) are added to the Navier–Stokes ~ c , ~z and zf proposed. As only three balance equations (respectively for Y system, this turbulent combustion model is straightforward to implement and is not expensive in terms of CPU time. After preliminary tests on 1-D planar flames, the proposed strategy is applied to a 3-D simulation of a swirled turbulent stratified flame. The agreements observed between numerical and measured data validate the modeling strategy. Finally, it is observed that the flame position cannot be captured without accounting for the subgrid scale mixture fraction heterogeneities that impact on the flame propagation. This result highlights the importance of the mixture fraction filtered density function on the filtered chemical reaction rate. Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Future aeronautical engines will introduce new combustion chambers based on lean-premixed, multi-injection, and/or prevaporized technologies. These systems promote combustion in lean premixed-like regimes to control the flame temperature and consequently the pollutant emissions such as nitrogen oxides [2]. For safety reasons, fuel and oxidizer are injected separately before entering the combustion chamber; therefore, in practice reactants are not perfectly mixed (non-spatially uniform equivalence ratio) and stratified combustion regimes are observed. The emergence of these complex advanced combustion devices introduces new challenges in terms of modeling [3]. A crucial issue is to capture the correct propagation speed of the flame front on grid coarser than the flame thickness. To overcome this difficulty, efforts have recently been made for fully premixed situations [1,4–8]. The popular thickened flame model (TFLES) initially developed by Colin et al. [8] is robust and can easily be applied to realistic combustion chambers [9,10]. However, a too large flame front thickening may affect the flame response to unsteady flow motions [11]. Another approach is to solve a smooth scalar field where a given iso-surface is identified to the instantaneous flame front posi⇑ Corresponding author at: Ecole Centrale Paris, Grande Voie des Vignes, 92290 Châtenay-Malabry, France. Fax: +33 1 47028035. E-mail address: benoit.fi[email protected] (B. Fiorina).
tion. In such technique, called G-equation model, the inner layer is tracked using a level-set technique. Initially developed in a RANS context [12], the G-equation has been reformulated for LES [4,7,13–15]. These methods are very efficient to capture the flame front propagation but, as pointed out in Refs. [16,17], the coupling between the flame front position and flow equations through heat expansion is very challenging. A recently-proposed alternative technique called Filtered Tabulated Chemistry for LES (F-TACLES) [1] consists in generating a chemical look-up database from filtered 1-D unstrained laminar premixed flames computed with detailed chemical schemes. This approach, straightforward to implement in a Navier–Stokes solver, has shown very promising results on turbulent premixed flames and provides a more accurate prediction of the DNS flame dynamics than the TFLES approach [11]. Unfortunately, in practical situations, the mixing between fuel and oxidizer is never perfect and the flame front propagates in an heterogeneous environment. In an LES context, the equivalence ratio fluctuations at the subgrid scale will affect the propagation of the filtered flame front. The G-equation formalism is well adapted to LES of premixed turbulent combustion because it is designed to accurately describe flame front displacements [7,15,18]. Interesting attempts have been made to perform LES of stratified systems with this approach [19,20], but modeling efforts are still required to apply this strategy in industrial combustors, especially when local exinctions occur.
0010-2180/$ - see front matter Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2012.03.006
2705
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Recent developments of the Linear Eddy Mixing (LEM) model combined with an Artificial Neural Network (ANN) enable LES of both premixed and non premixed systems with a skeletal chemical mechanism [21]. In the context of tabulated chemistry, the FPI-PCM [22] method accounts for mixture fraction subgrid scale fluctuations with the introduction of a presumed b filtered density function for the mixture fraction. However, this approach does not capture the correct propagation speed of the flame front. In particular, the laminar flame speed is not recovered when turbulence vanishes [1]. Finally, effects of subgrid scale stratification on the flame propagations are not predicted by the TFLES model that does not account for equivalence ratio fluctuations in chemical reaction rates closures. The objective of this paper is to draw a modeling route for LES of stratified flames, where the flame front propagates into nonperfectly mixed reactants. In the first section, the modeling issues for LES of stratified flames are presented. A tabulated chemistry formalism adapted to stratified combustion is introduced, and the propagation of a filtered flame front into a stratified environment is discussed. In a second section, the F-TACLES model initially developed for premixed combustion is extended to turbulent stratified flames. Numerical simulations of planar laminar flames that propagate in an heterogeneous equivalence ratio flow field are presented. Finally, 3-D LES of the experimental MOLECULES combustor where a swirled flame is stabilized behind a TURBOMECA injector are performed. Simulation results are compared with experimental data obtained by Janus et al. [23–25]. The mathematical derivation of the F-TACLES formalism is detailed in Appendix A while Appendix B is devoted to the practical implementation of tabulated chemistry. 2. Modeling issues for Large Eddy Simulation of stratified flames 2.1. Chemistry tabulation for stratified combustion A careful description of chemical reactions is required to accurately predict the flame structure and the pollutant formation. As hydrocarbon combustion generally involves hundreds of species, simulations of practical combustion systems including complex chemical schemes are not possible. Various studies have been performed during the last decades to account for detailed chemistry effects for a reasonable computational cost. One of the most popular techniques is the use of tabulated chemistry. Maas and Pope [26] have shown that in practical reactive systems chemical composition mainly evolves on a reduced subspace called manifold. Different strategies based on mathematical analysis such as Intrinsic Low Dimensional Manifold (ILDM) [26], or based on physical considerations like Flame Prolongation of ILDM (FPI) [27,28], or Flamelet Generated Manifold [29], or using artificial neural networks [21] have been suggested to identify these manifolds. For each of these approaches, thermochemical quantities are stored in a multi-dimensional look-up database. These techniques are able to predict regimes such as perfectly premixed or stratified combustion [30]. Recently, Nguyen et al. [31] have extended this formalism to diffusion flames. In this article, the FPI approach [27,28] is retained. The chemical subspaces accessed by a stratified flame in a complex geometry configuration is mapped by a collection of 1-D laminar premixed flames computed for various equivalence ratio within the flammability limits using detailed chemistry. The rich (lean) flammability limit is defined as the maximal (minimal) possible value of the equivalence ratio of a premixed flame. Above (bellow) that point, a freely propagating premixed planar laminar flame will extinguish. Note that this formalism excludes the possible combustion of reactant mixtures below or above flammability limits sustained by hot burnt gases issued from the combustion of flammable mix-
tures. In the present study, laminar thermochemical quantities are stored in a look-up table as a function of two coordinates: the progress variable Yc that evolves monotonically between fresh and burnt gases and a mixture fraction z equal to zero in pure oxidizer and unity in pure fuel. For methane–air combustion, preliminary studies have shown that an efficient definition of this progress variable is obtained from a linear combination of CO and CO2 mass fractions: Y c ¼ Y CO þ Y CO2 [28]. Such a definition of Yc ensures a one-to-one correspondence between species mass fractions and the progress variable. A normalized form of the progress variable c is introduced:
c¼
Y c Y fc
ð1Þ
Y bc Y fc
where Y fc and Y bc are the values of progress variable in fresh and burnt gases, respectively. The mixture fraction z is defined from Yz, an inert species or a linear combination of atomic elements representative of fuel and air mixing:
z¼
Y z Y ox z
ð2Þ
Y fz Y ox z
where Y fz and Y ox z denote values of Yz in fuel stream and oxidizer stream, respectively. 2.2. Propagation of a filtered flame front in a stratified mixture 2.2.1. Premixed flame front (uniform equivalence ratio) A premixed wrinkled reaction layer, propagating in perfectly mixed reactants (spatially uniform equivalence ratio), is first considered and illustrated in Fig. 1a. The solid black line represents the interface between fresh and burnt gases before applying the LES filter. The wide gray area shows the corresponding filtered flame front. As demonstrated in Appendix A, when the flame structure is assumed to be similar to the structure of a planar filtered flame, the progress variable reaction rate reads [1]:
e_ ¼ N q x Yc
Z
þ1
1
q x_ Y c x0n GD ðxn x0n Þdx0n
ð3Þ
where the subgrid scale flame surface wrinkling factor N ¼ jrcj=jrcj measures the ratio of the total flame surface to the resolved flame surface in the filter volume. The ⁄ superscript denotes quantities issued from computations of 1-D laminar unstrained premixed flames computed as function of the direction normal to the flame front, xn, for a constant equivalence ratio (or mixture fraction z). GD is a 1-D filter function of size D. The propagation speed SD of the filtered flame front is obtained by integrating the filtered reaction rate given in Eq. (3) along the normal direction to the reaction layer:
q0 SD ¼
Z
1 Y bc Y fc
þ1
Z
þ1
N
1
1
q x_ Y c x0n GD xn x0n dx0n dxn
ð4Þ
where q0 is the fresh gas density. This equation leads to, assuming that the flame wrinkling factor N is uniform over the flame brush:
SD ¼ NSl
ð5Þ
where Sl is the flame speed of a planar laminar premixed flame front:
q0 Sl ¼
Z
1 Y bc
Y fc
þ1
1
q x_ Y c ðxn Þdxn
ð6Þ
2.2.2. Stratified flame front (non-uniform equivalence ratio) The propagation of a flame front in a stratified environment where the equivalence ratio is not spatially uniform is illustrated in Fig. 1b where the dashed lines represent the iso-mixture fraction
2706
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Fig. 1. Propagation of a filtered flame front.
lines. It has been shown in Refs. [30,31] that the orientation of progress variable isosurfaces with mixture fraction isosurfaces, related to the cross-scalar dissipation rate vzc = qDrc rz, may affect the chemical flame structure. However, to take into account this effect in practical situations is not straightforward and considerably complicates the chemical look-up table generation. In order to obtain an easy-to-implement modeling strategy relevant for practical applications, the orientation of iso-c versus iso-z surfaces is assumed constant and perpendicular: vzc = 0. As indicated in Appendix B, this simplification supposes that the normalized form of the progress variable, c, and the mixture fraction, z, are independent variables at the subgrid scale level (unlike to Yc and z which are not independent, as the equilibrium value of Yc depends on the mixture fraction z [32]). Under these assumptions, Appendix A shows that the filtered reaction rate may be written as follows:
e_ ¼ N q x Yc
Z
1
0
_ Y c jz ¼ z0 iPðz0 Þdz hqx
ð7Þ
0
where P(z0 ) is the filtered density function of the mixture fraction _ c jz ¼ z0 i, the conditional filtered value of the reaction rate and hqx at z = z0 , is modeled by:
_ Y c jz ¼ z0 i ¼ hqx
Z
þ1
1
q x_ Y c
0 x0n ; z0 GD ðxn x0n Þdxn
ð8Þ
ð9Þ
where
1 e Sl ¼ q 0
Z
1 0
q0 ðz0 ÞSl ðz0 ÞPðz0 Þdz0
where q is the density, u the velocity vector, and D the diffusivity. x_ Y c is the progress variable reaction rate expressed in s1. Closure models have been proposed in [1] for the subgrid scale transport ~c q fc Þ, the filtered laminar diffusion term ~Y uY u terms r ðq e_ Y . This methodology is here qDrY c and the filtered source term x c extended to account for mixture fraction heterogeneities at the subgrid scale level. e_ Y 3.1. Filtered chemical reaction source term x c e_ Y reads as As shown in Appendix A the expression for x c follows:
1
e_ ¼ N x Yc q
Z
1
ð10Þ
where q0(z), the fresh gas density, is function of the mixture fracR 0 ¼ 01 q0 ðz0 ÞPðz0 Þdz0 . tion and q0 is then expressed by q 3. F-TACLES for stratified flames The balance equation for the filtered normalized progress variable ~c contains several cross-terms introduced by the dependency of Y bc and Y fc on the mixture fraction and is difficult to close in practice. Solving a transport equation for the filtered non-normalized e c is then preferred: progress variable Y
0
_ Y c jz ¼ z0 iPðz0 Þdz hqx
ð12Þ
0
_ Y c jz ¼ z0 i is estimated by filtering a 1-D laminar prewhere hqx mixed flamelet computed at a given mixture fraction z:
_ Y c jz ¼ z0 i ¼ hqx
Then, under these modeling assumptions, the propagation speed of a filtered flame front in a stratified mixture is given by:
SD ¼ N e Sl
fc Y @q e_ gc q fc Þ ¼ r ðqDrY c Þ r ðq fc Þ þ q ~Y ~Y u uY u x þ r ðq Yc @t ð11Þ
Z
þ1
1
q x_ Y c x0n ; z0 GD xn x0n dx0n
ð13Þ
In this article, P(z) is modeled by a b-function parametrized by the 002 . first- and second-order moments of the mixture fraction: ~z and zf e _ The filtered reaction rate x Y c is in practice pre-computed and e_ Y ¼ stored in a four-dimensional look-up database, x c e f f 002 x_ Y c ½ Y c ; ~z; Sz ; D where Sz ¼ z =ð~zð1 ~zÞÞ is the unmixedness factor, 002 ¼ 0) to unity going from zero in perfectly mixed situations ( zf when reactants are fully segregated (bimodal distribution z = 0 or z = 1). 3.2. Filtered laminar diffusion terms A model for the filtered molecular diffusion term has been discussed in [1] for premixed flames. The extension toward stratified combustion reads as follows:
fc Þ r ðqDrY c Þ ¼ r ðaY c q0 D0 r Y
ð14Þ
where q0 and D0 are reference values for the density and the molecular diffusion coefficient, respectively, for example, corresponding to fresh gas values. The correction factor aY c is estimated from 1D filtered premixed flames:
2707
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
aY c ½ Yfc ; ez ; Sz ; D ¼
c q D @Y @xn
ð15Þ
q0 D0 @@xYecn
e c =Y fc eq in Fig. 2 for Sz = 0 (perfect mixing) and a function of ~z and Y Sz = 0.02, respectively. 3.3. Unresolved transport term
The laminar diffusivity is computed assuming unity Lewis number, D ¼ k =q C p . In Eq. (15), q D @Y c =@xn is computed as follows:
q
@Y c
D
@xn
¼
Z
1
hq D
0
@Y c
@xn
0
0
0
jz ¼ z iPðz Þdz
ð16Þ
fc =@xn presents in the RHS denominator of The computation of @ Y Eq. (15), is decomposed as follows:
fc @Y @ q Y c ¼ @xn @xn q
!
The unresolved transport term is modeled as [1]:
1 @ q Y c f @ q Yc ¼ q @xn @xn
! ð17Þ
where
Z 1 @ q Y c @Y ¼ q c jz ¼ z0 Pðz0 Þdz0 @xn @xn 0 Z 1 @ q @q 0 ¼ jz ¼ z0 Pðz0 Þdz @xn @xn 0
ð18Þ ð19Þ
fc ; ~z; Sz ; D is As for the progress variable reaction rate, aY c ¼ aY c ½ Y stored in the look-up database. This correction factor is plotted as
gc q fc Þ fc ¼ NXY þ r ðaY ðN 1Þq D0 r Y ~Y uY u r q c c 0
ð20Þ
where the first RHS term, NXY c , corresponds to the thermal expansion and is estimated from 1-D detailed chemistry laminar premixed flame solutions as [1]:
XY c ¼ q0 ðzÞSl ðzÞ
fc @Y c @Y q0 ðzÞSl ðzÞ @xn @xn
ð21Þ
while the second RHS term in Eq. (20) is a gradient closure of the unresolved transport due to turbulence motions. Eqs. (12), (14), and (20) are designed to propagate a flame front at the turbulent flame speed SD ¼ Ne S l , according to Eq. (9). The first RHS term in Eq. (21) is estimated according to the procedure described in Appendix B.2:
q0 ðzÞSl ðzÞ
@Y c ¼ @xn
Z 0
1
q0 ðz0 ÞSl ðz0 Þ
@Y c 0 jz ¼ z0 Pðz0 Þdz @xn
ð22Þ
e c =Y fc eq and D ¼ 20dst for two values of mixture fraction unmixedness factor (dst is the flame thickness determined as the Fig. 2. Correction factor aY c as a function of ~z and Y r r full width at half maximum of the reaction rate in the stoichiometric conditions).
Fig. 3. 1-D simulations of stratified flame with F-TACLES for ~z ¼ zst . (a) Temperature profiles. (b) Flame speed; j: effective propagation speed, computed from Eq. (24); –: theoretical propagation speed given by Eq. (10).
2708
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
The second term is computed as follows:
q0 ðzÞSl ðzÞ
fc @ Y fc Z 1 @Y ¼ q ðz0 ÞSl ðz0 ÞPðz0 Þdz0 @xn @xn 0 0
5. LES of a 3-D turbulent stratified swirled flame
ð23Þ
fc =xn is estimated from Eq. (17). In practice, the subgrid with @ Y scale convection term of the progress variable is also stored in a fc ; ~z; Sz ; D. 4-D look-up database: XY c ¼ XY c ½ Y
4. Construction and validation of the filtered look-up table 4.1. Chemical database generation A chemical look-up table for methane/air combustion, representative of the conditions operated in the MOLECULES experimental setup [23–25], described in Section 5.1, is constructed. For that purpose, a collection of 1-D laminar premixed flames is computed using the GRI 3.0 mechanism [33] and the PREMIX [34] solver for equivalence ratio within the flammability limits. The unity Lewis assumption is retained for this study. The progress variable Yc is defined as Y c ¼ Y CO þ Y CO2 as in [28]. The mixture fraction z is defined from nitrogen assumed as inert as the NOx chemistry is not taken into account (Y z ¼ Y N2 ). This coordinate is normalized using Eq. (2). Y fz and Y ox denote values of z Y N2 in fuel stream (here pure methane at 388 K, Y fz ¼ 0) and oxidizer stream (here air at 623 K, Y ox z ¼ 0:77), respectively. From this flamelet library, the procedure proposed in Appendix B is followed to compute the filtered thermochemical quantities. The grid resolution of the filtered look-up table is 200 150 10 in e c , ~z, Sz) subspace. the ( Y
5.1. MOLECULES configuration The combustion chamber MOLECULES, designed to be representative of an actual aeronautical combustor, has been experimentally investigated by Janus et al. [23–25] at TU Darmstadt. The geometry shown in Fig. 4 consists of a plenum, a swirl-injector, and a combustion chamber. Preheated air is introduced in the swirler via the plenum. Fuel is directly injected through the center of the swirler. The combustion chamber pressure is 2 bar at the operating conditions. Air is injected at 623 K at a flow rate of 30 g/s. The global equivalence ratio is 0.8. The fuel used for the experiment is natural gas identified to pure methane in numerical simulations. The flow field has been characterized by Laser Doppler Anemometry for both non-reactive and reactive cases. The mixing has been
4.2. 1-D flame computations As explained in Section 2.2, the present modeling approach describes a flame front evolving in a stratified and turbulent environment at a velocity SD ¼ N Sel , where N is the subgrid flame wrinkling factor. As iso-c surfaces are assumed perpendicular to iso-z surfaces at the subgrid scale level, the propagation speed Sel of the filtered flame front is given by Eq. (10). To validate the model implementation, a series of planar flames evolving in a stratified mixture characterized by a constant subgrid scale variance are computed. The database described in Section 4.1 is used st with a filter size equal to D ¼ 20dst r (dr is the full width at half maximum of the progress variable reaction rate of the stoichiometric flame). A premixed flame under stoichiometric conditions without mixture fraction fluctuations at subgrid scale is computed as a reference. Then, 1-D stratified flames are computed with ~z ¼ zst and 0 < Sz < 1. Temperature profiles are plotted in Fig. 3a and the corresponding flame burning velocities are shown in Fig. 3b. Increasing the subgrid scale mixture fraction variance decreases the burnt gas temperature and increases the flame thickness as shown in Fig. 3a. The predicted flame propagation speed is computed after post-processing 1-D solution profiles data
1 Sel ¼
q0
Z
1
e_ ðxÞdx qx c
ð24Þ
0
and is compared with the flame speed given by Eq. (10). Fig. 3b shows that the effective flame propagation matches the one predicted by the model theory. The model implementation procedure is therefore validated. These tests also demonstrate that the filtered look-up table resolution is sufficient to capture the correct flame propagation speed.
Fig. 4. Schematic views of the MOLECULES burner.
Table 1 Computational parameters used for the simulations. MOLECULES combustion chamber Pressure Air temperature Air mass flow rate Global equivalence ratio Fuel temperature Fuel mass flow rate Reynolds number for air Laminar flame thicknesses Methane injector Swirler exit Chamber dimension
2 bar 623 K 30 g/s 0.8 388 K 1.4 g/s 46,000 5104 m < dc <1.5103 m 4 mm diameter (40 nodes in the radial direction) 28 mm diameter (112 nodes in the radial direction) 60 60 80 mm3
2709
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
measured by Planar Laser Induced Fluorescence (PLIF) technique using acetone as a fuel tracer. In addition, the flame front has been identified with PLIF of the OH radical. As fuel and oxidizer are not premixed before entering the combustion chamber, the modeling approach retained has to account for the mixture fraction stratification. The numerical setup is described in the following section. 5.2. Numerical setup The F-TACLES approach is implemented in the compressible LES code AVBP [35] using the Tabulated Thermochemistry for Compressible flow (TTC) formalism [36]. The filtered balance equation of the progress variable implemented in AVBP reads as follows: Fig. 5. Line positions where numerical simulations are compared with experimental data.
fc Y @q fc Þ fc ; ~z; Sz ; Dq D0 r Y fc Þ ¼ r ðNaY ½ Y ~Y u þ r ðq c 0 @t e_ ½ Y fc ; ~z; Sz ; D þ Nq fc ; ~z; Sz ; D ð25Þ x NXY c ½ Y Yc where the wrinkling factor N is estimated by the analytical model proposed by Charlette et al. [37]:
140
"
120
N¼
Axial velocity
100 80 60 40 20 0 -20 0
20
x [mm]
40
60
Fig. 6. Mean axial velocity component along the centerline for the non-reactive case. Symbols: measurements. Lines: simulation with F-TACLES.
D u0 1 þ min 0 ; C 0D dl Sl
#!0:5 ð26Þ
In this expression, d0l denotes the thermal flame thickness, u0D the subgrid-scale turbulent velocity, and C the spectral efficiency function (see [37]). In order to identify an adequate computational grid, LES of the cold flow has been initially performed on a coarse mesh. Then, the grid resolution has been increased until the convergence of the mean flow field. An unstructured 40 million cells mesh has then been identified and retained for the reactive flow simulation. The fuel nozzle is meshed with 40 cells along the diameter direction. st The filter size D ¼ 20dst r , where dr is the full width at half maximum of the reaction rate of the stoichiometric 1-D premixed flamelet, is retained to generate the filtered look-up table. It corresponds to a filter width 8 times larger that the laminar thermal layer. This choice ensures a computation of the filtered flame front over 5 nodes, which is the minimal resolution required to
Non reactive case 25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-20
-20
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
y [mm]
25
y [mm]
25
0
100
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
Fig. 7. Mean (left) and RMS (right) of axial velocity for the non- reactive case. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
2710
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Non reactive case x = 1mm
x = 5mm
x = 10mm
x = 15mm
x = 20mm
x = 5mm
x = 1mm
x = 10mm
x = 15mm
x = 20mm
25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-20
-20
-25
0 0.5
1
-25
0 0.5
1
-25
0
0.5
-25
0
0.5
-25
y [mm]
25
y [mm]
25
0
0.5
-25
0
0.3
-25
0
0.3
-25
0
0.3
-25
0
0.3
-25
0
0.3
Fig. 8. Mean (Left) and RMS (Right) of CH4 mass fraction for the non-reactive case. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
avoid numerical artifacts. More details and justifications may be found in Ref. [11]. The third-order finite element TTGC scheme [38] is used. Boundary conditions are prescribed using Navier–Stokes Characteristic Boundary Conditions [39]. The WALE turbulence model [40] closes the unresolved shear stress tensor. Boundary parameters are summarized in the Table 1 for both reactive and non-reactive simulations. For velocity field and fuel mass fraction, numerical simulations are compared with experimental data along the centerline and in the plane f = 0 at the locations indicated in Fig. 5.
5.3. Non-reactive case The non-reactive case is first computed. Mean and resolved Root Mean Square (RMS) quantities are computed by time averag140 120
5.4. Reactive case
100
Axial velocity
ing LES solutions over a physical time corresponding to 3 flow through times based on the air inlet velocity (14 ms). This observation time ensures sufficient convergence of the first-order statistical moments on the planes x = 1 mm to x = 20 mm, experimentally investigated by Janus et al. [23–25]. Mean axial velocity component is first compared with the experimental data along the centerline in Fig. 6. The fair agreement between measurements and simulation shows that the recirculation zone is well captured by the LES. Mean radial profiles of axial velocity and CH4 mass fraction are then plotted on Fig. 7 (left) and Fig. 8 (left), respectively. A correct agreement is observed between experimental and numerical profiles, showing that the mixing and the flow field are correctly predicted. RMS of axial velocity and CH4 mass fraction are plotted on Fig. 7 (right) and Fig. 8 (right), respectively. The general trends are well reproduced but an overestimation of the fluctuation amplitudes of CH4 by the simulations is observed at x = 5 mm and x = 10 mm (the relative differences between measured and computed data are about 30 %).
80 60 40 20 0 -20 0
20
x [mm]
40
60
Fig. 9. Mean axial velocity component along the centerline for the reactive case. Symbols: measurements. Lines: simulation with F-TACLES.
For the reactive case, mean and resolved root mean square (RMS) quantities are computed by time averaging LES solutions over a physical time corresponding to 8 flow through times (38 ms) based on the air inlet velocity. The performances of the F-TACLES model are assessed by comparing the predicted velocity and fuel mass fraction to the experimental data. Mean axial velocity component is compared with the experimental data along the centerline in Fig. 9. As for the cold case, the LES results match very well the velocity measurements. Mean axial and radial velocities are plotted in Fig. 10 (top) and Fig. 11 (top), respectively. A fair agreement is observed between experimental and numerical profiles, showing that the flow field is very well predicted in the reactive case. The corresponding RMS are plotted in Fig. 10 (bottom) and Fig. 11 (bottom). As these resolved RMS do not include the subgrid scale RMS’s conclusions regarding the model performance
2711
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
y [mm]
Reactive case 25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
y [mm]
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
0
100
-25
25
25
25
25
25
25
25
25
20
20
20
20
20
20
20
20
15
15
15
15
15
15
15
15
10
10
10
10
10
10
10
10
5
5
5
5
5
5
5
5
0
0
0
0
0
0
0
0
-5
-5
-5
-5
-5
-5
-5
-5
-10
-10
-10
-10
-10
-10
-10
-10
-15
-15
-15
-15
-15
-15
-15
-15
-20
-20
-20
-20
-20
-20
-20
-20
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
50
-25 0
0
100
50
Fig. 10. Mean (top) and RMS (bottom) of axial velocity in m s1. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at locations indicated in Fig. 5.
in terms of flame turbulence interactions are not straightforward. However, it is observed that LES RMS remains lower than measured RMS, as expected from theory, and that the evolution of the predicted RMS is very similar to the experimental results. Mean and RMS CH4 mass fraction are plotted in Fig. 12. The comparison shows that the mixture field is predicted in the reactive case excepted at x = 5 mm and x = 10 mm where deviations of about 30% are observed. A comparison of the predicted mean OH mass fraction versus the experimental PLIF signal is shown in Fig. 13. Both the lift-off height and the flame envelope are captured by the simulation.
The present modeling strategy relies on the assumption that the stratified flame is composed of a collection of flamelets that have a premixed-like structure. To investigate this hypothesis, scatter plots of filtered temperature versus the filtered mixture fraction, ~z, are shown in Fig. 14a where each computational node is represented by one dot. Statistics shown here are estimated from a single snapshot but the same conclusions are drawn when other time instants are analyzed. It shows that all the possible phase subspace is covered but no preferential flame structure is evidenced. Fig. 14b e ) colored by computational shows the reduced phase subspace (e z; T nodes population. This post-processing highlights two major
2712
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Reactive case 20
20
20
15
15
15
10
10
10
20
25 20
15
15
y [mm]
5
5
0
0
0
-5
-5
-5
-10
-10
-10
5
5
0
0 -5
-5
-10
-15
-20
-50 0
50
-15
-20
-50 0
50
20
20
15
15
15
10
10
10
-15
-20
-20
20
-50 0
-25
50
-50
0
50
-20
-50
0
50
20
25 20
15
15
y [mm]
5
5
0
0
0
-5
-5
-5
-10
-10
-10
5
5
0
0 -5
-5
-10
-20
-15
0
30
-20
-15
0
30
-20
20
20
15
15
15
10
10
10
5
5
5
0
0
0
-5
-5
-5
-10
-10
-10
-15
-15
-15
-20
-20
-20
-25
-50
0
50
-25
-50
0
50
-25
-50
25
25
25
20
20
20
15
15
15
10
10
10
5
5
5
0
0
0
-5
-5
-5
-10
-10
-10
-15
-15
-15
-20
-20
-20
0
50
-10 -15
-15
20
10 10
5
25
-10 -15
-15
25
10 10
5
25
-15
-20
0
-25
30
Fig. 11. Mean (top) and RMS (bottom) of radial velocity in m s locations indicated in Fig. 5.
0
30
-20
0
30
-25
0
30
-25
0
30
-25
0
30
1
. Symbols: measurements. Lines: simulation with F-TACLES. Profiles are extracted in the f = 0 plane at
trajectories: a line (1) and a narrow vertical strip (2). The line (1) shows that fuel and oxidizer are mixed together prior the chemical reactions. Indeed the transition between fresh and burnt gases occurs in the region (2) for mixture fractions that belong to the flammability limits. The same conclusion is observed in Fig. 15 where the joint filtered density function of temperature and mixture fraction is plotted. Statistical events are clustered in the vicinity of z = 0.05 plane. It indicates that the transition between fresh and burned gases follows premixed flames trajectories. Fig. 16 shows a scatter plot of the subgrid scale flame wrinkling factor N encountered in the flame front reaction zone. As N ranges
from 1 to 3, the impact of the subgrid turbulence on the flame propagation is important and therefore the flame front displacement speed is sensitive to the subgrid scale wrinkling factor model given in Eq. (26). A second reactive simulation has been performed without accounting for subgrid scale mixture fraction heterogeneities, that is, by modeling P(z) by a Dirac d function centered on ~z ðPðzÞ ¼ dðz ~zÞÞ. Figure 17 shows two 2-D snapshots of the instantaneous temperature iso-color field with P(z) modeled by a b-function (left) or by a Dirac function (right). The temperature, computed from the filtered chemical look-up table, is strongly
2713
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Reactive case
y [mm]
y [mm]
x = 1mm
x = 5mm
x = 10mm x = 15mm x = 20mm x = 30mm x = 40mm
20
20
20
20
20
20
20
10
10
10
10
10
10
10
0
0
0
0
0
0
0
-10
-10
-10
-10
-10
-10
-10
-20
-20
-20
-20
-20
-20
-20
0 0.5 1
0 0.5 1
x = 1mm
x = 5mm
0
0.6
0
0.4
0
0.3
0 0.03
0 0.1
x = 10mm x = 15mm x = 20mm x = 30mm x = 40mm
20
20
20
20
20
20
20
10
10
10
10
10
10
10
0
0
0
0
0
0
0
-10
-10
-10
-10
-10
-10
-10
-20
-20
-20
-20
-20
-20
-20
0 0.2
0 0.2
0
0.3
0
0.3
0
0.3
0
0.1
0 0.04
Fig. 12. Mean (top) and RMS (bottom) of CH4 mass fraction. Symbols: measurements. Lines: simulation with F-TACLES.
grid scale mixture fraction effects, the front displacement speed is too high. The flame is then stabilized at the burner lips and is not lifted anymore. Then, because of a misprediction of the flame anchoring mechanism, accurate statistics on temperature and velocity cannot be obtained.
5.5. Comparison with other modeling strategies applied in the MOLECULES combustor
Fig. 13. Iso-color of simulated (left) and measured (right) OH concentrations in the combustor. OH is made dimensionless by its maximal value.
dependent of the mixture fraction variance. As expected, peaks of temperature are smoothed when the impact of subgrid mixture fraction variance is considered. Also a strong influence of P(z) on the flame dynamics is observed. Indeed, without considering sub-
Numerical investigations of the MOLECULES combustion chamber have been performed in the past by other groups following different modeling strategies. U-RANS (Unsteady Reynolds Average Navier-Stokes) simulations have been conducted by Schneider et al. [41] with two turbulent combustion models: the G-equation [12] and the Bray-Moss-Libby approaches [42]. Results do not show specific advantage of one or other premixed combustion model in a RANS context and fair comparisons versus experiments has been achieved. However, the comparison between measured
2714
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Fig. 14. Left: Scatter plot of the filtered temperature versus the filtered mixture fraction in the MOLECULES case. The red line corresponds to the global equivalence ratio. Right: Combustion trajectories in the MOLECULES burner. Red color corresponds to a number of grid points larger than 10,000 and the blue lower than 50.
Fig. 15. Joint filtered density function of the mixture fraction and temperature.
e c =Y e eq Fig. 16. Scatter plot of the subgrid flame wrinkling for 0:1 < Y c < 0:9.
Fig. 17. 2-D snapshots of the temperature field. The Filtered Density Function P(z) is either modeled by a b function (left) or by a Dirac d function (right). The flame is then stabilized at the burner lips when the subgrid scale mixture fraction heterogeneities are not accounted for (right). The solid black line shows the stoichiometric line.
and simulated 2-D field of OH radical shows that the present LES results obtained with F-TACLES significantly improves the prediction of the flame front position. A possible explanation of the difference is the representation of the turbulence. Indeed, analyzing velocity data [41], the authors highlight the difficulty to accurately describe the turbulence with U-RANS. As shown both in the cold flow LES of the combustor performed by Sadiki et al. [43] and in the present non-reactive case simulation, a better representation of the flow dynamic is obtained with LES. LES of the reactive flow in the MOLECULES combustor has also been performed by Sadiki et al. [43] using two different turbulent combustion models: the RPV (Reaction Progress Variable) model and the steady diffusion flamelet model [12]. The diffusion flamelet model failed to reproduce the flame lift-off because the chemistry tabulation is not adapted to the stratified flame structure. In the RPV approach, the chemistry is tabulated from 1-D premixed flamelets as in FPI (Flame Prolongation of ILDM) [27] or FGM (Flamelet Generated Manifold) [29] approaches. However, unlike to F-TACLES, the filtered density function of the progress variable is simply
2715
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
modeled by a Dirac function. Subgrid scale flame front contributions are not accounted for and may cause errors in the flame front propagation. It may explain why predictions are closer to the experimental data with F-TACLES than with RPV.
6. Conclusions Issues relative to LES of stratified combustion are discussed in this article. An extension of the F-TACLES model, initially developed for premixed combustion, is proposed to tackle turbulent stratified flames. Detailed chemistry effects are accounted for in an LES context by storing in a look-up table a set of 1-D laminar filtered premixed flames. A mathematical formalism, including the impact of subgrid scale mixture fraction heterogeneities and subgrid scale flame wrinkling on the flame front propagation, is proposed to couple this filtered thermo-chemical table to the LES flow governing equations. The model is implemented in the compressible AVBP solver [35] to perform the 3-D simulations of a stratified swirled turbulent flame representative of an aeronautical combustor. Both reactive and non-reactive simulations are compared with experimental data. In general, a very good agreement between numerical simulations and measurements is observed. A comparison of the simulations performed with other techniques published in the literature demonstrates the performances of the F-TACLES approach.
e_ ¼ q x Yc
Appendix A. Modeling assumptions in F-TACLES formalism
jrcj
h/i ¼
Z
1
R
1
The filtered reaction rate of the progress variable is given by:
e_ ¼ q x Yc
Z
qðx0 Þx_ Y c ðx0 ÞGD ðx x0 Þdx0
ðA:1Þ
X
_ Y c ðx0 Þ is the progress variable reaction rate at location x0 where x and GD is the LES filter. This reaction rate may be rewritten as follows:
e_ ¼ q x Yc
Z
1
Z
qðx0 Þx_ Y c ðx0 Þ
0
jrcjdðcðx0 Þ c0 ÞGD ðx x0 Þdx0 dc jrcj 0 X Z 1 qx_ Y c ¼ Rðc0 Þdc0 ðA:2Þ jrcj c0 0
where d(c) is the Dirac delta-function and Rðc0 Þ ¼ jrcjdðc c0 Þ the density of the c = c0 -isosurface, that is, the available c = c0 -surface per unit filter volume. h/ic0 denotes the average of the quantity / over the c = c0 -surface in the filtering volume, defined as [44]:
R /jrcjdðcðx0 Þ c0 ÞGD ðx x0 Þdx0 RX / ¼ h ic0 0 0 0 0 X jrcjdðcðx Þ c ÞGD ðx x Þdx
ðA:3Þ
Following [44], a generalized surface density, R, is then introduced:
R¼
Z 0
1
Rðc0 Þdc0 ¼
Z
1
0
jrcjdðc c0 Þdc ¼ jrcj
0
and the filtered reaction rate expressed as follows:
ðA:4Þ
Rðc0 Þdc0 ¼
qx_ Y c
c0
jrcj
R¼
qx_ Y c jrcj
Njrcj ðA:5Þ
0
R /jrcjGD ðx x0 Þdx0 0 h/ic0 Rðc0 Þdc ¼ RX 0 0 X jrcjGD ðx x Þdx
ðA:6Þ
_ Y c =jrcji measures the mean reaction rate per unit of In Eq. (A.5), hqx flame area. The above relations are exact, according to the corresponding definitions. F-TACLES model assumptions are now discussed for two cases successively: (i) perfectly premixed combustion (i.e., the original F-TACLES model [1]); (ii) stratified combustion, an extension proposed in this paper. A.2. Perfectly premixed combustion By construction, the F-TACLES model is a flamelet model: the flame surface is identified to a collection of laminar flame elements. Surface densities R(c0 ) are assumed to be independent on c0 and the progress variable reaction rate per unit of flame surface, _ Y c =jrcji, is computed from one-dimensional laminar flame, hqx using Eq. (A.6):
q x0n x_ Y c x0n GD ðxn x0n Þdx0n
1
R þ1
jrcj
R þ1
0
j@c =@xn jGD ðxn x0n Þdxn Z þ1 1 q x0n x_ Y c x0n GD ðxn x0n Þdx0n j@c =@xn j 1 1
ðA:7Þ
where the ⁄ superscript denotes quantities issued from a 1-D laminar unstretched premixed flame computed in the xn direction. According to Eq. (A.5), the filtered progress variable reaction rate in turbulent flames is modeled as follows:
e_ Y q x c
Z
þ1
1
A.1. Definitions and general relations
where N ¼ jrcj=jrcj is the surface wrinkling factor (i.e., the ratio of the generalized surface divided by the resolved surface) and h/i the generalized surface average defined by:
qx_ Y c
We thank Prof. Andreas Dreizler and Dr. Simone Schoenfelder for providing us the burner geometry and the experimental data. This work was supported by the ANR-07-CIS7–008-04 Grant of the French Ministry of Research and accessed the high-performance computing resources of IDRIS and CINES under the allocation 2009-i2009020164 made by GENCI (Grand Equipement National de Calcul Intensif).
qx_ Y c
0
Acknowledgments
Z 1
q x0n x_ Y c x0n GD ðxn x0n Þdx0n N
ðA:8Þ
assuming that the resolved flame thickness are similar in 1-D laminar flames and in turbulent flames (i.e., j@c =@xn j jrcj). Eq. (A.8) is the F-TACLES model expression for the filtered progress variable reaction rate [1]. The quantity inside brackets is computed from 1-D laminar flames and tabulated as function of the filtered progress e c and the filter size D. variable Y A.3. Stratified combustion In stratified combustion regime, the reaction rate of the filtered progress variable is influenced by the mixture heterogeneities at the subgrid scale. The challenge is now to estimate the reaction _ Y c jrcji rate of the filtered progress variable per unit surface, hqx (Eq. (A.5)), in this situation. By definition:
qx_ Y c jrcj
R
¼
X
qðx0 Þx_ Y c ðx0 ÞGD ðx x0 Þdx0 R
R1 R ¼
0
R1
¼ R 10 0
X
jrcjGD ðx x0 Þdx0 0
0 _ 0 0 0 0 0 X qðx ÞxY c ðx Þdðzðx Þ z ÞGD ðx x Þdx dz R1 R 0 0 0 0 0 0 X jrcjdðzðx Þ z ÞGD ðx x Þdx dz
_ Y c jz ¼ z0 ÞPðz0 Þdz0 ðqx 0
ðjrcjjz ¼ z0 Þ Pðz0 Þ dz
ðA:9Þ
where z is the mixture fraction going from z = 0 for pure oxidizer to z = 1 for pure fuel. P(z0 ) is the subgrid scale filtered density function of the mixture fraction, defined by [45]:
2716
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
Pðz0 Þ ¼
Z
dðzðx0 Þ z0 ÞGD ðx x0 Þ dx0
ðA:10Þ
X
The factor aY c correcting the filtered laminar diffusion term (Eq. (15)):
ð/jz ¼ z0 Þ is the filtered value of / conditioned on the iso-surface z = z0 :
R
ð/jz ¼
z0 Þ
¼
Þdðzðx0 Þ z0 ÞGD ðx x0 Þdx0 0 Þ z0 ÞG ðx x0 Þdx0 dðzðx D X
/ðx XR
0
_ Y c jz ¼ ðqx
Z
All these terms can be expressed in the form:
¼ Z ðjrcjjz ¼ z0 Þ ¼
þ1
1 þ1
1
ðx0n ; z0 Þ
0 0 Y c ðxn ; z ÞGD ðxn
q x_
@c 0 0 0 0
@x ðxn ; z Þ GD xn xn dxn n
x0n Þ
0 dxn
ðA:12Þ ðA:13Þ
Eq. (A.9) becomes R 1 R þ1 0 0 0 0 q ðx ; z Þx_ Y c ðxn ; z ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0 qx_ Y c
hR n
i 0 R1
1 þ1 @c 0 0 jrcj
ðx0 ; z0 Þ G ðx x0 Þ dx Pðz0 Þ dz 0
1 jrcj
1
Z
1
Z
0
@xn
þ1
1
D
n
n
n
q ðx0n ; z0 Þx_ Y c ðx0n ; z0 ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0
Assuming that jrcj extracted from one-dimensional laminar flame computations combined with presumed filtered density functions P(z0 ) is similar to jrcj in turbulent flame simulations (i.e., jrcj jrcj) allow to express the filtered reaction rate of the progress variable (Eq. (A.5)) as follows:
e_ N q x Yc
Z
0
þ1
1
Z
1
0
Z
q ðx0n ; z0 Þx_ Y c ðx0n ; z0 ÞGD ðxn x0n Þ dx0n Pðz0 Þ dz0
ðA:15Þ The quantity inside brackets is computed from 1-D laminar flames while the filtered density function P(z0 ) of the mixture fraction is modeled as a b-function and depends on the mean, ~z, and the 002 =~ unmixedness factor, Sz ¼ zf zð1 ~zÞ, of the mixture fraction z. The quantity inside parenthesis is then tabulated as function of e c , the filter size D, ~z and Sz, while the filtered progress variable Y the flame surface wrinkling factor, N, quantifies the flame/turbulence interactions.
B.1. Introduction This Appendix describes the practical construction of the filtered database used in the F-TACLES model. According to the model derivation, three terms are estimated from one-dimensional flame computations and then stored in the database: The contribution to the filtered reaction rate (Eq. (A.15)): 1
0
hqxY c jz ¼ z0 iPðz0 Þdz
0
ðB:1Þ
The contribution of the thermal expansion to the unresolved transport (Eqs. (21)–(23)):
XY c ¼
Z
A collection of 1-D premixed laminar flame is first computed with a detailed chemical scheme for equivalence ratios covering the mixture flammability limits. Thermo-chemical data u, such as mass fractions and reaction rates and solutions of the 1-D flam elets, are mapped in the 2-D x0n ; z0 subspace where x0n denotes the direction normal to the flamelet and z0 is the mixture fraction:
u ¼ u x0n ; z0
0
q
0
¼ z0 iPðz0 Þdz
@xn fc Z @Y @xn
ðB:6Þ
Outside the flammability range, chemical reaction rates are set to zero while species mass fractions are linearly interpolated between the point of the last tabulated flame on the lean (respectively rich) side and the air (respectively fuel) stream. For a given value of the filter size D, conditional filtered quantities are then computed for constant mixture fraction value z = z0 according to:
hujz ¼ z0 iðxn ; z0 ; DÞ ¼
Z
þ1
1
uðx0n ; z ¼ z0 ÞGD xn x0n dx0n
ðB:7Þ
where GD denotes the one-dimensional filter function, chosen here as Gaussian. The practical integration of the RHS in Eq. (B.7) is detailed in Appendix B.3. According to Eq. (1), conditional values of the progress variable hYcjz = z0 i and of the normalized progress variable hcjz = z0 i are related through the relation:
hY c jz ¼ z0 i Y fc ðz0 Þ Y bc ðz0 Þ Y fc ðz0 Þ
ðB:8Þ
which is used to map the conditional variables huj z = z0 i(xn, z0 , D) into the subspace (hcjz = z0 i,z0 , D) where hcjz = z0 i replaces the space coordinate xn. P(z0 ) is now presumed by a b-function P b ðz0 ; ~z; Sz Þ parametrized by the filtered mixture fraction ~z and the corresponding unmixed002 =~ ness factor Sz ¼ zf zð1 ~zÞ. Eq. (B.5) is then recast as follows: 002 ; DÞ ¼ ðc; ~z; zf u
Z
1
0 hujz ¼ z0 iðhcjz ¼ z0 i; z0 ; DÞPb ðz0 ; ~z; Sz Þ dz
ðB:9Þ
0
where
Z
1
0 hcjz ¼ z0 iPb ðz0 ; ~z; Sz Þ dz
ðB:10Þ
0
@Y c
0 0 0 ðz ÞSl ðz Þh
ðB:5Þ
and are easily extracted from one-dimensional laminar flames. When the mixture fraction is constant: P(z0 ) = d(z0 z0), where d is the Dirac delta-function. Then, only one flame is considered for a given mixture fraction value z0, as in the original F-TACLES model presented in Appendix A.2. For variable mixture fractions, several 1-D flames are combined through the filtered density function P(z0 ) as described in the following.
c ¼ 1
hujz¼z0 iðxn Þ
hcjz ¼ z0 i ¼
Z
u x0n ; z0 GD ðxn x0n Þ dx0n Pðz0 Þ dz0
1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
Appendix B. Construction of the filtered database
q x_ Y c ¼
þ1
B.2. General procedure
ðA:14Þ
1
¼ u
n
Z
ðB:3Þ
q0 D0 @@xYecn
ðA:11Þ
_ Y c jz ¼ z0 Þ, as well as The filtered progress variable reaction rate, ðqx the gradient of the progress variable, ðjrcjjz ¼ z0 Þ, both conditioned on the z = z0 isosurface, should now be determined (see Eq. (A.9)). F-TACLES assumes that these terms can be estimated from a one-dimensional laminar flame at constant mixture fraction z = z0 . Then,
z0 Þ
c q D @Y @xn
aY c ¼
jz
0
1
q0 ðz0 ÞSl ðz0 ÞPðz0 Þdz0
ðB:2Þ
The value of the conditional normalized progress variable hcjz = z0 i should now be specified for each mixture fraction value z0 to integrate Eqs. (B.9) and (B.10). Assuming that iso-c and iso-z surfaces are orthogonal at the subgrid scale level, the conditional filtered
P. Auzillon et al. / Combustion and Flame 159 (2012) 2704–2717
GD ðxÞ ¼
6 p D2
1=2
6x2 exp 2 D
2717
ðB:18Þ
References
Fig. 18. Scheme for Gaussian convolution.
normalized progress variables hcjz = z0 i is kept constant during the integration. Note that this procedure corresponds to the following assumption: P(c, z) = P(c)P(z), where P(c, z) is the joint progress variable/mixture fraction FDF. Following the same methodology, Favre filtered thermochemical quantities are computed from the definiu ~ ¼ qu. In particular, the Favre filtered progress variable is tion q estimated by:
Z
fc ¼ 1 Y q
1
0 hqY c jz ¼ z0 iPb ðz0 ; ~z; Sz Þ dz
ðB:11Þ
0
e c , ~z, zf 002 , and u ~ are then mapped in the 4-D ( Y Filtered quantities u D) subspace to terminate the filtered chemical table generation. A constant value of D is chosen in the present article and balance 02 are solved to directly couple the filtered e c , ~z and zf equations for Y look-up table to the LES flow governing equations. B.3. Numerical integration of Eqs. (B.7) and (B.9) The F-TACLES model predictions are very sensitive to the quality of the numerical computations of (B.7) and (B.9). The integration of the b-FDF (filtered density function) P(z) in Eq. (B.9) is performed following the methodology proposed in [46]. The progress variable c(x) is a monotonic function that evolves between 0 (fresh gases) and 1 (burnt gases) as shown in Fig. 18. x0 and x1 are introduced such as c(x0) = 0 and c(x1) = 1. Thus, c is a bijection defined on [x0, x1] to [0, 1]. In order to compute the integral in Eq. (B.7), the flame is divided in three zones as shown in Fig. 18:
hujz ¼ z0 i ¼ I0 þ I þ I1
ðB:12Þ
where I0, I and I1 are as follows:
Z x0 I0 ¼ uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0 1 Z x1 uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0 I¼
ðB:13Þ ðB:14Þ
x0
I1 ¼
Z
þ1
uðz ¼ z0 ; x0 ÞGD ðx x0 Þdx0
ðB:15Þ
x1
I is numerically computed with a second-order integration scheme. The range ]1, x0] (resp. [x1, +1[) corresponds to fresh gases (resp. burnt gases); therefore, u(z = z0 , x0 ) is constant and equal to u(z = z0 , c = 0) (resp. u(z = z0 , c = 1)). I0 and I1 are then computed directly from the analytical solution:
!# pffiffiffi 6 1 erf I0 ¼ ðx x0 Þ 2 D " !# pffiffiffi 6 uðz ¼ z0 ; c ¼ 1Þ ðx x1 Þ I1 ¼ 1 erf 2 D
uðz ¼ z0 ; c ¼ 0Þ
"
the 1-D filter GD being defined by a Gaussian function:
ðB:16Þ ðB:17Þ
[1] B. Fiorina, R. Vicquelin, P. Auzillon, N. Darabiha, O. Gicquel, D. Veynante, Combust. Flame 157 (2010) 465–475. [2] A. Unal, Y. Hu, M. Chang, M.T. Odman, A. Russell, Atmos. Environ. 39 (3) (2005) 5787–5798. [3] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, R.T. Edwards, Inc., 2005. [4] H. Pitsch, Combust. Flame 143 (4) (2005) 587–598. [5] J. Galpin, A. Naudin, L. Vervisch, C. Angelberger, O. Colin, P. Domingo, Combust. Flame 155 (2008) 247–266. [6] C. Duwig, Combust. Theor. Modell. 13 (2) (2009) 251–268. [7] S. Menon, W. Jou, Combust. Sci. Technol. 75 (1991) 53–72. [8] O. Colin, F. Ducros, D. Veynante, T. Poinsot, Phys. Fluids 12 (7) (2000) 1843– 1863. [9] M. Boileau, G. Staffelbach, B. Cuenot, T. Poinsot, C. Berat, Combust. Flame 154 (1–2) (2008) 2–22. [10] J. Janicka, A. Sadiki, Proc. Combust. Inst. 30 (1) (2005) 537–547. [11] P. Auzillon, B. Fiorina, R. Vicquelin, N. Darabiha, O. Gicquel, D. Veynante, Proc. Combust. Inst. 33 (1) (2011) 1331–1338. [12] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [13] V.-K. Chakravarthy, S. Menon, Flow Turbul. Combust. 65 (2000) 133–161. [14] E. Knudsen, H. Pitsch, Combust. Flame 154 (4) (2008) 740–760. [15] V. Moureau, B. Fiorina, H. Pitsch, Combust. Flame 156 (4) (2009) 801–812. [16] J. Piana, D. Veynante, S. Candel, T. Poinsot, Direct numerical simulation analysis of the g-equation in premixed combustion, in: J. Chollet, P. Voke, L. Kleiser (Eds.), Direct and Large Eddy Simulation II, Kluwer Academic Publishers, 1997, pp. 321–330. [17] V. Moureau, P. Minot, H. Pitsch, C. Berat, J. Comput. Phys. 221 (2) (2007) 600– 614. [18] H. Pitsch, Annu. Rev. Fluid Mech. 38 (2006) 453–482. [19] B. Li, E. Baudoin, R. Yu, Z. Sun, Z. Li, X. Bai, M. Alden, M. Mansour, Proc. Combust. Inst. 32 (2) (2009) 1811–1818. [20] A. Roux, H. Pitsch, Large-eddy simulation of stratified combustion, in: Annual Research Briefs, Center For Turbulence Research, 2010, pp. 275–289. [21] B.A. Sen, S. Menon, Combust. Flame 157 (1) (2010) 62–74. [22] P. Domingo, L. Vervisch, D. Veynante, Combust. Flame 152 (3) (2008) 415–432. [23] B. Janus, A. Dreizler, J. Janicka, ASME Turbo Expo, Vienna (2004) 293–315. [24] B. Janus, A. Dreizler, J. Janicka, Flow Turbul. Combust. (75) (2004) 293–315. [25] B. Janus, A. Dreizler, J. Janicka, Proc. Combust. Inst. 31 (2) (2007) 3091–3098. [26] U. Maas, S. Pope, Combust. Flame 88 (1992) 239–264. [27] O. Gicquel, N. Darabiha, D. Thévenin, Proc. Combust. Inst. (2000) 1901–1908. [28] B. Fiorina, R. Baron, O. Gicquel, D. Thévenin, S. Carpentier, N. Darabiha, Combust. Theor. Modell. 7 (2003) 449–470. [29] J. van Oijen, L.P.H. de Goey, Combust. Theor. Modell. 8 (2004) 141–163. [30] B. Fiorina, O. Gicquel, L. Vervisch, S. Carpentier, N. Darabiha, Combust. Flame 140 (3) (2005) 147–160. [31] P. Nguyen, L. Vervisch, V. Subramanian, P. Domingo, Combust. Flame 157 (1) (2010) 43–61. [32] L. Vervisch, R. Haugel, P. Domingo, M. Rullaud, Flow Turbul. Combust. 5 (4) (2004) 1–36. [33] G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K. Hanson, S. Song, W.C. Gardiner, V.V. Lissianski, Z. Qin,