University of Birmingham Numerical Investigation of an Oxyfuel non-premixed combustion using a Hybrid Eulerian Stochastic Field/Flamelet Progress Variable approach: Effects of H2 / CO2 enrichment and Reynolds number

In the present paper, the behaviour of an oxy-fuel non-premixed jet flame is numerically investigated by using a novel approach which 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 detailed chemical reaction mechanism. This hybrid ESF/FPV approach overcomes the limitations of the presumedprobability density function (P-PDF) based FPV modelling along with the solving of associated additional modelled transport equations while rendering the T-PDF computationally less demanding. In Reynolds Averaged Navier-Stokes (RANS) context, the suggested approach is first validated by assessing its general prediction capability in reproducing the flame and flow properties of a simple piloted jet flame configuration known as Sandia Flame D. Second, its feasibility in capturing CO2 addition effect on the flame behaviour is demonstrated while studying a non-premixed oxy-flame configuration. This consists of an oxy-methane flame characterized by a high CO2 amount in the oxidizer and a significant content of H2 in the fuel stream, making it challenging for combustion modelling. Comparisons of numerical results with experimental data show that the complete model reproduces the major properties of the flame cases investigated and allows achieving the best agreement for the temperature and different species mass fractions once compared to the classical presumed PDF approach.


Introduction
High temperature processes are an essential element of advanced industrial societies. They dominate energy, metal and materials production processing. In the most of these industries, the huge amount of energy needed is currently covered by the combustion of fossil fuels. However, this classical combustion In the Eulerian framework, two approaches are being followed. Fox in reference [18] suggested the so-called multi-environment PDF method in which the joint PDF is solved by a set of deterministic Eulerian transported equations, in contrast to the stochastic particle method. Another working frame is the so-called Eulerian stochastic fields (ESF) method proposed by Valińo [16] and Sabel'nikov et al. [19].
Here, the solution of the stochastic differential equation in turbulent reacting flow is obtained by using smooth different stochastic fields for scalars to describe the PDF undergoing diffusion, turbulent convection and chemical reaction. Recently, the ESF method has been successfully used in a wide range of combustion problems, including laboratory gas jet flames like in references [3,13,20,21]. In these studies, the ESF method is essentially coupled to reduced chemical mechanisms. Recently, Jangi et al. [22][23][24], applied the ESF method in conjunction with detailed chemistry for some jet flames cases.
Only very few studies focused on developing a framework for the ESF method in coupling with a tabulated chemistry. In their large eddy simulation (LES) study, Amer et al. [25] reported a coupling of the ESF approach with a FGM tabulated chemistry using the FASTEST code.
To the knowledge of the authors, no numerical studies have addressed the FPV tabulated chemistry strategy coupled to the ESF method. As pointed out by [26] who compared the two different flamelet reduced order manifolds (FGM and FPV) for non-premixed flames, both FPV and FGM operate differently. The generated tables for FGM and FPV show similar behaviour in the mixture fraction space but different behaviour on the progress variable space. For more details the reader is referred to [26].
In the present work, the Eulerian Stochastic Field method is coupled to the FPV approach. The test and the validation of the method will be carried out by applying, RANS, the famous Sandia flame D experimentally investigated in reference [27]. The method will be then assessed in capturing the flame properties along with the CO 2 dilution effect on the flame behaviour of a non-premixed oxy-methane jet flame asexperimentally studied in reference [28]. It is characterized by a high CO 2 amount in the oxidizer and a significant content of H 2 in the fuel stream. In this environment, the flame becomes more prone to extinction, so that extra stabilization mechanisms for example by enriching the fuel stream with H 2 , appear mandatory. This results in a complex flow and mixing structure behaviour making it challenging for combustion modelling. Three aspects are especially investigated in the present paper, namely the impact of the addition of H 2 by comparing its enrichment in flames A1 and A3, the effect of high level dilution of CO 2 in oxidizer stream for all cases and the influence of the Reynolds number by comparing cases A1 and B3 from [28]. In fact, based on the experimental study by Sevault et al., and different studies reporting the carbon capture and storage technique, using O 2 /CO 2 mixtures instead of air for fuel combustion ideally produces water in exhaust gases, which can be easily separated by condensation and pure CO 2 that can be captured and stored. Following previous studies in references [29,30], it has been reported "that the molar percentage of oxygen in the oxidant should be around 30% to reach air flame stability." Thereby in the experimental study on which the numerical work is based, for the oxidizer mixture, the oxygen percentage is set to 3% and consequently 68% of CO 2 is diluted.
Garmory and Mastorakos have used, in Reference [31], the Conditional Moment Closure (CMC) to numerically study the oxyflame configuration A from reference [28], where LES with the Smagorinsky sub grid scale model was employed. In the present paper, rather a RANS framework usually applied in industrial environments is favoured.
The paper is organized as follows: Section 2 outlines the modelling approach including the RANS model used along with the hybrid ESF/FPV method. Section 3 introduces the experimental configurations for both the flame D and the Oxyfuel jet flame with different numerical setups. The achieved results are presented and discussed in Section 4. Concluding remarks are summarized in the last Section.

