Extension of the Spectral Difference Method to Premixed Laminar and Turbulent Combustion

A Spectral Difference (SD) algorithm on tensor-product elements which solves the reacting compressible Navier–Stokes equations is presented. The classical SD algorithm is shown to be unstable for test cases involving a multi-species gas whose thermodynamic properties depend on temperature and species mass fractions. It is observed that the extrapolation of conservative variables from the solution points to the flux points is the origin of the instability produced by the classic SD algorithm. An alternative scheme, where primitive variables are computed at the solution points and then extrapolated to the flux points, is shown to make the SD method stable for such cases. In addition, a new, more stable methodology for computing the diffusive fluxes at an interface is detailed. To allow the simulation of combustion, characteristic and wall boundary conditions for multi-species flows in the SD framework are introduced, and the thickened flame turbulent combustion model is adapted to the SD context. Validation test cases from one-dimensional and two-dimensional laminar flames to a three-dimensional turbulent flame are performed showing excellent agreement with each of the reference solutions. To the authors’ knowledge, this is the first time that a 3D turbulent flame is simulated using the SD method. The interest of employing a high-order method, such as the SD, in combustion is demonstrated: high values of the polynomial degree improve the quality of the results which advocates for the use of p-refinement when simulating reacting flows.


Introduction
Combustion phenomena are quite complex to understand and the use of numerical simulations, in concert with experimental studies, is essential for the design of combustion engines. These simulations must be able to accurately reproduce the multi-scale physical phenomena inside a combustion chamber, including flow structures, complex chemistry, and strong variations of thermodynamic properties. During the last decades, Large Eddy Simulations (LES) have appeared to be a good compromise between Reynolds Average Navier-Stokes (RANS) equations and Direct Numerical Simulations (DNS) in the combustion field. However, LES require accurate spatial discretizations with low dissipation and low dispersion errors.
With finite difference (FD) or finite volume (FV) techniques, increasing spatial accuracy is done by increasing the size of the stencil of the scheme, which by construction becomes very expensive when using unstructured meshes that are needed for complex geometries. In contrast, recently developed high-order discontinuous methods give access to high-order accuracy on unstructured meshes with efficient parallelization. All of these methods follow the same principle: they define a high-order representation (typically a polynomial of degree p) of the variables within each mesh element using a high-order interpolation procedure. A major advantage is to exploit the possibility of increasing the spatial resolution by either locally increasing the polynomial degree p, called p-refinement, or by doing local mesh refinement, called h-refinement. The first one reduces dissipation and dispersion errors in smooth flow regions whereas the second one is able to isolate regions with geometrical and physical discontinuities (Chapelier et al. 2014;de la Llave et al. 2018).
The most well-known high-order discontinuous approach is the Discontinous Galerkin (DG) method developed by Reed and Hill (1973) for the neutron transport equation. Almost 20 years later, the DG approach was applied to conservation laws and more specifically to the Navier-Stokes equations (NSE) (Cockburn et al. 1990; Shu 1997, 2001). Since then, it has been widely employed to perform LES of multiple problems such as turbulent jets (Anghan et al. 2019), laminar to turbulent transition (Beck et al. 2014), shock waves (Renac et al. 2015) and combustion applications (Billet et al. 2014;Ihme 2014, 2017;Johnson and Kercher 2020;Du and Yang 2022). However, the computational cost of the DG approach grows rapidly when the polynomial degree is increased. The use of the weak integral form of the NSE by the DG method requires integrals to be evaluated using quadrature rules, which become very costly as the scheme order is increased (Yu et al. 2014).
Following this observation, new high-order discontinuous methods built to solve the strong form of the NSE have emerged and they are typically splitted into two families. The first one is the Flux Reconstruction (FR) technique, also called the Correction Procedure for Reconstruction (CPR), developed in Huynh (2007). In this method the solution and the flux polynomials are collocated which implies that the flux divergence is no longer a polynomial of degree p but also that the scheme is not conservative because the flux is discontinuous at element interfaces. Flux correction functions taken as polynomials of degree p + 1 are then used to tackle these issues and make the FR method conservative. It has been applied to multiple configurations of flows over the past 10 years (Gao and Wang 2009;Haga et al. 2011Haga et al. , 2015Singh and Frankel 2019).
The second family of strong discontinuous approaches is the Spectral Difference (SD) method, originally introduced from the work of Kopriva and Kolias (1996) as the staggered-grid Chebyshev multidomain method. They applied it to structured quadrilateral elements using a tensor-product formulation. Later on, the method was extended to triangular elements and conservation laws (Liu et al. 2006). Then, Wang et al. (2007) used it for the Euler equations and the extension to the NSE was done by May and Jameson (2006) for triangular grids and by Sun et al. (2007) for hexahedral elements. For tensor-product cells such as quadrilaterals or hexahedra, the standard SD method builds a polynomial of degree p for the solution vector using solution values at what are called solution points (SP) as well as a polynomial of degree p + 1 for the flux vector using flux values at another set of points called flux points (FP). This leads to a scheme order of p + 1 . The SD method is becoming increasingly mature, having been successfully used to simulate wall boundedturbulent flows , bypass transition (Pinto and Lodato 2019;Lodato et al. 2021); shocks and detonations (Gupta et al. 2018;Lodato 2019;Tofaili et al. 2021) and with moving meshes (Yu et al. 2011;Zhang et al. 2016).
Some studies have tried to compare DG, FR and SD techniques but there is no final conclusion on which high-order method should be used for a given application. The comparative study of Liang et al. (2013) showed that FR is the fastest but least accurate of the three, whereas DG is the slowest but most accurate approach. In terms of speed and accuracy, the SD method occupies an intermediate position, with the recent study of Cox et al. (2021) confirming that the SD method is more robust and accurate than FR. This is very encouraging for pursuing the development of the SD method to simulate complex flows. Fiévet et al. (2020) have recently developed non-reflecting boundary conditions that are specially adapted to SD and FR algorithms on tensor-product cells, making it possible to simulate combustion in confined media using these approaches. To the authors' knowledge, only Gupta et al. (2018) and Tofaili et al. (2021) have published combustion simulations using the SD approach, and these studies were limited to one-dimensional (1D) detonations. Consequently, there is a need to explore the ability of such methods to compute combustion applications involving multi-species gases, combustion models and stiff reactive source terms.
In this paper, the SD code JAGUAR (proJect of an Aerodynamic solver using General Unstructured grids And high ordeR schemes) (Cassagne et al. 2015;Veilleux et al. 2021Veilleux et al. , 2022Fiévet et al. 2020), jointly developed by CERFACS and ONERA, is extended to reacting flows and validated on several academic test cases of increasing difficulty. The p-refinement recently implemented in JAGUAR and tested on the convection of an isentropic vortex and then on laminar flows over a cylinder (Hirsch et al. 2021) allowed a gain in reduction of the number of grid points from 38% to more than 50% compared to a uniform p calculation to achieve the same level of error. The gain in computational time is then significant, between 25 and 50% compared to a uniform p calculation. Similar observations were made when using the DG discretization on such academic test cases (Naddei et al. 2018). Efficient p-refinement is very appealing for combustion applications where the flame zone requiring fine discretization is very localized. Moreover, it does not need remeshing compared to the current state of the art in combustion which is using a FV discretization coupled with an Adaptive Mesh Refinement (AMR) technique (Bell 2005). Therefore, it motivates the present work, focusing on the resolution of combustion problems with the SD method and, as a first step, using a uniform polynomial degree p.
The objectives of the paper are twofold. First, the development of a stable SD algorithm on hexahedral elements to simulate multi-species reacting flows has to be considered. To reach this, a modification of the usual SD method coupled with characteristic boundary conditions for a multi-species gas and a new way to evaluate the interface gradients will be presented. The use of a positivity-preserving limiter introduced in Tofaili et al. (2021) was also employed in some cases. The second objective is to demonstrate the ability of the SD method to simulate several combustion cases of increasing complexity, with a particular emphasis on investigating the impact of changing the polynomial degree.

