Quasi-DNS Dataset of a Piloted Flame with Inhomogeneous Inlet Conditions

A quasi-DNS of the partially premixed turbulent Sydney flame in configuration FJ200-5GP-Lr75-57 has been conducted using detailed molecular diffusion for multi-component mixtures and complex reaction mechanisms. In order to study flame dynamics like regime transition in this flame for the development of new combustion models and to directly compare the quasi-DNS to different LES models, the simulation results are compiled into a data base. Because the simulation was performed with OpenFOAM, we demonstrate the quasi-DNS capabilities of OpenFOAM by performing canonical test cases. They attest that OpenFOAM’s cubic discretization has lower numerical diffusion compared to classical central difference schemes and can reach higher than second order convergence rate in some cases. The quasi-DNS of the Sydney flame is conducted with a self-developed reacting flow solver which is able to accurately compute molecular diffusion coefficients from kinetic gas theory and employs a fast implementation for detailed reaction mechanisms. The computational mesh is shown to be able to resolve the flow as well as the flame front sufficiently for the quasi-DNS. Comparisons with experimental data also show that the simulation can quantitatively reproduce measured time-mean and time-RMS statistics.

These partially premixed or mixed combustion mode flames exhibit characteristics of both premixed and nonpremixed flames. By deliberately creating inhomogeneous mixing conditions, the stability of flames can be improved and pollutant emissions decreased [1]. Classical combustion models often fail to correctly simulate this type of flame because they are derived for either perfectly premixed or nonpremixed flames. Mixture fraction based models, e.g. the flamelet model [2], are originally derived for nonpremixed flames, while other models such as level-set methods are only valid for premixed flames [3]. Although different modeling approaches have been developed and tested specifically for partially premixed flames [4][5][6], the complex nature of these inhomogeneous flames still requires further research and a deeper understanding is a prerequisite for designing more efficient combustion devices in the future [7,8].
An experimental setup for flames with inhomogeneous inlet conditions [9,10] is the Sydney/Sandia burner [1], which is actively being investigated. It consists of an inner fuel tube annularly surrounded by an air duct. The inner tube can be retracted to vary the degree of premixing of fuel and oxidizer for the developing flame. This setup is well suited for numerical simulations due to its simple geometry and availability of ample experimental data [1,7,11,12]. Therefore, many simulations of this setup have been performed in the last few years, especially using large eddy simulations (LES): Kleinheinz et al. [13] performed LES with a multi regime flamelet model to study the contribution of premixed and nonpremixed regimes to heat release rates and flame stability; Perry et al. [14,15] conducted simulations utilizing a combustion model with two mixture fractions to take the inhomogeneous inlet conditions into account; Ji et al. [16] applied a multi-environment probability density function model to the Sydney burner while Tian et al. [17] performed LES with a modified PDF transport approach to investigate the statistical distribution of mixture fraction and reaction progress variable.
In order to gain a deeper understanding of the flame dynamics in the Sydney flame and to develop improved combustion models for partially premixed flames in general, a model-free quasi direct numerical simulation (quasi-DNS) of that flame is performed in this work. Model free means that no subgrid combustion models or turbulence models are used and that all length and time scales are mostly resolved. Additionally, a complex reaction mechanism is used for the combustion of methane in air together with detailed molecular transport models for multi-component mixtures to account for preferential diffusion. The goal of this work is to provide highly resolved and detailed datasets for the mixed-mode Sydney flame configuration for studying partially premixed flames in general and to demonstrate the validity of the quasi-DNS by means of various validation cases.
The quasi-DNS of the Sydney flame has been performed in this work with the opensource CFD toolkit OpenFOAM [18], which solves the Navier-Stokes equations using the Finite Volume Method (FVM). OpenFOAM has been used by many groups to perform quasi-DNS in the past: Komen et al. [19,20] showed that DNS quality results can be achieved for channel flows and pipe flows if structured hexahedral meshes are used. Channel flows were investigated by Jin et al. [21] and compared to Lattice Boltzmann codes. Habchi et al. [22] performed DNS of electromagnetically forced turbulent flows. Addad et al. [23] investigated natural convection in cylindrical annuli and compared channel flow results to those from other DNS codes. Lecrivain et al. [24] tested the quasi-DNS capabilities of Open-FOAM for particle deposition and found good agreement with results from literature. Chu et al. [25][26][27] performed DNS of heated turbulent pipe flows. Zheng et al. [28,29] used OpenFOAM to perform DNS of turbulent pipe flows with non-Newtonian fluids and compared the results with a spectral element-Fourier DNS code, where they found good agreement for meshes with hexahedral cells. Bricteux [30] applied OpenFOAM to various canonical test cases including the Taylor-Green vortex. OpenFOAM has also been used to perform quasi-DNS in the context of combustion. Zhong et al. [31] computed turbulent premixed methane flames. Tufano et al. [32,33] performed DNS of coal combustion and Wang et al. [34,35] performed DNS of droplet combustion. Vo et al. [36] performed DNS of nanoparticle synthesis. Zhang et al. [37] performed resolved simulations of oscillating methane and hydrogen slot burner flames. The combustion DNS performed by Vo et al. [38] showed good agreement with another DNS code for nonpremixed combustion in a double shear layer. In this work, the quasi-DNS capabilities of OpenFOAM are demonstrated using different canonical test cases with a focus on the discretization schemes offered by OpenFOAM.
In contrast to numerous other CFD-codes based on finite differences (see Section 2) OpenFOAM does not offer discretization schemes of higher orders. Further, the computational mesh employed in this work achieves DNS quality resolution generally only in the upstream region near the nozzle exit where results from experiments are available (see Section 4.2). Therefore, the term quasi DNS is chosen for the model-free numerical simulations performed in this work.
The paper is structured as follows: Section 2 shows the general quasi-DNS capabilities of OpenFOAM by discussing its discretization schemes. Section 3 introduces the selfdeveloped solver for reacting flows. The numerical setup for the quasi-DNS of the Sydney flame is given in Section 4 together with validation of the mesh resolution. Section 5 compares the simulation results with two experimental data sets. Section 6 describes the contents of the database and gives a brief outlook of how the provided database will be used in future work.

