G Wissocq

T Coratger

G Farag

赵崧 S Zhao

P Boivin

P Sagaut

Restoring the conservativity of characteristic-based segregated models: application to the hybrid lattice Boltzmann method

A general methodology is introduced to build conservative numerical models for fluid simulations based on segregated schemes, where mass, momentum and energy equations are solved by different methods. It is here especially designed for developing new numerical discretizations of the total energy equation, adapted to a thermal coupling with the lattice Boltzmann method (LBM). The proposed methodology is based on a linear equivalence with standard discretizations of the entropy equation, which, as a characteristic variable of the Euler system, allows efficiently decoupling the energy equation with the LBM. To this extent, any LBM scheme is equivalently written under a finite-volume formulation involving fluxes, which are further included in the total energy equation as numerical corrections. The viscous heat production is implicitly considered thanks to the knowledge of the LBM momentum flux. Three models are subsequently derived: a first-order upwind, a Lax-Wendroff and a third-order Godunov-type schemes. They are assessed on standard academic test cases: a Couette flow, entropy spot and vortex convections, a Sod shock tube, several two-dimensional Riemann problems and a shock-vortex interaction. Three key features are then exhibited: 1) the models are conservative by construction, recovering correct jump relations across shock waves, 2) the stability and accuracy of entropy modes can be explicitly controlled, 3) the low dissipation of the LBM for isentropic phenomena is preserved.

I. INTRODUCTION

The lattice Boltzmann method (LBM) has emerged as a powerful alternative tool for computational fluid dynamics during the last three decades [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] . Its low dissipation properties [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF] , together with a simple and easily parallelizable algorithm [START_REF] Schornbaum | Massively Parallel Algorithms for the Lattice Boltzmann Method on NonUniform Grids[END_REF] and an ability to handle complex geometries thanks to immersed boundary conditions on an automatically generated Cartesian mesh [START_REF] Touil | Direct and large-eddy simulation of turbulent flows on composite multi-resolution grids by the lattice Boltzmann method[END_REF] have made it competitive for both academic and industrial applications ranging from 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 latticeboltzmann 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] to magnetohydrodynamics [START_REF] Dellar | Moment equations for magnetohydrodynamics[END_REF][START_REF] Dellar | Electromagnetic waves in lattice boltzmann magnetohydrodynamics[END_REF] .

This original numerical method differs from conventional approaches as inherited from the fluid kinetic description at a mesoscopic scale. The Boltzmann equation, balancing inter particles' collisions with their statistical motion, is numerically modelled by a lattice of discrete velocities, further discretized in space and time [START_REF] He | A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit[END_REF][START_REF] Dellar | An interpretation and derivation of the lattice Boltzmann method using Strang splitting[END_REF] . The resulting scheme can be splitted into two simple steps: a local non-linear collision term followed by a node-to-node transport consisting in a mere memory shift [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] . Original lattice Boltzmann (LB) schemes, sometimes referred to as the standard LBM, are based on first-order lattices with low numbers of velocities such as the commonly used D3Q19 lattice [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] . The ensuing method is inherently restricted to isothermal and weakly compressible flows. This is due to two main limitations of these models: 1) the numerical stability of the method is compromised as the Mach number increases, especially in the inviscid limit [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] , 2) the low-order velocity discretization of the lattice induces a Mach-related error in the momentum equation, commonly referred to as Galilean or symmetry-breaking error [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] , and prevents the consideration of temperature variations.

The first limitation has been the topic of many researches in the LBM community these last decades [START_REF] Coreixas | Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations[END_REF] . Their main purpose is to build more stable collision models than the simple but low-dissipative Bhatnagar-Gross-Krook (BGK) approximation [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] . This model indeed suffers from a lack of robustness due to the presence of non-hydrodynamic modes, eventually source of instabilities in the high-Mach and zero-viscosity limits [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF][START_REF] Wissocq | An extended spectral analysis of the lattice Boltzmann method: modal interactions and stability issues[END_REF] . Several families of collision models have then emerged, including multiple-relaxation-time (MRT) [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF][START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF][START_REF] Ginzburg | Two-Relaxation-Time Lattice Boltzmann Scheme: About Parametrization,Velocity, Pressure and Mixed Boundary Conditions[END_REF][START_REF] Geier | The cumulant lattice Boltzmann equation in three dimensions: Theory and validation[END_REF] , regularized [START_REF] Latt | Lattice Boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann methodbased large eddy simulation[END_REF] and entropic models [START_REF] Karlin | Maximum entropy principle for lattice kinetic equations[END_REF][START_REF] Boghosian | Entropic lattice Boltzmann methods[END_REF][START_REF] Ansumali | Minimal entropic kinetic models for hydrodynamics[END_REF][START_REF] Chikatamarla | Entropic lattice Boltzmann models for hydrodynamics in three dimensions[END_REF] . These substantial efforts have successfully led to considerable improvements in the numerical stability of the LBM for high Reynolds and Mach number flows. In this regime, the second limitation, resulting from the velocity space discretization, then becomes the most stringent one.

A natural way to deal with the second issue is to increase the number of lattice velocities, yielding so-called multi-speed approaches [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] . This strategy has recently been employed to successfully build compressible LB schemes [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Frapolli | Entropic lattice Boltzmann model for compressible flows[END_REF][START_REF] Frapolli | Entropic lattice Boltzmann models for thermal and compressible flows[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] . However, it suffers from drawbacks induced by the large number of velocities, especially in three dimensions: it is memory consuming, the implementation of boundary conditions is delicate and the stability is made sensitive by the large number of non-hydrodynamic modes [START_REF] Wissocq | An extended spectral analysis of the lattice Boltzmann method: modal interactions and stability issues[END_REF][START_REF] Siebert | Lattice Boltzmann equation linear stability analysis: Thermal and athermal models[END_REF] . Furthermore, their extension to polyatomic gases is not straightforward, eventually prompting the need for an additional set of distribution functions [START_REF] Frapolli | Entropic lattice Boltzmann models for thermal and compressible flows[END_REF][START_REF] Dellar | Two routes from the Boltzmann equation to compressible flow of polyatomic gases[END_REF] . A common strategy to extend the LBM to thermal and compressible flows while preserving the low stencil of the lattice is to solve the energy equation separately. This is the purpose of so-called segregated methods, which can be grouped into two main classes. The first one, referred to as the double distribution function (DDF) approach, relies on the introduction of a second set of distribution functions, allowing for the recovery of one form of the energy equation [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-distributionfunction 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] . The stability of this method is, however, hardly controllable due to the strong coupling occurring between the two lattices [START_REF] Li | Coupled double-distributionfunction lattice Boltzmann method for the compressible Navier-Stokes equations[END_REF] . The second class is referred to as hybrid methods. In the latter, temperature fluctuations are explicitly considered by coupling the LBM to a finite-difference form of the energy equation, while the Galilean error is handled via appropriate forcing terms [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] . This approach has recently proven its worth for simulations of compressible flows [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][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] 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] .

The present article focuses on hybrid models, regardless the collision kernel adopted for stabilizing the LBM. A central point of this approach involves the coupling between a LB scheme and an explicit discrete energy equation, which can lead to new instabilities [START_REF] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF] . The question of a stable and accuracy-preserving coupling is an active research topic of compressible LBMs [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] . To prevent the coupling issue, one solution exists, based on the characteristic variables of the hyperbolic system formed by the Euler equations. Their evolution obeys a simple advection equation, in one dimension,

∂ L ∂t + Λ ∂ L ∂ x = 0, ( 1
)
where L is a characteristic variable, also referred to as Riemann invariant [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] , and Λ is an eigenvalue of the hyperbolic system. In the absence of discontinuities, the system formed by the characteristic equations is equivalent to the Euler equations. Note that under a linear assumption, Eq. ( 1) comes down to a system of fully decoupled advection equations at constant velocities. The interest of dealing with the characteristic variables is then clear: it allows for an explicit control of one single characteristic by using a well-known discrete scheme for its advection equation. In a sense, this is exactly the purpose of Riemann solvers in numerical fluid dynamics, which aim at building numerical methods for nonlinear systems allowing a correct resolution of the characteristic waves [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . In the context of hybrid LBM, this idea is reflected in models based on the entropy variable [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] , which is precisely one characteristic L of the Euler equations. The entropy advection equation is then linearly decoupled from the rest of the system, modelled by the LBM. Its smart discretization allows an independent control of entropic phenomena without degrading the stability and accuracy of the LBM, which entirely deals with isentropic phenomena such as acoustic and vorticity propagation. This technique has recently extended the applicability of hybrid LBM to transonic and supersonic flows [START_REF] Ribeiro | Buffet Simulations with a Lattice-Boltzmann based Transonic Solver[END_REF][START_REF] Fares | Validation of a Lattice-Boltzmann Approach for Transonic and Supersonic Flow Simulations[END_REF][START_REF] Coratger | Large-eddy lattice-boltzmann modeling of transonic flows[END_REF] . However, for these applications, a new concern can be raised: the conservativity of the model. An explicit discretization of a characteristic equation (1) does not lead to a conservative system, for two reasons: 1) the evolution of all theoretically conserved quantities formed by the density, the momentum and the total energy, is not controlled, 2) Eq. ( 1) is not written under conservative form since the temporal evolution is not exactly balanced by the gradient of a flux, as shown by the Λ pre-factor in front of the spatial derivative. As well identified for the Euler equations [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] and observed with previous LB schemes [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF][START_REF] Hou | Why nonconservative schemes converge to wrong solutions: error analysis[END_REF] , a non-conservative formulation of the energy equation leads to incorrect jump conditions across shock waves. Restoring the conservativity of hybrid LBM is therefore a crucial subject for compressible flows.

