The Eulerian Stochastic Fields Method Applied to Large Eddy Simulations of a Piloted Flame with Inhomogeneous Inlet

Large Eddy Simulations of the Sydney mixed-mode flame with inhomogeneous inlet (FJ200-5GP-Lr75-57) are performed using the Eulerian Stochastic Fields (ESF) transported probability functions method to account for the sub-grid scale turbulence–chemistry interaction, to demonstrate the suitability of the ESF method for mixed-mode combustion. An analytically reduced 19-species methane mechanism is used for the description of the chemical reactions. Prior to the reactive case, simulation results of the non-reactive setup with cold and hot pilot stream are presented, which show differences in the jet breakup and radial species mass fluxes. The reactive case simulations are compared to experimental data and a recently conducted model free quasi-DNS (qDNS), showing very good agreement with the qDNS in terms of scatter data and radial mean values of temperature and species distribution, as well as mixture fraction conditional statistics. Further analysis is dedicated to sub-grid scale statistics, showing that mixture fraction and reaction progress variable are strongly correlated in this flame. The impact of the number of stochastic fields on the filtered temperature and species distribution is investigated; it reveals that the ESF method in conjunction with finite-rate chemistry is very insensitive to the number of employed fields to obtain highly accurate simulation results.


Introduction
Turbulent flames are usually classified into non-premixed and premixed combustion regimes. In non-premixed flames the fuel consumption is dominated by diffusion while premixed flames can propagate into a flammable mixture (Peters 2000). Based on this classification turbulent combustion models have been developed and extensively validated, but their application is limited to their specific regime. In practical combustion devices, such as internal combustion engines or gas turbines, the fuel-oxidizer mixture is rarely perfect and events like autoignition, recirculation of hot combustion products, and extinction/reignition can be important.
The effect of such inhomogeneities in the fresh gas composition on the stabilization of a piloted turbulent methane/air flame has recently been investigated experimentally by Barlow et al. (2015) as well as by Meares and Masri (2014). In these experiments, the main pipe of the burner consists of two concentric tubes that separately supply methane and air. The inner tube is retractable which allows to adjust the degree of inhomogeneity at the burner inflow plane. Hence, depending on the setup, the flame is either in a stratified-premixed or a diffusion-dominated mode. This configuration has gained attention in the combustion modeling community over the last years and is the subject of intensive investigation and model validation. Maio et al. (2017) conducted LES using a premixed flamelet approach to compare filtered tabulated chemistry with a thickened flame model. Good agreement with the experimental data could be obtained only close to the burner exit. Another approach is followed by Galindo et al. (2017), who employ the sparse-Lagrangian multiple mapping conditioning (MMC) method to simulate the homogeneous (FJ200-5GP-Lr300-59) and inhomogeneous (FJ200-5GP-Lr75-57) setups of this flame. While the homogeneous case is represented well they report difficulties to adequately reproduce the flame structure for the inhomogeneous case. Perry et al. (2017) combined the flamelet progress-variable (FPV) model with a second mixture fraction, showing some improvements over the single mixture fraction FPV model. Kleinheinz et al. (2017) proposed methods to combine premixed and non-premixed manifold representations to obtain a multi-regime model. Wang and Zhang (2017) performed a model assessment for different combustion regimes, including the Sandia partially premixed flame. Tian and Lindstedt (2019) performed RANS simulations with a particle tPDF method to investigate the statistical distribution of mixture fraction and reaction progress variable in the Sydney flames. Hansinger et al. (2017) investigated the effect of premixed and non-premixed manifold representations for the FPV model, both for the homogeneous and inhomogeneous case. Recently, Perry and Mueller (2019) used the inhomogeneous case to study the effect of different subfilter PDF models for the FPV model with multiple mixture fractions. However, as pointed out by Wu and Ihme (2016), a particular difficulty in the simulations of multi-mode flames using tabulated chemistry is the adequate representation of the thermochemical state space in terms of a reduced set of variables. Despite the continuous efforts to model partially premixed flames (Fiorina et al. 2015;Hu and Kurose 2019) there is further research required to understand the complex combustion process in these flames.
Transported PDF models on the other hand may be a suitable tool towards improved simulations of partially premixed flames as they a-priori do not presume any kind of combustion regime (Pope 1981). They allow to employ finite-rate chemistry in the framework of LES, the joint sub-grid probability density function of the reactive scalars is constructed via stochastic Monte-Carlo methods. Initially derived in the framework of Reynolds averaged Navier-Stokes (RANS) simulations and solved by coupling Lagrangian particles with the finite-volume flow solver, as described in the seminal work of Pope (1985), the PDF method has been reformulated independently by Valiño (1998) and also by Sabel'nikov and Soulard (2005) in an Eulerian framework, known as the Eulerian Stochastic Fields (ESF) method. The idea is that the joint sub-grid PDF is constructed from an ensemble of N s statistically independent fields which evolve according to a stochastic partial differential equation (SPDE) which itself is a stochastic equivalent to the transported PDF equation. There are two main advantages with this method: The first is its Eulerian formulation which allows a straightforward implementation into a finite volume CFD code. Second, the fields are smooth and differentiable in space (though not in time, given their stochastic nature) and statistical moments are easy to compute without spatially varying sampling errors. In the framework of LES it has been successfully applied in the simulation of reacting flows, as demonstrated by Mustata et al. (2006), Jones and Prasad (2010), Avdic et al. (2017), Hodzic et al. (2017) and Fredrich et al. (2019); an extensive review on the applications of PDF methods including a chapter on the ESF method is given by Haworth (2010).
So far not much research has been done on the application of ESF to mixed-mode combustion. In the present work we therefore perform LES of Barlow's inhomogeneous test case (FJ200-5GP-Lr75-57) using the ESF method. This work has close connections to a recent publication by Zirwes et al. (2020) who performed a highly resolved model free simulation of the same test case. They have named this type a quasi direct numerical simulation (qDNS) as it uses overall second order discretization schemes and achieves DNS resolution only in the upstream region of the flame. While the focus of Zirwes et al. is to gain a deeper understanding of the mixed-mode flame dynamics and provide a detailed data set of the Sydney flame our intention is to explore the suitability of the ESF method for mixed-mode combustion, its mesh dependency and sensitivity to the number of stochastic fields. Additionally, sub-filter statistics, extracted from the individual fields, are presented. To validate our results against experimental and qDNS results we use the same analytically reduced 19-species reaction mechanism by Lu and Law (2008) for methane/air combustion, employ topologically similar but coarser meshes, use the same inlet boundary conditions and apply both similar pressure based solvers implemented into OpenFOAM. Apart from the reactive flow simulations we also present simulation results of the non-reactive case with a cold and hot pilot stream, in order to gain further insights into the mixing process between fuel and oxidizer.
The paper is structured as follows: In Sect. 2 we present the mathematical background on LES and the ESF combustion model. The test case and the numerical setup are presented in Sect. 3. An analysis and discussion of the simulation results is given in Sect. 4. Finally, the paper is closed with a conclusion in Sect. 5.