Quasi-DNS Capabilities of OpenFOAM
Accurate numerical discretization schemes are a requirement for quasi-DNS. OpenFOAM is a finite volume method code that operates on unstructured numerical grids. Therefore, by default, it is not possible to use more than one layer of neighbor cells in the discretization schemes. OpenFOAM offers many different schemes for spatial discretization. The most accurate one uses a cubic interpolation, which is achieved by utilizing explicitly computed gradients at the cell centers. In this section, different canonical test cases are computed with standard OpenFOAM solvers and the accuracy and rate of convergence of different discretization schemes are tested. In this work, time is always discretized using an implicit second order backward Euler method, except for the cases which deliberately only use first order schemes.
The most accurate spatial discretization scheme available in OpenFOAM is the aforementioned cubic scheme. It is based on a third-order polynomial fit for interpolation which requires four coefficients. Assume the center of the current cell C is located at x = 0 and the center of a neighboring cell N at x = 1. The face f between the two cells is located at x = x f . The value of a property φ at the cell face can be interpolated by using a polynomial of third degree: The four coefficients of this polynomial are determined from the following four known values: where the gradients at the cell centers for Eqs. 4 and 5 are computed explicitly with a second order central difference scheme. Solving for the polynomial coefficients and introducing an interpolation The resulting polynomial is reordered so that the first two terms represent a standard linear interpolation which is treated implicitly by OpenFOAM. The other terms represent an explicitly computed correction term. This correction term increases the numerical accuracy compared to the standard linear interpolation leading to reduced numerical dissipation but makes this discretization scheme less stable and prone to oscillations due to the explicit correction. Therefore, low CFL numbers (CFL ≤ 0.2) and well resolved solution fields of all variables are required for the simulation.
In the appendix in Section A.1, four test cases are described to demonstrate Open-FOAM's quasi-DNS capabilities: 1D steady-state convection-diffusion equation, 2D and 3D Taylor-Green vortex and homogeneous isotropic turbulent decay. There, the numerical dissipation and convergence order of OpenFOAM's cubic discretization scheme are discussed in detail.

Reactive Flow Solver
The reactive flow solver [37,39,40] used for the simulation of the Sydney flame is implemented in the OpenFOAM framework. It is based on the standard rhoReac-tingBuoyantFoam solver and is coupled to the open-source library Cantera [41] for computation of thermo-physical and transport properties. Therefore, the solver extends OpenFOAM's capabilities by providing detailed molecular diffusion coefficients for momentum, energy and each chemical species. Additionally, it provides performance optimized routines for computing chemical reaction rates from detailed reaction mechanisms, which requires large computational resources [42]. The Navier-Stokes equations are solved in their fully compressible formulation. The equations are summarized below.