The natural way to address the conservativity issue is to introduce a coupling of the LBM with a conserved discrete equation for the total energy, which, unlike the entropy, is a quantity for which the flux is conserved across shock waves.

Unfortunately, the use of classical finite-difference schemes for the energy equation leads to strong coupling instabilities, as mentioned above and in recent publications [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 explanation is quite simple: when discretizing the total energy equation, the characteristic waves are implicitly solved in a non-controlled way by advection schemes whose numerical stability is not a priori ensured. To the authors' knowledge, the first promising step towards a conservative hybrid LB model was recently done by Zhao et al. [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] . In this work, conservative transport schemes are built in a general way by using the knowledge of LB distribution functions in the construction of the total energy flux. The purpose of this approach is to restore some consistency between the LBM fluxes and the total energy one. The idea of using gas distribution functions at cell interfaces to evaluate macroscopic fluxes was the key point of the gas kinetic schemes (GKS), where macroscopic Navier-Stokes variables are recovered from the Boltzmann equation through a Chapman-Enskog expansion [START_REF] Xu | Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes[END_REF] . Unfortunately, although promising, the methodology proposed by Zhao et al. [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] does not ensure the numerical stability of independent characteristics, which notably leads to instabilities for long simulations such as vortex convection. Besides, with the use of a completely new total energy scheme, some recent efforts done to stabilize the hybrid LBM on entropy [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 | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF] would have been done in vain.

The objective of the present work is to develop a general methodology for the construction of conservative total energy models by using the knowledge of a discrete scheme on entropy, which is a characteristic of the Euler system. Ensuring the linear equivalence between a conservative scheme and its non-conservative counterpart, the main advantages of the entropy scheme, linearly decoupled with the LBM, can be preserved. In particular, this allows for an explicit control of the entropy modes, which can be discretized by well-understood advection schemes, without compromising the high accuracy of the LBM for isentropic phenomena 2 . To this regard, the methodology proposed in this work can be connected to Riemann solvers, whose purpose is to construct conservative numerical schemes ensuring a correct transport of the characteristic waves. One big advantage of the proposed methodology is its generality: it can be applied to any LB model, and more generally to any segregated model. The resulting schemes being conservative, it will be shown that correct jump relations are recovered across discontinuities, paving the way for new compressible simulations with hybrid LBM.

The article is organized as follows. Section II describes the theoretical construction of a conservative scheme in a general way. To that end, any LB scheme is re-written as a finitevolume scheme involving fluxes. The linear equivalence between the conservative scheme and its non-conservative entropy counterpart, key point of the methodology, is demonstrated. Section III summarizes the discretization schemes adopted in the present work: a first-order upwind scheme, a Lax-Wendroff scheme and a third-order Godunov type method. They are further compared and validated with the entropy schemes on several classical compressible test cases in Section IV, including discontinuities. The LB model adopted in the present work is described in App. A, and some math-ematical development involved in the theoretical work is detailed in App. B and C.

II. CONSERVATIVE HYBRID LBM SCHEME BASED ON TOTAL ENERGY

In this section, the derivation of the proposed total energy scheme is introduced. It is first done considering the Euler equations for inviscid fluids and then extended to viscous flows. The proposed methodology being completely independent of the scheme used for solving the mass and momentum equations, the lattice Boltzmann algorithm used in the present work is detailed in App. A.

Unless otherwise stated, a dummy field Φ evaluated at (x,t), where x is the position and t is the time, is simply denoted as Φ. Also note that all vector quantities are written in bold.

A. Targeted continuous equations

The targeted continuous conservative system can be written in the inviscid case as

∂ t U + ∂ α F U,c α = 0, (2)
where an implicit summation is done over the index α ∈ {x, y, z}, U = [ρ, ρu x , ρu y , ρu z , ρE] T is the vector of conservative variables, ρ is the density, u α is the α-component of the velocity, E is the total energy and F U,c α is the vector of fluxes in the direction α. The latter reads

F U,c α = [ρu α , ρu x u α + pδ αx , ρu y u α + pδ αy , ρu z u α + pδ αz , ρHu α ] T , ( 3
)
where p is the pressure, H is the total enthalpy and δ αβ is the Kronecker symbol. For an ideal gas, p = ρRT where R is the gas constant and T is the temperature. E and H can be related to T as

E = RT γ -1 + u 2 α 2 , H = E + p ρ = γRT γ -1 + u 2 α 2 , ( 4
)
where γ is the adiabatic exponent. This system can be equivalently written under a nonconservative form as

∂ t V + A V α ∂ α F V,c α = 0. (5)
In the hybrid lattice Boltzmann literature, this is sometimes done by replacing the total energy equation by an entropy equation for stability purposes [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Nie | A Lattice-Boltzmann / Finite-Difference Hybrid Simulation of Transonic Flow[END_REF][START_REF] Ribeiro | Buffet Simulations with a Lattice-Boltzmann based Transonic Solver[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] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF][START_REF] Farag | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF][START_REF] Guo | An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice[END_REF] . In this case, one has V = [ρ, ρu x , ρu y , ρu z , s] T , where

s = R γ -1 ln(p/ρ γ ) (6)
is the entropy, and

A V α =      1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 u α      , (7)
F V,c α = [ρu α , ρu x u α + pδ αx , ρu y u α + pδ αy , ρu z u α + pδ αz , s] T . (8)
Note that this choice of non-conservative system is not commonly encountered in the Euler literature, where the primitive variables [ρ, u α , p] T may be preferred [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . The main advantage of the adopted non-conservative formulation is that the entropy is a characteristic of the system [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . This means that, in the inviscid case, it obeys a simple advection equation:

∂ t s + u α ∂ α s = 0. ( 9
)
Assuming linearity, the velocity appearing in the above equation has to be replaced by a constant, time-and spaceaveraged, velocity. As a consequence, this form of the energy equation is, under a linear assumption, completely decoupled from the mass and momentum equations. This is why coupling an isothermal LB scheme with an energy equation based on the entropy does not affect the linear stability of the LB system, as shown in a previous work [START_REF] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF] . Unfortunately, the main drawback of this formulation is that it is not conservative, hence not adapted to handle discontinuities such as shock waves.

B. Discrete hybrid system based on the entropy equation

In what follows, the mass and momentum equations are solved by a lattice Boltzmann scheme. As a reminder, the purpose of the LBM is to describe a fluid flow using a set of distribution functions f i , designed to mimic the kinetic description of the Boltzmann equation. The standard method can be divided into two parts: 1) a collision step, which intends to model inter-particle collisions, 2) a streaming step, where the distribution functions are simply transported on a lattice of Q discrete velocities denoted as c i . The scheme can generally be written as

∀i ∈ 0, Q -1 , f i (x,t + ∆t) = f coll i (x -c i ∆t,t), ( 10
)
where ∆t is the time step and f coll i are so-called post-collision distribution functions. They depend on the collision model under consideration. More details regarding the lattice and the collision model adopted in the present work are provided in App. A.

By definition of mass and momentum in the LB formalism in absence of body-force term affecting the mass and momen-tum equations, one has

ρ(x,t + ∆t) = Q-1 ∑ i=0 f i (x,t + ∆t) = Q-1 ∑ i=0 f coll i (x -c i ∆t,t), (11)
ρu α (x,t + ∆t) = Q-1 ∑ i=0 c i,α f i (x,t + ∆t) = Q-1 ∑ i=0 c i,α f coll i (x -c i ∆t,t). (12)
Furthermore, ρ and ρu α are collision invariant so that

ρ(x,t) = Q-1 ∑ i=0 f coll i (x,t), ρu α (x,t) = Q-1 ∑ i=0 c i,α f coll i (x,t). (13)
Hence, the discrete time evolution of mass and momentum can be written as

ρ(x,t + ∆t) -ρ(x,t) ∆t = 1 ∆t Q-1 ∑ i=0 f coll i (x -c i ∆t,t) -f coll i (x,t) , (14)
ρu α (x,t + ∆t) -ρu α (x,t) ∆t = 1 ∆t Q-1 ∑ i=0 c i,α f coll i (x -c i ∆t,t) -f coll i (x,t) . ( 15
)
As shown in App. B, this can be systematically re-written under a conservative formulation involving fluxes as

