Large eddy simulation of turbulent combustion by a dynamic second-order moment closure model

Fuel 187 (2017) 457–467

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Large eddy simulation of turbulent combustion by a dynamic second-order moment closure model Kun Luo, Jianshan Yang, Yun Bai, Jianren Fan ⇑ State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou 310027, PR China

h i g h l i g h t s
A dynamic second-order moment closure combustion model is developed.
The model is applicable for all combustion regimes.
Two experimental flames are predicted by this model.

a r t i c l e

i n f o

Article history: Received 26 April 2016 Received in revised form 13 August 2016 Accepted 26 September 2016 Available online 7 October 2016 Keywords: Turbulent combustion Large eddy simulation Second-order moment closure Sandia flame D Swirling non-premixed flame

a b s t r a c t A dynamic second-order moment closure model is developed for turbulent combustion in the form of large eddy simulation. The filtered reaction rate is directly closed in the form of Arrhenius law, and the whole temperature exponential function is treated as a single variable to avoid the traditional Taylor series expansion. The sub-grid unresolved reaction rate is modeled with a second-order moment closure model. All the coefficients in the sub-grid models are evaluated by the dynamic procedures. To validate and evaluate this model, a priori validation using a DNS database and posteriori validation by LES of the Sandia piloted jet flame (Flame D) and the Sydney bluff-body swirling flame (SM1) are performed. The results demonstrate that the dynamic second-order moment closure model coupled with LES is able to reasonably predict turbulent combustion even with simple chemistry, and has potential to predict more complex combustion with detailed reaction mechanism and acceptable computational cost. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction In most industrial equipment, such as coal-fired furnaces, combustion gas turbines and internal combustion engines, turbulent combustion has been encountered [1]. Predictive capacity for turbulent combustion is thus of paramount importance to those applications in which experimental measurements are difficult, limited and even not accessible. In the recent years, much more attention has been paid to large eddy simulation (LES) of turbulent combustion, since LES has been demonstrated to be able to provide more accurate information and useful qualitative understanding of turbulent combustion phenomena, compared to the traditional ways on the basis of the Reynolds Averaged Navier-Stokes [2–4]. The most challenging issue in turbulence combustion modeling is proper treatment of the multi-scale coupling process of chemical reaction and molecular diffusion. Many of the current approaches

⇑ Corresponding author. E-mail address: [email protected] (J. Fan). http://dx.doi.org/10.1016/j.fuel.2016.09.074 0016-2361/Ó 2016 Elsevier Ltd. All rights reserved.

could be categorized as PDF-like or flamelet-like based on how to deal with the coupling between chemical reactions and molecular diffusion [5]. The steady flamelet model (SFM) put forward by Peters [6,7] is one of the flamelet-like models. It is relatively simple and has been improved in different methods, such as the unsteady flamelet model (UFM) [8], the flamelet/progress variable model (FPV) [9] and so on. The advantage of the flamelet model is that detailed chemistry can be considered with reasonable computational cost. However, it was basically formulated and developed separately for premixed flames and non-premixed flames. The probability density function (PDF) models have been well developed for combustion by Pope [10,11]. Some PDF-like models have also been developed such as the multiple mapping conditioning (MMC) [12] and the one-dimensional turbulence (ODT) model [13]. One advantage of the PDF method lies in the fact that it can be applied to all combustion regimes, even challenging flames, like extinction, ignition flames and bluff-body swirling flames [14]. But the relatively simple mixing models and high computational cost are still the weaknesses in the current application of PDF models.

458

K. Luo et al. / Fuel 187 (2017) 457–467