Reynolds Averaged Transport Equations
A turbulent reacting Newtonian fluid flow is investigated. For its description in the frame of the FPV approach (see Section 2.2), the required averaged equations are expressed in a conservative form for mass, momentum and scalars in RANS context. The turbulence is modelled by means of the standard two equations k-epsilon model. Thereby, the following set of equations is then solved: In Equations (1)-(5), (.) stands for a Favre weighted quantity and (.) for mean quantity. ρ is the density, U i the velocity component of the flow field in different directions (i = 1, 2, 3), p the pressure, µ the dynamic molecular viscosity and µ t = ρC µ k 2 / ε the dynamic turbulent viscosity. The coefficient σ t represents the turbulent Schmidt number and G k stands for the production of turbulence kinetic energy (G k = 2µ t S ij S ij , with S ij being the strain rate tensor).
According to [32] the turbulence model coefficients are set as: The remaining quantity . ω Y k stands for the species source term which will be determined from the FPV tables.

Chemistry Modelling Using Flamelet/Progress Variable FPV Approach
Dealing with turbulent reacting flow, the problem of turbulent combustion model and simulation dwells upon finding a suitable way to determine the mean chemical source term . ω Y k using appropriate chemical reaction model or mechanism and then solving Equation (3). According to the FPV formulation, a tabulation based on two controlling variables, namely the mixture fraction and a reaction progress variable, is generated instead of using the scalar dissipation rate (strain rate) in the classical Steady Laminar Flamelet model (see in references [33][34][35][36]). A two-scalar representation emerges, which results in look-up tables with coordinates associated with the controlling variables. All relevant information (e.g., the chemical composition, chemical rates of formation/destruction, temperature, density . . . ) is therefore stored in the tables as function of the mixture fraction f and the reaction progress variable (PV) only according to [37,38]. 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. The FPV tables are created after solving in physical space laminar counter-flow flames using the flamelet generator Flame Master code [39]. In the present study, two definitions of the progress variable are used. In the case of the Sandia Flame-D: and for the Oxyfuel flame cases: PV = Y H2O W H2O . The chemical tables created for the two configurations, Sandia flame-D and Oxyfuel flames (A1, A3 and B3), are generated from steady non-premixed flamelets at different strain rate, from very small to extinguishing values and the state of the unsteady flamelet at the extinguishing strain rate. Figure 1 shows some features from the FPV tables for different flames cases. It displays the distribution of the production rate of heat release coloured by the source term of the progress variable in the mixture fraction space. As both fuel and oxidizer compositions are different from one case to another, the maximum value of heat release production rate is varying. In Figure 1a referring to Flame A1 and B3 with higher H 2 % enrichment, the maximum value of heat release locally reaches 7.10 9 J/m 3 s. It is reduced to 5.10 9 J/m 3 s once the H 2 % enrichment is decreased to 37% in Flame A3 as shown in Figure 1b. Assuming a unity Lewis number even in the Oxyfuel case, the effect of hydrogen high diffusivity is not considered. The detailed chemical mechanism used consists in 325 reactions and 53 species available in GRI-MECH 3.0 [40].

Modelling of Turbulence-Chemistry Interaction
Since a complete statistical modelling of turbulent flames includes both turbulence closure and a chemical reaction model which interact, it is of relevance to describe how the turbulence chemistry interaction is treated. Two approaches will be applied in this study, namely the T-PDF according to the Eulerian Stochastic Field approach (ESF) and the presumed probability density function based on the beta-function (presumed β-PDF).