Large Eddy Simulation
In the LES approach only the large scale contributions of velocity and scalar fluctuations are directly computed, whereas the effects of the unresolved sub-grid (SGS) motion has to be modeled. The separation of resolved and unresolved scales is mathematically achieved through the application of a low-pass filter. An arbitrary quantity f = f ( ) is then filtered through the spatial convolution with the filter operator G on the domain as: denotes the filter width, which itself can be a function of the spatial location. In LES, filtering is not done explicitly but is implicitly achieved through the choice of the numerical grid, which does not resolve all scales of motion. Hence, the computed quantities are already the filtered values and resemble the spatial mean on the given control volume, i.e. the computational cell. The effect of sub-grid motions needs to be modeled. In reactive flows, the often high variations of density can be treated through the use of a Favre-filter, defined as f = f ∕ . Applying both filter operations, the continuity, momentum, scalar, and enthalpy equations of motion read: Continuity Momentum with

Sensible enthalpy
In this set of equations ũ i denotes the filtered velocity in i-direction, is the cell averaged fluid density, is the dynamic viscosity and ij denotes the Kronecker delta. Ỹ k is the mass fraction for species k , h s the sensible enthalpy, ̇k is the filtered reaction source term for species k and ̇T is the filtered heat release due to chemical reaction (Poinsot and Veynante 2005) and is evaluated as ̇T = − ∑̇k h 0 k , with h 0 k being the standard enthalpy of formation for species k. Within this study Soret and Dufour effects, as well as differential diffusion are assumed to be negligible, thus the individual diffusion coefficients D k are approximated as the mixture averaged diffusion coefficient D. Since we assume unity Lewis number ( Le = 1 ) the laminar Prandtl and Schmidt numbers are equal ( Pr = Sc ). With Pr = ∕ C p the diffusion coefficient can be retrieved as D = ∕Sc. (1) In Eqs.
(3)-(6) unclosed terms ( ũ i u j , ũ j Y k and ũ j h s ) appear due to the filtering operation. These terms cannot be resolved on the grid anymore and require modeling. This sub-grid activity mainly accounts for an additional diffusive flux, i.e. it enhances scalar mixing and dissipates turbulence kinetic energy. In the Boussinesq approximation the sub-grid momentum flux sgs ij can be expressed with the turbulent eddy viscosity t and S ij as the resolved symmetric strain rate tensor: Introducing the turbulent Prandtl and Schmidt numbers ( Pr t = Sc t = 0.7 ) the unresolved SGS species and enthalpy fluxes can be modeled using the gradient flux hypothesis: The turbulent viscosity itself has to be modeled based on the underlying turbulence model. In this work we use the WALE sub-grid turbulence model (Nicoud and Ducros 1999).

