A mesh-free framework for high-order direct numerical simulations of combustion in complex geometries

The multiscale nature of turbulent combustion necessitates accurate and computationally efficient methods for direct numerical simulations (DNS). The field has long been dominated by high-order finite differences, which lack the flexibility and adaptivity for simulations of complex geometries and flame-turbulence-structure interactions in realistic settings. In this work we introduce a new approach to DNS of premixed combustion, based on a high-order mesh-free discretisation in combination with finite differences, enabling high-order simulations in non-trivial geometries. The approach is validated against a range of two- and three-dimensional flows, both laminar and turbulent, and reacting and inert. The present method a) has the resolving power for DNS of both laminar flames and inert turbulence with comparable accuracy to high-order finite differences, b) can capture the dynamics of unsteady bluff body stabilised flames, and c) is capable of simulating flame-turbulence interactions, with results comparing qualitatively well with published data. This work paves the way for DNS of combustion in complex geometries, offering an alternative approach to methods based on structured grids with immersed boundaries, or unstructured meshes. Further studies with the present method are proposed, which will aid understanding of fundamental flame dynamics in non-trivial geometries. Planned developments in adaptivity and extension of the mesh-free construction to all three dimensions will increase the value of the method, and support the push towards DNS of real geometries.


I. INTRODUCTION
Turbulent combustion is multiscale in nature, with chemical time scales of the order of nanoseconds, and submillimetre flame fronts, interacting highly non-linearly with much larger scale hydrodynamic, turbulent, and acoustic structures.Consequently, any direct numerical simulation (DNS) of combustion needs to be both highly accurate and computationally efficient.The need for high-order schemes is widely acknowledged in the literature (e.g.[1][2][3][4]), and has led to the field being dominated by high-order structured-grid methods, in particular finite differences.For a comprehensive review of the state-of-the-art in DNS of combustion, I refer the reader to [4].
As we move towards hydrogen fuels, DNS will continue to play a central role in improving our understanding of fundamental flame dynamics, providing insight impossible to obtain from experiments.Flame-turbulence-structure interactions are one key area of relevance, however current methods lack the ability to handle complex geometries.In the arena of high-order finite differences, codes such as SENGA+ [5], S3D [6] and KARFS [7] have had significant impact on combustion DNS.These are based on tenth-order (SENGA+) and eighth-order (S3D and KARFS) central finite difference schemes, and numerically solve the fully compressible reacting Navier-Stokes equations using explicit Runge-Kutta schemes for time advancement.Whilst extremely fast and accurate, these codes are limited to simple geometries, and lack the flexibility to simulate complex flame-structure interactions, whilst their construction largely precludes adaptivity.Body-fitted or block-structured meshes can be implemented for finite difference codes, but only in geometries of limited complexity.Another relevant structured-mesh code is the finite volume code HAMISH [8], which uses unstructured refinement on a Cartesian grid, along with fourth-order flux reconstruction to give an adaptive code with accuracy levels sufficient for DNS of turbulent combustion, though this remains limited at present to simple geometries.Complex geometries may be included in structured-mesh codes via immersed or embedded boundary methods, and in recent years several leading codes have benefitted from development in this area (e.g.S3D [9]), although these immersed boundary methods are generally limited to second order accuracy in practice.Simulations of turbulent combustion involve very fine flame fronts which dynamically evolve, but generally occupy only a small portion of the domain.This is a setting where adaptivity can bring significant benefits in reduced computational expense, and adaptivity is a highly desirable property for the next generation of combustion DNS codes [8].However, finite difference codes face challenges in this regard: variable resolution is limited to grid or coordinate stretching techniques, and adaptivity is generally not possible [10].Despite the advantages in accuracy and efficiency, the need for the next generation of combustion codes to include adaptivity and complex geometry capabilities has resulted in a move towards more flexible alternatives.
Unstructured-mesh methods provide an alternative with increased geometric flexibility, as a body fitted mesh may more easily be constructed to accurately represent the desired geometry.An example of a leading unstructured-mesh code in combustion DNS is AVBP, which uses the cell-vertex finite volume method, and although it is limited to low order spatial discretisation, it has been widely employed for simulations of combustion in non-trivial geometries (e.g.[11]).Codes based on the spectral difference method [12], in which high-order polynomial reconstruction of both solution and fluxes in each grid-element result in a highly accurate method on an unstructured mesh, show significant promise for the future of combustion DNS [4].
Where acoustics are not of interest, the variable density low-Mach approximation allows simulations with a larger time-step, potentially providing computational speedup.In the structured-mesh field, a key example of this approach is the adaptive-mesh code PELE [13] (and its subvariants), which solves the variable-density low-Mach multispecies Navier-Stokes equations with a second-order scheme, with recently developed embedded boundary capabilities [14].The elimination of the acoustic time-step constraint in the low-Mach approximation makes weak-form methods (which require solution of large, sparse linear systems every time-step) feasible.In this area, the leading Spectral Element method (SEM) code Nek5000 [15] has been adapted to low-Mach DNS of combustion in complex geometries (e.g.[16]).The crucial advantage of the SEM over finite difference and finite volume schemes is that it delivers extremely high accuracy on unstructured meshes and in complex geometries.However, where acoustics are of interest (as in the present case, and many practical combustion applications), weak-form methods such as SEM are prohibitively expensive.A second limitation, which applies to SEM and all unstructured-mesh methods, is the challenge of constructing a highquality body-fitted mesh.In complex geometries the construction of a high-quality mesh, where there is no significant loss of accuracy due to skewed mesh elements, can be an extremely resource-intensive task, and in some cases (e.g.porous media), can take longer than the simulation itself [17].This leads us to consider a fundamentally different approach: mesh-free methods.
Mesh-free methods are an alternative to grid-based or mesh-based methods, with Smoothed Particle Hydrodynamics (SPH) probably the most widely used (see [18] for a recent review).In the mesh-free literature, computational points are usually referred to as particles: an intuitive description as mesh-free methods are often constructed in a Lagrangian framework, where the collocation points move with the fluid.However, even in an Eulerian framework, mesh-free methods can have advantages over mesh-based methods, particularly for simulations on complex geometries.In mesh-based methods, the discretisation requires both the positions of computational points, and their topological connectivity.In mesh-free methods, we rely only on the relative positions of computational points.With no information on node connectivity required, discretising a geometry with an unstructured node-set (for example via propagating front algorithms [19]) is straightforward and easily automated.A limitation of many mesh-free methods is accuracy: for those in a Lagrangian frame, the construction of high-order interpolants at every time-step is prohibitively expensive, and hence the interpolation used to construct derivative approximations is typically low-order.For example, the standard SPH discretisation is formally zeroth order, or even divergent [20], and the resolving power falls far short of that required for DNS at reasonable resolutions [21].Even many of the proposed consistency corrections are only first, second or at most third order (e.g.[22][23][24][25][26][27][28][29]).
With these limitations in mind, I recently developed the Local Anisotropic Basis Function method (LABFM) [30], and demonstrated its potential with simulations of inert isothermal flows in [31].LABFM is a method for obtaining finite-difference style weights on a stencil of arbitrarily distributed unstructured nodes, with no information on node connectivity required.The weights are obtained by solution of local linear systems constructed to ensure that consistency is satisfied up to the desired order, with consistency between 4 th and 10 th order demonstrated [31].In [30] we assessed the eigenvalues of operators constructed with LABFM, which show very low levels of numerical dissipation, whilst in [31], we showed that the resolving power of LABFM is comparable, and in some cases outperforms, classical central finite differences of equivalent order.Compared with other consistent mesh-free methods (for a review of which, I refer the reader to the introductions of [30,31]), LABFM may be formulated at higher order, and has a lower computational cost for a given level of accuracy.LABFM has significant potential in the field of combustion DNS, with accuracy and stability properties akin to high-order finite differences, and excellent scalability, whilst retaining the geometric flexibility common to mesh-free methods.
In this work I demonstrate the potential for mesh-free DNS of combustion in non-trivial geometries.With LABFM and high-order finite differences having a similar construction, I show how they can be combined for three-dimensional simulations where the geometry is homogeneous in one dimension (e.g.flow past a cylinder).With this combination of LABFM and high-order finite differences, I vary the order of the spatial discretisation is between 4 th and 10 th , whilst using an explicit 3 rd Runge-Kutta scheme for time integration.The method is validated against a range of canonical inert and reacting flow problems, showing the framework has the resolving power required for both combustion and turbulence DNS.
The remainder of the paper is set out as follows.In Section II I introduce the governing equations.In Section III I describe the numerical implementation of the governing equations within the present framework, including details the LABFM/FD based discretisation, filtering procedure, and boundary conditions.Section IV contains the results of numerical tests on two-and three-dimensional reacting and turbulent flows.In Section V I discuss the computational performance and parallel scaling of the method.Section VI is a summary of conclusions.Before proceding further, I make a brief comment on my notation.The letter m is used to denote the polynomial consistency of the numerical scheme.Subscripts α and β are used to denote chemical species.Of the Latin subscripts, i, j, and k are reserved for vector indices in Einstein notation (and where these indices are repeated, summation is implied), whilst a and b are reserved for node (collocation point) indices.All other symbols are defined at the point of first use.

