A Generalised Series Model for LES of Premixed and Non-Premixed Turbulent Combustion

1. Introduction

A challenge in the development of turbulent reacting models is the broad spectrum of length and temporal scales introduced by large Reynolds (Re) and Damköhler (Da) numbers in modern combustors [1,2,3]. Compromising between the practicality of Reynolds Averaged Navier-Stokes (RANS) for industrial scale applications and the accuracy of Direct Numerical Simulations (DNS), Large Eddy Simulation (LES) has proven itself a promising technique for the simulation of propulsion and energy systems, if appropriate sub-grid-scale (SGS) models are employed. The vast majority of established reactive SGS models are “borrowed” from those in RANS, and can be classified into three categories [2,4]: a) geometric approaches, such as flamelet [5], flame wrinkling [6], and thickened flame [7], b) statistical methods, including presumed [8] and transported PDF [9] and Conditional Moment Closures [10]/Multiple Mapping Closures [11], c) models based on turbulent mixing description, like, eddy break up [12], and eddy dissipation assumption [13]. Extensive reviews can be found in [14,15,16].

While these models offer satisfactory predictions, if combined with fine grid resolutions, they rely heavily on assumptions related to specific combustion regimes or they are too expensive (such as the PDF transport model). Formerly, mathematical approaches [2,17,18] based on series expansion were attempted within the RANS framework, and were limited to conditions like supersonic combustion [17] and atmospheric boundary layers [2,18]. In RANS, the high-order series of the chemical source may partake larger contributions beyond the first-order one and thus nonnegligible truncation errors, due to the substantial temporal oscillations of the temperature and specie fraction. This gives rise to the challenging task of modelling accurately the involved scalar gradient ($\overline{{Y_F}^{″}{Y_O}^{″}},\overline{{Y_F}^{″}{T}^{″}}$). Nevertheless, the same challenge does not exist in LES since filtering is operated in physical space (rather than temporal averaging of RANS), so that the higher-order contribution becomes less given a relatively well-resolved mesh and short time-stepping. This is evident in the research [19,20] of LES of non-reactive flows, in which a Taylor expansion on the filtered velocity field operates to close the SGS Reynolds stress in the Navier-Stokes equations. The simulation is validated against DNS and obtain satisfactory precision [21]. The model is later applied to modelling a neutral boundary layer flowing over a rough wall and the velocities coincide with theoretic solutions very well, ameliorating against standard eddy-viscosity approaches [19].

Similar mathematical methods [22,23] have also been proven feasible in the premixed turbulent combustion context, where series expansion is acted upon the reactive scalar field to determine the F factor in an artificially thickened flame (ATF) model. LES of a turbulent Bunsen flame is conducted, then accurate predictions regarding velocities and scalars are achieved [22,23]. For turbulent scalar signals, fourth order terms are found not to play a major role [22].

Based on the above work, we have presented a novel series model for turbulent combustion in the previous research [24]. Taylor series expansion is performed in scalar space to avoid the highly spatial non-linearity of the chemical source term; meanwhile, the method is generalised for non-premixed and premixed combustion, owing to its mathematical derivation. This model was previously validated on premixed bluff-body stablised flames [24]. Due to the experimental limitation of species information [25], we only briefly examined the predictive capability of combustion dynamics, velocity and temperature profiles (but not species) in the premixed context, in addition to the grid resolution dependence [24].

In this paper we build upon our previous work in order to provide a more detailed and well-rounded model evaluation. We extend the analysis of the model in the premixed regime with a new case (Bunsen jet flame F3 [26]) which benefits from a wider range of experimental data for validation. We also present results for non-premixed cases (Sandia Flame D [27]). The suitability of the model in non-premixed combustion has not been presented before. The objectives of this study are twofold: exploring the model’s generality for both premixed and non-premixed combustion simulation, and assessing the predicting accuracy for reactive species distribution simulations. In the next sections, the derivation of the series model and the numerical implementation into OpenFOAM is demonstrated, followed by the error analysis. Subsequently, LESs of the two widely-used validation benchmarks are presented. The results are analysed along with experimental data and data from previous relevant publications. Finally, major conclusions and suggestions for future work are presented.

2. Methodology

2.1. Mathematical formulation

The Mathematical formulation is identical to that in the study [24]. In the LES governing equations [2], the SGS stress tensor is computed by the dynamic one equation eddy model [28], although other approaches can be in principle be used. For the closure of reaction rates, a Taylor series expansion is operated in scalar domain, averting the high nonlinearity in physical space (regarding the source term as a function of spatial coordinates). The chemical source term is first expanded in a single scalar form.

