Gauthier Wissocq

Said Taileb

Song Zhao

Pierre Boivin

A hybrid lattice Boltzmann method for gaseous detonations

This article is dedicated to the construction of a robust and accurate numerical scheme based on the lattice Boltzmann method (LBM) for simulations of gaseous detonations. This objective is achieved through careful construction of a fully conservative hybrid lattice Boltzmann scheme tailored for multi-species reactive flows. The core concept is to retain LBM low dissipation properties for acoustic and vortical modes by using the collide and stream algorithm for the particle distribution function, while transporting entropic and species modes via a specifically designed finite-volume scheme. The proposed method is first evaluated on common academic cases, demonstrating its ability to accurately simulate multi-species compressible and reactive flows with discontinuities: the convection of inert species, a Sod shock tube with two ideal gases and a steady one-dimensional inviscid detonation wave. Subsequently, the potential of this novel approach is demonstrated in one-and two-dimensional inviscid unsteady gaseous detonations, highlighting its ability to accurately recover detonation structures and associated instabilities for high activation energies. To the authors' knowledge, this study is the first successful simulation of detonation cellular structures capitalizing on the LBM collide and stream algorithm.

Introduction

The aeronautics and aerospace industries are currently experiencing significant transformations driven by the demand for groundbreaking technologies that offer improved efficiency and a reduced environmental impact. Developing new engines based on high-speed combustion, such as detonation, holds great promise due to the potential for increased pressure gain compared to traditional combustion methods. This concept is inspired by the recognition that the thermodynamic cycle known as the Fickett-Jacobs cycle surpasses the efficiency of the classical Brayton cycle [START_REF] Kailasanath | Recent developments in the research on pulse detonation engines[END_REF][START_REF] Schwer | Numerical investigation of rotating detonation engines[END_REF]. Various applications, including rotating detonation engines (RDE), have garnered significant attention in recent years. [START_REF] Anand | Rotating detonation combustors and their similarities to rocket instabilities[END_REF]. Attaining a comprehensive grasp of the physical phenomena taking place within such engines is imperative for optimizing the design process. However, gaseous detonation waves encompass reactive and compressible chaotic gas-dynamical events. Their observation reveals intricate patterns involving multidimensional and unsteady structures [START_REF] Soloukhin | Multiheaded structure of gaseous detonation[END_REF][START_REF] Soloukhin | Nonstationary phenomena in gaseous detonation[END_REF], thereby rendering precise experimental studies extremely challenging.

In this regard, computational fluid dynamics (CFD) presents an alternative approach with the aim of accurately capturing the flow physics. However, significant challenges still have to be addressed, further narrowing the selection of suitable numerical methods. Detonation waves encompass characteristic scales that span several orders of magnitude: chemical length scales are approximately one micrometer, the reactive fronts display inherent instability and form cellular structures spanning from millimeters to centimeters in length, while the applications involving detonation typically extend to the meter scale. Furthermore, the detonation front propagates at velocities on the order of several kilometers per second. Incorporating a detailed chemical reaction description compounds the complexity of the application, making it even more demanding. The intricate interplay of shocks, chemical reactions, and turbulence therefore requires the use of very robust and low dissipative numerical methods [START_REF] Lele | Shock-turbulence interaction: what we know and what we can learn from peta-scale simulations[END_REF]. To tackle these challenges, various strategies have been embraced in the literature, ranging from the use of high-order finite volume schemes based on HLLC (Harten-Lax-van Leer-Contact) solvers [START_REF] Sow | Mean structure of one-dimensional unstable detonations with friction[END_REF][START_REF] Hu | A modified fifth-order wenoz method for hyperbolic conservation laws[END_REF][START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF][START_REF] Sato | Detailed chemical kinetics based simulation of detonation-containing flows[END_REF] to Galerkin discontinuous methods [START_REF] Wang | Robust high order discontinuous galerkin schemes for two-dimensional gaseous detonations[END_REF][START_REF] Du | High-order bound-preserving discontinuous galerkin methods for stiff multispecies detonation[END_REF]. Despite significant enhancements in modeling strategies through adaptive mesh refinement [START_REF] Quirk | Amrita-a computational facility(for cfd modelling)[END_REF][START_REF] Houim | A low-dissipation and time-accurate method for compressible multi-component flow with variable specific heat ratios[END_REF][START_REF] Cai | Adaptive simulations of viscous detonations initiated by a hot jet using a high-order hybrid weno-cd scheme[END_REF], simulating an entire detonation engine remains computationally challenging with these high-order approaches.

From this perspective, the lattice Boltzmann method (LBM) represents a competitive alternative approach that has gained increasing attention from the CFD community over the last three decades. Its primary strengths stem from a simple numerical scheme, achieved by splitting the Boltzmann equation with discrete velocities into two parts: a local collision term mimicking inter particles interactions and a node-to-node streaming step consisting in a mere memory shift [START_REF] Krüger | The Lattice Boltzmann Method[END_REF]. Although second-order accurate in space and time, it has been shown that the standard LBM can be a low dissipative, efficient and accurate numerical method for solving fluid mechanics problems [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF]. These features have made it appealing for academic and industrial applications alike, encompassing fields such as turbulent flows [START_REF] Yu | Lattice boltzmann simulations of decaying homogeneous isotropic turbulence[END_REF][START_REF] Wilhelm | A new explicit algebraic wall model for les of turbulent flows under adverse pressure gradient[END_REF], combustion [START_REF] Lin | Double-distribution-function discrete boltzmann model for combustion[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF], multiphase flows [START_REF] Mazloomi | Entropic lattice boltzmann method for multiphase flows: Fluid-solid interfaces[END_REF][START_REF] Lafarge | Improved color-gradient method for lattice boltzmann modeling of two-phase flows[END_REF] and magneto-hydrodynamics [START_REF] Dellar | Moment equations for magnetohydrodynamics[END_REF]. Nevertheless, conventional lattice Boltzmann methods, rooted strictly in a kinetic model derived from the Boltzmann equation, still face difficulties when simulating the compressible and reactive flows inherent to detonation waves. One difficulty arises from the need to model the energy equation for compressible flows, which necessitates an increase in the number of microscopic velocities to be considered [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]. This, in turn, expands the spatial stencil of the method, leading to heightened complexity in the numerical scheme, a significant rise in memory requirements, and makes numerical stability more challenging [START_REF] Wissocq | An extended spectral analysis of the lattice Boltzmann method: modal interactions and stability issues[END_REF]. Admittedly, recent work have successfully made it possible to simulate detonation waves with kinetic approaches grounded in the Boltzmann equation. To the best of the authors' knowledge, Yan et al. were pioneers in proposing a lattice Boltzmann model for detonations, employing a finite difference lattice Boltzmann (FDLB) scheme [START_REF] Yan | Lattice Boltzmann model for combustion and detonation[END_REF]. More recently, Ji et al. introduced a discrete Boltzmann model (DBM) for compressible reactive flows, allowing them to simulate steady and unsteady detonations [START_REF] Ji | A three-dimensional discrete Boltzmann model for steady and unsteady detonation[END_REF]. Sawant et al. achieved successful simulations of detonation waves with the particles on demand (PonD) method [START_REF] Sawant | Detonation modeling with the particles on demand method[END_REF]. Despite their promise, these recent approaches do not fully leverage the advantages of the LBM, which lie in the simplicity and low dissipative feature of the collide and stream algorithm.

An alternative avenue for expanding the capabilities of LBM to encompass thermal, compressible, and reactive flows is through the implementation of a hybrid LBM. This concept stems from the observation that increasing the number of lattice velocities can lead to computational expenses, especially when incorporating a supplementary scalar equation like temperature evolution. The purpose of the hybrid LBM is, therefore, to explicitly solve a suitable form of the energy equation using a finite-difference (FD) or finite-volume (FV) scheme, sharing the same mesh and time step as the LBM. The appeal of this approach unfolds in three key aspects. First, it enables the consideration of compressible flows without increasing the spatial stencil and the memory cost due to the add of distributions functions in multi-speed [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF][START_REF] Philippi | From the continuous to the lattice Boltzmann equation: The discretization problem and thermal models[END_REF][START_REF] Frapolli | Multispeed entropic lattice Boltzmann model for thermal flows[END_REF][START_REF] Frapolli | Entropic lattice Boltzmann model for compressible flows[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF][START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF] or double distribution functions (DDF) [START_REF] He | A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit[END_REF][START_REF] Guo | Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model[END_REF][START_REF] Li | Coupled double-distribution-function lattice Boltzmann method for the compressible Navier-Stokes equations[END_REF][START_REF] Feng | A compressible lattice Boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF] approaches. Second, it allows for explicit specification of the gas heat capacity ratio and, in a broader sense, its equation of state. This is achievable as the internal energy is no longer restricted solely to translational energy. With fully kinetic approaches, the simple extension to polyatomic gases is acknowledged to be quite challenging [START_REF] Dellar | Two routes from the Boltzmann equation to compressible flow of polyatomic gases[END_REF]. Third, in the context of reactive flows, this strategy can be used to explicitly solve the mass fraction equation for each species, thereby sidestepping the additional memory cost of using multiple distribution functions, as done in the aforementioned methods. The hybrid LBM also opens the door to well-controlled chemical kinetics. Owing to these advantages, it has recently proven its worth for simulations of compressible [START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF][START_REF] Li | Application of a Higher Order Lattice Boltzmann/ Hybrid Method for Simulation of Compressible Viscous Flows with Curved Boundary[END_REF][START_REF] Ribeiro | Buffet Simulations with a Lattice-Boltzmann based Transonic Solver[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Renard | Improved compressible hybrid lattice Boltzmann method on standard lattice for subsonic and supersonic flows[END_REF][START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF][START_REF] Zhao | Toward fully conservative hybrid lattice Boltzmann methods for compressible flows[END_REF][START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF] and reactive [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF] flows.

An inquiry naturally emerges regarding the hybrid LBM: since part of the system (commonly energy and mass fraction equations) is numerically solved by a FD/FV scheme, is it possible to preserve the low dissipation property of the LBM even when it is not coupled with a high-order scheme? One can indeed expect the coupling to deteriorate the numerical properties of the LBM alone. Addressing this query involves the potential separation of the physical waves, commonly referred to as Kováznay modes [START_REF] Chu | Non-linear interactions in a viscous heat-conducting compressible gas[END_REF], which can be achieved when the hybrid LBM is based on the characteristics of the hyperbolic Euler system. Indeed, as discussed in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], when the LBM is supplemented by a characteristic wave equation solved through a FD/FV scheme, the two systems become linearly decoupled, so that the independent numerical properties of each scheme can be preserved. In practice, this leads to the choice of the entropy equation (which is precisely a characteristic one) as initially proposed by Nie et al. [START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF], to consider temperature fluctuations. Recent endeavors have successfully demonstrated that, when this strategy is adopted, the remarkable accuracy of the LBM can be preserved for any isentropic flow (e.g. acoustic propagation or vortical fields), while the numerical scheme adopted for the entropy equation only affects the behavior of entropy waves [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF].

The entropy-based hybrid LBM has therefore successfully extended the scope of the method to compressible and reactive flows, as discussed in [START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF]. Unfortunately, this model is not suitable for compressible reactive flows in presence of discontinuities such as shock waves, which is required for the simulation of detonations. The issue arises from the lack of conservativity of this approach, which fails to yield correct jump relations across a discontinuity [START_REF] Zhao | Toward fully conservative hybrid lattice Boltzmann methods for compressible flows[END_REF]. A conventional way to ensure the global conservativity of a numerical scheme is to deal with a FV discrete formulation of the total energy equation [START_REF] Hou | Why Nonconservative Schemes Converge to Wrong Solutions: Error Analysis[END_REF]. However, integrating this formalism with the hybrid LBM proves challenging, since it is prone to triggering severe numerical instabilities, as discussed in earlier references [START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF][START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF][START_REF] Guo | On the use of conservative formulation of energy equation in hybrid compressible lattice Boltzmann method[END_REF]. Additionally, the aforementioned benefits of dealing with a characteristic equation such as the entropy one may be lost when working with a total energy formulation. Recently, a solution to this predicament was proposed by Wissocq et al. [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF] for non-reactive flows. The idea is to construct FV total energy schemes maintaining all the linear properties of a nonconservative entropy-based hybrid system. Through this, new numerical schemes could be successfully built for the hybrid LBM, satisfying the following properties: 1) conservativity is ensured by construction, 2) linear stability can be controlled by the choice of FV scheme, 3) the low dissipation of the LBM for isentropic phenomena can be preserved. The aim of the present work is to extend this methodology to simulations of multi-species reactive flows. Our objective is to build conservative, robust and low-dissipative numerical schemes based on the hybrid LBM for compressible reactive flows in order to accurately simulate detonation waves.