Filtered Density Function
Due to the high nonlinearity of the reaction source term in Eq. 5, sub-rid fluctuations in the species composition cannot be neglected when evaluating the filtered reaction source term ̇k . These sub-grid fluctuations can be described using PDFs. The theoretical background is described hereafter, following Pope (1985) and Jaberi et al. (1999). The fine-grained probability density function P of any scalar quantity , i.e. species or enthalpy, at any time t and any point in space can be described by a Dirac delta function where is the 'phase space' for the scalar quantity , i.e. it describes the accessible sample space. The joint PDF P joint of all N scalar quantities involved is obtained as the product of the fine-grained probabilities P for each scalar quantity : with as the composition space of the N scalar quantities involved. Applying the spatial filter operator (Eq. 1) and Favre-filtering on Eq. 12 one obtains the joint sub-grid PDF P sgs for N scalar quantities as: (11) P ( ; , t) = ( − ( , t)) , 1 3 P sgs describes the probability of ≤ ≤ + d within the LES filter size at a given point in space and time. It is important to note that P sgs is an instantaneous quantity, which describes only the SGS statistics in terms of a PDF for a given cell volume at a given time. It does not contain any information on two point correlation or temporal evolution on the scalar composition. However, an exact evolution equation for this one-point, one-time PDF can be derived from the appropriate conservation equations. Following Gao and O'Brien (1993), Jaberi et al. (1999) and Haworth (2010) it reads: Equation (14) describes the temporal evolution (Term I) of P sgs according to convective transport on the resolved scales (Term II), chemical reaction (Term III), scalar conditioned velocity fluctuations on the sub-grid scale (Term IV), and SGS micro-mixing due to molecular diffusion (Term V). All terms on the left hand side (I-III) appear in closed form, this includes also the chemical source term, which does not require any further treatment as the filtered value is obtained through the convolution with the PDF. However, the last two terms (IV, V) do require modeling strategies as they are unclosed. Term IV can be interpreted as the sub-grid flux of the PDF due to SGS velocity fluctuations and is commonly approximated with the gradient flux hypothesis (Pope 1994;Mustata et al. 2006). Term V resembles the effect of molecular diffusion on the sub-grid scale. It counteracts scalar SGS variances, i.e. dampens the fluctuations caused by (IV). It has no impact on the filtered first moment of the PDF, and in the absence of chemical reactions it only changes the shape of the PDF, which eventually should relax to a Dirac delta function when perfect mixing on SGS is reached (in the limit of infinite time or an infinitely small filter width ). There are several models available to model this term, an overview is given, e.g. in Subramaniam (1998), Mitarai et al. (2005) and Meyer (2010). In this work we will apply the linear mean square estimate (LMSE) formulated by Villermaux and Devillon (1972), it is also known as interaction by exchange with the mean (IEM) model, as it is described by Dopazo (1979). Although there are more sophisticated models available, Jaberi et al. (1999) showed that the LMSE achieves good results in LES applications and it is relatively simple to implement. With the given modeling assumptions Eq. (14) can be recast into: � P sgs ( )