However, it has been demonstrated that both non-premixed flames and premixed flames may coexist in either traditionally defined premixed turbulent combustion or non-premixed combustion because of turbulence-chemistry interactions [15]. To predict both premixed and non-premixed flames, Knudsen and Pitsch [16] proposed a multi-regime flamelet model to consider different combustion regimes. This model was further extended to study IC engines by Mittal et al. [17]. In this multi-regime flamelet model, models for both premixed and non-premixed flames are combined to account for turbulent combustion with a regime indicator. As we mentioned above, the current combustion models, which could be accounted for different combustion regimes, are very difficult for a freshman in CFD to make use of. In this paper, a dynamic second-order moment closure (DSMC) model is developed to account for not only premixed but also non-premixed combustion in the context of large eddy simulation. Also for a beginner in LES of turbulent combustion, there is no difficulty to know well about the model. The chemical reaction rate term is directly closed in the form of Arrhenius law, while the whole temperature exponential function is treated as a variable so that the Taylor series expansion is avoided. The coefficients in the sub-grid model are determined by dynamic procedures. This model is then validated and evaluated by a priori validation using a direct numerical simulation (DNS) database and posteriori validation by LES of the Sandia piloted flame (Flame D) and the Sydney bluff-body swirling flame (SM1), as these flames have been widely selected as benchmark tests to validate various combustion models

3. The DSMC combustion model Considering a chemical system with M reactions, if the filtered reaction rates of all the reactions are evaluated accurately, the reaction rates of the scalars could be calculated. To highlight the effects of the filtering, we consider a rate expression of the form m  n   x_ ðq; Y F ; Y O ; TÞ ¼ A exp  RTE q WY FF q WY OO . Here m and n are the concentration exponents depending on the chemical kinetics. W denotes molecular weight. Y stands for the mass fraction. A and E denote the pre-exponential factor and the activation energy respectively. R is the universal gas constant. The subscripts F and O represent the fuel and oxidizer respectively. The filtered reaction rate term could be expressed as:

m  n  YF YO
q
q WF WO  E
Rn1
Rn2 ¼ Bexp  RT 

x_ ¼ Aexp 

E RT

ð5Þ

n where R1 ¼ ðqY F Þ; R2 ¼ ðqY O Þ; B ¼ AW m F W O . In the previous study, the temperature function is usually treated by the Taylor series expansion, and the filtered source term becomes so complicated that significant error appears [21]. To avoid this problem, the whole temperature function is treated as a variable K in the present study,  E namely K ¼ exp  RT . For Ri ði ¼ 1; 2Þ, the variance R0i is assumed to be much smaller compared with Ri , then the expression below can be regarded feasible:

2. Governing equations For low-Mach number and variable-density flow, the filtering continuity, momentum, species and energy equations for LES of turbulent combustion read [18]:

  m m1 0 m Rm  ðRi Þ þ mðRi Þ R0i : i ¼ Ri þ Ri

ð6Þ

Ignoring the third-order unresolved terms, the reaction rate turns into:



x_ ¼ BKRm1 Rn2  B KR1 m R2 n þ mnR1 m1 R2 n1 KR01 R02 þ mR1 m1 R2 n K 0 R01

~j Þ @ðqÞ @ðqu þ ¼ 0; @t @xj

ð1Þ



 ~ i Þ @ðqu ~i u ~j Þ @ðqu @p @ rij @  ~i u ~j Þ ; þ ¼ þ  qð ug i uj  u @t @xj @xi @xj @xj  @ðqY~k Þ @ðqY~k uj Þ @ @ Y~ þ ¼ qDk k @t @xj @xj @xj

!

ð2Þ

@ðq Te Þ @ðq Te uj Þ @ @ Te ¼ qDT þ @xj @xj @xj @t

ð7Þ

Neglecting the third-order variance correlation term, the mean reaction rate reads:   x_ ¼ B KR1 m R2 n þ mnR1 m1 R2 n1 KR01 R02 þ mR1 m1 R2 n K 0 R01 þ nR1 m R2 n1 K 0 R02 : ð8Þ

i  @ h g  q Y k uj  uj Y~k þ x_ k ; @xj ð3Þ



 þ nR1 m R2 n1 K 0 R02 þ mnR1 m1 R2 n1 K 0 R01 R02 :

! 

@ h g e i qð Tuj  uj T Þ þ x_ T ; @xj

ð4Þ

where q and u are the density and velocity of mixture gas respectively. p is the pressure of the mixture gas. Y k represents the mass _ k and x _ T are the filfraction of specie k, and T is the temperature. x tered chemical reaction and heat release source terms. The sub-grid momentum and scalar transport terms in the above equations are modeled by dynamic approach [19,20]. The turbulent stresses are calculated with an eddy viscosity assumption [20]. The residual scalar fluxes are modeled with the gradient-diffusion assumption [1]. The coefficients in the subgrid models are all calculated by using the dynamic procedures [19]. The closure for chemical reaction source term in Eq. (3) and heat release source term in Eq. (4) will be introduced in the next section in detail.