The article is organized as follows. Section 1 details the theoretical derivation of the conservative scheme for multi-species reactive flows in a general way. Mathematical details for this work are provided in Appendix A, Appendix B and Appendix C. In Section 2, the models retained in the present work for the simulation of detonation waves are introduced. They are based on the so-called unified lattice Boltzmann collision model, initially proposed by Farag et al. [START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF], and further detailed in Appendix D. Three conservative discretizations of the energy equation are proposed. They are further compared and validated on standard academic cases of growing complexity in Section 3: the convection of inert species, a shock tube with two gases and a stable one-dimensional detonation wave. In Section 4, the potential of this new approach is demonstrated on strongly unsteady inviscid detonations in one and two dimensions.

Construction of conservative hybrid LB schemes for multi-species reactive flows

The aim of this section is to detail the construction of stable conservative schemes for the hybrid LBM in presence of multi-species reactive flows. It is an extension of a previous work dedicated to the non-reactive case [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF].

In the present work, the inviscid case, modelled by the Euler equations, is considered. Unless otherwise stated, a dummy field Φ evaluated at (x,t), where x is the position and t is the time, is simply denoted as Φ. In order to make the demonstration as general as possible, no assumption is made on the spatial dimension which is represented by Greek indices obeying Einstein summation convention. For instance, the product u α δ α F α reads, in two dimensions: u x δ x F x + u y δ y F y . Also note that all vector quantities are written in bold.

Non-conservative model and definitions

We first consider a hybrid model where mass and momentum equations are numerically solved in their conservative form by a LB scheme, and coupled with non-conservative formulations of the entropy and mass fraction equations. This description is particularly suitable for previous hybrid LB models of the literature [START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Renard | Improved compressible hybrid lattice Boltzmann method on standard lattice for subsonic and supersonic flows[END_REF][START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Guo | Improved standard thermal lattice Boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows[END_REF]]. In the one-dimensional case, the system reads:

δ t ρ + F ρ,LB +∆x/2 -F ρ,LB -∆x/2 /∆x = 0, δ t (ρu) + F ρu,LB +∆x/2 -F ρu,LB -∆x/2 /∆x = 0, δ t s + u F s,FD +∆x/2 -F s,FD -∆x/2 /∆x = ωs , δ t Y k + u F Y k ,FD +∆x/2 -F Y k ,FD -∆x/2 /∆x = ωk , ∀k ∈ 2, N s , (1)
where (ρ, u, s,Y k ) are respectively the mass density, the fluid velocity, the specific entropy and the volume fraction of the species k, N s is the number of species, ωs and ωk respectively stand for the net productions of entropy and mass fraction of the species k due to chemical reactions. In practice, the mass fraction of the most inert species (here the first one) will not be transported and may be evaluated from

Y 1 = 1 -∑ N s 2 Y k if necessary.
The operator δ t is a numerical discretization of time derivatives, defined for a dummy field Φ as

δ t Φ = Φ(x,t + ∆t) -Φ(x,t) ∆t , ( 2
)
where ∆t is the time step and ∆x is the mesh size. Finally, F Φ are numerical fluxes of a quantity Φ. They can be either computed by a LB scheme ("LB" superscript) or by a finite-difference scheme ("FD" superscript). The index +∆x/2 indicates a flux estimated at the right boundary of a given cell, while the index -∆x/2 indicates the left boundary. A conservative computation of the fluxes yields

F Φ -∆x/2 (x,t) = F Φ +∆x/2 (x -∆x,t), (3)
so that the conservative scheme can be equivalently written involving the fluxes computed at the right boundary (F +∆x/2 ) only. The complete system can then be generalized to any space dimension under the following matrix form:

δ t V + A α δ α F V,d α = S V , (4)
where an implicit summation is assumed on the index α representing the spatial dimensions,

V = [ρ, ρu, s,Y k ] T is a vector of non-conserved quantities, F V,d α = [F ρ,LB +∆α/2 , F ρu,LB +∆α/2 , F s,FD +∆α/2 , F Y k ,FD +∆α/2 ]
T is the vector of discrete fluxes in the direction α, S V = [0, 0, ωs , ωk ] T is the vector of source terms, δ α is a discrete spatial gradient operator defined as

δ α Φ = Φ(x,t) -Φ(x -e α ∆x,t) ∆x , (5)
e α is the unity vector in the direction α and A α is a diagonal matrix whose coefficients are 1 for the mass and momentum equations, and u α for the non-conservative entropy and volume fraction equations. For instance, in two dimensions with N s = 2, this reads

A α =       1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 u α 0 0 0 0 0 u α       . (6)
for α ∈ {x, y}. The presence of this non-identity matrix A α is a clear indication of the non-conservativity of the system.

With the FV formulation of Eq. ( 4), the choice of an explicit numerical scheme for each equation boils down to the choice of a discrete flux.It is worth noting that in the context of the entropy and mass fraction equations, the term "flux" is not appropriately used as they are not expressed in a conservative form. Nevertheless, these fluxes are expected to be consistent with their continuous counterpart so that the consistency of the whole scheme can be ensured. This reads

F V,d α = F V,c α + O(∆x, ∆t), (7)
where F V,c α are the continuous Euler fluxes defined as

F V,c α =     ρu α ρu α u β + pδ αβ s Y k     , ( 8
)
where p is the thermodynamic pressure and δ αβ is the Kronecker symbol. In order to close the system, an equation of state has to be prescribed between p and other known quantities, i.e. p = p(ρ, s,Y k ). For instance, for an ideal single species it reads

s = c v ln p ρ γ , (9)
where c v and γ are respectively the heat capacity at constant volume and the heat capacity ratio of the fluid.

In practice, dealing with a system like Eq. ( 4) is not common for reactive flows owing to two difficulties stemming from the chosen entropy equation: calculating entropy might not be straightforward for real gas mixtures and the presence of the source term ωs increases the complexity of the ensuing numerical scheme. Consequently, reactive flows typically favor using enthalpy or total energy equations [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. However, the present work aims to initiate from a generalized form of the hybrid system based on the entropy equation, as it offers a distinct advantage over other equations: the evolution equations for entropy and mass fraction assume a characteristic form, manifesting as straightforward advection equations for quantities s and Y k at velocity u. As stated in previous work [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF][START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF][START_REF] Guo | On the use of conservative formulation of energy equation in hybrid compressible lattice Boltzmann method[END_REF], this characteristic nature yields two key benefits: 1) the characteristic equations remain linearly decoupled with each other, affording enhanced stability and accuracy control for the entire system, and 2) a substantial literature exists on numerical discretizations of advection equations [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. This is why coupling the LBM with a non-conservative entropy equation has been the favored approach for compressible non-reactive flows [START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Renard | Improved compressible hybrid lattice Boltzmann method on standard lattice for subsonic and supersonic flows[END_REF][START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Guo | Improved standard thermal lattice Boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows[END_REF]. However, this non-conservative formulation is plagued by a significant issue: it leads to incorrect solutions close to flow discontinuities [START_REF] Zhao | Toward fully conservative hybrid lattice Boltzmann methods for compressible flows[END_REF]. This flaw stems from the fact that solely fluxes of conserved quantities -namely mass, momentum, total energy and mass of each speciesremain continuous across a discontinuity. To address this issue, it is imperative to deal with conserved formulations, which are unfortunately likely to lead to coupling instabilities as shown in a previous work [START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF]. Two solutions have recently emerged to ensure the stability of conservative hybrid LB schemes: Zhao et al. developed a new construction of the fluxes of conserved quantities drawing inspiration from the way mass and momentum fluxes are computed by the LBM [START_REF] Zhao | Toward fully conservative hybrid lattice Boltzmann methods for compressible flows[END_REF], and Wissocq et al. derived a conserved total energy flux preserving the linear properties of the nonconservative system [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF]. In the following, an extension of the latter approach is proposed for multi-species reactive flows starting from Eq. ( 4).

Mass and momentum equations being solved, in the present work, by a LB scheme, the expression of their fluxes F ρ,LB +∆α/2 and F ρu α ,LB +∆α/2 are provided in Sec. 2.2 and Appendix C depending on the collision model adopted for the LBM. Regarding the entropy and mass fraction fluxes, we will assume they can be written under the following form:

F s,FD +∆α/2 = F α (s), F Y k ,FD +∆α/2 = G α (Y k ), (10)
where F α and G α are functions of the vector fields (in space) s and Y k . The most straightforward example of such an operator in one dimension is to consider F x (s) = s(x,t), which would lead to

δ x F s,FD +∆x/2 = s(x,t) -s(x -∆x,t) ∆x , (11)
so that the numerical scheme for the entropy is a first-order left-sided Euler scheme 1 . More advanced (stable and second-order) choices for these functions are provided in Sec. 2.3.

Construction of a conservative system

Following [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], the objective of the present section is to find a linearly equivalent system of Eq. ( 4) based on the conserved variables U = [ρ, ρu, ρE, ρY k ] T , where E is the total energy by unit of mass defined as

E ≡ H - p ρ , with H ≡ h + κ, h ≡ N s ∑ k=1 Y k h k , h k ≡ T T 0 c p,k (θ )dθ + ∆h 0 k , ( 12
)
where H is the total enthalpy, h is the specific enthalpy of the mixture, κ = ||u|| 2 /2 is the kinetic energy by unit of volume, h k is the enthalpy of each species, c p,k is the heat capacity at constant pressure of each species, T is the gas temperature and ∆h 0 k is the enthalpy of formation of the species k at temperature T 0 . Note that with this definition, total energy naturally includes the chemical potential of the gas mixture, so that it is a conserved quantity for multi-species reactive flows.

The knowledge of the Jacobian matrix M = ∂ U/∂ V is a key point in the passage from non-conserved quantities V to conserved ones U. The expression of this matrix depends on the equation of state under consideration. Notably, a particular expression of this matrix can be obtain from the first law of thermodynamics as shown in Appendix A. Multiplying the discretized form of the non-conservative system (4) by M yields, formally,

δ t U + MA α δ α F V,d α = MS V ≡ S U . (13)
A key point of the proposed methodology is that Eq. ( 13) is linearly equivalent to Eq. ( 4). Even though this system describes the discrete evolution of conserved variables U, it is not a conservative formulation because of the presence of the non-identity matrix MA α in front of the discrete flux. According to Appendix B, it is possible to transform this system without affecting its linear properties so as to obtain

δ t U + δ α F U,d α = S U , (14)
where

F U,d α = F ρ,LB +∆α/2 , F ρu,LB +∆α/2 , F ρE,FV +∆α/2 , F ρY k ,FV +∆α/2
T is a vector of conserved fluxes.

F ρE,FV +∆α/2 , F ρY k ,FV
+∆α/2 are the conserved numerical fluxes of total energy and mass fraction of each species. Their general expression is provided in Appendix B involving two arbitrary operators K 1,α and K 2,α . In what follows, for the sake of simplicity, one sets K 1,α = K 2,α = G α = F α , meaning that same spatial discretizations for the total energy and the mass fraction are assumed. Together with the expression of M provided in Appendix A, this leads to:

F ρE,FV +∆α/2 = F α (ρHu α ) (i) + (h -κ) F ρ,LB +∆α/2 -F α (ρu α ) (ii) + u β F ρu β ,LB +∆α/2 -F α (ρu α u β + pδ αβ ) (iii) , (15)
F ρY k ,FV +∆α/2 = F α (ρY k u α ) (iv) +Y k F ρ,LB +∆α/2 -F α (ρu α ) (v) , (16)
where an implicit summation is done over β . Total energy and mass fraction evolutions are now written under a conservative form, so that conservativity of the whole system can be ensured by construction. We now discuss on the role played by each term in the above equations.

(i) This term is nothing but the expected total energy flux ρHu α discretized by the numerical scheme denoted as F α , initially designed to compute the entropy gradients. As shown in a previous study [START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF], it is very likely that this term alone would lead to linearly unstable numerical schemes. The other terms in the total energy flux ((ii) and (iii)) can therefore be considered as numerical corrections to ensure linear stability.

(ii) Thanks to Eq. ( 7), one has F ρ,LB +∆α/2 -F α (ρu α ) = O(∆x, ∆t). As a consequence, this term does not affect the consistency of the total energy flux. Its role is to correct numerical errors induced by the use of different schemes for the mass and the total energy equation.

(iii) In the inviscid case, Eq. ( 7) leads to F ρu β ,LB +∆α/2 -F α (ρu α u β + pδ αβ ) = O(∆x, ∆t). The role of this term is then to correct numerical errors induced by the use of different schemes for the momentum and the energy equation. In the viscid case, this term is consistent with u β σ αβ , where σ αβ is the shear stress tensor. As shown in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], this leads to an implicit computation of the viscous heating term of the energy equation.

(iv) This term is the expected mass fraction flux ρY k u α discretized by the numerical scheme denoted with F α . As for (i), it is likely that this term alone would lead to a linearly unstable numerical scheme.

(v) Eq. ( 7) leads to

F ρ,LB +∆α/2 -F α (ρu α ) = O(∆x, ∆t).
This term is then a numerical correction due to the use of different schemes for the mass fraction and the total mass equations.

Finally, the location where pre-factors (hκ), u β and Y k are evaluated can be open to a discussion. Noticing that they are all multiplied by a O(∆x, ∆t) term, this choice does not affect the consistency of the scheme, nor its linear properties. To systematically ensure the symmetry of the scheme, these quantites will be extrapolated at the cell interface by averaging their value between x and (x + e α ∆x).

The crucial feature of the proposed construction is that all the steps leading to Eq. ( 14) are designed so as to preserve the linear properties of the initial scheme (4), which allows decoupling the characteristics of the system. As shown in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], this leads to two important properties shared by the conservative schemes:

• the accuracy in the transport of entropy and mass fraction can be explicitly and independently controlled by the choice of operator F α ,

• the linear properties of the LBM are preserved for acoustic and vortical phenomena.

Especially thanks to the second point, the low dissipation of the LBM for pressure waves can be preserved [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF][START_REF] Wissocq | Linear stability and isotropy properties of athermal regularized lattice Boltzmann methods[END_REF].

Adding chemical source terms

We now focus on the source term S U appearing in the conservative formulation [START_REF] Houim | A low-dissipation and time-accurate method for compressible multi-component flow with variable specific heat ratios[END_REF]. It is noteworthy that, given the definition provided by Eq. ( 12), sensible enthalpy variations induced by the chemical reactions are naturally absorbed in the total energy. In other words, this total energy is a conserved quantity and no source term related to chemical reactions are expected to appear in its evolution equation. As a consequence, the source term reads

S U = [0, 0, 0, ρ ωk ] T , (17)
and only the mass fraction equations are affected by a a net production due to chemical reactions. Without loss of generality, it can be written as a function of the conserved variables as ωk (U). In the present work, this source term is numerically implemented thanks to the following splitting approach [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]:

(ρY k ) * (x,t + ∆t) = (ρY k )(x,t) -∆tδ α F ρY k ,FV +∆α/2 , (ρY k )(x,t + ∆t) = (ρY k ) * (x,t) + ∆tρ(x,t + ∆t) ωk (U * (x,t + ∆t)) , (18)
where

U * = [ρ, ρu, ρE, (ρY k ) * ] T .
Even though the computation of the source term involves quantities evaluated at time t + ∆t, the overall scheme can be made fully explicit by performing the transport of mass density ρ and total energy ρE before the update of species. A detailed example of the inclusion of source term ωk is illustrated in Sec. 2 for a single-step chemistry involving two species, which is the topic of the rest of the article.

A model for inviscid gaseous detonations

This section aims at detailing the particular choices adopted in the rest of the article for: 1) the composition of the gas mixture, 2) the lattice Boltzmann scheme used for mass and momentum equations, 3) the choice of discretization operator F α appearing in Eqs. ( 15)-( 16).