Joint Probability Density Function and Eulerian Stochastic Field Method
The novel turbulence-chemistry interaction model referring to as the hybrid ESF/FPV method is employed and dynamically linked to a customized solver for combustion. For a specific instance of time t and a fixed point x i , the species Φ α involved in the reaction process can be represented by a marginal probability density function P α as explained in references [14,15,41], which can be expressed as follows: where δ represents the Dirac delta function and ψ represents the composition space of species. The joint PDF can then read: The joint PDF in Equation (7) is determined as the product of the one-point fine-grained PDF P α of each species Φ a . According to [42,43], the explicit transport equation for the density weighted PDF reads: where P ψ; x i , t , stands for the density weighted joint PDF defined as: P α ψ ≡ ρP α ψ /ρ. The first term {1} in Equation (8) represents the temporal evolution of the equation in the physical space, the convection due to the mean velocity is represented by the second term {2} and the last term in the left side of the {3} expresses the closed chemical production source term in the sample space. All these three terms in the left side of the equation are in closed form; however the two terms in the right side are unclosed and need modelling. The fourth term {4} in Equation (6) describes the turbulent transport of the PDF for which a gradient approach is used [41]. Usually for the term {5} representing the micro-mixing part of the PDF, an interaction by Exchange with the Mean model (IEM), reported in references [44,45], is adopted as closure. Thereby and after applying the assumption of equal diffusivity, the modelled expression for the filtered joint probability density function P can be obtained in the following form with τ = k/ ε and C Φ is the model coefficient which, according to [25], is set to C Φ = 2. Toward solving the joint probability density function represented in Equation (9) and based on a recent Monte Carlo formulation determined by the Eulerian stochastic field (ESF) method, proposed independently by Valińo in reference [16] and applied in reference [21], the P ψ; x i , t in Equation (9) is approximated for N α scalars using an ensemble of N fields ξ n α (x i , t) defined in the entire domain for 1≤ n ≤ N and 1 ≤ α ≤ N α . Thus, the modelled stochastic differential equation (SDE) can be obtained according to [16,21] in the following form Instead of solving Equation (3), the ESF method is employed in combination with the FPV chemistry reduction technique. Thereby, the ensemble of N Eulerian stochastic fields ξ n α (x i , t), are defined for each controlling variable α ≡ {PV, f } which spans the FPV chemistry table. Hence, in Equation (10), the chemical source term . ω α , is determined based on the evolution of the controlling variables. It is important also to mention here, that solving the different SDE (Equation (10)) for N stochastic fields is equivalent to solve a set of Equation (9), since they have the same one point-PDF and that the use of N stochastic fields for each controlling variable can be used to reconstruct the filtered density function equations.
The filtered mean of the controlling variable obtained by solving the SDE and the mean chemical source term defined from the table for different stochastic fields, are obtained respectively as: .
The last term in Equation (10) accounts for the stochastic component where dW n j stands for the increments of a vectorial Gaussian process referred to as vector Wiener process which is spatially independent, varying in time and different for each stochastic field as reported in reference [46]. For the different stochastic fields N and by employing the approximation used in reference [16], the Wiener term is determined by the discrete time-step ∆ t between the instances t p and t p+1 , multiplied by the dichotomic vector and that is: The Wiener process is essentially normally distributed with zero mean and variance of the time increment ∆ t reading: In the Eulerian stochastic method (see Valińoin reference [16]), the different stochastic fields are continuous in time due to the random walk of Wiener process where its components are uncorrelated for different time levels. Thereby the fields are smooth with the cell mesh volume and then discretized at the length scale of the grid size. These assumptions are highly required in this method in order to fully solve all stochastic differential equations (SDE) on the grid size level.