II. GOVERNING EQUATIONS
I consider the compressible Navier-Stokes equations in conservative form for a mixture of N S miscible reacting chemical species where ρ is the density, u i the i-th component of velocity, p is the pressure, E is the total energy, τ ki are components of the viscous stress tensor, f i is the i-th component of body force, q k the k-th component of the heat flux vector, Y α is the mass fraction of species α ∈ [1, N S ], and V α,k is the k-th component of the molecular diffusion velocity of species α. ωα is the production rate of species α.The total energy is related to the other thermodynamic quantities by where h α is the enthalpy of species α.The system is closed with an equation of state where R mix is the gas constant for the mixture, R 0 is the universal gas constant, W α is the molar mass of species α, and W mix is the molar mass of the mixture.In ( 2) and (3), and hereafter, sums over α are taken to be over all species α ∈ [1, N S ].The viscous stress is defined as where µ is the dynamic viscosity and δ ik is the Kronecker delta.The molecular diffusion term in (1d) is modelled as a Fickian process with in which ρ Vα,k Y α is obtained from the diffusion model (either constant Lewis numbers or mixtured-averaged, details below), and the final term is a correction to ensure that the compatibility condition is satisfied.The heat flux vector is given by where λ is the thermal conductivity of the mixture.The temperature dependencies of thermodynamic quantities (heat capacity c p,α , enthalpy h α ) for each species take polynomial form fitting the standard NASA polynomials [32].The energy in (2) can be expressed as a nonlinear function of T and the conserved variables, which is solved using a Newton-Raphson method to obtain T .In the specific case where c p,α is independent of temperature, the expression for temperature becomes linear and the Newton-Raphson method converges exactly in a single iteration.I take the last known value of T as a suitable guess to initialise the Newton-Raphson method, and for typical flows (e.g.premixed H 2 -air flames), this results in convergence to a tolerance of 10 −10 within two or three iterations.
I consider two models for transport properties.The constant Lewis number approximation, and a mixture-averaged transport model.For the case of constant Lewis numbers, I model the diffusive flux as in which D α is the molecular diffusivity of species α.I prescribe a temperature dependence of the viscosity of the form with T ref a reference temperature and r T a constant exponent typically set to r T = 0.7, except where otherwise specified.The thermal conductivity λ and molecular diffusivities D α follow from the definitions of Prandtl P r and Lewis Le α numbers.
For the case of mixture-averaged transport, I model the diffusive flux by In this case, the transport properties and binary diffusion coefficients for each species are evaluated using polynomials (polynomial in ln T ) following [33], and the mixture-averaged transport properties are obtained from combination rules following [34,35] with the Hirschfelder-Curtiss approximation for diffusion [36].Soret and Dufour effects are neglected.I consider a general reaction mechanism consisting of M steps, written in the form with the inclusion of third bodies where appropriate.Here M α symbolises the chemical symbol of species α, and ν ′ α,l and ν ′′ α,l are the molar stoichiometric coefficients of species α in reaction l.The reaction rate ω α of species α is given by the sum of the molar production rate over all reactions.The system of equations which describes the evaluation of molar production rates is not unique to this work, wherein it closely resembles the reaction frameworks used in other codes (e.g.[5,8,13,37]).I omit details for brevity, but note that third bodies, backwards reaction rates and Lindemann forms are included where required.For more details I refer the interested reader to [38].Although I have tested my method on a range of reaction mechanisms, in this work I focus only on two: a single-step Arrhenius mechanism designed to match the flame speed of stoichiometric methane-air combustion, and the 21 step, 9 species hydrogen-air mechanism of [39].For hydrogen combustion, I use the mixture-averaged transport model, whilst for the single-step methane mechanism, I assume constant Lewis numbers.

III. NUMERICAL IMPLEMENTATION
The spatial discretisation is based on the Local Anisotropic Basis Function Method (LABFM) in the first two dimensions, and high-order finite differences in the third dimension.This allows us to simulate complex geometries, although I am limited to those in which the geometry is homogeneous in the third dimension.I could implement a fully mesh-free formulation, with LABFM used for derivatives in all three dimensions, at the cost of increased stencil sizes.This is an area of development for the author.LABFM has been detailed and extensively analysed in [30,31], and I refer the reader to these works for a complete description, where the congergence properties of the method are comprehensively demonstrated.For convenience of exposition, I denote as A 12 any plane in which x 3 is constant, and A 12 (C) the specific realisation of A 12 in which x 3 = C.
The domain is discretised in A 12 (0) with an unstructured point cloud of N 12 nodes.The discretisation in the third dimension is uniform: there are N 3 copies of the A 12 (0) discretisation, equispaced along the x 3 axis.Hence, there are a total of N = N 12 N 3 nodes in the domain.Each node a ∈ [1, N ] has position x a,i , a distribution lengthscale s a , and a computational stencil lengthscale h a .s a is the average node spacing in A 12 (x a,3 ) around node a, and is analagous to the grid spacing in a finite difference scheme.The spacing in x 3 is denoted s 3 , and is uniform throughout the domain, whilst s a may vary as a function of x 1 and x 2 .
The left panel of Figure 1 shows an example two-dimensional discretisation in the present framework for an aerofoil.Wall, inflow and outflow boundaries are discretised with a locally structured node distribution (detailed in Section III B), whilst the remaining space is filled with an unstructured node-set.The resolution is non-uniform, and in the present example s a varies by a factor of approximately 30 through the domain.The node distribution is generated using a variation of the propagating front algorithm originally developed by [19], following [31].I set s 3 to the weighted mean of s a taken over the domain as where n s takes the value 1 for the arithmetic mean, and 2 for the area-weighted mean, with the choice of n s case dependent.Unless explicitly specified, I set n s = 1 in this work.In some cases, setting s 3 to the average of s a may result in the discretisation in the third dimension having insufficient resolution to capture the smallest scales of the flow (either turbulent structures or flame fronts).In all cases in the present work, a manual check sufficed to reassure that this issue did not arise.Ultimately in the present framework there is a compromise that must be made: larger s 3 may under-resolve flames and turbulent structures in some regions, whilst smaller s 3 will increase computational costs, over-resolving regions of little interest.Including adaptivity, and extending the framework to be mesh-free in all three dimensions, is a key avenue I intend to explore in future.Each node holds the conservative variables ρ a , (ρu i ) a , (ρE) a and (ρY α ) a .The governing equations are solved on the set of N nodes.I denote the difference between any property (•) at two nodes by (•) ba = (•) b − (•) a .The computational stencil for each node a is denoted N a , and is constructed to contain all nodes b in the same plane A 12 (x 3,a ) such that |x i,ba x i,ba | ≤ 4h 2 a , (with the repeated subscript i implying summation), alongside all the nodes b with x 1,ba = x 2,ba = 0 required to build a finite difference stencil in the x 3 dimension.The right panel of Figure 1 shows a typical computational stencil (for the central node in blue) for a three dimensional simulation, with the nodes in the A 12 plane in yellow, and the nodes used to construct the finite difference stencil in x 3 coloured red.