δ t ρ + δ β F ρ +∆β /2 = 0, (16)
δ t (ρu α ) + δ β F ρu α +∆β /2 = 0, (17)
where an implicit summation is done on β ∈ {x, y, z}. The δ t and δ β operators are defined as

δ t Φ = Φ(x,t + ∆t) -Φ(x,t) ∆t , (18)
δ β Φ = Φ(x,t) -Φ(x -e β ∆x,t) ∆x , ( 19
)
e β is the unity vector in the direction β , ∆x is the mesh size and F Φ +∆β /2 (x,t) is the intercell numerical flux of a quantitiy Φ between two cells centered about x and x + e β ∆x. The latter is specified in App. B as function of f coll i for several lattices of interest. Note that, in the present context, mesh size and time step are related to each other through an acoustic scaling [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] , which reads for the D1Q3, D2Q9 and D3Q19 lattices [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] ,

c s ≡ RT re f = 1 √ 3 ∆x ∆t , (20)
where T re f is an arbitrary reference temperature effectively used to tune the CFL number [START_REF] Renard | Improved compressible Hybrid Lattice Boltzmann Method on standard lattice for subsonic and supersonic flows[END_REF][START_REF] Courant | On the Partial Difference Equations of Mathematical Physics[END_REF] , defined as

CFL ≡ ||u|| + γRT ∆t ∆x = ||u|| 3RT re f + γT 3T re f . (21)
As discussed above, the LB scheme can be coupled with an entropy equation to account for the temperature variations. Let us here assume that Eq. ( 9) is discretized by the following explicit scheme:

δ t s + u α (x,t)δ * α s = 0, (22)
where δ * α is a spatial gradient operator associated to the selected numerical scheme for the entropy equation. Without loss of generality, such a spatial discretization can also be rewritten involving δ α of Eq. ( 19) as

δ * α s = δ α F * +∆α/2 (s), (23)
where F * +∆α/2 is a linear function of the scalar field s. Examples of numerical discretizations δ * α and their corresponding numerical fluxes F * +∆α/2 adopted in the present work are specified in Sec. III. It is noteworthy that, as a consequence of the consistency of the scheme,

F * +∆α/2 (s)(x,t) = s(x,t) + O(∆x). (24)
The complete discrete system can then be re-written as

δ t V + A V α δ α F V,d α = 0, (25)
with

F V,d α = [F ρ +∆α/2 , F ρu x +∆α/2 , F ρu y +∆α/2 , F ρu z +∆α/2 , F * +∆α/2 (s)] T .
The system of Eq. ( 25) is not written under conservative form for two reasons : (1) it does not describe the time evolution of conserved quantities U because the entropy s is considered instead of the total energy ρE, (2) a non-identity matrix A V α appears in front of the spatial gradient operator δ α . Furthermore, note that Eq. ( 25) is not closed because the numerical fluxes of the LB solver are likely to involve all the distributions f i of the LB scheme, so that the macroscopic quantities of V may not be sufficient to describe it. Nevertheless, it will be shown that this rewriting is sufficient to derive a conservative hybrid scheme in the next section.

C. Derivation of the conservative hybrid scheme

The purpose of the present section is to derive a conservative total energy scheme based on the numerical scheme adopted for the advection of entropy. The latter being a characteristic of the system, preserving its correct discretization would ensure that the linear stability of the energy equation does not affect its mass and momentum counterparts. To this end, let us multiply Eq. ( 25) by the following Jacobian matrix:

M = ∂ U ∂ V =      1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 h -κ u x u y u z ρT      , ( 26
)
where h = γRT /(γ -1) is the enthalpy and κ = u 2 α /2 is the kinetic energy. This Jacobian matrix allows the passage from non-conserved variables V to conserved ones U, so that one obtains

δ t U + MA V α δ α F V,d α = 0. ( 27
)
Even if this system now involves the conserved variables U, this formulation is still not conservative because of the matrix MA V,d α appearing in front of the space gradient operator. Note that the mass and momentum equations are not affected by this matrix (they still read as in Eqs. ( 16)- (17) which is in conservative form), but the energy equation now reads:

δ t (ρE) + (h -κ)δ α F ρ +∆α/2 + u β δ α F ρu β +∆α/2 + ρTu α δ * α s = 0. ( 28
)
This equation is clearly not conservative. To solve this issue, the continuous equations of Eqs. ( 2)-( 5) can be used. Multiplying Eq. ( 5) by M yields:

∂ t U + MA V α ∂ α F V,c α = 0. ( 29
)
Comparing it with the conservative continuous equation of Eq. ( 2), one can relate the differential form of conserved and non-conserved fluxes:

dF U,c α = MA V α dF V,c α . (30)
The last line of this matrix equality reads

d(ρHu α ) = (h -κ)d(ρu α ) + u β d(ρu α u β + pδ αβ ) + ρTu α ds. (31)
Replacing the differential form with the discrete spatial derivation δ * α and combining it with Eq. ( 28) yields

δ t (ρE) + δ * α (ρHu α ) + (h -κ) δ α F ρ +∆α/2 -δ * α (ρu α ) + u β δ α F ρu β +∆α/2 -δ * α (ρu α u β + pδ αβ ) = 0, (32)
which can equivalently be written as

δ t (ρE) + δ * α (ρHu α ) + (h -κ)δ α F ρ +∆α/2 -F * +∆α/2 (ρu α ) + u β δ α F ρu β +∆α/2 -F * +∆α/2 (ρu α u β + pδ αβ ) = 0. ( 33
) Finally, to obtain a conservative scheme, one can slightly modify the above equation by including the (hκ) and u β pre-factors inside the operator δ α , noticing that it does not affect the consistency and linear properties of the system (as shown in Sec. II E below). It yields

δ t (ρE) + δ α F ρE +∆α/2 = 0, ( 34
)
where

F ρE +∆α/2 = F * +∆α/2 (ρHu α ) (i) + (h -κ) F ρ +∆α/2 -F * +∆α/2 (ρu α ) (ii) + u β F ρu β +∆α/2 -F * +∆α/2 (ρu α u β + pδ αβ ) (iii) . ( 35
)
Thanks to these modifications, Eq. ( 34) is now written under a conservative form. Note that, as shown in App. C, the linear behavior of the scheme is not affected by the location where the pre-factors (hκ) and u β of Eq. ( 35) are evaluated. However, to systematically ensure the symmetry of the scheme, they will be extrapolated at the cell interface by averaging their value between x and (x + e α ∆x).

Let us now look in more details at each term of the resulting flux (35). The first term (i) is nothing but a total energy flux ρHu α discretized by the linear scheme F * +∆α/2 , initially designed to compute the entropy gradients. Regarding the second term (ii), by consistency of the LB scheme with the mass equation [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] and using the fact that

F * +∆α/2 (ρu α ) = ρu α + O(∆x), one has F ρ +∆α/2 -F * +∆α/2 (ρu α ) = O(∆x), (36)
so that this term does not affect the consistency of the total energy scheme: it vanishes as ∆t → 0. As will be shown below, its role is to restore the consistency between the total energy scheme and the mass equation, numerically discretized by another scheme. Finally, to interpret the last term (iii), one may use the fact that, by consistency of the LB scheme [START_REF] Farag | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF] ,

F ρu β +∆α/2 = Π αβ + O(∆x) (37)
= ρu α u β + pδ αβ + Π neq αβ + O(∆x), (38)
where

Π αβ = ∑ i c i,α c i,β f i is the second-order moment of f i , and Π neq αβ = Π αβ -(ρu α u β + pδ αβ )
is its so-called offequilibrium part. Furthermore, by looking at the hydrodynamic limits of the kinetic model using e.g. a Chapman-Enskog expansion [START_REF] Chapman | The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases[END_REF] or a Taylor expansion [START_REF] Wissocq | Hydrodynamic limits and numerical errors of isothermal lattice boltzmann schemes[END_REF] , it can be shown that in the low-Knudsen limit, Π neq αβ ∼ -Π NS αβ , where Π NS αβ is the Navier-Stokes viscous stress tensor defined as

Π NS αβ = µ ∂ α u β + ∂ β u α - 2 D δ αβ ∂ γ u γ + µ b ∂ γ u γ δ αβ , (39)
where µ and µ b are respectively the shear and bulk viscosities, and D is the spatial dimension. As a consequence, the last term (iii) can be re-written as

(iii) ∼ -u β Π NS αβ + O(∆x). ( 40
)
Its role is then multiple. Like the second term (ii), it aims at restoring the consistency between the total energy flux and the discrete momentum equation. It is also remarkable that when considering viscous flows, this term implicitly contributes for the viscous heat expected by the Navier-Stokes equations. Hence, the extension of the derived total energy schemes to viscous flows only requires the explicit computation of the conduction term, which is proposed in the next section.