Presumed Probability Density Function Approach
Setting up the composition space by the two controlling variables within combined Presumed-PDF/FPV model, an integration of a scalar field ϕ α with the PDF over the composition space yieldsall statistical moments, for example the mean and the variance, as given by where P is unknown and has to be prescribed. Assuming a statistical independence between the mixture fraction, f and the reaction progress variable, PV, in this study as usually within the classical FPV approach [47], a joint PDF based on presumed β-function is used for the mixture fraction and a Dirac function is applied for the reaction PV. Because the shape of the β-PDF depends on the weighted mean and variance of the mixture fraction f , f "2 , transport equations for these variables in Equations (16) and (17) are solved in addition to the transport equation of PV (Equation (18)). For practical convenience, the controlling variable PV is normalized by its maximum value PVmax, which is a function of the mixture fraction f. Therefore, together Equations (1), (2), (4) and (5) and the following Equations (e.g., [48]) where Equation (3) is replaced by Equation (18), will be solved in the CFD code. To determine the chemical source term in Equation (18), we set ϕ = .
ω and read . ω from the FPV table. The needed mean chemical source term is then computed via the integration described in Equation (15).

Presumed PDF/FPV
Subsequent to the implementation of the FPV chemical model within the PDF approach, the mean chemical source term and the mean density will be fed back (interpolated) to the CFD code. In practice, the numerical procedure consists of an iterative exchange of the turbulent flow and mixing field, the mean density and the chemical source term between the CFD code and the PDF sub-code. The mean values and variances of the controlling variables gained from the solution of their transport Equations (16)- (18) are utilized as input parameters for the look-up table interpolation.
This table interpolation provides, in turn, the quantities required in the transport equations of the controlling variables to be used in the next iteration step. The whole PDF integration is carried out as a pre-processing step. By using 401 equidistant discretisation points for the mixture fraction, 401 for the progress variable and10 for the variance of the mixture fraction and once the discretisation is accomplished a large look-up table that requires a memory of approximately 407 Mbytes is created for each flame.

Eulerian Stochastic Field/FPV
New libraries were created in the open source code OpenFOAM, [49] in order to cope the Eulerian Monte Carlo stochastic field method with the FPV chemistry reduction technique. The solution procedure is obtained through various numerical steps. Heretofore, 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 f : f n with n = {1, . . . , N} and N fields for 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. In this stage, the momentum and continuity transport equations are calculated where all necessary parameters such as density, viscosity and the chemical source term for each stochastic field, are extracted from the chemical table in 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 N = 128. Hence, the stochastic differential equations from ESF/FPV approach read: Within the solution procedure, the chemical reacting source term and different scalars are determined after updating the averaging of the different source terms and scalars from the different calculated stochastic fields ξ n α by accessing N times the chemistry FPV table according to Equations (11) and (12).
In order to solve the governing Equations (1), (5) and (19), (20), which are spatially discretized with a Finite Volume Method (FVM), momentum predictor, pressure solver and momentum corrector are utilized sequentially. Second order central differencing scheme is applied for the convection term of the velocity field. In the case of passive scalar flux, a Minmod differencing scheme is applied to make the solution total variation diminishing (TVD). A second order conservative scheme is applied for the Laplacian terms. The first order Euler integration method is used for the time derivative terms for RANS. For further details about the discretization procedure and the numerical schemes, the reader is referred to Open-FOAM programmer's guide [49].

Experimental Configurations and Case Setups
In this work, two different configurations are numerically studied by applying both the hybrid ESF/FPV approach and the presumed PDF method, respectively. The first configuration consists of the piloted CH4/Air diffusion Sandia flame-D [27] which is studied in order to validate the hybrid ESF/FPV approach. The second features an Oxyfuel configuration with cases series A and B [28] among them, three cases (A1, A3, B3) are numerically considered to point out the effect of O 2 /CO 2 dilution in the oxidizer stream, of H 2 % enrichment in the fuel nozzle and of Reynolds numbers.

Validation case: Sandia Flame D
The piloted coaxial methane-air jet flame (Sandia Flame D) [27] is characterized by 3 inlet streams and Re = 22,400. The main jet consists of fuel with a mixture of 25% methane and 75% air by volume. Its inside diameter is d fuel = 0.0072 m and the bulk velocity equal to 49.6 m/s. The pilot jet is composed of lean (phi = 0.77) mixture of C 2 H 2 , H 2 , air, CO 2 and N 2 with the same nominal enthalpy and equilibrium composition as methane/air at this equivalence ratio. The experimental configuration is represented in Figure 2 where also a 2D computational domain with 32,000 control volumes used is shown. Thereby, the symmetry property of the configuration is exploited in order to save computational costs. The reader is referred to [27] for more details about the experimental set up. Table 1 summarizes the operating conditions in accordance with experiments.
Note that this number of control volumes has been found sufficient to achieve grid-independent solutions. The numerical outflow conditions are imposed at x = 70 d. The outlet plane is set as waveTransmissive condition with P = 101.325 kPa, while all other variables have a zero gradient boundary condition. The time step ∆t used is equal to 1.1 × 10 −6 to ensure CFL-number below one. The convergence of the iterative procedure is assumed if all normalized residuals are smaller than 10 −6 . The calculations are achieved on 16 to 64 processors depending on the number of stochastic fields used (1 to 128SFi).