$$\dot{\mathsf{\omega }}\left(\mathrm{c}\right)=\dot{\mathsf{\omega }}\left(\overline{\mathrm{c}}\right)+{\left.\frac{\partial \dot{\mathsf{\omega }}}{\partial \mathrm{c}}\right|}_{\mathrm{c}=\overline{\mathrm{c}}}\frac{\partial \mathrm{c}}{\partial {\mathrm{x}}_{\mathrm{i}}}\mathsf{\delta }{\mathrm{x}}_{\mathrm{i}}+\frac{1}{2}{\left.\frac{{\partial }^2\dot{\mathsf{\omega }}}{\partial {\mathrm{c}}^2}\right|}_{\mathrm{c}=\overline{\mathrm{c}}}{\left(\frac{\partial \mathrm{c}}{\partial {\mathrm{x}}_{\mathrm{i}}}\mathsf{\delta }{\mathrm{x}}_{\mathrm{i}}\right)}^2+\dots$$

(1)

Note that index notation is employed. Then an isotropic filter is practiced to relation (1). All terms with odd powers of x, y, and z vanish by means of symmetry elimination [24]. For detailed process, one can refer to the study [24] and the similar manipulations in the publications [19,20,21,22,23].

$$\overline{\dot{\mathsf{\omega }}\left(\mathrm{c}\right)}=\dot{\mathsf{\omega }}\left(\overline{\mathrm{c}}\right)+\frac{{\mathsf{\Delta }}^2}{24}{\left.\frac{{\partial }^2\dot{\mathsf{\omega }}}{\partial {\mathrm{c}}^2}\right|}_{\mathrm{c}=\overline{\mathrm{c}}}\overline{{\left(\frac{\partial \mathrm{c}}{\partial {\mathrm{x}}_{\mathrm{i}}}\right)}^2}+\mathrm{O}({\mathsf{\Delta }}^4)$$

(2)

The same processes are executed on the chemical source, a function of multiple scalars ${\mathsf{\phi }}_{\mathrm{k}}$, including chemical species, temperature and pressure:

$$\begin{gathered} \overline{\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{k}}\right)}=\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left(\overline{{\mathsf{\phi }}_1},\overline{{\mathsf{\phi }}_2},\dots ,\overline{{\mathsf{\phi }}_{\mathrm{k}}}\right)+ \\ \frac{{\mathsf{\Delta }}^2}{24}{\left.\frac{{\partial }^2\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{k}}\right)}{\partial {\mathsf{\phi }}_{\mathrm{m}}\partial {\mathsf{\phi }}_{\mathrm{n}}}\right|}_{\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{k}}\right)=\left(\overline{{\mathsf{\phi }}_1},\overline{{\mathsf{\phi }}_2},\dots ,\overline{{\mathsf{\phi }}_{\mathrm{k}}}\right)}\overline{\frac{\partial {\mathsf{\phi }}_{\mathrm{m}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\frac{\partial {\mathsf{\phi }}_{\mathrm{n}}}{\partial {\mathrm{x}}_{\mathrm{i}}}}+\mathrm{O}({\mathsf{\Delta }}^4) \end{gathered}$$

(3)

The unclosed scalar gradient term within, treated as the pseudo-scalar-dissipation form[2,17,18,24], is approximated by an algebraic expression approach [29]:

$$\overline{\frac{\partial {\mathsf{\phi }}_{\mathrm{m}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\frac{\partial {\mathsf{\phi }}_{\mathrm{n}}}{\partial {\mathrm{x}}_{\mathrm{i}}}}=\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{m}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{n}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}+{\mathsf{\chi }}_{\mathrm{sgs}}\text{≈}{(1+\mathrm{C}}_{\mathrm{sgs}})\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{m}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{n}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}$$

(4)

For non-premixed combustion, ${\mathrm{C}}_{\mathrm{sgs}}$ is adopted as 0.1 following [30,31]. Otherwise, ${\mathrm{C}}_{\mathrm{sgs}}$ can be provided by transport equation models at the cost of high computational demands, which remains to be explored in future.

Finally, the model explicitly turns out to be:

$$\begin{gathered} \overline{\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{k}}\right)}=\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left(\overline{{\mathsf{\phi }}_1},\overline{{\mathsf{\phi }}_2},\dots ,\overline{{\mathsf{\phi }}_{\mathrm{k}}}\right)+ \\ \frac{{\mathsf{\Delta }}^2}{24}{\left.\frac{{\partial }^2\dot{{\mathsf{\omega }}_{\mathsf{\alpha }}}\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{n}}\right)}{\partial {\mathsf{\phi }}_{\mathrm{m}}\partial {\mathsf{\phi }}_{\mathrm{n}}}\right|}_{\left({\mathsf{\phi }}_1,{\mathsf{\phi }}_2,\dots ,{\mathsf{\phi }}_{\mathrm{k}}\right)=\left(\overline{{\mathsf{\phi }}_1},\overline{{\mathsf{\phi }}_2},\dots ,\overline{{\mathsf{\phi }}_{\mathrm{k}}}\right)}(1+{\mathrm{C}}_{\mathrm{sgs}})\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{m}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\frac{\partial \overline{{\mathsf{\phi }}_{\mathrm{n}}}}{\partial {\mathrm{x}}_{\mathrm{i}}}+\mathrm{O}\left({\mathsf{\Delta }}^4\right) \end{gathered}$$