As Ri can be calculated from the transport Eq. (3), the closure problem of chemical source term is changed into how to model the second-order variance correlation terms R01 R02 , K 0 R01 , K 0 R02 and the temperature function K. The highly non-linearity exponential function K is modeled using the fluctuation averaged closure method as:

8 3 39 2 2 > > > > < 7 7= 6 6 1 E E 7 7 6   þ exp  exp 6  K¼ 1 1 5 5>: 4 4   2> 00 00 2 00 00 2 > > : ; T T R T þ Tg R T  Tg

ð9Þ

00 00 The variance of temperature Tg T is modeled using the sub-grid variance model with the dynamic procedure to determine the coefficient [20]:

00 00 q Tg T ¼ C T qD2 jr Te j2 :

ð10Þ

The transport equation of R01 R02 can be derived following the similar method of deriving the transport equation of Reynolds stress as [22]

459

K. Luo et al. / Fuel 187 (2017) 457–467 0

0

þ R2 R01

@u0j

2200

_ þ tO W O R01 x _ : þ tF W F R02 x 0

0

ð11Þ

Multiplying the temperature Eq. (4) by expðE=RTÞ RE

1 , T2

one

could get the transport equation of K. Following similar procedure, the transport equation for the variance correlation K 0 R01 can be derived as: 0 0 0 @u0j @K 0 R01 @uj K 0 R01 @K uj R1 @K @R1 þ þ u0j R01 þ R1 K 0 þ K 0 u0j þ @xj @t @xj @xj @xj @xj

 0 @ @Y ¼ K0 qD1 1 þ tF W F K 0 w_ 0 þ R01 A01 þ R01 A02 þ R01 A03 ; @xj @xj

ð12Þ

ð13Þ

A2 ¼

 1 @K @ k ; q @xj @xj C p

ð14Þ

1

expðE=RTÞ

E 1 x_ T : R T2

2000 1800

0.1

1600 0.05 1400 1200 0

  k @ @K E 1 @K @T 2 @K @T  ; þ 2 R T @xj @xj T @xj @xj qC p @xj @xj

q

0.15

0

A1 ¼

A3 ¼

2400

þ R01 R02

mass fraction

@u0j

0.2

@u0j @uj þ R01 R02 @xj @xj @xj @xj

 0

 0 @ @Y @ @Y ¼ R02 qD1 1 þ R01 qD2 2 @xj @xj @xj @xj

þ R1 R02

temperature C7H16 O2

temperature (K)

0

@ðR01 R02 Þ @uj R01 R02 @R2 @R1 @uj R1 R2 þ þ u0j R01 þ u0j R02 þ @t @xj @xj @xj @xj

0.2

0.4

x/L

0.6

0.8

1

Fig. 1. The distribution of the section-averaged temperature and mass fractions in the x direction.

ð15Þ

To further simplify the transport equations, unity Lewis number is assumed. Then

 0 

 0  @ @R1 @ @R2 þ R02 D1 R01 D2 @xj @xj @xj @xj

 0 0  0 @ @R R @R @R01  2D 2 D 1 2 ; ¼ @xj @xj @xj @xj

ð16Þ

  _ 0 ¼ R1 w _  R1 w _ R01 w   0 0 n1 m n 0 0 m n1 0 0 0 ¼ B nRm K R1 R2 : 1 R2 KR1 R2 þ R1 R2 K R1 þ nR1 R2

ð17Þ

Keeping only the leading order term, we have

tF W F R02 x_ 0 þ tO W O R01 x_ 0 ¼ C 1 R01 R02 :

ð18Þ

n1 where C 1 ¼ tF W F BmRm1 Rn2 K þ tO W O BnRm 1 1 R2 K. Casting the little

terms, the transport equation of R01 R02 is reduced to:

  @ R01 R02 @t

þ

@uj R01 R02 @xj

@R01 R02

Following similar procedure, the transport equation of K 0 R01 can be simplified as:

!