Assumptions: ideal gas mixture of two species

In the following, we consider a mixture of two calorifically perfect gases denoted with indices 1 and 2, subject to a chemical reaction 1 → 2. Index 1 is then referred to the fresh gas, and index 2 to the burnt gas. Unless otherwise stated, both species have equal and constant heat capacities, but different formation enthalpies, such as

p = ρRT h = Y 1 h 1 +Y 2 h 2 = c p T -Y 2 Q, ( 19
)
where Q is the heat of combustion corresponding to the reaction. The chemical source term corresponding to the net production of the burnt gas is modelled with the Arrhenius law [START_REF] Arrhenius | Über die Dissociationswärme und den Einfluss der Temperatur auf den Dissociationsgrad der Elektrolyte[END_REF],

ω2 = kY 1 exp - E a RT , ( 20
)
where k is the pre-exponential factor and E a the activation energy. In the following, the mass fraction and production rate of the burnt gases will be simply be denoted Y and ω.

Mass and momentum: lattice Boltzmann scheme

Mass and momentum equations are resolved using a lattice Boltzmann scheme, which can be written under the following general (collide and stream) formulation:

f i (x,t + ∆t) = f coll i (x -c i ∆t,t), ∀i ∈ 0,V -1 , (21)
where

( f i ) is a set of V distribution functions, ( f coll i
) are so-called post-collision distribution functions and c i are the discrete velocities of a lattice. In the following, the focus is put on two-dimensional simulations using the D2Q9 lattice [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] for which V = 9 and

c i,x = (0, 1, 1, 0, -1, -1, -1, 0, 1)∆x/∆t, c i,y = (0, 0, 1, 1, 1, 0, -1, -1, -1)∆x/∆t. ( 22
)
Mass density ρ and momentum ρu α are then defined as moments of the distribution functions,

ρ(x,t) = ∑ i f i (x,t), ρu α (x,t) = ∑ i c i,α f i (x,t). (23)
As previously shown by Wissocq et al. [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], the discrete evolutions of mass and momentum can be written under the conservative form of Eq. ( 4), where the expressions of mass fluxes (F

ρ +∆x/2 , F ρ +∆y/2
) and momentum fluxes (F

ρu x +∆x/2 , F ρu x +∆y/2 , F ρu y +∆x/2 , F ρu y +∆y/2 ) are recalled in Appendix C. Regarding the post-collision distribution functions ( f coll i
), their computation depends on the collision model under consideration. In the present work and following [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], the unified density-based model of Farag et al. [START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] is adopted with a recursive regularized collision operator [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF]. It is detailed in Appendix D in the case of the D2Q9 lattice. It is noteworthy that performing a Chapman-Enskog expansion [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF] or a Taylor expansion [START_REF] Farag | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF][START_REF] Wissocq | Hydrodynamic limits and numerical errors of isothermal lattice Boltzmann schemes[END_REF], the relaxation time of the collision model τ can be related to the physical dynamic viscosity µ and an eventual artificial kinematic viscosity ν sc as

µ + ρν sc = τρc 2 s , (24)
where c s = 1 √ 3 ∆x ∆t is the so-called lattice sound speed and ∆t is the time step. In the present study, the shock sensor proposed in [START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF][START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF] is slightly modified as

ν sc = s c θ p(x -∆x,t) -2p(x,t) + p(x + ∆x,t) p(x -∆x,t) + 2p(x,t) + p(x + ∆x,t) , (25)
in the one-dimensional case, where θ = RT /c 2 s and s c is a free parameter whose value is specified in each case of Secs. 3-4. In multi-dimensions, the strategy proposed in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF] is adopted: ν sc is first computed as the average value of the above expression in each direction, then its local value is retained as the maximal value of all its direct neighbors. Additionnaly in the present work, a clipping is applied on the relaxation time as

τ = min µ + ρν sc ρc 2 s , ∆t 2 . ( 26
)
The idea behind this clipping is is rooted in the insight that, for τ = ∆t/2, most of the existing collision models (spanning from the Bhatnagar-Gross-Krook collision [START_REF] He | A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit[END_REF][START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] to the recursive regularization detailed in Appendix D) simplify to f coll i = f eq i +F E i ∆t/2, F E i being a general body-force term. As demonstrated in previous work [START_REF] Wissocq | Linear stability and isotropy properties of athermal regularized lattice Boltzmann methods[END_REF][START_REF] Wissocq | Hydrodynamic limits and numerical errors of isothermal lattice Boltzmann schemes[END_REF], this specific value results in a first-order accurate scheme which offers the advantage of filtering out any non-hydrodynamic mode from the flow physics, thereby enhancing stability properties. Therefore, Eq. ( 26) can be interpreted as a switch between a second-order scheme (when ν sc = 0) to a first-order one (when τ = ∆t/2), that can be used to better address pressure discontinuities. The free parameter s c then acts as a stiffness parameter for this switch. This particular feature will be assessed in the academic cases of Sec. 3.

Finite volume flux operators

This section details the choice of discretization operators F α involved in Eqs. ( 15)-( 16). Three schemes, referred here to as MHM 1D, MHM 2D and TVD Heun, will be assessed on academic cases. They are detailed below and their spatial stencil in two dimensions is displayed in Fig. 1.

MHM 1D

The first scheme, referred to as MHM (MUSCL-Hancock-Method) 1D in the present work, follows the work of van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow[END_REF] and Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. It was initially adopted for the compressible hybrid LBM by Farag et al. [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]. With the notations of the present article and following [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], it reads for a dummy field Φ,

F α (Φ) = Φ +∆α/2 (x,t) if u α (x,t) + u α (x + e α ∆x,t) ≥ 0, Φ -∆α/2 (x + e α ∆x,t) else, (27)
where

Φ +∆α/2 = Φ + ∆ α 2 -u α ∆t ∆x ∆ α 2 , Φ -∆α/2 = Φ - ∆ α 2 -u α ∆t ∆x ∆ α 2 , (28)
and where ∆ α is an approximation of the slope of Φ in the direction α, which can be given by

∆ α = 1 2 {(1 + η α ) [Φ(x,t) -Φ(x -e α ∆x,t)] + (1 -η α ) [Φ(x + e α ∆x,t) -Φ(x,t)]} . (29)
The parameter η α is set to η α = 1 3 2u α ∆t ∆xsign(u α ) as suggested in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] to obtain a third-order accurate convection scheme in space and time. Yet, as noted in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], the present algorithm is not effectively third-order accurate because of two reasons: (1) the time evolution of u should be considered to build a high-order scheme [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF], (2) in two dimensions, fluxes in the xand ydirections are computed separately, leading to a cross-shaped stencil (cf. Fig. 1) that can reduce the order of accuracy to the first-order for multi-dimensional flows. The second point is the reason why this scheme is referred to as MHM 1D. An improvement for multi-dimensions is proposed below.

MHM 2D

Noticing that the aforementionned MHM scheme can be unadapted to non-aligned flows, Yoo et al. [START_REF] Yoo | Compressible lattice boltzmann method with rotating overset grids[END_REF] recently proposed a slight improvement inspired by the work of Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. Compared to the model of Sec. 2.3.1, the only modification reads

Φ +∆α/2 = Φ + ∆ α 2 -u β ∆t ∆x ∆ β 2 , Φ -∆α/2 = Φ - ∆ α 2 -u β ∆t ∆x ∆ β 2 , (30)
where an implicit summation is done over the index β . For instance, in the two-dimensional case, this reads

Φ +∆x/2 = Φ + ∆ x 2 - ∆t 2∆x (u x ∆ x + u y ∆ y ) , Φ -∆x/2 = Φ - ∆ x 2 - ∆t 2∆x (u x ∆ x + u y ∆ y ) , (31)
Φ +∆y/2 = Φ + ∆ y 2 - ∆t 2∆x (u x ∆ x + u y ∆ y ) , Φ -∆y/2 = Φ - ∆ y 2 - ∆t 2∆x (u x ∆ x + u y ∆ y ) , (32)
where ∆ x and ∆ y are given by Eq. ( 29). Thanks to this modification, any flux computed in the xor y-direction involves both the slopes ∆ x and ∆ y , so that diagonal points are considered in this scheme as shown in Fig. 1. For this reason, it will be referred to as MHM 2D in the following.

TVD Heun