(5)

The final closure in expression (5) takes no presumption with respect to the combustion regime and, a-priori, it is applicable for premixed, non-premixed, and partially premixed combustion. Besides, accuracy is explicitly determined by the series truncation order. Finally, the approach preserves the reasonable restricting properties, where the SGS contribution decreases with the square of the filter width, see Relation (5), approaching DNS as  reaches Kolmogorov scales. For more details, one can refer to the study [24].

2.2. Numerical Implementation and Error Analysis

The series model is implemented on OpenFOAM platform [32]. In practice, the approach is incorporated in the reactingFoam solver through the manipulation on the chemical source in scalar balance equations. Besides, the library combustionModel is customised to import the sub-grid part of the series model into the source term interface. In order to utilise the information of Arrhenius reaction rate, the library chemicalModel is called. Then, the series model is coupled through the scalar transport equation with LES equations in the formation of a new reactingFoam-based solver. The PIMPLE algorithm is employed to deal with iterative procedures for coupling equations of momentum and mass conservation.

Despite the series model is derived mathematically, computational errors inevitably will occur, like any other sub-grid combustion models. The sources of error can potentially stem from:

First, high-order terms. The accuracy of the model is mainly determined by the order of Taylor series. In application, terms of the fourth and higher order are neglected, as they are not discovered to act a major role [23]. Nevertheless, under some extreme conditions, for example in modelling deflagration and detonation transition, these higher order contributions can significantly affect the flow field and lead to inaccuracy of predicting the sub-grid influence.

Second, the chemical mechanism. The series model requires the 2nd-derivative of the chemical source, which may become noisy under certain conditions. The choice of a well-established chemical scheme is crucial.

Thirdly, the scalar gradient term. The approximation in Equation (4) could be inaccurate in poorly resolved flames, where scalar gradients are large in the periphery of the reaction zones. This can be improved by employing a better refined/adaptive mesh.

3. Results and discussion

3.1. Non-premixed combustion simulation

The case under consideration is Sandia Flame D, an ideal experimental benchmark for non-premixed combustion, broadly utilised for model verification in past studies due to the large set of experimental data available. Experimentally, this jet flame has been studied on a piloted burner at Sydney University [33] by Barlow [34] who measured the temperature and species field, and Schneider [27] who provided LDV velocity details. The burner has an inner nozzle with a diameter D = 7.2 mm, where a mixture of 25% methane and 75% air by volume is injected at different bulk velocities 49.5m/s, resulting in Re = 22400. The pilot nozzle expands to a diameter of 18.2mm, coaxial to the main jet. From it exits a lean (Φ = 0.77, T = 1880K) mixture of C₂H₂, H₂, air, CO₂, and N₂. Surrounding the pilot nozzle air streams at a velocity of 0.9m/s. The burner exit is positioned approximately 15 cm above the exit of the vertical wind tunnel [34]