Governing equations
The governing equations are formulated as follows [43]: Conservation of total mass ∂ρ ∂t + ∇ · (ρu) = 0 where ρ is the density and u the fluid velocity; conservation of momentum ∂ (ρu) ∂t + ∇ · (ρuu) = −∇p + ∇ · τ + ρg (8) where p is the pressure and g the gravitational acceleration. The viscous stress tensor τ is computed for a Newtonian fluid using the Stokes assumption: where I is the identity tensor and μ the dynamic viscosity. The balance equation for the species masses is where Y k is the mass fraction of species k. The reaction rateω k is computed from detailed reaction mechanisms employing an operator splitting approach and Sundials [44] as ODE integrator. Different reaction types, like Arrhenius, Falloff or Third-Body reactions, are supported. N is the total number of species and u c the correction velocity [43]. Lastly, the transport of energy is formulated in terms of the total sensible enthalpy h s : where the species enthalpies h s,k and standard enthalpies of formation h • k are computed from JANAF polynomials. The energy diffusion is Thermal radiation is not included in the quasi-DNS of the Sydney flame. This allows easier comparison with LES since most simulations of this flame are performed without radiation models. Using a radiation model assuming an optically thin gas for CO 2 and H 2 O [45] showed that peak heat losses due to radiation reach approximately 0.1 % of the peak chemical heat release rates on average. The mass diffusion coefficients for each species are computed from rigorous kinetic gas theory and are based on the Chapman-Enskog solution of the Boltzmann equations [46]. The implementation uses routines from the open-source library Cantera [41]. The transport properties of the mixture are determined by mixing rules so that the diffusive flux of species k becomes j k = −ρD m,k ∇Y k (13) taking differential diffusion into account. Validation cases showing comparisons of the quasi-DNS solver implemented into Open-FOAM with reference solutions from Cantera for two canonical flame setups can be found in the appendix in Section A.2.

Numerical Setup for the Sydney Burner
After attesting the quasi-DNS capabilities of OpenFOAM in Section 2 and the suitability of the reactive flow solver for accurately reproducing the inner structure of flames in Section 3, this section presents the numerical setup for the quasi-DNS of the mixed-mode turbulent piloted Sydney/Sandia flame in configuration FJ200-5GP-Lr75-57.

Domain and numerical settings
In total, three consecutive simulations have been conducted where all solution variables at the outlet of the previous simulation are linearly interpolated to the inlet of the consecutive simulation: -a precursor LES to generate initial flow profiles (P) -a non-reactive quasi-DNS of the pipe flow and mixing (A) -a reactive quasi-DNS of the flame (B) This is illustrated schematically in Fig. 1. The domain of the precursor LES (P) on the left in Fig. 1 consists of two pipes: an inner pipe with a diameter d = 4 mm and a length of L/d = 6.25, providing methane with a bulk velocityū = 67 m/s; and an annular pipe with an inner diameter of 4.5 mm, outer diameter of D = 7.5 mm and a length of L/D = 3.3, providing air withū = 59.5 m/s. The time step is set to 1 · 10 −7 s and the mesh for the inner pipe consists of 2 million cells while the annular pipe consists of 1 million cells, both made of purely hexahedral cells. The LES subgrid model for the turbulence is the WALE model [47]. The simulation is run with OpenFOAM using second order discretization schemes for spatial and temporal derivatives. The boundaries at the inlet and outlet are periodic. In order to get a fully developed turbulent flow profile, this simulation was run for 120 flow-through times.
The second domain (A) for the nonreactive mixing of the methane and air flows has a total length of 10 D. The inner pipe ends after a length of 1 D from the inlet, so that the two flows can develop on the quasi-DNS mesh for a length larger than one integral length scale before methane and air flows start to mix. The inner pipe is retracted by 7.5 cm with respect to the combustion chamber in (B) in order to generate inhomogeneous mixing conditions. At the inlets of domain (A), velocity profiles are sampled at a rate of 1 · 10 −7 s from the LES solution in domain (P). The computational grid for domain (A) is block structured and consists of 150 million purely hexahedral cells. The mesh is refined radially with a smallest resolution of 5 μm at the walls. The third domain (B) has a total axial length of 68 D and a diameter of 6 D and its inlet plane is connected to the outlet plane of domain (A). The flame in domain (B) is simulated with a reactive quasi-DNS using the reactive flow solver described in Section 3 and a reaction mechanism by Lu et al. [48] further described below. The computational grid of domain (B) consists of 150 million purely hexahedral cells which are refined in the upstream direction, toward the walls of the pipe and the shear layer, with a smallest radial resolution of 10 μm. For a detailed discussion on the grid resolution see Section 4.2.
The simulation time step is set to 0.1 μs, which corresponds to a convective Courant number of about 0.18. Although the simulation is fully compressible in the sense that pressure and density are directly coupled, OpenFOAM solves the governing equations implicitly together with the PISO algorithm for pressure correction, so that the acoustic CFL number is not a stability criterion. A second order implicit Euler backward time discretization is employed together with the cubic discretization for all spatial derivatives (see also Section 2). The system of linearized equations from the governing equations are solved using preconditioned (bi-)conjugate gradient methods. The boundary conditions at the outlet are zero gradients for temperature T , velocity u and mass fractions Y k . For the pressure, a non-reflecting boundary condition [49] is used with an additional sponge layer near the outlet to filter out unphysical high frequency pressure waves. At the wall, the boundary conditions for T , u and pressure p are zero gradient and no-slip for u. The velocity at the inlet of the pilot (26.6 m/s) and the co-flow (15 m/s) is a block profile due to the small influence of the co-flow on the flame and the installations near the pilot exit in the experimental setup [1]. This is a commonly made assumption in numerical studies of this burner [16] and has been chosen here to ensure better comparability between different simulations. The composition of the co-flow is pure air at T = 300 K and burnt methane/air at equivalence ratio one for the pilot. The transient inlet condition for the inhomogeneous unburnt methane/air flow from the central pipe is sampled from the solution of the quasi-DNS in domain (A) on a similar mesh, which captures the mixing of methane and air in full detail. This is important because the properties of this mixed-mode flame strongly depend on the degree of mixing. Note that the coupling between domain (A) and (B) is one-way. This allows to run the two simulations consecutively and use the inlet data even for different numerical setups like LES. This is justified since pressure waves that would be reflected in domain (B) and travel upstream to domain (A) are filtered out at the sponge region near the outlet of domain (B).
The simulations have been performed on the German national supercomputer Hazel Hen at the High Performance Computing Center Stuttgart utilizing up to 28 800 CPU cores in parallel. For further information about parallel efficiency and numerical performance of the code see [42]. OpenFOAM was used in version 5.0 and v1712, which were the most recent versions when the simulations were performed.
The reaction mechanism used in the reactive quasi-DNS in domain (B) (see Fig. 1) is an analytically reduced mechanism based on computational singular perturbation (CSP) for the combustion of methane in air developed by Lu et al. [48] based on the GRI 3.0 [50] reaction mechanism. It contains 19 transported chemical species and 11 additional species which are assumed to be in steady state (QSS) and therefore their concentrations are computed algebraically. Internally, the rate coefficients of 184 reactions are computed based on the Arrhenius rate law. The mechanism can be used over a wide range of pressure and equivalence ratios and is able to capture different properties like ignition delay times and flame speed with a worst-case error of 10 % [48]. It is therefore suitable for the mixed-mode Sydney flame. A direct comparison with the GRI 3.0 reaction mechanism is given in the appendix in Section A.3.