D. Extension to viscous flows

In this section, the total energy schemes adopted in the present work are extended to the Navier-Stokes case. Let us recall that the continuous total energy equation reads, in presence of viscosity and heat conduction,

∂ t (ρE) + ∂ α (ρHu α ) = ∂ α (u β Π NS αβ ) + ∂ α (λ ∂ α T ), (41)
where λ is the heat conductivity. As mentioned above, the total energy schemes of Eqs. ( 34)-( 35) implicitly include the viscous heat production. Hence, only the last right-hand-side term (conduction term) has to be explicitly discretized. In the present work, this is done by replacing the "Euler" flux

F ρE +∆α/2
appearing in Eq. ( 34) by its Navier-Stokes counterpart F ρE,NS +∆α/2 , defined as

F ρE,NS +∆α/2 = F ρE +∆α/2 -λ δ α T (x + e α ∆x,t). ( 42
)
Note that for a constant conductivity λ , this scheme comes down to a second-order centered scheme for the heat conduction term, since

δ α δ α T (x + e α ∆x,t) = T (x -e α ∆x,t) -2T (x,t) + T (x + e α ∆x,t) ∆x 2 . ( 43
)
However, for a non-constant conductivity, Eq. ( 42) should be preferred since it naturally ensures the conservativity of the equations. The Euler part of the total energy flux F ρE +∆α/2 is provided by Eq. (35), where the linear function F * +∆α/2 depends on the entropy scheme one wants to mimic, and especially its spatial derivative operator δ * α . Their expression is provided in Sec. III for several schemes of interest.

E. Linear equivalence with an entropy-based hybrid scheme

The total energy schemes derived in Sec. II C have been obtained in a systematic way starting from a given discretization of the entropy equation. However, this derivation does not theoretically ensure that the numerical properties of the adopted entropy scheme are conserved. Especially, some assumptions result in neglecting the numerical errors in the construction of the scheme, e.g. replacing continuous derivatives (∂ t and ∂ α ) by their discrete counterparts (δ t and δ α ), as well as the inclusion of (hκ) and u β in spatial derivatives to obtain Eq. (34). In this section, it is shown that, even if the derived total energy schemes are essentially different from the entropy ones, they share equivalent properties in the linear approximation. Hence, provided that the targeted entropy scheme is linearly stable and sufficiently accurate, such properties are recovered in the total energy scheme.

To this end, let us assume that a dummy field Φ can be decomposed between a mean base flow, time-and spaceaveraged, denoted as Φ, and a small perturbation, denoted as Φ . In particular, when applied to the total energy, one assumes that ρE = ρE + (ρE) , where (ρE) ρE. A linearized form of the total energy scheme in the viscid case can then be written as

δ t (ρE) + δ α F ρE,NS +∆α/2 = 0. ( 44
)
Using linear transforms from the conserved perturbations U and their non conserved counterparts V , it is theoretically shown in App. C that, when coupled to the desired mass and momentum schemes, Eq. ( 44) is equivalent to the following linearized scheme on entropy:

δ t s + u α δ * α s = λ ρT δ α T (x + e α ∆x,t). ( 45
)
This equation is nothing but a linearization of the discrete entropy equation (22), supplemented by the heat diffusion. Note that the viscous heating is a purely non linear effect, this is why it does not appear in the above equation. The linear equivalence between the total energy scheme and the entropy one is then demonstrated. As a consequence, the conservative scheme derived in Sec. II C shares all its properties in the linear approximation. In particular:

• the linear stability of the total energy scheme can be directly controled by the choice of advection scheme through

F * +∆α/2 ,
• coupling with a LB scheme does not affect the linear behavior of the latter.

The first point is very interesting for stability purpose. The second point is important for the accuracy of the coupled scheme, since it ensures that the total energy coupling does not degrade the very low dissipation properties of the LBM 2 .

III. TOTAL ENERGY SCHEMES ADOPTED IN THE PRESENT WORK

This section aims at summarizing the schemes adopted in the present work.

Regarding the mass and momentum equations, they are solved by a LB solver based on the unified formulation introduced in a previous work [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] . The latter relies on a regularized collision model [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 methodbased large eddy simulation[END_REF] . The main steps of this algorithm, especially the way f coll i is computed, are recalled in App. A. Regarding the energy equation, it is discretized thanks to the methodology introduced in the previous section. The general form of the scheme is provided in Eq. (34), where the total energy flux is given by Eq. ( 35) in the inviscid case, and by Eq. ( 42) in the viscid one. These fluxes rely on the choice of entropy scheme one wants to mimic, through the appearance of the linear operator F * +∆α/2 , which is related to the numerical gradient δ * α . Three choices are adopted in the present work: (1) a first-order upwind scheme, (2) an isotropic Lax-Wendroff scheme, (3) a monotonic upstream-centered scheme for conservation laws (MUSCL). These standard schemes can be found in several textbooks, e.g. Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . They are recalled below using the non-conventional notations of the present article based on the F * +∆α/2 operator.

A. First-order upwind scheme

In the case of a first-order upwind scheme for an advection equation at constant velocity u, the operator δ * α is defined such that

δ * α Φ = [Φ(x,t) -Φ(x -e α ∆x,t)] /∆x if u α ≥ 0, [Φ(x + e α ∆x,t) -Φ(x,t)] /∆x else. ( 46
)
This can be equivalently re-written as δ α F * +∆α/2 (Φ) with

F * +∆α/2 (Φ) = Φ(x,t) if u α ≥ 0, Φ(x + e α ∆x,t) else. ( 47
)
For a non-constant velocity u, as it is the case in the present non-linear system of equations, a discussion has to be raised regarding the location where the sign of the velocity is evaluated. First note that it has no impact on the linear stability of the proposed algorithm, since a linearization only makes the mean base value u α appear. To systematically ensure the symmetry of the scheme, it is here proposed to evaluate u α appearing in Eq. ( 47) at the cell interface by averaging the values at x and (x + e α ∆x).

B. Lax-Wendroff scheme

In order to increase the order of accuracy of an advection equation discretization, Lax and Wendroff proposed a secondorder scheme preserving a narrow spatial stencil [START_REF] Lax | Systems of conservation laws[END_REF] . As proposed by Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] , its numerical flux can be written as

F * +∆α/2 (Φ) = 1 2 1 -u α ∆t ∆x Φ(x + e α ∆x,t) + 1 2 1 + u α ∆t ∆x Φ(x,t) + (δ αβ -1) u β 2 ∆t ∆x Φ(x + e β ∆x,t) -Φ(x,t) . (48)
One can easily show by Taylor expansions that F * +∆α/2 (Φ) = Φ + O(∆x). However, it can be noted that setting this flux function in Eqs. ( 22)- (23) does not exactly lead to a secondorder accurate scheme in space and time, because the time evolution of u α has not been considered in the construction of the advection scheme. In order to obtain a full second-order accuracy, Eqs. ( 22)-( 23) could be replaced by

δ t s + u α (x,t + ∆t/2)δ α F * +∆α/2 (s) = 0, (49)
which requires an estimation of u α after half a time step. This strategy is not proposed in the present work. Hence, the resulting scheme is expected to be linearly second-order accurate.

C. MUSCL-Hancock scheme

The order of accuracy of the advection equation discretization can be further increased with a so-called MUSCL-Hancock scheme [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF][START_REF] Van Leer | On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe[END_REF] . Such a strategy has been successfully employed to couple an entropy scheme with the LBM [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] Coratger | Large-eddy lattice-boltzmann modeling of transonic flows[END_REF] . In this scheme, the inter-cell flux can be constructed as:

F * +∆α/2 (Φ) = Φ +∆α/2 (x,t) if ũα ≥ 0, Φ -∆α/2 (x + e α ∆x,t) else, (50)
where, to ensure the symmetry of the algorithm, ũα = [u α (x,t) + u α (x + e α ∆x,t)]/2,

Φ +∆α/2 = Φ +∆α/2 + u α 2 ∆t ∆x Φ -∆α/2 -Φ +∆α/2 , (51)
Φ -∆α/2 = Φ -∆α/2 + u α 2 ∆t ∆x Φ -∆α/2 -Φ +∆α/2 , (52)
and

Φ +∆α/2 = Φ + ∆ α 2 , Φ -∆α/2 = Φ - ∆ α 2 . ( 53
)
∆ α 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)] , (54)
where η α = 1 3 2u α ∆t ∆x -sign(u α ) as suggested in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] to obtain a third-order accurate convection scheme in space and time. Yet, note that in the present work, this algorithm is not exactly third-order accurate for two reasons : (1) in the case of a coupled system, the time evolution of u α should be considered to build a high-order scheme, as discussed above with the Lax-Wendroff scheme, (2) in multi-dimensions, a high-order scheme would involve diagonal points, which are not considered here. To recover a high-order scheme even in multidimensions, one could use a directional splitting [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF][START_REF] Strang | On the construction and comparison of difference schemes[END_REF] , which is not adopted in the present work for the sake of simplicity.