The details of modelling parameters are exhibited in Table 1, along with comparisons with previous studies. According to it, the radial domain widely chosen ranged from 8 to 40D in the radial direction, while the axial part from 35 to 150D. In our work, the domain extends 70D in axial directions. Previous studies [35,36] suggest that clapping the grids to 50D does not affect the results, and in the radial direction it has an increase from 15 to 30D to capture the downstream flame/temperature expansion. The mesh is aligned with polar coordinates, but in the centre, a square section (or o-grid) is applied to avoid very fine meshes in the centreline. Two resolutions are employed under the Pope criterion [37]: a fine grid of 71 points in the tangential direction, 48 points in the azimuthal direction and 210 points in the axial direction, with an o-grid area of 12×12; a coarse one, which is approximately downscaled with a factor of 1.5 in each direction (52, 36, 139), with an o-grid zone of 9×9. Both grids stretch along the axial and tangential directions to handle the inlet variance and the strong gradients in shear layers. The sketch of simulation domain and grid allocation is depicted in Figure 1.

Snapshots of the simulated instantaneous fully-developed fields are shown in Figure 2. From the temperature and CO₂ distribution, the predicted flame structure displays the conventional characteristics of non-premixed combustion [2] as expected: a preliminary area close to the inlet nozzle where flames are thin and a subsequent zone beside which hot products fill the majority of the downstream realm.

Statistical collection is performed over 16 burner flow-through times based on the jet bulk velocity. In addition, 10 flow-through times are initially computed to guarantee the flows fully developed. To compare, published Flame D predictions by [43] in Table 1 are imported in the following result plots. The main reason is that it adopts the same 4-step chemistry mechanism and uses the sophisticated Eulerian stochastic field method for combustion modelling. This will allow us to evaluate the performance of the series method against a well-established and reliable model avoiding complexities associated with the chemical mechanism. The results (fine-grid set) presented are carried with a similar resolution to the fine mesh employed in this study.

Figure 3 displays the mean and root-mean-square (RMS) trends of mixture fraction, temperature, axial velocity and species along the centreline, and Figure 4 show the radial distributions. The mixture fraction is defined following Bilger [34]. Along the centreline, the series model shows good prediction with mixture fraction experimental data in the fine grid, with slight over-prediction from z/D=10 to 40, illustrating that the flame shape is replicated reasonably. For the coarse resolution, the result shows a reasonable trend, but diverges off the fine grid prediction from z/D = 15. This is largely due to the grid insufficiency of the coarse mesh in the downstream part. However, the situation is opposite in the reference case [43], where over-prediction is observed in the immediate vicinity of the inlet. The authors attributed it for the limited grid resolution of the mixing layers evolving between the jet and pilot close to the nozzle. Switching to the radial distributions in Figure 4, the mixture fraction profiles are well reproduced by the series model, although small over-predictions of the peak values are found at z/D =30 and 45 in the coarse grid simulation. It is consistent with the centreline observations. The results from the stochastic fields method [37] are comparable in terms of accuracy with the series results although some under-prediction in the inner (fuel-rich side) regions at z/D =7.5 and 15 is noticed, which is not observed in the series model.

The axial profile of mean and RMS temperature is also displayed in Figure 3. Along the centreline, the mean temperature, the flame location accordingly, is well reproduced for the series model, and the fine grid prediction behaves better than the coarse one as expected. The temperature RMS keeps the same level as the experimental data, although small over-prediction appears in the upstream locations before z/D =25. This evidences that the temperature fluctuation is sufficiently resolved by the series model with a simplified chemical mechanism. Besides, the RMS trend duplicates the observations in the experiments that the minimum arises around the location of the maximum mean profile. For the radial profiles in Figure 4, the statistical distributions agree well with the experiments, albeit with slight over-predictions of the mean values at z/D =15 and 30. The over-predictions are also spotted in the reference case [43] and the discrepancies may originate from the thermal radiation that is not considered in both LES simulations. Besides, it should be noticed that the series model performs relatively better than the PDF method [43] at predicting the location of the mean peak at the near-nozzle positions z/D=7.5 and 15. This is related to the under-prediction of mixture fraction profiles in the stochastic field simulation[43]. The centreline distribution of velocities is displayed in Figure 3 as well. Clearly, those quantities are decently reproduced by the series model, in spite of some acceleration after 15D in the coarse grid simulation. For the stochastic field model [43], it predicts a relatively better trend.