The two aforementioned schemes offer the advantage of a low numerical dissipation for the considered spatial stencil (maximum five points in each direction). Especially, they are linearly third-order accurate for flows aligned with the Cartesian mesh. Unfortunately, these schemes are not total variation diminishing (TVD), meaning that they are likely to generate spurious oscillations close to flow discontinuities, which may alter the detonation physics driven by shock waves. For this reason, it is also proposed to use a TVD scheme following the work of Dubois [START_REF] Dubois | Nonlinear Interpolation and Total Variation Diminishing Schemes (jun 2010)[END_REF], who adapted the MUSCL scheme of Heun with slope limiters, so that both TVD properties and second-order accuracy can be ensured. To the authors' knowledge, the Heun scheme is the only five-point-stencil scheme that enables the recovery of such properties. It is similar to the MHM 2D scheme of Sec. 2.3.2 with the following modification of the slope ∆ α :

∆ α = φα (r α ) (Φ(x + e α ∆x,t) -Φ(x,t)) , (33)
where φα , the slope limiter function, has to satisfy some conditions to ensure the TVD property and preserve the second-order accuracy [START_REF] Dubois | Nonlinear Interpolation and Total Variation Diminishing Schemes (jun 2010)[END_REF]. In the present work, it is proposed to use the absolute minmod limiter [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF], defined as

φα (r α ) = φ (r α ) if u α (x,t) + u α (x + e α ,t) ≥ 0, φ 1 r α else, (34)
φ (r) = |r| if |r| ≤ 1, 1 if |r| > 1, r α = Φ(x,t) -Φ(x -e α ∆x,t) Φ(x + e α ∆x,t) -Φ(x,t) . ( 35
)
Although not presented in the present article, some results obtained with other limiter functions tested in this work (in particular the Lagrange limiter of [START_REF] Dubois | Nonlinear Interpolation and Total Variation Diminishing Schemes (jun 2010)[END_REF]) will also be mentioned in the numerical validations of Sec. 3 below.

Summary of the hybrid scheme

The proposed fully conservative hybrid scheme for compressible reactive flows can be summarized as follows. Starting from the knowledge of the conserved fields ( f i ) i∈ 0,V -1 , ρ, ρu, ρE, ρY at time t, as well as (ρ, T ) a time (t -∆t):

1. Compute T from the knowledge of ρ, ρu, ρE and ρY at time t by inverting Eq. ( 12). 2. Compute ( f coll i ) at time t thanks to the recursive regularized collision model detailed in Appendix D. The knowledge of (ρ, T ) at time (t -∆t) is required in the computation of the correction term F E i . 3. Apply the streaming step to obtain ( f i ) at time (t + ∆t) (Eq. ( 21)). 4. Compute the LB mass and momentum fluxes from the knowledge of ( f coll i ) using Appendix C. 5. Compute ρ and ρu at time (t + ∆t) either as moments of f i (Eq. ( 23)) or using the LB fluxes (both methods are equivalent). 6. Total energy and mass fraction fluxes can then be computed using Eqs. ( 15)- [START_REF] Krüger | The Lattice Boltzmann Method[END_REF], where F α is a discretization operator proposed in Sec. 2.3. Remember that the pre-factors are evaluated at the mid-cell interface to ensure the symmetry of the scheme.

7. Total energy and mass fraction are transported as

(ρE)(x,t + ∆t) = (ρE)(x,t) - ∆t ∆x F ρE,FV +∆α/2 (x,t) -F ρE,FV +∆α/2 (x -e α ∆x,t) , (36)
(ρY ) * (x,t + ∆t) = (ρY )(x,t) - ∆t ∆x F ρY,FV +∆α/2 (x,t) -F ρY,FV +∆α/2 (x -e α ∆x,t) . ( 37
)
8. The source term of Eq. ( 20) is then introduced following the splitting approach of Sec. 1.3:

(ρY )(x,t + ∆t) = (ρY ) * (x,t + ∆t) + ∆tkρ(x,t + ∆t)(1 -Y * (x,t + ∆t)) exp - E a RT * (x,t + ∆t) , (38)
where Y * (x,t +∆t) = (ρY ) * (x,t +∆t)/ρ(x,t +∆t) and T * (x,t) is obtained in a similar way by inverting Eq. ( 12) with the knowledge of U * (x,t + ∆t).

All fields ( f i ) i∈ 0,V -1 , ρ, ρu α , ρE, ρY being then known at time (t + ∆t) as well as (ρ, T ) a time t, the procedure can be repeated until the last iteration. This scheme has been implemented in a massively parallel two-dimensional Fortran code and evaluated in the test cases of Sections 3 and 4 below.

Validation on academic cases

In this section, the proposed conservative scheme for compressible reactive flows is assessed on several one-and two-dimensional academic cases: the convection of inert species, a shock tube of two ideal gases and a stable onedimensional detonation. These test cases aim at numerically validating the stability and the conservativity of this approach, as well as guiding the choice between the schemes introduced in Sec. 2.3 for the simulation of gaseous detonations. In this regard, some models will be progressively discarded from the study.

All the test cases are simulated on a uniform mesh with an acoustic scaling [START_REF] Krüger | The Lattice Boltzmann Method[END_REF], meaning that the ratio ∆x/∆t is kept constant in the mesh convergence studies. For a given mesh size ∆x, the only numerical parameters are the shock sensor coefficient s c and the constant ratio ∆x/∆t, whose values will be specified in each case. The Courant-Friedrichs-Lewy (CFL) number [START_REF] Courant | On the Partial Difference Equations of Mathematical Physics[END_REF], defined as

CFL ≡ ||u|| + γRT ∆t ∆x , (39)
is not a constant prescribed parameter but an observed variable whose maximal value will be provided for each test case.

Convection of inert species

The first test case investigated in this work is the convection of inert (non reactive) species at a constant velocity. Its main purpose is to compare the numerical behavior of the schemes proposed in Sec. 2.3 for the transport of the mass fraction ρY , which behaves as a passive scalar. For this case, the gas mixture is characterized by γ = 1.4 and R = 1 and no chemical source term is considered ( ω = 0). A (L × L) fully periodic domain with L = 1 is discretized by a (200 × 200) Cartesian mesh. The simulation is initialized by density, pressure, velocity components and mass fraction as

ρ = 1, p = 1, u = (Ma x , Ma y ) √ γ, Y = exp - (x -x c ) 2 + (y -y c ) 2 R 2 c , (40)
where x c = y c = 0.5 are the coordinates of the center of the domain, R c = 0.1 is the radius of the Gaussian and Ma x and Ma y are the Mach numbers in the horizontal and vertical directions. Two cases with an absolute Mach number equal to unity are investigated:

1. horizontal convection: Ma x = 1, Ma y = 0, 2. diagonal convection: Ma x = 1/ √ 2, Ma y = 1/ √ 2.
Characteristic times corresponding to the flow-through-time period of the species can be respectively defined as t c = 1/ √ γ and t c = 2/γ. In both cases, we set ∆x/∆t = 7.9 so that max(CFL) ≈ 0.3. The shock sensor is not activated for this smooth case (s c = 0) and the fluid is considered inviscid (µ = 0). The exact solution is then supposed to be a convection of the initial Gaussian in mass fraction Y at the Mach number Ma. Fig. 2 displays isocontours of the mass fraction Y obtained with the three models introduced in Sec. 2.3 for the horizontal and diagonal convections at t = 20t c . For the horizontal convection, MHM1D and MHM2D lead to accurate results with a remarkable ability to transport species over long distances without transforming the Gaussian shape. In contrast, the Heun (TVD) scheme exhibits a dispersion error responsible for a distortion of the mass fraction, as well as more numerical dissipation. The observation that both MHM1D and MHM2D schemes lead to similar results is in agreement with the good properties of the MHM1D along the Cartesian directions. The diagonal points introduced by the MHM2D scheme have no contribution in this case. On the contrary, when the convection is performed along the diagonal, MHM1D rapidly results in numerical instabilities, which is in agreement with previous studies of an entropy spot convection [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF]. For this case, the add of diagonal points in the MHM2D scheme restores the stability and even more: the shape and amplitude of the Gaussian are especially well preserved. Regarding the Heun (TVD) scheme, a much larger dissipation is observed, even if the shape of the Gaussian is well conserved. Quantitative results regarding the amplitude loss for each case are provided in Table 1. Note that horizontal and diagonal simulations cannot be directly compared with each other because the characteristic time of the latter is larger in this case.

MHM1D MHM2D Heun (TVD)

Horizontal convection 0.980 0.980 0.688 Diagonal convection ✗ 0.957 0.465

Table 1: Amplitude of the mass fraction Y after 20t c along the horizontal and diagonal directions. To better explain these results, a convergence study is performed on similar cases for different meshes ranging from (100 × 100) to (3200 × 3200) nodes. Note that this test case is a pure advection at constant velocity without source term, for which the linear assumption is satisfied. The mass fraction is the only variable that is supposed to evolve with time. We therefore expect to recover the orders of convergence of the MHM 1D, MHM 2D and Heun (TVD) schemes as if they were not coupled with the LBM. At t = t c , the L 2 Y error on the mass fraction is computed as

L 2 Y = ∑ x,y [Y (x, y,t = t c ) -Y (x, y,t = 0)] 2 ∑ x,y Y (x, y,t = 0) 2 . ( 41
)
These errors are displayed in Fig. 3 with a log-log scale. Looking at the slope of each case leads to the following conclusions regarding the orders of accuracy:

• For the horizontal direction: MHM1D and MHM2D are both third-order accurate, Heun (TVD) is second-order accurate.

• For the diagonal direction: MHM1D is first-order accurate (and blows up with the finest mesh), MHM2D and Heun (TVD) are both second-order accurate.

In particular, the first-order accuracy of the MHM1D scheme for non-aligned flows leads to a viscosity-related numerical error which is responsible for the instability observed. For this reason, this model will be discarded from the rest of the study. Regarding the Heun (TVD) scheme, further studies indicate that the predominant numerical dissipation is due to the choice of the limiter function, namely the absolute minmod model. For the sake of completeness, other slope limiters have been tested during this work. In particular, the Lagrange limiter of Dubois [START_REF] Dubois | Nonlinear Interpolation and Total Variation Diminishing Schemes (jun 2010)[END_REF] leads to a much lower numerical dissipation, providing properties similar than the Heun scheme without limiter for this case. Unfortunately, another problem then occurs: the Heun scheme with Lagrange limiter proved to be too dispersive to be used in the following detonation simulations. In fact, the numerical dissipation introduced by the absolute minmod limiter has the advantage of reducing the dispersive behavior of the Heun scheme. This is why this model is retained for the following, together with the MHM2D scheme.

Sod shock tube with two ideal gases

The next test case aims to exhibit the behavior of the considered numerical schemes across flow discontinuities, as well as to validate the conservativity of the methods. To this end, the one-dimensional shock tube of Sod [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] considered, where a fluid flow is initialized with two constant (left and right) states. Since in absence of chemistry, this case does not allow us to validate the transport of species (which behave as passive scalars), it is proposed here to extend it to the study of two ideal gases with different heat capacity ratios γ, as suggested by Abgrall and Karni [START_REF] Abgrall | Computations of Compressible Multifluids[END_REF]. A (L = 1)-long domain, discretized with N = 400 points, is initialized with left and right states separated by a discontinuity localized at x = 0.5 as

(ρ, u, p, γ,Y ) L = (1, 0, 1, 1.4, 0), ( 42
) (ρ, u, p, γ,Y ) R = (0.125, 0, 0.1, 1.6, 1). ( 43
)
The varying heat capacity ratio is considered with a mass fraction model through the following relation