@R2 @R1 @ ¼ u0j R01  u0j R02 þ D @xj @xj @xj @xj |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

! @R01 @R02 þ C 1 R01 R02 : |fflfflfflffl{zfflfflfflffl} @xj @xj E4 |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}

E1

@K 0 R01 @uj K 0 R01 @K @R1 @ @K 0 R01 ¼ C4 þ D þ @xj @xj @xj @t @xj @xj

E2

 2D

ð19Þ

E3

The right side of the above equation consists of production term (E1 ), diffusion term (E2 ), dissipation term (E3 ), and reaction term (E4 ). Introducing the similar modeling strategies for the production term and dissipation term as those in Reynolds stress model [23], the transport equation of R01 R02 is finally reduced to:

@ðR01 R02 Þ @uj R01 R02 @R1 @R2 @ @R0 R0 ¼ C2 þ D 1 2 þ @t @xj @xj @xj @xj @xj

!

þ C 3 R01 R02 :

Fig. 2. Comparison of the filtered reaction rate calculated using Eq. (8) when ignoring the third-order unresolved correlation term with the filtered DNS chemical reaction rate.

ð20Þ

!

þ C 5 K 0 R01 :

ð21Þ

It should be noted that the terms C 1 , C 2 , C 3 , C 4 , C 5 in Eqs. (18)– (21) are all dimensional terms. In general, the transport equations of second-order variance correlation terms could be closed in the form as:

@/0 u0 @uj /0 u0 @ @/0 u0 @/ @ u ¼ C g1 þ C g2 þ C g3 /0 u0 ; þ @xj @xj @xj @xj @t @xj

ð22Þ

where / and u stand for R1 , R2 and K. C g1 , C g2 and C g3 are also dimensional terms. In principle, the above transport equations can be solved to get the second-order variance correction terms with closure models. However, we seek simple algebraic expression in the present work. In analogy to the method of the algebraic stress model (ASM)

460

K. Luo et al. / Fuel 187 (2017) 457–467

developed from the Reynolds stress model [22], the modeled equations can be further simplified to the algebraic forms by assuming a steady state and the convection and diffusion terms are in balance:

/0 u0 ¼ 

C g2 @/ @ u @/ @ u ¼ C / u D2 : @xj @xj C g3 @xj @xj

ð23Þ

The non-dimensional coefficient C /u is determined by a dynamic procedure as

q/0 u0 ¼ qð/u  /uÞ ¼ C /u D2 q

@/ @ u ; @xj @xj

ð24Þ

and

C /u ¼

^

hLMi ; hMMi

^ b / u; L ¼ qd /u  q

In order to validate the combustion model, a priori validation using a DNS database and posteriori validation by LES of the Sandia piloted jet flame (Flame D) and the Sydney bluff-body swirling flame (SM1) are performed. In LES, to generate turbulent inflow conditions for the Sandia Flame D and the Sydney SM1, separate large eddy simulations of pipe flow are conducted, and the data from the pipe flow is read in as the inflow boundary conditions for the combustion LES. The governing equations are solved by using a staggered space-time, conservative discretization scheme [24]. A second-order semi-implicit iterative technique is used for time integration, while a second-order finite difference scheme, based on kinetic energy conserving is used for spatial discretization. 4.1. A-priori validation

^

^ d b b 2 @ / @ u  qD2 o/ ou : M¼q oxj oxj @xj @xj

4. Results and discussions

ð25Þ

To sum up, with these second-order variance correlation terms available, the filtered reaction rate can be obtained. The chemical source term and heat release source term in the transport Eqs. (3) and (4) are then available, and the equations are closed. In theory, the chemical source term is directly closed without any premixed or non-premixed assumptions and this model can be applicable to predict turbulent combustion. However, some assumptions in the model should be carefully checked. In the next section, a-priori validation using a DNS database and posterior validation by LES are presented to validate this model.

In the present DSMC model, the assumptions of neglecting the third-order correlation term and the algebraic form closure for the second-order variance correlation terms are applied. To validate these assumptions a-priori, a DNS database for turbulent combustion is utilized in this section. The DNS is carried out by Wang et al. [25] using a one-step reaction global reaction for heptane droplet evaporation and combustion in an isotropic turbulent flow. The computational domain is 2p mm in each direction. The grid spacing D0 is 49 lm and the total grid number is 1283. The droplets are initially distributed homogeneously in the slab between 3/8 L and 5/8 L in the x direction, and then evaporate, with gaseous fuel ignition and combustion. For more information, please refer to the previous studies [25,26]. To eliminate the effect of droplets, the