D. Summary

The schemes adopted in the present work are summarized in Table I. Note that the order of accuracy is provided in the linearized one-dimensional case. TABLE I: Summary of the adopted discretization schemes.

IV. VALIDATION

In this section, the proposed conservative formulation is assessed on different compressible test cases aiming to exhibit its features. To this end, reference results will be compared with numerical solutions from

• Classical upwind, Lax-Wendroff and MUSCL-Hancock schemes (respectively denoted as s-upwind/LW/MHM) based on the entropy equation ( 9) in non-conservative form. These schemes are recalled in Sec. III.

• Proposed upwind, Lax-Wendroff and MUSCL-Hancock schemes (ρE-upwind/LW/MHM) based on the total energy equation in conservative form, following the proposed method Sec. II.

• The upwind scheme Econs-Up proposed by Zhao et al. [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] .

The specific features we intend to validate are the accuracy, robustness and the linear equivalence between classical sschemes and their ρE-counterparts.

In all the following simulations, the adiabatic exponent is set to γ = 1.4 and the gas constant is set to R = 1, which only affects the dimensionalization of the problems. Given the lattice Boltzmann model described in App. A, the only numerical parameters are the shock sensor coefficient s c and the CFL number, whose values are detailed for each of the forthcoming simulations. Regarding the CFL number, it is, for each case, prescribed to its maximum achievable value. The adopted lattice is the D3Q19 one [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] , and only one cell is considered in the invariant directions of one-and two-dimensional cases.

A. Thermal Couette flow

The first test case we present is the thermal Couette flow. An identical wall temperature T 0 is imposed on two parallel walls. The left wall remains at rest while the right one is moving with a constant velocity u R . When the temperature increase induced by the viscous heat is exactly compensated by heat diffusion, the velocity and temperature profiles reach a steady state [START_REF] Liepmann | Elements of gasdynamics[END_REF] characterized by

T th T 0 = 1 + x H Pr γ -1 2 Ma 2 1 - x H , (55)
u th = x H Ma γRT 0 . (56)
Here, Pr stands for the Prandtl number defined as Pr = µγR/[(γ -1)λ ], Ma = u R / √ γRT 0 is the Mach number of the right wall and H = 1 is the distance between the walls.

The numerical setup is the following: the one-dimensional domain is discretized by 100 nodes. The first and last nodes are located at a distance ∆x/2 from the left and right walls. On the left boundary, a no slip condition is imposed using the classical bounce-back approach [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] on the first node. The right boundary condition also uses a bounce-back scheme and an additional term [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] meant to prescribe a given wall velocity -corresponding here to Ma = 1.3 -and a prescribed wall density enforcing a zero pressure gradient normal to the right wall. The wall temperature T 0 = 1 is set by enforcing the value of the analytical profile of Eq. ( 55) on the first and last points (x = ∆x/2 and x = H -∆x/2). Note that the MHM scheme having a larger stencil, it is replaced by a LW scheme near the boundaries. Other parameters are Pr = 1, µ = 10 -4 , T re f = 10 so that max(CFL) ≈ 0.5. The solution being smooth, the shock sensor is turned off during this test case, s c = 0. Fig. 1 displays the time evolution of the L 2 T error on temperature, defined as

L 2 T = ∑(T -T th ) 2 ∑ T 2 th , (57)
where the sums are performed over the whole simulation domain. As can be seen in Fig. 1 the simulation is started from the expected steady solution. Eventually, due to numerical errors, a short unsteady regime is observed before reaching a plateau. The final steady state L 2 T errors can be found in Table II, they remain of the order of 10 -4 or below. Fig. 2 illustrates the normalized temperature profile using the upwind, LW and MHM schemes. Results are in good agreement with the theoretical solution for all tested schemes, as can be seen from the L 2

T errors previously discussed. The viscous heat term appearing in the total energy equation (41) is inherited from the momentum conservation, see Sec. (II D). The ability to accurately reach the steady state solution of the thermal Couette flow means that the equilibrium between viscous heat and heat diffusion is properly captured by our schemes, thereby validating the implicit computation of the viscous heat by the proposed schemes.

B. Entropy spot convection

The second test case is among the key configurations we need to check. Assuming a constant flow with perturbations of small amplitude, any fluctuation can be decomposed as a sum of entropy, vorticity and pressure fluctuations. These modes are known as the Kovásznay modes [START_REF] Chu | Non-linear interactions in a viscous heatconducting compressible gas[END_REF][START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF][START_REF] George | On chu's disturbance energy[END_REF] . When a LBM solver is coupled to a discretized entropy equation, Kovásznay modes are separated from each other,

• the entropy scheme carries the entropy mode,

• the LBM solver carries vorticity and pressure modes. This means that in the linear limit, entropy fluctuations should behave as if they were uncoupled from the LBM scheme, and the numerical errors of the entropy scheme are likely to impact the entropy mode only. The proposed conservative ρEupwind/LW/MHM schemes are designed to be linearly identical to s-upwind/LW/MHM schemes. This means that for small entropy fluctuations, a given entropy scheme and its proposed linearly equivalent total energy scheme are supposed to give identical results.

To illustrate this particular behavior, a (L × L) fully periodic domain with L = 1 is discretized by a (200 × 200) mesh. The simulation is initialized by density, pressure and velocity components defined as

ρ = ρ 0 1 + e -((x-x c ) 2 +(y-y c ) 2 )/R 2 c , (58)
p = 1, u x = Ma γRT 0 , u y = u z = 0, (59)
where a cold spot of amplitude = 0.01 and radius R c = 0.1 is centered in the middle of the domain (x c = y c = 0.5), superimposed on a constant flow characterized by its Mach number Ma = 2, its density ρ 0 = 1 and its temperature T 0 = 1. The Euler limit is modeled by µ = 0, the shock sensor is turned off s c = 0 and the reference temperature is set to T re f ≈ 104, so that the mean CFL number is 0.2. The Kovásznay entropy mode [START_REF] Chu | Non-linear interactions in a viscous heatconducting compressible gas[END_REF] obeying a passive scalar advection, the reference solution is simply advected horizontally as a frozen pattern by the mean flow. After 20t c corresponding to 20 flow through time periods t c = 1/(Ma √ γ), numerical solutions are compared. Note that an equivalent test case is also performed on the diagonal direction using a slightly modified velocity field, (u

x , u y ) = ( √ 2, √ 2)Ma √ γ/2,
t c = √ 2/(Ma √ γ).
For the horizontal case, entropy schemes s-upwind, s-LW and s-MHM are compared to their linearly equivalent total energy ρE schemes in Fig. 3. For the diagonal case, MHM schemes led to unstable solutions and the upwind ones led to extremely damped amplitudes. Hence, only the results from the LW schemes are reported on Fig. 4. Density amplitude ratios, defined as (max(ρ) -1)/ , are summarized in Table III. We can see that:

• Upwind schemes lead to robust but extremely damped results. This is a consequence of the low order (firstorder accuracy) of this scheme, inducing a large numerical viscosity with the adopted mesh.

• LW schemes lead to slightly damped amplitudes, with about 10% loss in both cases but also exhibit a strong dispersion error.

• MHM schemes lead to accurate results for horizontal but unstable ones for diagonal advection.

Importantly, note the almost point-by-point agreement observed in Figs. 34between a given entropy scheme and its total energy counterpart, as theoretically expected.

Order of accuracy

In order to explain these results, a convergence study is performed using a simpler configuration. A (1t c )-long simulation with Ma = 1 and CFL = 0.2 is performed for different grid resolutions from (100 × 100) to (1600 × 1600) points. Again, both the horizontal and diagonal cases can be found in Fig. 5. For the horizontal case, the expected accuracy orders are observed,

• upwind schemes and Econs-Up are 1 st -order accurate,

• LW schemes are 2 nd -order accurate,

• MHM schemes are 3 rd -order accurate, except for very fine grids where the 2 nd -order accuracy of the LBM dominates the error. For the diagonal case, a slightly different result is observed: the accuracy of the MHM schemes is reduced to the 1 st order. On the contrary, LW schemes preserve a 2 nd order of accuracy whatever the direction of the advection flow.

Stability and precision

In order to summarize the properties of the different schemes, we plot the function (max(ρ)ρ 0 )/ on Fig. 6 for both horizontal and diagonal advections in the aforementioned Ma = 2 case. Note that the analytical solution of the entropy spot is a frozen pattern passively advected over time, such that its amplitude should remains constant. The main observations for the different schemes are,

• Econs -Up 52 solutions crash before 4t c for both the horizontal and diagonal cases.

• Upwind solutions are stable but they also lead to very high damping of the amplitude. For both horizontal and diagonal transport, 50% of initial amplitudes are lost before reaching 5t c . This is a consequence of the numerical viscosity induced by the first-order accuracy of the scheme.