Oxy-Fuel Jet Flame Series
The oxy-fuel jet flame series were investigated experimentally at the Sandia National Laboratories [28]. They consist of two main flame sets, oxy-flames A (1-3) and B (1-3), where the differences between the two series are the Reynolds number and the different compositions of both fuel and oxidizer jets. According to [28], the burner consists of a fuel main jet with inside diameter of 5 mm and 0.5 mm wall thickness. The fuel nozzle is surrounded by a laminar coflow of diameter equal to 96.5 mm. Experimentally the H 2 content in the fuel jet helps the flame to remain attached to the nozzle. The fuel jet has its tip 40 mm above the coflow so that the mixed flows are considered fully developed once they reach the tip of the nozzle. Regarding the inlet compositions of the Oxyfuel flame, Table 2 summarizes the different compositions of the fuel for flames A1, A3 and B3 under investigation in the present paper. For the oxidizer stream, a molar percentage of 32% of O 2 and 68% of CO 2 -diluted are used rather than N 2 . As shown in Figure 3 the two jets are surrounded by a third co-flow which is considered as wind tunnel from where fresh air is flowing in order to accompany the flow of interest and prevent early mixing with ambient air. The third flow is for purely experimental reasons as clearly described in reference [28]. According to [28,50], using the O 2 /CO 2 mixtures instead of air for fuel combustion, the flame temperature can be reduced and NOx emissions are expected to be much lower than in air-diluted conditions. In particular, the composition of the fuel and oxidizer streams generates density fields very different to those found in methane-air flames. The numerical outflow conditions are imposed at x = 80 d, where the outlet plane is set as waveTransmissive condition with P = 101.325 kPa, while all other variables have a zero gradient boundary condition. The temperature is initially uniformly distributed with 300 K. Taking benefit of the configuration symmetry, for the numerical calculations, a simple 2D computational grid with 28,000 cells, is designed. This was found sufficient for grid independent solutions. The time step ∆t used is equal to 2.1 × 10 −6 . Also in these cases, the convergence of the iterative procedure is assumed if all normalized residuals are smaller than 10 −6 .The number of CPU's applied to carry out the simulations varies from 16 to 64 depending on the number of SFi employed (varying between 1 and 128SFi).

Validation Case: Sandia Flame D
Since this configuration is used to validate the hybrid ESF/FPV approach, the first task was to generate the FPV tables. The FPV tables used are based on non-premixed flamelets and generated according to the procedure in Section 2.2. To examine the convergence of the number of stochastic fields, SFi, the simulations have been done for different numbers of stochastic fields (8, 16, 48 and 128). To assess the combustion/turbulent model, the obtained results are compared to the experimental data in Figure 4 where different radial profiles of mean value of mixture fraction f, the root mean square (rms) of mixture fraction, temperature T and velocity U, are reported at three different axial locations,  Even though the Sandia Flame D is considered as non-complex case, the satisfactory results in Figure 4 qualify the developed hybrid ESF/FPV approach for further investigations.

