1 Introduction

Combustion is in a position of great importance for energy structures at the present time. Although the development of new energy conversion technology has achieved great progress in recent years, it can be inferred that combustion will remain irreplaceable for an extended period^[1]. Among different combustion classes, the premixed mode has been widely used in gas turbine and aero-engine, because of its uniform temperature field and low emission^[2]. However, higher demands for premixed combustor optimal design are required, as there usually exist some adverse phenomena especially for lean condition, such as backfire, blowoff, and thermal acoustic instability^[3].

Large eddy simulation (LES), which solves the large scale structures accurately and models small scale fluctuations using some specific methods through a filter operation, is a promising approach. However, a thorny issue is that the flame thickness in premixed combustion is usually too thin to be captured under a typical LES mesh density^[4]. The artificial thickened flame (ATF) model^[5], which increases the diffusivity and reduces the reaction source term accordingly by a thickening factor, can preserve the flame propagation characteristic, has been widely used for turbulent premixed flame simulations^[6-7]. Coupled with the flamelet generated manifold (FGM) method, the ATF-FGM model can resolve the turbulent premixed flame structures on an LES grid under a reasonable computational cost for considering the detailed chemical kinetic mechanisms^[8-9]. Recently, a new ATF-FGM method which includes the effects of relevant variances has been proposed by Zhang et al.^[10]. Its performances in the simulations of the Cambridge (SwB) flames are acceptable. One should note that the entire flow field is artificially thickened in the original ATF model so that the mixing description is modified seriously by changing the molecular diffusion coefficients. Therefore, a quantity to mark the flame front, namely, the flame sensor, is required to avoid excessive thickening^[11]. However, the transformations of temperature and species are mainly controlled by the diffusion and convection processes in the pre-heat zone and the diffusion and reaction processes in the reaction zone^[12]. Therefore, for the thickened flame model, it is not appropriate to reflect the diffusion and reaction processes which have different spatial scales using a single flame sensor.

Most of the presented sub-grid scale (SGS) models are proposed with the help of equilibrium hypothesis between turbulence and reaction. The transport approach, which evades the equilibrium hypothesis by solving a transported equation, is closer to the real physics. However, it is barely used in the LES as there are so many unclosed items in the relevant transport equation. Compared with the transported model, the dynamic model can take the nonequilibrium effect into account more conveniently and efficiently^[13]. It determines the desired quantities dynamically through a "Germano-like" procedure, and has been popularized in recent studies^[14-16].

The turbulent Bunsen burner designed and experimented by Chen et al.^[17] is a well benchmark for testing new combustion models^[18-21]. In the present research, a new coupled ATF-FGM SGS combustion model is proposed and applied to the turbulent premixed flame. Its dynamic formulation is also tested. Two self-contained flame sensors are employed to track the diffusion and reaction processes in the flame front, respectively. The article is structured as follows. Brief introductions of laminar flame structures and experimental setups are shown in Section 2. In Section 3, detailed descriptions of the presently used SGS models are presented. The numerical methods and example settings are shown in Section 4. The numerical results and the comparisons of different SGS models are included in Section 5. Finally, the conclusions are highlighted in Section 6.

2 Laminar flame structures and experimental setups

The ATF model is put forward from the analysis of laminar flame structures. Therefore, the laminar flame structures are discussed firstly. We all know that combustion is a complex physical phenomenon, which contains many spatial and temporal scale processes. Typically, the laminar flame structures are mainly made up of a pre-heat zone and a reaction zone, as shown in Fig. 1. The flame front propagates from the burnt side to the unburned side at the laminar flame speed S_L. The transformations of temperature and species are mainly controlled by the diffusion and convection processes in the pre-heat zone and the diffusion and reaction processes in the reaction zone^[12]. Each process has its own spatial scale. Compared with the pre-heat zone, the reaction zone is much thinner.

The basic sketch of the Bunsen burner is shown in Fig. 2. The central jet is stoichiometric CH₄/air mixture with the mean velocity U_j, and the inner diameter is D = 12 mm. The pilot stream is the high temperature combustion products of stoichiometric CH₄/air mixture, the mean velocity is U_p = 1.5 m/s, and the inner diameter is d_p = 68 mm. The heat loss caused by the water cooling makes the temperature of the pilot stream lower than that of the adiabatic flame. To avoid the effects of the environment on the flame, the coflow is normal temperature air, and its velocity is U_co = 0.22 m/s.