Figure 3 also exhibits the axial profiles of reactive species. Observing the series model results, the reactants CH₄ and O₂ tend to consume faster than the experimental data, consistent with the slightly over-predicted temperature trend. The predictions of CO₂ and H₂O agree well with the measurements, though coarse grid results are less accurate. Nevertheless, the prediction of CO and H₂ is significantly over-predicted. The discrepancies as well as the accelerated CH₄ decay rate are like those observed in the simulation [43] and it can be attributed to the limitations of the simplified mechanism [42]. The findings[43,59] suggest that the C1 scheme in the mechanism[42] gives rise to an over-prediction of reaction rates on the fuel-rich side of the non-premixed flames. It is evidenced in the radial distribution of CH₄ in Figure 4. Besides, the reduced reaction mechanism suffers from the shortcomings of predicting intermediates like H₂ and CO as it is susceptible to diffusive transport [42]. This accounts for the inferior distributions of the same intermediate species in Figure 3. In spite, good radial agreement is achieved for both the series model and stochastic fields method [43], regarding the mean and RMS of CO₂ and H₂O. The peak locations near the nozzle exit are slightly under-predicted in the reference case [43], owing to the same reason of the temperature distribution.

Overall, both the mean and RMS behaviours of the non-premixed flame are well captured by the series model. Improvements are obtained with grid refinement, especially in terms of reproducing the temperature, mixture fraction and major species distributions. In contrast to the stochastic field approach using the same mechanism [43], the series model demonstrates comparatively good prediction in general, and even better at some occasions: like correctly capturing the radial distributions of mixture fractions and the peak locations of mean temperature and chemical species. This gives the series model a advantage to some extent since its formulation is simpler, the cost cheaper, and can be incorporated directly to the species transported equation.

3.2. Premixed combustion simulation

In this section, the methane/air turbulent Bunsen flame F3 experimented by Chen [26] is chosen. The configuration is a typical turbulent premixed flame with a wide range of available velocity, temperature and species experimental data provided in [25], and has been extensively applied for model validation in the combustion community. The burner has an inner nozzle with a diameter H = 12 mm, from which a stoichiometric mixture of methane and air is injected at different bulk velocities 30 m/s, leading to Re = 23,000. Surrounding the main fuel injection, a laminar stream (the burning products of stoichiometric methane and air) is piloted via a perforated plate for stabilising the turbulent flame [60], and the outer diameter is 68mm. Outward, the fresh air is entrained as a form of a low-velocity co-flow. In the diagram of the premixed regimes, the F3 flame lies in the thin reaction zone, more precisely near the flamelet regime.

Simulation details are presented in Table 2, along with the comparison with previous research. The simulated domain extends 30H downstream the nozzle, and 12H in the radial direction so as to capture the downstream flame/temperature expansion. As Table 2 summarises, the chosen domain is sufficient for the flame propagation. The grid used is aligned with polar coordinates, but in the centre, a square section (o-grid) is applied to avoid very fine meshes in the centreline as in Sandia Flame D. Two resolutions are employed following the Pope criterion [37]: a fine grid of 69 points in the tangential direction, 48 points in the azimuthal direction and 200 nodes in the axial direction; and a coarse one, which is approximately downscaled with a factor of 1.5 in each direction (49, 36, 134). Grids are stretched along the axial and tangential directions to seize the strong gradients stemming close to the nozzle and shear layers. The sketch of simulation domain and grid allocation is detailed in Figure 5.

On the inlet boundary, velocity is prescribed with the detailed flow field experimental data [26] above the nozzle exit. The inlet turbulence intensities are imposed using the same approach as the Sandia D flame. For the main jet species, flat profiles are specified using the measurements[26]. The composition of the pilot stream is following [60]. The pilot temperature is presumed as 1785K [60,73,74], corresponding to 20% heat loss. Free slip is used for the lateral boundary, while a non-reflecting outflow condition is used on the outflow plane. The numerical schemes, the time set-up and the chemistry mechanism are the same as in Sandia Flame D simulations.

Statistical collection is performed over 15 burner flow-through times based on the jet bulk velocity. Before it, 8 flow-through times are computed to guarantee the flows fully developed. As a comparison, a set of F3 flame predictions published by [22] in Table 2 are selected, since it adopts a similar mathematical approach based on Taylor series expansion on the scalar field to determine the F factor in the ATF model. The augmented reduced version of GRI 3.0 schemes is introduced to denote chemistry. Besides, the mesh resolution is higher than the fine grid in this work.