Fig. 3. The correlations of unresolved correlation terms and algebraic modeling terms for (a) R01 R02 , (b) K 0 R01 , (c) K 0 R02 . Top: filtered by 2D0 ; bottom: filtered by 4D0 .

K. Luo et al. / Fuel 187 (2017) 457–467

dataset at the time when the heptane droplets disappear and the flame becomes a purely gaseous flame are utilized to validate the DSMC combustion model. In the DNS, the following one-step chemistry [27] is used:

11O2 þ C7 H16 ! 7CO2 þ 8H2 O:

ð26Þ

The chemical reaction rate is



x_ ¼ A expðE=RTÞ
q

Y1 W1

m  n Y2
q W2

ð27Þ

where subscripts 1 and 2 stand for C7H16 and O2 respectively, and m and n are 0.25 and 1.5 respectively. The pre-exponential factor A is 5:1  1011 cm3 =mol  s and the activation energy E is 30:0 kcal=mol. Fig. 1 presents the distribution of the y-z plane-averaged temperature and mass fractions in the x direction. It is obvious that chemical reaction appears in the middle slab area with high temperature and high concentrations of reactants. There are flame

461

fronts which separate the fresh air and the burnt gas on the left and right sides of the reaction area. Fig. 2 compares the averaged _ calculated using Eq. (8) with the averaged chemical reaction rate x DNS chemical reaction rate. It can be found that the second-order closure model is able to reasonably represent the chemical reaction rate, and the influence of the third-order variance correlation term can be neglected although some deviations are observed in the middle strong reaction area. Thus the assumption of ignoring the third-order variance correlation term in the DSMC combustion model is believed to be acceptable. The algebraic form closure for the second-order variance correlation terms is validated by testing the correlations of variance correlation terms and their algebraic modeling terms. To be consistent with the LES method, the values in the coarse grids are obtained by filtering the DNS data. The DNS data are filtered by two different gird sizes, 2D0 and 4D0 (D0 is the DNS gird size) to get the mean values and variances with a box filter. Fig. 3 illustrates the correlations

00 00 Fig. 4. The correlations of unresolved correlation terms and modeling terms for (a) K and (b) Tg T . Top: filtered by 2D0 ; bottom: filtered by 4D0 .

462

K. Luo et al. / Fuel 187 (2017) 457–467

Co-flow 144mm Pilot Jet

18.2mm

7.2mm

(a)

(b)

(c) Fig. 5. The sketch map of the geometry configuration and computational mesh for Sandia flame D.

Table 1 The detailed setup for four cases. Case Case Case Case Case

1 2 3 4

Grid Points (x, r, h)

Cl

240  120  32 = 921 600 300  160  64 = 3 072 000 360  160  128 = 7 372 800 300  160  64 = 3 072 000

Dynamic Smagorinsky Dynamic Smagorinsky Dynamic Smagorinsky 0.05

of the DNS second-order variance correlation terms and their algebraic closure modeling terms. In the figure, the SMC /0 u0 term is cal@/ culated by @x j

@u @xj

and the DNS /0 u0 term is given as:

/0 u0 ¼ /u  /u;

ð28Þ

where / and u are the filtered mean values. Fig. 3 shows that the variance correlation terms are really proportional to the products of the two filtered mean values’ gradients, although scatters appear in coarser mesh. This justifies the assumption of the algebraic closure modeling of the second-order variance correlation terms in certain range of grid size. Furthermore, in the present DSMC combustion model, the highly non-linearity exponential function K is closed using the fluctuation averaged closure method as shown in Eq. (9), and the 00 00 T is modeled using Eq. (10). To examvariance of temperature, Tg ine these treatments, the correlations of DNS and SMC parts are 00 00 e j2 , T is calculated by jr T shown in Fig. 4. In Fig. 4(a), the SMC Tg 00 00 00 00 g g fT eT e , where T e is the and the DNS T T is given as T T ¼ TT filtered Favre average for T. In Fig. 4(b), the SMC K is calculated 8 2 3 2 39 > > > > < 6 7 6 7= 1 E E 6 7 6 7   , and the DNS K þ exp 4 by 2 exp 4 1 1 5 5 > >   2 2 > > ; : R T þg R T g T 00 T 00 T 00 T 00  E is the filtered mean value K ¼ exp  RT . Fig. 4 shows that the