In this paper, the low Reynolds number flame F3 is chosen as a target case. The specific experimental conditions are shown in Table 1.

3 Mathematical and physical models 3.1 FGM method and presumed probability density function (PPDF) model

In the FGM method, the detailed chemical reaction processes are reflected by a small number of tabulated scalars. The detailed derivation can be found in Oijen and Goey^[22]. The database for the chemistry tabulation is obtained by calculating a series of one-dimensional freely propagating premixed flames under different equivalent ratios using the code FlameMaster.

The tabulated scalars are generally set as the mixture fraction Z and the normalized progress variable c for the laminar FGM table, i.e., φ =φ (Z, c). The two scalars are defined below. According to the experimental conditions of the Bunsen burner, the equivalent ratio and mixture fraction are defined to ensure that for the central stoichiometric CH₄/air jet, Z=1, while for the air coflow, Z=0. The normalized progress variable is defined as c=(Y_c-Y_cmin)/(Y_ceq-Y_cmin), where Y_c is the progress variable defined as a linear combination of the main species, i.e., CO+CO₂+H₂O+H₂, and Y_ceq and Y_cmin are its maximum and minimum values under a given mixture fraction, respectively.

Modeling the SGS flame structures, the turbulent FGM table is obtained by the integration of the laminar FGM table as

(1)

where is the presumed joint probability density function (PDF) for mixture fraction and normalized progress variable. With an assumption that the normalized progress variable is statistical independent of the mixture fraction, the joint PDF can be simplified as the products of two relative marginal PDFs, .

The PDFs for mixture fraction and normalized progress variable are both presumed using the beta function. Therefore, the turbulent FGM table is expressed as from Eq. (1).

Notably, in order to avoid too many unclosed items in the transport equation, the progress variable is solved rather than the normalized one. Therefore, the lookup process needs to be transformed as^[1]

(2)

(3)

3.2 ATF-FGM model

In the ATF-FGM model, the flame is thickened to be resolved under the general LES mesh density on the condition that the laminar flame speed remains unchanged. The Favre filtered transport equation for the progress variable using the ATF-FGM model is expressed as

(4)

where Ξ, F, and Ω are the wrinkling factor, the thickening factor, and the flame sensor, respectively. The three parameters are discussed below. The filtered density ρ, the diffusion coefficient D, and the chemical source term ω_Yc are pre-computed and determined from the turbulent FGM table. It should be highlighted that the source terms in the turbulent FGM table are un-integrated and not affected by the PPDF, which is the same as the previous work^[7].

The wrinkling factor is used to compensate the reduction of flame front area change caused by the thickening process. A power law formulation proposed by Wang et al.^[14] based on the work of Charlette et al.^[6] is used,

(5)

where β is the model coefficient and needs modulation for different cases. The SGS velocity fluctuation u'_Δ can be expressed using the resolved velocity with a similarity assumption as

(6)

It can also be estimated using a simple model based on a shear stress related closure proposed by Lilly^[23] as

(7)

where ν_t is the turbulent viscosity modeled by the dynamic Smagorinsky model, and the model constant C_L = 0.094. The meanings of the rest quantities in Eq. (5) can be referred to Zhang et al.^[10].

The flame sensor Ω is defined to represent the flame region, where the artificially thickening procedure is carried out. The thickening factor F is expressed as a function of flame sensor Ω,

(8)

(9)

where n is the number of grid point to resolve the flame structure, which is set to be 5 according to Charlette et al.^[6]. Δ_mesh is the mesh scale. δ_l⁰ is the laminar flame thermal thickness, and , where T_b and T_u are the burnt and unburned temperature, respectively.

In this paper, two independent flame sensors are used to represent the diffusion and reaction processes, respectively. Same with above, the flame sensor Ω₁, which represents the diffusion process, is determined as a function of gradient of progress variable^[9],

(10)

where the subscript 1-D stands for the quantity determined from the turbulent FGM table.

The reaction process is represented using the flame sensor Ω₂ as^[11]

(11)