γ(Y ) = (1 -Y )γ L /(γ L -1) +Y γ R /(γ R -1) (1 -Y )/(γ L -1) +Y /(γ R -1) . ( 44
)
Thanks to this variable heat capacity ratio, the mass fraction Y is not a passive scalar in absence of chemistry, which allows us to validate the behavior of the whole system in presence of flow discontinuities. The gas constant is R = 1 and the time step is such that ∆x/∆t = 5, leading to max(CFL) ≈ 0.45. Left and right boundaries are imposed as Neumann conditions, i.e. a zero-gradient for every variable. Figure 4 displays the density and temperature profiles at t = 0.2 √ RT L /L obtained with the four following schemes: (1) MHM with s c = 0, (2) Heun (TVD) with s c = 0, (3) MHM with s c = 1, (4) Heun (TVD) with s c = 1. The first observation is that the jump relations observed between the rarefaction wave, the contact discontinuity and the shock wave are in very good agreement with the reference provided by an exact Riemann solver [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. This allows us to validate the conservativity of the coupled scheme. The main difference between the four models lies in the appearance of numerical oscillations of the temperature profile close to flow discontinuities. More precisely, for s c = 0, the MHM model gives rise to oscillations in the shock and in the contact discontinuity. The latter are removed by the use of the Heun (TVD) scheme. However, although theoretically TVD, numerical oscillations are still present in the shock with this scheme. This observation is correlated with a particular feature of the proposed hybrid schemes, for which entropy waves are decoupled from acoustic ones. As previously shown in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], the contact discontinuity is related to an entropy and mass fraction wave, which is here numerically solved by a TVD scheme. This is why no oscillations are observed in the contact discontinuity. However, the shock wave is related to pressure fluctuations, which are purely solved by the LBM. For this reason, the only way to damp its oscillations is to apply modifications to the LB scheme itself. In Fig. 4, the LB scheme is modified by adding an artificial viscosity through the shock sensor with s c = 1. This results in the successful suppression of oscillations in the shock, regardless of the finite volume scheme. The decoupled behavior of the contact and shock discontinuities are summarized in Table 2.

MHM

Heun (TVD) Given these conclusions, and since non-oscillating flow discontinuities are of primary importance for the numerical simulation of detonations, it is decided to retain the Heun (TVD) scheme only for the rest of this work. This choice may seem suprising in view of the dissipative behavior of this scheme exhibited in Sec. 3.1. In fact, for most of the applications involving inert mixtures, the MHM2D scheme is likely to outperform the Heun (TVD) one, precisely thanks to its low dissipation in every direction. However, the numerical oscillations observed across contact discontinuities turned out to be problematic for the simulations of detonations, where the flow physics is strongly driven by the discontinuities. Temperature oscillations can result in spurious chemical reactions that pollute the solution. Therefore, even though the Heun (TVD) scheme can be more dissipative than the MHM2D one, it proved to be a very good scheme for the studies that follow.

s c = 0 contact: ✗ shock: ✗ contact: ✓ shock: ✗ s c = 1 contact: ✗ shock: ✓ contact: ✓ shock: ✓

Stable one-dimensional detonation

The third test case aims at validating the ability of the proposed method to simulate a single detonation wave. For this purpose, the one-dimensional steady model initially developed by Zel'dovich, von Neumann and Döring (ZND) [START_REF] Zeldovich | On the theory of the propagation of detonation in gaseous systems[END_REF][START_REF] Neuman | Theory of Detonation Waves[END_REF][START_REF] Döring | Über den Detonationsvorgang in Gasen[END_REF] will be considered as a reference. The principle of this approach is to numerically integrate the steady one-dimensional conservation equations in order to have access to the structure of the detonation wave, starting from the knowledge of the physical properties of the mixture and the fresh gas steady state. The detonation is then especially characterized by a half reaction length l 1/2 corresponding to the distance between the shock and a fresh gas mass fraction Y = 0.5, by a von Neumann state corresponding to a pressure peak behind the shock and by the Chapman-Jouguet detonation velocity D CJ . When the mixture is composed of two ideal gases with an identical heat capacity ratio γ and with a heat release Q, the detonation velocity reads

D CJ = γRT 0 + (γ 2 -1)Q 2 + (γ 2 -1)Q 2 , (45)
T 0 being the temperature of the fresh gas [START_REF] Zhang | Shock Waves Science and Technology Library[END_REF]. In the present case, such a mixture will be considered with the following properties:

γ = 1.2, R = 692.83 J.kg -1 .K -1 , Q = 50RT 0 , E a = 24RT 0 , k = 2.006 × 10 7 s -1 , p 0 = 101325 Pa, T 0 = 295 K. ( 46
)
With these parameters, note that the half reaction length of the ZND profile is l 1/2 = 6.92 × 10 -4 m. This particular test case is commonly investigated in the literature for the validation of numerical methods [START_REF] Zhang | Shock Waves Science and Technology Library[END_REF][START_REF] He | The dynamical limit of one-dimensional detonations[END_REF][START_REF] Sharpe | Numerical simulations of pulsating detonations: I. Nonlinear stability of steady detonations[END_REF][START_REF] Kasimov | On the dynamics of self-sustained one-dimensional detonations: A numerical study in the shock-attached frame[END_REF][START_REF] Ng | Nonlinear dynamics and chaos analysis of one-dimensional pulsating detonations[END_REF][START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF]. The reason is that thanks to linear stability analyses, it can be theoretically demonstrated that the ZND profile is stable under these conditions [START_REF] Lee | Calculation of linear detonation instability: one-dimensional instability of plane detonation[END_REF][START_REF] Sharpe | Linear stability of idealized detonations[END_REF][START_REF] Short | Cellular detonation stability. Part 1. A normal-mode linear analysis[END_REF]. Therefore, the expected solution over time is a mere propagation of the ZND structure at the detonation velocity D CJ . In this section, the density, velocity, temperature, and mass fraction profiles obtained by the ZND theory are considered for the initialization of a detonation wave propagated to the right. Compared to other possible ways of initiating a detonation, the use of the ZND solution offers several advantages: one can make sure that a nonoverdriven detonation takes place and the number of points in l 1/2 , referred to as N 1/2 , can be explicitly prescribed. The simulation setup is then the following. A one-dimensional domain of L = 0.025 m is considered, where Neumann conditions are imposed on the left and right boundaries. In order to make the detonation wave freely propagate over long distances (larger than the domain length L), a stategy referred to as recycling is adopted: when the detonation wave is detected close to the right boundary, the solution is shifted to the left and fresh gas is filled between the shock and the boundary [START_REF] Sow | On the viscous boundary layer of weakly unstable detonations in narrow channels[END_REF]. Furthermore, in order to avoid spurious reflexions due to imperfect absorption properties of the Neumann conditions, a viscosity sponge zone is introduced close to the left boundary. It consists in a continuous increase of the LB relaxation time τ so as to reach τ = ∆t/2 close to the boundary, which, as discussed in Sec. 2.2, is a particular value leading to a first-order LBM filtering any non-hydrodynamic modes. An inviscid fluid being considered, the relaxation time is then computed as a superposition of the sponge zone with the artificial viscosity of the shock sensor as

τ = min ν sc c 2 s + 1 4 (1 -tanh ((x -s x )/s w )), ∆t 2 , ( 47
)
where s x = 0.01 m is the x-position of the sponge zone and s w = 0.001 m its characteristic width. These parameters are adopted so that the sponge zone always remains far away (more than 10l 1/2 ) from the reaction zone. Furthermore, the ratio between the mesh size and the time step is set to ∆x/∆t = 2 × 10 4 , leading to max(CFL) ≈ 0.2. Figure 5 displays the dimensionless pressure and temperature profiles obtained at time t = 4000t c , where t c = l 1/2 / √ RT 0 , for a mesh resolution N 1/2 = 32 in two cases: s c = 0 (no shock sensor) and s c = 0.1. Note that the detonation wave is represented here as propagated to the left direction, which is a conventional representation of ZND profiles. The pressure is non-dimensionalized by the von Neumann value p N ≈ 4.26 × 10 6 Pa. One can notice a very good agreement between the profiles obtained by the HLBM approach with and without shock sensor, and the ZND theory. The reaction zone is accurately captured, indicating the ability of the model to reproduce the complex coupling occurring between the shock wave and the chemical reaction. Only two differences are noticeable between the s c = 0 and the s c = 0.1 cases:

• For s c = 0, a pressure peak is observed in the shock, leading to an incorrect von Neumann pressure.

• For s c = 0.1, the shock discontinuity is slightly smoother than without shock sensor.

The numerical oscillations observed in absence of shock sensor are likely to be problematic for the forthcoming numerical simulations of detonations, especially since the pressure behind the shock is one of the most important observed quantities. For this reason, the use of the shock sensor with the parameter s c = 0.1 will be retained in the following. This is all the more motivated by the fact that, as shown in Fig. 5, the sensor only increases the relaxation time (hence the viscosity) in the shock region, while τ = 0 is recovered everywhere else and especially in the reaction zone. Thus, this strategy effectively reduces the numerical oscillations without increasing the viscosity away from the shock. 45)).

Resolution

In the case s c = 0.1, the measured detonation velocities are displayed for several meshes and compared to the Chapman-Jouguet theory (Eq. ( 45)) in Table 3. Note that the exactly same velocity D is obtained for N 1/2 = 16 and N 1/2 = 32, which is due to the fact that, over the considered time, the detonation wave can only flow through an integer number of mesh sizes ∆x. In any case, one can notice a very good agreement between the HLBM simulations and the expected solution with a relative error in the order of 10 -6 . This result is highly accurate compared to the numerical simulations based on the Euler equations of Kasimov et al. [START_REF] Kasimov | On the dynamics of self-sustained one-dimensional detonations: A numerical study in the shock-attached frame[END_REF], who obtained a relative error in the order of 10 -3 with a finer resolution (N 1/2 = 200). According to Henrick et al. [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF], their lack of convergence was due to a first-order accuracy of the numerical scheme. More accurate results could be obtained with a true fifth-order shock-fitting scheme, leading to a relative error in the order of 10 -7 for N 1/2 = 20. With regards to these results, the hybrid LBM is competitive with high-order standard numerical schemes for the simulation of one-dimensional stable detonation waves.

Physical instabilities of inviscid gaseous detonations

Given the results of the previous section on academic test cases, only a single model is retained here: the Heun (TVD) scheme of Sec. 2.3 with a shock sensor parameter s c = 0.1, unless otherwise stated. In this section, the physical instabilities of one-and two-dimensional gaseous detonations are investigated.

One-dimensional detonations

According to linear stability analyses, the one-dimensional detonation investigated in Sec. 3.3 becomes unstable when the activation energy is increased above a critical value E a ≳ 25.3RT 0 , all others parameters being unchanged [START_REF] Zhang | Shock Waves Science and Technology Library[END_REF][START_REF] Lee | Calculation of linear detonation instability: one-dimensional instability of plane detonation[END_REF][START_REF] Sharpe | Linear stability of idealized detonations[END_REF][START_REF] Short | Cellular detonation stability. Part 1. A normal-mode linear analysis[END_REF]. This means that the ZND profile is not the unsteady solution any more. This topic has been the purpose of many non-linear analyses based on unsteady numerical simulations [7, 79-83, 88, 89]. In this section, it is proposed to reproduce the study of Ng [START_REF] Ng | Nonlinear dynamics and chaos analysis of one-dimensional pulsating detonations[END_REF] who developped a finite-volume high-order SLIC scheme to investigate the non-linear behavior of four unsteady detonations characterized by different activation energies. In order to preserve the same half reaction length as in Sec. 3.3, i.e. l 1/2 = 6.92 × 10 -4 m, the exponential pre-factor is adapted for each case. Compared to Eq. ( 46), the simulation parameters are modified as: These cases are especially selected since they are expected to give birth to different instability modes according to the bifurcation diagram [START_REF] Zhang | Shock Waves Science and Technology Library[END_REF][START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF]. Except the changes in the activation energy and the exponential pre-factor, the setup of the numerical simulations is the same as that described in Sec. 3.3. Notably, the shock sensor parameter is set to s c = 0.1 and we set ∆x/∆t = 2 × 10 4 , which allows us to have a maximal CFL number lower than 0.3 for all the cases. Figure 6 displays the leading shock pressure histories as well as the corresponding phase plots obtained for each case. The results are in very good agreement with Ng [START_REF] Ng | Nonlinear dynamics and chaos analysis of one-dimensional pulsating detonations[END_REF] of the flow behavior with the activation energy. In particular, Cases 1, 2, 3, 4 respectively lead to 1, 2, 4, 8 modes of instabilities. It is also noteworthy to point out the long simulation time required to reach a non-linear limit cycle with a converged structure. This observation is in agreement with the conclusions of Sharpe and Falle [START_REF] Sharpe | Numerical simulations of pulsating detonations: I. Nonlinear stability of steady detonations[END_REF] who run simulations for thousands of half-reaction times to reach a convergence. Moreover, by performing a mesh convergence study for each case, we notice that when the mesh is coarse (N 1/2 = 16), the higher the activation energy and the more sensitive the computation is to small perturbations of the flow. Therefore, it is sometimes difficult to obtain a temporal convergence of the simulation and instability bifurcations from one mode to another can be observed on a single computation. However, this behavior vanishes when the mesh is sufficiently refined, so that a converged structure can be finally observed in any case. It is also noticed that, even when a limit cycle can be obtained, the number of modes appearing may not be the expected one. Here again, refining the mesh finally allows reaching the correct unsteady structure. This mesh convergence is illustrated on Fig. 7 where phase plots are displayed for several meshes in each case. One may especially note that the eight unstable modes expected by Case 4 can be obtained with no less than N 1/2 = 512 points. Also note that when E a increases, a finer mesh may be required to converge in the pressure amplitude even when the expected mode is reached. All these observations explain why each case of Fig. 6 is presented with a different mesh resolution N 1/2 .