Grid resolution
Since the simulation of the Sydney flame is a quasi-DNS without any subgrid turbulence or combustion models, both the flow and the reaction zone of the flame have to be resolved sufficiently. The Reynolds number at the inlet of domain (B) based on mean velocity and the diameter of the central pipe with D = 7.5 mm is Re = 2.4 · 10 5 . Although this Reynolds number is quite high, it should be noted that a hot pilot flow with T = 2231 K is present at the very beginning of the inlet plane, which has a Reynolds number of well below 1000 due to the high temperature and therefore high viscosity. This hot pilot confines the region where the flow is still cold and unburnt and therefore exhibits a small Kolmogorov length only near the symmetry axis. The surrounding co-flow itself has a low velocity of 15 m/s so that there is little need to resolve this outer region finely. The central methane-air flow leaving the nozzle is also not wall-bound anymore (see Fig. 1) and develops further as a free jet. Additionally, although the whole domain for the flame simulation has a length of 68 D, sufficient resolution is provided in the range of 0 ≤ x/D ≤ 20 where most of the experimental measurements are performed. The wall-near flow is resolved with y + = 0.58 for the non-reactive premixing simulation in domain (A). Due to the use of a purely hexahedral mesh, this resolution also extends downstream in the shear layer. The circumferential resolution φ + near the wall in wall units is about 8, and it is 6 in axial direction x + . The wall-near flow of partially premixed methane and air upstream of domain (B) in the central pipe is resolved with y + = 1, x + = 5, φ + = 10 and the pilot flow in the outer pipe is resolved with y + = 0.32, x + = 1.5, φ + = 1.9.
The mesh consists of purely hexahedral cells which are refined towards the flame reaction zone and the core region of the flow where the gas is cold and unburnt. Since the highest gradients are expected in radial directions, especially in terms of reactive scalars, the size of the cells in radial direction is generally smaller than in axial direction.
A direct comparison of grid resolution versus Kolmogorov length is given in Fig. 2. The radial distribution of the Kolmogorov length (solid lines) is shown at different axial positions. It is approximated from where the overline symbol indicates Reynolds averaging. The grid resolution is given in terms of the cube root of the cell volumes V c and shown as dashed lines. Near the symmetry line at the center of the domain (r/D < 1) and close to the inlet (x/D < 5), both Kolmogorov length and cell size are smallest due to the high velocity and low temperatures. The cube root of the cell volumes gives a length that is slightly larger than the Kolmogorov length but still within a factor of two ( 3 √ V c /l K < 2). It should be noted that the finest resolution is in the radial direction, which is lower than 3 √ V c and therefore is able to resolve the Kolmogorov length at the center in radial direction. At r/D > 0.5 the cold methane-air flow from the central pipe comes into contact with the hot pilot and chemical reactions are initiated, leading to a rapid increase of the Kolmogorov length with temperature. The computational mesh however is still relatively fine in order to resolve the reaction zone of the flame. At the outer regions r/D > 2, the flow is dominated by the co-flow which is cold

