1 Introduction

Combustion plays a significant role in energy conversion nowadays. Its efficient utilization attracts much attention with the purpose to reduce energy consumption and protect the environment. Premixed combustion, especially the lean premixed mode, is widely used in gas turbine, because of its uniformly distributed temperature field and low emissions^[1]. Generally, ideal full premixed condition is difficult to achieve. In most cases, the reactants are non-uniform. Hence, stratified combustion is more common in pratical situations. Compared with premixed combustion, stratified combustion catches more attention lately, as it could extend the combustion limits and improve combustion stability^[2]. Large eddy simulation (LES) which decomposes the large scale structures and small scale motions through a filtering process has been widely used in turbulent flow and combustion simulations^[3-5].

Recently, some fundamental experiments for premixed and stratified flames have been carried out. Among them, the experiment carried on Cambridge stratified swirl burner designed and operated by Sweeney et al.^[6] is a well benchmark owing to its substantial data of the flow and scalar fields. In order to develop more suitable models for premixed and stratified combustion, massive efforts have been invested^[7-10]. It should be mentioned that for premixed or stratified combustion, the flame thickness is usually too thin to be resolved using an LES mesh scale. Hence, some sub-grid scale (SGS) combustion models should be proposed for these two combustion to handle the above problem. The flame surface density (FSD) model^[11], which describes the chemistry occurring in the flame front through a monotonically increasing characteristic scalar, namely, progress variable, is a reasonable choice for the premixed and stratified combustion simulations, and could address the flame propagation. Ma et al.^[12-13] performed detailed reviews for both the algebraic flame surface density model and the transported flame surface density model. Though the transported flame surface density model provides more detailed description of actual combustion process, the numerous unclosed terms in the equation prevent its application. Being simple and efficient, the algebraic flame surface density model has been employed widely^[14-16]. However, in the FSD model, the lack of the species information across the flame front holds back analyzing flame structures meticulously. The chemistry tabulation combined with presumed probability density function (PPDF) closure is a very efficient way for turbulent combustion simulations by decoupling flow and reaction, with an assumption that the compositions change through low-dimensional manifolds, but it fails to preserve the flame propagation property^[17]. Lecocq et al.^[18] took advantage of the two models to propose a new coupled model and applied it in the LES of a quasi-steady burner. The results verified that the new coupled model was capable of locating the reaction zone and predicting the distributions of species in the flame front. However, more works need to be conducted to improve the coupled model and to make it more suitable for premixed and stratified combustion.

In the present research, an SGS combustion model by coupling the FSD model and the FGM method is developed, and LESs are performed to investigate Cambridge stratified flames, SwB1 and SwB5. The characteristics of the improved coupled SGS model in premixed and stratified flames are evaluated versatilely. The present paper is composed of the following sections. In Section 2, brief introduction to the Cambridge stratified swirl burner and descriptions of the present used SGS models, including the FSD model, the FSD-FGM model, and the PPDF-FGM model, are presented. In Section 3, the numerical methods are described in detail. The numerical results and flame structures are discussed in Section 4. Finally, the paper ends up with conclusions in Section 5.

2 Experimental setup and physical model 2.1 Cambridge stratified swirl burner

The geometry of the Cambridge stratified swirl burner is shown in Fig. 1. It consists of two coaxial round pipes and a central bluff body. The inner circular passage supplies the CH₄/air jet under the normal pressure and temperature with the bulk axial velocity U_j. The outer one also supplies the CH₄/air swirling flow at the same conditions, with the bulk axial velocity U_s and the bulk circumferential velocity W_s. The experiments are carried out in a wind tunnel supplying the air at the room temperature and atmospheric pressure with a diameter of 196.8 mm. Its bulk axial velocity is U_e. In the present research, the non-swirling premixed flame SwB1 and the non-swirling stratified flame SwB5 are studied by the LES. The operating conditions are given in Table 1, and ϕ_j, ϕ_s, and ϕ_e represent the equivalence ratios of three inlets, respectively.

According to Zhou et al.^[19], the studied flames lie in the thin reaction zone (TRZ) regime on the modified Borghi diagram. Hence, the present SGS combustion model based on the FSD model and the FGM method are applicable.