Governing Equations
In this paper, the three-dimensional reacting compressible NSE for a multi-species gas composed of N s species are considered (Poinsot and Veynante 2005): are respectively the sum of convective and diffusive fluxes of U along x, y and z directions and S is a source term vector due to chemical reactions. They read as: In Eqs. (2-3), is the gas density, u = (u, v, w) T is the velocity vector along physical coordinates x = (x, y, z) , E is the total non-chemical energy per unit mass, Y k is the mass fraction of species k ∈ [[1, N s ]] and P is the static pressure. Furthermore, , q = q x , q y , q z T and M k = M kx , M ky , M kz T are respectively the viscous stress tensor, the energy and species diffusion flux vectors defined as (Poinsot and Veynante 2005): where I is the identity matrix, is the thermal conductivity, is the dynamic viscosity, T is the temperature, W is the molar mass of the mixture and h sk , X k and W k are respectively the mass sensible enthalpy (tabulated every 100 K from Joint Army Navy Air Force (JANAF) thermochemical tables (Stull and Prophet 1971)), the mole fraction and the molar mass of species k. In Eq. (4), V k and V c are the diffusion velocities and their associated correction velocity under the Hirschfelder and Curtiss approximation for k = 1, N s (Poinsot and Veynante 2005): with D k the diffusion coefficient of species k into the rest of the mixture, computed here assuming a constant Schmidt number Sc k for each species whose expression is written in Table 1. The dynamic viscosity is obtained using a power law for as a function of T as stated in Table 1 where T ref , ref and m are respectively a reference temperature, the dynamic viscosity at this reference temperature and the power law exponent. For the thermal conductivity, the Prandtl's number Pr of the mixture is assumed to be constant so that is computed following the expression shown in Table 1 with C p the heat capacity at constant pressure of the mixture. The net production rate of each species ̇k is computed using a classical Arrhenius's law (Poinsot and Veynante 2005) and the heat release rate ̇T is deduced from Eq. (6): where Δh 0 f ,k is the formation enthalpy of species k also obtained using JANAF thermochemical tables (Stull and Prophet 1971). Finally, the equations are closed assuming that the mixture and each species behave as a thermally perfect gas where the static pressure is the sum of partial pressures of each species (Dalton's law): with R = 8.314 J.mol −1 .K −1 the ideal gas constant. In Eq. (1), one equation is redundant since for a multi-species gas: Consequently, the mass conservation is not solved by the numerical scheme and the density is recomputed using Eq. (8) with the transported Y k values.

The Spectral Difference Method
In this section, the SD discretization process is introduced for hexahedral elements only.

Isoparametric Transformation for Hexahedral Elements
Let's consider a computational domain Ω divided into N e non-overlapping hexahedral elements inside which Eq. (1) is to be solved. Each element Ω e of Ω is transformed into a standard hexahedron H e = {( , , ), 0 ≤ , , ≤ 1} following a so-called isoparametric transformation characterized by a Jacobian matrix J along with its inverse (assuming a non-singular transformation) representing the reverse transformation: where the components of J −1 are the grid metrics whose expressions can be found in Sun et al. (2007). Equation (1) is solved in the reference domain H e so that it will be written in this domain to get (Sun et al. 2007): with where |J| is the determinant of J.

General Principle
The SD discretization will be described for an order of accuracy of p + 1 inside H e . First, two sets of points are defined inside H e : a set of solution points (SP) and a set of flux points (FP). The solution is the vector of conservative variables whereas fluxes are the convective and diffusive terms in the NSE. Building a polynomial degree p (respectively p + 1 ) for the solution (respectively flux) along each direction , and , the solution (respectively flux) has to be known at p + 1 (respectively p + 2 ) points in each of these directions. In this work, the p + 1 SP in one dimension are the Gauss-Chebyshev points of the first kind defined in Sun et al. (2007) and the p + 2 FP are taken as the Gauss-Legendre quadrature points for the p interior FP and the two remaining FP are located at the element boundaries 0 and 1 (Sun et al. 2007). Solution and flux polynomials are then built using Lagrange interpolation bases respectively at SP and FP. Fluxes at FP are obtained from the solution polynomial built at SP and evaluated at FP. These fluxes are continuous inside each standard element but discontinuous at their interfaces. At interface FP, a numerical flux is used for both convective and diffusive fluxes to ensure conservativity across the whole domain. For the convective fluxes, a Riemann solver is used to obtain a unique flux value that replaces the fluxes at each interface FP. The Riemann solver employed during this work is the Harten Lax and van Leer Contact (HLLC) scheme (Batten et al. 1997). For the diffusive fluxes at interface FP, two different options are used: 1. A centered scheme as in Sun et al. (2007) where viscous fluxes at the interface FP are computed using the arithmetic averages of the left and right solutions and gradients at this interface FP. 2. A new methodology developed in this work, called the SDLIFT approach detailed in Sect. 3.4, where correction functions as in the FR framework are introduced. It allows a reduction of stencil, compared to the centered scheme, and shows a gain in stability illustrated in Appendix 4.
Once internal and interface fluxes have been computed at all FP, flux derivatives are evaluated at SP using derivatives of Lagrange polynomials built at FP. Finally, equations can be marched in time using any explicit temporal scheme at each SP inside all elements. In this work, either the three-stage third-order in time Strong Stability Preserving (SSP) Runge-Kutta (RK) scheme of Gottlieb and Shu (1998) or the five-stage fourth-order in time SSP-RK scheme of Spiteri and Ruuth (2002) are employed.

Extrapolation of Primitive Variables from SP to FP
As described in Sect. 3.2, the usual discretization process in SD considers an extrapolation of conservative variables from SP to FP. This approach will be referred as CONS. Nevertheless, another approach can be considered which is to compute T, u , P and Y k at SP from conservative variables at SP and then extrapolate T, u , P and Y k at FP. This approach can be referred as TUPY. Both methodologies can be summed up as: where Q = T, u, P, Y k T is the vector of primitive variables. It was found that the CONS approach is unstable, whereas the TUPY approach is not, when considering a multi-species gas while using the SD method. To highlight this non-equivalence, a 1D domain with air ( Y O 2 = 0.23 and Y N 2 = 0.77 ) on the left side at T = 300 K and the same air on the right side at T = 2010 K is considered. Initial pressure and velocity are set constant respectively at P = 101325 Pa and u = 0.2815 m s −1 . This situation is illustrated on Fig. 1 where values of T, u and P at SP are shown. Characteristic boundary conditions for a subsonic inlet and for a subsonic outlet are used. The way to impose these boundary conditions within the SD context is described in Sect. 3.6. Then, the simulation is run with diffusion fluxes and combustion source terms set to zero to only see convective fluxes effects. This is the situation of a contact discontinuity with only the temperature (and so the density) that is changing in the transition zone between cold and hot air. In this case, pressure and velocity must remain constant during the computation. Unfortunately, this is not the case for the CONS approach where wiggles appear on the pressure (and consequently on velocity too) in the temperature transition zone as depicted in Fig. 2. On the contrary, it can be seen that these wiggles are not present for the TUPY approach. It is also shown on the same figure that there are no wiggles as well when the CONS approach is used for a mono-species gas with constant thermodynamic properties in the same situation (so with a temperature and density gradient).
To better understand what is happening, values of T, u and P can be represented at the initialization (when no time iterations were done) on ten points equally spaced per mesh element, called output points (OP), for both the mono-species case and the multi-species case. This is done to see values of T, u and P on points other than SP from which their continuous polynomials are built. Figure 3 shows these values at initialization, where it can be seen that the pressure field in the transition zone between cold and hot air shows wiggles when the CONS approach is used. These instabilities are not present for the other two cases. Therefore, for the multi-species case, it seems that building polynomials of conservative variables, evaluate them on another points than SP and compute P at theses points create pressure oscillations. On the other hand, building polynomials of P and T at SP and evaluating them on points other than SP does not create pressure oscillations. Thus, the TUPY method is retained in this work. These differences of results between CONS and TUPY approaches stem from the fact that for a multi-species gas, pressure and total energy are not connected by a linear relation, which has been already studied in the literature (Karni 1994;Abgrall 1996;Abgrall and Karni 2001).
Remark: Another approach where is extrapolated at FP instead of T, called UPY was also tested. It was also much more stable than CONS approach and results were similar to the ones of TUPY.

The SDLIFT Formulation
The centered scheme proposed in Sun et al. (2007) has a five elements stencil for the gradient evaluation at interface FP (Huynh 2009). It was found to be very unstable in the case of shear layers where the mesh discretization was not fine enough as observed in the case presented in Sect. 5.2. Although other approaches can be found in the SD literature, mostly  (Moreira et al. 2016), the local-DG like algorithm (Sun et al. 2007) or a Bassi and Rebay 2 (BR2) scheme adapted to SD (Van den Abeele et al. 2009), they did not give satisfying results compared to the centered scheme. Consequently, another methodology was developed in this work using principles from FR discretization. It is illustrated here in the 1D case. The situation considered is represented in Fig. 4 where an element H e is shown with its two interfaces between respectively H e−1 on its left and H e+1 on its right. The objective is to explain how common states and common gradients are computed on these interfaces using a new approach. Starting from the solution polynomial built at SP with conservative variables in the physical domain, a (p − 1) -degree polynomial can be built inside element H e : with U e j =Û e j ∕|J| e j and l SP j the j-th 1D polynomial of Lagrange basis of degree p built at SP j using all the other SP s with s ≠ j . If U e−1∕2 and U e+1∕2 are the common state values (given by the chosen diffusion scheme as in Eq. (17)), then U e h ∕ can be corrected into a p-degree polynomial as in the FR approach: where g L e (respectively g R e ) is the left (respectively right) correction function at the right (respectively left) interface of H e . The superscripts L and R thus refer to the sides relative to a given interface and not to the sides of the element. The adopted convention then differs to the one used in Huynh's original papers (Huynh 2007(Huynh , 2009). For instance, g L e has superscript L because it is used for the left side of the right interface in e + 1∕2 and g R e is defined equivalently for the right side of the left interface in e − 1∕2 . These correction functions are polynomials of degree p + 1 that approximate zero in some sense (Huynh 2007(Huynh , 2009 and have the following constraints: Thanks to these correction functions, Ũ e h ∕ is of degree p but is discontinuous across elements. Consequently, a new correction step is done by introducing common gradient values at element interfaces (also given by the chosen diffusion scheme as in Eq. (16)), and Ũ e h ∕ e+1∕2 , to have: Equation (15) defines a (p + 1)-degree polynomial for the gradient of U which is continuous across the computational domain and is evaluated at all FP to get Finally, values of ( U∕ ) C i are multiplied by x at each FP i to get the gradient of U along x direction. This new methodology is named the SDLIFT approach in this work since it uses the SD method coupled with lifting/correction functions as in the FR formalism. It offers the possibility to make different kinds of choices for the common states and gradients at element interfaces. Each of these choices is considered as a numerical diffusion scheme since it is a user choice on the way to set these common states and gradients. In this work, the I-continuous (IC) diffusion scheme, originally introduced in Huynh (2009) in the FR context, was considered. It computes the common state value at an element interface, for instance between H e and H e+1 , such that: The same calculation holds for the interface between H e−1 and H e . Once the common state has been computed with Eq. (17), a common gradient is computed based on Eq. (16). From Eqs. (13) and (16), it can be seen that the IC scheme has a three element stencil, in contrast to the five element stencil of the centered scheme. The correction function that is employed in this work is: where R R,p+1 is the Right Radau polynomial of degree p + 1 on ∈ [0, 1] , whose definition can be found in Huynh (2007Huynh ( , 2009 and g L e is given by: g L e ( ) = g R e (1 − ) . The reasoning presented here was done with the conservative variables U but it can be applied in the same way to primitive variables Q = T, u, P, Y k T .

Positivity-Preserving Limiter for Mass Fractions
The condition 0 ≤ Y k ≤ 1 is ensured by a limiting procedure that was first introduced by Zhang and Shu (2010) for the compressible Euler equations, and then extended to the SD framework for shock-capturing by Lodato (2019) and for multi-species Euler equations by Tofaili et al. (2021). In the present work, the enforcement of this condition differs slightly from how it is described in the current literature: the limitation on species mass fraction is almost identical, whereas it was found that limiting the 1 3 temperature field, and not the pressure field as it is commonly done (Zhang and Shu 2010), offers a gain in stability. Without loss of generality, the process is described for a given element Ω e which was flagged for limiting process. First, the spatial average of conservative variables over SP in Ω e is calculated denoted as U e . Then, e k (x) is replaced by: with: where e k,min = min SP∕FP k for k = 1, N s . Density at each SP is then recomputed using Eq. (8) to get ̃ e (x) . Then, the following vector in Ω e is defined: in order to solve the following equation, for a variable t ∈ [0, 1] , at each SP or FP noted where: < T , avoiding temperature to go below T , otherwise t = 1 . Moreover, Eq. (22) is a second-order algebraic equation in t which can be solved either analytically or using an iterative procedure such as Ridders' iterative algorithm (Ridders 1979) as suggested in Lodato (2019). Then, new corrected conservative variables at SP in the flagged element are obtained: where 2 = min t . This methodology is applied at each RK stage. The present limiter differs from those found within the SD literature (Lodato 2019;Tofaili et al. 2021), since: • There is no limiting procedure applied at FP to modify conservative variables at these points. The values at FP are only used to spot the minimum values of k and T in an element and check if they are respectively above 0 and T . Only conservative variables at SP are modified using Eq. (23) but the 2 value employed can come from a t in Eq. (22) solved at a FP. • The numerator in Eq. (20) is not set to e k − Y , with Y ≈ 10 −13 , since it is possible to have Y k = 0 for some species and then e k is really zero.

NSCBC in the Multi-species Case for the SD Method
Characteristic boundary conditions are widely employed to simulate combustion applications, but their use with the SD method was only recently formulated by Fiévet et al. (2020). In their work, a methodology to employ Navier-Stokes characteristic boundary conditions (NSCBC) for an SD discretization with a mono-species and non-reacting gas was developed. The present work extends the approach of Fiévet et al. to multi-species gases. Chemical reactions are not taken into account here, since the flame is not positioned close to a boundary.

Mathematical Formulation
Without loss of generality a boundary located at a -normal face at = 1 as in Fig. 5 is considered. The idea of the NSCBC treatment is to impose a value of Ê ∕ at the boundary FP at = 1 according to the characteristic waves crossing this boundary. The determination of this flux derivative at the FP in = 1 starts with the diagonalization of Ê ∕ in Eq. (10) which gives the following wave equation (Fiévet et al. 2020): Here, W , N and S are respectively the vector of characteristic variables, the normal and tangential strengths of the characteristics given by: with A c ( ) and A d ( ) accounting for the mesh non-orthogonality whose expressions can be found in Fiévet et al. (2020). In Eqs. (25-27), P U is the transformation matrix from conservative to characteristic variables usually expressed as: where matrices P Q and Q∕ U have been extended (from the work of Fiévet et al. (2020) to the multi-species case in this work and are given respectively in Eq. (41) and Eq. (40). (24) Illustration of the NSCBC treatment at a boundary FP where to impose a value of Ê ∕ . The polynomial degree is set to p = 3 with five FP along direction Vector W is composed of 3 + N s entropy waves W 1 , W 2 , W 3 , W 5+k , propagating at u n = u.n u ( n u is defined in Appendix 1.2), and of two acoustic waves W + and W − propagating respectively at u n + c and u n − c . In Eqs. (28-40), Q = , u, v, w, P, Y k T should not to be confused with the set of primitive variables used for the TUPY approach described in Sect. 3.3. Consequently, the main steps to compute characteristic boundary conditions in generalized coordinates for the boundary FP in = 1 in Fig. 5 are as follows: 1. Evaluate initial guesses for N and S at NSCBC FP in = 1 using Eq. (26)

Types of NSCBC
The objective of this paragraph is to briefly explain how some values of N are modified in order to impose a given NSCBC. Two types of NSCBC will be used in this work: a subsonic outflow imposing a constant pressure and a subsonic inflow imposing velocity, temperature and species mass fractions. For the subsonic outflow, to have a numerically well-behaved boundary condition (Poinsot and Lele 1992), only N − has to be modified to get (Fiévet et al. 2020): where K P s −1 is the pressure relaxation rate, P t and P are respectively the target and current pressure at the boundary and is a relaxation parameter usually taken as the averaged bulk Mach number over the whole boundary (Fosso et al. 2012). At this point, only convective fluxes have been treated. According to Poinsot and Lele (1992), for a subsonic outflow, diffusive fluxes have to satisfy (in the case of a -normal boundary): where t u 1 and t u 2 are two unit tangential vectors in the exit plane of unit normal n u . Concerning the subsonic inflow imposing u , T and Y k : N 1 , N 2 , N 3 , N + and N 5+k have to be modified. It is done by solving the system given by Eqs. (44-48) and Eq. (49) derived in Appendix 2. In the case of time-constant target values at the inlet boundary, a relaxation procedure is applied as in the subsonic outflow condition, so that time derivatives in Eqs. (44-48) and in Eq. (49) are replaced by: where K u , K T and K Y k are respectively the relaxation rates for velocity components, temperature and species mass fractions. Note that here no condition on diffusive fluxes is imposed as stated in Poinsot and Lele (1992) and Sutherland and Kennedy (2003).

Thickened Flame Model for the SD Method
For LES of turbulent premixed flames, the spatial resolution is not sufficient to resolve the flame front which entails the use of a turbulent combustion model. Several approaches may be found in the literature such as the G-equation (Kerstein et al. 1988), the use of filtered progress variables (Boger et al. 1998) or the thickened flame model for LES (TFLES) . Within the present work, the TFLES model is chosen. Artificially thickening the flame front makes it resolvable on a relatively coarse LES mesh, and a subgrid efficiency function is applied to describe subgrid-scale flame turbulence interactions. The aim of this section is to briefly present how the TF model can be implemented within the SD context.

General Principle
The TFLES model is based on a modification of the local flame structure which is done in two consecutive steps: 1. The flame is thickened by a factor F .
2. An efficiency function E is applied to increase transport and kinetics in the thickened flame region to compensate the loss in wrinkling due to the thickening step.
In practice, the use of this model consists in multiplying thermal and species diffusion coefficients, and D k , by EF and the species net production rates ̇k by E∕F . It is usually done only in regions where a flame sensor S is activated otherwise E and F tend to one. In this work: • S is computed following the procedure described in Jaravel (2016) which is based on the transport of a fictive species. • E is obtained using Charlette's model (Charlette et al. 2002) modified by Wang et al. (2011) under its saturated version as advised in Veynante and Moureau (2015) and with a constant parameter set to 0.4 in the simulations considered here.

Practical Application in the SD Context
For an SD discretization, and D k are evaluated at FP, because they are needed for flux computation, whereas ̇k are evaluated at SP. Thus, there will be values of E and F at FP noted E FP and F FP and also values of E and F at SP noted E SP and F SP . The main steps for implementing the TFLES model when using the SD method are summed up in Algorithm 1.

Results and Discussion
To validate the implementation of reacting flow models in JAGUAR, simulations from 1D to 3D flames have been performed. For some simulations, JAGUAR results are compared with results obtained with the reference solver AVBP (Schonfeld and Rudgyard 1999;Gourdain et al. 2009) developed by CERFACS. AVBP solves exactly the same reacting NSE for a multi-species gas with the same transport model. All AVBP simulations were carried out using the TTGC scheme (Colin and Rudgyard 2000) for convective fluxes without artificial viscosity for a fair comparison with JAGUAR. For diffusive fluxes, AVBP uses a finite element scheme of order 2 (Colin and Rudgyard 2000). 1D cases are also compared to the reference chemistry code CANTERA (Goodwin et al. 2017). In this section, only a 1D premixed flame and a 3D turbulent premixed flame will be presented whereas a 2D laminar burner has also been considered whose results are gathered in Appendix 3 for completeness.

One-Dimensional Flame Using a Two-Reaction Scheme
The objectives of this test case are to validate in JAGUAR: • The combustion source terms and species transport computations with a reduced tworeaction scheme. • The implementation of JANAF thermochemical tables. • The inlet and outlet NSCBC on a 1D multi-species case.

Presentation of the Case and Numerical Setup
A 1D methane-air premixed flame is considered. The chemical scheme is the two-reactions CH4/Air-2S-BFER developed by Franzelli et al. (2010). The characteristics of this flame are given in Table 2 where is the global equivalence ratio, T f is the fresh gases temperature, T b and S 0 L are respectively the burnt gases temperature and the laminar flame speed estimated by CANTERA and 0 L is the laminar flame thickness obtained from the temperature profile of the CANTERA solution. The computational domain is a 1D segment of length L x = 0.02 m discretized with N e elements (whose values can be found in Table 3). The left boundary condition is a NSCBC subsonic inflow imposing u in = S 0 L , T in = T f and Y in k such that in = 0.8 at the inlet for the 6 species involved in the CH4/Air-2S-BFER scheme. The right boundary is a NSCBC subsonic outflow imposing a static pressure P out = 101325 Pa. To allow a fair comparison between JAGUAR and AVBP, all calculations must be done with a similar number of points in the 1D domain. These points are called degrees of freedom (DOF): they correspond to cell nodes for AVBP and to SP for JAGUAR. In space dimension d, a standard hexahedral element where Û varies as a p-degree polynomial contains N SP = (p + 1) d SP so that if all these elements have the same degree p, the number of DOF inside the computational domain (total number of SP) is: Here it was chosen to keep the number of DOF around 400 in order to have a number of points in 0 L around nine which is sufficient to resolve the flame front well. While this imposes a fixed number of 400 nodes in AVBP, JAGUAR has the possibility to use different values of N e depending on the polynomial degree p to keep DOF SD close to 400 as summed up in Table 3. In this test case for all JAGUAR simulations, the centered diffusion scheme introduced in Sun et al. (2007) was employed along with the third-order in time SSP RK scheme of Gottlieb and Shu (1998) for time integration.

Results
Both JAGUAR and AVBP simulations are initialized using a CANTERA solution which contains approximately 30 points in the flame region. After a transient phase due to the transition from a constant pressure (CANTERA solution) to a compressible solution, the solution converges and final profiles of density, velocity, pressure and temperature are represented in Fig. 6 for CANTERA, AVBP and JAGUAR with p = 4 and p = 6 . Major mass fractions are plotted in Fig. 7 for CANTERA and JAGUAR at p = 4 only since JAGUAR at p = 6 and AVBP gave the same results.
All profiles are in excellent agreement. The pressure jump through the flame front is captured by both JAGUAR and AVBP while CANTERA runs at constant pressure. Theoretically, this pressure drop is given by Poinsot and Veynante (2005): where P b and P f are respectively the pressure in burnt and fresh gases and f is density of the fresh gases. Using values of Table 2, Eq. (33) leads to a pressure drop of − 0.511 Pa i.e., very close to the value of −0.512 Pa measured from JAGUAR and AVBP solutions. Finally, flame speeds, based on consumption speed, can be compared: they are summed up in Table 4 for both AVBP and JAGUAR solutions where rel is the relative error compared to a CANTERA reference value given in Table 2. Flame speeds estimated by JAGUAR are in very good agreement with the value given by CANTERA which shows its capability to well-capture the flame front.

Simulation of the Cambridge Burner
The main objective of this test case is to validate the implementation of the TFLES model, described in Sect. 4, on a fully turbulent combustion case.

Presentation of the Case
The Cambridge burner (Sweeney et al. 2011a) is an academic configuration which is composed of two concentric tubes surrounding a central bluff-body as shown on Fig. 8 where burner dimensions and injected streams are illustrated. Inner (indexed i) and outer (indexed o) inlet streams contain CH 4 /Air mixtures that are controlled independently in terms of bulk velocity U b,i∕o and equivalence ratio i∕o . The last stream is a co-flow (indexed cf) of air ( cf = 0 ) at U b,cf = 0.4 m s −1 used to isolate the flame from ambient perturbations. Finally, all injected streams work at T 0 = 298 K. This burner was widely studied experimentally by changing the different levels of CH 4 /Air stratification (variations of i and o ) and swirl providing an abundant experimental database of fifteen different configurations denoted from SwB1 to SwB15 (Sweeney et al. 2011a(Sweeney et al. , b, 2012aZhou et al. 2013). Consequently, it was also investigated numerically in LES (Nambully et al. 2014a, b;Mercier et al. 2015;Proch and Kempf 2014) and DNS (Proch et al. 2017a, b) contexts. The simulation that is considered in this work is the non-swirling premixed case SwB1 with the following operating conditions for inner and outer injected streams:

Computational Domain
The computational domain has the following characteristics: • The injector streams start at z = − 65.4 mm as illustrated on Fig. 8 and are of length L inj = 6D bb , where D bb = 12.7 mm is the diameter of the bluff-body, which extends to z = 0 . This length was shown to be sufficient to have an established turbulent flow inside the injectors. • The co-flow stream starts at z = − 20 mm. • The domain is extended 250 mm downstream of the injector exits in the z-direction so that the length of the chamber is L ch = 250 mm. The radius of the chamber is set to R ch = 140 mm.

Chemistry and Combustion Model
The chemical scheme that is employed is the two-reaction CH4/Air-2S-BFER (Franzelli et al. 2010) scheme previously detailed in Sect. 5.1. Previous numerical studies have employed more detailed mechanisms such as the Lindstedt (Lindstedt 1997) scheme used in Mercier et al. (2013Mercier et al. ( , 2015 or the GRI − 3.0 scheme (Smith et al. 2000) considered in Proch and Kempf (2014); Proch et al. (2017aProch et al. ( , 2017b. In these cases, detailed mechanisms were affordable because the LES or DNS were performed using tabulated chemistry where only the mixture fraction Z and the progress variable C are transported by the flow and the combustion source terms were read in a (Z, C) table. However, the CH4/Air-2S-BFER mechanism should be accurate enough to reproduce the correct evolution of lean mixture temperature, laminar flame speed and major species which is sufficient for validating the simulations in this work. The TFLES combustion model introduced in Sect. 4.1 is considered with a maximum thickening computed to have at least seven points in the flame front and was set to 0.4 in Charlette's efficiency function based on several preliminary tests with different values ranging from 0.3 to 0.5. It should be mentioned that no subgrid scale (SGS) turbulent models was employed compared to previous numerical studies (Mercier et al. 2013Proch and Kempf 2014) which used the -model (Nicoud et al. 2011). Although there were recent investigations in SD to see the impact of using a SGS turbulent model Lodato et al. 2021), SD approaches commonly use implicit turbulence subgrid modeling, as is done within the present work. It can be attributed to the upwind treatment of the convective fluxes at interfaces (Riemann solvers), which improves the robustness of calculations (Jameson and Lodato 2014).

Boundary Conditions
Inlet and outlet boundary conditions are imposed through the methodology described in Sect. 3.6. For the inlet of the injection tubes at z = −65.4 mm, an analytic power-law profile is imposed for the mean velocity of the form: where r 1,i∕o and r 2,i∕o are respectively the radius of the inner and outer walls of the given injected stream and U cl,i∕o is the centerline velocity in a given injected stream. Values of U cl,i∕o are computed to have the expected bulk velocity in each injector and in this work they are set to: Homogeneous and isotropic turbulence (HIT) is superimposed to these mean velocity profiles using a synthetic random Fourier method (Bechara et al. 1994;Bailly and Juve 1999) with a von Karman-Pao energy spectrum (Bechara et al. 1994) generated with 10% of turbulent intensity and an integral length scale L e = 1 mm. These parameters are the result of a sensitivity analysis performed on a reduced domain to find the values that best fit measurements directly downstream of injector streams. Species mass fractions in these injector streams are set to have the equivalence ratios defined in Eq. (34): For the inlet co-flow at z = − 20 mm, a constant axial flow at U b,cf is imposed with species mass fractions set to Y O 2 ,cf = 0.2226 and Y N 2 ,cf = 0.7774 . Outlet boundary conditions imposing P = 101325 Pa are prescribed on the sides of the chamber and at its outlet. Finally, all the walls are considered using adiabatic no-slip boundary conditions imposed following the same methodology used in the 2D burner case shown in Appendix 3.

Mesh Discretization and Other Numerical Parameters
Two different discretizations, both based on a mesh size Δ x = 0.3 mm in the flame region, were considered. The first mesh sets p = 2 (order 3) in all elements and is composed of N e = 378424 second-order hexahedral elements ending up with approximately 10.2 millions of SP. The second mesh sets p = 3 (order 4) in all elements and is composed of N e = 205,070 second-order hexahedral elements ending up with approximately 13.1 millions of SP. The objective is to see if there is an impact of the scheme order on the results even if the characteristic mesh size in the flame region is the same for both discretizations. For this configuration, the SDLIFT formulation described in Sect. 3.4 was needed to stabilize the numerical solution in shear layers compared to the centered scheme which showed high numerical instabilities highlighted in Appendix 4. The use of the positivity-preserving limiter introduced in Sect. 3.5 also helped to avoid undesired numerical behavior. Finally, time integration is carried out by the fourth-order in time SSP-RK scheme of Spiteri and Ruuth (2002).

Results Analysis
As suggested in Proch and Kempf (2014), simulations are run for a total of fifteen flowthrough times based on the inner stream velocity, with statistics recorded over the last eight flow-through times. Contour plots in the burner mid-section of temperature, heat release rate, equivalence ratio and axial velocity are shown in Fig. 9. The flame is anchored at the burner bluff-body because of the adiabatic conditions as already observed in Mercier (2015) while experiments predict a lifted flame due to the non-adiabaticity of the bluffbody surface. It can be seen that the velocity behind the bluff-body is very stable without turbulent fluctuations, attributed to re-laminarisation effects of combustion, which was already highlighted in previous works on the Cambridge burner (Zhou et al. 2013;Mercier et al. 2013;Nambully et al. 2014b;Proch et al. 2017a).
Radial profiles at different downstream locations from the burner exit of the statistics of axial and radial velocities are illustrated in Fig. 10 for the cases p = 2 and p = 3 . They are compared with Laser Doppler Anemometry (LDA) experimental data of Zhou et al. (2013) as well as numerical results from the LES of Proch and Kempf (2014). Proch and Kempf used a mesh spacing of Δ x = 0.25 mm in the flame region, which is very close to the spacing of Δ x = 0.30 mm employed here. A general good agreement is observed between the JAGUAR results, the experimental data and the LES of Proch and Kempf. The current simulation and the LES of Proch and Kempf produce small, nearly identical discrepancies with the experimental data. Looking at the mean axial velocity at z = 10 mm, it seems that JAGUAR slightly overpredicts the recirculation zone behind the bluff-body compared to both the experiment and the simulation of Proch and Kempf. This overprediction was said to be originated from the lifted character of the flame which cannot be reproduced using adiabatic boundary conditions on the bluff-body (Proch and Kempf 2014;Nambully et al. 2014b). The turbulent injection, necessary to reproduce experiments, has also probably some influence on this observation although it was employed with the best set of parameters following preliminary tests. Finally, radial profiles of temperature statistics and of mean equivalence ratio, at the same downstream distances from the burner exit, are presented in Fig. 11. The different simulations are compared with the Rayleigh-measurements provided by Sweeney et al. (2012a). Mean temperature and equivalence ratio obtained with JAGUAR are in good agreement with both experimental and numerical results. However, JAGUAR results underestimate temperature fluctuations compare to Proch and Kempf's results. This probably originates from the fact that the combustion model used in the LES of Proch's is the TFLES model but combined with a premixed flamelet generated manifolds (PFGM) tabulation method which very often ends up with more wrinkled flames compared to the classic TFLES approach (Philip et al. 2015). This effect is less visible on the p = 3 simulation, especially at z = 10 mm and z = 30 mm. It shows that, for a given mesh size in the flame region, increasing the spatial order of the scheme improves the prediction of the wrinkling even if the thickening factor is almost identical. Consequently, the limitations of the TFLES model can be compensated by using higher order schemes, making the use of high-order methods in combustion attractive.
This test case has shown that the SD method developed in this work for reacting flows is able to simulate three dimensional turbulent flames with satisfying results close to the ones obtained with classical numerical methods doing LES of combustion. It also highlights the fact that coarser mesh discretizations can be employed using high-order methods for the same results: Δ x was set to 0.3 mm in the flame region for JAGUAR with a maximum of 13.1 × 10 6 DOF using p = 3 whereas Proch and Kempf (2014) had a Δ x of 0.25 mm in this region but with more than 103.2 × 10 6 DOF. This high difference is due to the use of equidistant orthogonal cartesian grids in the LES of Proch and Kempf whereas in this work the mesh is refined only in the flame region. Other LES of the Cambridge burner, made with unstructured FV meshes refined only in the flame region, were also considered Gruhlke et al. 2021) and had a number of DOF between 6 × 10 6 to 17 × 10 6 for a minimum Δ x value of 0.25 mm. Therefore, the discretizations employed in this work are also in this range of DOF values. For future works, it will be very interesting to use local polynomial adaptation in order to keep high p values in the flame region and use low p values elsewhere. It will keep accuracy at a lower computational cost by reducing the total number of DOF in the dom ain.

Conclusion and Perspectives
The implementation of combustion source terms, the transport of a multi-species gas, the treatment of boundary conditions in the multi-species case and a new way to evaluate interface diffusion fluxes have been presented within a SD context. For multi-species flow simulations, it was shown that temperature and pressure must be computed at SP first and then interpolated on FP.
The extended SD method has been validated on 1D and 2D (see Appendix 3) laminar premixed flames showing perfect agreement with already existing solvers. Finally, simulations of the Cambridge burner (Sweeney et al. 2011a) were carried out using either p = 2 or p = 3 within all elements. It was shown that increasing p better predicts the wrinkling of the flame that is typically lost due to the thickening of the flame when using the TFLES model. Further investigations, out of the scope of such paper, must be considered to see if there is a need to take into account the value of p when computing the efficiency function E since the resolved wrinkling seems to vary with p. It is inline with the work of  who developed a SGS turbulent model which adapts the sub-grid dissipation locally to an element based on the polynomial energy spectrum within each element. Finally, to the authors' knowledge, this is the first 3D turbulent flame simulated with the SD discretization.
To conclude, for all the cases studied here, the capability of the SD method to simulate multi-species reacting flows is demonstrated. The present work is a starting point for developing the SD method on combustion applications. Future work will focus on the development of p-refinement techniques since local high values of p improve the results very well. Comparisons between p-refinement and the AMR technique on several reacting flows configuration is the next challenge.

Transformation Matrices from Conservative to Primitive Variables
Let's denote by U = , u, v, w, E, Y k T the vector of conservative variables and Q = , u, v, w, P, Y k T the vector of primitive variables with density as the first variable and pressure as the fifth variable. For a multi-species thermally perfect gas, E is the sum of sensible and kinetic energies per unit of volume (Poinsot and Veynante 2005): C pk T � dT � , C pk being the heat capacity at constant pressure of species k and T 0 a reference temperature and ||u|| 2 is the L 2 -norm of the velocity vector. Moreover, the multi-species gas is also assumed to behave as an ideal gas so that Eq. (7) holds true. The passage matrix from primitive to conservative variables U∕ Q for a multi-species thermally perfect gas is given by: Matrix Q∕ U is the inverse matrix of U∕ Q introduced in Eq. (39): Compare to a calorically perfect gas, the four first lines of ( U∕ Q) are exactly the same: the differences lie in the fifth line because the expression of E is different for a thermally perfect gas.

Transformation Matrix from Primitive to Characteristic Variables
From a wave analysis of Eq. (10) written for a non-viscous and non-reacting multi-species gas, a transformation matrix from primitive to characteristic variables can be found (Hirsch 1988): where c is the local speed of sound and n u = n u x , n u y , n u z T is the face unit normal vector (the normal at a given FP for instance) taken as ∇ ∕||∇ || for a boundary located at a -normal face. Its inverse can also be employed, for instance when computing the system of equations for a NSCBC inlet derived in Appendix B, so it is stated here for clarity: n u x 0 n u z ∕c − n u y ∕c − n u x ∕c 2 0 n u y − n u z ∕c 0 n u x ∕c − n u y ∕c 2 0 n u z n u y ∕c − n u

Presentation of the Case and Numerical Setup
To go a step further in the validation, a 2D laminar methane-air premixed flame is now computed. The geometry is a 2D burner shown in Fig. 12. Fresh gases enter the burner axially (x-axis) through a NSCBC subsonic inflow imposing a parabolic profile given by Eq. (50): where u 0 = 4 m s −1 and l 0 = 0.65 mm. The two-reaction mechanism is again the CH4/ Air-2S-BFER (Franzelli et al. 2010) with equivalence ratio and fresh gas temperature also set respectively to 0.8 and 300 K at the inlet so that flame characteristics are the ones in Table 2. At the outlet, a subsonic outflow imposing P = 101325 Pa is employed. Side walls are modeled as adiabatic no-slip walls. In SD, there are several methods to specify a wall boundary condition at an interface I due to the fact that only one state is available at this interface. This condition is enforced weakly by imposing values of convective and diffusive fluxes at interface FP on the wall. Assuming that the boundary left state Q L FP ( Q = T, u, P, Y k T here) is the known state from the interior domain, the ghost state Q R FP is set in this work as: where T w , P w and Y k,w are wall temperature, pressure and species mass fractions computed such that normal temperature, pressure and species gradients are zero at the interface FP. Once left and right states are known, the convective and diffusive fluxes at interface I are computed as follows: with (∇Q) R FP the gradient of each variable in Q FP computed using the values of Q R FP on the wall. For the JAGUAR simulation, the domain is discretized using 1216 uniform quadrilateral elements of characteristic size Δ e = 2.27 × 10 −4 m. The polynomial degree within each element is set to p = 4 ending up with 30,400 DOF inside the SD mesh and eight points inside the flame front. The AVBP solution domain is still discretized with almost the same number of DOF (31,574 nodes in that case) to have a fair comparison. Finally as for the 1D flame, the centered diffusion scheme (Sun et al. 2007) and the third-order in time SSP RK scheme of Gottlieb and Shu (1998) were used respectively for interface diffusion fluxes and for time integration. Figure 13 shows the 2D heat release rate field obtained with the JAGUAR simulation. The flame structure of a 2D burner is retrieved. The theoretical flame length can be calculated using Eq. (53) (Nambully et al. 2014b): In Eq. (53), U b = 2.667 m s −1 , = 0.106 (computed with S 0 L = 0.2815 m s −1 ) so that L th f = 6.11 mm which is very close to the value of 6 mm for the JAGUAR simulation which can be deduced from Fig. 13. Vertical profiles of both JAGUAR and AVBP solutions at x = 12 mm in Fig. 14 are very similar. Horizontal profiles along the centerline y = 0 are also represented in Fig. 15 for completeness showing also good agreement between the two solvers. Finally, Table 5 shows that both codes predict the same burning temperature T b and almost the same maximum of heat release rate in the whole domain. CANTERA results obtained for the equivalent 1D flame are also added as reference values. In terms of computational cost in that case, JAGUAR has an efficiency (cost per iteration and per DOF) around 7.5 μs∕ite∕DOF whereas AVBP using TTGC scheme has an efficiency around 6.0 μs∕ite∕DOF . This is quite encouraging since JAGUAR is a very recent code where no real optimization was done compared to AVBP, which has been under continuous development for more than 20 years.

Appendix 4: Comparison Between the Centred Diffusion Scheme and the SDLIFT Diffusion Scheme on a 1D Flame with a Coarse Discretization
The same 1D flame introduced in Sect. 5.1, whose characteristics are given in Table 2, is considered but this time with 185 DOF, instead of 400, set with N e = 37 elements at p = 4 . The number of points in 0 L using this discretization is around 4 which is close to the usual acceptable limit for discretizing a flame front (Poinsot and Veynante 2005). The simulation is run with either the centered diffusion scheme (Sun et al. 2007) or the SDLIFT formulation described in Sect. 3.4 using the IC with g DG correction function. It should be mentioned that the positivity-preserving limiter for mass fractions was deactivated to focus on the impact of the diffusion scheme on the results. Figure 16 shows the mass fraction   (b) Heat release rate.   Comparison between JAGUAR and AVBP of temperature, heat release rate (zoomed between x min = 1.54 cm and x min = 1.66 cm), axial and vertical velocity profiles at y = 0 mm along x-axis for the 2D burner case of CO obtained with both diffusion schemes showing that the centered diffusion scheme produces severe oscillations across the flame front and predict a much higher value of Y CO in the burnt gases. On the contrary, the SDLIFT formulation is very stable and predicts the correct value of Y CO in the burnt gases. This simple test case demonstrates that the centered diffusion scheme is not adapted when the flame front is not sufficiently discretized which can happen in the case of three-dimensional turbulent flames, even if a combustion model is used. The SDLIFT formulation seems more appropriate to gain in stability for such configurations with coarser meshes.
Author Contributions TM writing-original draft, methodology, software, investigation, validation. HD writing-original draft, methodology, software, investigation, validation. J-FB methodology, review, editing, supervision. BC methodology, review, editing, supervision. RM methodology, review, editing, supervision Funding This work was conducted under a Cifre contract funded by Safran and the Association nationale de la recherche et de la technologie (ANRT) with the research Project Number 2019/0794.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethical Approval This material is the authors' own original work, which has not been previously published elsewhere. It reflects the authors' own research and analysis in a truthful and complete manner.
Informed Consent All authors agreed with the content and that all gave explicit consent to submit and that they obtained consent from the responsible authorities before the work was submitted.  Fig. 16 Comparison of mass fractions of CO obtained with the centred diffusion scheme (Sun et al. 2007) and the SDLIFT formulation described in Sect. 3.4 using the IC with g DG correction function on a 1D flame with four points inside the flame front