A. Evaluation of spatial derivatives
In present framework, all spatial derivative operators take the form where d identifies the derivative of interest (e.g.d = x 1 indicates the operator L d approximates ∂ (•) /∂x 1 , and d = x i x i (summation implied) indicates the operator approximates the Laplacian), and w d ba are a set of inter-node weights for the operator.Note that the construction of (13) takes the same form for both LABFM-based and finite difference-based derivative approximations.For derivatives in the x 3 direction, the weights w d ba are simply central finite difference weights of order m.For derivatives in A 12 , LABFM is used to set w d ba , which is constructed from a weighted sum of anisotropic basis functions centred on node a. Note, derivatives at boundaries require a slightly different treatment, detailed in Section III B and Appendix A. The basis functions used herein are formed from bivariate Hermite polynomials multiplied by a radial basis function, and the basis function weights are set such that the operator achieves a specified polynomial consistency m.For polynomial consistency m, the errors in the approximation of first derivatives are O (s m ), whilst the errors in the approximation of second derivatives are O s m−1 .Calculation of these weights is done as a pre-processing step, and described in detail in [31].
The order of consistency of the scheme is specified between m = 4 and m = 10, and although there is capability for this to be spatially (or temporally) varying, in this work I set m = 8 uniformly away from boundaries (except where explicitly stated in my investigations on the effects of changing m).This provides an appropriate compromise between accuracy and computational cost.At non-periodic boundaries, the consistency of the LABFM reconstruction is smoothly reduced to m = 4, as detailed in Section III B. The stencil scale h a is initialised to h a = 2.7s a in the bulk of the domain, and h a = 2.4s a near boundaries (choices based on experience from previous studies as being large enough to ensure stability [30]).In bulk of the domain, the stencil scale is then reduced (effectively reducing the support radius) following the optimisation procedure described in [31].This has the effect of both reducing computational costs, and increasing the resolving power of LABFM, and h a subsequently takes a value slightly larger than the smallest value for which the discretisation remains stable with the non-uniformity in h a /s a accounting for the differences in local node distribution within each stencil.
The evaluation of derivatives using (13) constitutes the bulk of the computational cost, and hence it is desirable to use analytic expressions to relate derivatives of secondary properties to those evaluated via (13) where possible.First derivatives of ρ, u i , ρE and Y α , alongside the temperature T and pressure p, are evaluated directly using (13).Convective terms are then constructed from combinations thereof: for example Laplacians and direct second (i.e.∂ 2 (•) /∂x i ∂x i ) derivatives are also evaluated with (13).To avoid the explicit evaluation of cross-second derivatives of velocity, the viscous stress divergence is evaluated as where ∂ 2 u k /∂x i ∂x k is the gradient of the velocity divergence, which is calculated through two iterations of first derivative evalations with (13).The benefit of this approach is that it avoids the need to explicitly evaluate the derivatives ∂ 2 /∂x 1 ∂x 3 and ∂ 2 /∂x 2 ∂x 3 .This allows for a smaller stencil size, and provides a significant reduction in the complexity associated with parallelising the scheme.Gradients of secondary properties are evaluated from analytic expressions in terms of gradients T and Y α .The gradient of the enthalpy of species α is given by and the gradient of the mixture specific heat capacity is When constant Lewis numbers are assumed, additional relations are available to evaluate the spatial derivatives of transport properties: For the mixture averaged model, gradients of µ, λ and ρD α are obtained directly using (13).Even with these substitutions, for three-dimensional simulations the constant Lewis number transport model requires 3 (8 + N S ) first derivative evaluations, and 4 + N S Laplacian evaluations to construct the RHS of the governing equations.For the mixture-averaged transport model, this increases to 3 (10 + 2N S ) gradients evaluations and 6 + N S Laplacian evaluations.Hence any reduction in stencil size whilst retaining accuracy and numerical stability has the potential to significantly reduce computational costs, and the optimisation of stencils is an active area of research for the author.

B. Boundaries
The boundary framework extends the approach described for isothermal flows in [31] in two ways.Firstly, it is extended to accommodate thermal and reacting flows, utilising the formulation of [40].Secondly, for wall boundaries, a LABFM-based interpolation is used to partially replace the finite difference stencil, allowing for the band of uniformly distributed nodes along boundaries to be reduced from 5 rows to 3, permitting greater geometric flexibility whilst retaining the 4 th order consistency of the derivative operators at the boundaries.
The following discretisation procedure applies to derivatives in the plane A 12 only.All derivatives in the third dimension use high-order finite differences.Panel a) of Figure 2 illustrates the node distribution and computational stencil used near wall boundaries.Nodes are placed along the boundary (black nodes).For every boundary node a 0 , and additional node a 1 (yellow, in the left panel of Figure 2) is placed with position x a1,i = x a0,i + s a0 n a0,i for where n a0,i is the i-th component of the unit vector normal to the boundary directed towards the fluid, and a node a 2 (shown in dark blue) is placed with position x a1,i = x a0,i + 2s a0 n a0,i .For inflow and outflow boundaries, I place a further two nodes (denoted a 3 and a 4 ) along the normal vector to complete the 5 point finite difference stencil.This formulation for inflow and outflow boundaries is as detailed in [31], and is illustrated in Figure 2 b).For wall boundaries, I don't complete the finite difference stencil, instead discretising the remainder of the fluid with an unstructured node set (nodes coloured light blue near the boundary, and white a distance of more than 4.5s a from the boundary, in panel a) of Figure 2).To construct the finite difference approximation to derivatives normal to the wall boundary on nodes a 0 and a 1 , I interpolate fluid properties onto a set of two fictitious nodes (again denoted a 3 and a 4 , coloured red in panel a) of Figure 2) to complete the stencil.For clarity, the governing equations are not evolved on these fictitious nodes -these nodes are only used for interpolation to construct the finite difference stencil.Properties are interpolated from regular nodes onto the fictitious nodes a 3 and a 4 .I refer to this scheme as an interpolated finite difference scheme.Full details are given in Appendix A. Details of the discretisation scheme used for different nodes near the boundary are provided in Table I, where the distance between a node and the boundary is denoted d a .Beyond the structured nodes a 0 and a 1 , I use LABFM for all derivatives in the plane A 12 .For nodes a 2 I set m = 4, and for nodes with 2 < d a /s a ≤ 4.5 I impose the limit m ≤ 6.For derivatives tangential to the boundary at nodes a 0 and a 1 , I use the one-dimensional form of LABFM described in [31], to which I refer the reader for complete details.For nodes a 0 , the stencil for tangential derivatives involves only other nodes on the boundary (i.e. the solid black nodes in Figure 2), and for nodes a 1 , the stencil contains only nodes a distance s a from the boundary (i.e. the yellow nodes in Figure 2).With these stencils, difference operators are constructed using a one-dimensional formulation of LABFM as detailed in the appendix of [31].On boundaries, derivative operators are constructed in the boundary frame of reference, whilst in the first row, a rotation is applied to transform them into the Cartesian frame, as described in [31].Where the resolution is non-uniform along the boundary (as illustrated in panel b) of Figure 2), an additional rotation is applied to the first derivatives evaluated on node a 1 .
In terms of the numerical boundary conditions, no-slip walls (either isothermal or adiabatic), inflows (either nonreflecting or "hard" with specified u i and T ), and non-reflecting outflows are imposed through the Navier Stokes characteristic boundary condition (NSCBC) formalism (see, e.g.[38]).For non-reflecting inflows, I include the transverse relaxation terms following [41].As is standard in the NSCBC formalism, non-reflecting outflows are constructed to track a desired pressure p out , introducing a small degree of acoustic reflection in return for ensuring the outflow pressure does not drift over time.For outflows, I follow the formulation of [40].The improved outflow conditions of [41] have been implemented, but were found to be less stable than the method of [40], particularly when handling passage of flames at angles oblique to the boundary.Further investigations are required in this area.

C. Temporal integration and filtering
With the right hand sides of (1) evaluated using the derivative operators described in Section III A, and chemical source terms evaluated explicitly, the governing equations are integrated in time using an explicit Runge-Kutta scheme.I use a four-step, third-order, low-storage scheme, with embedded second order error estimation, denoted RK3(2)4[2R+]C in the classification system of Kennedy et al. [42].The value of the time step is set adaptively using a PID controller to keep time-integration errors below a specified threshold (set to 10 −4 in throughout the present work).As is common in high-order collocated methods, the low levels of numerical dissipation in the spatial derivative operators result in a scheme in which high-wavenumber modes of instability can develop.To avoid this, the solution is de-aliased after each time-step by high-order spatial filtering using the filters described in [31] combined with a finite difference filter following [43].The combined filter is defined for even m by where κ m is a filter coefficient, specific to each node, calculated as in [31], and ∇ m A12 is the two-dimensional Laplacian raised to the power of m/2 (as used in the filtering in [31]), defined in index notation as The term within the square brackets in (19) is the two-dimensional filter in A 12 , whilst the operator outside the square brackets is the filter in the x 3 dimension.The filters are applied sequentially in this manner to ensure high wavenumber modes of all orientations are suppressed, rather than just those in and orthogonal to A 12 .This approach follows that used in finite difference codes, where the three coordinates are filtered sequentially.With the LABFM filter, the entirety of A 12 is filtered in one pass, because the filter includes cross-derivative terms, and is hence able to suppress all high wavenumber modes in A 12 simultaneously, regardless of their in-plane orientation.For density, momenta and energy, I simply apply the filter as defined by (19), such that the filtered variables passed to the next time step are As with the molecular diffusion terms, the filtering procedure does not guarantee satisfaction of the compatibility criteria (6).Hence, I introduce a correction to the filtering procedure for the species mass fractions, which takes inspiration from the diffusion correction velocity.For the species mass fractions, I modify the procedure defining where the filter correction F c is given by This ensures that where Y α = ρY α / ρ, is satisfied by construction.The weighting ρY α / ρ ensures that the correction is distributed across the species in proportion to their filtered mass fractions.Hence, as Y α approaches zero, the correction for species α also approaches zero, ensuring that a locally negative F c doesn't push Y α < 0 locally.Equation ( 22) can be straightforwardly re-arranged to give the correction in multiplicative form as where the correction term F ⋆ c is simply I use the correction procedure defined by ( 25) and ( 26) throughout this work, for which the error in satisfaction of ( 24) remains approximately 1 × 10 −18 at all times.