where η is a parameter controlling the thickness of the thickened zone, and the hyperbolic tangent function is used to make the flame sensor characterize the source term more reasonably^[11]. The distributions of different flame sensors for one-dimensional laminar CH₄/air flames under two equivalence ratios are shown in Fig. 3, where ω^* is the normalized source term (ω_Yc(x)/max(ω_Yc(x))). Gradient represents the flame sensor calculated by Eq. (10), Tanh3.0, Tanh5.0, and Tanh10.0 represent the flame sensor calculated by Eq. (11) with η equal to 3.0, 5.0, and 10.0, respectively. It can be seen from Fig. 3 that the distribution of the source term is represented using Ω₂ logically. In this paper, η is set to be 3.0 to avoid excessive thickening.

In order to thicken the flame enough and meanwhile maintain the laminar flame speed, the thickening factor for progress variable is set as the maximum of F₁ and F₂^[7],

(12)

(13)

The distributions of the reaction source term in the turbulent FGM table colored by different flame sensors are shown in Fig. 4. As shown in the figure, the source term may not be thickened enough using the flame sensor Ω₁, while most of its higher values can be covered by the flame sensor Ω₃.

The transport equations for mixture fraction and its variance are not affected by the thickening process using double flame sensors, as there is no chemical source term and the diffusion process is predominant. According to the method proposed by Domingo et al.^[24], the related transport equations are shown below.

3.3 Dynamic formulation

In order to achieve a reasonable description for the turbulent flame structures at all stages, in this paper, a "Germano-like" dynamic procedure is used to calculate the model coefficient β in Eq. (5). The dynamic procedure depends on the comparisons of the reaction source term at different scales. The reaction source term can be written as a function of flame surface density. Assume that the reaction source term for the test filter progress variable equals the one at the test filter scale at a given range 〈·〉. Then, the "Germano-like" equation is obtained as

(14)

where the superscript represents the test filter operation with the test filter scale , c is the Reynolds averaged normalized progress variable, γ is the effective filter scale coefficient, and .

With an assumption that the wrinkling factor remains unchanged at the given range 〈·〉, the model coefficient β can be calculated as

(15)

In this paper, the filter scale for the ATF model is set to be Δ = 1.4F₃δ_l⁰, the test filter scale is set to be , and the given range is a cubic with a length of 4Δ. The test filter operation is realized using the second-order Gaussian filter as^[16]

(16)

where x represents space coordinates.

3.4 LES governing equations

In the present research, the Favre filtered Navier-Stokes equations are solved with the SGS stress modeling by the dynamic Smagorinsky model. The transport equations for Favre filtered tabulated scalars are also required for a complete description for the combustion field. Meanwhile, a transport equation for filtered total enthalpy , the sum of sensible enthalpy and chemical enthalpy, is also solved to consider the heat loss. The related equations are as follows:

(17)

(18)

(19)

(20)

(21)

It needs to be emphasized that there are a lot of unclosed items in Eqs. (17)-(21). The SGS scalar flux are modeled using the eddy diffusion model. The SGS scalar dissipation rates and are modeled through the linear relaxation hypothesis, i.e., , where C_χψ = 2.0, C_ε/C_u = 2.0, and Sc_t=0.4.

In this paper, the ATF-FGM model using double flame sensors is called as the double-flame-sensor ATF-FGM (DATF-FGM) model, and its dynamic formulation (i.e., the model coefficient is calculated using Eq. (15)) is called as the DDTF-FGM model. The ATF-FGM model using single flame sensor (i.e., replace F₃ and Ω₃ in Eqs. (19)-(20) with F₁ and Ω₁, respectively) is called as the single-flame-sensor ATF-FGM (SATF-FGM) model.

As for the PPDF-FGM model, the wrinkling factor Ξ, the thickening factor F, and the flame sensor Ω in Eqs. (17)-(21) are not required, and can be set as 1.0 for Ξ and F and 0 for Ω. The reaction source terms in Eqs. (20)-(21) are determined from the turbulent FGM table through the PPDF by Eq. (1).

Considering the effect of heat loss caused by the water cooling, the temperature is determined as

(22)

where T₀, , and are the reference temperature, the enthalpy of formation for mixture, and the constant-pressure specific heat for mixture, respectively, and the reference temperature is set to be T₀ = 298.5 K. The filtered density is determined through the ideal gas state equation,

(23)

where p, , and R₀ are the filtered pressure, the mixture molecular weight, and the universal gas constant, respectively. , and are determined from the turbulent FGM table.