2.2 FSD model

A normalized quantity called progress variable is required for the FSD model and generally expressed using the mass fraction of fuel Y_F as

(1)

where the superscripts u and b indicate unburned and burnt states, respectively. Z is the mixture fraction representing the distribution of reactants.

With low Mach number and unity Lewis number assumptions, the Favre filtered transport equation for the progress variable can be simplified as^[15]

(2)

Notably, the last term on the right-hand side of Eq. (2) is the cross scalar dissipation rate, which is caused by the spatial inhomogeneity for fuel. As its order of magnitude is much smaller than that of the reaction source term , the last term can be expressed using the quantities at the resolved scale as^[15]

(3)

where is the turbulent diffusion coefficient and calculated by Germano dynamic process^[20]. The filtered reaction source term is modeled by the following FSD model.

2.2.1 Modeling strategy using Charlette2 wrinkling factor model

The diffusion term and reaction source term in Eq. (2) can be integrated as a function of general flame surface density according to the methodology proposed by Boger et al.^[11],

(4)

where is the Favre surface averaged flame displacement speed and approximated as , in which ρ₀ is the unburned gas density, and S_L is the unstretched laminar flame speed. A modified laminar flame speed can also be used to consider the effect of curvature^[21]. Σ_gen is the generalised flame surface density and can be expressed as^[12]

(5)

where c is the Reynolds averaged progress variable. Ξ is the wrinkling factor, which represents the effect of turbulent eddies on the flame surface area density. In this paper, the Charlette2 wrinkling factor model is chosen and expressed as

(6)

Each parameter in Eq. (6) was determined by Ma et al.^[12].

By using the Weller model to compute the SGS scalar flux term in Eq. (2), the diffusion, filtered reaction source, and turbulent transport terms can be combined as^[16]

(7)

Then, the Favre averaged transport equation for the progress variable can be expressed as

(8)

2.2.2 Modeling strategy using Muppala wrinkling factor model

Moreover, the reaction source term in Eq. (2) itself can be approximated as a function of the generalised flame surface density as proposed by Muppala et al.^[22] and Ma et al.^[12],

(9)

where is the Favre averaged progress variable. In the present study, a wrinkling factor model based on the analysis from the experimental data, namely, the Muppala model, is also used. Its expression is shown as

(10)

where Le is the Lewis number and set to be 1.0. p, p₀, Δ, and ν are the filtered pressure, the reference pressure, the mesh scale, and the kinematic viscosity, respectively. The SGS fluctuating velocity u'_Δ is modeled using the Smagorinsky model with the model coefficient determined by the Germano dynamic process.

According to the above formulae, the Favre averaged transport equation for the progress variable is also expressed as

(11)

Compared with Eq. (8), the diffusion term is still retained in Eq. (11). As mentioned before, the FSD model needs to be coupled with some specific approaches to determine the detailed composition structures across the flame front.

2.3 FGM method and PPDF model

In the FGM method, the laminar flame equations are solved, and a small number of tabulated scalars are employed to describe the detailed chemical reaction processes based on the method proposed by Oijen and Goey^[23]. In this paper, a series of one-dimensional freely propagating laminar premixed flames are calculated by the code Flame Master to obtain the original flame database. The lean and rich flammability of equivalence ratio is set to be 0.4 and 1.8, respectively. The laminar FGM table can be tabulated based on mixture fraction and progress variable, namely, φ = φ(Z, c)(c = (Y_c - Y_cmin) /(Y_ceq - Y_cmin), where Y_ceq and Y_cmin are the maximum and minimum values for the unnormalized progress variable Y_c under a given mixture fraction, respectively).

Two expressions for the progress variables c₁ and c₂ are discussed in restoring scalars under the lean (ϕ = 0.5) and stoichiometric ratio (ϕ= 1.0) conditions covered by the flames SwB1 and SwB5. c₁ and c₂ are obtained from Y_c1 = Y_CH4 and Y_c2 = Y_CO + Y_CO2 + Y_H2O + Y_H2, respectively. Obviously, Y_c1 is identified with Y_F. Therefore, the expression for c₁ is the same as Eq. (1).