Eulerian Stochastic Fields Method
The PDF equation (15) is usually solved via stochastic methods. As an alternative to standard Lagrangian Monte-Carlo Particle methods (e.g., Pope 1981;Raman et al. 2005), Valiño (1998) proposed the Eulerian Stochastic Fields method (ESF). The joint PDF P sgs is here constructed from N s continuous scalar fields. Each of these fields contains a composition of the N scalars n ( , t) for 1 ≤ n ≤ N s and 1 ≤ ≤ N . Valiño et al. (2016) modified the ESF method to be consistent at low Reynolds numbers. This formulation is used here. Using the Itô formalism the stochastic evolution equation of the Eulerian fields can be written : The first term on the RHS represents the stochastic Wiener term which models the effect of turbulent subgrid diffusion. As we follow Valiño et al. (2016) only the subgrid diffusivity in the stochastic term is used, while Valiño (1998) originally included also the laminar dif- In the newer formulation, laminar flows with scalar gradients do not exhibit unphysical fluctuations from the stochastic term. Except for this detail, our implementation follows Jones and Prasad (2010) with the Wiener term dW n i = n i √ dt evaluated with n i being a dichotomic random value of {−1; + 1} which has to satisfy ∑ N s n=1 n i = 0 . The micro mixing time scale sgs , which already appeared in Eq. (14), is modeled as −1 Valiño et al. (2016). It should be noted that each single realization of a stochastic field n is not necessarily a physical realisation of the particular field (Valiño 1998;Valiño et al. 2016;Prasad et al. 2011;Jones and Prasad 2010) in contrast to the Favre filtered mean value evaluated as ̃ = 1 Besides the Itô formalism used here, Eq. (16) may also be written by means of a Stratonovich formalism, leading to a different interpretation and numerical implementation. This approach is used in Sabel'nikov and Soulard (2005) and Sabel'nikov (2006b, 2006a). A general overview of stochastic methods and the differences between Itô and Stratonovich integration can be found in Gardiner and Gardiner (2009).
Recently Wang et al. (2018) reported on the mathematical inconsistency of the ESF method and provided correction terms for the RANS context for the case of a single scalar field. They (Wang et al. 2018) argue that the second and any higher moments of scalars obtained from the ESF method are not consistent with the actual PDF equation (15). According to them, the original ESF formulation (Valiño 1998 spurious production of scalar variance, whereas the modified formulation (Valiño et al. 2016) (ESF) introduces a spurious dissipation effect to the scalar variance. In the present work, no correction terms have yet been considered due to various reasons. First, Wang et al. (2018) clearly state that the first moments obtained with the ESF method are correct, which is eventually the quantity we are interested in. Second, Wang et al. (2018) based their analysis on RANS where the second moments are more important and of higher magnitude as the cells are bigger. Third, they (Wang et al. 2018) state that the inconsistency does especially matter at small Reynolds numbers where turbulent diffusivity D t is considerably small, e.g. Wang et al. (2018) used test cases at Re = 20 . Since we investigate a highly turbulent test case we consider the error to be rather small. Moreover, the correction terms have been derived for the single scalar case only. A straight forward extension to the multi-species/scalar case is not obvious. However, we directly compare the SGS fluctuations of the ESF-O with the ESF method in Sect. 4.2.5 and compare radial averaged and mixture fraction conditional quantities in the Appendix.

Chemical Reaction Mechanism
A detailed chemical reaction mechanism of the methane/air combustion (GRI-3.0) 1 comprises 53 species and 325 single reaction steps. Due to the large number of species and the stiffness of this mechanism, an analytically reduced methane mechanism developed by Lu and Law (2008) is employed. The mechanism has been derived from GRI-3.0 using directed relation graph (DRG) (Lu and Law 2005) and computational singular perturbations (CSP) (Lu et al. 2001) methods. DRG resulted in a 30-species mechanism, with 11 quasi-steady-state (QSS) species identified by CSP, resulting in 19 species ( H 2 , H, O, O 2 , OH, H 2 O, HO 2 , H 2 O 2 , CH 3 , CH 4 , CO, CO 2 , CH 2 O, CH 3 OH, C 2 H 2 , C 2 H 4 , C 2 H 6 , CH 2 CO, N 2 ) and 184 reactions. Source term calculations involve additional algebraic equations to account for the QSS species, which are solved analytically. A detailed validation of this mechanism against the GRI-3.0 can be found in Zirwes et al. (2020) and shows excellent agreement for species and temperature distribution in the case of a freely propagating flame.

Experimental Setup
The main pipe of the burner consists of an outer tube with diameter D = 7.5 mm providing air (oxidizer) and a retractable inner tube with an inner diameter of d = 4 mm providing the methane (fuel). The burner is surrounded by a concentric pilot which is 18 mm in diameter. The axial position of the inner tube defines the degree of homogeneity between fuel and oxidizer at the jet outlet. If the tube is completely recessed ( L r = 300 mm), oxidizer and fuel are perfectly mixed. If the inner tube ends at the pilot plane ( L r = 0 mm), fuel and oxidizer are completely separated. From the experimental series by Barlow et al. (2015) we simulate the test case with a recess length of L r = 75 mm namely case FJ200-5GP-Lr75-57. Here the mixture leaves the burner in an inhomogeneously mixed state before it gets ignited by the surrounding hot pilot stream consisting of burnt products at stoichiometric conditions ( Z st = 0.055 at T = 2226 K). The bulk flow velocity of the fuel/air mixture is at U bulk = 57 m/s, corresponding to a Reynolds number ( Re = UD∕ ) of Re = 26, 800 , the pilot and surrounding co-flow streams have mean velocities of U P = 26.6 m/s, and U co = 15 m/s, respectively.

Numerical Setup and Mesh
The LES code solves the filtered compressible Navier-Stokes equations with a pressure based Low Mach algorithm and is implemented in OpenFOAM 4.1. 2 Temporal integration is performed by a second order implicit scheme. Second order accuracy is given for the momentum equations, due to the stochastic term the scalar transport equations are only first order accurate. Spatial discretization of the convective momentum terms is done by a second order central differences scheme (CDS) with filtering for high-frequency modes. The convective terms in scalar transport equations are discretized by a second order CDS with a van Leer limiter (1974). The time step is limited to fulfill the criterion CFL max < 0.3 . The turbulent Schmidt and Prandtl numbers, appearing in Eqs. (9) and (10), are set to Pr t = Sc t = 0.7 . The ESF method is implemented by a fractional step method, meaning the stochastic term ( dW n i ) is evaluated in a first step and fixed for each stochastic field prior to the solution of the species transport equation (Eq. 16). This procedure is equivalent to the first order Euler-Maruyama approximation of an Itô process (Garmory 2008).
Block structured meshes of different resolution have been used. The base mesh consists of approximately 4.6 million cells, the coarse version comprises 1 million cells. Compared to the mesh with 150 million grid points employed by Zirwes et al. (2020) for the qDNS our base mesh is coarsened by a factor of 33, or 3.17 in each spatial direction. In addition a highly resolved fine mesh with circa 60 million cells has been used for the nonreactive simulations only to validate our flow solver. All meshes are topologically identical and have a thorough mesh refinement in regions of shear layers and where the flame front is expected. Figure 1 shows exemplary the filter width over the radius r at the downstream position x∕D = 5 for the three different meshes. The Kolmogorov length l k , as it is  Zirwes et al. (2020), is presented as black dashed line for comparison. In the region of the hot pilot even the coarse mesh is able to resolve the turbulent length scales. The fine mesh resolves all length scales, except the eddies in the cold fuel jet.
The simulation domain starts at an axial position of x = − D in the fuel pipe and extends 66 D in stream wise direction; in radial direction it extends to a radius of 6 D. The inert mixing process in the pipe was not part of the simulation domain and the same boundary conditions as in Zirwes et al. (2020) have been used. There, a wall resolved qDNS on 150 million cells has been conducted to simulate the inert mixing process between methane and air in the fuel pipe. The velocity and species profiles have been sampled at a rate of 1 × 10 −7 s at a distance of one burner diameter D before the pipe outlet. This data is filtered to the LES grid and used as transient inflow boundary condition of the central pipe.
For the pilot and the surrounding co-flow velocity, block profiles with U P = 26.6 m/s and U co = 15 m/s have been applied. The composition of the pilot is assumed as burnt methane/ air at stoichiometric conditions with a temperature of T P = 2226 K, the pure air co-flow is set to T co = 300 K.

Results and Discussion
In this section we present at first a numerical study of the non-reactive FJ200-5GP-Lr75-57 configuration on different meshes; two setups, one with a hot, the other one with a cold pilot stream, are considered. Table 1 gives an overview.
The second part compares the simulation results of the reactive case FJ200-5GP-Lr75-57 with the experimental data and the qDNS results of Zirwes et al. (2020). It focuses on the impact and importance of the number of stochastic fields N s and provides insights into the sub-grid statistics and joint correlations of mixture fraction Z and a reaction progress variable c. The mixture fraction used for the following analysis follows Bilger's definition (Bilger et al. 1990): where Y s indicates the elemental mass fraction of element s and W k is the atomic weight of species k. Subscripts 1 and 2 denote the fuel and air streams, respectively. The reaction progress variable, which is considered in Sect. 4.2.5 is defined as:

Non-reactive Cases
The investigated test case has been one of the target flames of the 14th international workshop on turbulent non-premixed flames (TNF). 3 One of the workshop's findings was that this configuration challenges not only the combustion models, but also the turbulence modeling in terms of an accurate reproduction of the flow velocities and the turbulent diffusion of mixture fraction. To further investigate this we performed simulations of two nonreactive configurations. The first configuration uses a cold pilot stream at T P = 300 K. To keep the pilot's Reynolds number identical to the hot pilot case the inflow velocity is set to U P cold = 3.26 m/s. Experimental data are available in terms of velocity measurements. 4 The second non-reactive case uses the identical setup as the reactive case, i.e. a hot pilot and the same inflow velocities as described in Sect. 3, but the chemical source terms have been deliberately switched off. Although this may not seem physically correct, it allows to study the turbulent diffusion of mixture fraction under the influence of the hot pilot stream. Both of these non-reactive cases have been simulated on the fine, base, and coarse mesh in order to first have a highly resolved reference solution (fine) for further model validation, and second to demonstrate mesh convergence of the base and coarse mesh, which are used in the combustion simulations. Figure 2 compares the instantaneous (a) and time averaged (b) axial velocity U x of these two cases close to the jet inlet between −1 < x∕D < 10 on the fine mesh. The upper half shows the case with the hot pilot, the lower half the cold pilot, respectively. Figure 3 shows the same situation for the instantaneous (a) and time averaged (b) mixture fraction. The comparison of the instantaneous velocities shows that the velocity core of the cold case breaks up further upstream than in the case of a hot pilot. This is even more evident in the comparison of the time averaged velocity fields where the black lines depict the iso-contour of ⟨U x ⟩ = 65 m/s ( ⟨⋅⟩ denots temporal averaging). While the iso-contour of the mean velocity at 65 m/s of the hot case reaches up to x∕D ≈ 9 the iso-contour has its tip already at x∕D ≈ 5 in the cold pilot case. We attribute this to the adjustment of the cold pilot inlet velocity, which is with U P cold = 3.26 m/s significantly lower than the hot pilot velocity at U P = 26.6 m/s. This leads to higher shear velocities between jet and pilot in the cold case and higher momentum transfer in radial direction, eventually leading to a shear layer breakup and turbulent dissipation of the jet's kinetic energy. The effect of the radial momentum transport can also be seen in the mixture fraction distribution. In the comparison of the time averaged mixture fraction (Fig. 3b) the iso-contour of ⟨Z⟩ = 0.5 is depicted as a black line. In the hot case the tip of the iso-contour reaches up to x∕D ≈ 6 and surpasses the cold case by about 1 D. Figure 4 shows the time and circumferentially averaged mean and RMS values of axial velocity U x , and mixture fraction Z of the non-reactive cold pilot case at different axial positions, obtained by averaging the flow simulation over 10 ms (approx. 15 flow through times, based on the experimentally investigated range of 0 ≤ x∕D ≤ 30 ). The black line refers to the fine, red to the base, and blue to the coarse grid, solid lines show the mean, dashed lines the RMS values. Experimental data is only available for the velocities at x/D = 5, 10, and 20. The velocities from the fine and base mesh are in good agreement with the experiments and both of them are almost indistinguishable. On the contrary, the velocities from the coarse simulation show a tendency to underpredict the experimental mean in the center region of the jet by approximately 10%. The mixture fraction results in the second row show a high discrepancy between base, coarse and fine simulation, especially close to the jet exit at x/D = 5. While the base results slightly overpredict those of the fine case, the coarse one underpredicts these by about 20% close to the jet exit.  Figure 5 shows the same results on all three meshes for the non-reactive case with hot pilot, including the temperature which can be considered as passive scalar, as there are no reactions considered here. While there are significant differences in the velocity and mixture fraction fields between base and coarse simulation in the cold case, there are only minor differences between all three simulations in the hot pilot case. The simulated velocities are identical and the mixture fraction mean values show only differences in the core region of the jet up to x/D = 10, with the base solution being closer to the fine result. The same findings are true for the temperature where the results are very similar and only show small differences in the peaks of the mean values. Generally there is good agreement between all three simulations. From the previous results of the non-reactive cases with cold and hot pilot two major conclusions can be drawn: 1. The hot pilot confines the fuel stream and dampens turbulent fluctuations in radial direction, weakening a turbulent mixing process between fuel and pilot stream, which is a result of the low densities and higher viscosities (compared to the jet) in the hot pilot stream, while the cold pilot with lower flow velocity triggers a shear layer breakup between jet and pilot stream and thus fosters a turbulent mixing process between these two streams.

Although identical numerical setups and boundary conditions have been used for both
cases the mesh resolution significantly affects the velocity and mixture fraction results of the cold case, while for the hot case even the results from the coarse mesh are in generally good agreement with the reference solution of the fine mesh. There is a slight tendency on the coarser meshes to predict higher mixture fraction peak values on the center line of the jet. However, as there is no reaction taking place in this region it is expected that the simulations on the coarse mesh yield adequate results in the combustion simulations.

Reactive Cases
This

Time Averaged Flow Quantities
Time averaged mean and RMS values have been obtained by averaging over 10 ms which corresponds to approximately 15 flow-through times within the experimentally investigated range between 0 ≤ x∕D ≤ 30 . Figure  The temperature results (1st row) from the coarse and base simulation are almost identical. Near the pipe exit plane ( x∕D ≤ 10 ) they are also in very good agreement with experimental data and the qDNS results. Further downstream temperatures on the outer flame region ( r ≥ 6 mm) are overpredicted, when compared to the experiments, however the mean values are still in good agreement with the high-fidelity qDNS data. The qDNS RMS values are slightly higher at the outer flame region, indicating higher turbulent fluctuations of the temperature field in the shear layer between flame and co-flow.
The simulated mixture fractions (2nd row) agree very well close to the jet inlet ( x∕D = 5 ), although all three simulations overpredict the experimental mixture fraction in the core region of the jet. Further downstream, the qDNS predicts well the experiments, also in the center of the jet, while coarse and base results continue to show the tendency of predicting too high mixture fractions in the jet center but do not show significant differences among the two of them. Here it is speculated that the chosen turbulence modeling approach in LES, where the turbulent viscosity is computed from the eddy viscosity model, underestimates the production of turbulent diffusivity which results in scalar fluxes that are too small on the sub-grid scale. However, the accurate prediction of Z is one of the  difficulties of the investigated test case, as it has been reported and discussed on the TNF-14 workshop 5 among several other groups. The next rows show the results for selected major and minor species. For the major species ( O 2 , H 2 O, CO 2 ) they follow a similar tendency, coarse and base results generally agree very well and are also in good agreement with experimental and qDNS results in the first section of the domain, while they differ further downstream on the outer flame region. However, no significant mesh dependency is evident. For the minor species (CO, H 2 ) the situation is different. Apparently the resolution of the mesh has an influcence on the accuracy of the results. The best agreement with experimental data is seen for the qDNS, followed by the results on the base and coarse mesh.
The last row depicts and compares the mean heat release from the simulations (without experimental data as this quantity cannot be measured). The heat release Q output from OpenFOAM is usually in the unit watt per volume (W/m 3 ). The cell volumes of all three meshes differ significantly and so does the heat release in (W/m 3 ). For a better comparison Q has been normalized with the respective cell volume and is presented in watt (W) only.
While qDNS and base results agree very well and their peaks gradually increase in downstream direction, the heat release from the coarse mesh has its peak at x∕D = 5 , moreover, it is approximately more than eight times higher. It then decreases at x∕D = 10 and increases again until it is finally in the same order of magnitude as qDNS and base results at x∕D = 20 . It is assumed that the base mesh is already fine enough to resolve reactions zones, which explains why it coincides so well with the qDNS, while the coarse mesh isdespite using the ESF combustion model-too coarse to adequately predict the heat release in the front region where the flame is in a predominantly premixed combustion mode. This assumption is supported by a recent study from Picciani et al. (2018) who showed that the ESF method significantly overpredicts the fuel consumption rate in premixed combustion when the grid spacing is not in the order of magnitude of the laminar flame thickness. Fuel consumption rate and heat release rate are both a result of the species reaction rates and tend to be overpredicted on the coarse mesh. This could explain the Q peak at x∕D = 5 , which decreases and is in better agreement further downstream, where the flame transitions to a non-premixed combustion mode.

Instantaneous Scatter Data
From the experiments instantaneous scatter data is available. Figure 8 compares the instantaneous temperature in mixture fraction space from the experiments (1st column), and the results obtained from the base (2nd column) and coarse (3rd column) mesh at the different axial locations. The fourth column presents the mean value of temperature conditional on mixture fraction ⟨T�Z⟩ , here are also the qDNS results included. In order to guarantee comparability among the plots the samples are taken at the same radial locations as in the experiments and always the same number of scatter points is displayed. Figure 9 shows scatter data for the CO mass fraction.
The experimental temperature data show the transition from the partially premixed/ auto-ignition state close to the inlet ( x∕D = 5 ) to a non-premixed diffusion flame further downstream. On both meshes this transition is reproduced very well. Visually both results are also in good qualitative agreement with the experiments, demonstrating the 1 3 suitability of the ESF method to cope with different combustion regimes. Minor differences can be seen at x∕D = 20 where the experiments show some extinction and re-ignition events between 0 ≤ Z ≤ 0.1 , which are not represented in the simulations. However, this is not attributed to the combustion model itself, it is rather speculated that these events are a result of the turbulent breakup of the flame in the rear part, which is not reproduced well in the LES. The scatter data in Fig. 9 indicate a similar behavior for the CO mass fraction. Although this species is usually difficult to model the scatter plots of both cases are in very good agreement with the experiments and qDNS. The experimental peak values of CO mass fraction around stoichiometry are slightly underpredicted.
Conditional mean values have been obtained from the scatter data by averaging the scatter data over the mixture fraction in bins of Z = 0.001 . The quantitative comparison of the conditional temperatures in the fourth column of Fig. 8 shows very good agreement of the two LES with the experiments and the qDNS. Further downstream the LES start to slightly differ from the experiments on the fuel rich side, but are still close to the qDNS data. Both LES coincide very well, no significant difference can be seen between them. The conditional values of CO (Fig. 9, 4th column) show very good agreement for both LES cases on the fuel lean side at all downstream locations. Around stoichiometry and up to Z = 0.1 the base mesh shows better agreement with the qDNS. Towards the fuel rich side both LES tend to overpredict CO mass fractions.

Probability of Local Extinction
While scatter plots provide a useful visual insight into the flame structure and the degree of local extinction at each measurement position, conditional statistics of temperature and reactive scalars might be more useful for a quantitative comparison between simulation and experiments and to assess the quality of the employed combustion model. Figure 10 shows the PDFs of temperature, as well as CO 2 , CO, and H 2 mass fractions conditional on mixture fraction within the interval 0.005 ≤ Z ≤ 0.006 around stoichiometry, i.e. where the highest temperatures are expected at different axial planes. All PDFs have been normalized, so that they integrate to one. Experimental data is presented in the form of (grey) histograms, while LES and qDNS data are shown in the form of line plots. As expected from previous analysis the temperature PDFs of LES and qDNS generally agree well. At x∕D = 10 the experiments show the highest probability for temperature at around 2000 K, while all three simulations have a shifted peak at 2100 K, i.e. predict higher temperatures with a higher probability. Further downstream, the probability of predicting temperatures at 2000 K is slightly overpredicted by both LES. A possible explanation for temperature deviances between experiments and simulations might be the fact that simulation models do not account for radiative heat losses, generally leading to higher temperatures in the simulations. CO 2 mass fractions are predicted in both LES with the same probabilities and coincide with the experiments close to the jet exit ( x∕D ≤ 10). Discrepancies can be observed further downstream where qDNS and experiment agree well and show a wide spread range of mass fraction probabilities between 0.1 ≤ CO 2 ≤ 0.14 . It is speculated that a lack of turbulent fluctuations in the LES leads to a higher and more confined PDF profile. Regarding the probabilities of CO and H 2 mass fractions all simulations are in good agreements with the measurements, except for H 2 at the plane x∕D = 5 . Here, the experiments indicate a normal distribution centered around 0.001, whereas the simulations all consistently exhibit an exponentially-shaped distribution with maxima at H 2 = 0 . This is possibly attributable to the use of different pilot flame compositions. While the simulations assume a fully burnt CH 4 /air mixture at an equivalence ratio of = 1 it is a five component gas mixture in the experiment. Besides CH 4 one of the components is H 2 which eventually contributes the large deviations in the H 2 PDFs.
Trends in the probability of local extinction can be compared further by calculating the burning index BI T based on temperature. BI T is computed from the individual temperatures T i of all samples N within the mixture fraction interval 0.05 < Z < 0.06 around stoichiometry: Here 300 K is the inlet temperature for both fuel and coflow inlet, the burnt condition has been chosen as 2100 K, which was the maximum temperature which has been measured in the experiments (Barlow et al. 2015). Figure 11 compares the burning indices of the experimental data with the results from base, coarse, and qDNS (Zirwes et al. 2020) simulation. All BI T reproduce well the trend of local extinction from the measurements.
At x∕D = 5 the burning index of the coarse result matches the experimental index with BI T = 1 , while the base simulation indicates slight local extinction. Further downstream, both LES and qDNS results show a very similar BI T and follow the trend of local extinction observed in the experiments, though, they all predict a slightly higher burning state. Compared with the index values obtained from the qDNS both LES results are in very good agreement.

Effect of Number of Stochastic Fields
In this section it is investigated if and how the number of stochastic fields N s affects the accuracy of the filtered results. Based on previous studies with the ESF method in turbulent combustion simulation (e.g. Mustata et al. 2006;Jones et al. 2007;Jones and Navarro-Martinez 2008;Jones and Prasad 2010) a number of N s = 8 stochastic fields has been established to be sufficient to describe the SGS scalar fluctuations. We therefore performed additional simulations with different numbers of N s = 1, 8, 16 and 64 stochastic fields and investigate its impact. All these simulations have been conducted on the coarse mesh only, due to limited computational resources and because the influence of the number of fields is expected to increase on coarser meshes. Figure 12 compares the radial mean and RMS results (only presented: T, Z, CO, CO 2 ) of the four simulations with different number of stochastic fields; qDNS and experimental data are not shown in this plot as the focus is on a direct comparison between the simulation results. However, as it can be seen there are only small differences in the radial mean and the RMS results. Figure 13 depicts the mean values of temperature and CO mass fraction conditional on mixture fraction. Again, independent of the number of stochastic fields all simulations exhibit very similar results. Minor differences can be seen for N s = 1 in the CO mass fraction at x∕D = 5 and 10 where the conditional values are slightly lower. However, based on the underlying results it cannot be confirmed that a high number of N s fields improves the overall accuracy significantly. As a matter of fact using one stochastic field (which is the approximation of ̇k ≈̇k ) already accurately represents the present configuration.  Fig. 13 Comparison of mixture fraction conditional mean values of temperature and CO mass fraction using different numbers of N s stochastic fields

SGS Statistics
In this section the SGS statistics, obtained from the individual fields, will be examined in more detail. Figure 14 shows the temporal evolution of the temperature at a fixed point ( x∕D = 5 , r = D ) in the shear layer and reaction zone between fuel jet and hot pilot for the first millisecond of simulation time (approx. 2200 time steps) for the coarse simulation with N s = 8 . The thick black solid line represents the filtered temperature T , the temperatures of the individual stochastic fields are depicted in grey, the red shaded regions illustrate the RMS around the mean. Two aspects are worth mentioning here. First, the individual fields need some time to evolve until they decorrelate properly and constitute a PDF with a standard deviation around the mean. In this particular simulation this is the case after 0.15 ms (approx. 200-300 time steps). Second, although at the given location the filter width is with x ≈ 0.6 mm still very fine the individual fields may decorrelate strongly. This can be seen after 0.45 ms where individual fields predict temperatures between almost 2000 K at the maximum and as low as 1200 K. However, the SGS temperature fluctuations are confined to a narrow region in the reaction zone between x∕D = 0 and 20. This is expected since the stochastic velocity term is  Figure 16 shows the time averaged temperature SGS fluctuations at different axial positions in a more quantitative fashion. In addition to the ESF solution with 8 and 64 fields also the SGS fluctuations from a simulation with the original ESF-O formulation (Valiño 1998) and 8 fields are presented for comparison. In the reactive shear layer ( r ≈ 4 mm) the results are comparable in terms of magnitude, although the fluctuations of the ESF-O simulation are the highest. However, the differences between ESF-O and ESF with N s = 64 are smaller as between the two ESF results with N s = 8 and N s = 64 fields. A more evident difference can be noted in the shear layer between pilot stream and coflow ( 8 ≤ r ≤ 12 mm) where the ESF-O simulation introduces high SGS activity, which is only slightly present in the modified ESF simulations. Since there is no reaction in the outer shear layer these fluctuations have no impact on the combustion process. A direct comparison of the mean and RMS results of the ESF and ESF-O simulation on the coarse grid and a brief discussion is given in the Appendix.
From the simulation with N s = 64 fields we further investigate the common assumption of statistical independence between the mixture fraction Z and a normalized reaction progress variable c n , which is based on the mass fraction of reaction products and normalized by the maximum of the reaction progress conditional on the mixture fraction c max|Z : The common assumption states that the joint sub-filter probability distribution P(Z, c n ) can be approximated by a convolution of the marginal probability distributions P(Z, c n ) ≈ P(Z)P(c n ) , assuming statistical independence between both variables (e.g. Pierce and Moin 2004). Statistical independence demands that both variables are not correlated, so they need to have a correlation coefficient of zero ( R Zc n = 0 ). The converse, that uncorrelated variables are automatically independent is also not true. In order to check if this assumption is valid in the given test case the correlation coefficient R Zc (Eq. 26) has been computed from the individual stochastic fields based on the covariance Z′′ c ′′ n (Eq. 25) of mixture fraction and progress variable and the individual variances Z ′′ 2 (Eq. 23) and c ′′  Figure 17 shows the scatter data of R Zc n over the mixture fraction at different axial positions. It can be seen that, independent of the axial position, the samples are positively correlated ( R Zc n ≈ 1 ) on the fuel rich and negatively correlated on the fuel lean side ( R Zc n ≈ −1 ), with a transition point around stoichiometry. It appears that close to the exit, where the flame is in a partially premixed state, the correlation between Z and c n is slightly suppressed; further downstream, where the flame evolves towards a non-premixed flame, the correlation (both negative and positive) is more evident. As mentioned before, for Z and c n to be uncorrelated it is necessary that |R Zc n | << 1 and this assumption clearly does not hold. The joint sub-filter statistics from the ESF method clearly indicates a strong correlation between Z and c n and that the assumption of statistical independence is implausible for the given flame. The findings are also in very good agreement with the results of Tian and Lindstedt (2019) who did a similar analysis of the same flame using RANS with a particle PDF method.  Fig. 17 Scatter data of correlation coefficient R Zc over mixutre fraction Z 1 3

Conclusion
The Sydney mixed-mode flame configuration FJ200-5GP-Lr75-57 has been successfully simulated using LES with the ESF transported PDF method to account for turbulence-chemistry interaction and a 19-species analytically reduced methane/air mechanism to describe the chemical reactions. In addition, highly resolved simulations of the non-reactive case with cold and hot pilot stream are presented to shed light on the spatial evolution of the mixing process between fuel and oxidizer and are used to ensure mesh convergence of the base and coarse meshes, which are used in the combustion LES. Moreover, the comparison shows that the hot pilot significantly dampens turbulent fluctuations in radial direction and thus suppresses the turbulent mixing process and thereby increasing the importance of laminar transport. The combustion LES performed on the base and coarse mesh with N s = 8 stochastic fields show good agreement with the experiments and very good agreement with the qDNS results of Zirwes et al. (2020) in terms of scatter data and the radial distribution of mean and RMS values of temperature and species mass fractions. The mesh resolution has only very little effect on the filtered results of temperature and major species. The mixture fraction is also in very good agreement with the qDNS in the regions where the chemical reactions take place. It is, however, highly overpredicted on the center line of the jet. Many researchers who numerically investigate the present test case have found similar results. Future research could address a comparative study with different LES turbulence models and a variation of the turbulent Prandtl and Schmidt numbers to close this gap.
Minor species (CO, H 2 ) which were simulated on the base mesh are in slightly better agreement with the qDNS results. A significant difference in the performance of the meshes can be seen in the direct comparison of the mean heat release Q , where the result from the base mesh coincides very well with the qDNS, whereas the result on the coarse mesh significantly overpredicts Q . The reason why the prediction accuracy for the radial distribution of CO, and H 2 mass fractions, and Q deteriorates with decreasing mesh resolution, while major species and temperature are in good agreement, is not fully clear. It might be attributed to the ESF method, which has difficulties to integrate underresolved scalar quantities, e.g. the heat release in the thin reaction layer at x∕D = 5 where the flame is in a predominantly premixed mode, on very coarse meshes.
However, analysis of the temperature and CO scatter plots, as well as statistical analysis of the probability of local extinction further prove the high quality of the results, since extinction and re-ignition events and the burning index BI T are again in very good agreement with qDNS and experimental data.
The impact of the number of stochastic fields on the resolved flow has been investigated on the coarse mesh and is proven to show only very little effect, indicating that even simulations with N s = 1 stochastic fields would yield good results in the present case, though this would mean deliberately neglecting turbulence-chemistry interaction by approximating the filtered reaction rates with the resolved ones ( ̇k ≈̇k).
Finally, insights into the sub-filter statistics are presented and the correlation coefficient R Zc between mixture fraction and progress variable are computed. It reveals a high correlation between the two variables indicating that the commonly used assumption in flamelet/progress variable methods of statistical independence between Z and c n P(Z, c n ) ≈ P(Z)P(c n ) is not valid in the studied flame. The results of this work demonstrate that LES in combination with a complex reaction mechanism and the stochastic fields method are a valuable tool for the simulation of complex combustion problems arising in mixed-mode combustion and partially premixed flames.
1 3 interesting to apply the corrections terms proposed by Wang et al. (2018) if their model can be generalized to the case of reacting multi-species chemistry.