1 Introduction

Turbulent combustion is common in many combustion devices such as gas turbine combustors. The accurate prediction of species and other thermo-physical parameters is crucial for the designs of these combustion devices. Due to the complexities of the intrinsic physical and geometric configurations, the accurate measurements for the detailed flow information are extremely difficult and expensive. However, computational tools can provide designers with abundant flow details, playing an important role in the design process.

Owing to the existence of a large number of species, chemical stiffness, and scale separation, a direct numerical solution (DNS) of combustion systems becomes prohibitive for practical applications. The vast majority of computational modeling work so far has relied on simplified or reduced chemical mechanisms. Among various techniques in chemistry reduction, the flamelet generated manifold (FGM)^[1] approach allows the computation of the chemistry to be performed independently of the flow simulation and stored in a tabulated database as a function of a limited number of scalars. The mixture fraction (Z), its variance (Z^"2), as well as the progress variable and its variance are typical control variables in the FGM approach. During actual simulations, the quantities of interest are read and interpolated from the database, thus dramatically decreasing the computational cost and allowing the use of detailed chemical mechanisms. The FGM approach has been successfully applied to large eddy simulation (LES) studies on gaseous flames^[2-4] and spray flames^[5-8].

Although FGM combustion models are available in commercial computational fluid dynamics (CFD) softwares, they have different limitations. For example, the FLUENT software cannot deal with unsteady flamelet combustion and the effect of heat loss resulting from evaporation and radiation. In this paper, an FGM model capable to overcome the limitations of the commercial FGM models is implemented with the OpenFOAM solver^[9]. The model uses the mixture fraction, the progress variable, and their variances to parameterize the flame behavior, and the resulting database is accordingly referred to as four-dimensional (4D) flamelet manifolds. The model is tested against an LES of a piloted methane/air jet diffusion flame (Sandia flame D^[10]). The numerical predictions are compared with the experimental results available in the literature.

Besides, since the formation of NO is a relatively slow chemical process, which cannot be mapped onto the same manifold, it therefore needs a special treatment. In this work, we solve a transport equation for NO, in which the source term is determined from the look-up database. Moreover, to extend the applicability of the developed solver for its future application in gas turbine combustors, we introduce here a different formulation of the FGM approach, where the temperature is not obtained from the look-up table, but is calculated from the total energy solved by a transport equation and the species mass fractions, in order to allow for the heat loss effect due to spray evaporation and radiation.

The paper is structured in the following way. In Section 2, The LES numerical formulation is first presented, which is followed by the FGM-CFD coupling technique, the solution algorithm, and the numerical schemes. In Section 3, the Sandia flame D configuration and the corresponding FGM database are described. In Section 4, the numerical results are analyzed and compared with the experimental data in terms of the flow field and species profiles and the improved prediction of NO. The capability to handle the effect of heat loss is demonstrated. Finally, the conclusions are summarized in Section 5.

2 Numerical formulation

In this section, the compressible transport equations are first presented, and then the FGM-CFD coupling approach, the solution algorithm, and the numerical schemes are introduced.

2.1 Compressible transport equations

In the present work, the reactive flow is simulated using the LES approach. In the LES, the turbulent flow is decomposed into coherent large scales, which are directly resolved, and small unresolved scales. Applying a filter operation to the governing equations, the filtered governing equations for mass (1) and momentum (2) read

(1)

(2)

where ρ, u_i, and p denote the density, the velocity, and the pressure, respectively. The superscript is used for the filtering operation, and the superscript is used for the Favre filtering operation. The dynamic viscosity μ is a function of temperature according to Sutherland's three-coefficient law. The subgrid-scale viscosity μ_t is modeled with the dynamic Smagorinsky model^[11]. The rate of strain τ_ij is defined by

(3)

The mixture fraction Z and the progress variable are used to describe the chemical state of the the mixture. The transport equations for Z and can be expressed as follows:

(4)

(5)

where D_Z and D_Y denote the molecular diffusion coefficients. By assuming the Lewis numbers of unity for all species, ρD_Z and ρD_Y in Eq. (4) are the same, and are equal to

(6)

where λ is the thermal conductivity, and c_p the specific heat. The turbulent diffusivity D_t is calculated according to

where Sc_t is the turbulent Schmidt number and has a constant value here, i.e., 0.4.

Here, the progress variable is a linear combination of the reaction product species as follows:

(7)

where Y_i and M_i represent the mass fraction and the molar mass of species i, respectively. Its definition is chosen to keep a monotonically increase in both steady and unsteady flamelet solutions. The relevant progress variable source term is a linear combination of the chemical source term for the mass fraction Y_i of species i.

The element-averaged mixture fraction, i.e., the Bilger mixture fraction^[12], is adopted:

(8)

where the subscripts C and H correspond to the elements C and H, respectively, and M refers to the corresponding molar mass. In a mixing system of fuel and oxidizer streams, the values of the mixture fractions based on different elements may be different because of differential diffusion effects. This formulation preserves the stoichiometric value of the mixture fraction.

2.2 FGM-CFD coupling

Turbulent combustion is modeled with the FGM approach based on the assumption that a high-dimensional state space of chemical species can be approximated by a low-dimensional manifold (LDM). How the FGM model couples with the online CFD solver is illustrated in Fig. 1, where PDF is the abbreviation of probability density function. In the literature, the LDM is usually obtained from canonical cases, i.e., homogeneous reactors and one-dimensional (1D) laminar premixed flame and laminar counter-flow diffusion (CD) flame depending on the simulated combustion modes. Therefore, the first critical question is to choose a canonical test case which has a similar chemistry to our CFD simulation illustrated as the "raw data" part in Fig. 1. Since the Sandia flame D is a diffusion flame, steady and unsteady laminar counter-flow diffusion flames are used to create the LDM in this work.

The set of 1D equations for the CD flames, which are also referred to as flamelet equations, are^[13]

(9)

where m is the mass burning rate. The subscript i stands for the ith specie. Y is the mass fraction of the species, and h is the enthalpy. All the equations are expressed in terms of the arc-length normal to the flame front x. λ is the mixture conductivity, and c_p is the mixture heat capacity at constant pressure. Le is the Lewis number. In the present work, Le is equal to 1, which is suitable for most species in our simulations. is the reaction rate, and N_s is the number of the involved species.

The boundary conditions for the nonpremixed flamelets are the Dirichlet boundary conditions. Fuel enters at one side,

(10)

and oxidizer enters at the other side,

(11)

In the present work, Eq. 9 is solved with the Ember code. First, different steady flamelet equations are calculated with the scalar dissipation rate increasing from a very small value to the extinguishing value. Then, the extinguishing unsteady flamelets are calculated. Now, we have all species, enthalpy, and production rates in the physical space, i.e., the (x, t)-space. The state space is then parameterized by the mixture fraction and the progress variable to obtain an LDM, which is the "coordinate transformation" part as shown in Fig. 1. Figure 2 shows the LDM space, where the red lines stand for steady flames, the blue lines are for extinguish flames, and the black line is for the pure mixed solutions. The unsteady solutions allow the time-accurate transition from ignition to steady state to be captured. Now, we have a 2D FGM database to store the thermal physical properties for laminar flames, i.e., 2D FGM as shown in Fig. 1.

The influence of turbulence fluctuations on the local flame structure is accounted through the joint PDF for the mixture fraction and reaction progress variable. A widely used approach is based on the presumed joint PDFs. This process is illustrated as the "β-PDF integration" in Fig. 1. In the LES, all thermal physical properties are formulated as the Favre-filtered quantities:

(12)

where is the joint PDF, and β-PDF is used in this work. Since depends on Z, a normalized progress variable C is defined by

(13)

where denote the minimum and maximum values of at a given Z, respectively. The definition of C allows the assumption of statistical independence for Z and C^[14]:

(14)

The β-PDF equation for a scalar ϕ is

(15)

where Γ is the Γ-function, and a and b are the β-PDF parameters defining its shape depending on the mean value and its variance.

(16)

The subgrid mixture fraction and progress variable variance are calculated by

(17)

where α=1 is a reasonable choice^[4], and Δ denotes the LES filter width. An external numerical integration routines, i.e., Cuhre, in the CUBA library^[15] is used for the 2D integrals of Eq. (12). When is zero, Eq. (12) degenerates to 1D integration and the gsl_integration_qags routine in the gsl library is used. To avoid modeling the source term in Eq. (5), a transport equation for the unscaled progress variable is solved and C is only used for the value lookup from the FGM database. Finally, we obtain a 4D FGM for the turbulent combustion.

The coupling of the FGM database with the fluid solver is achieved via a lookup process as shown in Fig. 1. The species mass fractions, sensible enthalpy, and progress variable source term are stored as a function of and in the FGM database. During a single Pimple loop, the LES solver send these variables (old time state) to the lookup process subroutine firstly to get the source term in Eq. (5). Then, Eqs. (4) and (5) are solved to get the field variables of Z and at a new time step. After that, the subgrid variable variance is updated as shown in Eqs. (13) and (17). These variables at a new time step are sent to the lookup process to update the species mass fraction and enthalpy. Finally, the thermophysical and transport properties, i.e., species and transport coefficients, are subsequently updated in the LES solver. We note that since the specie mass fraction and enthalpy are looked up from the FGM database during the simulations instead of the solutions of extra transport equations, the computational cost decreases significantly.

2.3 Solution algorithm and numerical schemes

The adopted LES flow solver is based on a pressure-based compressible solver from the open-source OpenFOAM^[10]. The pressure-velocity-density solution algorithm is based on the PIMPLE algorithm^[16], a combination of PISO and SIMPLE. The spatial accuracy of the solver is formally second-order, and an implicit second-order time integration is used. The detailed information on the flow solver can refer to Ref. [17].

3 Case description

The validation case is a piloted methane jet flame (Sandia flame D)^[10-17]. The fuel is a mixture of 25% methane and 75% air by volume. The fuel nozzle has a diameter of D = 7.2 mm, and is enclosed by a pilot flow and an outer air co-flow. The pilot flow is a pre-burned mixture, which has been adjusted to match the species composition and temperature of the main fuel/oxidizer flow with a mixture fraction Z = 0.271. The fuel bulk velocity is 49.6 m/s, which leads to a fuel stream based Reynolds number of Re=22 400. The bulk inflow velocities of the pilot nozzle are 11.4 m/s, and the air coflow velocity is 0.9 m/s.

The computational domain is a cylinder with a diameter of 40D (see Fig. 3). The streamwise coordinate x is defined to be zero at the nozzle exit, and the outflow flow boundary is imposed at x = 80D. The grid clusters around the cylinder axis, and has a resolution of 164 ×120 × 300 and a minimum grid spacing of D/60 in the x-direction and D/80 in the y- and z-directions. This configuration is the result of a grid convergence analysis performed on a non-reactive flow field. A convective outflow boundary is applied with zero normal derivatives for all variables except for the pressure. The pressure is equal to the ambient pressure on the lateral boundary, and a wave transmissive pressure boundary condition is applied to the outlet boundary in the streamwise direction. The inflow conditions are critical to reproduce the experimental results. All the required mean quantities are specified according to the experimental data given in Ref. [10]. Perturbations are added on the prescribed mean velocity profiles via a vortex method^[18] with a decaying turbulence inflow generator implemented in the LEMOS extensions.

The simulations start from an initial state, where all velocities and scalars are initialized to zero, except that the streamwise velocity is initialized to the coflow velocity of 0.9 m/s and the temperature is set to the ambient condition. The Courant-Friedrichs-Lewy (CFL) number is set to be 0.8. The time-averaging statistics starts at 0.11 s, and the simulations run until approximately 0.5 s.

The chemistry data in the FGM table is generated with steady and unsteady laminar counterflow diffusion flames. The steady flames at a broad range of strain rates (0.7 < a < 390) are included. Additionally, the extinction flames are solved to improve the interpolation precision. All laminar CD-flame simulations are performed with the GRI 3.0 mechanism. A linear interpolation scheme is used to fill the manifold in the (Z, C)-space. The resolution of the mainfolds is 161 points for both variables, equally spaced between 0 and 1. Figure 4 shows the 2D laminar manifolds for the source term and the temperature T. Bilinear interpolation is used to retrieve the necessary quantities from the database.

The subgrid chemistry effects are modeled with the presumed β-PDF approach. In this approach, the 2D database is extended by another two variance entries, i.e.,

The reason why we use ζ and η instead of is to ensure that the γ-function in Eq. (15) is in the range of (0, +∞). The range of each variance is discretized with 21 points, equally distributed between 0 and 1. Finally, the FGM table resolution is set to Z × C × ζ × η = 161× 161 × 21 × 21 points with equidistant spacing in all directions. A linear interpolation is used when this table is accessed in the large-eddy simulations to return necessary values for the local values of the control variables.

4 Results and discussion

In this section, the developed solver is validated by comparing the obtained predictions with the experimental data. A quantitative comparison with the measurements is performed for the averaged and root mean square (RMS) profiles in Subsection 4.1. In Subsection 4.2, a transport equation for NO is solved to improve the prediction of NO. The ability of solving the non-adiabatic combustion system is demonstrated in Subsection 4.3 by solving a scalar transport equation for the absolute enthalpy.

4.1 Prediction of the flow field and species

The mean axial velocity and its RMS along the axial line are given in Fig. 5. The predictions are found to be sensitive to the imposed inlet velocity boundary condition. Both quantities are in good agreement with the experimental data.

The statistical results for the mixture fraction and progress variable are then analyzed. It is vital to predict these two quantities with accuracy, since they are highly related to the chemistry process. Figure 6 shows the mean and fluctuation results of the mixture fraction and the unscaled progress variable on the axial line. The mean quantities are well predicted by the current simulations until z/D=30. Thereafter, the experimental data are slightly underpredicted in the range of 30 < z/D < 60, which might be attributed to the interpolation error of the source term during the lookup process because of the sharp gradient (see Fig. 4(a)). The mixture fraction RMS is slightly overpredicted in the near field region. The underprediction when z/D < 10 is due to the imposed turbulent inlet boundary condition, which is not a fully developed turbulent condition.

The radial mixture fraction and unscaled progress variable profiles are shown in Fig. 7 for four different axial locations. The comparison with the experimental data shows that the general agreement for these variables is reasonable at all these positions. Small discrepancies of the mixture fraction at r/D=0 can also be observed very clearly in the axial profile in Fig. 6. The mixture fraction and unscaled progress variable for z=15D is slightly underpredicted in the range of 1 < r/D < 3, which is mainly attributed to the interpolation error of the source term.

Having verified that the mixture fraction and unscaled progress variable are reasonably predicted, we turn to the presentation of quantities, which are directly interpolated from the database. The calculated mean temperature and its variance along the axial line are shown in Fig. 8. Both the mean results and the RMS results agree well with the experimental data. Even the decrease in the temperature RMS around the maximum mean temperature, which can be clearly observed in the Fig. 8, is well represented by the simulation. In the fuel-rich part of the flame in the range of 20 < z/D < 30, the temperature RMS is slightly overpredicted.

The axial profiles for the most important chemical species in the flow are shown in Fig. 9. Generally, all profiles except the NO profile are reasonably well predicted. NO is predicted less satisfactorily, as shown in Fig. 9(h), which will be discussed in Subsection 4.2. Almost all discrepancies are located in the range of 30 < z/D < 65, where H₂O, CO, and H₂ are underpredicted, while OH is overpredicted. The reason for this is attributed to the underprediction of the unscaled progress variable and that different species have different sensitivities to the unscaled progress variable. CO₂ is underpredicted, and O₂ is overpredicted in the far field. The reason for this is that the manifolds of these two species are very sensitive for both low mixture fraction and high unscaled progress variables.

4.2 Prediction of NO

The formation of NO is a slow process compared with that of other major species. Therefore, it cannot be assumed in a quasi-steady state. Figure 10(a) shows the NO mass fraction mapped onto the space. It can be seen that the main variation of Y_NO occurs close to .

A scalar transport equation for the mass fraction of NO is solved to improve its prediction accuracy:

(18)

where ω_NO denotes the source term, which is retrieved from the manifolds as shown in Fig. 10(b). A significant improvement is obtained (see Fig. 9(h)). The small discrepancies can mainly be attributed to the underprediction of the NO source term. It is interesting to note that both the mass fraction of NO and the unscaled progress variable, whose governing equations have the source term, are underpredicted due to the interpolation error. The sharp gradient with respect to the mixture fraction and unscaled progress variable, similar to the progress variable source term manifold, is shown in Fig. 4(a).

4.3 Inclusion of the additional energy equation

The combustion in industrial combustion devices is usually accompanied with heat losses due to the evaporation of liquid fuel and radiation. To extend the capability of the FGM approach to handle the heat loss effect, a transport equation for absolute enthalpy is solved. The temperature is not looked up from the chemistry table, but is calculated from the absolute enthalpy and species mass fractions. This approach is suitable for the cases with an insignificant enthalpy deficit effect, and can deal with a slight compressible effect.

(19)

where h is the absolute enthalpy, and S_h is the source term introduced by evaporation and radiation. In the present study, the enthalpy deficit is ignored.

The predictions of the axial mean and fluctuation results of temperature obtained from the energy equation are shown in Fig. 11. Since no evaporation and radiation is included in the simulation, no significant discrepancy is observed.

5 Conclusions

In this paper, a turbulent jet diffusion flame (Sandia flame D) is studied numerically. The simulation is carried out with a customized OpenFoam solver. In this solver, the FGM method is formulated as a combustion model for the LES. A transport equation is solved to improve the prediction of NO. The source term is retrieved from manifolds by interpolation, and the ability to handle the enthalpy loss effect due to evaporation or radiation is also demonstrated.

The results are compared with the experimental data for the flow field and the species. The agreement is very reasonable for all quantities. The remaining differences are discussed. In the analysis of the computational results, it has been found that the interpolation errors due to the sharp gradients in the manifolds lead to lower accuracy of some components.