Fig. 6. Axial profiles of the mean and rms temperatures in the centerline for all of the cases.

K. Luo et al. / Fuel 187 (2017) 457–467

463

Fig. 7. Radial profiles of mean and rms for velocities and O2 mass fraction at three axial locations for case 2.

Fig. 8. Instantaneous distribution of temperature at the plane of z = 0. (a) Experiment, (b) reduced two-step mechanism, (c) one-step mechanism. The black line represents the iso-line of stoichiometric mixture fraction.

00 00 gradient closure model for Tg T is acceptable with linear correlations, and the fluctuation averaged closure model for K is very good with excellent accuracy. The above a-priori validations demonstrate that the assumptions in the present DSMC combustion model are reasonable. To further validate the model, large eddy simulations of the Sandia flame D and Sydney bluff-body stabilized swirling flame are performed a-posteriori in the next sections.

4.2. Sandia Flame D The geometry configuration and computational mesh for the Sandia Flame D are depicted in Fig. 5. The central jet’s diameter is D = 7.2 mm while the diameter of outer is of 18.2 mm. The central jet is injected with the mixture gas of fuel and air (75% air and 25% CH4 by volume) in the velocity of 49.6 m/s. The velocity of the pilot flame is 11.4 m/s. The velocity of co-flow air is 0.9 m/s. In this

464

K. Luo et al. / Fuel 187 (2017) 457–467

Fig. 9. Radial profiles of the mean and rms axial velocities at different axial locations. (Symbols) experiment, (black solid line) reduced two-step mechanism and (grey Long Dash) one-step mechanism.

simulation, hydrogen–air combustion products are set to be the pilot flame with a temperature of 1880 K. The computational domain is taken as a cylindrical domain with an axial length of 0.6 m (about 83D) and a radius of 0.072 m (10D). To check the grid independence, three cases with different grid resolutions are performed. A constant C l is also used in case 4 to check the influence of the dynamic Smagorinsky model to the flow field computation. The detailed setup for these four cases is listed in Table 1. The chemistry in the simulation is described by a one-step irreversible overall reaction for methane [23]. Fig. 6 presents the axial distributions of mean and root mean square (rms) temperatures along the centerline. It is obvious that both the mean and rms temperatures are in good agreement with the experiment data except the case 1 with coarser grid. Cases with finer grids show few differences which means that the computational results obtained are independent of the grid resolution. Thus we take the medium grids in case 2 in our further simulations. Meanwhile, the predicted results in case 4 where the flow filed are calculated without a dynamic procedure also show a reasonable agreement with the data from experiment. This indicates that

combustion model plays a vital role in the prediction of the flame no matter the dynamic procedure is used for flow calculation. Fig. 7 presents the radial profiles of mean and rms of axial velocities and the O2 mass fraction for case 2. It is evident that the predicted mean and rms values agree well with the relevant measured data at different axial locations. This suggests that the proposed LES DSMC model is able to reasonably reproduce Flame D even with a simple chemistry. 4.3. Bluff-body stabilized swirling non-premixed flame (SM1) To further examine the predictive capability of the LES DSMC model for turbulent combustion, the bluff-body stabilized swirling non-premixed flame (SM1) is studied in this section. The central jet’s diameter is D = 3.6 mm. The pure fuel (100% CH4) is injected through the central jet with a velocity of 32.7 m/s. In the 5 mm wide annulus, which is surrounded the bluff-body, whose diameter is 50 mm, swirl air is injected. The bulk axial velocity and bulk tangential velocity are 38.2 m/s and 19.1 m/s. The geometrical swirl numbers is 0.5. Co-flow air is injected with the velocity of

K. Luo et al. / Fuel 187 (2017) 457–467

465