Oxyflame B3 and Combustion Modelling Comparison
As pointed out above, three flames A1, A3 and B3 from the oxyflame case series (see Table 2) are considered in order to highlight the effect of O 2 /CO 2 dilution and ofH 2 % content as well as the impact of the Reynolds number on the flow field and the flame properties using the novel hybrid ESF/FPV approach.
Some of the properties of the oxyflame case B3 are presented in Figures 5 and 6, where both show contour plots of different quantities distribution at mid-plane using the ESF/FPV approach. In Figure 5 both the mean temperature field distribution in (a) and the mean CO 2 mass fraction in (b) are smooth and continuous. In Figure 6, the mean mixture fraction distribution and the mean axial velocity profile are displayed. From the contour plots, the flame seems to be attached to the nozzle of the burner by reason of important amount of diluted H 2 in the fuel mixture as it is mentioned in the experimental paper [28]. But unfortunately clear experimental qualitative images from the experiments are not available to be compared with the modelled contour-plots in Figures 5 and 6.  The prediction results of the properties of the flame are then presented by comparing them with the experimental data and with the results obtained by means of the presumed ß-PDF approach. Especially for the hybrid ESF/PFV approach, different number of stochastic fields (SFi = 1, 16,48,128) are considered in order to assess the convergence of the ESF. In Figure 7, a comparison of the mean mass fraction of O 2 and H 2 from the ESF/FPV calculation using different stochastic fields numbers, with experimental data is reported at 3 different axial positions (z/d = 3, z/d = 5 and z/d = 10). For the H 2 species, an acceptable agreement is clearly observed between the data at different positions for all SFi numbers. However, for the mean O 2 mass fraction prediction, the ESF/FPV calculations using high numbers of SFi are matching the experimental data, whilst there is an under-prediction of O 2 mass fraction with 1SFi results at z/d = 5 and z/d = 10. The deviation of the 1SFi results with respect to measurements is also visible in Figure 8 for the prediction of mean mass fraction of CO and H 2 O going from z/d = 3 to z/d = 10 for both species. These results are expected since using 1SFi means applying simple laminar FPV chemistry without any sub-grid model which is the similar behaviour of a perfectly stirred reactor. Increasing the number of SFi to 16 SFi leads to an under-prediction of both CO and H 2 O mass fraction at z/d = 10 while the cases using 48 and 128 SFi seem to be in very good agreement with experimental data further from the nozzle. This means that, the use of very high number of stochastic fields in the order of 128 is reproducing closely similar results to the case with 48 stochastic fields. This makes clear that the 48 SFi emerge as acomprise number between better prediction and computational costs. To note is that 128 SFi necessitate 64 CPUs while the 48 SFi only requires 32 CPUs.
Regarding the calculations using different combustion sub-models, the two different approaches employed, the hybrid ESF/FPV and the presumed ß-PDF, deliver different prediction results for the CO and H 2 O mass fraction once compared to experimental data. While the hybrid ESF/FPV reproduce closely the reference experimental data of the species, the assumed β-PDF approach suffers from some limitations due to intrinsic assumptions made, like the consideration of statistical independence between single PDF, along with the modelling applied to the source term during its calculation. This prediction failure of the presumed PDF method is also confirmed by the evolution of the mean temperature profile as function of the mixture fraction at z/d = 3 close to the fuel nozzle in Figure 9. In reference [28], the stoichiometric mixture fraction is reported as 0.056 with maximum adiabatic temperature around 1750 K. The hybrid ESF/FPV calculations reproduce a value of 1700 K for both simulations with 48 and 128 SFi. In contrast, the ß-PDF results clearly under-estimate the temperature evolution with a maximum value of 1300K. Table 2 shows that, with the same Reynolds number equal to 15,000, the two jet flame cases A1 and A3 feature different percentage of H 2 % enrichment in the fuel side, with 55% and 37%, respectively. As reported in reference [28], the extinction level increases from flame A1 to A3 and its effect is reported together with the reduction of the mean temperature values around the stoichiometric mixture fraction as reproduced numerically at the axial positions z/d = 3 and z/d = 5 in Figure 10.

H 2 % Enrichment in Fuel Side
Although the increase in temperature with H 2 % enrichment is observed close to the nozzle from A1 to A3, this temperature difference disappears at positions far from the nozzle at z/d = 10 where the numerical prediction of the temperature profile of both flames leads to similar results. Unfortunately, experimental results of temperature at further positions, z/d ≥ 10, are not available as can be seen in Figure 10. Not only is the temperature's peak value changing while reducing the H 2 % enrichment in the fuel for flame A3 at z/d = 3 but also the shifting of the maximum adiabatic temperature from stoichiometry toward rich side is visible. Further, Figure 11 presents the radial profiles of mean mass fraction of both CO and H 2 O.
For the flame A1 with the highest CH 4 /H 2 ratio, both CO and H 2 O productions are higher than in flame A3 at different axial positions in accordance with experimental findings. This is of interest despite that the Lewis number effect has not been included at the present stage of the FPV combustion model assuming Le = 1.