Figure 2 shows the comparisons of temperature and mass fraction of CH₄ and CO₂ versus different progress variables. Compared with c₂, the gradient of scalars obtained using c₁ is larger near c=1, especially for the stoichiometric ratio condition, which may cause a slightly larger error for the interpolation during the lookup procedure. It needs to be emphasized that there are more unclosed terms in the transport equation for c₂ which are hardly modeled. Therefore, c₁ is chosen for the tabulation to keep pace with the progress variable used in the FSD model.

The PPDF is employed to model the SGS flame structures as follows:

(12)

where (Z, c) is the filtered joint probability density function (PDF) of mixture fraction and progress variable. Assume that the two scalars are statistically independent. The joint PDF can be approximated as the product of the marginal PDF of those two scalars, i.e., (Z, c) =(Z) × (c).

In this paper, the PDFs of mixture fraction and progress variable are both modeled by the beta PDF^[5]. Then, the turbulent FGM table can be expressed as .

2.4 Combination of FSD and FGM models

In the FSD-FGM model, the turbulent scalars are determined from the turbulent FGM table by solving the transport equations for each characteristic scalar. Note that the transport equation for progress variable is expressed as either Eq. (8) or Eq. (11) based on the FSD model concept. The unstretched laminar flame speed S_L is also determined from the turbulent FGM table. To preserve the flame propagation characteristic, S_L is un-integrated and not affected by the PDF. According to the method proposed by Domingo et al.^[24], the other three related governing equations are expressed as follows:

(13)

(14)

(15)

The eddy diffusivity model is adopted for the SGS scalar fluxes , with the model coefficients determined by the dynamic procedure. The SGS scalar dissipation rate of mixture fraction is modeled using a linear relaxation hypothesis, i.e., . The SGS scalar dissipation rate of progress variable is modeled using a coupled approach proposed by Domingo et al.^[25], i.e., , where Sc is the normalized progress variable variance, and the parameters are set as C_χψ = 2.0, C_ε /C_u= 2.0, and S_ct=0.4. The difference between the present FSD-FGM model and the PPDF-FGM model is the approach to consider the turbulence/flame interaction, i.e., the source term in Eq. (2). The detailed description of the PPDF-FGM model can be found in our previous work^[10]. The reaction source term in the transport equation for progress variable variance is pre-computed and stored in the turbulent FGM table. The filtered density and diffusion coefficient are also determined from the turbulent FGM table.

In the present study, the widely used GRI3.0 mechanism (53 species, 325 reactions) is employed for tabulation. The turbulent FGM table contains 150, 25, 150, and 25 grid points in the corresponding coordinate direction. The distributions of grid points along each direction are referred to Zhang et al.^[10].

3 Numerical methods

The LES calculations are performed using an in-house FORTRAN code. The Favre filtered governing equations are discretized using the finite difference method in a cylindrical coordinate system. The detailed numerical schemes are described by Zhang et al.^[1].

The mesh and boundary conditions are shown schematically in Fig. 3. The computational domain is a cylinder of length L_x=240 mm and radius L_r=80 mm. The mesh used has 323 and 256 non-uniform grid points in the axial and radial directions, respectively, and 64 uniform nodes in the circumferential direction. The mesh is refined both at the exit and in the shear layer region. The total number of grid points is almost 5.29 million, which is reasonable for this LES computation, supported by our previous work^[10]. The no-slip boundary condition is adopted on the wall, and the convective boundary condition is implemented at the outlet boundary. The instantaneous inlet velocity is determined by the mean velocity superimposing the white noise. The mean inlet velocity profiles for the central jet and swirling flow are appointed as 1/7 power-law distributions. In coflow, the mean inlet velocity profile is set as the experimentally reported bulk-flow velocity. The time step is adjusted dynamically to guarantee the Courant number to be less than 0.5. The statistic ensembling average lasts for 10τ (τ=L_x/U_j, where L_x=240 mm and U_j=8.31 m/s) after running over 10τ to ensure the effectiveness of statistical results.

4 Results and discussion 4.1 Statistical results