Fig. 10. Radial profiles of the mean and rms azimuthal velocities at different axial locations. (Symbols) experiment, (black solid line) reduced two-step mechanism and (grey Long Dash) one-step mechanism.

20 m/s. In the simulation, the computational domain is taken as a cylindrical domain with an axial length of 300 mm (about 83D) and a radius of 75 mm (about 21D). The mesh is stretched in the axial and radial directions. The mesh is consisted of 240  160  64 in the axial, radial and azimuthal directions after grid-independency study. Both a one-step irreversible reaction and a reduced twostep reaction rate mechanism [27] are used here to highlight the influence of chemistry. Fig. 8 presents the distribution of the instantaneous distribution of temperature in the central plane of z = 0. The black iso-lines are the stoichiometric mixture fraction, which are superimposed in the figure. It can be found that a fully developed turbulent flame is generated with an axial length of 0.15 mm reaction region. One can observe a necking region around the axial location of X = 50 mm, where the flame cross-sectional area is reduced, and then increase in the tail region for both cases. These predicted flame structures are similar to those observed in the experiment. Figs. 9 and 10 present the radial profiles of the mean and rms axial and azimuthal velocities at different axial locations. It can be found that both mean and rms velocities are in reasonable agreement with the measurement data, though the axial velocities

near the centerline are slightly over-predicted in the positions of X = 20 mm and X = 40 mm, which may be caused by the improper inflow boundary conditions. The two-step mechanism performs a little better than the one-step chemistry. The radial profiles of the mean temperature, mixture fraction and H2O mass fraction are depicted in Fig. 11. It’s found the LES DSMC model can correctly predict the mean mixture fraction, temperature and H2O mass fraction in general, except a slight overestimation of the mixture fraction in the jet central region at X = 20 mm. The mean temperature is also overestimated, but the results from the reduced two-step mechanism are better than those from the one-step mechanism, indicating that reduced chemistry is able to improve the predictive capability of the present model. In addition, the over-prediction of temperature is associated with the absence of the radiation effect in the current simulations. The above results demonstrate that the proposed LES DSMC model has reasonable predictive capability for turbulent combustion even with simple chemical mechanism. The predictive capability can further improved by more detailed chemistry. In the present work, the averaged computation time per step is about

466

K. Luo et al. / Fuel 187 (2017) 457–467

Fig. 11. Radial profiles of mean temperature, mixture fraction and H2O mass fraction at various axial locations. (Symbols) experiment, (black solid line) reduced two-step reaction rate mechanism and (grey Long Dash) one-step finite-rate mechanism.

K. Luo et al. / Fuel 187 (2017) 457–467

9.21 s and 10.75 s for the one-step mechanism case and the reduced two-step mechanism case. Thus the present model has the potential to be extended to detailed chemistry with reasonable computational cost for turbulent combustion. 5. Conclusions A dynamic second order moment turbulent combustion model has been developed in the context of large eddy simulation for turbulent combustion. The filtered reaction rate is directly calculated in the form of Arrhenius law. Of which, the whole temperature function is treated as a single variable so that the large truncation errors induced by the temperature function series expansion can be avoided. The unresolved reaction rate is modeled with the second order moment closure. All the coefficients in the sub-grid models are estimated by dynamical procedures. To validate and evaluate this model, a-priori validation using a DNS database and LES of the Sandia Flame D and the Sydney SM1 flame with simple chemistry are performed. The assumptions in the model are checked, and the general agreement with experimental data demonstrates the capability and efficiency of the proposed LES DSMC turbulent combustion model. As the chemical reaction source terms are directly closed, this model has the potential to predict premixed, non-premixed and partially premixed flames. Its predictive capability can be further improved by using complex chemistry with acceptable computational cost. Acknowledgement This work was supported by the National Nature Science Foundation of China (No. 91541202) and the National Key Research & Development Plan (No. 2016YFB0600102). We appreciate the inspired and useful discussions on the model with Dr. Andreas Kempf at the Universität Duisburg-Essen, Dr. Benoît Fiorina at the Laboratoire EM2C, Dr. Nils Haugen at the SINTEF Energy Research, and Professor Lixing Zhou at Tsinghua University. References [1] Poinsot T, Veynante D. Theoretical and numerical combustion. 2nd ed. Philadelphia: Edwards; 2005. [2] Pitsch H. Improved pollutant predictions in large-eddy simulations of turbulent non-premixed combustion by considering scalar dissipation rate fluctuations. Proc Combust Inst 2002;29:1971–8.