D. Parallelisation
The present implementation of the numerical scheme is accelerated with a shared-distributed framework, using both OpenMP and MPI.Non-uniform block-structured three-dimensional domain decomposition is used, with the physical sizes of subdomains set such that each MPI rank contains the same number of nodes.The decomposition takes as input the desired number of processors in each dimension N P 1 , N P 2 , and N P 3 , giving a total of N P = N P 1 N P 2 N P 3 processors.The allocation of nodes to processors then forms a preprocessing step, following the algorithm: 1. Generate the node discretisation in A 12 (0).
2. Re-arrange the nodes in order of increasing x 1 .
3. Split the discretisation into N P 1 slices each containing the same number of nodes.(c) For each block, create the x 3 discretisation by creating N 3 copies equispaced along the x 3 axis, and distribute into N P 3 processor domains.
This approach results in each processor containing the same number of nodes, and the variation in processor load arises only from variations in the number of neighbours N .Providing the number of nodes in each MPI rank is large enough that the time for MPI communications is small relative to the time for all other calculations (which scale almost perfectly with increasing N P ), the overall parallel efficiency of the scheme is excellent.In strong scaling tests, with N = 2 × 10 6 nodes, I found the scheme has a parallel efficiency of 96.5% for N P = 1024 (relative to a baseline of N P = 4).I note that for smaller N/N P , the limitation on the parallel efficiency is the load balancing, rather than the cost of MPI communications.Although each MPI rank contains the same number of nodes, the size of computational stencils N is not uniform (in particular, N is larger where resolution gradients are larger).The decrease in parallel efficiency at higher N P is because the load balancing is less effective as N/N P is reduced.For typical simulations, such as those presented in Section IV of two-dimensional bluff body flames and three-dimensional flame-turbulence interactions, I use N/N P ≈ 5000.Although the decomposition and load balancing are sufficient for the present work, the implementation of a dynamic domain decomposition framework, with on-the-fly load balancing, is a key area of interest for the author.

IV. NUMERICAL TESTS
With the convergence and stability properties of LABFM extensively investigated in [30], and the ability of the method to accurately simulate two-dimensional isothermal flows established in [31], I focus here on tests which demonstrate the performance of the framework for reacting and turbulent flows.

A. Laminar hydrogen-air flames
First I present the results of two-dimensional simulations of planar laminar freely-propagating stoichiometric hydrogen-air flames.The computational domain is rectangular, with length L x1 = L = 10mm and a narrow aspect ratio with L x2 = L/40.The domain is periodic in x 2 , with inflow and outflow boundaries in x 1 .In all cases, the temperature T in = 300K, velocity U = 2.05m/s, and mass fractions for a stoichiometric hydrogen-air mixture are imposed at the inflow boundary, and a partially non-reflecting outflow condition tracks an outflow pressure of p out = 10 5 P a.I use the 9 species, 21 step reaction mechanism of [39].The initial temperature, density, velocity and mass fraction fields are defined by an error function flame profile with width L/40.The simulation is run until a steady flame profile is obtained.Figure 3 shows an example of the computational domain and node distribution for this test, with the bottom panel showing the resolution s, the middle panel showing the temperature, and the top panel showing a close up view of the node distribution around the flame front, coloured by the mass fraction of Y H2O2 .H 2 O 2 is chosen for visualisation, as it has a very narrow distribution profile within the flame front.For this discretisation, I set s/L = 0.003 at the inflow and outflow, reducing linearly over a region of 0.1L to s/L = 0.001 in the central 0.1L of the domain, which corresponds to a resolution of s = 10µm.As can be seen in the top panel of Figure 3, this results in the structure of the flame being resolved by at least 10 nodes.
Figure 4 shows the spatial profiles of temperature (top left) and species mass fractions through a freely propagating hydrogen-air flame, obtained with the present method (black lines), and two reference codes.For the present method, the one-dimensional profiles are obtained by interpolating (using LABFM with m = 8) values from the two-dimensional set of nodes to a uniform set of points along x 2 = 0.The first reference code is SENGA+ [5], one of the combustion community's leading high-order DNS codes.SENGA+ is based on a tenth-order finite difference scheme, with a fourth-order explicit Runge-Kutta scheme for time integration.The set of equations solved by SENGA+ is very similar to those solved in my mesh-free framework, with only slight differences in the choice of combination rules for mixture-averaged transport properties.I also note there are differences between the present work and SENGA+ in the construction used for convective terms (SENGA+ uses a skew-symmetric form, the present work does not), and viscous terms.The simulations performed with SENGA+ use a uniform resolution of 10µm.The second reference code is Cantera [44], an open-source flame code, which in the present work is used to implicitly obtain steady one-dimensional flame solutions with an adaptive grid, and chemical reactions represented with the San Diego mechanism [45].To aid comparison, the data in Figure 4 have been translated in x 1 such that the flame fronts (taken as the location where Y H2 is half its stoichiometric value) from all three codes align.
In Figure 4 there is a very close match between all three codes for the profiles of temperature and mass fractions of the main reactant and product species (top panels).For the O and H radicals, there is a very close match.However, for the mass fractions of OH, HO 2 , and H 2 O 2 , the results from Cantera over-predict peak concentrations relative to both the present method and SENGA+.This may be due to the reaction mechanism used; with both my mesh-free method and SENGA+ I use the 9 species mechanism of [39], whilst for Cantera, the San Diego mechanism [45] is used.However, it could also be due to different formulations for molecular diffusion.The mixture-averaged formulation is used in all three codes, although there are differences in the combination rules used to evaluate D α , µ and λ between codes.For example, in the present work the power-weighted average descibed in [34] is used for viscosity, whilst in SENGA+ the combination rule of [46] is used.Additionally, in the present work, Soret and Dufour terms are neglected, whilst they are included in SENGA+.Whilst these differences in formulation appear to yield negligible differences in results for the stoichiometric hydrogen-air flames modelled here, for lean flames, they may be more significant.Notwithstanding these minor discrepancies with the solutions obtained from Cantera, with SENGA+ widely established as one of the leading codes for combustion DNS (see e.g.[4]), the excellent match between the present method and SENGA+ in Figure 4 provide a demonstration of the comparable ability of LABFM to simulate accurately simulate laminar flames.
I briefly explore the effect of the order of accuracy of my scheme, by simulating the same case, but changing m. Figure 5 shows the effects of changing m on the profiles of T and Y OH .In all cases, the resolution is uniform at s = 10µm.I test four orders of LABFM: m = 4, 6, 8, 10.For both the temperature (left panel) and OH mass fraction (right panel) profiles, the results are indistinguishable when looking at the entire flame (inset in both panels).However, there are localised differences just upstream of the flame front (main axes).For smaller m the resolving power of the scheme is lower, and the resolution is insufficient to accurately calculate the advection terms.This results in slight oscillatory behaviour just upstream of the flame front, with local minima in temperature, and a region of (unphysical) negative mass fraction of OH.For m = 4, these are clearly visible in Figure 5, whilst for m = 6 the magnitude of these oscillations is much reduced.For m ≥ 8, these oscillations are not present -the temperature increases monotonically through the flame, and Y OH ≥ 0 everywhere.I conducted simulations with an in-house one-dimensional finite difference code, constructed to solve the same set of equations as my mesh-free framework, but using high-order central finite differences instead of LABFM, and found the same behaviour with changing order of the scheme.Note that whilst negative mass fractions are unphysical, that the numerical method admits them is not a fundamental problem, and is common to this and other methods (e.g.central finite differences) which are not strictly monotonicity preserving.In such schemes (including the present method), negative mass fractions are avoided by using high-order and sufficient resolution.
Next I keep m = 8 and vary the resolution.In this test I set s to be uniform throughout the domain.The left panel of Figure 6 shows the profile of Y OH for a range of resolutions, from s = L/250 to s = L/2000, with the full flame shown in the inset, and a close up of the preheat zone in the main axes.The right panel of Figure 6  FIG. 4. The spatial profiles of temperature and species mass fractions through a freely propagating stoichiometric hydrogen-air flame, obtained with the present method (black lines), SENGA+ [5] (red lines), and Cantera [44] (blue lines).For the present method and SENGA+, the 21 step mechanism of [39] was used.For Cantera, the San Diego mechanism [45] was used.
the profile of the pressure field through the flame, for the same range of resolutions, alongside data obtained with SENGA+ (blue line).For s ≤ L/1000 the solution is converged, and discrepancies are between these finer resolutions are negligible.Conducting a simulation with a very fine resolution of s = L/4000 as a reference, I find discrepancies in the volume integrated heat release rate and the volume integrated kinetic energy are less thant 5.1 × 10 −4 for all s ≤ L/750.For larger s, small errors appear just upstream of the flame front, in the form of localised oscillations, similar to those observed in Figure 5 for m = 4, with Y OH becoming locally negative.The magnitude of these errors increases with increasing s (decreasing resolution).This observation is expected, and common to codes based on high-order central finite differences.It is a consequence of the discretisation scheme, resulting from the treatment of the convective terms: the scheme is not total variation diminishing -there is no upwinding or similar monotonicity preserving treatment.As in my investigation of changing m, I have performed additional simulations with an in-house finite difference code, which yielding the same oscillatory behaviour ahead of the flame for equivalent resolutions.As with the mass fraction profiles, the pressure field shows good agreement with SENGA+ for resolutions s ≤ L/1000, both in terms of the pressure jump across the flame (to within 0.5P a), and the profile through the flame.For coarser resolutions oscillatory behaviour is observed in the pressure field, and the magnitude of the jump across the flame changes too.These pressure oscillations are a consequence of the oscillations in the evolved mass variables (ρ, ρu i , ρE and ρY α ): the the temperature is obtained via the non-linear system described in Section II, and the pressure follows from the equation of state.I note that even for the coarse resolution of s = L/500, the maximum deviation in pressure is of the order of 3P a, which corresponds to a relative deviation (given atmospheric pressure of 10 5 P a) of 3 × 10 −5 .For all values of m, the simulation was run with a uniform resolution of s = 10µm.Note YOH has been scaled by a factor of 10 3 .
- The results of Figures 5 and 6 taken together demonstrate the benefits of high-order schemes for flame simulations.With a high-order scheme, equivalent accuracy can be obtained at lower resolution (and hence cost).
Finally, in addition to the results for hydrogen-air flames presented here, I have tested the method for stoichiometric methane-air flames, using the 16 species, 25 step reaction mechanism of [47].In the interests of brevity I omit detailed results for the methane-air flame, but note that I obtain qualitatively similar results: i.e. close agreement with Cantera for the major species, and very close agreement with results from SENGA+ for all species.