CO 2 Dilution in Oxidizer Side
The CO 2 amount diluted within the oxidizer is constant in flames A and B series and is about 68% which is quite a high level. This results in the high level of production of CO and H 2 O as it can be seen in Figure 11. For flame cases A1 and B3 which include both the same H 2 % enrichment in fuel side and same CO 2 dilution, the CO mass fraction locally reaches at z/d = 5 an amount of 0.14 and increases at z/d = 10 to reach 0.16 which is not a regular value for cases with air-diluted flames. This confirms previous finding by Masri et al in reference [51] who presumed that the CO 2 diluted in oxidizer is not inert and CO high level formation is the result of the reaction of CO 2 with H to form CO species. Therefore, the CO production level in flames A1 and B3 is manifestly higher than in A3 at all positions. Only the results with the hybrid ESF/FPV could reproduce this trend except for H 2 . Furthermore, the temperature is slightly reduced once compared to the air/methane flame.

Reynolds Number Effect
Both jet flame cases A1 and B3 share the same CH 4 /H 2 ratio in the fuel side but are characterized by different Reynolds numbers. From Figure 10, one can observe that the maximum adiabatic temperature location toward the mixture fraction space remains nearly the same for both flames at 3 and 10 diameters above the nozzle. However in Figure 11, there is a clear augmentation in the production of mean mass fraction of both CO and H 2 O species at the axial position z/d = 3 for flame A1 with lower Re-number. The difference of CO and H 2 O formation is reduced far from the nozzle. It turns out that the mixing state near the nozzle with lower jet Reynolds number likely leads to higher CO formation level. It turns out, that only the hybrid ESF/FPV is able to reproduce satisfactorily this trend.

Conclusions
Numerical investigations of an oxy-fuel non-premixed jet flame, with high CO 2 diluted level in the oxidizer and different CH 4 /H 2 ratios, have been carried out using a novel hybrid Eulerian Stochastic Field (ESF)/Flamelet Progress Variable (FPV) combustion model within the RANS modelling framework. After a successful validation of the combustion model in the piloted CH4/Air jet flame, the combustion behaviour of an oxyfuel jet configuration featuring the flame series A1, A3 and B3 has been studied. These flames exhibit different CH 4 /H 2 and O 2 /CO 2 ratios in the fuel and oxidizer streams, respectively and are characterized by different Reynolds numbers. This study allowed for tracing the impact of these properties on the temperature profiles and the CO and H 2 O formation. The obtained results were compared to available experimental data and to achievements accomplished by using a presumed β-PDF combustion model. Following important conclusions can be drawn down:

1.
A good prediction of different experimental and flow field variables, is reported by using the novel ESF/FPV approach.

2.
Related to the convergence of the stochastic field number (SFi), it turned out that starting the calculations with 48 SFi emerged as the compromise between accurate prediction and computational costs.

3.
Comparing the two different PDF-based combustion models, it turned out that the hybrid ESF/FPV clearly showed superiority in better predicting the temperature, H 2 O mass fraction and specially CO mass fraction, unlike the presumed β-PDF model which under-estimates the maximum adiabatic temperature and over-predicts the CO formation level at different positions above the nozzle downstream.

4.
With lower H 2 % enrichment in fuel side, fixed CO 2 /O 2 ratio and constant Reynolds number, the maximum adiabatic temperature value decreases in a significant manner near the fuel nozzle and its location in the mixture fraction space is shifted toward the reach side of the fuel.

5.
With lower Reynolds number, constant CO 2 /O 2 and CH 4 /H 2 ratios, the CO formation is considerably intensified near the nozzle.
Despite encouraging results provided by the hybrid ESF/FPV approach, it is still important to consider the H 2 differential diffusion effect which may affect the prediction of other minor species. The results reported in the present paper may suffer from the limitations of the RANS turbulence model used. The task of coupling the hybrid ESF/FPV approach with Large Eddy Simulation technique for more accurate prediction of the turbulence, along with the turbulence-chemistry interaction, is left for future research work.