The GRI3.0 chemical reaction mechanism (53 species and 325 elementary reactions) is chosen to construct the turbulent FGM table . The table contains 150, 25, 150, 25 grid points at the corresponding coordinate directions, respectively. The grid points distribute uniformly along the and , while parabolically at the rest directions.

4 Numerical methods and example settings

The LES calculations are performed using an in-house FORTRAN code. The Favre filtered governing equations are discretized using the finite difference method with staggered grids. The detailed numerical methods can be found in our previous work^[1]. The Courant number is guaranteed to be less than 0.5 using dynamically adjusted time step.

The cylindrical computational domain has 20D and 6D in length at axial and radial directions, respectively. There are 350 and 150 non-uniform grid points in axial and radial directions, and 32 uniform grid points in the circumferential direction, respectively. The mesh is refined around jet exit and shear layer region. The total number of grid points is almost 1.68 million, and the one is set as Mesh 1. According to the mesh density in Dodoulas and Navarro-Martinez^[18], the one in this study is enough for LES of Bunsen flame F3. A grid convergence study is also conducted. Mesh 2 has the same grid points in axial and radial directions and 64 uniform grid points in the circumferential direction. The no-slip boundary condition is chosen on the wall, and the convective boundary condition is used for the domain exit. The inlet velocity profile in the central jet is defined by the instantaneous velocity from a solution of pre-computed fully developed pipe flow using the bulk-flow velocity given by experiments and periodic boundary condition. The inlet velocity profiles in pilot stream and coflow are defined by the bulk-flow velocity reported from experiments. The enthalpy value in the central jet is set to be the one of stoichiometric CH₄/air mixture at T = 298.5 K, and the enthalpy value in coflow is set to be the one of normal temperature air. The distribution of the temperature in the pilot stream is not uniform and varies between 1 785 K and 2 248 K. In this paper, the temperature in the pilot stream is set to be T = 1 950 K, which is the same as that in Ref. [19]. Langella and Swaminathan^[20] indicated that the different inlet temperatures for the pilot stream only influence the temperature calculation at upstream and have no significant effect on the prediction of the overall temperature field distribution.

According to Wang et al.^[14] and Volpiani et al.^[15], when considering the DATF-FGM model and SATF-FGM model, the model coefficient β is set to be 0.25 for flame F3. For the DDTF-FGM model, the model coefficient β is determined dynamically. Table 2 makes a list of all simulation cases in the present study. Each case is running over 10τ (τ = L_x/U_j, where L_x = 20D) before beginning to take statistic ensemble average to ensure the statistical stationary, and the statistic process also lasts over 10τ. It should be mentioned that the DATF-FGM model is referred to Case 1 if there is no special declaration.

5 Results and discussion 5.1 Grid convergence study

In this section, the results of normalized temperature and axial velocity calculated on different meshes shown in Figs. 5 and 6 are compared to conduct the grid convergence study. The normalized temperature is defined as , where is the statistical filtered temperature. The normalized axial velocity is defined as , where is the statistical filtered axial velocity. One can see from the two pictures that the distributions of normalized temperature and velocity are captured well by both meshes and the two results are similar. Considering the large calculated quantities in combustion simulations, the following cases are carried on Mesh 1.

5.2 Statistical results based on static models

The results of normalized temperature, axial velocity, and radial velocity are given to show the differences between different static models. The normalized radial velocity is defined as , where is the statistical filtered radial velocity. Figure 7 shows the radial distributions of normalized temperature obtained using different static models. As shown in the figure, all static models fail to predict the distributions of normalized temperature at x = 2.5D, because of the effects of the inlet conditions for the pilot stream, which is similar with the results in Langella and Swaminathan^[20] and Kolla and Swaminathan^[21]. At downstream locations, the DATF-FGM model is able to predict the distributions of non-dimensional values of temperature more accurately compared with the other models. It can be deduced that the DATF-FGM model can preserve the premixed flame propagation well.

Figure 8 shows the radial distributions of normalized axial and radial velocity using different static models. As seen from Fig. 8, the DATF-FGM model shows ability in predicting the distributions of axial velocity reasonably for flame F3. The results are slightly underestimated at downstream, which may be due to the large dissipation caused by the present model there. For the radial velocity, compared with the PPDF-FGM model, relative good results are achieved using the DATF-FGM model. The predicted normalized radial velocity is slightly underestimated at downstream x = 8.5D. Probably, it also has something to do with the large dissipation caused by the present model.