B. Two-dimensional bluff-body stabilised flame
With the ability of the present method to accurately simulate planar laminar flames established, I now move on to a more challenging test.The key benefit of the present method over many other high-order combustion DNS codes is the ability to handle complex geometries, and as a demonstration of this, I simulate a two-dimensional bluff-body stabilised flame.Although comparatively simple (a single obstacle in a channel), this geometry is beyond the reach of most high-order combustion DNS codes.Similar configurations have been previously studied by a number of authors, both in the context of low-order codes (e.g.Ansys Fluent [48,49]), low-Mach number codes (e.g.[50]), and high-order finite difference codes [51][52][53][54][55].I note that although the simulations in [51][52][53][54][55] were performed with a high-order finite difference code, the geometric flexibility remains limited: in those works, the domain decomposition was constructed such that the bluff body boundary lay on the boundary between the domains of different processors.A bluff body with different shape (e.g.circular), could not be simulated with such an approach.
For this test the configuration consists of a square two-dimensional domain with side length L x1 = L x2 = L = 10mm.At the left and right boundaries non-reflecting inflow and outflow (respectively) conditions are imposed.The domain is bounded at the top and bottom by no-slip adiabatic walls.A bluff body is located with its leading edge a distance L/5 = 2mm downstream of the inlet, and positioned centrally in the cross-channel direction.I consider two shapes of bluff body.First, a square body of side length d = L/20 = 0.5mm, on which the corners have been rounded to have a radius of curvature of d/8, and second, a circular body of diamter d = 0.5mm.A no-slip adiabatic condition is imposed on the surface of the bluff body.The corners of the square bluff body are fully resolved with a resolution of s min = 5µm, whilst for the circular bluff body I set s min = 10µm.The domain is discretised with a non-uniform resolution with s = s min at the bluff body surface, increasing to s = s max = 30µm away from the bluff body, in regions where combustion is not expected to occur.In regions combustion is anticipated to occur (i.e.near the bluff body, and directly downstream), the resolution is restricted to s ≤ 10µm.This resolution matches that used in [53], and ensures the flame front is resolved by at least 15 nodes.Although no closed analytic expression is available to describe the variation of s a with x 1 and x 2 , input files containing lists of nodes and resolutions are available on request.The entire domain is discretised with approximately 2 × 10 5 nodes, a significant reduction compared to the 10 6 nodes required for a uniform resolution of s = 10µm, and 4 × 10 6 required for a uniform resolution of s = 5µm.The simulations are run on 200 cores spread across 100 MPI ranks.The domain, node distribution, MPI decomposition, and resolution are illustrated in Figure 7.
Following [53], the inflow composition is a hydrogen-air mixture with an equivalence ratio of Φ = 0.5, and the inflow temperature is set at T in = 300K.The partially non-reflecting outflow boundary is set to track a pressure of p out = 10 5 P a.The simulations is initialised with zero velocity, and ignition triggered by superposition of a Gaussian hot spot immediately downstream of the bluff body, with size 0.25mm, and peak temperature 2500K.The inflow velocity is then set to a fully developed (i.e.parabolic) profile with mean velocity which ramps linearly from zero to U = 30ms −1 over a period of T ramp = L/2U .For U = 30ms −1 , the Reynolds number based on the bluff body diameter of Re = 810.Simulations at sequentially higher Re are performed, each starting from the previous simulation, and again increasing the inflow velocity linearly over a period of T ramp = L/2U .Figure 9 shows the geometry of these flames, as defined by the isocontour along which Y H2O is half its equilibrium value, alongside the same measure taken from [54].For the square bluff body (solid black lines) I see an approximate match with the data from [54] (red lines), but with a slightly reduced rate of expansion of the flame with increasing streamwise distance, and small differences around the flame attachement point.my method yields an attachment point slightly further upstream than in [54], with a shallower (more swept back) profile near the bluff body.I postulate that the differences in attachement are due to the slight differences in geometry: my bluff body has rounded corners, whilst that in [54] is a square.The rounded corners on my geometry result in reduced flow separation at the sides of the bluff body, and in a higher energy boundary layer and hence the more swept back flame profile.The downstream differences in flame shape are likely to follow from the differences in the attachment.For the circular bluff body (dashed black lines), the flame attaches to the bluff body midway between the leading and trailing edges, and the expansion profile is shallower still.Although these differences are minor, they may have an impact on the dynamics of unsteady flames at higher Re, and serve to highlight the importance of simulations which accurately account for bluff body geometry.An in depth study of the effects of bluff body geometry on flame dynamics is planned for future work.Panels c) and d) of Figure 8 show the symmetric vortex shedding regime (c) and asymmetric vortex shedding (d), both at Re = 2160 (U = 80ms −1 ), for the square bluff body.I see qualitatively similar dynamics to those presented in [52][53][54].As the inflow velocity is increased, the steady flame develops mild fluctuations, which develop into symmetrically shed vortices, and then asymmetrically shed vortices.Regions of local extinction occur, and eventually these lead to global extinction.The mechanism appears to be the entrainment of cold gases into the wake of the bluff body due to the vortex shedding (near the point marked 1 in panel d) of Figure 8), as observed in [53].Whilst my simulations show flames with similar dynamics to those of [53], I do not match their values of U for the transition between regimes.In my work, I generally find the transition from steady to unsteady flames occurs at lower U than that found in [53].I believe this is due to the finite size of the domain.Whilst for the symmetric vortex shedding (and mild fluctuation) mode the fluctuations propagate convectively downstream, the instability appears to be an absolutely unstable one: that is, disturbances downstream can propagate upstream.The transition from a steady flame to the mildly fluctuating mode starts with small deviations from the steady flame far downstream from the bluff body, and the streamwise extent of these deviations then migrates upstream as they increase in amplitude.It is not clear whether this is due to the confinement (the flame is acoustically confined, even though hydrodynamically unconfined), or whether disturbances could propagate upstream even for a truly unconfined flame.Previous stability analysis of bluff-body stabilised flames (e.g.[56,57]) has relied on assumptions of low Mach number, unconfined flames, and parallel flow, and open questions remain on the nature of this instability.Even for an acoustically unconfined flame, acoustic waves are partially reflected from the flame itself, and disturbances could potentially propagate upstream between the two flame brushes.Figure 10 shows the pressure field around the flame in the asymmetric vortex shedding regime.Whilst the high pressure at the leading edge stagnation point, and the low pressure in the recirculating vortex just behind the bluff body are visible, there is also significant large scale structure in the developed pressure field due to the interaction between the unsteady flame and the acoustic field.
A further complication in this regard is the outflow boundary condition.Although able to handle the passage of an obliquely angled flame, the outflow boundary conditions (both in the present work, based on [40], and in e.g.[53], based on [41]) provide imperfect (though very close) approximations to the incoming acoustic waves.The passage of a flame through the boundary generates a small amount of acoustic noise, and it is possible that this triggers the transition from a steady to an unsteady flame.This noise is expected to be greater for flames at an oblique angle than those perpendicular to the boundary, as the boundary condition is constructed to relax towards a state with zero tangential velocity.In the present case, the steady flame solution has non-zero u x2 within the flame, and without knowing this a priori, designing the boundary condition to relax towards u x2 = 0 on the boundary is, although not strictly correct, the best available option.Furthermore, the outflow boundary conditions cannot predict the flame dynamics in the far field.It is possible that in reality an instability can develop far downstream and propagate upstream, but this would not be captured in the present simulations with a finite domain.
I have conducted additional simulations with the circular bluff body, in which the length of the domain is increased in the x 1 direction by a factor of both 1.5 and 2. I find a dependence of the flame dynamics on the domain length.Figure 11 shows isocontours of the heat release rate for flames attached to the circular bluff body at Re = 810, with the upper a), middle b) and lower c) panels showing the domain lengths of L x1 , 1.5L x1 , and 2L x1 respectively.All three simulations are started with zero velocity, uniform composition, and a Gaussian hot spot just downstream of the bluff body to trigger ignition.In all cases the inflow velocity is increased linearly to reach U = 30 at t = 0.25ms.The results in Figure 11 are taken for all three cases at t = 3.75ms.
For the shorter domain the flame is stable: the flame profile shown in panel a) of Figure 11 remains unchanged when the simulation is run for a further 5ms.For the longest domain, the flame is unstable, developing mild fluctuations which move convectively downstream.For the domain with length 1.5L x1 , the flame is close to the stability limit; perturbations grow in time, but very slowly.This can be seen quantitatively in Figure 12, which shows the time evolution of the volume averaged heat release rate for the three domain lengths.In all cases, the heat release rate starts at zero, and ramps up as ignition occurs and the flame develops.The rate of increase is different for the three domain lengths, because of the normalisation by the domain volume.As the flame reaches the outflow boundary, some perturbations are triggered by the boundary conditions (visible in the lower left inset in Figure 12) for all three domain lengths.The behaviour of these perturbations at later times is shown in the upper right inset of Figure 12.For the short domain with length L x1 , these perturbations are damped, and the volume averaged heat release rate reaches a steady value.For the longer domains these perturbations are less damped.For the domain with length 1.5L x1 , the perturbations grow, but extremely slowly, whilst the growth of the perturbations is clear for the domain with length 2L x1 .
None of the above discussion implies that the unsteady dynamics shown in Figure 8 (or in [52][53][54]) are unphysical or necessarily incorrect.Indeed, once the unsteady flame modes have been triggered, the boundary conditions are sufficient to capture the dynamics, which evolve convectively.It is only critical values of U at which the dynamics transitions from steady to unsteady which are uncertain.It is not clear whether this is due to the construction of the boundary conditions, or the finite domain length, or both.Further work is needed in this area, and is ongoing within the author's research group.