These observations are in agreement with previous numerical simulations who concluded that a fine mesh, with more than 100 points in l 1/2 , may be required to observe an expected instability when the activation energy increases [START_REF] Sharpe | Numerical simulations of pulsating detonations: I. Nonlinear stability of steady detonations[END_REF]. It can be explained by two phenomena. First, due to the strong oscillations observed in the detonation front, the effective number of points in l 1/2 varies with time so that, in practice, the reaction zone may be less resolved than expected with a steady ZND profile. Secondly, one should mention that the unsteadiness is a result of acoustic waves propagating in the reaction zone, between the shock discontinuity and the Chapman-Jouguet state. Their wavelength is likely to be smaller as the number of modes increases, thus requiring a finer mesh to be captured. The mesh convergence analysis can be used to accurately estimate the time period of the oscillations in each case. The procedure adopted here is as follows. For each activation energy, the period of the pressure amplitude is measured with different meshes, ranging from N 1/2 = 16 to N 1/2 = 512 points. Then, noticing that a first-order convergence is systematically obtained, a linear regression is performed with an inverse function. This makes it possible to obtain estimations of the time period in the limit of an infinitely fine mesh, referred to as (Period) l . The results are summarized in Table 4. In the first case, the measured period is very close to that of Henrick et al. [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF] obtained with a fifth-order accurate scheme.

E a (
The convergence analysis displayed in Fig. 8 confirms the first-order convergence of the measured periods toward the asymptotic limits given in Table 4. This first-order accuracy may seem surprising while the LBM, which is secondorder accurate, is coupled with a MUSCL scheme. In fact, two properties of the adopted numerical scheme can explain this effective first-order accuracy: 1) the implementation of the chemical source terms using a splitting approach, 2) the add of artificial viscosity through the use of a shock sensor. These conclusions can be retained for future work aiming to increase the effective accuracy order of the scheme.

Although it is clearly first-order accurate, the proposed HLBM model remains competitive with regards to the numerical schemes commonly encountered in the literature, even when compared with high-order methods. The detonation instabilities can be accurately captured with mesh resolutions of the same magnitude as in previous work, and the detonation velocity could be measured with great accuracy in Sec. 3.3. The good accuracy obtained with this HLBM scheme has two explanations. First, the artificial viscosity introduced by the shock sensor, which is one source of first-order accuracy, is only triggered in a very thin zone in the shock region as shown in Fig. 5. Knowing that the oscillations of the shock wave result from physical modes propagating in the reaction zone, it is worth noting that this region is not, or very slightly, affected by the shock sensor. The second reason lies in an important feature of the proposed HLBM solver together with the physical nature of the detonation instabilities. Precisely, one should mention that the physically unstable modes propagating in the reaction zone are acoustic waves. One the other hand, the proposed hybrid scheme has been especially designed to ensure a linear decoupling of the flow characteristics. As such, it exactly behaves as if acoustic phenomena were addressed by the LBM only, as shown by Wissocq et al. [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF]. The very good numerical properties of the LBM for acoustic waves propagation, which is comparable to high-order Navier-Stokes schemes [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF], are therefore preserved. This is why the physical modes responsible for the oscillations of the detonation wave can be accurately captured with a moderate number of points per half reaction zone, further explaining the very good accuracy of this approach for this case. This discussion highlights an important feature of the proposed hybrid method, preserving the main interest of the LBM even when coupled with finite volume schemes.

Time-averaged detonation velocities obtained for each activation energy with the finer mesh (N 1/2 = 512 points) are displayed in Figure 9 in comparison with D CJ . Here again, a very good agreement is observed with the numerical results of Henrick et al. [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF], further validating the ability of the HLBM to accurately simulate unsteady detonations.

It is also worth noting that the relative deviations with the steady Chapman-Jouguet theory is of very small order of magnitude (∼ 10 -4 ). Having an accurate estimation of the detonation velocity with more than five digits is therefore crucial to obtain this result. Notably, the first-order numerical scheme of Kasimov et al. [START_REF] Kasimov | On the dynamics of self-sustained one-dimensional detonations: A numerical study in the shock-attached frame[END_REF] could not provide sufficiently reliable measurements. Similarly, the detonation velocities measured by Ng [START_REF] Ng | The effect of chemical reaction kinetics on the structure of gaseous detonations[END_REF], showing a chaotic evolution with E a and relative deviations with D CJ larger than 10 -3 , were likely to be not accurate enough. To the authors' knowledge, only the fifth-order method proposed by Henrick et al. [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF] could provide satisfactory results so far. In this regard, the confrontation carried out in Fig. 9 is not only a numerical validation of the present HLBM model, but also a confirmation of the physical trend observed in [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF].

A final discussion can be raised regarding the role played by the numerical parameters s c and ∆x/∆t in this study. We recall that s c controls the add of artificial viscosity in the shock region. Since a clipping is done in Eq. ( 26), above a given value, s c mainly acts as a stiffness parameter for the switch between a second-order scheme and a first-order one. Regarding the ratio ∆x/∆t, it is closely related to the CFL number, playing an important role on the numerical stability of the method. To better exhibit their impact, case 2 is reproduced with N 1/2 = 64 and three different choices for these parameters:

• s c = 0, ∆x/∆t = 2 × 10 4 ,

• s c = 0, ∆x/∆t = 10 4 ,

• s c = 0.1, ∆x/∆t = 10 4 , and compared with the results provided in Fig. 6. Note that setting ∆x/∆t = 10 4 leads to an increase of the CFL number up to max(CFL) ≈ 0.48, which nonetheless did not harm the numerical stability of the simulations. Fig. 10 displays the time evolution of the pressure peaks (left) and the pressure profiles measured when the maximal amplitude is reached (right) in each case. This figure illustrates the fact that removing the artificial viscosity in the shock (s c = 0) and increasing the CFL number have a similar effect on the pressure profile: important oscillations can be generated close to the discontinuity, which, however, has very little impact on the reaction zone. As expected, two modes of instabilities are observed in any case, but the numerical oscillations can lead to a large overestimation of their amplitude, especially in the worst-case scenario tested here (s c = 0, ∆x/∆t = 10 4 ). This observation motivates the choices adopted in this section: s c = 0.1, ∆x/∆t = 2 × 10 4 leading to max(CFL) ≤ 0.3 and no overestimation of the pressure peak. ) simulated different values the shock sensor parameter s c and the ratio ∆x/∆t, with N 1/2 = 64 points. Note that ∆x/∆t = 2 × 10 4 leads to max(CFL) ≈ 0.23, while ∆x/∆t = 4 leads to max(CFL) ≈ 0.48. Left: leading shock pressure history, right: pressure profile when the maximal pressure amplitude is reached.

Two-dimensional detonations

In this section, two-dimensional simulations are performed in order to assess the ability of the proposed numerical method to accurately capture the physical properties of multidimensional detonations. We consider the thermodynamics parameters from Reynaud et al. [START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF][START_REF] Reynaud | Computation of the mean hydrodynamic structure of gaseous detonations with losses[END_REF], where a stochiometric H 2 -O 2 mixture at atmospheric pressure p 0 = 101.25 kPa and temperature T 0 = 295 K was investigated. A single-step chemical model is used with a reduced heat release Q/RT 0 = 28.3. The ratio of specific heats is γ = 1.333 to accurately recover the post-choc properties and the CJ speed provided in Table 5. The mixture sensitivity is varied by considering different reduced activation energies E a /RT 0 = 20, 30, 38.23, 48, so that a wide range of detonation instabilities can be covered. As in Sec. 4.1, the preexponential factor is modified accordingly in order to keep the same half-reaction length for all the reduced activation energies, in this case l 1/2 = 9.079 µm. This length scale is calibrated in order to fit the experimental two-dimensional cell size. The pre-exponential factors and the stability criteria χ for the different activation energies are gathered in Table 6. The factor χ is computed from the definition of Ng et al. [START_REF] Ng | Numerical investigation of the instability for one-dimensional chapman-jouguet detonations with chain-branching kinetics[END_REF] adapted to a single-step chemistry. In this study, the authors proceeded to a classification of various mixtures according to their sensitivity by experimental observations of the cellular structure. A limit allowing them to distinguish stable from unstable mixtures was found for χ = 10. Beyond this critical value, detonations are considered unstable. Based on the value obtained from this analysis, mixtures with the activation energies E a /RT 0 = 38.23 and 48 are expected to be unstable. To further attest this observation, the reduced activation energies as function of the CJ Mach number are compared in Fig. 11, where a dashed line indicates the stability limit [START_REF] Eckett | Numerical and analytical studies of the dynamics of gaseous detonations[END_REF]. Experimental H 2 -O 2 mixtures are displayed with triangles for comparison with different degrees of dilution with Argon. Stable detonation are expected for E a /RT 0 = 20 since it is below the stability limit while E a /RT 0 = 30 is just above it. On the other hand, E a /RT 0 = 38.23 and E a /RT 0 = 48 are in the unstable zone. Experimental H 2 -O 2 mixtures show that dilution makes it possible to switch from an unstable zone to the stable one. E a /RT 0 = 48 is close to the representative points of the CH 4 -2O 2 and C 3 H 8 -5O 2 mixtures, which can be considered as "very" unstable ones with very irregular cell structures. Simulating such mixtures would allow us to cover the large area of physical properties and multidimensional dynamics that can be met in detonation applications.

Two-dimensional simulations have been performed to check if the physics of the detonation can be correctly reproduced for each activation energy. Analyzing the flow pattern allows us to estimate the size of the cells and their regularity. A Cartesian computational domain of dimensions L x × L y = 15 mm ×5 mm = 165l 1/2 × 55l 1/2 is considered. The domain is chosen so that enough cellular structures can be exhibited to investigate their regularity. Upstream and downstream boundary conditions are similar as in Sec. 3.3: Neumann conditions are considered together with a viscosity sponge zone on the left boundary. Top and bottom boundaries are slip walls imposed by ensuring a symmetry thanks to the use of ghost nodes. The numerical resolution is ∆x = l 1/2 /16 which is expected to be sufficient to capture main physical patterns of the gas dynamics [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF].