Figures 4 and 5 show the predicted radial distributions of mean and root mean square temperature at different streamwise locations for the two flames, SwB1 and SwB5. The results obtained by using the FSD-FGM model with two different wrinkling factor models are compared with the experiments and the previous results obtained with the PPDF-FGM model^[10]. Firstly, from Fig. 4, for the premixed flame SwB1, almost even better results can be observed for the FSD-FGM model compared with the PPDF-FGM model. At x=10 mm, the predictions of temperature using different models are similar. At the downstream locations, the FSD-FGM model using both the wrinkling factor models shows the ability in describing the flame propagation. Secondly, from Fig. 5, for the stratified flame, much better results are achieved using the FSD-FGM model. At x=10 mm, the calculated temperature is overestimated due to the ignorance of heat loss on the bluff body. The discrepancies at x=50 mm reflect the inaccurate prediction for the location of flame front. The larger deviations using the PPDF-FGM model are probably related to the inappropriate PPDF profile for the stratified flame. Similar to the previous work^[8], the predictions of root mean square fluctuations of temperature by all the SGS models are less accurate which are likely to be caused by the large filter size in the downstream region and the incomplete boundary conditions in the upstream.

The ability of the FSD-FGM model to restore the distributions of major quantities is assessed. Figures 6 and 7 show the radial distributions of mean mixture fraction and mass fraction of major species for flames SwB1 and SwB5, respectively, using the FSD-FGM model with different wrinkling factor models. For the premixed flame SwB1, at locations x=10 mm and x=30 mm, the predicted results are in good agreement with the experiment data. The deviations at x=50 mm may be due to the slight underestimation for the turbulent flame speed. For the stratified flame SwB5, at x=10 mm, the generally satisfied predicted results using two different wrinkling factor models are achieved. In the experimental results, there are three steps in the profile of mixture fraction. As discussed by Nambully et al.^[8], the first step at r < 6 mm is possibly caused by the effect of differential diffusion around the bluff-body, which is not captured by the current LES. At x=30 mm, the deviations for the distributions of major species reflect the inaccurate predictions for the turbulent flame front thickness. At x=50 mm, this effect becomes more obvious. Generally, the predicted results show that the SGS models are suitable for the calculation of upstream, while need further improvement in downstream. Compared with the Muppala wrinkling factor model, the Charlette2 wrinkling factor model is a better choice to consider the flame surface area change caused by the motions of turbulent eddies.

4.2 Flow field structures

In the following sections, the results using the data obtained by the Charlette2 wrinkling factor model are given and discussed.

Figures 8 and 9 show the radial distributions of mean and root mean square values of axial velocity at different streamwise locations for SwB1 and SwB5. In general, the LES results are in good agreement with the experiments. At x=10 mm, the velocity in the recirculation zone is underestimated for both flames, due to the ignorance of heat loss. No recirculation zone is observed at x=30 mm and x=50 mm. Furthermore, the peaks of root mean square values of axial velocity locate in the shear layers.

The Instantaneous axial velocity fields for premixed and stratified flames are shown in Fig. 10. The recirculation zone formed behind the bluff body is closed to be stable owing to the re-laminarisation caused by heat release of combustion. Through momentum exchange and thermal expansion, both recirculation zones are dissipated gradually at downstream. For stratified combustion, a relative shorter recirculation zone due to the higher heat release can be observed.

4.3 Combustion characteristics

Figure 11 shows the snapshots of instantaneous mixture fraction and temperature for both flames SwB1 and SwB5. Seen in the figure, the effect of turbulence on flame is strengthened at downstream locations, where the flame front wrinkles more obviously. Besides, the flame brush of the stratified combustion SwB5 is wider than the premixed combustion SwB1, resulting from the higher heat release rate. The intersections of flame and mixing layer are also observed after x=30 mm for the stratified flame SwB5. This phenomenon is illustrated by a scatter diagram shown in Fig. 12. The expression is used to reflect the intersections. Seen from the figure, the values are almost equal to zero at both x=10 mm and x=30 mm, while nonzero values are predominant at x=50 mm, which indicates that the intersections mainly happen in the downstream region.

Compared with premixed combustion, stratified combustion is more complicated. It can be regarded as a special kind of partial premixed combustion, and the flame front travels along the mixtures of a series of equivalent ratios^[26]. The combustion mechanisms for stratified flames can be deduced by a flame index^[7],

(16)