Resolution of the flame front
Because the focus of this work is the creation of a database for investigating flame dynamics of flames with mixed combustion mode such as flame stabilization and flame regime transition, the resolution of the flame front by the computational mesh is important. Figure 3 shows the volumetric heat release rateQ on a 2D cutting plane. The highest values ofQ and therewith thinnest flame fronts appear upstream near the inlet. Therefore, a zoomed view into flame front near the inlet is given together with the grid lines showing the mesh resolution. One way to estimate the resolution of the flame front is comparing the cell sizes with the thermal thickness δ of laminar stretched flames [43]: A worst case estimate for the flame thickness can be computed by taking the highest stretch rates K = 1/A dA/dt, where A is the surface area of an infinitesimal area element on the flame front, of about 5 · 10 3 s −1 , that appear in the quasi-DNS, and compute laminar premixed counterflow twin-flames or diffusion flames with the same stretch rate. Both setups yield similar laminar flame thicknesses of about 330 μm at this stretch rate. The computational grid has cells with a radial resolution of about 10 μm to 50 μm in the regions where reaction zones can appear, leading to a worst case resolution of the flame front with 8 to 30 cells. Since this is based on the highest stretch rates, the flame front will be resolved with more cells on average. Additionally, the reaction zones become thicker downstream as the flame shows strong characteristics of purely premixed flames only near the inlet. Another example is given in Fig. 4, which shows the profile of instantaneous heat release rate along a radial line at x/D = 1, where the flame front is thinnest. Each point represents the intersection of the line with a cell face and therefore illustrates the resolution of the flame front.

Validation with Experimental Data
After having demonstrated the adequacy of the numerical settings and computational mesh for the quasi-DNS of the Sydney flame, the simulation results are compared with experimental data in this section. There are two different sets of measured data available for the Sydney flame in configuration FJ200-5GP-Lr75-57: One from 2013 (I2013-5GP), which will be denoted as Exp. 2013 in the following, and one from 2015 (I2015-5GP), which will be denoted as Exp. 2015 [7]. The measured data comprise of time averaged and instantaneous data for temperature and species mass fractions. The data from Exp. 2015 also contains velocity measurements.

Mean and RMS values
Time averaged values and their root-mean-square (RMS) have been obtained from the simulation by Reynolds averaging quantities of interest over 40 ms or 60 flow-through times. During the last 10 ms, averaged profiles changed by less than 2 % on average. Figure 5  shows a comparison between the simulation and the two experimental data sets for radial profiles of the time averaged temperature and species mass fractions as well as their RMS at different axial locations. At x/D ≥ 5, the simulation results for both mean and RMS profiles mostly lie on or between the two measured data sets which illustrate the experimental uncertainties. Additionally, the flame in the experiments is not perfectly symmetrical. This can be seen at the inner flame region (low values of r/D) where two values of measured data appear for the same radius, contributing to the experimental uncertainties (for example the green circles representing the O 2 profile at x/D = 12). The biggest differences between measurement and simulation are upstream near the nozzle at x/D = 1 (left column in Fig. 5). On the one hand, the Bilger mixture fraction Z [51] is higher in the simulation in the core region at r/D = 0 and x/D = 1 compared to the experiments. These differences however vanish further downstream. On the other hand, temperatures at x/D = 1 are higher and H 2 and CO mass fractions are much lower in the simulation. The reason for this deviation comes from the hot pilot flow: in the simulation, the pilot is assumed to be in chemical equilibrium for a methane-air mixture at φ = 1 and an initial temperature of 300 K. This however is not the case in the experiment. Figure 6 shows the measured mixture fraction profile Z exp and the measured temperature profile T exp near the inlet at x/D = 1. Based on the measured mixture fractions Z exp , the adiabatic flame temperature can be computed from the chemical equilibrium of that mixture. This is shown as T ad . For r/D < 0.5 the flame cannot ignite due to the lack of a heat source and the very high mixture fractions outside the ignition range which explains the difference between T exp and T ad . At r/D ≈ 1 however, the pilot flow should be relatively undisturbed and approximately reach the adiabatic flame temperature. Because the simulation assumes perfectly burned methane-air for the pilot while the experiments show a still not fully burned pilot, intermediate species profiles and temperature near the pilot exit at x/D = 1 in Fig. 5 deviate. Additionally, the pilot flames in the experiment consists of a five component mixture with the same elemental composition and adiabatic flame temperature as methane-air at φ = 1. These five components include not only methane as fuel but Fig. 6 Reynolds averaged profiles of temperature and mixture fraction from the experiment near the inlet at x/D = 1. The measured temperature profile T exp is compared to the adiabatic flame temperature T ad , which is computed from the measured mixture fraction profile Z exp and T 0 = 300 K at constant enthalpy and pressure also hydrogen, leading also to the large deviations in the H 2 profiles. Further downstream (x/D ≥ 5) these differences disappear as the pilot flow approaches chemical equilibrium in the experiments and the influence of the combustion of the inhomogeneous methaneair from the central pipe flow increases. Although the correct pilot composition could have been included in the quasi-DNS, the simplified composition was chosen to make comparisons with different LES combustion models, which may not be able to include this, easier in the future. For the same reason, heat conduction in the pipes is not considered in the quasi-DNS.
The experimental data from Exp. 2015 also include sparse data for the velocity fields, which are compared with the simulation results in Fig. 7. Both mean and RMS values show again very good agreement, except at the central location at x/D = 1. This may be attributed to the velocity block profile used for the pilot flow, causing stronger dissipation of the central jet.