The radial profiles of mean temperature and reactive species at different streamwise locations are first depicted in Figure 6. In overall, the mean temperature is over-produced, especially near the nozzle. However, different researchers [60,64,66,67,69,75] reported the similar discrepancies related to experimental measurements (largely extending 10–25 % experimental uncertainties) using either reduced or detailed mechanisms. This discrepancy is largely due to the uncertainties of temperature fields in the pilot stream, where a large and unidentified portion (not reported in experiments [26]) of heat is lost to the burner exteriors [64,67]. Past studies in Table 2 indicate that the loss proportion varies between 0[64,65] and 34%[60,70]. On the other hand, the simulation domain starts just downstream the nozzle exit for computational cost reduction, and wall heat transfer is not considered.

The radial trends of species mass fraction are displayed in Figure 6. The profiles of CH4 and O2 are reasonably replicated; In spite of some minor under-predictions on the fuel-lean sides (r/H>0.5), the calculated reactants meet the measurements well at the fuel-rich side and keep the descending trend from upstream to downstream. This indicates that the consummation rates of CH₄ and O₂ are well reproduced along the centreline. On contrast, the reference results [22] are under-predicted on both the fuel-lean and fuel-rich sites.

The radial distributions of CO₂ and H₂O by the series model are in a reasonably decent consistence with experimental observations, despite of some under-predictions at z/H = 6.5. Meanwhile, improvements are obvious with the increase of grid resolutions. Compared with the reference case [22], the series model performs better on the fuel-rich sides. It is related with the situation that the profiles of the reactants CH₄ and O₂ are not well resolved in these areas [22]. Note that in ATF models, although the thickened flames are solved, the species transport equation is modified, the interaction between combustion and turbulence is transformed from a transport-dominant combustion regime to a chemistry-dominant one, and the impact of the heat release upon the flow field is not represented sufficiently. In contrast, the series model operates on the reaction rate term directly, without altering the formation of species balance equations. In terms of CO, the series model also obtains a good prediction, while the profiles are slightly under-predicted in the reference case[22]. On the other hand, the computed CO products are less sensitive to grid resolutions than major species, as less difference is found with grid refinement. The measurement error of major species varies from 8% to 15% and that of intermediate species reaches 25% [26]. The major and intermediate species predicted by the series model principally cater to the accuracy in both grid resolutions.

To further inspect the impact of the temperature discrepancies, a new progress variable is established by the definition of a linear combination of species concentration [76]:

$$c=\frac{Y_{{\mathit{CO}}_2}+Y_{\mathit{CO}}}{({Y_{{\mathit{CO}}_2}+Y_{\mathit{CO}})}_{\mathit{pilot}}}$$

(6)

The radial distribution is depicted in Figure 7 and compared with experimental measurements. In this regard, the progress variable based on species demonstrates a reasonable and much better agreement than the temperature profile. The trend is refined even in the vicinity of the nozzle exit, where substantial discrepancies in the mean temperature are detected.

Usually in non-premixed combustion where the flame is thick, one would eventually start resolving it as the mesh is refined. However, in premixed combustion, the reaction zones are very thin and they would always remain under the grid size in LES. In general, with the refinement of grid resolution, the series model improves at predicting the premixed Bunsen flame. Comparing the series model with the simulations of [22], which adopted a similar series approach to determine the F factor of the artificially thickened flame model, the predictions are as good, if not better. Considering that the reference simulation[22] employs a finer mesh and also more detailed chemistry, the series model shows the potential improvement if finer meshes are used or/and more detailed chemistry.

4. Conclusions

The generality and prediction accuracy of the generalised series model for LES of premixed and non-premixed turbulent combustion is presented. It is based on Taylor series expansion of the chemical source term in scalar space around the filtered value and implemented into OpenFOAM. The mathematical model does not depend on the combustion regimes, has the correct limiting behaviour, and determines the order of accuracy explicitly.

The model is first applied to a piloted methane/air non-premixed jet flame (Sandia Flame D) with two mesh resolutions. Overall, the statistical (time-averaged and RMS) results agree well with experimental measurements, particularly with regard to the mixture fraction, velocity, temperature, and the concentrations of major reactive species CH₄, CO₂, H₂O, and O₂. The concentrations of intermediates CO and H₂ are over-predicted, due to the limitations of the reduced reaction mechanism employed. Then, a Bunsen piloted premixed flame is simulated. Most of the statistical properties of both reactive species and progress variable are also well reproduced by the series model. The only major discrepancy evident is in the temperature and it is attributed to the experimental uncertainties of temperature fields in the pilot stream.

Overall, the results have demonstrated the model’s generality for both premixed and non-premixed combustion simulation, and its accuracy of prediction of reactive species distribution. Future research will extend the present work to partially premixed flames with detailed chemistry.