C. Inert isotropic decaying turbulence
I next demonstrate the ability of the method to handle inert turbulence via the benchmark problem of the threedimensional Taylor-Green vortex.This problem has been widely used as a benchmark for high-order numerical methods [8,[58][59][60], and reference data associated with [59] is available online [61].The computational domain is cubic, periodic in all dimensions, with side length 2π.The domain is discretised with a uniform resolution, with unstructured nodes A 12 (a section of the node distribution in A 12 is shown in the left of Figure 1), and a regular distribution along the x 3 axis.
The initial conditions are those given in [59] for compressible flows (sinusoidal and divergence free, but details omitted here for brevity), with an initial Reynolds number of Re = 1600, and a Mach number of M a = 0.1.For this inert flow, I remove the temperature dependence of the transport properties by setting r T = 0, set c p to a constant independent of temperature, and P r = 0.71.Denoting the resolution by the number of nodes (i.e.N = (2π/s) 3 ), I systematically vary the resolution from 64 3 to 512 3 , and the order of the scheme from m = 4 to m = 10. Figure 13 shows the evolution of volume averaged enstrophy with time for this case, for this range of resolutions and orders.Figure 14 shows the time evolution of kinetic energy with m = 0.8, for four resolutions (black lines).In both figures, the red lines correspond to the reference data accompanying [59].
My results approach the reference solution with both increasing resolution, and increasing polynomial consistency m.For m = 8 and m = 10, the maximum errors in enstrophy are small (below 3.6% and 1.8% respectively) at a resolution of 512 3 .For m = 4, the error in the enstrophy deviates significantly from the reference solution even at the finest resolution.This is because the filter used for de-aliasing is only 4 th order, and the scales on which it acts are too close to the length scales of the flow.In terms of the kinetic energy evolution with m = 8 (Figure 14), the match is closer still, with only the coarsest resolution showing significant deviation from the reference data.The size of the computational stencil (and hence the cost) increases with m, and this has guided the choice in this work to generally opt for m = 8 as a compromise between accuracy and efficiency.I note that a promising approach is to use a m = 4 or 6 for the spatial discretisation, but combined with m = 8 for the filter.This will provide the cost benefit of reduced stencil sizes for the bulk of the calculation, whilst the higher order filter will not be overly dissipative at the physical scales of the solution.Furthermore, such an approach could easily be implemented adaptively, with m = 8 in regions of interest (where fine flow structures or flames exist), and lower (m = 6 or m = 4) elsewhere.
These results demonstrate that the LABFM discretisation possesses sufficient accuracy for simulations of turbulent flows.The discrepancies between the present method and the reference solution are comparable to those of other leading methods for DNS, including the filtered finite difference schemes in [62] and the high-order flux reconstruction scheme in [8].The discrepency compared with the reference solution is expected.The reference solution uses a pseudo-spectral method de-aliased by Fourier smoothing (see references within [58]), which exhibits negligible numerical dissipation up to about 80% of the Nyquist wavenumber of the discretisation, whilst the LABFM operators have dissipation profiles similar to high-order finite differences, and for m = 8, have negligible numerical dissipation only up to about 50% of the Nyquist wavenumber, as shown in [31].Furthermore, given the highly three-dimensional nature of this flow, they provide a validation of the coupling between the LABFM discretisation in A 12 , and high-order finite differences in x 3 .I have performed additional simulations with the initial velocity field rotated by 90 degrees, and obtain almost identical results in terms of the time evolution of enstrophy and kinetic energy, with the maximum discrepancies being 0.35% and 0.057% respectively.

D. Flame-turbulence interactions
With the ability of the framework to simulate both flames and inert turbulence demonstrated, I briefly turn my attention to flame-turbulence interactions.I simulate the evolution of a statistically planar 3D flame under decaying isotropic turbulence.This problem has been studied extensively (see e.g.[8,[63][64][65][66][67][68][69]) although often with slight differences in the detail, for example in the choices of turbulence intensity relative to laminar flame speed u ′ /S L , the turbulence lengthscale relative to flame thickness l/δ th , or the method by which the initial turbulence field is generated.
The computational domain has size L × L/2 × L/2, with L = 10mm.The mean flow is in the x 1 direction, in which partially non-reflecting inflow and outflow conditions are applied, whilst the boundaries in the transverse directions are periodic.The resolution is uniform in the transverse directions, and varies from s = L/125 = 80µm at the inflow and outflow, to s = L/500 = 20µm at the centre of the domain where x 1 = L/2.The resolution in x 3 is uniform, FIG. 13.The variation of volume averaged enstrophy with dimensionless time for the three-dimensional Taylor-Green vortices at Re = 1600, for different orders of accuracy (top left -4 th order, top right -6 th order, bottom left -8 th order, bottom right -10 th order).In all plots the red line shows the reference data accompanying [59].
and set at the number-weighted mean value of s, giving s 3 = L/292 = 34.2µm.In the interests of computational expedience, and to match existing simulations in the literature, the chemistry is represented by a single-step Arrhenius model tuned to match the laminar flame speed of a stoichiometric methane-air mixture, following [8].Both reactions and products are assumed to have unity Lewis number, I set P r = 0.7, and take the temperature dependent exponent for the viscosity to be r T = 0.7.The flow is initialised from a one-dimensional laminar flame solution superimposed on a turbulent velocity field, with the unburnt gases having temperature T u = 300K, and the outflow tracking pressure p out = 10 5 P a.The turbulent velocity field is obtained by applying diffusion to a field of uncorrelated, uniformly distributed random numbers, following the method of [70].This approach is well suited to the present method where nodes are not arranged on a structured grid, and as the diffusion process is fully explicit, it is easily incorporated into the numerical framework.Furthermore, it is well suited to complex geometries, and the study of flame-turbulence interactions in non-trivial geometries is a key focus of future work.There are two key limitations of the method of [70].Firstly, the resulting velocity field is not divergence free [71], which can result in non-negligible acoustic energy which must be dissipated in the early stages of the simulation.Secondly, it offers little control over the energy spectrum.By construction, the spectra of the velocity components are flat at small wavenumbers (2π/l), with sharp decay at higher wavenumbers.Despite these drawbacks, the method of [70] is still used in turbulent combustion simulations, and I follow [8] in using it to generate the initial velocity field.The turbulent velocity field is homogeneous and isotropic, and I perform simulations at three different turbulent lengthscales l/δ th ∈ [1, 2, 3] (where δ th is the (laminar) thermal Figure 15 shows isosurfaces of the reactant mass fraction (in blue), over a background slice of the vorticity magnitude (in red), for two values of turbulence length scale: a) l/δ th = 3 (left panel) and b) l/δ th = 1 (right panel).In both cases, the results correspond to time t = 2l/u ′ .For l/δ th = 3, the configuration matches that in [8], and the resulting flame is qualitatively similar, whilst for l/δ th = 1 the wrinkling of the flame is correspondingly smaller, both in typical lengthscale, and amplitude.For both cases, the vorticity magnitude downstream of the flame is suppressed at all scales due to the gas expansion across the flame, as predicted by theory [72] and previous numerical experiments [68].
The turbulent flame speed is calculated by integrating the heat release rate over the domain, and the time evolution of the ratio of turbulent to laminar flame speeds S T /S L is shown in Figure 16, alongside data from similar simulations taken from [64] (solid red line) and [65] (dashed red line).For l/δ th = 2 (solid black line) I see a similar trend between the present simulations and the results of [64], with S T /S L increasing by a factor of approximately 2 over the first two eddy turnover times.For l/δ th = 1 (dashed black line), the turbulent flame speed increases more slowly, and for l/δ th = 3 (dot-dash black line) the turbulent flame speed increases more quickly, as expected.The discrepancy between the data of [64] and [65], despite similar configurations, serves to highlight the challenge in accurate quantitative comparisons for this case.The evolution of the turbulent flame speed is strongly dependent and l/δ th = 1 (dashed black line).The solid red line corresponds to data taken from [64], and the dashed red line corresponds to data taken from [65].
on initial turbulent energy spectrum, as this has a significant bearing on how the turbulence develops over time.
Whilst [64] generated a divergence free initial velocity field in Fourier space, with a Batchelor-Townsend energy spectrum, my initial velocity field is not divergence free, and has a flat energy spectrum, although the lengthscale below which the energy rapidly drops off is similar.The results in [65], and my own numerical experiments, show some variability between individual DNS runs, so quantitative comparisons should be made on ensemble averages rather than the results of a single simulation.Nevertheless, the results shown in Figures 15 and 16, especially when taken in the context of the results in Sections IV A and IV C indicate the present numerical framework is capable of accurately simulating flame-turbulence interactions.

E. Hydrogen ignition in an array of cylinders
Finally, I present the results of a simulation designed to show the potential of the numerical method for more complex geometries.I offer no experimental data or comparison with other numerical methods, and include this case simply as an example of the geometries for which the present method is suited.Whilst the fundamentals of the numerical framework have been extensively analysed in previous papers [30,31], and the ability of the method to simulate flames, turbulence, and flame turbulence interactions has been assessed above, this test introduces a step change in the level of complexity, and further assessment and validation will be conducted in an in-depth study in future.
The computational domain has size 13.53D×2.5D×D,with D = 1mm.Non-reflecting inflow and outflow boundary conditions are imposed in x 1 , with periodic conditions in x 2 and x 3 .An array of cylinders each with diameter D is located with the centre of the central cylinder a distance 2.25D downstream of the inlet.The cylinders are arranged in an isometric configuration, with spacing 1.25D.The domain contains a stoichiometric hydrogen-air mixture, with inflow velocity U = 10ms −1 , and inflow temperature T in = 300K.For the results reported herein, I use dimensionless time t ⋆ = tU/D.The outflow boundary tracks a pressure of p out = 10 5 P a.The domain is discretised with a nonuniform resolution, varying from s min = (5/6) × 10 −5 m at the cylinder surfaces, to s out = 2s min at the outflow, and s in = 6s min at the inflow.No closed form analytical expression is available to describe the spatial variation of s through the domain, but input files containing lists of node positions and resolutions are available on request.This discretisation ensures that in regions where the flame is expected to reside, the resolution is approximately the minimum sufficient for laminar flame simulations, as identified in Section IV A. i conduct both two-and three-dimensional simulations.In each case, an inert simulation is run for one flow-through time, to obtain an initial velocity field.The simulation is then restarted with a Gaussian hot spot (with peak temperature T hot = 2500K and characteristic width 0.2mm for two-dimensional simulations, and 1mm for three-dimensional simulations), located just upstream of the cylinder array, a distance 2mm upstream of the central cylinder, and centered in the x 2 and x 3 directions.The hot spot triggers ignition of the mixture, and I simulate the flame propagating through the cylinder array and

downstream.
Figure 17 shows the development of the flame kernel through the cylinder array for the two-dimensional simulation.The left panels a) to d) show isocontours of the heat release rate.The flame develops from the Gaussian hot spot and is initially drawn in to the array as can be seen in panels a) and b).A later times, the flame kernel begins to propagate through the cylinder array wake (panels c) and d) of Figure 17), at which point the flame speed increases due to the vortical structures in the flow.The right panels e) to h) show isocontours of different fields at dimensionless time t ⋆ = 1.25.In panel f) the increase in secondary and unsteady vortical structures in the downstream half of the cylinder array is apparent.In panels g) and h), which show the mass fractions Y O and Y HO2 respectively, the assymetry of the flame near the leading edge of the cylinder array is visible.This can also be seen in the isocontours of the heat release rate in panel d).At this stage the flame doesn't connect with the upper and lower cylinders in the central column.I also see in panel g) a very non-uniform distribution of Y O within the flame in the cylinder array, with significant accumulations at certain locations on the cylinder walls, likely due to the small scale vortical structures within the array.
Figure 18 shows the flame kernel at time t ⋆ = 0.6 for the three-dimensional simulation.Note, the hot spot used for ignition is larger for the three-dimensional simulation, and the flame develops faster, whilst the upstream flame-front has not yet been advected towards the cylinder array.The cold-flow velocity field is largely two-dimensional, with fluctuations in x 3 of u ′ x3 ≈ 5 × 10 −3 U .As the flame front develops, it triggers significant three-dimensional flow, especially in the vortical structures downstream of the cylinder array, and the flame kernel itself becomes highly three-dimensional in this region.At later times (not shown), the velocity fluctuations in the third dimension reduce and the flow becomes nearly two-dimensional again.
Although I offer only observations of this test (I leave detailed analysis and further study of such problems to a future paper), it highlights the type of problem to which the present method is well suited: DNS of combustion in complex geometries, discretised to a high-order, and able to capture flame-vortex-structure (or flame-turbulencestructure) interactions.With such problems beyond the the capabilities of leading high-order finite difference codes for combustion DNS, the present method offers a complementary approach which extends the state-of-the-art for high-order methods in non-trivial geometries.