• LW solutions are stable and lead to similar results whatever the direction, where roughly 90% of the initial amplitude remains at the end of the simulations.

• MHM solutions outperform the other ones for the horizontal case with around 3% loss on the maximum amplitude. However, for a diagonal advection, the maximum amplitude grows over time such that the solution rapidly diverges. Again, note the point-by-point agreement for this linear case between entropy schemes and their total energy counterparts proposed in the present work. In summary, for streamlines aligned with an axis of the Cartesian mesh, the MHM scheme highly outperforms the other schemes. However, for more general flows where streamlines are not aligned with the grid, the LW scheme leads to more robust results with in return a small loss of accuracy.

10 -3 log(L 2 T ) -1 slope -2 slope -3 slope Econs -Up s -Up s -LW s -MHM ρE -Up ρE -LW ρE -MHM (a)
For the diagonal advection, the accuracy decrease of the MHM schemes from 3 rd -to 1 st -order and the corresponding unstable results are the consequence of the cross-shaped stencil used in [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] in which the diagonal points are not considered. Note that the high order of accuracy of the MHM could be restored using a directional splitting [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF][START_REF] Strang | On the construction and comparison of difference schemes[END_REF] , but this would lead to a multi-step algorithm that would hardly be implemented in the LB solver used in the present work.

C. Isentropic vortex convection

In this section, a classical isentropic vortex advection is performed. The entropy being constant in this case, it involves the entropy scheme only through numerical errors. Therefore, while the entropy spot advection was primarily meant to assess the stability of the finite difference part of the algorithm, the isentropic vortex is primarily meant to validate the robust- ness of the LBM, as well as a correct decoupling with the entropy characteristic. Compared to the entropy spot, only the initial conditions differ as

(max(ρ) -ρ 0 )/ǫ Econs -Up s -Up s -LW s -MHM ρE -Up ρE -LW ρE -MHM (b)
ρ = 1 - (γ -1) 2 M v 2 e -((x-x c ) 2 +(y-y c ) 2 )/R 2 c 1/(γ-1) , p = ρ γ , u x = Ma √ γ -M v √ γ y -y c R c e -((x-x c ) 2 +(y-y c ) 2 )/(2R 2 c ) , (60)
u y = M v √ γ x -x c R c e -((x-x c ) 2 +(y-y c ) 2 )/(2R 2 c ) .
The vortex and mean Mach numbers are M v = 1/(4π √ γ) and Ma = 2, similarly to previous studies [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] . Other parameters are identical to the entropy spot and the total simulation time is 20t c , again corresponding to 20 flow-through-time periods.

Qualitative results for the classical s scheme and the proposed linearly equivalent ρE-upwind, LW and MHM schemes can be found in Fig. 7, where a very close agreement between s-and ρE-schemes is observed. All entropy and total energy variants of the upwind, LW and MHM schemes are able to preserve the circular shape of the vortex. Only the ρE-LW scheme exhibits little discrepancy with its entropy counterpart, which can be attributed to a non-linear dispersion. The LW scheme is indeed known to suffer from numerical dispersion [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] .

Note that the test case is not reproduced in the diagonal direction because the results are, this time, very similar to the horizontal case, expect for the MHM schemes for which the instability of the diagonal entropy mode, exhibited in the previous section, leads to an instability.

More quantitative results can also be found in Fig. 8, where the normalized amplitudes of the density, temperature and pressure are plotted over time for the different investigated schemes. Results highlight that:

• The Econs-Up scheme is almost immediately unstable.

• s-schemes along with the ρE-upwind scheme lead to a very small damping of the maximum amplitude of less than 1%.

• However, the ρE-LW and ρE-MHM schemes exhibit an unexpected -although very small -behavior. This is attributed to weak non-linear effects that seem to create density and temperature fluctuations with constant pressure. Therefore, this is interpreted as a non-linear shearto-entropy production phenomenon similar to what was pointed out by Renard et al. in a linear framework [START_REF] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF] . The generated spurious entropy wave is then slightly dispersed by the LW scheme as shown in Fig. 7.

Except for these negligible effects observed with the total energy LW and MHM schemes, all results are in good agreement between a given entropy scheme and its proposed total energy counterpart. Additionally, the results presented in Figs. 7-8 are also in line with the expected passive advection of the vortex. After 20t c , the amplitude errors remain below 4% for the three entropy schemes and their linearly equivalent total energy schemes, which is of the same order as reported by previous models [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] .

D. Sod shock tube

The Sod shock tube [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] is a classical test case for compressible solvers. A (L = 1)-long domain is discretized by 400 points. Left and right constant states separated by a discontinuity localized at x = 0.5 are initialized as (ρ L , u L , p L ) = (1, 0, 1) and (ρ R , u R , p R ) = (0.125, 0, 0.1). The CFL number is set using T re f = 9T L , corresponding to max(CFL) ≈ 0.45.

As previously shown by Zhao et al. [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] , the use of the entropy equation leads to inaccurate temperature jumps when the shock Mach number is sufficiently large. The natural way to circumvent this problem is to use a conservative form of the system based on a total energy equation. Therefore, while a point-by-point agreement between entropy and corresponding total energy schemes was assessed on the previous test cases, at least in the linear approximation, substantial discrepancies are now expected.

Results are shown in Fig. 9, where a gap is easily noticeable on the temperature profile between entropy and total energy schemes. While the Econs-Up and proposed total energy schemes are able to accurately capture jump conditions, all entropy schemes underestimate the temperature behind the right facing shock. Comparing entropy and total energy schemes, jump contact discontinuities are different from each other but remarkably close Gibbs oscillations are observed between lin- early equivalent formulations of the same numerical scheme. These oscillations are due to our choice to set µ = 0, s c = 0 and to remove the artificial bulk viscosity introduced by Eq. (A15) for this specific test case. By doing so, we want to emphasize the robustness of all the investigated schemes. Of course, Gibbs phenomenon could have been removed by using a shock sensor or a finite viscosity, similarly to other studies [START_REF] Latt | Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria[END_REF][START_REF] Coreixas | Compressible lattice boltzmann methods with adaptive velocity stencils: An interpolation-free formulation[END_REF][START_REF] Saadat | Extended lattice boltzmann model for gas dynamics[END_REF] . Another possible solution to damp the Gibbs oscillations located close to the contact discontinuity, which is related to an entropy wave, is to rely on a total variation diminushing (TVD) formulation of the energy equation [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . The principle of TVD schemes is to locally decrease the order of accuracy of the numerical scheme by redefining the flux as a combination of a high-order and a first-order flux. When applied to the total energy scheme, it reads

F ρE,NS +∆α/2 = φ α F high +∆α/2 + (1 -φ α )F low +∆α/2 , (61)
where F high +∆α/2 is a high-order flux, computed here either with the LW or the MHM scheme, F low +∆α/2 is a low-order flux, com-puted with the upwind scheme and φ α is a slope limiter function in the direction α. For instance, the van Albada slope limiter [START_REF] Van Albada | A comparative study of computational methods in cosmic gas dynamics[END_REF] reads

φ α = r 2 + r r 2 + 1 , r = ρE(x,t) -ρE(x -e α ∆x,t) ρE(x + e α ∆x,t) -ρE(x,t) . (62)
As illustrated in Fig. 10, this strategy efficiently removes the Gibbs oscillations located close to the contact discontinuity.

E. 2D Riemann problems

From the above test cases, LW and MHM seem to outperform the upwind and Econs-Up schemes. The MHM exhibits an advantageous low numerical dissipation whereas the LW scheme is subject to a spurious numerical dispersion. Except the unstable diagonal advection of fluctuations, which would require an advanced multi-dimensional model, all results are in favor of the use of the MHM scheme. Therefore, the last two test cases are presented with the MHM scheme only. The Gibbs oscillations, of low amplitude with the MHM model, can be addressed thanks to the shock sensor Eq. (A4), as shown below.

In order to assess the ability of the proposed method to handle complex 2D patterns of interacting discontinuities and to demonstrate benefits of using a total energy conservative scheme, configurations 3, 4 and 12 of Lax & Liu 75 are reproduced with s-MHM and ρE-MHM schemes. For these configurations, a 2D (L × L) domain with L = 1 is discretized by a (400 × 400) mesh. The time step is set as ∆t/∆x = 0.15, 0.25 and 0.25 respectively. Initial conditions of the four quadrants can be found in Table IV.

Note that we set µ = 0 in order to mimic an Euler solver, leading to Gibbs oscillations that are removed by setting s c = 1.5, 0.5 and 0.5 for configurations 3, 4 and 12.

Results are displayed in Fig. 11, where the middle and right columns correspond to the s-MHM and ρE-MHM schemes, respectively. Configurations 4 and 12 involving relatively weak discontinuities, both entropy and total energy schemes are able to provide satisfactory results, although a slightly better agreement is observed on ρE-MHM than s-MHM isocontours.

As pointed out by Zhao et al. [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] , stronger discontinuities are involved in configuration 3, where the entropy equation is unable to capture the interaction of the four strong shock waves, as can be seen on Fig. 11. However, the proposed ρE-MHM scheme is perfectly able to capture this configuration, and details are in remarkable agreement with the reference solution. Also note the appearance of intermediate plateaux between the zones at constant density with the entropy scheme, which are due to its inability to recover the Rankine-Hugoniot relations imposed at the initial state. This phenomenon is greatly attenuated with the total energy scheme.

F. Shock-vortex interaction

Inoue & Hattori 76 investigated the shock-vortex interaction for seven different sets of parameters listed from A to G. Previous works [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] Guo | An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice[END_REF] successfully managed to reproduce case C using an entropy equation. Here, our purpose is to show that case G, having stronger mean and vortex Mach numbers Ma = 1.29 and M v = 0.39, can only be accurately reproduced using a total energy scheme. To this end, case G, corresponding to the experimental results by [START_REF] Dosanjh | Interaction of a starting vortex as well as a vortex street with a traveling shock wave[END_REF] , is simulated with both s-MHM and ρE-MHM schemes and then compared to the reference solution [START_REF] Inoue | Sound generation by shock-vortex interactions[END_REF] .

The computational domain is a rectangle of size [0, 28] × [0, 24] discretized by a 1120 × 960 mesh. A Ma = 1.29 planar shock is located at x = 0. Upstream of the shock, a M v = 0.39 vortex is initialized following Eq. ( 60) with (x c , y c ) = (6, 12) and R c = 1, where we implicitly consider an arbitrary system of units where ρ = 1 and RT =1 in the unshocked region. Other parameters such as heat capacity ratio γ = 1.4, Re = ρ √ γR c /µ = 800 and Pr = 0.75 are similar to case C [START_REF] Inoue | Sound generation by shock-vortex interactions[END_REF] . The reference temperature is set as T re f ≈ 8.4 so that Results can be found in Figs. 12-13. Again, a better agreement is found with the total energy ρE-MHM than with entropy s-MHM. While Fig. 13 shows an excellent agreement with reference of the circumferential distributions of sound pressure around the vortex, small defects between entropy and total energy can be observed on the density isocontours Fig. 12 : • The vortex wake is accurately reproduced by ρE-MHM only.

• Rankine-Hugoniot jump relations across the shock are not a solution of the entropy equation. This creates the spurious vertical line located at x ≈ -9.

Again, this test case clearly favors the use of the ρE-MHM as it leads to more accurate results for large values of the Mach number than s-MHM.

V. CONCLUSION

Based on the fact that the entropy equation offers several advantages for the stability and accuracy of segregated models, an original methodology has been introduced to build conservative schemes. It is especially designed for the hybrid LB method, so as to extend its applicability to transonic and supersonic flows.

The construction of the new hybrid schemes relies on two steps: 1) the LB scheme is equivalently re-written as a (conservative) finite-volume scheme, inspired by a previous work [START_REF] Zhao | Toward fully conservative hybrid lattice boltzmann methods for compressible flows[END_REF] , 2) a total energy scheme involving the LB fluxes is systematically derived. Subsequent schemes can then be in- terpreted as classical discretizations of the total energy equation, to which some corrections are added in order to address the coupling instabilities. The key point of this strategy is to end up with a conservative scheme that is linearly equivalent to its non-conservative entropy counterpart, thus preserving most of its advantages. The latter are:

• Entropic phenomena can explicitly be controlled by the choice of an appropriate scheme, notably ensuring their linear stability.

• The low dissipation of the LBM for isentropic phenomena, such as acoustic propagation and vorticity, is preserved.

• The implementation of the viscous heating is made implicit and consistent with the viscous stress tensor mod-elled by the LBM.

The methodology has been used to develop three different schemes, namely a first-order upwind, a Lax-Wendroff and a MUSCL-Hancock schemes. They have been further assessed on classical test cases for compressible flows: a Couette flow, entropy spot and vortex convections, a Sod shock tube, several two-dimensional Riemann problems and a shock-vortex interaction. The simulations allowed validating two important properties of the models: they are conservative and linearly equivalent to their entropy counterpart.

The algorithms proposed in the present work are very similar to previous hybrid LB models coupled with an entropy equation, especially those introduced by Farag et al. [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] . To this extent, they share the ability of the hybrid LBM to deal with complex geometries thanks to the use of immersed-like boundary conditions [START_REF] Coratger | Large-eddy lattice-boltzmann modeling of transonic flows[END_REF] , and the easy parallelization is preserved since the spatial stencil is not modified. As a matter of fact, the only drawback of this new conservative approach, compared to former non-conservative hybrid methods, relies in the additional computational cost induced by the correction terms appearing in the total energy scheme. In this regard, a further study and optimisation of the computational cost of the present method on increasing complexity simulations is planned for a future work.

The present article can be viewed as a first step towards the adaptation of Riemann solvers to the LB formalism. The proposed conservative schemes have indeed been built by ensuring the correct behaviour of a characteristic variable, namely the entropy, which is precisely the purpose of Riemann solvers [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] . Following this adaption would require modifications of the LB scheme in order to properly address acoustic and vorticity waves. Such a task may be the purpose of future work, together with the development of TVD schemes for a better consideration of discontinuities. Notably, this topic has been partially addressed in Sec. IV D, where the effect of a slope limiter on the energy equation has been demonstrated on the Sod test case. However, a complete TVD scheme remains necessary to reduce the Gibbs oscillations induced by shock waves. Furthermore, improvements of the proposed MHM scheme for non-aligned flows will have to be considered. Finally, one may take the advantage of the original finite-volume interpretation of the LB scheme to develop new models such as conservative boundary conditions. and H i are discrete Hermite polynomials in the velocity space, defined as:

H (0) i = 1, (A7) H (1)
i,α = c iα , (A8) H

i,αβ = c iα c iβc 2 s δ αβ , (A9) H

i,αβ γ = c iα c iβ c iγc 2 s c iα δ β γ + c iβ δ γα + c iγ δ αβ . (A10) The forcing term F E i is introduced in the collision step Eq. (A2) to correct the isotropy defect of the lattice [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] . For the equilibrium distribution function of Eq. (A5), it reads

F E i = ω i 2c 4 s H (2)
i,αβ a F,(2)

αβ , (A11)
where, as proposed in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF]

αβ γ , (A13)
where

∂ γ D f eq , (3) αβ γ
= δ αβ ∂ α (ρu 3 α ) + (1δ αβ )∂ ψ (ρu x u y u z ). (A14) In the above equation, no summation is done over the index α and ψ is chosen as ψ = α and ψ = β . In Eq. (A12), the spatial derivative operator is discretized using a second-order centered scheme and the time derivative using a temporal upwind scheme. In Eqs. (A13)-(A14), spatial gradient operators are discretized using a first-order upwind scheme.

The body-force term F

µ b i reads F µ b i = - ω i 2c 4 s H (2)
i,αβ δ αβ ρc 2 s ∂ γ u γ 0.07Ma 2 , (A15)