As seen from Figs. 7 and 8, the results obtained by all ATF-FGM based models are similar. The main differences between Case 1 and Case 4, observed for the radial velocity at upstream, indicate that the flame sensor used in the DATF-FGM model can thicken the flame front more reasonably. Besides, the results of Case 1 are better than those of Case 5 because the SGS velocity fluctuation expressed using Eq. (6) can reduce the contribution of thermal expansion.

The distributions of mass fraction of CH₄, CO₂, and CO are chosen to evaluate the ability in restoring the species for the DATF-FGM model. Its performance is compared with the PPDF-FGM model. Figure 9 shows the radial distributions of mass fraction of major species obtained by Case 1 and Case 3, respectively. One can see that the predicted radial distributions of mass fraction of CH₄ obtained using the PPDF-FGM model are better at upstream x = 2.5D and x = 4.5D, while these are in reverse at the rest positions. The results of the radial distributions of mass fraction of CO₂ obtained using the PPDF-FGM are slightly better than those using the DATF-FGM model, which indicates that the more rigorous coupling method between the ATF model and the FGM method remains attention. Both models can hardly predict the distributions of mass fraction of CO accurately. At the upstream x = 2.5D, the DATF-FGM model shows a reasonable prediction for the peak value. The peak values are underestimated at the rest positions.

5.3 Combustion characteristics

In this section, the combustion characteristics are demonstrated using the Q criterion, the instantaneous and mean normalized progress variable, and the flame brush thickness. All the analyses are based on the results obtained from Case 1.

The Q criterion is expressed as the quadratic invariant of the velocity gradient tensor, of which the positive value can identify the vortex. Figure 10 shows the snapshots of the isosurface of Q (Q = 1.5 × 10⁷ s^-2) colored by temperature. As shown in Fig. 10, the large scale vortices are formed due to the Kelvin-Helmholtz instability at upstream, while the vortex breaking phenomenon is more obvious at downstream.

The snapshots of instantaneous normalized progress variable are shown in Fig. 11. The flame height is almost 14.3D observed from the mean snapshot, which is close to the curved-fitted one reported in the experiment^[17]. The quite high level flame wrinkling can be observed from the instantaneous snapshots, and the flame wrinkles more severely at the flame tip.

The mean flame brush thickness at different axial locations is shown in Fig. 12. Although the differences between the LES results and the experiments are obvious, the trend is consistent. As suggested by Chen et al.^[17], the turbulent flame brush thickness increases linearly along the axial direction determined by turbulence.

5.4 Analysis of the dynamic procedure

The differences between the dynamic model (DDTF-FGM) and the static model (DATF-FGM) are revealed by the distributions of normalized temperature and radial velocity shown in Figs. 13 and 14. As seen from the two pictures, both results are similar, due to the approximate steady-state for flame F3. The DDTF-FGM model has certain advantages in predicting the distributions of temperature in the flame front at downstream and the distributions of radial velocity at upstream, which indicates that the dynamic model can preserve the premixed flame propagation characteristics better. Compared with Case 1, the CPU time is almost 8.1% higher for dynamic formulation, which is an acceptable increase.

Figure 15 shows the snapshots of instantaneous source term obtained from Case 1 and Case 6. As seen in the figure, both results are similar, and the source term is larger at the more wrinkling zone. As discussed by Volpiani et al.^[15], the instantaneous source term length of Case 6 (the dynamic model) is shorter because of the higher wrinkling factor at downstream.

6 Conclusions

In this paper, a new ATF-FGM SGS model is proposed and applied to LESs of Bunsen flame F3. In the new SGS model, two flame sensors are used to trace the diffusion process and chemical reaction process with different space scales, respectively. The differences between the dynamic model (DDTF-FGM) and the static model (DATF-FGM) are also analyzed. The conclusions are as follows:

(ⅰ) The DATF-FGM model shows a marvelous ability in predicting the distributions of the temperature and velocity, while the distributions of the species need to be improved.

(ⅱ) The vortex breaking phenomenon is more obvious at downstream. The flame wrinkles remarkably and more severely at the flame tip for flame F3.

(ⅲ) The DDTF-FGM model has certain advantages in predicting the distributions of temperature in the flame front at downstream and the distributions of radial velocity at upstream with an acceptable increasing CPU time, which indicates that the dynamic model can preserve the premixed flame propagation characteristics better.

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