V. COMPUTATIONAL PERFORMANCE AND SCALING
As discussed in Section I, DNS of turbulent combustion is computationally demanding.The requirement to resolve the Kolmogorov scales gives rise to the number of collocation points scaling with Re 9/4 [4], and the requirement to resolve the flame structure is often even more restrictive.With the present method having similar resolving power to high-order finite differences, the resolution requirements are similar to other leading combustion codes (e.g.SENGA+ [5], S3D [6] and KARFS [7].Simulations at larger Re, although possible, require signicant increases in computational resources, and with the present work focussed on demonstrating the potential of the method, I have limited my investigations to flows up to Re ≈ 2000. For the quasi-one-dimensional laminar flame tests conducted in Section IV A, a comparison with the two reference codes is not possible as for this test SENGA+ is used to perform truly one-dimensional simulations, whilst the present method is performing two-dimensional simulations, and Cantera (also truly one-dimensional) is solving for the steady state, rather than the unsteady system.However, comparisons can be made against other high-order finite difference codes.Broadly, the construction of the numerical framework is similar to other high-order finite difference codes for combustion (e.g.SENGA+), and the primary difference in computational costs arise from the difference in stencil sizes.For a finite difference scheme of order m, the computational stencil contains N = dm + 1 elements, where d is the number of dimensions.For the present framework, the computational stencils are larger: for m = 6, N ≈ 55, and for m = 8, N ≈ 75 in two dimensions (and on average, as there is spatial variation in this due to resolution gradients).For three-dimensional simulations the values of N for the present method increase by +m to include the finite difference stencil in x 3 .Hence, for a two-dimensional simulation of order m = 8, the computational stencil is approximately 4.5 times larger for the present method than an equivalent finite difference scheme.This is the cost of geometric flexibility.The benefit, however, is that variable resolution can be easily implemented, and in a more flexible and localised manner than is possible with grid-stretching techniques for finite difference schemes.I take as an example the two-dimensional bluff body simulations in Section IV B. The simulations of [53] utilised a uniform grid with approximately 10 6 collocation points.In the present work, the same geometry is discretised with approximately 2 × 10 5 collocation points, with a coarser resolution in regions where the flame is unlikely to occur (i.e.upstream of the bluff body, and away from the channel centre).For this case, the benefits of variable resolution (smaller N ) and the costs of a mesh-free scheme (larger N ) approximately balance (to an order of magnitude), with [53] reporting costs of approximately 800 CPU-hours per millisecond simulated, compared with 1132 CPU-hours per millisecond for the present work.A caveat to this comparison is that the two sets of simulations use codes of differing maturity (a new code here, vs. one subject to extensive development and optimisation in [53]), run on different HPC systems, and with different approaches for certain aspects (e.g.use of Chemkin libraries, vs. in-house libraries in the present work).Furthermore, this is just one case, and the balance will be different in different scenarios: there will be cases where a uniform grid-based scheme is always faster, and cases where the present method performs better.There will also be cases which simply cannot be simulated by existing high-order finite-difference combustion codes, such as the final test case in this work, and these are the problems on which I intend to focus this method in future.
The memory requirements of the present method are greater than finite difference schemes (and more akin to highorder flux-reconstruction (FR) schemes), as a set of weights (or interpolant coefficients for FR schemes) must be stored for each collocation point (element/cell in FR schemes).This requirement is only prohibitive when performing scaling tests (e.g.three-dimensional Taylor-Green vortices at a resolution of N = 512 3 running on 8 processors).Although the number of elements in a computational stencil is larger than for a finite difference scheme, the stencil lengthscale h/s is comparable.Hence, the amount of information which must be passed between subdomains is comparable with an equivalent order finite difference scheme (i.e.similar sized halos).As with high-order finite difference schemes, provided the number of nodes per MPI rank is large enough (typically > 6000), the MPI communication does not present a bottleneck to simulation, with the cost of communication remaining small compared to the cost of calculation.
Figure 19 shows the strong scaling performance of the method for the Taylor-Green vortex test case from Section IV C, for both two-dimensional (circles) simulations at a resolution of 2048 2 , and three-dimensional simulations at the finest resolution of 512 3 .In both cases, simulations were performed with a single thread per MPI rank, and so this is a test of the scaling performance of the distributed memory parallelism.In both cases there is almost perect scaling from the smallest number of cores (4 and 8 in two and three dimensions respectively) to the largest available on the system (1024 and 1000, for two and three dimensions respectively).For the two-dimensional case, the parallel efficiency from 4 to 1024 cores is 96.49%.Deviations above ideal scaling for the three-dimensional case at low numbers of cores (< 128) are due to slight differences in configuration (e.g.not utilising every core on a compute node, to give each MPI rank enough memory.).For the three-dimensional simulations, the scaling efficiency is ideal.For the three-dimensional simulations run on 1000 cores, there are approximately 1.342 × 10 5 collocation points per MPI rank, and the cost of MPI communication is small relative to the cost of calculations.

VI. CONCLUSIONS
Whilst direct numerical simulations (DNS) of combustion have long been restricted to simple geometries, the development of methods which provide the accuracy, speed and flexibility to enable the computation of flows in non-trivial geometries is an active area of research for the community.In this work I have presented a framework which allows for high-order discretisations of complex geometries.The discretisation is mesh-free in two-dimensions, utilising the Local Anistropic Basis Function method (LABFM) to evaluate spatial derivatives, combined with highorder finite differences for the third dimension.The method is validated against several two-and three-dimensional test cases, yielding excellent agreement with a leading high-order finite difference code for laminar flame simulations.The ability of the method to handle turbulent flows is demonstrated with convergence towards reference solutions for simulations of a three-dimensional Taylor-Green vortex, whilst the invariance of the numerical results to rotations of the discretisation frame show the compatibility of the coupling of a mesh-free method in two dimensions to high-order finite differences in the third.My results show that increasing the polynomial consistency (and hence the resolving power) allows equivalent levels of accuracy to be obtained at coarser resolutions, and I find that using the method with consistency of order m = 8 provides sufficient accuracy for resolving both flames and turbulent structures at reasonable resolutions.For flame-turbulence interactions, results are comparable with those in the literature, whilst close agreement with published simulations is obtained for bluff body stabilised flames, demonstrating the potential of the method for DNS of combustion in non-trivial geometries.
This new mesh-free method has sufficient resolving power for combustion DNS, and offers an alternative route to high-order combustion DNS in non-trivial geometries.Whilst the present method can handle spatially varying resolutions, the discretisation is constant.A key advance which is under development is a framework for resolution adaptivity, such that the discretisation is dynamically refined in regions where the length-scales of flow features are small.Additionally, in depth studies are planned for bluff body stablised flames, and the effects of bluff body geometry on flame dynamics.
In the present work, the geometric complexity of the domain is limited to two dimensions, with homogeneity a requirement for the third due to the use of finite differences in that dimension.I note that there is nothing precluding the use of LABFM for all three dimensions, although the size of the computational stencil will be significant, increasing costs.Developments are under way create optimised stencils, which still yield high-order interpolants, but at reduced costs.This will render a three-dimensional implementation of LABFM more feasible, allowing truly complex geometries to be simulated.

FIG. 1 .
FIG. 1. Left panel: A typical node distribution in the x1 − x2 plane A12 around an aerofoil.The color represents the local resolution lengthscale s normalised by the chord length C. Right panel: A typical stencil for a three-dimensional simulation.The nodes shown are those in the stencil of the central node in blue.The unstructured nodes in yellow lie in the same plane A12 as the central node.The nodes in red are those in the finite difference stencil for the x3 direction.

FIG. 2 .
FIG. 2. A schematic of the discretisation near boundaries, for a) walls, and b) inflows and outflows.The solid black lines indicate a boundary, and the small circles indicate computational nodes.The straight dashed black lines indicate the boundary normal direction.For walls, the fictitious nodes required to complete the finite difference stencil are shown in red.Nearboundary nodes, at which I impose the limit m = 6 are shown in light blue.Internal nodes are shown with a white fill.For wall boundaries, the red dashed lines indicate the regions at which different orders are imposed, and the dashed black circle shows the computational stencil used for the fictitious node interpolation.