E a /RT 0 k (s -
The detonation is initiated by setting a ZND profile at the left-hand-side of the domain, enabling a propagation into a fresh reactive mixture in the right-hand-side. Similarly to the one-dimensional cases, when the detonation reaches the right boundary, a recycling of the solution is performed by adding extra fresh reactive mixture downstream of the shock. More detailed description of the recycling methodology can be found in Sow et al. [START_REF] Sow | On the viscous boundary layer of weakly unstable detonations in narrow channels[END_REF].

Figure 12 displays the resulting cellular structure for all the activation energies once a steady state propagation is reached. An excellent cellular regularity is observed for the case of E a /RT 0 = 20. With the increase of the activation energy, the cellular structure tends to loose its regularity and becomes highly irregular for E a /RT 0 = 48. Following some classifications of the literature [START_REF] Libouton | Cinétique, structure et entretien des ondes de détonation[END_REF][START_REF] Strehlow | The structure of marginal detonation waves[END_REF], structures of excellent regularity are encountered in highly diluted detonations with monatomic gases such as argon and helium. The associated reactive mixtures are generally light fuels (H 2 , C 2 H 2 ...) reacting with oxygen. Structures of good regularity are recovered with E a /RT 0 = 30. They are generally found in low dilution and high pressure mixtures. The size of the cells is relatively constant and can be measured with a reasonable accuracy. By increasing the activation energy to E a /RT 0 = 38.23, the resulting structures show a rather large distribution in cell sizes. Highly irregular structures as with E a /RT 0 = 48 are usually the signature of a heavy fuel reacting with O 2 or CH 4 reacting with air. This irregularity can be related to the properties of the transverse waves. Indeed, for stable mixtures, transverse waves are very weak, while in unstable ones, transverse waves are strong and can, under certain conditions, give birth to transverse detonations [START_REF] Gamezo | Fine cellular structures produced by marginal detonations[END_REF][START_REF] Kessler | Multilevel detonation cell structures in methane-air mixtures[END_REF][START_REF] Taylor | Numerical simulations of hydrogen detonations with detailed chemical kinetics[END_REF].

The general characteristics of the detonation are then revealed through an analysis of instantaneous fields. For simplicity, a qualitative comparison is drawn between two activation energies: E a /RT 0 = 20 and E a /RT 0 = 48. In this Increasing the activation energy results into a more complex detonation front. Figure 14 displays the resulting temperature flow field. Flow symmetry is broken and Mach stem/incident shock are difficult to distinguish. Unburned pockets are detached from the front and convected downstream into the burnt gas. This consumption mechanism is due to several phenomena such as surface burning caused by instabilities at the interface with burnt gas and volume burning by the global increase of the temperature inside the pocket [START_REF] Radulescu | A detonation paradox: Why inviscid detonation simulations predict the incorrect trend for the role of instability in gaseous cellular detonations?[END_REF]. Transverse waves are weaker than with E a = 20RT 0 [START_REF] Lee | The detonation phenomenon[END_REF].

To further highlight the effect of the shock sensor on the flow, the dimensionless relaxation time τ/∆t = 1/2+τ/∆t is displayed in Fig. 15. As expected, the relaxation time is increased in the shock region, corresponding to the add of artificial viscosity. As discussed in Sec. 2.2, the particular value τ = ∆t/2 (or τ = ∆t) reached in the shock can be interpreted as the local use of a first-order accurate scheme allowing us to damp the numerical oscillations of the LBM close to discontinuities. Away from discontinuities, a value of τ = 0 (τ = ∆t/2), is recovered in agreement with an inviscid description of the flow. For stable detonation, the artificial relaxation remains activated in the transverse shock wave downstream of the triple point over a long distance. On the contrary, it vanishes rapidly in the case of high activation energy case due to the strong instability of the transverse waves. This result confirms the previous observation for high activation energies.

Conclusion

In this work, new finite-volume schemes have been proposed for the simulation of compressible reactive flows with a hybrid LBM solver in the aim of numerically modelling inviscid gaseous detonations. For this purpose, the methodology initially proposed in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF] has been extended to multi-species reactive flows. The key point in the strategy is to end up with a fully conservative numerical scheme that is linearly equivalent to a non-conservative one, where the LBM is supplemented by the characteristic equations of entropy s and mass fractions Y k . As a consequence, the Kovásznay modes remain linearly decoupled from each other, preserving the following advantages:

• the low dissipation of the LBM for isentropic phenomena such as pressure waves is preserved,

• entropy and mass fraction transport can be explicitly controlled, in terms of stability and dissipation, by the choice of an appropriate finite-volume scheme.

Thanks to this strategy, stable and conservative numerical discretizations of the total energy equation have been derived, which a challenging task when the hybrid LBM is considered [START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF][START_REF] Guo | On the use of conservative formulation of energy equation in hybrid compressible lattice Boltzmann method[END_REF]. The introduction of production terms modelled by the Arrhenius law and numerically integrated by a splitting approach has then made it possible to consider the chemical reactions involved in gaseous detonations. Subsequently, three numerical schemes have been proposed for the the total energy and mass fraction equations. The first one is directly derived from the MUSCL-Hancock method previously adopted by Farag et al. [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF] in their entropy-based hybrid LBM. The second one is an extension of it for non-aligned flows, considering diagonal points in the spatial stencil as recently proposed by Yoo et al. [START_REF] Yoo | Compressible lattice boltzmann method with rotating overset grids[END_REF]. The third scheme is based on a Heun version of the MUSCL-Hancock method supplemented by a slope limiter in order to ensure TVD properties. A first numerical validation based on standard academic cases (transport of inert species, Sod shock tube and stable one-dimensional detonation) allowed us to draw comparisons between these numerical schemes. When supplemented with artificial viscosity to address numerical oscillations of the shock wave, the Heun (TVD) scheme turned out to provide the best results for compressible cases in presence of flow discontinuities. For this reason, and despite its rather dissipative behavior exhibited in the transport of inert species, this scheme was retained for the numerical simulations of gaseous detonations performed in the present work. The reader should however note that this conclusion may not remain true in all cases. It is very likely that, as long as numerical oscillations remain admissible close to flow discontinuities, the second scheme would lead to the best results. Regarding this point, the development of a second-order low-dissipative and non-oscillating scheme for the total energy equation may be the purpose of future work.

The retained methodology led to very satisfying results for one-and two-dimensional gaseous detonations. In particular, the detonation velocity as well as the time period of unsteady structures could be recovered with more precision than standard low-order methods, for example that of Kasimov et al. [START_REF] Kasimov | On the dynamics of self-sustained one-dimensional detonations: A numerical study in the shock-attached frame[END_REF]. The observed accuracy is not so far from the one obtained by Henrick et al. [START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF] with a fifth-order numerical scheme. This interesting property can be explained by the low dissipation of the LBM for the acoustics [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF], which is preserved with the present hybrid model. It is indeed known that most of the detonation instabilities are consequent of pressure waves developing in the reaction zone, which can therefore be accurately simulated with the proposed method.

These promising results are an important step toward the use of the hybrid LBM for industrial applications involving compressible reactive flows. It is worth noting that the primary focus of the present study is centered on the formulation and validation of numerical fluxes. This is why only the most simple introduction of chemical source term, based on a splitting approach, was considered. Since the proposed method is not dependent on the specific integration of the source term, it can readily be replaced by any numerical scheme. For instance, the use of a predictor-corrector scheme can be more adapted to handle stiff chemistry [START_REF] Mott | Chemeq2: A solver for the stiff ordinary differential equations of chemical kinetics[END_REF][START_REF] Blanchard | Chemical modeling for methane oxy-combustion in liquid rocket engines[END_REF]. Also note that while the present study is restricted to inviscid detonations, the introduction of nonzero dynamic viscosity and thermal diffusion is straightforward with the proposed method, as shown in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF]. Its extension to unsteady detonations in the viscid case, as done in [START_REF] Sow | Mean structure of one-dimensional unstable detonations with friction[END_REF][START_REF] Romick | The effect of diffusion on the dynamics of unsteady detonations[END_REF], would require the introduction of a diffusion term in the mass fraction equation, which may be carried out in further studies. The validation of this numerical method for more realistic simulations including rotating detonation engines can be the purpose of future work. Its extension to even more challenging compressible flows including high pressure ratio may require further developments of: 1) more robust LBM schemes for higher Mach numbers, 2) new strategies to deal with the flow discontinuities avoiding the need for a shock sensor. Last but not least, an in-depth investigation of the competitiveness of this method with regard to finite volume schemes could be conducted in the future.

Acknowledgement

Alois Sengissen and Gabriel Farag are gratefully acknowledged for fruitful discussions. Centre de Calcul Intensif d'Aix-Marseille and GENCI (Grant 2022 -A0112B11951) are acknowledged for granting access to their high performance computing resources. This research was supported by the ANR projects ALBUMS (grant ANR-18-CHIN-0003-01) and MALBEC (grant ANR-20-CE05-0009).

Appendix A. Jacobian matrix between conserved and non conserved variables

The aim of this appendix is to provide an expression for the Jacobian matrix M = ∂ U/∂ V which remains valid for a large variety of gas mixtures, where U = [ρ, ρu, ρE, ρY k ] T and V = [ρ, ρu, s,Y k ] T . To build this matrix, let us use the first law of thermodynamics, de = T ds + p

ρ 2 dρ + N s ∑ k=1 g k dY k = T ds + p ρ 2 dρ + N s ∑ k=2 (g k -g 1 )dY k , (A.1)
where g k = h k -T s k is the Gibbs free energy and where s k is the entropy of the specy k. We assume here that the entropy of the mixture can be obtained by s = ∑Yksk. The aim of this appendix is to provide the mathematical details for the construction of the conservative scheme obtained in Eq. ( 14). To this end, starting from Eq. ( 13), we only have to apply modifications to the non-conservative transport term MA α δ α F V,d α . This is done here by ensuring that an equivalence remains true in the linear approximation. We therefore define the operator = representing a linear equality between left and right hand sides of an equation.

The first step of the present demonstration starts by noticing that, away from flow discontinuities, a strict equivalence exists between the continuous non-conservative system

∂ V ∂t + A α ∂ F V,c
MA α δ α K α (F V,c α ) =δ α K α (F U,c α ), (B.5)
which is one of the key relations that can be used to transform the transport term MA α δ α F V,d α in Eq. ( 13). We then define the projection matrices P LB , P E and P Y such that

P LB U =          ρ ρu 0 0 . . . 0          , P E U =          0 0 ρE 0 . . . 0          , P Y U =          0 0 0 ρY 2 . . . ρY N s          . (B.6)
Note that necessarily, P LB + P E + P Y = I. Then, one has, thanks to Eq. (B.5) and for two arbitrary functions K 1,α and K 2,α , the following linear equivalence: the linear equivalence can be written as

MA α δ α F V,d α =MA α δ α F V,d α + P E δ α K 1,α (F U,c α ) -MA α δ α K 1,α (F V,c α ) + P Y δ α K 2,α (F U,c α ) -MA α δ α K 2,α (F V,
MA α δ α F V,d α =P LB MA α δ α F V,d α + P E δ α K 1,α (F U,c α ) + MA α δ α F V,d α -K 1,α (F V,c α ) + P Y δ α K 2,α (F U,c α ) + MA α δ α F V,d α -K 2,α (F V,c α ) . (B.9)
This relation can then be transformed noticing two facts. The first thing is that, since both the original and conservative operators use LBM for mass and momentum transport, the matrix MA β has an identity block corresponding to the mass and momentum variables. This reads

P LB MA β = P LB . (B.10)
The second thing is that, due to the consistency of the discrete scheme and the operator K 1,α , we have

F V,d α -K 1,α (F V,c α ) = O(∆x, ∆t), (B.11)
and similarly for K 2,α . As a consequence, the MA α pre-factors can be included inside the operator δ α without modifying the linear equivalence. This leads to MA α δ α F V,d α = P LB δ α F V,d α + P E δ α K 1,α (F U,c α ) + MA α F V,d α -K 1,α (F V,c α )

+ P Y δ α K 2,α (F U,c α ) + MA α F V,d α -K 2,α (F V,c α ) . (B.12)

Finally, one can build the conservative scheme of Eq. ( 14) with

F U,d α = P LB F V,d α + P E K 1,α (F U,c α ) + MA α F V,d α -K 1,α (F V,c α ) + P Y K 2,α (F U,c α ) + MA α F V,d α -K 2,α (F V,c α ) . (B.13)
We note that the mass and momentum fluxes are unchanged from F V,d α . Therefore, the discrete conserved fluxes can be written as F U,d α = F To illustrate this, the choice of the Jacobian matrix M of Eq. (A.5) leads to In what follows, this model will be summarized with ζ = 0 and κ = 0. Regarding the collision operator, it can be viewed as the hybrid recursive regularized (HRR) model of Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF] with a hybridization parameter σ = 1 (no finite differences in the regularization procedure). It reads where f eq i is the equilibrium distribution function, f neq i is the regularized off-equilibrium part, F E i is a body-force term designed to restore the Galilean invariance at the Navier-Stokes level and τ is the relaxation time. For the D2Q9 lattice, the equilibrium distribution function reads [START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] f eq i = w i ρ 1 + where θ = RT /c 2 s , c s is the so-called lattice sound speed related to the mesh size and time step (acoustic scaling [START_REF] Krüger | The Lattice Boltzmann Method[END_REF]) as

f coll i (x,
c s = 1 √ 3 ∆x ∆t , (D.3)
and w i = [4/9, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36] are the Gaussian weights of the lattice [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]. The offequilibrium distribution functions are computed thanks to a recursive regularization [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF] as Finally, as proposed in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF], the correction term F E i reads

