University of Birmingham Combustion characteristics of a non-premixed oxy- flame applying a hybrid filtered Eulerian stochastic field/flamelet progress variable approach

The oxidation of methane under oxy-fuel combustion conditions with carbon capture is attractive and deserves huge interest towards reducing CO2 and NOx emissions. The current paper reports on the predictions and analysis of combustion characteristics of a turbulent oxy-methane non-premixed flame operating under highly diluted conditions of CO2 and H2 in oxidizer and fuel streams, respectively. These are achieved by applying a novel, well-designed numerical combustion model. The latter consists of a large eddy simulation (LES) extension of a recently suggested hybrid model in Reynolds averaging-based numerical simulation (RANS) context by the authors. It combines a transported joint scalar probability density function (T-PDF) following the Eulerian Stochastic Field methodology (ESF) on the one hand, and a flamelet progress variable (FPV) turbulent combustion model under consideration of detailed chemical reaction mechanism on the other hand. This novel hybrid ESF/FPV approach removes the weaknesses of the presumed-probability density function (P-PDF)-based FPV modeling, along with the solving of associated additional modeled transport equations while rendering the T-PDF computationally less affordable. First, the prediction capability of the LES hybrid ESF/FPV was appraised on the well-known air-piloted methane jet flame (Sandia Flame D). Then, it was assessed in analyzing the combustion properties of a non-premixed oxy-flame and in capturing the CO2 dilution effect on the oxy-fuel flame behavior. To this end, the so-called oxy-flame B3, already numerically investigated in a RANS context, was analyzed. Comparisons with experimental data in terms of temperature, scalar distributions, and scatter plots agree satisfactorily. Finally, the impact of generating the FPV chemistry table under condition of unity Lewis number, even with CO2 dilution, was investigated on the general prediction of the oxy-fuel flame structure, stability and emissions. In particular, it turns out that 68% molar percentage of CO2 leads to 0.39% of CO formation near the burner fuel nozzle and 0.62% at 10 dfuel above the nozzle.