Scatter data
Not only time averaged profiles are available from the measurements but also instantaneous scatter data, presented in Fig. 8: Each point in the columns on the left and in the middle show instantaneous temperature conditioned on mixture fraction at a single location on an axial plane. In order to ensure comparability, the same number of points from the simulation (left column) as in the measurements (middle column) sampled at the same radial locations as in the experiment which are not uniformly distributed, are depicted. Again, the distributions of instantaneous temperature over mixture fraction values have their largest differences at x/D = 1 while the distributions further downstream become very similar. This shows that the simulation yields the correct correlation between chemical scalars and control variables across the different flame regimes (premixed, nonpremixed) as well as quenching and ignition events. The right column of Fig. 8 shows the same values of temperature as from the scatter plots conditioned on mixture fraction: All temperature values are sorted into mixture fraction bins of Z = 0.005 and their average and RMS values in each bin are computed. Again, x/D = 1 shows the largest differences as explained in the previous section while the agreement further downstream is very good. In general, mean values show a better agreement with the measurements from Exp. 2015 (blue dots) than Exp. 2013 (black dots). At x/D = 1, a sharp bend in the conditioned mean temperature profile at Z ≈ 0.05 can be seen. This is the radial position where the central pipe wall at the inlet plane is located.  The simulation results from the quasi-DNS of the mixed combustion mode Sydney flame have been compiled into a database to allow further investigation in future work. Figure 9 shows an example of simulation results in the database. Depicted are 2D cuts through the full 3D dataset at a single time step showing solution fields of the simulation like temperature and velocity as well as post-processing data like Bilger mixture fraction Z, reaction progress variable c (16) and flame index FI [52]  Since the aim of this work is to present the quasi-DNS for the generation of the database and to show its validity, a detailed analysis of the database with respect to flame stabilization and flame regime transition phenomena is left as future work. For instance, the data may be used to investigate the regime transition between premixed and nonpremixed combustion mode and how this affects the statistical distribution of common control variables like mixture fraction and reaction progress. Figure 10 on the left shows the joint probability density function (JPDF) of the mixture fraction Z and reaction progress c at a single location which shows a distinct correlation. Therefore, Z and c cannot be regarded as independent parameters, which is often assumed to construct joint PDF models [53]. This data can for example be further conditioned on the flame index F I , in order to separate the contributions of premixed and nonpremixed combustion modes to this distribution. Other examples would be to investigate systematically different definitions of the reaction progress variable and how this affects the JPDF. Another goal for future works is to compare the results from the quasi-DNS to LES with different combustion models to test whether common assumptions from LES combustion models hold for this type of flame and to develop improved combustion models for mixed combustion regimes.

Conclusion
A model-free quasi-DNS has been performed of the turbulent Sydney flame in configuration FJ200-5GP-Lr75-57, which is actively being investigated, and the results have been compiled into a database. Because the simulation was conducted using finite rate chemistry from a complex reaction mechanism as well as detailed molecular diffusion models for multispecies mixtures which take preferential diffusion into account, the simulation results allow studying the flame dynamics in this mixed combustion mode flame in great detail, particularly due to the very high resolution of the reaction zone. The database provided with this paper will allow to test the validity of different modeling approaches for inhomogeneous flames and study regime transition events.
The simulation was performed with a self-developed reactive flow solver in OpenFOAM. The quasi-DNS capabilities of OpenFOAM are discussed with a focus on spatial discretization schemes. Canonical test cases included in this paper like the Taylor-Green vortex show that OpenFOAM's cubic discretization scheme yields better accuracy and lower numerical diffusion compared to classical central difference schemes and can achieve higher-than second order convergence rate in some cases. On the same mesh as a spectral DNS, Open-FOAM was able to compute dissipation rates within 1 % maximum deviation of the spectral DNS solution. Furthermore, the reactive flow solver implementation is validated by demonstrating its ability to compute the correct flame profiles in canonical 0D and 1D setups.
The simulation of the Sydney flame itself is validated by demonstrating that the grid resolution is sufficient for a model-free simulation, both in terms of the resolution of the flow as well as the flame front. This is also reflected in comparisons with two data sets from experiments which are in quantitatively good agreement with the simulation for scatter data as well as averaged and RMS values for temperature, concentrations and flow velocity. Deviations near the inlet are explained by the different pilot flow compositions in the simulation and experiments.
The numerical setups for the test cases from the Appendix as well as the database itself can be downloaded here: http://vbt.ebi.kit.edu/index.pl/specialtopic/DNS-Links. The numerical setup consists of a one-dimensional domain with equidistant cells of size x and a total length L. The setup was computed using OpenFOAM's standard solver scalarTransportFoam in version v1812 for varying cell sizes (or Péclet numbers Pe) and discretization schemes. Figure 11 shows results of the relative L 1 error after the results have converged up to a relative residual of 10 −7 .
Here, N is the number of grid points. Every point in Fig. 11 represents a 1D simulation. The rate of convergence can be estimated from the slope of the relative L 1 error in Fig. 11 from d log(L 1,rel ) d log( x/L) . As expected, the rate of convergence of the first order upwind scheme is 1, while the second order central difference (CD) scheme is second order. The cubic scheme also exhibits a second order convergence rate but with much lower numerical dissipation, hence the smaller relative L 1 error compared with the central difference scheme.

