- Email: [email protected]
Numerical simulation of a turbulent partially premixed flame with inhomogeneous equivalence ratio
Fuel 124 (2014) 57–65
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Numerical simulation of a turbulent partially premixed flame with inhomogeneous equivalence ratio A. Benarous a,b,⇑, D. Karmed c, A. Liazid d, M. Champion c a
Department of Mechanical Engineering, University of Chlef, PO Box 151, Hay Essalem, Chlef, Algeria LCEMSM Laboratory, University of Chlef, PO Box 151, Chlef, Algeria c Department of Fluid, Thermal and Combustion Sciences, Pprime Institute (UPR3346), PO Box 30179, Futuroscope Chasseneuil Cedex, France d LTE Laboratory, National Polytechnic School, PO Box 1523, Oran, Algeria b
h i g h l i g h t s The model uses a two-scalar four delta presumed Pdf to describe mixing and reactions. The LW–P model is validated against the ICARE lean premixed flame. The model successfully predicts flow dynamics and turbulent flame characteristics.
a r t i c l e
i n f o
Article history: Received 9 March 2012 Received in revised form 20 January 2014 Accepted 21 January 2014 Available online 2 February 2014 Keywords: Partially premixed combustion LW–P model Scalar dissipation Bunsen burner
a b s t r a c t Reynolds averaged numerical simulation (RANS) of a turbulent partially premixed flame is performed using a presumed probability density function (Pdf) approach. The mixture fraction and the fuel mass fraction are respectively chosen to describe the local composition and the chemical reaction progress. The so called LW–P model was employed for the calculation of a lean methane–air flame stabilized by a stoichiometric pilot-flame. Numerical predictions of the jet-core as well as the mean quantities of the flame, namely the modified progress variable and the static temperature are compared successfully with the experiments. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction For the purpose of theoretical and numerical analysis, combustion systems are often divided into two groups: premixed systems where the fuel and oxidizer are perfectly mixed at the molecular level before being consumed by the flame and non-premixed or diffusion systems where the fuel and oxidiser enter the flame separately, mixing within the reaction zone. However, in practical situations relevant to working conditions in energy conversion devices, turbulent mixing of fuel and oxidizer prior to combustion leads to a reactive mixture that is not homogeneous. The development of systems operating in this mixed combustion mode has been driven by the need to reduce pollutant emissions, especially nitric oxides (NOx). By operating systems under fuel lean conditions, the reduction in temperature has a beneficial effect on the NOx level. However, as a system is run ⇑ Corresponding author at: Department of Mechanical Engineering, University of Chlef, PO Box 151, Hay Essalem, Chlef, Algeria. Tel.: +213 6 62 85 92 51; fax: +213 27 72 17 94. E-mail address: [email protected] (A. Benarous). http://dx.doi.org/10.1016/j.fuel.2014.01.069 0016-2361/Ó 2014 Elsevier Ltd. All rights reserved.
increasingly lean it becomes susceptible to several problems, including extinction, ignitability and flame stability [1]. By carefully designing and controlling spatial variation of mixture strength or areas of rich mixture, it becomes possible to achieve a stable lean combustion. Accordingly the equivalence ratio of the mixture varies in space and time, and combustion occurs in the form of partially premixed flames. Various experimental and numerical studies have been carried out to evaluate the influence of spatial equivalence ratio variations for different geometrical and initial conditions. The most noticeable effects that have been evidenced can be summarized as follows; (i) extension of the flammability limits [2,3], (ii) modification of the inner structure of the flame [4], and (iii) strong dependence of the combustion efficiency on both turbulence and scalar length scales [5,6]. Most of the models devoted to turbulent partially premixed combustion are based on approaches such as presumed Pdfs [7,8], flame surface density [9] and coherent flame modelling [10]. The most attractive among the proposed models is the presumed Pdf variety (using the joint Pdf of the mixture fraction
58
A. Benarous et al. / Fuel 124 (2014) 57–65
Nomenclature B c Cl D DT I k LT p P(Y, n) r R S ScT T uk, U X Y
pre-exponential factor progress variable turbulence model constant burner diameter turbulent diffusion coefficient turbulence intensity (%) turbulent kinetic energy turbulent length scale static pressure probability density function radial coordinate turbulent mixing model constant segregation factor turbulent schmidt number static temperature velocity component, mean axial velocity axial coordinate fuel mass fraction
Greek symbols dF flame brush thickness e turbulent kinetic energy dissipation Y, n, nY scalar dissipation rates n mixture fraction lT turbulent viscosity
and a reactive scalar) which offers a combination of simplicity, reduced computational cost and reasonable accuracy [8]. There are, however, several issues that arise including the shape of the Pdf, the statistical dependence of the variables and the closure of the scalar terms which appear in the transport equations of second moment quantities [11]. The Libby–Williams (LW) approach has already demonstrated its ability to recover not only the flamelet regime of turbulent combustion but also the thickened flame regime, at least for fully premixed situation [7]. The generalized form of the LW model introduced by Robin et al. [8] for partially premixed conditions is used in the present numerical study and will be denoted LW–P (Libby–Williams–Poitiers) in the following. In the latter model, the closure relies on a presumed joint scalar Pdf performed with a four Dirac delta functions. The application of the LW–P model to various situations involving either partially premixed or non premixed combustion has been successfully carried out by Ribert et al. [12] for the description of a reactive shear layer stabilized by a parallel incoming stream of hot products. Robin et al. [8] used the generalized form of the model to perform a RANS simulation on the ORACLES rig made of a combustion chamber fed with two parallel streams of premixed reactants, having different equivalence ratios. Recently, these authors validated the LW-P model on a stratified V-shaped flame of methane and air [13]. As stated above, the LW–P model has been validated on various partially premixed configurations witnessing slight gradients of equivalence ratio. In the present investigation turbulent partially premixed combustion is studied in the special case where a strong mean gradient of equivalence ratio exists. In this respect, a new closure expression for the scalar variances is used to interpolate between thin flamelets and thickened flame regimes. The flow configuration consists of a turbulent lean premixed burner-flame stabilized by an annular pilot-flame at stoichiometric conditions, as already investigated experimentally by Pavé [14] and Lachaux [15].
x q sij
instantaneous fuel consumption rate density viscous stress tensor element
Subscripts b burner 0 centreline value st stoichiometric value Superscripts g Reynolds average of g g 00 Favre fluctuations of g Favre average of g g~ Acronyms CFD Computational Fluid Dynamics ICARE Institut Combustion Aerothermique, Reactivite Environnement (Research Institute) LDA Laser Doppler Anemometry ORACLES A Combustion Test-Bench (ENSMA) NOx Nitrogen Oxides Pdf Probability Density Function RANS Reynolds Averaged Numerical Simulation
2. Mathematical modelling Favre averaging form for each extensive quantity (except pressure and density) is considered leading to the following mean equations of continuity, momentum and energy:
~k @q u @q þ ¼0 @t @xk u ~k @ q u ~k u ~j @ p @q @ 00 00 þ þ ¼ quk uj skj @xj @xk @t @xj ! ~ ~ ~ k hÞ h @ðq u @q @ @h @p 00 þ ¼ qD qu00k h þ @t @xk @xk @xk @t
ð1Þ
h denotes the sensible enthalpy of the mixture, given by:
h¼
Z 0 Y D h þ a F; a a
X
T
T0
cpa ðTÞdT
ð2Þ
0
where Ya, DhF;a are respectively the mass fraction and heat of formation of the species a. In the previous set of equations, body forces and radiation fluxes have been neglected. Turbulent fluxes of velocity and enthalpy are closed using viscosity and diffusivity hypothesis, respectively as:
qu00k u00j ¼ lT
~k @ u ~j 2 @ u ~i @u l @ h~ 00 qu00k h ¼ T þ dkj @xj @xk 3 @xi PrT @xk
where PrT = 0.7 is the turbulent Prandtl number, and bulent viscosity:
lT ¼ C l q k~2 =~e
ð3Þ
lT is the turð4Þ
The previous system is completed by the transport equations of ~ and its dissipation rate ~e [1]. the turbulent kinetic energy k In the LW–P model, all species mass fractions are related to the fuel mass fraction Y and the mixture fraction n defined as: max min n ¼ ðY max N2 Y N2 Þ=ðY N2 Y N2 Þ
ð5Þ
59
A. Benarous et al. / Fuel 124 (2014) 57–65
~1 ; P ~ 2 are the two conditional Pdf at n ¼ n1 and n ¼ n2 , defined as: P
~ 1 ðYÞ ¼ bdðY Y 11 Þ þ ð1 bÞdðY Y 12 Þ P ~ 2 ðYÞ ¼ cdðY Y 21 Þ þ ð1 cÞdðY Y 22 Þ P
ð7Þ
where a, b, c, n1, n2, Yii are the nine Pdf parameters depending on the location inside the flow field, as depicted in Fig. 2. Eq. (6) describing the Pdf allows for the evaluation of the mean related to the mean fuel mass fraction: chemical rate x
ZZ
¼q x
D
xðn; YÞ ~ Pðn; YÞdndY qðn; YÞ
ð8Þ
where D is the Pdf definition domain in the composition space. In the present model, we only consider the single-step global reaction, so that the rate of fuel consumption x is written as:
xðn; YÞ ¼ qðn; YÞ:BðnÞðY Y min Þ expðT a =TÞ
Fig. 1. Definition domain of (n, Y).
where Y N2 is the mass fraction of nitrogen and the superscripts max and min correspond to pure air and fuel, respectively. As shown in Fig. 1, the coordinate plane ðn; YÞ illustrates three different straight lines, the pure mixing line Y ¼ n, and the equilibrium lines for lean and rich mixtures Y ¼ 0 and Y = (n nst)/ (1 nst), respectively. The extreme situation where the mixture fraction n varies in the whole range 1 P n P 0 corresponds to a non premixed flame case. In the present situation, the definition of mixture fraction becomes: max n ¼ 1 Y N2 =Y max N2 where Y N2 is the mass fraction of nitrogen in air.
where Ta is the activation temperature, B the pre-exponential factor and Y min the minimum value of ; Ymin = 0 for lean mixtures and Ymin = (n nst)/(1 nst) for rich mixtures. The pre-exponential factor is assumed to be depending on the mixture fraction n [16]. It is selected such as the laminar burning velocity obtained with a detailed kinetic scheme is recovered [13]. 2.2. Scalar transport equations Regarding the mean reaction rate evaluation, five balance equations are needed to provide a closed system. These are the first and second order scalar quantities: Mean mixture fraction:
@ ~ @ @ @n nÞ þ ðq q u~ k ~n ¼ qD qu00k n00 @t @xk @xk @xk
ð10Þ
Variance of the mixture fraction:
2.1. Expression of the mean reaction rate As described by Robin et al. [8], one introduces a mass-weighted joint Pdf of n and Y made of four Dirac delta functions:
~ YÞ ¼ aP ~ 1 ðYÞdðn n1 Þ þ ð1 aÞP ~ 2 ðYÞdðn n2 Þ Pðn;
ð9Þ
@ 002 @ @ @n002 ~ k qn002 ¼ qn þ qD qu00k n002 u @t @xk @xk @xk
ð6Þ
2qD
!
@n00 @n00 @ ~n 2qu00k n00 @xk @xk @xk
ð11Þ
Mean fuel mass fraction:
@ ~ @ @ @Y q Y þ q u~ k Y~ ¼ qD qu00k Y 00 @t @xk @xk @xk
! þx
ð12Þ
Variance of the fuel mass fraction:
@ 002 @ @ @Y 002 ~ k qY 002 ¼ qY þ qD qu00k Y 002 u @t @xk @xk @xk 2qD
!
@Y 00 @Y 00 @ Y~ 2qu00k Y 00 @xk @xk @xk
þ 2Y 00 x
ð13Þ
Cross correlation:
@ 00 00 @ @ @Y 00 n00 ~ k qY 00 n00 ¼ qY n þ qD qu00k Y 00 n00 u @t @xk @xk @xk
!
@Y 00 @n00 @ ~n qu00k Y 00 @xk @xk @xk @ Y~ qu00k n00 þ n00 x @xk
2qD
Fig. 2. Locations of the four delta functions [8].
ð14Þ
60
A. Benarous et al. / Fuel 124 (2014) 57–65
In these equations, turbulent scalar fluxes are modelled by a gradient approximation:
The reactive contribution to the mean dissipation rate is given by:
@ ~n @ Y~ DT ; qu00k Y 00 ¼ q @xk @xk ~n002 @ @ Y~ 002 002 002 DT qu00k n ¼ q DT ; qu00k Y ¼ q @xk @xk 00 ~00 n @Y 00 00 qu00k Y n ¼ q DT @xk
~ ~ AÞ ð1 þ AÞð1 ~n002 x 2Y 00 x 2Y Y~ max Y~ þ Y~ min x 2q Y~ max Y~ min ! @ ~n ~ ~ v ~ n þ qu00k n00 þ ð1 hÞAÞ ð1 þ AÞðh @xk ! ~ ~ ~ qu00 Y 00 @ n þ qu00 n00 @ Y þ A ~v ~n þ ð1 þ AÞ k k @xk @xk
qu00k n00 ¼ q DT
ð15Þ
2
where DT = Clk /ScTe and C l ¼ 0:09. k and e are respectively the turbulent kinetic energy and its dissipation rate, and ScT is a turbulent Schmidt number set to its usual value 0.7.
The scalar dissipation rates which appear in the transport Eqs. (11), (13), and (14) represent the mean mixing rate at small scales. In the case of partially premixed combustion, it has been shown that the reaction rate depends on three dissipation rates [17]; ~ n ¼ 2q2n ¼ 2qD @n the mixture fraction dissipation rate v @x ~ Y ¼ 2q2Y ¼ 2qD @Y fuel mass fraction dissipation rate v @x ~ nY ¼ 2q2nY ¼ 2qD @n the cross-dissipation rate v @x
00 k
@Y 00 . @xk
00 k
00 k
@n00 , @xk @Y 00 @xk
and
The modelling
~e v~ n ¼ R ~ qn002 k
ð16Þ
where the constant R is the timescale ratio of the turbulence to scalar fluctuation. Its current is about 2 [18]. This model is efficient for passive scalars. However, the dissipation rate of the reactive scalar fluctuations is influenced by turbulence, chemical reaction and their interactions. As the classical model accounts only for turbulence it will not describe correctly the reactive scalar dissipation rate. The model discussed in [13,19] was derived by considering the behaviour of reactive scalar fluctuations in the limit of thickened flamelets and infinitely thin flamelets. A segregation factor S is employed to interpolate between these two regimes such as:
~e v~ Y ¼ 2q S2Y þ ð1 SÞR ~ Y~ 002 k
~ 1A @ ~n ~ v ~n002 x ~ n þ 2qu00k n00 þ ðh þ ð1 hÞAÞ @xk Y~ max Y~ min
qu00k Y 00
@ ~n @ Y~ qu00k n00 @xk @xk
The stainless steel cylindrical high pressure combustion chamber of ICARE can sustain premixed flames up to 1 MPa [14]. The chamber consists of two cylindrical halves of 600 mm high (Fig. 3) and 300 mm inner diameter. Each halve is equipped with four windows of 100 mm diameter designed for optical diagnostics. The burner with a 25 mm diameter is placed in the centre of the chamber and can be moved along the vertical axis by means of a stepping motor. The upper part of the burner has a convergent profile. A porous grid is placed in the bottom part of this burner section to ensure chamber-inlet turbulence. The main flame consists of a lean methane–air mixture with an equivalent ratio of 0.6. It is surrounded by a 2 mm pilot flame at stoichiometric conditions. 4. Numerical results and discussion The LW–P model described in the previous sections has been implemented in the computational fluid dynamic software Code_Saturne developed by EDF, Archambeau et al. [20]. It is a parallel code based on a finite volume method. The set of the considered equations consists of the Favre averaged Navier–Stokes equations
ð17Þ ð18Þ
where S is defined as:
S¼
qY 002 ðqY 002 Þmin ðqY 002 Þmax ðqY 002 Þmin
ð19Þ
and
qY 002
max
~ 2 ~n002 Y~ max Y~ Y~ Y~ min þ h þ ð1 hÞA ¼q
ð20Þ
and
~ 2 Þ~n002 ðqY 002 Þmin ¼ ðh þ ð1 hÞA
ð21Þ
~ n and h is defined as h ¼ ðY ~ Y ~ min Þ= ðh þ ð1 hÞAÞ where qY n ¼ q ~ max Y ~ min Þ and Y ~ max and Y ~ min are defined by the bounds of the ðY ~ is the mean slope of the equilibrium line such composition space. A 00 00
ð23Þ
3. The ICARE rig
~e v~ nY ¼ 2q S2nY þ ð1 SÞR ~ Y 00 n00 k
!
the
of these three dissipation rates is one of the most challenging problems in the simulation of partially premixed systems. The classical model describing the dissipation rate of the mix~ n is based on the analogy with the turture fraction fluctuations v bulent mixing [18]:
and
2nY ¼ n00 x 2q
2.3. Scalar dissipation terms
ð22Þ
~002
~ ¼ 0 for lean mixture and A ~ ¼ 1=ð1 n Þ for rich mixture, that A st where nst stands for the stoichiometric mixture fraction.
Fig. 3. ICARE test-rig view [14].
61
A. Benarous et al. / Fuel 124 (2014) 57–65
for incompressible flows completed with sub-models for turbulence and equations describing the additional scalars (enthalpy, species concentration, and possible user-defined scalars). The operating scheme is based on a prediction of the velocity field followed by a pressure correction step. Equations for turbulence and scalars are solved separately afterward. The space discretization is based on the fully conservative unstructured finite-volume framework, with a fully collocated arrangement for all variables. In the present work, turbulent mixing is represented by the two-equations (k–e) model. An unstructured grid of approximately 76,000 cells has been generated to represent half of the physical domain. The retained mesh has been obtained from several numerical tests on finer grids which did not influence the centreline static temperature field by more than 2% (Fig. 4a).
4.1. Boundary conditions
1800
1500
1200 Mesh 1 : 29378 Cells Mesh 2: 49348 Cells
ð24Þ
has been imposed. In Eq. (24), rb denotes the burner radius and u0 is the centreline velocity. For such ambient pressure case, the velocity exponent is taken to be n ¼ 0:2818 [15]. According to LDA measurements of the turbulent kinetic energy k [15], a polynomial expression was imposed at the main inlet flow, with a minimum value of k0 ¼ 0:047 m2 =s2 (Fig. 5b). The mean dissipation rate ~e has been chosen in such a way that the measured grid turbulence decay is recovered. 1.00
0.80
0.60
Velocity Profile at burner exit
0.40
Experimental Measurements at X/D = 0.2 [14]
0.20
Mesh 3: 76321 Cells
900
0.00 0.00 600
0.10
0.20
0.30
0.40
0.50
DimensionlessRadialDistance(R/D) Fig. 5a. Radial velocity profile at X/D = 0.2.
300
0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D)
0.10
Fig. 4a. Mesh study based on the centreline temperature.
Turbulent Kinetic Energy (m**2/s**2)
Centreline Static Temperature (K)
uðrÞ ¼ u0 ð1 r=r b Þn
Dimensionless Radial velocity (U/U0)
Inlet boundary conditions must be handled with special care, since they may have a strong influence on combustion inside the computational domain. This is often complicated by several constraints: (i) the measurements of all variables at the exact location of the numerical inlet boundary are not always available; (ii) also, an additional constraint results from the difficulty to measure
some transported variables related to the turbulent combustion model. This is especially true for the turbulence mean dissipation rate ~e. Boundary conditions on the left side of the computational domain are given by the flow symmetry (Fig. 4b). The right and top sides are approximated by outlet conditions, whereas an adiabatic wall condition is specified on the burner thickness (burner lips). To preserve a fixed bulk velocity ud ¼ 2:3 m/s [14] within the burner (Fig. 5a), a radial velocity profile:
0.09
TKE Profile at burner exit Experimental Measurements at X/D = 0.2 [14]
0.08
0.07
0.06
0.05
0.04 0.00
0.10
0.20
0.30
0.40
0.50
DimensionlessRadialDistance (R/D) Fig. 4b. The computational domain.
Fig. 5b. Radial profile of the turbulent kinetic energy at X/D = 0.2.
A. Benarous et al. / Fuel 124 (2014) 57–65
In order to recover the experimental values of k and u at the first measurement station X/D = 0.2, we have adjusted the numerical boundary conditions at the domain inlet. The pilot flame is modelled as a co-flow of hot flue gases being burned at stoichiometric conditions. Accordingly, a uniform and low value of the axial velocity was imposed at the second inlet whereas turbulent quantities were expressed in terms of intensity and hydraulic diameter. The numerical values for the boundary conditions are depicted in Table 1.
4.2. The nonreactive flow As a first step, it is interesting to perform a non-reactive flow computation. The so-called ‘‘cold flow simulation’’ helps to validate the appropriate turbulence model describing the mixing and tune the ‘‘numerical’’ boundary conditions. Moreover, solving the scalar transport equations without chemical source terms provides a ‘‘guess’’ solution field for the reactive calculations. Figs. 6 and 7 show respectively the centreline profiles of the mean velocity com~ along the axial coordi~ and the turbulent kinetic energy k ponent u nate. A quasi-linear decrease is noticed on the velocity profile until X/D = 3.2 where a slope change occurs. As explained by Hinze [21], this behaviour is due to the training effect characterised by an important momentum transfer from the jet to its surrounding. As the jet develops axially, the mixing zone thickens and eventually reaches the axis while affecting the centreline velocity. Good agreement between numerical and experimental data is observed for the dynamic field (Fig. 6). The predicted jet-core size (X/ D 3.2) is quite close to the measured value (X/D = 3.0) [15]. As stated above, turbulence decay along the chamber axis (Fig. 7) was obtained as a polynomial distribution for the mean dissipation
Table 1 Boundary conditions used in the simulation.
Velocity (m/s) k (m2/s2) e (m2/s3) T (K) Equivalence ratio
Main flow
Co-flow
u0 = 2.56 r2 polynomial r3 polynomial 300 0.6
0.05 ud k = k (I = 5) e = e (I = 5%) T (nst) 1.0
Dimensionless Axial Velocity (U/U0)
1.00
0.95
0.90
0.85 Cold Flow Calculations Experimental Measurements [15]
0.80
0.75
0.0
1.0
2.0
3.0
4.0
5.0
Dimensionless Axial Distance (X / D) Fig. 6. Centreline profile of the axial velocity.
6.0
0.13 0.12 Experimental Measurements [15]
Turbulent Kinetic Energy (m**2/s**2)
62
0.11
Cold Flow Calculation
0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) Fig. 7. Axial profile of the turbulent kinetic energy.
rate ~e at the domain inlet. This was done by imposing a constant turbulent length scale LT ¼ 2:9 mm which is quite close to the turbulence grid size [14]. Predicted values of the mean turbulent kinetic energy at various axial stations (Fig. 7) seem to fit well the ~ is reached at experimental measurements. Minimum value for k the potential core region, exactly at the axial location X=D 2:48 where its exhibits an underestimation of 17% compared with the measurements. This discrepancy is certainly due to the turbulent stresses isotropy considered for the (k–e) model formulation, which underestimates the jet spreading in both axial and radial directions. Owing to the important chamber-to-burner diameter ratio (rch/ rb = 12), the wall effects seem to be unimportant which further justifies the free character of the jet within the chamber.
4.3. Partially premixed combustion The more general case of partially premixed combustion is now considered. The study is performed in the case of the stratified condition with an equivalence ratio varying from unit value for the hot co-flow to 0.6 at the main flow entrance. In order to allow comparison with the reactive measurements, we define a pseudo-progress ~ ~ variable ~c ¼ 1 Y= n which accounts for both fuel consumption and turbulent mixing [22]. The previous expression is only valid for lean mixture; therefore the corresponding variable is zero in the fresh mixture and takes a unit value in the burnt gases. Regarding the dynamic variables in the presence of the flame, ~ and the turcentreline profiles of the mean velocity component u ~ are depicted in Figs. 8a and 8b respectively. bulent kinetic energy k It is noticed that the slope of the axial velocity decay upstream the flame front is quite similar to that of the non reactive case. Momentum transfer still exists in this region. However, the isobaric character of the flame yields an acceleration of the flow before the potential core region (X/D 1.5). Globally, the calculated profile of the axial velocity is in good agreement with the experimental results. The acceleration behind the flame induced by the thermal expansion is well predicted (Fig. 8a). However, the profile of the numerical turbulent kinetic energy is less accurate. Our calculation shows a relative underestimation up to 26% downstream the flame brush, compared with experimental measurements (Fig. 8b).
63
A. Benarous et al. / Fuel 124 (2014) 57–65 1.0
0.8
1.04
0.6 1.00
Flame Calculations
0.4
Progress
Dimensionless Axial Velocity (U/U0)
1.08
Experimental Measurements [15]
0.96
Pseudo - Progress variable
0.2
Fig. 9a. Spatial contour of the static temperature. 0.92
0.0
1.0
2.0
3.0
4.0
5.0
6.0
0.0
2400
Dimensionless Axial Distance (X/D)
Calculations at r/D = 0.56
Fig. 8a. Centreline velocity distribution.
2100
Calculations at r/D = 0.5
1800
Experimental Measurements [15]
1.0
0.05
0.8
0.04
Flame Calculations Pseudo-Progress variable
0.6
Experimental Measurements [15]
0.4
Progress
Turbulent Kinetic Energy (m**2/s**2)
Static Temperature (K)
Centreline Calculations
1500
1200
900
600
300 0.0
0.03
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) 0.2
0.02 0.0
1.0
2.0
3.0
4.0
5.0
0.0 6.0
Dimensionless Axial Distance (X/D) Fig. 8b. Turbulent kinetic energy distribution along the chamber axis.
~ in the flame region is certainly due The calculated decrease of k to the use of the first order (k–e) model in the closure of the scalar turbulent fluxes (Eq. (15)). Indeed, the model do not account for the flame generated turbulence resulting from the counter-gradient mechanism [23]. Spatial contour of the static temperature seems to have a conical shape (Fig. 9a). It is also noticed that in the present case, the coflow entrance is the hottest zone of the combustion chamber. This was expected since the co-flow mixture was burned at stoichiometric conditions, yielding an adiabatic value for the static temperature (2318 K). However, the co-flow does not seem to affect the thermal level of the flue gases throughout the chamber. This is a consequence of the low flow rate and the laminar character imposed at the co-flow inlet. Within the chamber, the flue gases are still at a much lower temperature (1645 K) than that of the co-flow entrance. Fig. 9b shows the static temperature profiles along the axial distance of the chamber at three radial locations, the centreline (r/
Fig. 9b. Static temperature profiles for various radial locations.
D = 0), the burner lips (r/D = 0.5) and the mid-entrance of the coflow (r/D = 0.56). Along the axial location X ¼ 0, the temperature profile moves from the highest values (2318 K), characteristic of the hot co-flow, to a value T~ ¼ 380 K induced by the slight thermal diffusion, and then finally to the lowest value (300 K) which refers to the temperature of the fresh gases. However, the tendency is different downstream, where the three profiles asymptotically reach a limiting value (1645 K) corresponding to the temperature of the products. The axial distance where the flue gas temperature is reached (Fig. 9b) decreases with the radial location which confirms the conical shape of the flame. Compared with the chamber centreline measurements, the static temperature exhibits a global overestimation of 6%. This discrepancy is probably due to the global reaction modelling adopted for the chemistry description. Indeed, the model considers only carbon dioxide CO2 and water vapour H2O as final products in the methane–air reaction and does not account for enthalpy losses by dissociation. The limiting values for ~c are clearly shown in the spatial contour depicted in Fig. 10a. The turbulent flame seems to be stabilized over the burner lips and exhibits a typical conical shape. Axial distribution of the progress variable agrees the measurements (Fig. 10b), although the experimental curve is slightly steeper than
64
A. Benarous et al. / Fuel 124 (2014) 57–65 0.6
0.4
LW-P Calculations Linear Fitting Y = (-1/a) X + b
Z / Flame Brush
0.2
0.0
-0.2
-0.4
Fig. 10a. Spatial contour of the mean progress variable. -0.6 -2.50 -2.00 -1.50 -1.00 -0.50 0.00
1.00
0.50
1.00
1.50
2.00
2.50
Ln((1-)/) Fig. 11. Progress variable across the flame brush.
Pseudo Progress Variable
0.80
5. Conclusion 0.60
0.40
Experimental Measurements [15] Numerical (LWP) Calculations
0.20
0.00 0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) Fig. 10b. Progress variable distribution along the chamber axis.
the numerical one. Accordingly, the maximum thickness of the conical flame dF can be evaluated by means of the centreline distribution of the progress variable. It is defined as the axial distance for which ~c increases from 0.01 to 0.99 [15]. Hence, the calculation shows an apparent flame thickness dF ¼ 62:5 mm slightly greater (11%) than the measurements [14]. As it was found by Robin et al. [22], the mean flame is most probably thicker than the experimental flame. This discrepancy is due to the gradient assumption considered on the scalar turbulent fluxes qu00 Y 00 . As already stated by Lipatnikov et al. [24], the axial progress variable profiles may be represented by a general form when using a reduced spatial coordinate across the flame brush. If we note this quantity Z ¼ X X ~c¼0:5 , it is possible to recover an exponential tendency of the progress variable along the flame thickness. Indeed, a linear evolution of Z=dF with the quantity Lnðð1 ~cÞ=~cÞ is clearly visible in Fig. 11. In order to obtain an analytic expression, we perform a linear fitting using the least square method. The corresponding slope a ¼ 4:20 is found to agree the experimental value as pointed out by Pavé (a = 4.28) [14] for a Bunsen premixed flame and Shepherd (a = 4.02) [25] for a V-shaped flame. Hence, it appears that the evolution of ~c in the flame brush may be represented by a single curve which seems to be fairly universal.
We have studied a partially premixed pilot-stabilized methaneair flame. The LW–P model based on a two-scalar ðY; nÞ four delta Pdf was used to describe mixing and reaction progress within the burner-chamber device. As it was stated in the experiments, the hot co-flow has allowed a stabilization of the main flame without affecting the global thermal level within the chamber. Moreover, the computed static temperature overestimates the measured one by only 6%, despite the use of a one-step chemistry. Spatial contour of the progress variable has clearly shown the conical shape of the main flame. Owing to the bridging expression used for the cross dissipation terms, values of ~c agree quite well with the experiments. In the cold flow situation, a good qualitative agreement with the experimental data is found, especially for the jet-core size and the turbulence decay along the chamber axis. Moreover, the mean velocity profile downstream the flame brush fits well the experiments despite a slight underestimation in the fresh gas region. The discrepancy is due to the fact that the adopted numerical conditions at the burner exit were issued from the cold flow case. However, it is shown that the present version of the LW–P model does not improve the prediction of turbulent kinetic energy downstream of the flame brush. The numerical depletion observed on the mean turbulent kinetic energy is probably due to the use of a first order model (k–e), a model that does not take into account the combustion induced turbulence. Some efforts are still necessary to improve the turbulent mixing by a second-order model. Such a model will allow us to deal with the flame-generated turbulence and counter-gradient diffusion in partially premixed conditions. The extension of the present model to include detailed chemistry effects still remains a challenging task. References [1] Peters N. Turbulent Combustion. 1st ed. Cambridge: Cambridge University Press; 2000. [2] Ra Y, Cheng WK. Laminar flame propagation through a step-stratified charge. In: Proceedings of the Fifth Int. Symp on Diagnostics and Modeling of Combustion in Internal Engines, 2001, p. 251–57. [3] Marzouk YM, Ghoniem AF, Najm HN. Dynamic response of strained premixed flames to equivalence ratio gradients. Proc Combust Inst 2000;28:1859–66.
A. Benarous et al. / Fuel 124 (2014) 57–65 [4] Mura A, Galzin F, Borghi R. A unified PDF-flamelet model for turbulent premixed combustion. Combust Sci Technol 2003;175(7):1–37. [5] Renou B, Boukhalfa A, Puechberty D, Trinité M. Local scalar flame properties of freely propagating premixed turbulent flames at various Lewis numbers. Combust Flame 2000;123:l507–1521. [6] Jimenez C, Cuenot B, Poinsot T, Haworth DC. Numerical simulation and modeling for lean stratified propane–air flames. Combust Flame 2002;128:1–21. [7] Bray KNC, Champion M, Libby PA. Bray-Moss thermochemistry revisited. Combust Sci Technol 2002;174(7):167–74. [8] Robin V, Mura A, Champion M, Plion P. A multi-Dirac presumed PDF model for turbulent reactive flows with variable equivalence ratio. Combust Sci Technol 2006;178:1843–70. [9] Bondi S, Jones WP. A combustion model for premixed flames with varying stoichiometry. Proc Combust Inst 2002;29:2123–9. [10] Helie J, Trouve A. A modified coherent flame model to describe turbulent flame propagation in mixtures with variable composition. Proc Combust Inst 2000;28:193–201. [11] Kolla H, Rogerson J, Chakraborty N, Swaminathan N. Scalar dissipation rate modeling and its validation. Combust Sci Technol 2009;181:518–35. [12] Ribert G, Champion M, Plion P. Modeling turbulent reactive flows with variable equivalence ratio: application to the calculation of a reactive shear layer. Combust Sci Technol 2004;176:907–23. [13] Robin V, Mura A, Champion M, Degardin. O, Renou B, Boukhalfa M. Experimental and numerical analysis of stratified turbulent V-shaped flames. Combust Flame 2008;153:288–315. [14] Pavé D. Contribution to study the structure of a lean premixed turbulent methane–air flame. PhD thesis, Orléans University, France, 2002.
65
[15] Lachaux T. Pressure effects on the structure and dynamics of lean premixed turbulent methane–air flames. PhD thesis, Orléans University, France, 2004. [16] Fernandez-Tarrazo E, Sanchez A, Linan A, Williams FA. A simple one-step chemistry model for partially premixed hydrocarbon combustion. Combust Flame 2006;147(1–2):32–8. [17] Bray K, Domingo P, Vervich L. The role of progress variable in models for partially premixed turbulent combustion. Combust Flame 2005;141:431–7. [18] Spalding DB. Mixing and chemical reaction in confined turbulent flames. Chem Eng Sci 1971;26:95–107. [19] Mura A, Robin V, Champion M. A second-order model for turbulent reactive flows with variable equivalence ratio. Combust Flame 2007;149:217–24. [20] Archambeau F, Mehitoua N, Sakiz M. A finite volume code for the computation of turbulent flows. Int J Finite, vols. 1–62. 2004. [21] Hinze JO. Turbulence. 2nd ed. New York: McGraw Hill Series in Mechanical Engineering; 1975. [22] Robin V, Guilbert N, Mura A, Champion M. Turbulent combustion modeling of mixtures with inhomogeneous equivalence ratio: application to a flame stabilized by the sudden expansion of a 2-D channel. Comptes Rendus de Mécanique 2010;338:40–7. [23] Robin V, Mura N, Champion M. Algebraic models for turbulent transports in premixed flames. Combust Sci Technol 2012;184:1718–42. [24] Lipatnikov AN, Chomiak J. Turbulent flame speed and thickness: phenomenology, evaluation and application in multi-dimensional simulations. Prog Energ Com Sci 2002;28:1–73. [25] Shepherd I.G. Flame surface density and burning rate in premixed turbulent flames. In: 26th Symposium (International) on combustion; The Combustion Institute, 1996. p. 373–79.
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Numerical simulation of a turbulent partially premixed flame with inhomogeneous equivalence ratio A. Benarous a,b,⇑, D. Karmed c, A. Liazid d, M. Champion c a
Department of Mechanical Engineering, University of Chlef, PO Box 151, Hay Essalem, Chlef, Algeria LCEMSM Laboratory, University of Chlef, PO Box 151, Chlef, Algeria c Department of Fluid, Thermal and Combustion Sciences, Pprime Institute (UPR3346), PO Box 30179, Futuroscope Chasseneuil Cedex, France d LTE Laboratory, National Polytechnic School, PO Box 1523, Oran, Algeria b
h i g h l i g h t s The model uses a two-scalar four delta presumed Pdf to describe mixing and reactions. The LW–P model is validated against the ICARE lean premixed flame. The model successfully predicts flow dynamics and turbulent flame characteristics.
a r t i c l e
i n f o
Article history: Received 9 March 2012 Received in revised form 20 January 2014 Accepted 21 January 2014 Available online 2 February 2014 Keywords: Partially premixed combustion LW–P model Scalar dissipation Bunsen burner
a b s t r a c t Reynolds averaged numerical simulation (RANS) of a turbulent partially premixed flame is performed using a presumed probability density function (Pdf) approach. The mixture fraction and the fuel mass fraction are respectively chosen to describe the local composition and the chemical reaction progress. The so called LW–P model was employed for the calculation of a lean methane–air flame stabilized by a stoichiometric pilot-flame. Numerical predictions of the jet-core as well as the mean quantities of the flame, namely the modified progress variable and the static temperature are compared successfully with the experiments. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction For the purpose of theoretical and numerical analysis, combustion systems are often divided into two groups: premixed systems where the fuel and oxidizer are perfectly mixed at the molecular level before being consumed by the flame and non-premixed or diffusion systems where the fuel and oxidiser enter the flame separately, mixing within the reaction zone. However, in practical situations relevant to working conditions in energy conversion devices, turbulent mixing of fuel and oxidizer prior to combustion leads to a reactive mixture that is not homogeneous. The development of systems operating in this mixed combustion mode has been driven by the need to reduce pollutant emissions, especially nitric oxides (NOx). By operating systems under fuel lean conditions, the reduction in temperature has a beneficial effect on the NOx level. However, as a system is run ⇑ Corresponding author at: Department of Mechanical Engineering, University of Chlef, PO Box 151, Hay Essalem, Chlef, Algeria. Tel.: +213 6 62 85 92 51; fax: +213 27 72 17 94. E-mail address: [email protected] (A. Benarous). http://dx.doi.org/10.1016/j.fuel.2014.01.069 0016-2361/Ó 2014 Elsevier Ltd. All rights reserved.
increasingly lean it becomes susceptible to several problems, including extinction, ignitability and flame stability [1]. By carefully designing and controlling spatial variation of mixture strength or areas of rich mixture, it becomes possible to achieve a stable lean combustion. Accordingly the equivalence ratio of the mixture varies in space and time, and combustion occurs in the form of partially premixed flames. Various experimental and numerical studies have been carried out to evaluate the influence of spatial equivalence ratio variations for different geometrical and initial conditions. The most noticeable effects that have been evidenced can be summarized as follows; (i) extension of the flammability limits [2,3], (ii) modification of the inner structure of the flame [4], and (iii) strong dependence of the combustion efficiency on both turbulence and scalar length scales [5,6]. Most of the models devoted to turbulent partially premixed combustion are based on approaches such as presumed Pdfs [7,8], flame surface density [9] and coherent flame modelling [10]. The most attractive among the proposed models is the presumed Pdf variety (using the joint Pdf of the mixture fraction
58
A. Benarous et al. / Fuel 124 (2014) 57–65
Nomenclature B c Cl D DT I k LT p P(Y, n) r R S ScT T uk, U X Y
pre-exponential factor progress variable turbulence model constant burner diameter turbulent diffusion coefficient turbulence intensity (%) turbulent kinetic energy turbulent length scale static pressure probability density function radial coordinate turbulent mixing model constant segregation factor turbulent schmidt number static temperature velocity component, mean axial velocity axial coordinate fuel mass fraction
Greek symbols dF flame brush thickness e turbulent kinetic energy dissipation Y, n, nY scalar dissipation rates n mixture fraction lT turbulent viscosity
and a reactive scalar) which offers a combination of simplicity, reduced computational cost and reasonable accuracy [8]. There are, however, several issues that arise including the shape of the Pdf, the statistical dependence of the variables and the closure of the scalar terms which appear in the transport equations of second moment quantities [11]. The Libby–Williams (LW) approach has already demonstrated its ability to recover not only the flamelet regime of turbulent combustion but also the thickened flame regime, at least for fully premixed situation [7]. The generalized form of the LW model introduced by Robin et al. [8] for partially premixed conditions is used in the present numerical study and will be denoted LW–P (Libby–Williams–Poitiers) in the following. In the latter model, the closure relies on a presumed joint scalar Pdf performed with a four Dirac delta functions. The application of the LW–P model to various situations involving either partially premixed or non premixed combustion has been successfully carried out by Ribert et al. [12] for the description of a reactive shear layer stabilized by a parallel incoming stream of hot products. Robin et al. [8] used the generalized form of the model to perform a RANS simulation on the ORACLES rig made of a combustion chamber fed with two parallel streams of premixed reactants, having different equivalence ratios. Recently, these authors validated the LW-P model on a stratified V-shaped flame of methane and air [13]. As stated above, the LW–P model has been validated on various partially premixed configurations witnessing slight gradients of equivalence ratio. In the present investigation turbulent partially premixed combustion is studied in the special case where a strong mean gradient of equivalence ratio exists. In this respect, a new closure expression for the scalar variances is used to interpolate between thin flamelets and thickened flame regimes. The flow configuration consists of a turbulent lean premixed burner-flame stabilized by an annular pilot-flame at stoichiometric conditions, as already investigated experimentally by Pavé [14] and Lachaux [15].
x q sij
instantaneous fuel consumption rate density viscous stress tensor element
Subscripts b burner 0 centreline value st stoichiometric value Superscripts g Reynolds average of g g 00 Favre fluctuations of g Favre average of g g~ Acronyms CFD Computational Fluid Dynamics ICARE Institut Combustion Aerothermique, Reactivite Environnement (Research Institute) LDA Laser Doppler Anemometry ORACLES A Combustion Test-Bench (ENSMA) NOx Nitrogen Oxides Pdf Probability Density Function RANS Reynolds Averaged Numerical Simulation
2. Mathematical modelling Favre averaging form for each extensive quantity (except pressure and density) is considered leading to the following mean equations of continuity, momentum and energy:
~k @q u @q þ ¼0 @t @xk u ~k @ q u ~k u ~j @ p @q @ 00 00 þ þ ¼ quk uj skj @xj @xk @t @xj ! ~ ~ ~ k hÞ h @ðq u @q @ @h @p 00 þ ¼ qD qu00k h þ @t @xk @xk @xk @t
ð1Þ
h denotes the sensible enthalpy of the mixture, given by:
h¼
Z 0 Y D h þ a F; a a
X
T
T0
cpa ðTÞdT
ð2Þ
0
where Ya, DhF;a are respectively the mass fraction and heat of formation of the species a. In the previous set of equations, body forces and radiation fluxes have been neglected. Turbulent fluxes of velocity and enthalpy are closed using viscosity and diffusivity hypothesis, respectively as:
qu00k u00j ¼ lT
~k @ u ~j 2 @ u ~i @u l @ h~ 00 qu00k h ¼ T þ dkj @xj @xk 3 @xi PrT @xk
where PrT = 0.7 is the turbulent Prandtl number, and bulent viscosity:
lT ¼ C l q k~2 =~e
ð3Þ
lT is the turð4Þ
The previous system is completed by the transport equations of ~ and its dissipation rate ~e [1]. the turbulent kinetic energy k In the LW–P model, all species mass fractions are related to the fuel mass fraction Y and the mixture fraction n defined as: max min n ¼ ðY max N2 Y N2 Þ=ðY N2 Y N2 Þ
ð5Þ
59
A. Benarous et al. / Fuel 124 (2014) 57–65
~1 ; P ~ 2 are the two conditional Pdf at n ¼ n1 and n ¼ n2 , defined as: P
~ 1 ðYÞ ¼ bdðY Y 11 Þ þ ð1 bÞdðY Y 12 Þ P ~ 2 ðYÞ ¼ cdðY Y 21 Þ þ ð1 cÞdðY Y 22 Þ P
ð7Þ
where a, b, c, n1, n2, Yii are the nine Pdf parameters depending on the location inside the flow field, as depicted in Fig. 2. Eq. (6) describing the Pdf allows for the evaluation of the mean related to the mean fuel mass fraction: chemical rate x
ZZ
¼q x
D
xðn; YÞ ~ Pðn; YÞdndY qðn; YÞ
ð8Þ
where D is the Pdf definition domain in the composition space. In the present model, we only consider the single-step global reaction, so that the rate of fuel consumption x is written as:
xðn; YÞ ¼ qðn; YÞ:BðnÞðY Y min Þ expðT a =TÞ
Fig. 1. Definition domain of (n, Y).
where Y N2 is the mass fraction of nitrogen and the superscripts max and min correspond to pure air and fuel, respectively. As shown in Fig. 1, the coordinate plane ðn; YÞ illustrates three different straight lines, the pure mixing line Y ¼ n, and the equilibrium lines for lean and rich mixtures Y ¼ 0 and Y = (n nst)/ (1 nst), respectively. The extreme situation where the mixture fraction n varies in the whole range 1 P n P 0 corresponds to a non premixed flame case. In the present situation, the definition of mixture fraction becomes: max n ¼ 1 Y N2 =Y max N2 where Y N2 is the mass fraction of nitrogen in air.
where Ta is the activation temperature, B the pre-exponential factor and Y min the minimum value of ; Ymin = 0 for lean mixtures and Ymin = (n nst)/(1 nst) for rich mixtures. The pre-exponential factor is assumed to be depending on the mixture fraction n [16]. It is selected such as the laminar burning velocity obtained with a detailed kinetic scheme is recovered [13]. 2.2. Scalar transport equations Regarding the mean reaction rate evaluation, five balance equations are needed to provide a closed system. These are the first and second order scalar quantities: Mean mixture fraction:
@ ~ @ @ @n nÞ þ ðq q u~ k ~n ¼ qD qu00k n00 @t @xk @xk @xk
ð10Þ
Variance of the mixture fraction:
2.1. Expression of the mean reaction rate As described by Robin et al. [8], one introduces a mass-weighted joint Pdf of n and Y made of four Dirac delta functions:
~ YÞ ¼ aP ~ 1 ðYÞdðn n1 Þ þ ð1 aÞP ~ 2 ðYÞdðn n2 Þ Pðn;
ð9Þ
@ 002 @ @ @n002 ~ k qn002 ¼ qn þ qD qu00k n002 u @t @xk @xk @xk
ð6Þ
2qD
!
@n00 @n00 @ ~n 2qu00k n00 @xk @xk @xk
ð11Þ
Mean fuel mass fraction:
@ ~ @ @ @Y q Y þ q u~ k Y~ ¼ qD qu00k Y 00 @t @xk @xk @xk
! þx
ð12Þ
Variance of the fuel mass fraction:
@ 002 @ @ @Y 002 ~ k qY 002 ¼ qY þ qD qu00k Y 002 u @t @xk @xk @xk 2qD
!
@Y 00 @Y 00 @ Y~ 2qu00k Y 00 @xk @xk @xk
þ 2Y 00 x
ð13Þ
Cross correlation:
@ 00 00 @ @ @Y 00 n00 ~ k qY 00 n00 ¼ qY n þ qD qu00k Y 00 n00 u @t @xk @xk @xk
!
@Y 00 @n00 @ ~n qu00k Y 00 @xk @xk @xk @ Y~ qu00k n00 þ n00 x @xk
2qD
Fig. 2. Locations of the four delta functions [8].
ð14Þ
60
A. Benarous et al. / Fuel 124 (2014) 57–65
In these equations, turbulent scalar fluxes are modelled by a gradient approximation:
The reactive contribution to the mean dissipation rate is given by:
@ ~n @ Y~ DT ; qu00k Y 00 ¼ q @xk @xk ~n002 @ @ Y~ 002 002 002 DT qu00k n ¼ q DT ; qu00k Y ¼ q @xk @xk 00 ~00 n @Y 00 00 qu00k Y n ¼ q DT @xk
~ ~ AÞ ð1 þ AÞð1 ~n002 x 2Y 00 x 2Y Y~ max Y~ þ Y~ min x 2q Y~ max Y~ min ! @ ~n ~ ~ v ~ n þ qu00k n00 þ ð1 hÞAÞ ð1 þ AÞðh @xk ! ~ ~ ~ qu00 Y 00 @ n þ qu00 n00 @ Y þ A ~v ~n þ ð1 þ AÞ k k @xk @xk
qu00k n00 ¼ q DT
ð15Þ
2
where DT = Clk /ScTe and C l ¼ 0:09. k and e are respectively the turbulent kinetic energy and its dissipation rate, and ScT is a turbulent Schmidt number set to its usual value 0.7.
The scalar dissipation rates which appear in the transport Eqs. (11), (13), and (14) represent the mean mixing rate at small scales. In the case of partially premixed combustion, it has been shown that the reaction rate depends on three dissipation rates [17]; ~ n ¼ 2q2n ¼ 2qD @n the mixture fraction dissipation rate v @x ~ Y ¼ 2q2Y ¼ 2qD @Y fuel mass fraction dissipation rate v @x ~ nY ¼ 2q2nY ¼ 2qD @n the cross-dissipation rate v @x
00 k
@Y 00 . @xk
00 k
00 k
@n00 , @xk @Y 00 @xk
and
The modelling
~e v~ n ¼ R ~ qn002 k
ð16Þ
where the constant R is the timescale ratio of the turbulence to scalar fluctuation. Its current is about 2 [18]. This model is efficient for passive scalars. However, the dissipation rate of the reactive scalar fluctuations is influenced by turbulence, chemical reaction and their interactions. As the classical model accounts only for turbulence it will not describe correctly the reactive scalar dissipation rate. The model discussed in [13,19] was derived by considering the behaviour of reactive scalar fluctuations in the limit of thickened flamelets and infinitely thin flamelets. A segregation factor S is employed to interpolate between these two regimes such as:
~e v~ Y ¼ 2q S2Y þ ð1 SÞR ~ Y~ 002 k
~ 1A @ ~n ~ v ~n002 x ~ n þ 2qu00k n00 þ ðh þ ð1 hÞAÞ @xk Y~ max Y~ min
qu00k Y 00
@ ~n @ Y~ qu00k n00 @xk @xk
The stainless steel cylindrical high pressure combustion chamber of ICARE can sustain premixed flames up to 1 MPa [14]. The chamber consists of two cylindrical halves of 600 mm high (Fig. 3) and 300 mm inner diameter. Each halve is equipped with four windows of 100 mm diameter designed for optical diagnostics. The burner with a 25 mm diameter is placed in the centre of the chamber and can be moved along the vertical axis by means of a stepping motor. The upper part of the burner has a convergent profile. A porous grid is placed in the bottom part of this burner section to ensure chamber-inlet turbulence. The main flame consists of a lean methane–air mixture with an equivalent ratio of 0.6. It is surrounded by a 2 mm pilot flame at stoichiometric conditions. 4. Numerical results and discussion The LW–P model described in the previous sections has been implemented in the computational fluid dynamic software Code_Saturne developed by EDF, Archambeau et al. [20]. It is a parallel code based on a finite volume method. The set of the considered equations consists of the Favre averaged Navier–Stokes equations
ð17Þ ð18Þ
where S is defined as:
S¼
qY 002 ðqY 002 Þmin ðqY 002 Þmax ðqY 002 Þmin
ð19Þ
and
qY 002
max
~ 2 ~n002 Y~ max Y~ Y~ Y~ min þ h þ ð1 hÞA ¼q
ð20Þ
and
~ 2 Þ~n002 ðqY 002 Þmin ¼ ðh þ ð1 hÞA
ð21Þ
~ n and h is defined as h ¼ ðY ~ Y ~ min Þ= ðh þ ð1 hÞAÞ where qY n ¼ q ~ max Y ~ min Þ and Y ~ max and Y ~ min are defined by the bounds of the ðY ~ is the mean slope of the equilibrium line such composition space. A 00 00
ð23Þ
3. The ICARE rig
~e v~ nY ¼ 2q S2nY þ ð1 SÞR ~ Y 00 n00 k
!
the
of these three dissipation rates is one of the most challenging problems in the simulation of partially premixed systems. The classical model describing the dissipation rate of the mix~ n is based on the analogy with the turture fraction fluctuations v bulent mixing [18]:
and
2nY ¼ n00 x 2q
2.3. Scalar dissipation terms
ð22Þ
~002
~ ¼ 0 for lean mixture and A ~ ¼ 1=ð1 n Þ for rich mixture, that A st where nst stands for the stoichiometric mixture fraction.
Fig. 3. ICARE test-rig view [14].
61
A. Benarous et al. / Fuel 124 (2014) 57–65
for incompressible flows completed with sub-models for turbulence and equations describing the additional scalars (enthalpy, species concentration, and possible user-defined scalars). The operating scheme is based on a prediction of the velocity field followed by a pressure correction step. Equations for turbulence and scalars are solved separately afterward. The space discretization is based on the fully conservative unstructured finite-volume framework, with a fully collocated arrangement for all variables. In the present work, turbulent mixing is represented by the two-equations (k–e) model. An unstructured grid of approximately 76,000 cells has been generated to represent half of the physical domain. The retained mesh has been obtained from several numerical tests on finer grids which did not influence the centreline static temperature field by more than 2% (Fig. 4a).
4.1. Boundary conditions
1800
1500
1200 Mesh 1 : 29378 Cells Mesh 2: 49348 Cells
ð24Þ
has been imposed. In Eq. (24), rb denotes the burner radius and u0 is the centreline velocity. For such ambient pressure case, the velocity exponent is taken to be n ¼ 0:2818 [15]. According to LDA measurements of the turbulent kinetic energy k [15], a polynomial expression was imposed at the main inlet flow, with a minimum value of k0 ¼ 0:047 m2 =s2 (Fig. 5b). The mean dissipation rate ~e has been chosen in such a way that the measured grid turbulence decay is recovered. 1.00
0.80
0.60
Velocity Profile at burner exit
0.40
Experimental Measurements at X/D = 0.2 [14]
0.20
Mesh 3: 76321 Cells
900
0.00 0.00 600
0.10
0.20
0.30
0.40
0.50
DimensionlessRadialDistance(R/D) Fig. 5a. Radial velocity profile at X/D = 0.2.
300
0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D)
0.10
Fig. 4a. Mesh study based on the centreline temperature.
Turbulent Kinetic Energy (m**2/s**2)
Centreline Static Temperature (K)
uðrÞ ¼ u0 ð1 r=r b Þn
Dimensionless Radial velocity (U/U0)
Inlet boundary conditions must be handled with special care, since they may have a strong influence on combustion inside the computational domain. This is often complicated by several constraints: (i) the measurements of all variables at the exact location of the numerical inlet boundary are not always available; (ii) also, an additional constraint results from the difficulty to measure
some transported variables related to the turbulent combustion model. This is especially true for the turbulence mean dissipation rate ~e. Boundary conditions on the left side of the computational domain are given by the flow symmetry (Fig. 4b). The right and top sides are approximated by outlet conditions, whereas an adiabatic wall condition is specified on the burner thickness (burner lips). To preserve a fixed bulk velocity ud ¼ 2:3 m/s [14] within the burner (Fig. 5a), a radial velocity profile:
0.09
TKE Profile at burner exit Experimental Measurements at X/D = 0.2 [14]
0.08
0.07
0.06
0.05
0.04 0.00
0.10
0.20
0.30
0.40
0.50
DimensionlessRadialDistance (R/D) Fig. 4b. The computational domain.
Fig. 5b. Radial profile of the turbulent kinetic energy at X/D = 0.2.
A. Benarous et al. / Fuel 124 (2014) 57–65
In order to recover the experimental values of k and u at the first measurement station X/D = 0.2, we have adjusted the numerical boundary conditions at the domain inlet. The pilot flame is modelled as a co-flow of hot flue gases being burned at stoichiometric conditions. Accordingly, a uniform and low value of the axial velocity was imposed at the second inlet whereas turbulent quantities were expressed in terms of intensity and hydraulic diameter. The numerical values for the boundary conditions are depicted in Table 1.
4.2. The nonreactive flow As a first step, it is interesting to perform a non-reactive flow computation. The so-called ‘‘cold flow simulation’’ helps to validate the appropriate turbulence model describing the mixing and tune the ‘‘numerical’’ boundary conditions. Moreover, solving the scalar transport equations without chemical source terms provides a ‘‘guess’’ solution field for the reactive calculations. Figs. 6 and 7 show respectively the centreline profiles of the mean velocity com~ along the axial coordi~ and the turbulent kinetic energy k ponent u nate. A quasi-linear decrease is noticed on the velocity profile until X/D = 3.2 where a slope change occurs. As explained by Hinze [21], this behaviour is due to the training effect characterised by an important momentum transfer from the jet to its surrounding. As the jet develops axially, the mixing zone thickens and eventually reaches the axis while affecting the centreline velocity. Good agreement between numerical and experimental data is observed for the dynamic field (Fig. 6). The predicted jet-core size (X/ D 3.2) is quite close to the measured value (X/D = 3.0) [15]. As stated above, turbulence decay along the chamber axis (Fig. 7) was obtained as a polynomial distribution for the mean dissipation
Table 1 Boundary conditions used in the simulation.
Velocity (m/s) k (m2/s2) e (m2/s3) T (K) Equivalence ratio
Main flow
Co-flow
u0 = 2.56 r2 polynomial r3 polynomial 300 0.6
0.05 ud k = k (I = 5) e = e (I = 5%) T (nst) 1.0
Dimensionless Axial Velocity (U/U0)
1.00
0.95
0.90
0.85 Cold Flow Calculations Experimental Measurements [15]
0.80
0.75
0.0
1.0
2.0
3.0
4.0
5.0
Dimensionless Axial Distance (X / D) Fig. 6. Centreline profile of the axial velocity.
6.0
0.13 0.12 Experimental Measurements [15]
Turbulent Kinetic Energy (m**2/s**2)
62
0.11
Cold Flow Calculation
0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) Fig. 7. Axial profile of the turbulent kinetic energy.
rate ~e at the domain inlet. This was done by imposing a constant turbulent length scale LT ¼ 2:9 mm which is quite close to the turbulence grid size [14]. Predicted values of the mean turbulent kinetic energy at various axial stations (Fig. 7) seem to fit well the ~ is reached at experimental measurements. Minimum value for k the potential core region, exactly at the axial location X=D 2:48 where its exhibits an underestimation of 17% compared with the measurements. This discrepancy is certainly due to the turbulent stresses isotropy considered for the (k–e) model formulation, which underestimates the jet spreading in both axial and radial directions. Owing to the important chamber-to-burner diameter ratio (rch/ rb = 12), the wall effects seem to be unimportant which further justifies the free character of the jet within the chamber.
4.3. Partially premixed combustion The more general case of partially premixed combustion is now considered. The study is performed in the case of the stratified condition with an equivalence ratio varying from unit value for the hot co-flow to 0.6 at the main flow entrance. In order to allow comparison with the reactive measurements, we define a pseudo-progress ~ ~ variable ~c ¼ 1 Y= n which accounts for both fuel consumption and turbulent mixing [22]. The previous expression is only valid for lean mixture; therefore the corresponding variable is zero in the fresh mixture and takes a unit value in the burnt gases. Regarding the dynamic variables in the presence of the flame, ~ and the turcentreline profiles of the mean velocity component u ~ are depicted in Figs. 8a and 8b respectively. bulent kinetic energy k It is noticed that the slope of the axial velocity decay upstream the flame front is quite similar to that of the non reactive case. Momentum transfer still exists in this region. However, the isobaric character of the flame yields an acceleration of the flow before the potential core region (X/D 1.5). Globally, the calculated profile of the axial velocity is in good agreement with the experimental results. The acceleration behind the flame induced by the thermal expansion is well predicted (Fig. 8a). However, the profile of the numerical turbulent kinetic energy is less accurate. Our calculation shows a relative underestimation up to 26% downstream the flame brush, compared with experimental measurements (Fig. 8b).
63
A. Benarous et al. / Fuel 124 (2014) 57–65 1.0
0.8
1.04
0.6 1.00
Flame Calculations
0.4
Progress
Dimensionless Axial Velocity (U/U0)
1.08
Experimental Measurements [15]
0.96
Pseudo - Progress variable
0.2
Fig. 9a. Spatial contour of the static temperature. 0.92
0.0
1.0
2.0
3.0
4.0
5.0
6.0
0.0
2400
Dimensionless Axial Distance (X/D)
Calculations at r/D = 0.56
Fig. 8a. Centreline velocity distribution.
2100
Calculations at r/D = 0.5
1800
Experimental Measurements [15]
1.0
0.05
0.8
0.04
Flame Calculations Pseudo-Progress variable
0.6
Experimental Measurements [15]
0.4
Progress
Turbulent Kinetic Energy (m**2/s**2)
Static Temperature (K)
Centreline Calculations
1500
1200
900
600
300 0.0
0.03
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) 0.2
0.02 0.0
1.0
2.0
3.0
4.0
5.0
0.0 6.0
Dimensionless Axial Distance (X/D) Fig. 8b. Turbulent kinetic energy distribution along the chamber axis.
~ in the flame region is certainly due The calculated decrease of k to the use of the first order (k–e) model in the closure of the scalar turbulent fluxes (Eq. (15)). Indeed, the model do not account for the flame generated turbulence resulting from the counter-gradient mechanism [23]. Spatial contour of the static temperature seems to have a conical shape (Fig. 9a). It is also noticed that in the present case, the coflow entrance is the hottest zone of the combustion chamber. This was expected since the co-flow mixture was burned at stoichiometric conditions, yielding an adiabatic value for the static temperature (2318 K). However, the co-flow does not seem to affect the thermal level of the flue gases throughout the chamber. This is a consequence of the low flow rate and the laminar character imposed at the co-flow inlet. Within the chamber, the flue gases are still at a much lower temperature (1645 K) than that of the co-flow entrance. Fig. 9b shows the static temperature profiles along the axial distance of the chamber at three radial locations, the centreline (r/
Fig. 9b. Static temperature profiles for various radial locations.
D = 0), the burner lips (r/D = 0.5) and the mid-entrance of the coflow (r/D = 0.56). Along the axial location X ¼ 0, the temperature profile moves from the highest values (2318 K), characteristic of the hot co-flow, to a value T~ ¼ 380 K induced by the slight thermal diffusion, and then finally to the lowest value (300 K) which refers to the temperature of the fresh gases. However, the tendency is different downstream, where the three profiles asymptotically reach a limiting value (1645 K) corresponding to the temperature of the products. The axial distance where the flue gas temperature is reached (Fig. 9b) decreases with the radial location which confirms the conical shape of the flame. Compared with the chamber centreline measurements, the static temperature exhibits a global overestimation of 6%. This discrepancy is probably due to the global reaction modelling adopted for the chemistry description. Indeed, the model considers only carbon dioxide CO2 and water vapour H2O as final products in the methane–air reaction and does not account for enthalpy losses by dissociation. The limiting values for ~c are clearly shown in the spatial contour depicted in Fig. 10a. The turbulent flame seems to be stabilized over the burner lips and exhibits a typical conical shape. Axial distribution of the progress variable agrees the measurements (Fig. 10b), although the experimental curve is slightly steeper than
64
A. Benarous et al. / Fuel 124 (2014) 57–65 0.6
0.4
LW-P Calculations Linear Fitting Y = (-1/a) X + b
Z / Flame Brush
0.2
0.0
-0.2
-0.4
Fig. 10a. Spatial contour of the mean progress variable. -0.6 -2.50 -2.00 -1.50 -1.00 -0.50 0.00
1.00
0.50
1.00
1.50
2.00
2.50
Ln((1-
Pseudo Progress Variable
0.80
5. Conclusion 0.60
0.40
Experimental Measurements [15] Numerical (LWP) Calculations
0.20
0.00 0.0
1.0
2.0
3.0
4.0
5.0
6.0
Dimensionless Axial Distance (X/D) Fig. 10b. Progress variable distribution along the chamber axis.
the numerical one. Accordingly, the maximum thickness of the conical flame dF can be evaluated by means of the centreline distribution of the progress variable. It is defined as the axial distance for which ~c increases from 0.01 to 0.99 [15]. Hence, the calculation shows an apparent flame thickness dF ¼ 62:5 mm slightly greater (11%) than the measurements [14]. As it was found by Robin et al. [22], the mean flame is most probably thicker than the experimental flame. This discrepancy is due to the gradient assumption considered on the scalar turbulent fluxes qu00 Y 00 . As already stated by Lipatnikov et al. [24], the axial progress variable profiles may be represented by a general form when using a reduced spatial coordinate across the flame brush. If we note this quantity Z ¼ X X ~c¼0:5 , it is possible to recover an exponential tendency of the progress variable along the flame thickness. Indeed, a linear evolution of Z=dF with the quantity Lnðð1 ~cÞ=~cÞ is clearly visible in Fig. 11. In order to obtain an analytic expression, we perform a linear fitting using the least square method. The corresponding slope a ¼ 4:20 is found to agree the experimental value as pointed out by Pavé (a = 4.28) [14] for a Bunsen premixed flame and Shepherd (a = 4.02) [25] for a V-shaped flame. Hence, it appears that the evolution of ~c in the flame brush may be represented by a single curve which seems to be fairly universal.
We have studied a partially premixed pilot-stabilized methaneair flame. The LW–P model based on a two-scalar ðY; nÞ four delta Pdf was used to describe mixing and reaction progress within the burner-chamber device. As it was stated in the experiments, the hot co-flow has allowed a stabilization of the main flame without affecting the global thermal level within the chamber. Moreover, the computed static temperature overestimates the measured one by only 6%, despite the use of a one-step chemistry. Spatial contour of the progress variable has clearly shown the conical shape of the main flame. Owing to the bridging expression used for the cross dissipation terms, values of ~c agree quite well with the experiments. In the cold flow situation, a good qualitative agreement with the experimental data is found, especially for the jet-core size and the turbulence decay along the chamber axis. Moreover, the mean velocity profile downstream the flame brush fits well the experiments despite a slight underestimation in the fresh gas region. The discrepancy is due to the fact that the adopted numerical conditions at the burner exit were issued from the cold flow case. However, it is shown that the present version of the LW–P model does not improve the prediction of turbulent kinetic energy downstream of the flame brush. The numerical depletion observed on the mean turbulent kinetic energy is probably due to the use of a first order model (k–e), a model that does not take into account the combustion induced turbulence. Some efforts are still necessary to improve the turbulent mixing by a second-order model. Such a model will allow us to deal with the flame-generated turbulence and counter-gradient diffusion in partially premixed conditions. The extension of the present model to include detailed chemistry effects still remains a challenging task. References [1] Peters N. Turbulent Combustion. 1st ed. Cambridge: Cambridge University Press; 2000. [2] Ra Y, Cheng WK. Laminar flame propagation through a step-stratified charge. In: Proceedings of the Fifth Int. Symp on Diagnostics and Modeling of Combustion in Internal Engines, 2001, p. 251–57. [3] Marzouk YM, Ghoniem AF, Najm HN. Dynamic response of strained premixed flames to equivalence ratio gradients. Proc Combust Inst 2000;28:1859–66.
A. Benarous et al. / Fuel 124 (2014) 57–65 [4] Mura A, Galzin F, Borghi R. A unified PDF-flamelet model for turbulent premixed combustion. Combust Sci Technol 2003;175(7):1–37. [5] Renou B, Boukhalfa A, Puechberty D, Trinité M. Local scalar flame properties of freely propagating premixed turbulent flames at various Lewis numbers. Combust Flame 2000;123:l507–1521. [6] Jimenez C, Cuenot B, Poinsot T, Haworth DC. Numerical simulation and modeling for lean stratified propane–air flames. Combust Flame 2002;128:1–21. [7] Bray KNC, Champion M, Libby PA. Bray-Moss thermochemistry revisited. Combust Sci Technol 2002;174(7):167–74. [8] Robin V, Mura A, Champion M, Plion P. A multi-Dirac presumed PDF model for turbulent reactive flows with variable equivalence ratio. Combust Sci Technol 2006;178:1843–70. [9] Bondi S, Jones WP. A combustion model for premixed flames with varying stoichiometry. Proc Combust Inst 2002;29:2123–9. [10] Helie J, Trouve A. A modified coherent flame model to describe turbulent flame propagation in mixtures with variable composition. Proc Combust Inst 2000;28:193–201. [11] Kolla H, Rogerson J, Chakraborty N, Swaminathan N. Scalar dissipation rate modeling and its validation. Combust Sci Technol 2009;181:518–35. [12] Ribert G, Champion M, Plion P. Modeling turbulent reactive flows with variable equivalence ratio: application to the calculation of a reactive shear layer. Combust Sci Technol 2004;176:907–23. [13] Robin V, Mura A, Champion M, Degardin. O, Renou B, Boukhalfa M. Experimental and numerical analysis of stratified turbulent V-shaped flames. Combust Flame 2008;153:288–315. [14] Pavé D. Contribution to study the structure of a lean premixed turbulent methane–air flame. PhD thesis, Orléans University, France, 2002.
65
[15] Lachaux T. Pressure effects on the structure and dynamics of lean premixed turbulent methane–air flames. PhD thesis, Orléans University, France, 2004. [16] Fernandez-Tarrazo E, Sanchez A, Linan A, Williams FA. A simple one-step chemistry model for partially premixed hydrocarbon combustion. Combust Flame 2006;147(1–2):32–8. [17] Bray K, Domingo P, Vervich L. The role of progress variable in models for partially premixed turbulent combustion. Combust Flame 2005;141:431–7. [18] Spalding DB. Mixing and chemical reaction in confined turbulent flames. Chem Eng Sci 1971;26:95–107. [19] Mura A, Robin V, Champion M. A second-order model for turbulent reactive flows with variable equivalence ratio. Combust Flame 2007;149:217–24. [20] Archambeau F, Mehitoua N, Sakiz M. A finite volume code for the computation of turbulent flows. Int J Finite, vols. 1–62. 2004. [21] Hinze JO. Turbulence. 2nd ed. New York: McGraw Hill Series in Mechanical Engineering; 1975. [22] Robin V, Guilbert N, Mura A, Champion M. Turbulent combustion modeling of mixtures with inhomogeneous equivalence ratio: application to a flame stabilized by the sudden expansion of a 2-D channel. Comptes Rendus de Mécanique 2010;338:40–7. [23] Robin V, Mura N, Champion M. Algebraic models for turbulent transports in premixed flames. Combust Sci Technol 2012;184:1718–42. [24] Lipatnikov AN, Chomiak J. Turbulent flame speed and thickness: phenomenology, evaluation and application in multi-dimensional simulations. Prog Energ Com Sci 2002;28:1–73. [25] Shepherd I.G. Flame surface density and burning rate in premixed turbulent flames. In: 26th Symposium (International) on combustion; The Combustion Institute, 1996. p. 373–79.