4 .
For each slice (a) Re-arrange the nodes in order of increasing x 2 .(b)Split the slice into N P 2 blocks each containing the same number of nodes.

FIG. 3 .
FIG. 3. The node distribution for the planar laminar hydrogen-air flame simulations.Dimensionless resolution s/L (bottom), temperature profile (middle), and a magnified view of H2O2 mass fraction (top).

12 FIG. 5 .
FIG. 5.The spatial profiles of temperature (left) and OH mass fraction (right) through a freely propagating stoichiometric hydrogen-air flame, for different orders of LABFM: m = 4 (dash-dot), m = 6 (dotted), m = 8 (dashed), and m = 10 (solid).For all values of m, the simulation was run with a uniform resolution of s = 10µm.Note YOH has been scaled by a factor of 10 3 .

FIG. 7 .
FIG. 7. Spatial discretisation for the two dimensional square bluff body configuration.Left, upper half: the resolution, normalised by the domain lengthscale L. Left, lower half: the domain decomposition with each coloured block representing a different MPI rank.Right: a close view of the bluff body, showing isocontours of the heat release rate (the colour level is a maximum at 10 10 W/m 3 ).The fields shown in the top half are interpolated onto a background mesh for visualisation purposes.

Figure 8
shows the flame front obtained with the present method for several Reynolds numbers, illustrating the different dynamics.Panels a) and b) show a steady flame at Re = 810 for the circular (a), and square (b) bluff body.

FIG. 8 .FIG. 9 .
FIG. 8. Isocontours of the heat release rate for a hydrogen-air flame anchored to a bluff body.a) Circular bluff body with Re = 810 and a steady flame; b)-d) a square bluff body with three flame regimes: b) steady, at Re = 810, c) symmetric vortex shedding at Re = 2160, and d) asymmetric vortex shedding at Re = 2160.The colour level is a maximum at 2 × 10 10 W/m 3 ).

FIG. 10 .
FIG.10.Isocontours of pressure (relative to the outflow pressure) around a flame anchored to the square bluff body in the asymmetric vortex shedding regime.The contour at which YH 2 O is half its equilibrium value is shown in white.

FIG. 11 .
FIG. 11.Isocontours of the heat release rate for three domain lengths a) Lx 1 , b) 1.5Lx 1 , and c) 2Lx 1 .In all cases, the simulations are started from the same initial conditions, and run until t = 3.75ms.

32 FIG. 12 .
FIG.12.Time evolution of the volume averaged heat release rate for simulations of the circular bluff body flame, for three domain lengths.

FIG. 14 .
FIG. 14.The variation of volume averaged kinetic with dimensionless time for the three-dimensional Taylor-Green vortices at Re = 1600, for m = 8 and various resolutions.The red line shows the reference data accompanying [59].

FIG. 19 .
FIG.19.Strong scaling of the code for the benchmark problem of decaying Taylor-Green vortices, in two-(circles) and three-(crosses) dimensions.The dashed line shows the slope for perfect scaling.Note that times have been re-scaled to arbitrary units.

TABLE I .
Schemes for derivatives near wall boundaries.IFD indicates the interpolated finite difference scheme described in detail in Appendix A.Node Distance from boundary First derivatives Second derivatives Filtering operator a0 shows Left panel: the spatial profile of OH mass fraction through a freely propagating stoichiometric hydrogen-air flame, for different resolutions.Right panel: the spatial profile of pressure (relative to the outlet pressure) across the flame, for the same set of resolutions, with comparison to the results from SENGA+ (blue line).For all resolutions, the simulation was run with m = 8.Note YOH has been scaled by a factor of 10 3 .