A.1.2 2D Taylor-Green vortex
A more complex test case is the canonical 2D Taylor-Green vortex (TGV) [55]. Compared to the previous test case, this case is transient and two-dimensional and the full set of incompressible Navier-Stokes equations are solved with the standard OpenFOAM solver pimpleFoam in version v1812. The numerical setup consists of a 2D domain built-up from equidistant cells with a height and width of L × L, where all boundaries are periodic. The analytical solution for this is: p(x, y, t) = ρu 2 Flow, Turbulence and Combustion  Figure 12 on the right depicts the field of the x-component of the velocity. The simulations are performed with the first order upwind scheme, second order central difference scheme and cubic scheme for all spatial derivatives and with varying number of cells or cell sizes x, respectively. The CFL number is held constant at 0.03 (corresponding to t = 0.1 ms for 63 2 cells), which is a commonly used value for the 2D TGV [56].
The solution of u x at t = 10t c along the center line (dashed line from Fig. 12) for a grid consisting of 63 × 63 cells is depicted in Fig. 13 along with the analytical solution. The simulation using the first order upwind scheme shows large differences to the  analytical solution. The benefit of using the cubic scheme over the classical second order central difference scheme is also evident: The difference between the simulation and the analytical solution at x = 1 3 L is below 1 % with the cubic scheme but over 10 % with the central difference scheme. In order to achieve the same accuracy with the central difference scheme as the cubic scheme on 63×63 cells, over 500×500 cells would be necessary. As shown in Eq. 6, the cubic scheme is a strictly more accurate version of the central difference scheme with the drawback of less stability due to the explicit correction term.
The rate of convergence based on the x-component of velocity obtained from the cubic scheme is shown in terms of the relative L 1 error in Fig. 14 on the left and in terms of the L 2 error on the right.
The rate of convergence can again be estimated from d log(L 1,rel ) d log( x/L) or d log(L 2 ) d log( x/L) , respectively. Since this case is transient and includes the solution of the full set of Navier-Stokes equations, the convergence behavior is different from the simple 1D steady-state case in Fig. 11. For coarse grids, the rate of convergence of the cubic scheme is slightly better than second order while for fine grids and small time steps it may even reach fourth order.

A.1.3 3D Taylor-Green vortex
The benefit of using the cubic scheme compared to the central difference scheme in terms of accuracy is demonstrated in this section with the 3D Taylor-Green vortex case. The setup is similar to the 2D TGV case except that the domain is a cube with side length of 2πL where all six faces are set to cyclic boundary conditions. The counter rotating vortices in the corners start to decay and form a turbulent flow field. This is shown in Fig. 15. The simulations are performed with the standard pimpleFoam solver in version v1712, which solves the incompressible Navier-Stokes equations, using upwind, central difference and cubic schemes for the spatial discretization. The maximum CFL number during the turbulent with u 0 = 1 m/s, L = 1 m, p 0 =0, ρ = 1 kg/m 3 . Viscosity is set to 1/1600 m 2 /s. Since there is no analytical solution for the 3D TGV, simulation results from OpenFOAM and those from a spectral DNS code [57] are compared. The OpenFOAM simulations are performed on the same numerical grid as the spectral DNS, which consists of 512 3 equidistant cells. The comparison is shown in Fig. 16. On the left, the volume averaged kinetic Fig. 16 Normalized volume averaged kinetic energy (left) and normalized dissipation rate (right) during the turbulent decay of the vortices from the 3D Taylor-Green vortex case energy E k inside the domain using the different discretization schemes is presented and compared with the spectral DNS solution.
At this resolution, there is only little difference to the spectral DNS, except with the upwind scheme. The advantage of the cubic scheme becomes clearer when looking at the volume averaged dissipation rate ε: The dissipation rate demonstrates that the first order upwind scheme is not able to reproduce the turbulent decay correctly. The second order central difference scheme underestimates the peak dissipation rate by 5 %. The cubic scheme, however, is able to correctly reproduce the evolution of the dissipation rate and to capture the correct values for the peak dissipation rates (difference below 1 % compared with the spectral DNS). The test cases discussed in this section, from the simple 1D convection-diffusion equations to the 2D and 3D Taylor-Green vortex cases, clearly demonstrate that the cubic discretization scheme implemented in OpenFOAM offers more accurate solutions compared to the classical second order central difference scheme on hexahedral meshes and may even reach convergence rates higher than two. Therefore, the cubic scheme was chosen for all spatial discretizations for the quasi-DNS of the Sydney burner, which will be described in the following sections. One drawback of using the cubic scheme is a reduced numerical stability due to the explicit correction method from Eq. 6 so that low CFL numbers, high quality meshes and well resolved solution fields are necessary.