f neq i = w i (c i,α c i,β -c 2 s δ αβ ) 2c 4
F E i = w i 2c 4 s (c i,α c i,β -c 2 s δ αβ ) a F, (2)
αβ , (D.8)

where, for the D2Q9 lattice, a F,(2)

αβ = c 2 s δ αβ ρ ∂ u γ ∂ γ - ∂ (ρ(1 -θ )) ∂t + c 2 s u α ∂ (ρ(1 -θ )) ∂ β + u β ∂ (ρ(1 -θ )) ∂ α - ∂ (ρu 3 α ) ∂ α 3 δ αβ , (D.9)
where the first spatial derivative (velocity divergence) is estimated thanks to a second-order centered finite-difference scheme, the time derivative is discretized using a temporal upwind scheme and every other spatial derivatives are discretized using a first-order upwind scheme.

Figure 1 :

 1 Figure 1: Sketch of the spatial stencils of the finite volume schemes considered in this work. : node where the fluxes are computed, : points considered in the computation of the fluxes, : points not considered.

Figure 2 :

 2 Figure 2: Contours of the mass fraction Y for the convection of an inert species after 20 characteristic times. 25 contours are displayed for Y 1 ∈ [0, 1]. From left to right: MHM1D, MHM2D, Heun (TVD). Top: horizontal convection, bottom: diagonal convection.

Figure 3 :

 3 Figure 3: Convergence study of the specy advection along the horizontal and diagonal directions with the finite volumes schemes proposed in Sec. 2.3. Dotted line: -1 slope, dashed line: -2 slope, Dash-dotted line: -3 slope.

  is

Figure 4 :

 4 Figure 4: Density and temperature profiles at time t = 0.2 √ RT L /L for the shock tube test case with two γ. Solid line: reference, blue dashed line: HLBM simulations, from left to right: [MHM, s c = 0], [Heun (TVD), s c = 0], [MHM, s c = 1], [Heun (TVD), s c = 1].

Figure 5 :

 5 Figure 5: Pressure and temperature profiles of the stable one-dimensional detonation obtained at t = 4000t c for N 1/2 = 32 points with and without shock sensor.

Case 1 : 1 Case 2 : 1 , 3 : 1 , 4 :

 1121314 E a = 27.00 RT 0 , k = 3.230 × 10 7 s -E a = 27.40 RT 0 , k = 3.445 × 10 7 s -Case E a = 27.80 RT 0 , k = 3.672 × 10 7 s -Case E a = 27.82 RT 0 , k = 3.683 × 10 7 s -1 .

Figure 7 :

 7 Figure 7: Mesh convergence of the phase plots showing mode oscillations for different activation energies.

Figure 8 :

 8 Figure 8: Mesh convergence of the deviation obtained between the measured periods and the converged periods given by 4 for varying activation energies. Dotted line: -1 slope.

Figure 9 :

 9 Figure 9: Relative devation of the detonation velocity with D CJ as a function of the activation energy for the selected cases. Comparisons with Henrick et al. [83].

NFigure 10 :

 10 Figure 10: Detonation case 2 (E a = 27.40 RT 0 , k = 3.445 × 10 7 s -1 ) simulated different values the shock sensor parameter s c and the ratio ∆x/∆t, with N 1/2 = 64 points. Note that ∆x/∆t = 2 × 10 4 leads to max(CFL) ≈ 0.23, while ∆x/∆t = 4 leads to max(CFL) ≈ 0.48. Left: leading shock pressure history, right: pressure profile when the maximal pressure amplitude is reached.

DTable 5 :

 5 CJ P N (bar) P CJ (bar) T N (K) T CJ (Chapman-Jouguet and von Neumann parameters obtained for the thermo-chemical parameters retained in the two-dimensional study.

Figure 11 :

 11 Figure11: Representation of the reduced activation energies in the stability region. The stability limit is represented by dashed lines. The percentage shows the degree of dilution with Argon of the H2 -O2 mixture. From Ecket et al.[START_REF] Eckett | Numerical and analytical studies of the dynamics of gaseous detonations[END_REF].

Figure 12 :

 12 Figure 12: Detonation cellular structure for different reduced activation energies

Figure 13 :

 13 Figure 13: A sequence of instantaneous temperature fields showing the dynamics of unstable detonation for E a = 20RT 0 .

Figure 14 :

 14 Figure 14: A sequence of instantaneous temperature fields showing the dynamics of unstable detonation for E a = 48RT 0 .

Figure 15 :

 15 Figure 15: Dimensionless relaxation time τ/∆t = 0.5 + τ/∆t in the flow fields of regular (left) and irregular (right) detonations.

5 )

 5 The differential form of the total energy by unit of volume is given byd(ρE) = e dρ + ρde + d (ρu α ) α d(ρu α ) + ρT ds + ρ s,Y k = hκ, ∂ (ρE) ∂ (ρu α ) ρ,s,Y k = u α , ∂ (ρE) ∂ s ρ,ρu α ,Y k = ρT, ∂ (ρE) ∂Y k ρ,ρu α ,s = ρ(g kg 1 ), (A.3)where κ = ||u|| 2 /2 is the kinetic energy. Derivatives of the volume fraction can be obtained in a straightforward way,∂ (ρY ) ∂ ρ ρu α ,s,Y k = Y k , ∂ (ρY ) ∂ (ρu α ) ρ,s,Y k u x u y u z ρT ρ(g 2g 1 ) . . . ρ(g N sg 1 )Appendix B. Derivation of the conservative scheme

3 )

 3 α = ρu α , (ρu α u β + pδ αβ ) α∈{x,y,z} , ρHu α , (ρY k u α ) k∈ 2,N s T is the vector of conserved fluxes. Multiplying Eq. (B.1) by M = ∂ U/∂ V yields By identification with Eq. (B.2), it follows that MA α dF V,c α = dF U,c α . (B.4) Let us now consider an arbitrary function K α of a vector field (in space) Φ verifying K α (Φ) = Φ + O(∆x, ∆t). The operators proposed in Sec. 2.3 are particular examples verifying this property. Eq. (B.4) then yields

2 T.

 2 The last two fluxes (discrete total energy and mass fraction fluxes) are given by Eq. (B.13).

F ρE,FV +∆α/ 2 =

 2 K 1,α (ρHu α ) + (hκ) F ρ,LB +∆α/2 -K 1,α (ρu α ) + u β F ρu β ,LB +∆α/2 -K 1,α (ρu α u β + pδ αβ ) + ρTu α (F α (s) -K 1,α (s)) + ρu α N s ∑ k=2 (g kg 1 ) (G α (Y k ) -K 1,α (Y k )) , (B.14) F ρY k ,FV +∆α/2 = K 2,α (ρY k u α ) +Y k F ρ,LB +∆α/2 -K 2,α (ρu α ) + ρu α (G α (Y k ) -K 2,α (Y k )) . (B.[START_REF] Cai | Adaptive simulations of viscous detonations initiated by a hot jet using a high-order hybrid weno-cd scheme[END_REF] functions.

  ,α c i,βδ αβ c 2 s )u α u β

  ,α c i,βc 2 s δ αβ ) f i -

Table 2 :

 2 Summary of the ability of MHM and Heun (TVD) schemes to simulate flow discontinuities without numerical oscillations.

Table 3 :

 3 Detonation wave velocities D measured with the HLBM simulations with s c = 0.1 for several meshes compared to the Chapman-Jouguet theory D CJ / √ RT 0 ≈ 6.8094746 (Eq. (

	D/ √ RT 0	|D -D CJ |/D CJ
	N 1/2 = 16 6.8094177	8.4 × 10 -6
	N 1/2 = 32 6.8094177	8.4 × 10 -6
	N 1/2 = 64 6.8094552	2.8 × 10 -6

  and previous non-linear analyses indicating bifurcations

		2.0 0.8						0.8			0.8
													0.6
	d(p/p N )/d(t/tc)	3900 3925 3950 3975 4000 4025 4050 4075 4100 t/tc 1.0 1.5 0.4 p/p N 0.6 0.4 0.6 0.0 0.2 0.2 N )/d(t/tc) 0.0 d(p/p -0.2 N 16 -0.2	d(p/p N )/d(t/tc) N 16	0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 p/p N
					N 32				(a) E a = 27RT 0 (N 1/2 = 16)	N 32
		2.0 -0.4	N 64 N 128 N 256					-0.4	N 64 N 128 N 256	0.8
		1.0 1.5 0.6 -0.6 p/p N 0.8	0.8	1.0 (a) E a = 27RT 0 1.2 p/p N	1.4	1.6	1.8	0.6 -0.6 0.8	0.8 N )/d(t/tc) d(p/p	-0.2 0.0 0.2 0.4 0.6 1.0 (b) E a = 27.4RT 0 1.2 p/p N	1.4	1.6	1.8
		0.6	3900 3925 3950 3975 4000 4025 4050 4075 4100 0.6		0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.6 -0.4
		0.4				t/tc			0.4			p/p N
	d(p/p N )/d(t/tc)	(b) E a = 27.4RT 0 (N 1/2 = 32) 3900 3925 3950 3975 4000 4025 4050 4075 4100 1.0 1.5 2.0 p/p N 0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.6 p/p N 0.6 -0.6 -0.2 0.0 0.2 N 32 -0.2 0.2 N )/d(t/tc) 0.0 d(p/p N 64 -0.4 N 128 N 256 -0.4 N 512	0.6 0.8 1.0 1.2 1.4 1.6 1.8 p/p N -0.6 1.0 1.2 1.4 1.6 -0.4 -0.2 0.2 0.4 0.6 0.8 0.0 N 128 0.8 N )/d(t/tc) N 256 d(p/p N 512	1.8
							t/tc						p/p N
						(c) E a = 27.8RT 0							(d) E a = 27.82RT 0
								(c) E a = 27.8RT 0 (N 1/2 = 128)
			2.0									0.8
													0.6
		p/p N	1.0 1.5							d(p/p N )/d(t/tc)	-0.2 0.0 0.2 0.4
												-0.4
				3900 3925 3950 3975 4000 4025 4050 4075 4100		0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.6
							t/tc						p/p N
								(d) E a = 27.82RT 0 (N 1/2 = 512)

Figure 6: Leading shock pressure history (left) and corresponding phase plot (right) of instable one-dimensional detonations for varying activation energies.

Table 4 :

 4 Period) l /t c Asymptotic limits for the time periods of the pressure oscillations measured for each one-dimensional detonation according to the activation energy E a . The first case can be compared with the fifth-order method of Henrick et al.[START_REF] Henrick | Simulations of pulsating one-dimensional detonations with true fifth order accuracy[END_REF].

		Henrick et al. [83]
	27.00 RT 0 11.876650650	11.877201192
	27.40 RT 0 23.780806692	-
	27.80 RT 0 47.678392393	-
	27.82 RT 0 95.399963581	-

Table 6 :

 6 Pre-exponential factor and stability criteria for different reduced activation energies.

		1 )	χ
	20	8.121 ×10 7 2.95
	30	3.356 ×10 8 8.32
	38.23	1.1 ×10 9	14.49
	20	4.605 ×10 9 26.74

  MA β = P LB MA β + P E MA β + P Y MA β , (B.8)

	c α ) .	(B.7)
	By noticing that

Note that this very simple example is unconditionally unstable for u x < 0.

Appendix C. Lattice Boltzmann fluxes

As shown in [START_REF] Wissocq | Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method[END_REF], any LB scheme given by Eq. ( 21) leads to a conservative evolution of mass and momentum that can be written as

In this appendix, the expressions for the mass and momentum fluxes are recalled in the particular case of the D2Q9 lattice, with the convention of Eq. [START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF]. They read:

where x ± = x ± ∆x, y ± = y ± ∆x and where the conventions of Eq. ( 22) have been adopted for the lattice velocities.

Appendix D. Unified density-based model with recursive regularized collision operator

In the present work, the unified density-based model of Farag et al. [START_REF] Farag | A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] is adopted as a LB scheme, together with a recursive regularized collision operator [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation[END_REF]. The unified model is an extension of the so-called pressurebased [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF] and improved density-based [START_REF] Guo | Improved standard thermal lattice Boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows[END_REF] models involving two parameters ζ and κ in the equilibrium distribution