Introduction
Oxy-fuel combustion is currently known to have significant advantages over conventional combustion concepts where fossil fuel is oxidized with air, affecting the environment in terms of emission increase of carbon dioxide (CO 2 ), NO x formation through the oxidation at high temperature of nitrogen (N 2 ) in air. According to the Emission Database for Global Atmospheric Research [1], global emission of CO 2 increased to 35.8 billion tones in 2016, compared to 22.06 billion tones registered two decades ago. The concept of carbon capture and storage (CCS) is actually one accepted strategy towards cleaner combustion of fossil fuels. Three main approaches are currently emerging, namely, pre-combustion capture, post-combustion capture, and during-combustion capture or Oxy-fuel combustion (see in [2]). Thereby, the concept of oxy-fuel combustion has been especially promoted. It is mainly based on replacing the air containing nitrogen with enriched air or pure oxygen as oxidizer, where the resulting flue gas is composed primarily of CO 2 and H 2 O, from which CO 2 can be easily separated, enabling its capture, storage or recycling [2][3][4]. The oxy-combustion process is classified as a high thermochemical process with faster chemical reaction and burning velocity when compared to air combustion. Thus, in order to control the flame's temperature, structure and heat transfer, recycled CO 2 is often added to the oxidizer stream as a diluent [3].
With all these challenging changes with respect to the traditional combustion environment, it is essential to better understand the impact of the oxy-fuel conditions on the flow properties and flame characteristics. Relying on numerical investigations, appropriate numerical simulation tools are highly needed. Focusing on the description of complex turbulent flow, the large eddy simulation (LES) technique has proven to be able of to capture very complex turbulent reacting flows better than Reynolds averaging-based numerical simulation (RANS) [5]. Using a LES filter, large turbulent flow structures are classically separated from the unresolved scales. The former are explicitly computed while the unresolved structures need to be modeled by means of sub-grid scale models. The classical LES approach [5][6][7] is applied in this paper in conjunction with a tabulated chemistry-based combustion modeling.
As the chemistry has to be considered within the LES framework, the integration of the required chemistry equations needs special care, in order to reduce the overall computational cost related with detailed chemistry computations essential for fully describing the chemistry process evolving. That is why various chemical mechanism reductions and chemistry tabulation or storage/retrieval approaches are now available in the literature and widely utilized, as pointed out in [5,8,9]. The mathematical procedures do not exclude species. The species concentrations are rather calculated or by tabulating them as functions of a few preselected controlling progress variables. Dealing with tabulated chemistry techniques, the most common methods used in different studies include the intrinsic low dimensional manifold (ILDM) in [10], the reaction-diffusion manifold (REDIM) [11] flamelet/progress variable approach (FPV) in [4,8], flame prolongation of ILDM (FPI) or flamelet generated manifold (FGM) in [12], and filtered tabulated chemistry for LES (F-TACLES) [9,13]. Thereby, each chemical look-up table corresponds to one set of boundary conditions in terms of fuel and oxidizer compositions and is mapped by means of a small number of reaction progress variables (RPV). When dilution by burnt gases and heat losses or multiphase phase process effects are to be considered, additional RPV may dramatically increase the size of the tables to be stored, resulting in huge look-up tables. To overcome this disagreeability, novel techniques were recently proposed, like the so-called hybrid transported-tabulated chemistry (HTTC) approach, combining the transport of the main species in the flow with the tabulation of the radical species [14]. The tables resulting from these techniques are coupled either to RANS-or to LES-based Computational Fluid Dynamics (CFD). Focusing especially on the flamelet/progress variable (FPV) method as suggested by Pierce in [8], the mixture fraction and the progress variable, tracing the local reaction progress, are the controlling variables for the multidimensional chemistry. Within the FPV-based combustion model, they are solved by their own transport equations in addition to the continuity, momentum, energy/enthalpy equations in single phase reacting flow cases (e.g., [4,8,[15][16][17]). This method has been widely used, especially for CH4/Air Appl. Sci. 2019, 9,1320 3 of 23 jet flame-D by various authors (e.g., [4,15,18]) and for predicting extinction in Sandia piloted jet flames (E and F) (e.g., [16]).
Dealing with oxy-combustion modelling, several advanced approaches have been evaluated to study the turbulence-flame interaction phenomena. In the numerical work of Mayer et al. in [19], both the laminar steady flamelet model (SFM) and the eddy dissipation concept model (EDC) were used to study the influence of additional oxy-fuel jet burners in gas phase on the productivity of the investigated furnace. Also, Prieler et al. in [20], in their oxy-combustion study, compared the same methods in a RANS context with the eddy dissipation model (EDM) for various oxygen enrichments. Both studies employed different chemical mechanisms. However, Mayer et al. in [21] discussed the usability and limits of the SFM approach in oxy-fuel combustions. Hidouri et al. in [22] applied in the LES framework an alternative approach based on the joint PDF method coupled to FGM technique based on the mixture fraction and the reaction progress variable to investigate the behavior of two reacting separated oxy-fuel jets. Thereby the sub-grid PDF shape is described by the presumed beta-PDF assumption. However, according to [23], not all transported scalar properties could be appropriately captured using presumed-PDF methods. Nevertheless, the conditional moment closure (CMC) combustion model has been used for different turbulent studies, including the characteristics of turbulent combustion of natural gas flame in oxy-fuel combustion environments by Kim et al. in [24] and the oxy-fuel coaxial jets by Garmory et al. in [25]. However, very few contributions dealt with oxy-combustion using the FPV approach. One of the early turbulent simulations incorporating the FPV approach in the LES framework has been reported in [17], where it was shown that the LES-FPV approach possesses in turbulent regimes potential promise in addressing the issue of differential diffusion in such a flame. Similar to FPV, the FGM method has been used in [22] with Le = 1 to tackle oxy-fuel combustion process in twin-jets.
It is worth noting that both FPV and FGM parameterize all species and temperature by a mixture fraction and a reaction progress variable or reaction progress parameter. However, the two models treat the flamelet manifold between equilibrium and flame extinction in a different manner (see in [26]). As pointed out in [26], the FPV model solves flamelet equations in the (unstable) middle and lower branches of the flamelet S-curve, while the FGM model solves an unsteady extinguishing flamelet from the last stable burned steady solution before extinction. Nevertheless, the generated tables show similar behavior in the mixture fraction space while exhibiting different behaviors on the progress variable space. In the current study, the FPV approach is especially adopted due to its better features in both lean and rich flame conditions. It is then used at the first time in coupling with LES to investigate non-premixed oxy-flames.
It is worth mentioning that using such a tabulated chemistry technique calls for assumptions regarding the flame structure along with the statistical turbulence-chemistry interaction (TCI) phenomena. For that purpose, a probability density function (PDF) of the controlling parameters of the table is usually introduced. To determine this PDF, two main methodologies are suggested in the literature [4,5,7,[27][28][29][30][31][32][33][34][35][36][37][38][39][40]. The first method consists in assuming a presumed shape for the PDF as reported in the most contributions which apply the FPV approach [8,10,[16][17][18]. Such a shape assumption for the PDF for all scalar distributions calls usually for considering statistical independence between single PDF, which is not the case in reality. Very few considerations relax this assumption, as in [41]. It is noteworthy that the chemical source term within this first method is expressed in an unclosed form. It is usually modeled from the FPV table as function of the reaction progress variable along with the presumed PDF. According to previous LES studies in oxy-flames [17,22], the application of a presumed shape of PDF with Beta-function experienced shortcomings regarding the general prediction capability of the temperature distribution and minor species evolution. As pointed out in [17], the results obtained within the LES-FPV framework are affected by the pre-calculated flame structures through the tabulated transport properties and source terms from the manifolds, and all data are additionally PDF-integrated. Similar limitations were reported within the RANS context on oxy-fuel configurations in [4]. The second method to determine the PDF is therefore adopted in this study. It is a more efficient approach that describes the PDF by solving a transported probability density function (T-PDF) equation, as reported, for example, in [4,[27][28][29][30][31][32][33][34][38][39][40].
Two major formulations of the T-PDF method have been proposed in the literature, namely the Lagrangian Monte-Carlo approach as applied in [28,29] and the Eulerian based approach in [30][31][32][33][34][37][38][39][40]. Solving a filtered transported PDF equation means applying a modeled form of the equations of the PDF or filtered density function (FDF) in LES, where the filtered and the sub-grid scalar field can be formulated. One attractive approach to model the FDF equations consists of using the Eulerian stochastic field (ESF) approach [31][32][33] and [34].
Following Valińo et al. in [35] and Sabel'nikovet al. in [36], multiple stochastic fields for each scalar are utilized to describe the FDF undergoing diffusion, turbulent convection and chemical reaction. The resulting stochastic fields (SFi) are continuous and differentiable in space. Within this method, the chemical source term emerges in a closed form, allowing modeling of the combustion and the involved TCI in an accurate way. Many previous LES investigations have employed the ESF method coupled with different chemical mechanism approaches in both simple and complex scales to simulate various jet flames in an air environment. In particular, the ESF method was applied with reduced mechanisms in order to study premixed and partially-premixed flames in [31,32] and [37,42]. From these studies, it was concluded that 8SFi could well characterize the influence of the sub-grid fluctuations on the combustion properties. Jangi et al. in [35,38,39] employed the ESF method in both RANS and LES contexts in combination with detailed chemistry. Even though all these previous studies proved the efficiency of the ESF approach in a LES context, the choice of the chemical mechanism technique has an enormous effect in terms of computational requirements.
Only very few studies have focused on developing a LES framework for the ESF method in coupling with a tabulated chemistry. In their large eddy simulation (LES) study, Amer et al. [33] reported a coupling of the ESF approach with the FGM-based tabulated approach. To the knowledge of the authors, coupling the ESF method to the FPV tabulated chemistry strategy is new in the LES context. In an oxy-fuel environment, Rihab et al. [4] were the first to develop a hybrid ESF/FPV approach within the RANS framework and to validate it in three oxy-fuel jet flames (A1, A3 and B3) as experimentally investigated in [43]. They compared the results with those obtained by means of presumed PDF and with experimental data. In contrast, the oxy-fuel jet flames (A1 and A3) have been investigated in LES, including a Smagorinsky model and a conditional moment closure (CMC) combustion model by Garmory and Mastorakos [25]. Since a unity Lewis number diffusion model was used, the model suffers from real limitations in predicting the main flame characteristics and major reacting species as well as the temperature evolution at different positions.
Relying on the abovementioned limitations of the LES-FPV approach and on the attractive features of the hybrid ESF/FPV approach, it is promising to extend the hybrid ESF/FPV approach in which the TCI is well accounted for to LES. This task is attempted in the present paper, where the resulting LES hybrid ESF/FPV approach includes in a first step FPV table generated with Le = 1 in order to evaluate its effect on the prediction of the oxy-flame properties.
The current paper reports on the predictions and analysis of combustion characteristics of a turbulent oxy-methane non-premixed flame operating under highly diluted conditions of CO 2 and H 2 in oxidizer and fuel streams, respectively, by rather applying a novel, well-designed numerical combustion model. The latter consists of a LES extension of a recently suggested hybrid ESF/FPV approach in RANS context by the authors [4]. It combines a transported joint scalar probability density function (T-PDF) following the Eulerian stochastic field methodology (ESF) and a flamelet progress variable (FPV) turbulent combustion model which includes detailed chemical reaction mechanism. This hybrid ESF/FPV approach overcomes the limitations of the presumed-probability density function (P-PDF) by accurately characterizing the influence of the sub-grid fluctuations on the flame structure and on combustion properties, along with the TCI. The prediction capability of the LES hybrid ESF/FPV was first appraised on the well-known air-piloted methane jet flame (Sandia Flame D). Then, it was evaluated in analyzing the combustion properties of a non-premixed oxy-flame and in capturing the CO 2 dilution effect on the oxy-fuel flame behavior. To this end, the turbulent non-premixed methane oxy-flame, the so-called oxy-flame B3, as experimentally studied in [43] and already numerically investigated in RANS context in [4], was analyzed. Finally, the impact of generating the FPV chemistry table under condition of unity Lewis number even with CO 2 dilution was examined on the general prediction of the oxy-fuel flame structure, stability and emissions. In this configuration, a significant CO 2 dilution level was considered in the oxidizer and a high percentage of H 2 enrichment in the fuel stream. As pointed out in [44], the CO 2 diluted oxy-fuel flame evolves differently compared to a classical air-fuel flame and more stabilization mechanisms are required. These changes in the combustion environment of the oxy-fuel case make it more attractive and challenging for sub grid scale and combustion modeling.
The remainder of the paper is structured as follows. Section 2 outlines the numerical methodology. Section 3 presents the experimental configurations with corresponding numerical set-up details. Section 4 outlines and discusses the obtained results. Concluding remarks are summarized in the last section.

Large Eddy Simulation Transport Equations
In this section, the LES hybrid ESF/FPV approach is presented. Actually, it combines in the LES framework a transported joint scalar probability density function (T-PDF) following the Eulerian stochastic field methodology (ESF) and the flamelet progress variable (FPV) turbulent combustion model under consideration of detailed chemical reaction mechanism.
The filtered governing equations for describing a turbulent reacting Newtonian fluid flow are obtained after a filtering procedure in space from which the large-scale energetic motions are resolved and the unresolved sub-grid scale structures are accounted for by so-called sub-grid scale (SGS) models (e.g., [5][6][7]). The filtered transport equations for LES emerge then in the condensed form as where the quantity ϕ is defined as ϕ = ρϕ/ρ [5,45]. The symbols (.) and (.) express spatially filtered and density-weighted (Favre-filtered) quantities, respectively. In Equation (1), u i is the velocity component (with i = {1,2,3}) and x i the position vector in i-direction. The quantity ϕ may represent unity, the velocity component, thescalar mass fraction, Y, for the equation of mass, momentum, and progress variable transport, respectively. Q i represents the molecular flux vector quantity and Q i sgs a SGS flux counterpart, while π stands for a supply term (viscous/pressure term, possible radiation) and S C for the chemical source term. To close Equation (1), models for the unclosed SGS terms Q i sgs and S C are required. In the momentum and scalar transport equations this term is the SGS stress tensor and scalar flux vector, respectively, expressed as: Various models have been suggested in the literature (see [5,7,46]) to close these terms. In the present paper, the SGS stress tensor is postulated following the eddy viscosity approach as: where µ sgs = ρ c m ∆ grid 2 S ij S ij 1/2 and S ij = 1 2 ∂ u j ∂x i , according to Smagorinsky [6]. Thereby, C m is a model coefficient and ∆ grid = ∆ x ∆ y ∆ z 1 3 the grid filter width. To account for the stream jet's wall surface, the van Driest wall damping function is introduced [45], so that the grid filter width ∆ grid is replaced by the expression: where k r is the Karman constant and r + is the dimensionless wall distance. The SGS scalar flux Q sgs iϕ α is modelled as: wre α refers to species and σ sgs expresses the SGS Schmidt number. According to [47,48] the values of the model coefficients applied are set as: Restricted ourselves to the FPV approach in which the remaining quantity S C in Equation (1) for the mass density, momentum and species mass fraction, Y k , respectively. Thereby, P stands for the mean pressure, µ for the dynamic molecular viscosity, and . ω Y k = S C for the mean chemical source term of the species mass fraction.

Chemistry Reduction with Flamelet/Progress Variable Approach (FPV)
In the classical FPV approach, the problem of turbulent combustion model and simulation consists of determining the mean chemical source term, . ω Y k , using the appropriate chemical reaction model or mechanism to solve Equation (8). In the present paper, a tabulated chemistry method following the flamelet/progress variable (FPV) approach [7] was employed. Accordingly, a steady-state flamelet library was generated based on two controlling variables, namely the mixture fraction and a reaction progress variable. All relevant information (e.g., the chemical rates of formation/destruction, chemical composition, temperature, density, viscosity, etc.) was therefore stored in the tables as a function of just the two variables. In this way, other flame thermochemical states (i.e., fully burning state, transient solution between burned and unburned states and the unburned state) can be parameterized with the advantage that the reaction progress variable allows description of local extinction and reignition phenomena [15,16,18,49]. In order to build such FPV tables, a laminar counter-flow flame was first solved using the flamelet generator FlameMaster from [50], where the chemical mechanism in use consists of 325 reactions and 53 species available in GRI-MECH 3.0 [51]. The results obtained from 1D flamelets were then mapped into the mixture fraction and progress variable space. Finally, each instantaneous thermochemical state was retrieved by means of the two controlling variables. Given that different definitions of the reaction progress variable (PV) could be used [52], the following progress variables were employed in the current study for the Sandia flame D and the oxy-flame B3, respectively, in accordance with [4,25].
where [-] denotes the specie mass fraction divided by the specie molar mass. Figure 1 presents some features from the FPV tables created according to the progress variable for both Flame D and Flame B3 (see Equations (9) and (10)). Thereby the distribution of the heat release production rate is shown colored by the reaction rate of the progress variable in the mixture fraction space. Considering the high amount of CO 2 diluted in the oxidizer in Oxyflame B 3 instead of N 2 as it is in Flame D and since the heat transfer properties of CO 2 are different from those of N 2 , one can expect a significant variation of the production rate of the heat release. In Figure 1a referring to Flame D with 79% of N 2 in the oxidizer, the maximum value of heat release locally reaches 4.10 9 J/m 3 s. It is in the oxy-flame B 3 radically increased to 7.10 9 J/m 3 s where an amount of 68% of CO 2 is injected as diluent instead of N 2 (Figure 1b).

Turbulence Chemistry Interaction, Transported Filtered Joint Probability Density Function and Eulerian Stochastic Field Approach
To achieve complete statistical modelling of turbulent combustion it is essential to consider both sub-grid scale turbulent motions and chemical reactions along with their interaction in a reliable way. Approaches based on solving the transported probability density function offer compelling advantages in providing an effective resolution to the closure problems that arise from averaging or filtering of the highly nonlinear chemical source terms, from modeling the effects of convection and body forces. In [4] the predictive capability of the hybrid ESF/FPV method in RANS was very well established. This is due to the inherent property of the transported PDF method to include the complete statistical information about its variables. Instead of solving Equation (8) as in the classical FPV approach, for which a presumed PDF is usually adopted, the present paper adapts the RANSbased hybrid ESF/FPV approach to LES. This allows accurate description of the turbulent flow field and to reliably capture the turbulence chemistry interaction, while providing both the chemical source term without model and the SGS PDF shape without assumption.

Transported Filtered Joint Probability Density Function
According to Pope in [53] and Jones in [31,32,54], a transported equation of a density weighted sub-grid probability density function (PDF) for different scalar quantities required to describe the reaction (FPV controlling variables) and for given instance t and point , can be derived from Equation (1) as follows: where the assumption of equal diffusivity has been applied [55]. In Equation (11),

Turbulence Chemistry Interaction, Transported Filtered Joint Probability Density Function and Eulerian Stochastic Field Approach
To achieve complete statistical modelling of turbulent combustion it is essential to consider both sub-grid scale turbulent motions and chemical reactions along with their interaction in a reliable way. Approaches based on solving the transported probability density function offer compelling advantages in providing an effective resolution to the closure problems that arise from averaging or filtering of the highly nonlinear chemical source terms, from modeling the effects of convection and body forces. In [4] the predictive capability of the hybrid ESF/FPV method in RANS was very well established. This is due to the inherent property of the transported PDF method to include the complete statistical information about its variables. Instead of solving Equation (8) as in the classical FPV approach, for which a presumed PDF is usually adopted, the present paper adapts the RANS-based hybrid ESF/FPV approach to LES. This allows accurate description of the turbulent flow field and to reliably capture the turbulence chemistry interaction, while providing both the chemical source term without model and the SGS PDF shape without assumption.

Transported Filtered Joint Probability Density Function
According to Pope in [53] and Jones in [31,32,54], a transported equation of a density weighted sub-grid probability density function (PDF) for different scalar quantities required to describe the reaction (FPV controlling variables) and for given instance t and point x i , can be derived from Equation (1) as follows:

Transported Filtered Joint Probability Density Function
According to Pope in [53] and Jones in [31,32,54], a transported equation of a density weighted sub-grid probability density function (PDF) for different scalar quantities required to describe the reaction (FPV controlling variables) and for given instance t and point , can be derived from Equation (1) as follows: where the assumption of equal diffusivity has been applied [55]. In Equation (11), with being the joint PDF. The quantity stands for sample space for the variable ɸ and represents the Dirac delta function. The density weighted sub-grid PDF, P , then represents the probability to capture the variable ɸ between two succeeding phase spaces and + , for different realizations of the controlling variable within the filter volume described in Equation (1). On the left hand side of Equation (11), the first term ❶ represents the where the assumption of equal diffusivity has been applied [55]. In Equation (11), being the joint PDF. The quantity ψ stands for sample space for the variable F α and δ represents the Dirac delta function. The density weighted sub-grid PDF, P sgs , then represents the probability to capture the variable F a between two succeeding phase spaces ψ and ψ + dψ, for different realizations of the controlling variable within the filter volume described in Equation (1). On the left hand side of Equation (11), the first term represents the accumulation part describing the rate of change in the physical space, the second term is the convection term due to the filtered velocity, and the last term stands for the closed chemical production source term in the phase space. Unlike the three terms defined already which are in closed form, the terms in the right hand side of Equation (11) require modelling as they are unclosed. The fourth term describes the transport at the sub-grid level of the FDF. The last term represents the influence of the molecular mixing at scales smaller than the size of the filter width. The transport of the SGS FDF is postulated by a simple gradient ansatz in accordance to the Smagorinsky model used above. To close the micro-mixing term, the linear mean-square estimation model approach (known as interaction by exchange of mean, IEM) [56] is applied in this work for its simplicity. Correspondingly, the modeled expression for the filtered joint probability density function P sgs can be written as In Equation (13) the sub-grid mixing time scale τ sgs is modeled by using a mixing time model referring to Jones et al. [54] as: It is worth noting that various approaches for the micro-mixing modeling exist. Recent reviews can be provided in [57,58]. It turns out that the coefficient C F can take different values, for more details, see [57,58]. In order to solve the closed form of the SGS density weighted transported PDF equation (Equation (12)) the Eulerian stochastic field formulation according to Valińo [35], Jones and Martinez [54], and Mahmoud et al. [4] is applied. In conjunction with the FPV chemistry reduction technique, the temporal evolution of the SGS joint PDF is represented by an ensemble of Ns Eulerian stochastic fields ξ n α (x i , t) defined in the entire domain for each controlling variable α = {PV, f } with 1 ≤ n ≤ Ns and 1 ≤ α ≤ N α . These fields evolve following the equation: which has to be solved in addition to the governing equations of mass and momentum (Equations (6) and (7)) within the LES hybrid ESF/FPV approach. Thereby, the chemical source term . ω α is determined based on the evolution of the controlling variable ξ n α , and the last term, dW n j , which stands for the increments of a vectorial Gaussian process, is referred to as the vector Wiener process, which is spatially uniform, varying in time and different for each stochastic field [59]. Hence, for the different stochastic fields Ns and by employing the approximation used in [35], the Wiener term is determined by the discrete time-step ∆ t between the temporal sample t n and t n+1 , multiplied by the dichotomic vector following (see also in [4]): The Wiener process as a random walk is normally distributed for Ns stochastic fields with zero mean and variance of the time increment ∆ t reading: It is important to note that the different stochastic fields are continuous in time due to the Wiener process. The fields are assumed to be smooth with the grid control volume and then discretized at the length scale of the grid size. This assumption is highly relevant, since sub-grid scales are absent in the stochastic fields and stochastic differential equations (SDE) are fully resolved on the grid size level. Further, it is worth mentioning that solving the different SDE means solving an equivalent set of equations to Equation (12) since they have the same one point-PDF, and the use of Ns stochastic fields for each controlling variable yield to reconstruct the filtered density function equations. Both the filtered mean of the different controlling variable calculated by solving the SDE and the chemical source term can be expressed, respectively, by: .

Numerical Implementation
The reported simulations have been carried out by means of the open source code OpenFOAM [60], in which the hybrid ESF/FPV was implemented [4] and extended in this study to be coupled with LES turbulence models. The solution procedure is carried out through various numerical steps. According to [4], once the time loop starts, N number of stochastic Fields for the two controlling variables are created, which means 2*N new fields are created for both scalars: N fields for the mixture fraction, f : f n with n = {1, . . . , N} and N fields for the progress variable, PV:PV n with n = {1, . . . , N}. The random Wiener term is computed, at the same stage, separately for each individual stochastic field in all spatial directions: dW n j with j = {1, 2, 3} allowing the random procedure of these created fields for both chemical table controlling variables. At this stage, the momentum and continuity transport equations are calculated where all necessary parameters such as density, viscosity and chemical source term for each stochastic field, are extracted from the chemical table as function of calculated f n and PV n .
In the numerical set-up of the studied cases, a stochastic field properties file (in terms of OpenFOAM code) is adopted in the way that the number of stochastic fields can be specified, either N = 1, 2, . . . , or 128. Hence, the stochastic differential equations from ESF/FPV approach read: Since only two controlling variables are addressed for the chemical look-up table, Equation (8) is replaced by the SDE equations in the context of ESF/FPV. In Equation (20), the filtered chemical reaction source term is evaluated by averaging the source terms of different stochastic fields, and there with the sub-grid level fluctuations due to stochastic terms are taken into account.
Following the numerical set-up explained in [4] and keeping the main numerical scheme features, momentum predictor, pressure solver and momentum corrector are applied to solve the compressible Navier-Stokes governing equations (Equations (6) and (7)). A second order central differencing scheme is employed for the convection term of the velocity field, and a second order conservative scheme is applied for the Laplacian terms. A minmod differencing scheme is addressed allowing the solution total variation diminishing (TVD) for scalar flux and the second order backward integration method is used for the time derivative terms for LES calculations. Further details about the discretization procedure and the numerical schemes are explained in the Open-FOAM programmer's guide [60].

Experimental and Simulation Set-Up Details
Two configurations have been selected to assess the prediction capability of the proposed LES hybrid ESF/FPV approach in both air and oxy-enriched environments, and to analyze in the oxy-flame the impact of CO 2 dilution on the general prediction of the oxy-fuel flame structure, stability and emissions. The first configuration studied is the well-known reference case, the air-piloted jet flame, Sandia Flame-D [61,62]. The second configuration features a turbulent non-premixed methane oxy-flame B3 [43].

Air Piloted-Jet Flame (Sandia Flame D)
The burner geometry of the turbulent reacting configuration Sandia Flame-D [62] consists of a central fuel nozzle with diameter d fuel = 7.2 mm (the wall thickness is 0.25 mm) surrounded by a pilot jet burner with d pilot = 18.2 mm (the wall thickness is 0.35 mm). The main jet is constituted by a mixture of 25% methane and 75% air by volume. The pilot jet burns a premixture of acetylene (C 2 H 2 ), hydrogen (H 2 ), air, carbon dioxide (CO 2 ), and nitrogen (N 2 ). It is operated lean. This configuration was planted in an air-wind tunnel around the nozzle presenting the coflow stream with diameter d coflow = 165 mm. Based on the central fuel nozzle properties it exhibited a Reynolds number of 22,400. Various scalar data obtained from Raman/Rayleigh/LIF (Laser-Induced fluorescence) measurement techniques and velocity data from two-component LDV measurements are available [62]. Figure 2 features the computational configuration studied where a 3D structured grid with 40 blocks and 1,808,440 control volumes has been found sufficient to accomplish grid-insensitive solutions. The height of the numerical grid applied was 80 d fuel . Regarding the numerical set-up, the waveTransmissive condition was used for the outlet with fixed pressure value; P = 1 atm, and zero gradient boundary conditions were set for all other parameters. All the operating conditions are provided in Table 1. To assess the convergence of the number of stochastic fields (SFi), the LES was carried out for different numbers of stochastic fields (4,6,8,16 SFi) for 0.38 s of physical time, respectively, to ensure that the flow reaches the statistically steady state with time step ∆t = 3 × 10 −7 . The calculations were accomplished on 96 processors. techniques and velocity data from two-component LDV measurements are available [62]. Figure 2 features the computational configuration studied where a 3D structured grid with 40 blocks and 1,808,440 control volumes has been found sufficient to accomplish grid-insensitive solutions. The height of the numerical grid applied was 80 dfuel. Regarding the numerical set-up, the waveTransmissive condition was used for the outlet with fixed pressure value; P = 1 atm, and zero gradient boundary conditions were set for all other parameters. All the operating conditions are provided in Table 1. To assess the convergence of the number of stochastic fields (SFi), the LES was carried out for different numbers of stochastic fields (4,6,8,16 SFi) for 0.38 s of physical time, respectively, to ensure that the flow reaches the statistically steady state with time step Δt = 3 × 10 −7 . The calculations were accomplished on 96 processors.

Oxy-Fuel Jet Flame (Flame B3)
The oxy-fuel jet flame B3 is one of the flame series that were experimentally studied in [43]. It exhibits a Reynolds number Re = 18,000, making it interesting for LES investigations. The reacting turbulent burner is composed mainly of two jets; a central jet with an inner diameter dB3-fuel = 5 mm (the wall thickness is 0.5 mm) for the fuel stream with a mixture of 55% H2 and 45% CH4, and a coaxially arranged jet with diameter dB3-oxy = 96.5 mm for the premixed oxidizer composition of 68% CO2 and 32% O2. In fact, it has been reported in [63,64], that air flame stability is established once the O2% content is around 30% which explains the specie percentages used for the oxidizer stream in this study. Further, the existence of H2 in the fuel jet helps the flame to remain attached to the nozzle as indicated in [43]. Figure 3 presents the 3D configuration of flame B3 featuring the computational

Oxy-Fuel Jet Flame (Flame B3)
The oxy-fuel jet flame B3 is one of the flame series that were experimentally studied in [43]. It exhibits a Reynolds number Re = 18,000, making it interesting for LES investigations. The reacting turbulent burner is composed mainly of two jets; a central jet with an inner diameter d B3-fuel = 5 mm (the wall thickness is 0.5 mm) for the fuel stream with a mixture of 55% H 2 and 45% CH 4 , and a coaxially arranged jet with diameter d B3-oxy = 96.5 mm for the premixed oxidizer composition of 68% CO 2 and 32% O 2 . In fact, it has been reported in [63,64], that air flame stability is established once the O 2 % content is around 30% which explains the specie percentages used for the oxidizer stream in this study. Further, the existence of H 2 in the fuel jet helps the flame to remain attached to the nozzle as indicated in [43]. Figure 3 presents the 3D configuration of flame B3 featuring the computational domain in a conical shape assuming that the base and the sides are inlets for the oxidizer, which helps the numerical convergence of the calculations. The block-structured hexahedral grid applied consists of 89 blocks and 2,800,770 cells in an O-ring arrangement. 8450 cells are used for the fuel pipe jet, with minimum cell size equal to 5.31 × 10 −12 m 3 and maximum cell size located in the outlet region with 1.09 × 10 −10 m 3 . To assess the LES quality, various criteria are available, as reported in [46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65]. In the present paper we adopt the approach proposed by Hanjalic et al. in [66] and Pope in [7]. Thus, based on the ratio of the LES local mesh size ∆ to the Kolmogorov length scale η, it turns out that ∆/η ≤ 12 in almost the entire computational domain, which reflects that more than 80% of the turbulent kinetic energy are resolved, along with the Taylor scales. Thereby, the Kolmogorov length scale, that has to be estimated by RANS solution following the procedure in [67], has been evaluated by RANS calculations from [4]. The outlet plane is at x = 80 d B3-fuel , with waveTransmissive condition and fixed pressure P = 1 atm. Similar to the previous case, zero gradient boundary conditions are imposed for other variables and the temperature is initially uniformly distributed. More details for the operating conditions are summarized in Table 2. To evaluate the convergence of the number of stochastic fields, SFi, calculations for different cases (using 4,8,16 SFi) were also run in this case, however on 160 processors for 0.16 s of physical time with time step ∆t = 2.5 × 10 −7 . For the inlet turbulent flow field velocity of both studied cases, an in house turbulence inflow generator according to [68] was used. other variables and the temperature is initially uniformly distributed. More details for the operating conditions are summarized in Table 2. To evaluate the convergence of the number of stochastic fields, SFi, calculations for different cases (using 4,8,16 SFi) were also run in this case, however on 160 processors for 0.16 s of physical time with time step Δt = 2.5 × 10 −7 . For the inlet turbulent flow field velocity of both studied cases, an in house turbulence inflow generator according to [68] was used.

Validation on an Air-Piloted Jet Flame (Sandia Flame D)
In this section, the Sandia flame D is used to support the evaluation of the performance of the LES hybrid FESF/FPV. This is achieved in terms of comparison between numerically obtained results with experimental data for various radial profiles of some scalar quantities, scatter plots and probability density function at three axial positions, z/d = 7.5; 15 and 30. Thereby, different simulations using different numbers of Stochastic Fields (SFi = 4,6,8 and 16) are presented, in order to assess the convergence of the number of SFi. Time mean and root mean square (rms) of both axial and radial velocity are displayed in Figure 4. For the axial velocity, the obtained radial profile results feature a good overall agreement with experimental data at different positions for both mean (solid line) and rms quantities (dashed line). However, some under-estimations are noticed in reproducing the behavior of the radial velocity at z/d = 15 while an over-prediction is observed at z/d = 30. By applying various numbers of stochastic fields, it turns out that the mean quantities of the radial velocity obtained with 8SFi are predicted better at all axial positions.

Validation on an Air-Piloted Jet Flame (Sandia Flame D)
In this section, the Sandia flame D is used to support the evaluation of the performance of the LES hybrid FESF/FPV. This is achieved in terms of comparison between numerically obtained results with experimental data for various radial profiles of some scalar quantities, scatter plots and probability density function at three axial positions, z/d = 7.5; 15 and 30. Thereby, different simulations using different numbers of Stochastic Fields (SFi = 4, 6, 8 and 16) are presented, in order to assess the convergence of the number of SFi. Time mean and root mean square (rms) of both axial and radial velocity are displayed in Figure 4. For the axial velocity, the obtained radial profile results feature a good overall agreement with experimental data at different positions for both mean (solid line) and rms quantities (dashed line). However, some under-estimations are noticed in reproducing the behavior of the radial velocity at z/d = 15 while an over-prediction is observed at z/d = 30. By applying various numbers of stochastic fields, it turns out that the mean quantities of the radial velocity obtained with 8SFi are predicted better at all axial positions. Figure 5 depicts the radial profiles of temperature and mixture fraction at different axial positions. Regarding the mixture fraction evolution, numerical results indicate an excellent agreement for both mean and rms quantities at the axial position z/d = 7.5, and a slight under-estimation at the middle axial location further from the centerline without exhibiting any differences once using varied stochastic field numbers. However, the mean temperature at the centerline was slightly over-predicted at both positions z/d = 7.5 and 15 in Figure 5a, independently from the number of stochastic fields. This over-prediction was mostly pronounced at the centerline only at the further downstream positions, which is explained by the under-estimation of mixture fraction at same axial position. The influence of using different stochastic field numbers is revealed by comparing the rms quantities of the temperature. It appears clearly that the use of 8SFi allows for an overall better prediction. Such a SFi number agrees with that found by Jones et al. [31]. Further, the reported results proved the efficiency of the suggested novel method once compared to available measurements and to other previous numerical studies in which different flamelet models coupled to the FPV approach have been adopted. In fact, regarding the different computational setups employed in the case of Sandia Flame-D, Coclite et al. in [69] used the statistically most likely distribution model (SMLD) coupled to the FPV technique and found the mixture fraction profiles at different axial positions rather slightly under-predicted near the centerline. In [15], the mean temperature profiles were slightly over-predicted in some axial positions far from the centerline.   Figure 5 depicts the radial profiles of temperature and mixture fraction at different axial positions. Regarding the mixture fraction evolution, numerical results indicate an excellent agreement for both mean and rms quantities at the axial position z/d = 7.5, and a slight underestimation at the middle axial location further from the centerline without exhibiting any differences once using varied stochastic field numbers. However, the mean temperature at the centerline was slightly over-predicted at both positions z/d = 7.5 and 15 in Figure 5a, independently from the number of stochastic fields. This over-prediction was mostly pronounced at the centerline only at the further downstream positions, which is explained by the under-estimation of mixture fraction at same axial position. The influence of using different stochastic field numbers is revealed by comparing the rms quantities of the temperature. It appears clearly that the use of 8SFi allows for an overall better prediction. Such a SFi number agrees with that found by Jones et al. [31]. Further, the reported results proved the efficiency of the suggested novel method once compared to available measurements and to other previous numerical studies in which different flamelet models coupled to the FPV approach have been adopted. In fact, regarding the different computational setups employed in the case of Sandia Flame-D, Coclite et al. in [69] used the statistically most likely distribution model (SMLD) coupled to the FPV technique and found the mixture fraction profiles at different axial positions rather slightly under-predicted near the centerline. In [15], the mean temperature profiles were slightly over-predicted in some axial positions far from the centerline.  Figure 6 shows two series of instantaneous scatter plots of the temperature distribution at 3 axial positions. The first series exhibit the calculated instantaneous scatter plots plotted versus the instantaneous mixture fraction quantities, when 8 SFi are used (Figure 6, bottom). The second series present the reference instantaneous scatter plots (Figure 6, top). It turns out that the maximum  Figure 6 shows two series of instantaneous scatter plots of the temperature distribution at 3 axial positions. The first series exhibit the calculated instantaneous scatter plots plotted versus the instantaneous mixture fraction quantities, when 8 SFi are used (Figure 6, bottom). The second series present the reference instantaneous scatter plots (Figure 6, top). It turns out that the maximum temperature corresponding to the stochiometric mixture fraction is well predicted for all locations except the first axial position where it is very slightly overestimated with small deviations. Some local extinction values could not be predicted numerically in the region close to the centerline. However, one can note the presence of few unrealistic structures of temperature in the mixture fraction range 0.8 < f < 0.9 for both axial positions z/d = 7.5 and 15. This could shape up due to the standard mixing model used in this work. The Interaction by exchange of mean (IEM) is known, by nature, to be not local in composition space which may lead to unphysical events (see in [70]). As reported in [57] where two different mixing models (IEM and EMST (Euclidean minimum spanning tree)) have been tested on all piloted Sandia flames including the Flame-D, the IEM approach could not capture perfectly the conditional temperature distribution leading to a similar unphysical trend within the same mixture fraction range as observed also in [58]. Another quantity which is essential to evaluate numerical combustion models in predicting the flame characteristics is the probability density function (PDF). Figure 7 displays the PDF of the temperature at 2 axial locations further downstream from the flame nozzle obtained numerically using different SFi numbers in comparison with experimental data. This PDF was calculated for a narrow band of the mixture fraction around the stoichiometry where 0.30 < f < 0.40. The PDFs obtained calculated using 8 SFi are showing better results comparing to other numerical results simulated with less fields, where very similar behavior to the available experimental data is illustrated. In particular, the probability of the presence of temperature ranging from 1650 K to 2150 K is comparable to the measurements at both positions, z/d = 15 and 30. Next, the novel LES hybrid ESF/FPV model was applied to investigate the evolution of oxy-fuel flame characteristics. Another quantity which is essential to evaluate numerical combustion models in predicting the flame characteristics is the probability density function (PDF). Figure 7 displays the PDF of the temperature at 2 axial locations further downstream from the flame nozzle obtained numerically using different SFi numbers in comparison with experimental data. This PDF was calculated for a narrow band of the mixture fraction around the stoichiometry where 0.30 < f < 0.40. The PDFs obtained calculated using 8 SFi are showing better results comparing to other numerical results simulated with less fields, where very similar behavior to the available experimental data is illustrated. In particular, the probability of the presence of temperature ranging from 1650 K to 2150 K is comparable to the measurements at both positions, z/d = 15 and 30. Next, the novel LES hybrid ESF/FPV model was applied to investigate the evolution of oxy-fuel flame characteristics. obtained calculated using 8 SFi are showing better results comparing to other numerical results simulated with less fields, where very similar behavior to the available experimental data is illustrated. In particular, the probability of the presence of temperature ranging from 1650 K to 2150 K is comparable to the measurements at both positions, z/d = 15 and 30. Next, the novel LES hybrid ESF/FPV model was applied to investigate the evolution of oxy-fuel flame characteristics.

Application to Oxy-Fuel Jet Flame B3
The LES hybrid ESF/FPV was now appraised on the oxy-flame B 3 in terms of comparison with available measurements, and then applied to analyze the effect of CO 2 dilution on the oxy-flame behavior. First, Figure 8 shows qualitatively both experimental and numerical instantaneous temperature field for oxy-flame B3. Thereby, the region close to the nozzle of the burner is particularly zoomed in order to clearly display the measurement planes. The predicted flame length observed in Figure 8c agrees reasonably well with the experimental data, which reached approximately L flame = 100 d B3-fuel (see in Figure 8a). The LES hybrid ESF/FPV was now appraised on the oxy-flame B3 in terms of comparison with available measurements, and then applied to analyze the effect of CO2 dilution on the oxy-flame behavior. First, Figure 8 shows qualitatively both experimental and numerical instantaneous temperature field for oxy-flame B3. Thereby, the region close to the nozzle of the burner is particularly zoomed in order to clearly display the measurement planes. The predicted flame length observed in Figure 8c agrees reasonably well with the experimental data, which reached approximately Lflame = 100 dB3-fuel (see in Figure 8a).  Figure 8b shows the instantaneous temperature contour on radial planes at different axial positions; it can be seen that the temperature field spreads radially and toward the outlet of the burner. The flame is also shown attached to the burner nozzle as it is observed experimentally under consideration of the important amount of diluted H2 in the fuel mixture. In the zoomed region above the fuel nozzle (Figure 8c right), the temperature field reaches the maximum predicted value, which corresponds to 2250 K found in experiments.
The quantitative validation was achieved by comparisons with experimental data in Figures 9-12 of the evolution of main species and temperature, as well as of the scatter plots and probability density function of scalars at different axial positions. To ensure the convergence of the number of stochastic fields (SFi) various simulations are performed by applying different numbers of Stochastic  Figure 8b shows the instantaneous temperature contour on radial planes at different axial positions; it can be seen that the temperature field spreads radially and toward the outlet of the burner. The flame is also shown attached to the burner nozzle as it is observed experimentally under consideration of the important amount of diluted H 2 in the fuel mixture. In the zoomed region above the fuel nozzle (Figure 8c right), the temperature field reaches the maximum predicted value, which corresponds to 2250 K found in experiments.
The quantitative validation was achieved by comparisons with experimental data in Figures 9-12 of the evolution of main species and temperature, as well as of the scatter plots and probability density function of scalars at different axial positions. To ensure the convergence of the number of stochastic fields (SFi) various simulations are performed by applying different numbers of Stochastic Fields (4, 8 and 16). Figure 9 shows the evolution of the mean values of CO and H 2 O species mass fractions at different axial positions. The amounts of CO and H 2 O mass fractions were much higher than values obtained with air diluted flame cases. Especially, close to the nozzle, the CO and H2O mass fractions were equal to 0.11 and 0.18, respectively. However, these values increased to almost 0.2 for both species at higher axial positions. In this way, the present LES results confirm the RANS-based findings of Rihab et al. in [4], which agree with Masri et al. in [71], who reported that the CO 2 68% diluted with pure oxygen in the oxidizer inlet is not inert. This plays an important role in increasing CO and H 2 O gases formation since the additional CO 2 reacts with H to form CO. Figure 10 presents the mean quantities of O 2 and H 2 species mass fractions. In both figures, the results with 16 SFi achieved better agreements at radial profiles once compared to results with 4 SFi. This is not surprising, since four stochastic fields are not sufficient to reproduce the exact evolution of measured quantities (see also previous case: Sandia flame-D results using 4SFi).  In general, the results of CO and H2O mass fraction at axial positions close to the nozzle using 8 SFi are shifted in Figure 9 toward the axial centerline. Further downstream, the agreement with experimental data is excellent. Similar observations can be made in Figure 10 for the results of O2 and H2 species mass fractions with 8 SFi. While slight under-prediction of O2 at z/d = 10 was recorded, the use of different SFi numbers had a negligibly small impact on the H2 mass fraction evolution.
The instantaneous scatter plots of the temperature distribution at axial position near the fuel nozzle (z/d = 3) were compared with available experimental data in Figure 11. Thereby, the vertical black line drawed in both plots corresponds to the stochiometric mixture fraction value. The maximum reached value of the temperature distribution was 2100 K, which matches the experimental data. All values referring to the extinction locations were properly reproduced numerically. However, the lowest value predicted with the simulation reads 750 K, while experimentally it is 500 K. Furthermore, the mean temperature profile is lightly overestimated near the stoichiometry region comparing to experimental results.  Figure 12 compares the PDF calculated using different numbers of SFi with those obtained from measurements at axial positions z/d = 3; 5 and 10. For both sets, the PDF was estimated for a narrow band of the mixture fraction around the stoichiometry where ∆ = 0.04. In most locations, the shape of the PDF of the temperature could be reproduced only with simulation results using 16 SFi, but slight differences were noticed at both last axial positions where few values of temperature in the range of 1500 K and 1700 K were over-predicted numerically. Nevertheless, the probability of the presence of the temperature range from 1800 K to 2100 K is matching the values obtained experimentally for all positions in a good way. In general, the results of CO and H2O mass fraction at axial positions close to the nozzle using 8 SFi are shifted in Figure 9 toward the axial centerline. Further downstream, the agreement with experimental data is excellent. Similar observations can be made in Figure 10 for the results of O2 and H2 species mass fractions with 8 SFi. While slight under-prediction of O2 at z/d = 10 was recorded, the use of different SFi numbers had a negligibly small impact on the H2 mass fraction evolution.
The instantaneous scatter plots of the temperature distribution at axial position near the fuel nozzle (z/d = 3) were compared with available experimental data in Figure 11. Thereby, the vertical black line drawed in both plots corresponds to the stochiometric mixture fraction value. The maximum reached value of the temperature distribution was 2100 K, which matches the experimental data. All values referring to the extinction locations were properly reproduced numerically. However, the lowest value predicted with the simulation reads 750 K, while experimentally it is 500 K. Furthermore, the mean temperature profile is lightly overestimated near the stoichiometry region comparing to experimental results.  Figure 12 compares the PDF calculated using different numbers of SFi with those obtained from measurements at axial positions z/d = 3; 5 and 10. For both sets, the PDF was estimated for a narrow band of the mixture fraction around the stoichiometry where ∆ = 0.04. In most locations, the shape of the PDF of the temperature could be reproduced only with simulation results using 16 SFi, but slight differences were noticed at both last axial positions where few values of temperature in the range of 1500 K and 1700 K were over-predicted numerically. Nevertheless, the probability of the presence of the temperature range from 1800 K to 2100 K is matching the values obtained experimentally for all positions in a good way. In general, the results of CO and H 2 O mass fraction at axial positions close to the nozzle using 8 SFi are shifted in Figure 9 toward the axial centerline. Further downstream, the agreement with experimental data is excellent. Similar observations can be made in Figure 10 for the results of O 2 and H 2 species mass fractions with 8 SFi. While slight under-prediction of O 2 at z/d = 10 was recorded, the use of different SFi numbers had a negligibly small impact on the H 2 mass fraction evolution.
The instantaneous scatter plots of the temperature distribution at axial position near the fuel nozzle (z/d = 3) were compared with available experimental data in Figure 11. Thereby, the vertical black line drawed in both plots corresponds to the stochiometric mixture fraction value. The maximum reached value of the temperature distribution was 2100 K, which matches the experimental data. All values referring to the extinction locations were properly reproduced numerically. However, the lowest value predicted with the simulation reads 750 K, while experimentally it is 500 K. Furthermore, the mean temperature profile is lightly overestimated near the stoichiometry region comparing to experimental results. Figure 12 compares the PDF calculated using different numbers of SFi with those obtained from measurements at axial positions z/d = 3; 5 and 10. For both sets, the PDF was estimated for a narrow band of the mixture fraction around the stoichiometry where ∆ = 0.04. In most locations, the shape of the PDF of the temperature could be reproduced only with simulation results using 16 SFi, but slight differences were noticed at both last axial positions where few values of temperature in the range of 1500 K and 1700 K were over-predicted numerically. Nevertheless, the probability of the presence of the temperature range from 1800 K to 2100 K is matching the values obtained experimentally for all positions in a good way.

Conclusions
In this paper, it was demonstrated that the combustion characteristics of a turbulent oxy-methane non-premixed flame operating under highly diluted conditions of CO 2 and H 2 in oxidizer and fuel streams, respectively, can be predicted and analyzed using a novel hybrid ESF/FPV model designed for that purpose. It combines a transported joint scalar probability density function (T-PDF) following the Eulerian stochastic field methodology (ESF), and a flamelet progress variable (FPV) turbulent combustion model under consideration of a detailed chemical reaction mechanism. This novel hybrid ESF/FPV approach removes the weaknesses of the presumed-probability density function (P-PDF)-based FPV modeling along with the solving of associated additional modeled transport equations while rendering the T-PDF computationally less affordable.
The prediction capability of the LES hybrid ESF/FPV was first appraised successfully on an air-piloted methane jet flame (Sandia Flame D). Then, it was assessed in analyzing the combustion properties of an oxy-flame, the so-called oxy-flame B 3 under investigation and in capturing CO 2 effect on the behavior of the oxy-flame. These configurations were previously investigated by the authors using a RANS-based formulation of this model. Relying on the results achieved, the following conclusions can be listed:

•
Overall, the LES hybrid filtered ESF/FPV approach demonstrated its high capability in capturing the main flame characteristics and flow field variables not only in the Sandia flame D but also in an oxy-flame configuration with CO 2 and H 2 dilution in oxidizer and fuel streams, respectively, using FPV tables based on Le = 1. This indicates that the H 2 induced differential diffusion effect is not important in the investigated area of this oxy-flame B 3 . • Compared to our RANS results in our previous work [4], the LES hybrid ESF/FPV model provides more accurate predictions. • A good convergence of the optimal number of stochastic fields (SFi) strongly depends on the complexity of the combustion case. Even though more SFi could help to achieve fully complete convergence in this regard, in search of saving computational costs it turned out that a simulation with at least 8 stochastic fields leads to an accurate prediction in Sandia flame D as in [31] and at least 16 stochastic fields allow achieving better results in accordance with measurements once complex configuration like the oxy-fuel flame-B 3 is investigated.

•
Regarding the general prediction of the oxy-flame structure, stability and emissions, it turns out that 68% molar percentage of additional CO 2 enrichment in the oxidizer side leads to 0.39% of CO formation near the burner fuel nozzle and 0.62% at 10 d fuel above the nozzle. These amounts of CO-formed gases are clearly high compared to ordinary flame cases with air/fuel conditions with 0.35% near the burner nozzle and 0.62% at 10 d fuel above the jet tip.
The novel model needs to be further evaluated by using FPV tables with a non-unity Lewis number and by considering advanced micro-mixing models coupled to the LES hybrid ESF/FPV approach. All of these tasks are left for future work.

Conflicts of Interest:
The authors declare no conflict of interest.