where Ma = ||u||/(c s γθ ). As proposed in [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] , it is designed to increase the numerical stability of the LB scheme as the Mach number increases. Finally, the off-equilibrium distribution functions f neq i are computed thanks to a recursive regularization [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann methodbased large eddy simulation[END_REF] . It reads Note that the second-order coefficient is here replaced by its traceless counterpart as suggested in [START_REF] Farag | A pressure-based regularized lattice-Boltzmann method for the simulation of compressible flows[END_REF] to increase the robustness of the LB scheme. This filtering process can be interpreted [START_REF] Wissocq | Hydrodynamic limits and numerical errors of isothermal lattice boltzmann schemes[END_REF] as the regularization of the non-equilibrium moment a neq,(2) αα .

f neq i = ω i H (

Appendix B: Lattice Boltzmann fluxes

In this appendix, the following conventions are adopted for the numbering of the lattice velocities: (1) the index i = 0 is attributed to the static velocity, i.e. c 0 = 0, (2) velocities are numbered so that

∀i ∈ 1, Q/2 , c i+ Q/2 = -c i , ( B1
)
where • is the floor function of a real number. In what follows, the notation ī = i + Q/2 is adopted for referring to the opposite index of i. With these conventions, Eq. ( 14) can be re-written as

δ t ρ = 1 ∆t Q/2 ∑ i=1 f coll i (x -c i ∆t,t) -f coll i (x,t) + f coll ī (x + c i ∆t,t) -f coll ī (x,t) , (B2)
or equivalently 16)- (17). This is done below for several lattices of interest. Note that the proposed relationships are, in general, not unique, but the choices adopted below are made to enforce the symmetry of the scheme.

δ t ρ + 1 ∆x Q/2 ∑ i=1 F ρ i (x,t) -F ρ i (x -c i ∆t,t) = 0, ( B3
One can easily check that the mean values of these mass fluxes are

FIG. 1 :

 1 FIG. 1: Convergence study of the L 2 T error of the thermal Couette flow.

FIG. 3 :FIG. 4 :

 34 FIG. 3: Temperature iso-contours of the horizontal entropy spot convection at t = 20t c with t c = 1/(Ma √ γ), Ma = 2 and CFL = 0.2. Top, from left to right: s-upwind, s-LW, s-MHM. Bottom, from left to right: ρE-upwind, ρE-LW, ρE-MHM.

FIG. 5 :

 5 FIG. 5: Convergence study of the energy schemes on the entropy spot convection at Ma = 1. (a) Horizontal convection. (b) Diagonal convection.

FIG. 6 :

 6 FIG. 6: Density amplitude decay of the entropy spot convection. (a) Horizontal convection, with t c = L/U. (b) Diagonal convection, with t c = √ 2L/U.

FIG. 7 :

 7 FIG. 7: Temperature iso-contours of the horizontal vortex convection at Ma = 2, for t = 20t c with t c = 1/(Ma √ γ). Top, from left to right: s-upwind, s-LW, s-MHM. Bottom, from left to right: ρE-upwind, ρE-LW, ρE-MHM.

-

  FIG.8: Density, temperature and pressure amplitude of the vortex convection.

  FIG. 9: Density, Mach and temperature profiles at time t = 0.2 √ RT L /L in the Sod shock tube test case[START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]

FIG. 10 :

 10 FIG.10: Density, Mach and temperature profiles at time t = 0.2 √ RT L /L in the Sod shock tube test case[START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] , using the van Albada slope limiter for high-order LW and MHM schemes.

FIG. 11 :

 11 FIG. 11: Two-dimensional problems of Lax & Liu 75 . From left to right: reference 75 , s-MHM, ρE-MHM. From top to bottom: density fields for configurations 3, 4 and 12 at time t = 0.3, 0.25 and 0.25 respectively.

FIG. 12 :

 12 FIG. 12: Density isocontours of the shock-vortex interaction, case G of Inoue & Hattori 76 , from Dosanjh & Weeks 77 . Left: Navier-Stokes reference 76 , middle: s-MHM, right: ρE-MHM. Contours levels from 1.26 to 1.60 with an increment 0.0029.

2 D δ αβ ρc 2 s

 22 ∂ γ u γδ αβ c 2 s ∂ t (ρ(1θ )) + a C αβ , (A12)D is the spatial dimension and a C αβ is the equilibriumdependent component of the force term,a C αβ = c 2 s u α ∂ β (ρ(1θ )) + u β ∂ α (ρ(1θ )) -∂ γ D f eq ,(3)

f 1 -f 5 + f 2 -f 6 -f 4 + f 8 =f 3 -f 7 + f 2 -f 6 + f 4 -f 8 =

 152648372648 ρu x , (B17) ρu y . (B18) Finally, regarding the x-and y-momentum fluxes, it is sufficient to replace every f coll i appearing in F ρ +∆β /2 by c i,α f coll i to obtain the expressions of F ρu α +∆β /2 . It yields

TABLE II :

 II L 2 T errors after 4.10 6 iterations of the thermal Couette flow. Regarding the MHM model, an average is done between 3.10 6 and 4.10 6 iterations.

TABLE III :

 III Density amplitude ratios of the convected entropy spot. (a) Horizontal advection at t = 20t c , with t c = 1/(Ma √ γ). (b) Diagonal advection at t = 14t c , with

  (x + c i ∆t,t) . (B6) Eqs. (B3)-(B5) are now written under a conservative form involving fluxes. However, the latter involve all the possible directions of the lattice, including diagonal ones, which is not in agreement with the standard formalism involving Cartesian directions only. It is therefore necessary to find a relationship between the fluxes F

	with
	F i ρu α	(x,t) = c i,α	∆x ∆t	f coll i	(x,t) + f coll ī
	ρ i and F i ρu α +∆β /2 and F terparts F ρ ρu α +∆β /2 , which appear in Eqs. (, and their Cartesian coun-
							)
					with
					F i (x,t) = ρ	∆x ∆t	f coll i	(x,t) -f coll ī	(x + c i ∆t,t) .	(B4)
					Similarly, the momentum equation (15) can be written as
					δ t (ρu α ) +	1 ∆x	Q/2 ∑ i=1	F i ρu α	(x,t) -F i ρu α	(x -c i ∆t,t) = 0, (B5)

ACKNOWLEDGMENTS

Part of this research was supported by ANR, Renault, Airbus and SafranTech by the Industrial Chair Program ALBUMS (ANR-CHIND-18-ALBUMS) and by the French project BALBUZARD. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to their high performance computing resources. We also acknowledge support from Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir".

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Lattice Boltzmann collision model

The collision model adopted in the present work is based on the unified model recently 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] , relying on the hybrid recursive regularization (HRR) of Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann methodbased large eddy simulation[END_REF] . The free parameters of the unified model are set, for all the test cases introduced in Sec. IV, to the following constant values: σ = 1, κ = 0, ζ = 0 (see 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] for the definitions of these parameters). In particular, the fact that σ = 1 reduces the HRR model to the more common recursive regularization (RR) [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF] . Furthermore, the D3Q19 [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] lattice is adopted, whose discrete velocities (see Fig. 14) read c * i,x = {0, 1, 1, 0, -1, 0, 0, 0, 1, -1, -1, -1, 0, 1, 0, 0, 0, -1, 1}, c * i,y = {0, 0, 0, 0, 0, 1, 1, -1, 1, 1, 0, 0, 0, 0, -1, -1, 1, -1, -1}, c * i,z = {0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, -1, -1, -1, 0, -1, -1, 0, 0}, (A1) where c * i,α = c i,α ∆t/∆x. The present appendix aims at detailing the main steps of this LB scheme.

Post-collision distribution functions read

where f eq i is the equilibrium distribution function, f neq i is the off-equilibrium part, F E i is a body-force term designed to restore the Galilean invariance at the Navier-Stokes level, F µ b i aims at introducing a stabilizing bulk viscosity and τ is the relaxation time. The latter is related to the fluid dynamic viscosity µ as

where c s = RT re f = 1 √ 3 ∆x ∆t , T re f is a reference temperature and where ν sc is an artificial kinematic viscosity introduced

where s c is a free parameter whose value is specified in each case of Sec. IV. In multi-dimensions, ν sc is computed as the average value of the above expression in each direction. Additionally, in order to spatially filter the artificial viscosity, its local value is retained as the maximal value of all its direct neighbours.

For the D3Q19 lattice, the equilibrium distribution function is expanded onto the orthogonal basis of Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann methodbased large eddy simulation[END_REF] as

i,xxy + H

i,yzz )(ρu x u x u y + ρu y u z u z )

where θ = T /T re f , ω i are the weights associated to lattice velocities defined as With the D1Q3 lattice, composed of three velocities (Q = 3) in one dimension [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] , let us use the convention

i=1 is simply reduced to one element, which corresponds to the flux along the x-axis. Then, one simply has

D2Q9 lattice

With the D2Q9 lattice, composed of nine velocities (Q = 9) in two dimensions [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] , let us use the following convention for the numbering of velocities:

Then, for a given quantity Φ ∈ {ρ, ρu x , ρu y },

where x ± = x ± ∆x and y ± = y ± ∆x. This quantity can be decomposed as x-derivatives and y-derivatives:

x-derivative

This can finally be written as

with

This equality can be applied to both F ρ +∆β /2 and F ρu α +∆β /2 . Using Eq. (B4), this leads to

It can then be easily checked that F ρu α +∆β /2 = Π αβ .

D3Q19 lattice

Regarding the D3Q19 lattice, the similar strategy as for the D2Q9 lattice is employed to construct the Cartesian flux. With the convention adopted in App. A for the lattice velocities, this leads to the possible expression for the fluxes:

As mentioned with the D2Q9 lattice, the expressions of F In this appendix, any space-and time-averaged value of a dummy quantity Φ is denoted as Φ, and its small perturbations are referred to as Φ .

Let us write the complete discrete conservative system as

Then, using the fact that U = ∂U/∂V V = M V , one has

The last line of this matrix system reads

Let us now focus on the last left-hand-side term of this equation representing the fluctuating flux of total energy. As shown by Eq. ( 35) and Eq. ( 42), it reads:

where

This expression can then be considerably simplified. First, one can notice that, by consistency of the scheme adopted for the entropy equation in Eq. ( 22), one systematically has F * +∆α/2 (Φ) = Φ. This can be easily checked for all the adopted schemes of Sec. III. Moreover, by consistency of the mass and momentum schemes (e.g. solved by a LB approach), one has This can also be easily checked from the LB fluxes obtained in App. B. Thanks to these observations, the second lines of Eqs. (C7)-(C8) directly vanish. Furthermore, by using Eq. ( 31), the first term of the fluctuating flux (i) can be expressed as