A.1.4 Decaying homogeneous isotropic turbulence
The numerical setup for a decaying homogeneous isotropic turbulence case is a cubic mesh consisting of 256 3 cells with a side length of 0.9π m. All boundary conditions are periodic. The initial velocity field is set [58] according to a prescribed kinetic energy spectrum, which corresponds to the case t.U 0/M=42 from [59]. OpenFOAM's pimpleFoam Fig. 17 Volume averaged enstrophy over time for decaying homogeneous isotropic turbulence in a box. Simulation performed with and without viscosity solver is used, which simulates incompressible fluids, together with the cubic discretization for spatial derivatives. The simulation was performed two times: once with a viscosity of ν = 10 −5 m 2 /s and once with ν = 0. Figure 17 shows the volume averaged enstrophy over time. In the case with non-zero viscosity, enstrophy reduces over time due to the dissipation of kinetic energy. In the inviscid case, however, volume averaged enstrophy first increases as explained in [60]. As the turbulence decays into smaller length scales over time, numerical dissipation leads to a reduction in enstrophy as the Kolmogorov length goes to zero because ν = 0.

A.2 Validation cases for the reactive flow solver
In this section, simulations of two canonical flame setups are briefly presented and the results are compared with results from Cantera in order to validate the reactive flow solver implementation.
The first test case is the zero-dimensional isochoric auto-ignition of methane in air at 1 bar and an initial temperature of 1600 K. The mixture has an equivalence ratio of 1 with additional 1.7 mass-% H 2 and 0.2 mass-% H in order to speed up the ignition process. The reaction mechanism is GRI 3.0 [50]. The computational mesh in OpenFOAM consists of only one cell while Cantera's idealGasReactor model is used for the simulation. Since this is a 0D process, there are no effects of convection and diffusion so that this setup is suitable to validate the implementation of the procedure for calculating chemical reaction rates. The results from Cantera and the reactive flow solver in OpenFOAM are displayed in Fig. 18. On the left, the temporal evolution of the main species mass fractions is shown. Maximum deviation between the results from Cantera and OpenFOAM are below 0.01 %. The same holds for the right of Fig. 18, where some intermediate species mass fraction profiles during the ignition process are compared.
The second setup is the one-dimensional freely propagating premixed flame. In this case, convection and diffusion play an important role in addition to the chemical reactions. The computational domain in OpenFOAM consists of 8 000 cells with an inlet on the left and an outlet on the right. Unburnt gas at 300 K, 1 bar with a composition of methane-air at φ = 1 enters the domain at the inlet. Initially, the left half of the domain is set to the unburnt mixture conditions while the right half is set to the burnt mixture conditions. The reaction   Fig. 19, with the main species on the left and intermediate species on the right. The plots are aligned so that the maximum heat release rate is at x = 0 mm. The maximum deviation between the two solutions is about 0.5 %. This agreement cannot be achieved with standard OpenFOAM, since it lacks detailed diffusion models [61].

A.3 Reaction mechanism by Lu et al.
In this section, results from the reaction mechanism by Lu et al. are compared to GRI 3.0. A direct comparison in Fig. 20 with the GRI 3.0 mechanism for the 1D flame setup from Section 1 shows, that the species mass fraction profiles of the main species as well as intermediate species agree within 1 %. Therefore, the reaction mechanism by Lu et al. is able to reproduce the correct flame structure. The computing time is however reduced by a factor of 2.5 in comparison to the original GRI 3.0 reaction mechanism.

A.4 Scatter plots
Additional scatter plots comparing experimental measurements and simulation results for different quantities and at different locations (see Section 5.2). In order to ensure Fig. 20 1D premixed freely propagating methane-air flame at p = 1 bar, T = 300 K and φ = 1 comparing the species mass fraction profiles predicted by the GRI 3.0 [50] and Lu19 [48] mechanism (both computed with the reactive flow solver in OpenFOAM). Left: Main species. Right: Intermediate species comparability, the same number of points from the simulation (left column) and the measurements (middle column) are sampled at the same radial locations (which are not uniformly distributed).