The positive value of flame index represents the premixed mode, the negative value represents the non-premixed mode, and the zero indicates the unburned state.

Figure 13 shows the snapshots of instantaneous and mean flame index for the flame SwB5. The black long dash lines represent the isolines of normalized source term equaling 0.1 and 0.9, and the black solid line represents the isoline of c=0.5. Both are used to identify the reaction zone. Seen from the figure, the value of flame index in the reaction zone is positive, which demonstrates that the premixed mode is dominant.

In this paper, the extent of stratification is quantified using the probability distributions of equivalence ratio. Figure 14(a) shows the probability distributions of the instantaneous equivalence ratio p(ϕ) at different axial locations for the flame SwB5. The flame zone is assigned as 0.01 ≤ ≤0.99 to include enough valid samples. As shown in the figure, at x=30 mm, the maximum of p(ϕ) is located at ϕ = 1.0, while there also exists small amount of mixtures at ϕ < 1.0 in the flame zone. It can be inferred that the combustion mechanism is in a transition stage shifting from the premixed mode to the stratified mode. At x=40 mm and x= 50 mm, the maximum locations move towards the lean side. Meanwhile, p (ϕ = 1.0) approaches or equals zero. Therefore, one can conclude that the stratified mode is dominant at the two sections.

In order to achieve a deeper understanding of the structures of stratified flames, a quantity named the orientation angle is used to distinguish the combustion modes. It is defined as the angles between the gradient of progress variable and that of the mixture fraction^[8],

(17)

Generally, the local flame is in the back-supported mode when θ < 90°, the front-supported when θ > 90°, and the premixed mode when θ = 90°. Same as above, the flame zone is set as 0.01 ≤ ≤ 0.99. One can see from Fig. 14(b) that θ ≤ 90° makes sense at all three sections, which indicates that the leading mode is the premixed mode and the back-supported mode in the flame zone. As shown in the figure, at x=30 mm, the maximum of p (θ) is located at θ= 0° and p (θ =90°) ≠ 0, which means that the local flame is in the stage of transition. Further, downstream at x=40 mm and x=50 mm, the maximum locations are around θ = 30°. Meanwhile, p (θ =90°) = 0. Hence, it can be concluded that the flame at downstream locations is back-supported. The conclusions drawn from the probability distributions of equivalence ratio and orientation angle are in agreement with the previous work^[7-8], which indicates that the FSD-FGM model can predict the flame structures to a good extent.

5 Conclusions

LESs are performed to investigate Cambridge premixed and stratified flames SwB1 and SwB5 using an improved coupled SGS model. In the coupled SGS model, the FSD model with two different wrinkling factor models is used to describe the combustion/turbulence interaction, and the FGM method is employed to determine the major scalars. Some conclusions can be obtained as follows.

(ⅰ) Generally, the predicted results using different SGS models are similar in the premixed flame. The FSD-FGM model is superior to the PPDF-FGM model for the stratified flame in describing the flame propagation characteristics. The Charlette2 wrinkling factor model is a better choice to consider the flame surface area change caused by turbulent eddies. The discrepancies observed in the stratified flame using the PPDF-FGM model are probably caused by the failure of the PPDF to represent the turbulence/flame interaction in the stratified condition.

(ⅱ) For both premixed and stratified flames, the recirculation zones formed behind the bluff body are rather stable. The effect of turbulence on flame is strengthened at downstream, which makes the flame front wrinkle more obviously. The flame brush of stratified combustion is wider than the premixed combustion, resulting from the higher heat release rate. The intersections for the flame and mixing layer are also observed after x=30 mm for the stratified flame.

(ⅲ) The snapshots of instantaneous and mean flame index for the stratified flame demonstrate that the combustion mechanism is dominant by the premixed combustion. Judged from the probability distributions of equivalence ratio and orientation angle, the local stratified flame belongs to the premixed mode at upstream and back-supported mode at downstream.

In conclusion, the analyses of flow field structures and combustion characteristics indicate that the correct descriptions for the investigated premixed and stratified flames can be achieved by the FSD-FGM model.

Acknowledgements

All the numerical simulations have been done on the supercomputing system in the Supercomputing Center of the University of Science and Technology of China.