467

[3] Jones WP, Marquis AJ, Vogiatzaki K. Large-eddy simulation of spray combustion in a gas turbine combustor. Combust Flame 2014;161:222–39. [4] Knudsen E, Pitsch H. A dynamic model for the turbulent burning velocity for large eddy simulation of premixed combustion. Combust Flame 2008;154:740–60. [5] Pope SB. Small scales, many species and the manifold challenges of turbulent combustion. Proc Combust Inst 2013;34:1–31. [6] Peters N. Laminar diffusion flamelet models in non-premixed turbulent combustion. Prog Energy Combust Sci 1984;10:319–39. [7] Peters N. Laminar flamelet concepts in turbulent combustion. Symp (Int) Combust 1988;21:1231–50. [8] Haworth DC, Drake MC, Pope SB, Blint RJ. The importance of time-dependent flame structures in stretched laminar flamelet models for turbulent jet diffusion flames. Symp (Int) Combust 1989;21:589–97. [9] Pierce CD, Moin P. Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J Fluid Mech 2004;504:73–97. [10] Pope SB. The probability approach to the modelling of turbulent reacting flows. Combust Flame 1976;27:299–312. [11] Pope SB. PDF methods for turbulent reactive flows. Prog Energy Combust Sci 1985;11:119–92. [12] Klimenko AY, Pope SB. The modeling of turbulent reactive flows based on multiple mapping conditioning. Phys Fluids 2003;15:1907–25. [13] Kerstein AR. One-dimensional turbulence: Model formulation and application to homogeneous turbulence. J Fluid Mech 1999;392:277–334. [14] Masri AR, Pope SB, Dally BB. Probability density function computations of a strongly swirling nonpremixed flame stabilized on a new burner. Proc Combust Inst 2000;28:123–31. [15] Wang H, Luo K, Yi FX, Fan JR. Direct numerical simulation study of an experimental lifted H2/N2 flame, Part 2: flame stabilization. Energy Fuels 2012;26:4830–9. [16] Knudsen E, Pitsch H. Capabilities and limitations of multi-regime flamelet combustion models. Combust Flame 2012;159:242–64. [17] Mittal V, Cook DJ, Pitsch H. An extended multi-regime flamelet model for IC engines. Combust Flame 2012;159:2767–76. [18] Pierce CD. Progress-variable approach for large-eddy simulation of turbulent combustion Diss. Stanford University; 2001. [19] Germano M, Piomelli U, Moin P, Cabot WH. A dynamic subgrid-scale eddy viscosity model. Phys Fluids 1991;3(7):1760–5. [20] Pierce CD, Moin P. A dynamic model for subgrid-scale variance and dissipation rate of a conserved scalar. Phys Fluids 1998;10(12):3041–4. [21] Luo K, Bai Y, Yang JS, Wang H, Zhou LX, Fan JR. A-priori validation of a secondorder moment combustion model via DNS database. Int J Heat Mass Transfer 2015;86:415–25. [22] Launder BE, Reece GJ, Rodi W. Progress in the development of a Reynoldsstress turbulence closure. J Fluid Mech 1975;68:537–66. [23] Ghosal S, Lund T, Moin P, Akselvoll K. A dynamic localization model for largeeddy simulation of turbulent flows. J Fluid Mech 1995;286:229–55. [24] Desjardins O, Blanquart G, Balarac G, Pitsch H. High order conservative finite difference scheme for variable density low Mach number turbulent flows. J Comput Phys 2008;227:7125–59. [25] Wang H, Luo K, Fan JR. Direct numerical simulation and CMC (Conditional Moment Closure) sub-model validation of spray combustion. Energy 2012;46:606–17. [26] Wang H, Luo K, Fan JR. Effects of turbulent intensity and droplet diameter on spray combustion using direct numerical simulation. Fuel 2014;121:311–8. [27] Westbrook CK, Dryer FL. Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames. Combust Sci Technol 1981;27:31–43.