G Farag

赵崧 S Zhao

G Chiavassa

Pierre Boivin
email: pierre.boivin@univ-amu.fr

Consistency study of Lattice-Boltzmann schemes macroscopic limit

Owing to the lack of consensus about the way Chapman-Enskog should be performed, a new Taylor-Expansion of Lattice-Boltzmann models is proposed. Contrarily to the Chapman-Enskog expansion, recalled in this manuscript, the method only assumes an su ciently small time step. Based on the Taylor expansion, the collision kernel is reinterpreted as a closure for the stress-tensor equation. Numerical coupling of Lattice-Boltzmann models with other numerical schemes, also encompassed by the method, are shown to create error terms whose scalings are more complex than those obtained via Chapman-Enskog. An athermal model and two compressible models are carefully analyzed through this new scope, casting a new light on each model's consistency with the Navier-Stokes equations.

INTRODUCTION

The Navier-Stokes-Fourier (NSF) system of conservation equations is widely accepted, to study mass, momentum and energy conservation in uid systems. Yet, its derivation from a more general and purely atomistic point of view is one of the challenges of the 6 ℎ Hilbert problem [START_REF] Hilbert | Mathematical problems[END_REF] . Formal solutions of the Boltzmann equation (BE) [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] were obtained following perturbation theory [START_REF] Cercignani | The boltzmann equation[END_REF][START_REF] Saint-Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF][START_REF] Gorban | Hilbert's 6th problem: exact and approximate hydrodynamic manifolds for kinetic equations[END_REF] , but nding full thermo-hydrodynamic solutions of the BE remains an active research topic in mathematics [START_REF] Villani | Limites hydrodynamiques de l'équation de boltzmann[END_REF] .

Nonetheless, this lack of theoretical understanding on the link between the NSF and BE formalisms has not slowed down the rapid development of Lattice-Boltzmann methods (LBM), now an invaluable simulation tool widely used in the engineering and scienti c communities. LBM emerged in the 1980 and consist of a speci c BE discretization. (i) First a discrete set of velocities is used to represent the velocity space, leading to the discrete velocity Boltzmann equation (DVBE). (ii) Second, time and space are discretized, as in most computational uid dynamics (CFD) methods. Albeit initially limited to low-Mach athermal ows, the range of applicability has been steadily growing, to encompass compressible ows [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible ows[END_REF][START_REF] Latt | E cient supersonic ow simulations using lattice boltzmann methods based on numerical equilibria[END_REF][START_REF] Guo | An e cient lattice boltzmann method for compressible aerodynamics on d3q19 lattice[END_REF][START_REF] Saadat | Arbitrary lagrangian-eulerian formulation of lattice boltzmann model for compressible ows on unstructured moving meshes[END_REF] , multiphase ows [START_REF] Zhang | Interface and surface tension in incompressible lattice boltzmann multiphase model[END_REF][12][START_REF] He | A lattice boltzmann scheme for incompressible multiphase ow and its application in simulation of rayleigh-taylor instability[END_REF][START_REF] Zhang | Surface tension e ects on two-dimensional twophase kelvin-helmholtz instabilities[END_REF][START_REF] He | On the three-dimensional rayleigh-taylor instability[END_REF][START_REF] Inamuro | A lattice boltzmann method for incompressible two-phase ows with large density di erences[END_REF] and combustion [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF][START_REF] Hosseini | Weakly compressible lattice boltzmann simulations of reacting ows with detailed thermo-chemical models[END_REF] .

In understanding the link between the equations resolved by LBM and the macroscopic NSF system, the so-called Chapman-Enskog (CE) expansion [START_REF] Chapman | The mathematical theory of non-uniform gases, 3rd edn. cambridge math[END_REF] is the most popular method, provided as an Appendix to most LBM papers. Yet, the CE expansion can have limitations in understanding speci c aspects of modern LB methods. For instance, the aforementioned applications (beyond athermal ows) often correspond to Knudsen numbers too high for the underlying LBM theory to hold, but LBMs reportedly yields reasonable results nonetheless. The impact of the choice of collision kernel, central in the method's robustness [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-grosskrook collision model for lattice boltzmann method-based large eddy simulation[END_REF] , is also hard to study with the CE expansion, often carried out with a simpli ed Bhatnagar-Gross-Krook BGK collision model [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems[END_REF] . Last but not least, the CE expansion can not be easily performed for the wide variety of models in which a LB distribution is resolved coupled to another distribution or scalar (which can represent energy, species, or any transported scalar).

The purpose of the present study is threefold. First, we provide a review of the methods traditionally used to derive the macroscopic equations from a given LBM. Second, the implicit assumptions underlying the CE expansion are discussed. Third, we propose a rigorous and systematic method to analyze LBMs, based on a modi ed equation analysis [START_REF] Warming | The modi ed equation approach to the stability and accuracy analysis of nite-di erence methods[END_REF][24][START_REF] Toro | Riemann solvers and numerical methods for uid dynamics: a practical introduction[END_REF] using a Taylor expansion in time and space. Although use of Taylor expansion to that goal is already reported in the LB literature [START_REF] Dubois | Equivalent partial di erential equations of a lattice boltzmann scheme[END_REF][START_REF] Junk | A nite di erence interpretation of the lattice boltzmann method[END_REF][START_REF] Junk | Rigorous navier-stokes limit of the lattice boltzmann equation[END_REF] for athermal models, the presented method is the rst -to the authors knowledge -to encompass arbitrary LB numerical schemes with multi-physics coupling, arbitrary collision kernel, arbitrary force terms and arbitrary non-dimensional numbers. It will be shown that the method allows to identify error terms beyond the CE expansion, necessary to fully understand recent LB models.

The article is organized following these three goals. After a brief Section I introducing the NSF system of equations along with the necessary notations, Section II focuses on the continuous BE.

We recall two popular methods used to analyze it and derive a NSF system from the BE, namely the CE expansion [START_REF] Chapman | The mathematical theory of non-uniform gases, 3rd edn. cambridge math[END_REF] and the Grad moment system [START_REF] Grad | On the kinetic theory of rare ed gases[END_REF] . Section III discusses the application of the CE formalism to LBMs. In particular, we will point out the lack of consensus found in the literature around the CE expansion. Underlying assumptions and limitations are also discussed.

Section IV contains the principal novelty of the present work, following the arguments presented in Sections II-III, and proposes an alternative to the CE expansion formalism. The stepby-step algorithm to build and understand a LB scheme is thoroughly explained. Resting on a naive Taylor expansion of the numerical scheme this method is seen to be fully deductive and ansatz-free in the sense that its derivation automatically and unequivocally gives the conditions for the scheme to be consistent to an expected set of macroscopic equations in the small time-step limit Δ → 0, while keeping the so-called acoustic scaling coe cient Δ /Δ constant [START_REF] Krüger | The lattice boltzmann method[END_REF] . As a rst textbook example, the classical athermal BGK [START_REF] Krüger | The lattice boltzmann method[END_REF] is analyzed through the scope of the Taylor expansion in Section V. Then, in light of the proposed step-by-step algorithm, Section VI proposes a new interpretation of the LB collision kernel strictly based on macroscopic equations instead of the usual kinetic interpretation. Lastly, the scope of our new theoretical framework is illustrated for two advanced LB thermal models recently published by our group, namely the RR- [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible ows[END_REF][START_REF] Renard | Improved compressible hybrid lattice boltzmann method on standard lattice for subsonic and supersonic ows[END_REF] and RR-7 models, respectfully in Sections VII and VIII.

I. THE NAVIER-STOKES-FOURIER SYSTEM

Before carrying out any comparison between the BE and NSF systems, it is useful to introduce the NSF governing equations, along with appropriate de nitions.

A. Navier-Stokes de nitions

Mass and momentum conservation read

+ = , (1)
+ + -T = F , (2)
where is the volume mass, is the local velocity vector and is the pressure. In addition, and F are respectively any forcing term in the mass an momentum equations. These forces can model physical phenomena e.g. gravity or mass source, but they can also correspond to numerical terms such as sponge-zones [START_REF] Feng | Solid wall and open boundary conditions in hybrid recursive regularized lattice boltzmann method for compressible ows[END_REF] . Lastly, T is the stress tensor,

T = + - 2 3 , (3)
with the shear viscosity. The bulk viscosity is neglected in the framework of this paper, but can readily be included in the analysis.

Recombining Eqs. (1, 2) we obtain the kinetic tensor equation

+ + -T + -T = F + F - , (4)
not to be confused with the kinetic energy evolution equation, corresponding to half the trace of the tensor evolution Eq. ( 4). When the ow is assumed to be athermal, the system is fully closed by assuming, e.g.

= 2 , (5)
where is the constant sound speed.

B. Fourier system de nitions

When thermal e ects cannot be neglected, one needs to consider additionally the total energy density equation

+ ( + ) + -T = F + , (6)
with the total energy de ned as the sum of internal and kinetic energies, = + /2. In Eq. ( 6), corresponds to the heat ux, and is an energy source.

To close the NSF system, a thermodynamic closure is required, e.g.

= , = , (7)
with = / , the perfect gas constant, the molecular weight, the temperature and the mass heat capacity at constant volume. Laws for the heat transfer are also required, e.g.

= - / , (8)
where is the heat conductivity [START_REF] Liepmann | Elements of gasdynamics[END_REF] . Having introduced the NSF system, let us now move forward with the BE system.

II. CONTINUOUS BOLTZMANN EQUATION ANALYSIS

The Boltzmann equation (BE) corresponds to the kinetic description of a gas out of equilibrium. In the absence of external forces, and assuming the gas to consist of small hard spherical particles bouncing elastically against each other, the BE reads

+ = ( ) . (9)
where (ξ, x, ) is the density probability function of nding particles with velocity vector ξ = , at location x = and time . The left hand-side of (9) corresponds to the free streaming of the particles, and the right hand-side takes into account their collisions. Gas kinetic theory dictate that will relax through collisions towards a Maxwellian equilibrium distribution function [START_REF] Maxwell | V. illustrations of the dynamical theory of gases.-part i. on the motions and collisions of perfectly elastic spheres[END_REF] , e.g.

= (2 ) 3/2 -(||ξ-u|| 2 )/2 , (10)
for monoatomic ideal gases.

For mass, momentum and energy to be conserved through collision , they must relate to the density probability functions through 36

= ∫ ξ = ∫ ξ , (11)
= ∫ ξ = ∫ ξ , (12)
= 1 2 ∫ ξ = 1 2 ∫ ξ , (13)
where , and are the macroscopic quantities, respectively the density, velocity vector and total energy.

For the sake of simplicity, let us consider for now the most simple collision operator, namely the Bhatnagar-Gross-Krook (BGK) collision kernel [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems[END_REF] ,

( ) = - 1 ( - ) , (14)
where is a characteristic relaxation time related to the viscosity of the uid. This approximation holds reasonably when the Knudsen number Kn, ratio of the mean free path of the particles to a characteristic ow length, is small. Note that Kn can also be expressed in terms of Reynolds Re and Mach Ma numbers through the Von Kármán relation Kn ∝ Ma/Re.

The formal derivation of the continuum NSF equations from the BE is not straightforward and usually needs a perturbation analysis which is by construction only valid for some speci c asymptotic cases [START_REF] Saint-Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF][START_REF] Gorban | Hilbert's 6th problem: exact and approximate hydrodynamic manifolds for kinetic equations[END_REF][START_REF] Villani | Limites hydrodynamiques de l'équation de boltzmann[END_REF] . For the purpose of this article, we shall recall hereafter two historical methods allowing to link the kinetic BE to the continuum NSF system.

A. Hilbert and Chapman-Enskog expansions

The question of a systematic derivation of the NSF system from kinetic theory is owed to Hilbert [START_REF] Hilbert | Begründung der kinetischen gastheorie[END_REF][START_REF]incompressible two-phase ows at high density ratio[END_REF] . Its rst attempt is based on the assumption that Kn 1 and that the collision characteristic time = ˜ , where is a small parameter usually identi ed with Kn. The rescaled BE then reads

+ = - 1 ˜ ( - ) . (15)
A singular perturbation procedure is then performed by taking the limiting case → 0. Then, we can search for the solution as an in nite expansion,

= ∞ =0 ( ) = (0) + (1) + 2 (2) + ... , (16)
where all ( ) are O (1). Making the ansatz that Eq. ( 16) is convergent we insert it inside Eq. (15).

Then by assuming a scale separation between orders in and collecting terms by orders we end up with an in nite hierarchy of equations

(0) - = 0 , (17)
- 1 ˜ ( ) = ( -1) + ( -1) , (18)
with > 0. This hierarchy of equations, called Hilbert expansion, exhibits some remarkable properties. The rst equation con rms that the zeroth order distribution (0) should match the equilibrium distribution

. The second equation shows that the ℎ equation depends on the ( -1) ℎ distribution only in a sequential manner. Then, by truncating at any order in the in nite expansion it is possible to get an approximate solution of the BE.

The problem of this solution is that it fails to capture very steep gradients where Kn = O (1)

such as in boundary layers or shocks [START_REF] Cercignani | The boltzmann equation[END_REF][START_REF] Chapman | The mathematical theory of non-uniform gases, 3rd edn. cambridge math[END_REF]39 . In those regions the distribution rapidly changes on a time scale of order Kn meaning that / scales as ∼ /Kn. The Hilbert expansion is then ill-equipped to deal with such applications. A way to circumvent this problem is to use the popular Chapman-Enskog expansion [START_REF] Cercignani | The boltzmann equation[END_REF][START_REF] Chapman | The mathematical theory of non-uniform gases, 3rd edn. cambridge math[END_REF] instead of the Hilbert expansion.

The only di erence with the Hilbert expansion lies in the fact that the time derivative is now also expanded,

= ∞ =0 = 0 + 1 1 + 2 2 + ... , (19)
where / denotes the contribution from the ℎ order to the physical time derivative / .

Plugging Eqs. (16, 19) inside Eq. ( 15) leads to

∞ =0 + ∞ =0 ( ) = - 1 ˜ ∞ =0 ( ) - . (20)
Assuming a scale separation between orders in and collecting terms by orders we end up with a new in nite hierarchy of equations. The leading orders in now read

(0) - = 0 , (21)
- 1 ˜ (1) = (0) 0 + (0) , (22)
- 1 ˜ (2) = (1) 0
+ (0) 1 + (1)
.

Again, the zeroth order distribution satis es (0) = , but a slight di erence appears in higher orders, the ℎ equation now depending not only on the ( -1) ℎ but also on any ( -) ℎ order with < . This recursive behavior means that CE expansion only addresses low-Knudsen solutions with depending only implicitly on time via the macroscopic variables appearing inside the Maxwellian Eq. (10). In other words, the CE expansion only describes solutions ( ) with autonomous time dependencies ( ( ), ( ), ( ), ∇ ( ), ∇ ( ), ∇ ( )) with ∇ the ℎ order rank space derivative. More general solutions are simply out of the scope of the CE expansion 40 .

The next step is to take successive moments of this in nite hierarchy. Integrating the rst three order moments of Eq. ( 22) lead to the Euler equations

0 + = - 1 ˜ ∫ (1) ξ , (24) 0
+ + = - 1 ˜ ∫ (1) ξ , (25)
0

+ ( + ) = - 1 ˜ ∫ (1) ξ . (26)
Similarly, the rst three order moments of Eq. ( 23) lead to

1 + ∫ 0 + (1) ξ = - 1 ˜ ∫ (2) ξ , (27) 1
+ ∫ 0 + (1) ξ = - 1 ˜ ∫ (2) ξ , (28) 1 +
∫ 0 + (1) ξ = - 1 ˜ ∫ (2) ξ , (29)
which can be interpreted as a correction to the Euler equations Eqs. (24-26), leading to the NSF system of equations. Note that this system is not closed yet because we do not know how to evaluate ∫ ( ) ξ with > 0 and = [1, , ]. The CE expansion being a formal search of BE solution any constraint can be used to close the system. Thus, solvability conditions are applied,

∫ ( ) ξ = 0 , > 0 . (30)
This is an essential step of the CE expansion because it allows to close the system and prevents the in nite hierarchy to impact the (low) orders of interest. Using the solvability conditions, neglecting higher orders and collecting Eqs. (24-29) the NSF system reads

+ = O ( 2 ) , (31)
+ + + ∫ (1) ξ = O ( 2 ) , (32)
+ ( + ) + ∫ (1) ξ = O ( 2 ) , (33)
this system being closed by moments of Eq. ( 22) [START_REF] Krüger | The lattice boltzmann method[END_REF] .

To summarize, the main CE expansion assumptions are the convergent nature of the expansion (16), the scale separation between orders in and the solvability conditions (30).

B. Grad moment system

Aside the well known CE expansion, another attempt to link the BE to the continuum mechanics was performed by Grad [START_REF] Grad | On the kinetic theory of rare ed gases[END_REF] . It is of historical importance for LBMs as it introduced the use of Hermite polynomials H ( ) to analyze the BE. The main idea is to project the distribution onto the Hermite polynomials [START_REF] Grad | Note on n-dimensional hermite polynomials[END_REF] basis composed by ℎ rank symmetric tensors H ( ) of ℎ degree polynomials in ξ, leading to

(ξ) = (ξ) ∞ =0 1 ! a ( ) : H ( ) (ξ) , (34)
a ( ) = ∫ H ( ) ξ , (35)
(ξ) = 1 (2 ) 3/2 -ξ•ξ/2 , (36)
H ( ) (ξ) = (-1) (ξ) ∇ (ξ) . ( 37
)
Assuming that is su ciently well approximated by a ℎ order truncated version of Eq. ( 34),

(ξ) ≈ (ξ) =0 1 ! a ( ) : H ( ) (ξ) , (38)
it is possible to get a closed system of equations by injecting Eq. ( 38) into the BE (9). Taking moments of this system and using the orthogonality properties [START_REF] Krüger | The lattice boltzmann method[END_REF][START_REF] Grad | Note on n-dimensional hermite polynomials[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] of the Hermite polynomials lead to a closed set of macroscopic equations, the most famous one being the Grad-13 system of equations [START_REF] Grad | On the kinetic theory of rare ed gases[END_REF]39 that describes the evolution of 13 di erent moments ( (0) , [START_REF] Hilbert | Mathematical problems[END_REF] , [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] and [START_REF] Cercignani | The boltzmann equation[END_REF] ).

Because of the completeness of the Hermite basis, when → ∞ the Grad moment method explores every possible BE solutions.

C. The in nite hierarchy of equations and the closure problem

By performing a CE or a Grad moment analysis we are trying to reduce the in nite cascade of moments to a nite system. In this process we are purposely losing some information and we expect that this ansatz was relevant for the considered physical phenomena. In the context of CFD this usually means that we are trying to extract low-Knudsen solutions out of the discretized LB model.

Whether we use the CE or the Grad moment expansion, both methods rely on a truncated expansion. To the best of the authors' knowledge the range of validity of such nite expansions is usually not questioned nor recalled in the LB community.

For instance, the CE expansion rst assumes that the time derivative Eq. ( 19) and populations

Eq. ( 16) can be expanded, but none of these assumptions seems to be explicable, particularly when we also remark that this leads to a singular perturbation analysis, meaning that the smallness parameter appears in front of the highest order derivative in the equation. Such singular perturbation expansion is known to exhibit unphysical behaviors when truncated [START_REF] Hunter | Asymptotic analysis and singular perturbation theory[END_REF] , i.e. the CE expansion is asymptotic rather than convergent [START_REF] Krüger | The lattice boltzmann method[END_REF]39,40,[START_REF] Grad | Asymptotic theory of the boltzmann equation[END_REF][START_REF] Öttinger | Beyond equilibrium thermodynamics[END_REF] . Therefore, for a given physical application there is an optimal number of terms to keep in the in nite CE expansion because higherorder terms may introduce unphysical behaviors. For example, higher order approximations than NSF, namely the Burnett and Super-Burnett equations, can be derived from a truncated CE expansion of respectively the Grad-13 and Grad-26 systems [START_REF] Grad | Asymptotic theory of the boltzmann equation[END_REF][START_REF] Struchtrup | Macroscopic transport equations for rare ed gas ows[END_REF] but also directly from the BE itself.

However, negative viscosity for high gradients and short wave instabilities of the Burnett and Super-Burnett equations [START_REF] Gorban | Hilbert's 6th problem: exact and approximate hydrodynamic manifolds for kinetic equations[END_REF] are reported in the literature. This phenomenon was rst observed by Bobylev 47,[START_REF] Bobylev | Instabilities in the chapman-enskog expansion and hyperbolic burnett equations[END_REF] . This means that higher order approximations in the CE expansion may lead to less stable and less physical results [START_REF] Gorban | Hilbert's 6th problem: exact and approximate hydrodynamic manifolds for kinetic equations[END_REF][START_REF] Struchtrup | Macroscopic transport equations for rare ed gas ows[END_REF] , endorsing a non convergence of clipped in nite expansions such as the CE expansion.

However, when trying to link LBMs with some macroscopic equations we usually resort to a truncated CE expansion. This link is obtained with strong assumptions whose practical use can be reasonably questioned by the existence of Bobylev instabilities for Burnett and Super-Burnett models.

It can also be highlighted that the CE expansion is often performed with a very simple BGK collision operator, which is now hardly used in practical applications for its behavior in shear ows [START_REF] Coreixas | Recursive regularization step for high-order lattice boltzmann methods[END_REF] and porous ows [START_REF] Ginzburg | Consistent lattice boltzmann schemes for the brinkman model of porous ow and in nite chapman-enskog expansion[END_REF] . Although the formal CE expansion is possible with the quadratic Boltzmann collision operator [START_REF] Cercignani | The boltzmann equation[END_REF][START_REF] Chapman | The mathematical theory of non-uniform gases, 3rd edn. cambridge math[END_REF] , extension to complex collision kernels such as Regularized [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible ows[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible ows[END_REF][START_REF] Renard | Improved compressible hybrid lattice boltzmann method on standard lattice for subsonic and supersonic ows[END_REF] or Entropic [START_REF] Frapolli | Entropic lattice boltzmann model for compressible ows[END_REF][START_REF] Frapolli | Entropic lattice boltzmann model for gas dynamics: Theory, boundary conditions, and implementation[END_REF][START_REF] Frapolli | Lattice kinetic theory in a comoving galilean reference frame[END_REF] remains unknown, at least to the authors' knowledge. In those kernels at least a part of the non-equilibrium population is systematically ltered out and replaced by a reconstructed population. This is therefore out of the scope of the formal CE expansion because nding an explicit solution ( ) as a function of ( -1) , ( -2) , ... as simply as in Eqs. (21-23) would require to invert a complex collision kernel.

We have now presented the basis of the two most common analysis tools applied to the BE.

Let us now focus exclusively on Lattice-Boltzmann methods, and leave the continuous BE.

III. LATTICE-BOLTZMANN METHODS

A. From Boltzmann to Lattice-Boltzmann

Besides the time and space discretization on cartesian grids, the main di erence between the BE and LBMs is the velocity space discretization.

Assume a dimensional lattice with discrete velocities , referred to as the lattice [START_REF] Krüger | The lattice boltzmann method[END_REF] .

Through discretization, Eqs. (11-13) become

= =1 = =1 , (39)
= =1 = =1 , (40)
= 1 2 =1 = 1 2 =1 , (41)
where and are the discretized probability density functions corresponding to their BE counterpart from Eqs. (9, 10). Accordingly the BE is also discretized such that we should now solve BEs, one for each discrete velocities ,

+ = - 1 ( - ) . (42)
The last piece of the puzzle is to project the distributions onto the Hermite basis using Eq. (38).

Because of the boundedness of the lattice a natural clip appears inside the Hermite basis. Each lattice only being able to represent independent Hermite polynomials, any Hermite polynomial beyond is in the span of the base.

On the bright side, this implies that Eq. ( 38) no longer involves an in nite hierarchy of equations, so it is enough to search for the solution as a nite expansion.

B. The Chapman-Enskog expansion in LBM

In a context of growing interest on LBMs, new models are being published on a monthly basis, usually including a CE expansion appendix. That CE expansion is usually left unquestioned, yet there is little consensus on its application to LBM.

To highlight this lack of consensus, we included in Tab. I a summary of di erent CE expansions found in four reference LBM textbooks [START_REF] Krüger | The lattice boltzmann method[END_REF]39,[START_REF] Mohamad | Lattice Boltzmann Method[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF] . Table I shows that only [START_REF] Mohamad | Lattice Boltzmann Method[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF] agree with Reference textbook f = ... / = ... / = ...

Mohamad 54 ∞ =0 ( ) 1 + 2 2 1
Guo and Shu 55

∞ =0 ( ) 1 + 2 2 1 Krüger et al. 30 ∞ =0 ( ) ∞ =1 1 Succi 39 ∞ =0 ( ) ∞ =1 ∞ =1
each other. They chose to expand the time derivative not as an in nite expansion but as a linear composition of a fast convective time 1 and a slow di usive time 2 , in contrast to [START_REF] Krüger | The lattice boltzmann method[END_REF]39 that considers -as the historical CE -/ as mathematical derivatives. Under the assumption that each physical phenomenon takes place on the same scale, space derivative is almost never expanded, except by 39 .

Despite these di erences, and reassuringly, these expansions all lead to the same lowest orders corresponding to the expected Navier-Stokes equations. In analyzing error terms in the higher orders, however, these di erences will lead to di erent results.

IV. TAYLOR EXPANSION OF A GENERIC LATTICE-BOLTZMANN SCHEME

LBMs being extensively used by engineers and researchers as a thermo-hydrodynamic solver, it can be interesting to analyze it purely as a CFD method for macroscopic equations. In other words, knowing a given LB scheme we could simply expand it as a Taylor series in small parameter Δ → 0. Then, we recast the scheme whose unknowns are into schemes whose unknowns are macroscopic variables such as mass, momentum and stress tensor. Lastly, we can a priori deduce in terms of nondimensional numbers consistency conditions of this scheme. This naive formalism allows to bypass the CE expansion assumptions, now only needing Δ to be small, as is usual for consistency studies. This strategy will also allow to simultaneously tackle possible numerical couplings between LBM and other schemes.

Hermite moments 45

,( ) 1 ...

Raw moments 46

,( ) 1 ...

Lattice isotropy defect 47

A. Lattice-Boltzmann de nitions

To facilitate the reading, let us now introduce LB speci c quantities. For future reference, all quantities introduced in this Section are summarized in Table II.

The rst Hermite polynomials read

H (0) = 1 , H (1) = , H (2) = -2 , (43)
H (3) = -2 [ + + ] , (44)
where the lattice sound speed is = Δ /( √ 3 Δ ) for standard lattices. Higher order polynomials do not generally belong to the Hermite base of standard lattices and are not provided here [START_REF] Krüger | The lattice boltzmann method[END_REF] .

For a given Hermite basis and an arbitrary population we de ne its Hermite moments ,( ) 1 ...

and macroscopic (raw) moments

,( )

1 ... as ,( ) 1 ... = H ( ) 1 ... , (45)
,( )

1 ... = 1 ... . (46)
Since the number of discrete velocities is nite, there always exists an order involving a non-zero Table III: Third-order isotropy defects of standard lattices.

, ,( )

1 ... = ∫ 1 ... c - ,( )
1 ... . (47)
For nearest-neighbors lattices, this isotropy defect appears from the third order ( = 3) [START_REF] Krüger | The lattice boltzmann method[END_REF] . For future reference, the term is provided in Table III for the athermal LB model detailed in the next Section.

Next, it is convenient to de ne several populations besides , appearing at di erent stages of an LBM algorithm :

• : The total population is the most important population because mass, momentum, etc are macroscopic moments of this population.

• : The equilibrium population, here de ned as

≡ H (0) ,(0) + H (1) 2
,(1) + H (2) 2 4 ,(2) + H (3) 6 6 ,(3) + ... , (48)
which is usually a Maxwellian projected onto the Hermite basis and properly truncated [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] .

• : The non-equilibrium population.

≡ - . ( 49
)
To ensure 2 order accuracy [START_REF] Krüger | The lattice boltzmann method[END_REF] , o set distributions ( , ) are introduced

• : The modi ed total population, de ned by

≡ + 1 + Δ 2 - Δ 2 , (50)
≡ + - Δ 2 , (51)
≡ H (0) ,(0) + H [START_REF] Hilbert | Mathematical problems[END_REF] 2

,(1) + H (2)
2 4 ,(2) + H [START_REF] Cercignani | The boltzmann equation[END_REF] 6 6 ,(3) + ... .

Note that the forcing scheme considered here was presented by [START_REF] Guo | Discrete lattice e ects on the forcing term in the lattice boltzmann method[END_REF] .

• : The modi ed non-equilibrium population ensuring a 2nd order accurate BGK scheme,

≡ - + Δ 2 . (53)
• : The population at the end of a collide step,

≡ + 1 - Δ + Δ 2 , (54)
≡ + 1 - Δ 2 + Δ 2 , (55)
where we used the shorthand = + Δ /2. is the population that is streamed during the propagation step. Note that depending on the collision kernel may be replaced by a reconstructed non-equilibrium population ˜ .

B. Structure of a generic lattice-Boltzmann scheme

A LBM model time iteration typically consists of the steps listed in Tab. IV. Let us detail these steps one by one. Although these steps may vary depending on LBM, the general structure remains. At the center of the algorithm lies the collision

( , x) = ( , x) + 1 - Δ ( , x) + Δ 2 ( , x) , (56)
and streaming steps

( + Δ , x) = ( , x -c i Δ ) . (57)
They remain generally identical for all LBM as they are the main ingredients responsible for the method's computing e ciency [START_REF] Löhner | Towards overcoming the les crisis[END_REF] , and low dissipation [START_REF] Marié | Comparison between lattice boltzmann method and navierstokes high order schemes for computational aeroacoustics[END_REF] .

Steps 1 (resp. 7) corresponds to the link between the macroscopic quantities and the population (resp. and back), and vary depending on the model. The force update (step 2) is also model dependent, but is generally computed from the macroscopic variables.

Step 3 includes the de nition of the collision kernel. Computation of will be highly dependent on the collision model, e.g. BGK [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems[END_REF] , multi-relaxation times models [START_REF] Lallemand | Theory of the lattice boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF][START_REF] Di Ilio | Simulation of turbulent ows with the entropic multirelaxation time lattice boltzmann method on body-tted meshes[END_REF] , regularized models [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-grosskrook collision model for lattice boltzmann method-based large eddy simulation[END_REF] , and so on.

Step 6 is optional, and corresponds to the resolution of a set of quantities φ, not encompassed in the main LBM population. φ is a dynamic eld, and can correspond to scalars (as in hybrid methods) or additional probability functions (as in double distribution functions). They can represent, e.g.

-Energy [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible ows[END_REF] -Species, if the ow is multi-constituent [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF] -an equation for turbulence modeling (e.g. for Spalart-Almaras ormodels), -a liquid mass or volume fraction, for multiphase ows.

While usually presented di erently [START_REF] Krüger | The lattice boltzmann method[END_REF]39,[START_REF] Mohamad | Lattice Boltzmann Method[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF] this structure is shared by all LB schemes using Hermite polynomials. Note that the coupling with the FD/FV/DDF scheme, here identi ed as step 6, can be performed at any moment in the algorithm except after Step 7, because the macroscopic quantity update may depend on it.

C. Taylor expansion

Now that the basic numerical LB scheme structure has been recast, we can compute its Taylor expansion. We introduce the distribution's Taylor Expansion in space as

( , x -c Δ ) = ( , x) + ∞ =1 (-Δ ) ! ( , x) , (58)
where is a dummy index such that / = / + / + / . Taking the ℎ order macroscopic moment of the streaming Eq. ( 57) leads to the moment streaming equation

( +Δ ,x),( ) 1 ... = ( ,x-c i Δ ),( ) 1 ... . (59)
Similarly, Eq. ( 50) can be recast into a moment equation

( +Δ ,x),( ) 1 ... = ,( )
1 ... ( + Δ , x) + Δ 2 1 ( +Δ ,x),( ) 1 ... -( +Δ ,x),( ) 1 ... . (60)
Combining Eqs. ( 59) and ( 60) nally leads to ,( )

1 ... ( + Δ , x) = ( ,x-c i Δ ),( ) 1 ... - Δ 2 1 ( +Δ ,x),( ) 1 ... -( +Δ ,x),( ) 1 ... . ( 61
)
This equation is the update rule for ( + Δ ) moments. It is nothing but the LBM numerical scheme written explicitly for the ℎ order moment ( +Δ ,x),( ) 1 ...

. We shall now Taylor expand ( ,x-c i Δ ),( ) 1 ... , using Eq. ( 58),

( ,x-c i Δ ),( ) 1 ... = 1 ... 1 + ∞ =1 (-Δ ) ! + + ( , x) . (62)
Using the fact that the discrete velocities are xed leads to

( ,x-c i Δ ),( ) 1 ... = ( ,x),( ) 1 ... + ∞ =1 (-Δ ) ! + ( ,x),( + ) 1 ... + . ( 63
)
It is shown in Appendix A that inserting Eq. ( 63) into Eq. ( 61) eventually leads to a secondorder accurate Crank-Nicolson scheme for the continuous equation:

,( )

1 ... + ,( +1) 1 ... +1 +1 = - 1 ,( ) 1 ... + ,( ) 1 ... + O (Δ 2 ) . (64)
Note that this equation is only relevant for independent equations, the isotropy defect (47) leaving higher orders redundant (see Table III).

Because of the moments cascade,

( ,x),( +1) 1 ...
is always the ux of ( ,x),( ) 1 ...

, allowing moment

computation to be algorithmically explicit as long as --1

,( ) 1 ... + ,( ) 1 ... = 0. However, when the collision kernel --1 ,( ) 1 ... or a force term + ,( ) 1 ... is non-zero one has to solve implicitly Eq. ( 61) in order to get a second-order accurate scheme. This layered structure between orders shows that non-equilibrium moments follow their own evolution equation, they are not algebraically enslaved to lower order moments as suggested by the Chapman-Enskog expansion through the scale separation hypothesis. Mandatory conditions for the consistency of LB schemes as NSF solvers will be discussed later for speci c kernels.

LBM can be seen as a very smart change of variable from macroscopic moments to distribution functions. First the macroscopic information is stored inside the Hermite basis through ( rst change of variable from ( ) to ). Then the transport is performed in the distribution space during the streaming, followed by the macroscopic reconstruction that lters only the relevant information for each macroscopic moment ( ) (second change or variable from to ( ) ), prompting us to draw a parallel between classical CFD and LBM, which is now only seen as a numerical macroscopic solver. This shows that even if the CE expansion leads to the desired set of macroscopic equation it does not leads to the necessary result that would allow us to fully understand what solves LBMs, namely Eq. ( 64), a Grad-system of equations.

In the following we will illustrate the Taylor expansion on three di erent LBMs, namely the athermal, density-based recursive regularized (RR-) and pressure-based recursive regularized (RR-) models. We shall also demonstrate that although Eq. ( 64) could have been easily guessed from Eq. ( 42), spurious terms will appear during the coupling between LB algorithm and other numerical schemes. Because these terms are purely stemming from the numerical coupling they are out of the scope of the CE expansion.

V. ATHERMAL LBM

The last section showed that each macroscopic moment follows its own evolution equation (64), advocating a change of paradigm. Instead of considering LBM as a kinetic solver let us consider it as a Grad-solver for an extended set of thermo-hydrodynamic equations. Some of them are desired conservation laws such as mass and momentum conservation, others corresponds to higher order equations in the nite hierarchy of equations related to the lattice . Therefore, and all other populations previously de ned lose their kinetic meaning and are now merely seen as temporary variables in the macroscopic CFD solver known as "LBM".

A. Athermal numerical scheme

For the sake of clarity we rst apply the proposed Taylor expansion on the classical athermal LBM on standard lattices [START_REF] Krüger | The lattice boltzmann method[END_REF] with a force term speci cally designed to get rid of the well known O (Ma 3 ) error detailed in Table III. This model, traditionally said to be athermal in the LB literature, is often used in practice to solve isothermal ows. Following the new paradigm we de ne the initial solution simply by initial macroscopic elds ( , x), ( , x).From this initial condition we would like to nd a LB algorithm that predicts ( + Δ , x) and ( + Δ , x) following an approximate Navier-Stokes system, hopefully matching Eqs. (1,2). Let us now detail step-by-step the algorithm proposed in Table IV, applied to the classical athermal LBM.

• Step 1 : Equilibrium construction Because we restrict ourselves to standard lattices, some third-order Hermite polynomials H [START_REF] Cercignani | The boltzmann equation[END_REF] do not belong to the Hermite basis. For this reason we do not expand further than the third order, isotropy defects being corrected by an appropriate force term. The equilibrium reads = H (0) + H [START_REF] Hilbert | Mathematical problems[END_REF] 2 + H [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] 2 [START_REF] Saint-Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF] [ ] + H [START_REF] Cercignani | The boltzmann equation[END_REF] 6 [START_REF] Villani | Limites hydrodynamiques de l'équation de boltzmann[END_REF] [ ] .

•

Step 2 : Force construction The forcing population is extended to second order, ≡ H (0) ,(0) + H [START_REF] Hilbert | Mathematical problems[END_REF] 2

,(1) + H (2) 2 4 , (2)
,

with its Hermite moments de ned as

,(2) = - , (3)
+ 2 2 3 + F + F - , (67)
, [START_REF] Hilbert | Mathematical problems[END_REF] = F , (68)

,(0) = , (69)
with , [START_REF] Cercignani | The boltzmann equation[END_REF] the isotropy defect of the equilibrium population, related to the lattice and the particular equilibrium function. Remember that and F correspond respectively to a mass source and volume force in the macroscopic equations (1, 2). Now that macroscopic quantities, namely mass, velocity and stress-tensor ( + Δ , x), ( + Δ , x) and ,(2) ( + Δ , x) have been explicitly updated let us analyze the equivalent continuous equations of System (74, 75, 77) and compare it with the target set of equations Eqs. (1, 2) with = 2 (5).

Using the continuous limit Eq. ( 64) of the LB scheme leads to an extended Grad-system that conditionally approximates the athermal Navier-Stokes:

+ = (78) + + + ,(2) = F , (79)
,

+ ,(3) - ,(3) = - 1 ,(2) + ,(2) . ( (2)
) 80
where the lattice-dependent isotropy defect , [START_REF] Cercignani | The boltzmann equation[END_REF] can be found in Table III.

From this system we can infer that -, [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] is the e ective stress tensor. We also see that contrary to the usual Navier-Stokes equations this system has an evolution equation for the stresstensor, Eq. ( 80). This evolution equation involves ,(3) through , [START_REF] Cercignani | The boltzmann equation[END_REF] , which depends on higher order contributions and is assumed negligible by the CE expansion. The term ,(2) / hidden inside Eq. ( 80) can be replaced using Eq. ( 4). Using the second order moment of the equilibrium population (67,69) in Eq. ( 80) and the athermal equation of state ( 5) nally leads to

- ,(2) = 2 + - 2 3
+ , (2)
+ ,(3) - , (3)
- , (2)
+ , (2)
(81) with 2 = obtained by identi cation with the usual de nition of the stress tensor (3). Note

that the e ect of the collision kernel is entirely hidden inside , [START_REF] Cercignani | The boltzmann equation[END_REF] . When using a simple BGK collision step ,(3) is enslaved to higher order unphysical contributions because of the lack of isotropy that leads to an under-resolved nite hierarchy of equations in the velocity space.

This last equation is not algebraic as the truncated CE expansion asserts but rather an evolution equation for the unknown , [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] .

The next step is to demonstrate in which cases the Lattice Boltzmann stress-tensor equation is consistent with a Navier-Stokes stress tensor (80)

- ,(2) ≈ T = + - 2 3 . ( 82
)
To that end, let us nondimensionalize Eq. ( 81). First we need to identify what is the shortest physical time scale , corresponding to the fastest and dominant physical phenomenon. Depending on the situation, mainly two relevant candidates exist: the viscous timescale = 2 0 / and the convective timescale = 0 / 0 .

If the shortest timescale is , then the appropriate nondimensionalization reads

= 1 * , = 1 0 * , (83)
,(2) = 0 * , , (2)
, ,

= 0 * , , (3)
, (84)

= 0 * , = 0 * , = 0 * . ( 85
)
where * superscript quantities are O (1) and non-dimensional. Applying this change of variable, - * ,

,(2) = T * + 0 0 2 0 0 1 * * , , (3)
- * , ,

* + Ma 2 Re 1 * * * , , (3)
* + * * , , (2)
* + 0 2 1 * * , , (2)
* , ( (2)
) 86
where the Reynolds number Re = / and athermal Mach number

Ma = 0 / , (87)
have been used. This implicitly means that in the athermal case Ma is enslaved to the CFL number 24 because = Δ /( √ 3 Δ ), leading to

CFL = 0 + Δ /Δ = Ma + 1 √ 3 . ( 88
)
Note that the stability criterion CFL ≤ 1 boils down to the usual athermal Mach limit Ma ≤ √ 3 -1 ≈ 0.732, which is consistent with previous studies [START_REF] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF][START_REF] Wissocq | Linear stability of athermal regularized lattice boltzmann methods[END_REF] . If the convective scaling is chosen the stress-tensor becomes

- * , , (2)
= T * + O 0 0 2 0 0 + O Ma 2 Re . (89)
The condition of consistency on LBMs is that prefactors should be negligible quantities to verify our initial assumption Eq. ( 82). The last term of the right hand side can be neglected for the di usive and convective timescale respectively if Ma 2 /Re 2 and Ma 2 /Re are small enough. On the other hand the factor 0 / 0 2 0 0 corresponding to the scaling between the third order non-equilibrium and second order non-equilibrium is not generally trivial. Because the isotropy defect completely modi es the convective term in the evolution equation of , [START_REF] Cercignani | The boltzmann equation[END_REF] , it leads to an evolution equation whose physical meaning is unclear. More importantly, this unphysical evolution equation can be a source of instabilities [START_REF] Masset | Linear hydrodynamics and stability of the discrete velocity boltzmann equations[END_REF] . However, when the recursive regularized kernel [START_REF] Coreixas | Recursive regularization step for high-order lattice boltzmann methods[END_REF] is used, 0 is known and the stress-tensor becomes

- * , , (2)
= T * + O Ma 2 Re . (90)
One sees here that the usual "small Knudsen" assumption is not even su cient because Kn.Ma terms appeared in the scaling analysis. To get back the proper NSF stress-tensor one have to carefully analyze one by one each of these spurious terms. Additionally the scaling 0 is related to the heat ux, suggesting that the Prandtl number Pr should also intervene in the range of consistency of LBMs.

The Taylor expansion showed us that the consistency condition was not as simple as the CE expansion suggests. The small Knudsen assumption is not enough and both the choice of the lattice and the collision kernel changes the consistency defect in Eq. (90). More discrete velocities mean that the isotropy defect is pushed away from the Navier-Stokes and the stress-tensor equation, but it also means that more unphysical equations are also solved. Those equations are likely to modify or even undermine the validity of the solution, such as for Burnett-like systems.

VI. COLLISION KERNEL

Now that we analyzed a simple athermal BGK LB scheme let us discuss the collision kernel by reviewing a sample of techniques that can be used to increase robustness of LBMs. From the Taylor expansion we've seen that mass and momentum conservation equations were correctly discretized up to O (Δ 2 ) errors by the LB scheme. On the other hand, the system is not closed through an algebraic constitutive equation as in usual CFD solvers. Instead we inherit Eq (81)

from the hierarchy of moments. Rearranging its terms leads to ,

+ , (2)
= -

1 ,(2) + T , (91)
where T (3) is the physical stress tensor. Higher order contributions , [START_REF] Cercignani | The boltzmann equation[END_REF] , as already discussed, do not necessarily match a physical behavior, especially for standard lattices because of isotropy defects. Therefore, a question that could drive us towards the use of a particular kernel is a correct modeling of the stress tensor by Eq. ( 91). Being the only moment that does not appear in the hydrodynamic equations Eqs. (78-79) the non-equilibrium tensor ,(3) is our only degree of freedom to modify the modeling of Eq. ( 91) towards a physically meaningful stress tensor transport equation. For example with a BGK collision operator, ,

purely stems from higher order non-hydrodynamic equations. In this case, assuming that the lattice is large enough such that ,(3) = 0 (this property is only enforced on some non-diagonal components of [START_REF] Cercignani | The boltzmann equation[END_REF] for standard lattices) we can write for example the ,(2) evolution equation, ,

= -

1 ,(2) + T - , (3)
.

In the left hand side we recognize a transport term in the direction with an unexpected backward propagation -2 while the rst term in the right hand side is a relaxation term that steers the variable , [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF] towards the expected T . The second term in the right hand side is the coupling with higher-order non-hydrodynamic moments.

A. E ect of regularized kernels

A particularly e cient way to control the time evolution of the stress tensor was identi ed in the regularized and recursive regularized collision kernels. The simple regularization [START_REF] Latt | Lattice boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Mattila | High-order regularization in lattice-boltzmann equations[END_REF] simply discards , [START_REF] Cercignani | The boltzmann equation[END_REF] during the collision leading to a ltered non-equilibrium population neq =

H (2) 2 4 neq , (2)
,

= -eq + Δ

2

(94) allowing those ltered moments to be compartmentalized from the hydrodynamic moments, and e ectively canceling the last term of Eq. ( 92).

+ neq , (2)
,

leading to a new evolution equation, ,

+ , (2)
+ 2

,(2) = - 1 ,(2) + T , (96)
that now exhibits a forward transport, which may explain why high velocity ows are more stable [START_REF] Coreixas | Recursive regularization step for high-order lattice boltzmann methods[END_REF] with this regularized collision kernel.

B. E ect of the Hybrid Recursive Regularization

An extension of the recursive regularization was developed [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-grosskrook collision model for lattice boltzmann method-based large eddy simulation[END_REF] by introducing ∈ [0, 1] into the non-equilibrium reconstruction Eq. (95) as neq

,(2) = =1 - + Δ 2 -(1 -) 2 + - 2 3 , (97)
where the last term is evaluated from a nite di erence scheme. It has been shown [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-grosskrook collision model for lattice boltzmann method-based large eddy simulation[END_REF] that this modi cation leads, for < 1, to the introduction of a numerical hyperviscosity in the momentum equation.

In light of the previous Section, an alternative explanation for the enhanced stability is that neq ,(2) may deviate from its target value T . Using 0 < < 1 (resp. = 0) as a weighting parameter is equivalent to a partial (resp. total) reset of neq ,(2) to its fully relaxed value T at the end of each time step, leading to a stronger steering of neq ,(2) towards T by the resulting LB scheme.

C. Trace of the stress tensor

The pressure work is of paramount importance in compressible ows and was already emphasized as a major source of instabilities for thermal LBMs [START_REF] Renard | A linear stability analysis of compressible hybrid lattice boltzmann methods[END_REF] . Because the non-equilibrium moment neq , (2) is the e ective stress tensor in LBMs, any spurious term appearing on its trace will act as a spurious pressure in momentum equation. Therefore, arti cially enforcing a traceless stress

tensor 7 , neq ,(2) = - - eq + Δ 2 , (98)
during the collision allows to get rid of this spurious pressure.

Because and are conserved moments in LBMs, the collision kernel has no direct impact on their numerical schemes. But we've seen in this section that a choice of collision kernel impacts the closure for the neq ,(2) equation. Therefore, a choice of collision kernel is a choice of closure for non-equilibrium moments. For this reason, we believe the present method to be a good candidate to analyze the LBM closure, and design future collision kernels.

VII. THERMAL RR-

Now that we analyzed both the athermal model and the collision kernel, let us analyze two thermal models using the recursive regularized collision kernel. These models are of particular interest to highlight advantages of the Taylor expansion over CE as it has both a coupling with another algorithm and a complex collision kernel.

Let us start with the recursive regularized pressure-based model (RR-), presented earlier this year [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible ows[END_REF] , to simulate both compressible ows and reactive ows [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF] .

The main feature in this model is that the equilibrium population is chosen so that the 0th order moment of the equilibrium population is a normalized pressure, the rest of the equilibrium is similar to the athermal equilibrium distribution Eq. (65). It leads to an unphysical mass conservation equation that is mended by a correcting term that explicitly uses the Taylor expansion, which is by itself a su cient reason to analyze this model through the scope of the Taylor expansion instead of the CE expansion.

Since standard lattices are not able to recover accurately the second-order moment corresponding to energy, the energy conservation ( 6) is solved separately, as φ in Step 6 of Table IV. Also note that the athermal equation of state is replaced, from here on, by the perfect gas equation of state (7).

A. RR-numerical scheme let us assume that we have a consistent and convergent explicit scheme for Eq. ( 6), expressed

as

[ ] ( + Δ , x) = [ ] ( , x) + Δ ( , x) + O (Δ ) , (99)
with ≥ 1 and an operator that represents everything but the time derivative in Eq. ( 6). We shall now detail the step-by-step algorithm leading to the RR-model.

•

Step 1 : Equilibrium construction The equilibrium is expanded as

= H (0) + H (1) 2 + H (2)
2 [START_REF] Saint-Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF] [ ] + H [START_REF] Cercignani | The boltzmann equation[END_REF] 6 [START_REF] Villani | Limites hydrodynamiques de l'équation de boltzmann[END_REF] [ ] .

(100)

Note the slight modi cation in the 0th order, with this model the density is replaced by a normalized pressure

= / 2 . (101)
•

Step 2 : Force construction The forcing population is simply extended to second order, ≡ H (0) ,(0) + H [START_REF] Hilbert | Mathematical problems[END_REF] 2

,(1) + H (2) 2 4 , (2)
,

with its Hermite moments de ned as

,(2) = 2 (1 -) + 2 (1 -) + 2 2 3 -2 (1 -) - , (3)
+ F + F - , (103)
, [START_REF] Hilbert | Mathematical problems[END_REF] = F ,

,(0) = .

(105)

• Step 3 : Non-equilibrium construction Though this particular pressure-based coupling does not presuppose the use of any particular collision kernel we chose the recursive regularized collision operator [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible ows[END_REF] . The non-equilibrium population is then reconstructed following Eq. ( 95).

the equivalent continuous equations of the system Eqs. (109,110,111,99) and compare it with the target set of equations Eqs. (1,2,6). The Taylor expansion leads to the system of equations

+ = , (112)
+ + + , (2)
= F ,

,

+ , (2)
-,

= - 1 ,( ) + ,( ) , (3)
+ ( + ) + + , (114)
= F + .

We recognize mass, momentum and energy conservation along with an evolution equation for the stress tensor, let us analyze this equation. Using the kinetic tensor equation Eq. ( 4), the force terms Eqs. (103,105) along with the 2nd and 3rd order macroscopic moments of the equilibrium population lead to the stress tensor equation,

- , (2) = 2
+ - 2 3
+ , (2)
+ ,(3) - , (3)
-,

+ , (2)
,

with 2 = . This is the exact same evolution equation as in the athermal case. Because equilibrium moments are enslaved to low order macroscopic moments the only place where the higher order moments have an impact in thermo-hydrodynamic equations is through , [START_REF] Cercignani | The boltzmann equation[END_REF] .

Fortunately here we used a hybrid recursive regularized collision kernel so that even if higher order macroscopic equations are unphysical their e ect on low order macroscopic equations is explicitly ltered during non-equilibrium reconstruction [START_REF] Coreixas | Recursive regularization step for high-order lattice boltzmann methods[END_REF][START_REF] Wissocq | Linear stability of athermal regularized lattice boltzmann methods[END_REF] . To nd out in which cases the stress-tensor evolution equation is su ciently close to its expected value we use the same nondimensionalization as in the athermal case, except that 0 = 0 0 is known because of the RR collision kernel. Also note that is not anymore related to the physical sound speed, therefore changing the de nition of the Mach number to Ma = 0 / 0 . This leads to

- * , , (2)
= T * + Ma 2 Re 0 2 1 * * , ,(3) - * , , (3) *
- Ma 2 Re 0 2 1 * * * , , (2)
* + * * , , (2)
* + 0 2 1 * * , , (2)
* . ( 117
) CFL = 0 + 0 Δ /Δ , (118)
and reminding that in the acoustic scaling = Δ /( √ 3 Δ ) we can see that

0 2 = 3CFL 2 (Ma + 1) 2 . ( 119
)
By assuming that our fastest relevant phenomenon is the convection, such that = 0 / 0 , we can rearrange Eq. ( 117) into

- * , , (2)
= T * + O Ma 2 CFL 2 Re (Ma + 1) 2 , (120)
which is the consistency condition of the pressure-based model [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible ows[END_REF] . This very simple expression was obtainable only because we used a collision kernel in which the scaling of third order nonequilibrium moments is algebraically enslaved to variables whose scaling is known (see Eq. ( 95)).

This equation shows a very special feature of LBMs: the discretized continuous equations directly depends on the CFL number.

VIII. THERMAL RR-

This scheme have been thoroughly used in the literature to perform simulations on compressible [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible ows[END_REF][START_REF] Renard | Improved compressible hybrid lattice boltzmann method on standard lattice for subsonic and supersonic ows[END_REF] , combustion [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF]69 . Mass and momentum are solved by the LB algorithm while the energy and potential species conservation are solved by a FD scheme, as already presented in the last section.

A. RR-numerical scheme

Contrarily to the RR-model this scheme encompasses the thermal e ects coming from the FD energy equation inside the 2nd and 3rd order equilibrium moments, leading to a more complex equilibrium distribution function.

•

Step 1 : Equilibrium construction The equilibrium is expanded as

= H (0) + H (1) 2 + H (2) 2 4 [ + 2 ( -1) ] + H (3) 6 6 { + 2 ( -1) [ + + ]} , (121)
where keeps its de nition (101). Now that macroscopic quantities, namely mass, velocity, stress-tensor and temperature, ( +

Δ , x), ( +Δ , x), , (2)
( +Δ , x) and ( +Δ , x) have been explicitly updated let us analyze the equivalent continuous equations of the system Eqs. (128,129,130,99) and compare it with the target set of equations Eqs. (1,2,6). The Taylor expansion leads to the system of equations

+ = , (131)
+ + + , (2)
= F ,

,

+ ,(3) - , (2)
= - 1 ,( ) + ,( ) , (3)
+ ( + ) + + , (133)
= F + .

We shall now deduce the stress tensor evolution equation. First, using Eqs. (131,132) to isolate the pressure time derivative / inside Eq. ( 134) and assuming that the adiabatic exponent is a constant leads to the pressure equation,

( + ) + + + , (2)
= +

-,

+ , (2)
+ - - , (2)
(136) with = . Here we end up with a slightly di erent evolution equation for the stress tensor when compared to the pressure-based and athermal cases. Again the higher order unphysical moments are ltered out because , [START_REF] Cercignani | The boltzmann equation[END_REF] was reconstructed following Eq. (95).

To allow a fair comparison with the RR-model we use the previously de ned nondimensionalization. This leads to

- * , , (2)
= * , + Ma 2 Re

1 * * * , , (3)
- * , ,

IX. ILLUSTRATION: SPHERICAL FLAME

In order to illustrate that di erent consistency defects can lead to very di erent results a circular laminar ame is simulated with both the density-based (Sec. VIII) and the pressure-based (Sec. VII) models, as adapted for reactive ows [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF] . This test case was found to particularly aggravate spurious currents e ects because of the presence of heat-release along a curved density gradient. This two dimensional test is setup as follows: the domain size is chosen as 4cm×4cm, with a grid spacing = 0.2mm. The domain is initialized with fresh stoichiometric mixture of propane and air at 300K and 1atm and a spot of corresponding burnt gases in its center. The thermodynamic, transport and chemical closures applied in the current work are identical to those in [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF] , which correspond to a laminar ame speed ≈ 0.25m/s and a thermal ame thickness ≈0.6mm. It is worth noting that the grid is quite coarse, corresponding to around three points in one thermal ame thickness. The boundary conditions are set to characteristic in/out boundaries with a far-eld pressure ∞ = 1atm.

The ow elds of two very early time steps of simulation are presented in Fig. 1. It can be seen that the density-based scheme creates spurious currents near the ame front which leads to unstable simulation, whereas the pressure-based scheme presents the expected behavior.

The temperature, heat release rate (HR) and velocity elds from the simulation pressure-based are given at a later instant in Fig. 2. Pro les are seen to remain circular even when approaching the boundary, showing both the isotropy of the scheme as well as the e ectiveness of the characteristic boundary.

This comparison between density-based and pressure-based formulations illustrates how spurious currents can be triggered by defects of consistency of the stress-tensor. • bounded consistency error with increasing Mach number and vanishing error with decreasing CFL number, for the pressure-based model (iii).

The Taylor analysis presented here therefore clearly favors the pressure-based model, while the Chapman-Enskog expansion was unable to discriminate the two models.

rigorously speaking this 2nd order accuracy stands for each of the moments associated to the DdQq lattice only in the BGK case. In other words Eq. (A4) is strictly veri ed only for BGK kernel. Meaning that depending on the chosen collision kernel the schemes associated to higher order moments than (0) and [START_REF] Hilbert | Mathematical problems[END_REF] may not exhibit a 2nd order accuracy [START_REF] Wissocq | Linear stability of athermal regularized lattice boltzmann methods[END_REF] .

1 .

 1 .. between continuous and discrete moments,

  with / = 1/( -1). Then injecting Eqs. (4,135,123-125) into Eqs. (133) nally leads to -

  is dominant, = 0 / 0 , Eq. (137) can be recast into should be compared to Eq. (120). While the consistency defect of the densitybased model increase with Ma it remains bounded for the pressure-based. Another interesting di erence between density-based on one side and both athermal and pressure-based on the other side is that the former model does not have a decreasing consistency defect when the CFL is decreased, an evident limitation for simulating compressible ows.

Figure 1 :

 1 Figure 1: Streamlines of the 2D circular ame simulation colored by velocity magnitude (in m/s). Left column: pressure-based model 18 , right column: density-based model 17 . Note the very di erent ranges of velocity magnitude from the two methods. The yellow contour is the heat release rate peak indicating the ame front.

Table I :

 I Chapman-Enskog expansions in the lattice-Boltzmann literature.

Table II :

 II Lattice-Boltzmann notations

	Notation	Representation	Equation
		Total population	+
		Equilibrium population	48
		Non-equilibrium population	49
	,	Modi ed populations for 2 order accuracy	50, 53
		Population after collision	54
		Forcing term	52
	H ( )	Hermite polynomials	43

,( ) 1 ...

Table IV :

 IV Generic LBM structure.

	Step	Description	Input		Output	Equation(s)
	1	Equilibrium update	macroscopic		48
	2	Force update	macroscopic	( , x)	52
	3	Non-equilibrium update	, ( , x)		53 (BGK)
	4	Collision	,	, ( , x)		56
	5	Streaming	( , x -c i Δ )	( + Δ , x)	57
	6	Coupling with another scheme	φ( , x)		φ( + Δ , x)	optional
	7	Macroscopic quantities update			macroscopic	61

Table V :

 V Summary of consistency errors for the stress tensor, for the three models studied.Here Ma is the lattice-speed based Mach number, Eq. (87)

	Model	Stress tensor error	Equation
	Athermal model (Section V)	O	2 Re Ma			(90)
	Compressible pressure-based (Section VII)	O	Ma 2 CFL 2 Re (Ma + 1) 2	(120)
	Compressible density-based (Section VIII)	O	Ma 2 Re	+ O	1 Re Pr	(138)

ACKNOWLEDGEMENTS

We thank Gauthier Wissocq and Pierre Sagaut for inspiring discussions. We also acknowledge support from Labex MEC (ANR-10-LABX-0092), the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir" and the French Space agency (CNES) for supporting Song Zhao at M2P2.

DATA AVAILABILITY

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

ulation is obtained as

Otherwise, following a particular choice of collision kernel construct ˜ and set

•

Step 4 : Collision With , and built in the previous steps, compute the collided population as

•

Step 5 : Streaming Transport the populations according to

• Step 6 : Coupling update No coupling is necessary in the athermal case, but this stage can be used to transport passive scalars [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF][START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b uids[END_REF][START_REF] Küng | Comparison of passive scalar transport models coupled with the lattice boltzmann method[END_REF] .

• Step 7 : Update macroscopic variables Using the macroscopic update rule Eq. ( 61) for = 0, 1, 2 respectively leads to

, [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF]

By splitting ,(2) ( +Δ , x) into its equilibrium and non-equilibrium parts, the above leads to the stress-tensor scheme

.

(77)

•

Step 5 : Streaming Shift the populations according to

• Step 6 : Coupling update Here simply solve the explicit numerical scheme Eq. (99) to get [ ] ( + Δ , x).

• Step 7 : Update macroscopic When applied to = 0 the macroscopic update rule Eq. ( 61) usually gives us the updated density ( + Δ , x). But here because we modi ed the equilibrium distribution it leads to an updated pressure ( * / 2 )( + Δ , x) such that

which is equivalent to an unphysical equation for the pressure, ( * / 2 )/ + / = , see Eq. ( 64) for = 0. To recover the correct mass conservation equation it is mandatory to modify the 0th order update rule such as :

The update rule for = 1, 2 remains unchanged,

.

Then compute ( + Δ , x), from [ ] ( + Δ , x), ( + Δ , x) and ( + Δ , x).

B. Continuous equivalent equations

Now that macroscopic quantities, namely mass, velocity, stress-tensor and temperature, ( +

( +Δ , x) and ( +Δ , x) have been explicitly updated let us analyze

,

with its Hermite moments de ned as

,(0) = .

(125)

• Step 3 : Non-equilibrium construction To allow a fair comparison with the pressurebased model, let us use again the recursive regularized collision operator, Eq. (95).

•

Step 4 : Collision Thanks to previous steps , and have been built, compute the collided population such that

•

Step 5 : Streaming Shift the populations according to

• Step 6 : Coupling update Here simply solve the explicit numerical scheme Eq. (99) to get [ ] ( + Δ , x).

• Step 7 : Update macroscopic Using the macroscopic update rule Eq. ( 60) for = 0, 1, 2 respectively leads to

.

Then compute ( + Δ , x), from [ ] ( + Δ , x), ( + Δ , x) and ( + Δ , x).

X. CONCLUSION

We have presented a new systematic method to analyze Lattice-Boltzmann models. Based on a

Taylor-Expansion, it allows to derive the set of macroscopic equations consistent with the model.

The analysis only requires the assumption of a su ciently small time discretization Δ → 0, as is usual for consistency studies.

A major advantage of the method over the traditional Chapman-Enskog framework is that it allows to study more carefully numerical errors. In light of this Taylor-Expansion, collision in Tab. V show that, for the athermal model (i), the consistency error decreases with the CFL and increases with the Mach number. For compressible models, two distinct behaviors were found:

• increasing consistency error with increasing Mach number, and bounded error with decreasing CFL number, for the density-based model (ii),

DECLARATION OF INTERESTS

The authors report no con ict of interest.

To study the order of accuracy let us use Eq. ( 63) evaluated for the ( +1)-order moment along with Eq. ( 59) to get

which can be injected back into Eq. ( 63), then using Eq. ( 59) leads to

On the other hand, Eqs. (52,55,50) leads to

At this point it is curious to note that the collision forcing -1 ,( ) 1 ... and the external forcing ,( ) 1 ... have the exact same treatment in the algorithm. Therefore the collision discretization is nothing else than a particular Guo forcing applied to a collisionless discretized BE. Also note that (A5)

Injecting Eq. (A5) and Eq. (A4) respectively in the right hand side and left hand side of Eq. (A2) leads to the general second order numerical scheme :


M Tayyab

S Zhao

Pierre Boivin

Lattice-Boltzmann modeling of a turbulent bluff

come

Lattice-Boltzmann modelling of a turbulent bluff-body stabilized flame M. Tayyab, 1 S. Zhao, [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF][START_REF] Löhner | Towards overcoming the LES crisis[END_REF] and P. Boivin 1, a) This letter reports the first large eddy simulation of a turbulent flame using a Lattice-Boltzmann model. To that end, simulation of a bluff-body stabilized propane-air flame is carried out, showing an agreement similar to those available in the literature. Computational costs are also reported, indicating that Lattice-Boltzmann modelling of reactive flows is competitive, with around 1000cpuh required to simulate one residence time in the 1,5m burner.

Lattice-Boltzmann methods (LBM) are an active topic of investigation in the computational fluid dynamics field. In the past decade, they led to unprecedented results in the field of low-Mach, isothermal, aerodynamics and aeroacoustics. The two main ingredients of this success can be summarized as follows. (i) LBM is highly scalable and HPC friendly by nature, due to its stream and collide algorithm [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF][START_REF] Löhner | Towards overcoming the LES crisis[END_REF] . In practice, CPU costs are reportedly divided by 5 to 10 compared to conventional Navier-Stokes solvers. (ii) Besides being secondorder accurate in time and space [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] , LBM presents excellent dissipation properties for pressure waves, comparable with high-order methods [START_REF] Marié | Comparison between lattice boltzmann method and navier-stokes high order schemes for computational aeroacoustics[END_REF] . Owing to this success, there is today an unprecedented effort in extending those methods to tackle problems beyond low-Mach athermal flows: compressible flows [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Latt | Efficient supersonic flow simulations using lattice boltzmann methods based on numerical equilibria[END_REF] , reactive flows [START_REF] Hosseini | Hybrid lattice boltzmann-finite difference model for low mach number combustion simulation[END_REF][START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Feng | A lattice-boltzmann model for lowmach reactive flows[END_REF] , multiphase flows [START_REF] Li | Lattice boltzmann methods for multiphase flow and phase-change heat transfer[END_REF] .

Our group recently presented a pressure-based LB model designed to cope with both high Mach flows [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] and reactive flows [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF] , emerging from earlier studies [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Feng | A lattice-boltzmann model for lowmach reactive flows[END_REF] . These studies validated the behavior of the model on standard combustion test cases, including, e.g. laminar flame propagation [START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Feng | A lattice-boltzmann model for lowmach reactive flows[END_REF] , counterflow diffusion flame [START_REF] Feng | A lattice-boltzmann model for lowmach reactive flows[END_REF] , vortex/flame interactions [START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF] , diffusion flame stabilized in double shear layer [START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF] , and thermo-diffusive instabilities [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF] . This letter aims at proving the method to be also suitable for more complex configurations, through investigation of a premixed turbulent flame. To that end, we selected the so-called VOLVO burner, for which a rich experimental database exists [START_REF] Sjunnesson | Lda measurements of velocities and turbulence in a bluff body stabilized flame[END_REF][START_REF] Sjunnesson | Cars measurements and visualization of reacting flows in a bluff body stabilized flame[END_REF] . The test case, consisting of a premixed turbulent flame stabilized behind a bluff-body, has been extensively used to test and tune new turbulent combustion models [START_REF] Ma | A posteriori testing of algebraic flame surface density models for les[END_REF][START_REF] Maestro | Les reliability of the volvo bluff-body stabilized ame dynamics[END_REF] .

From the numerical modelling standpoint, many studies have been conducted on the VOLVO burner. Because of its dimensions (1,5m in length), Direct Numerical Simulation (DNS) is out of reach. Reynolds Averaged Navier-Stokes (RANS) simulations were carried out by many researchers [START_REF] Bai | Modelling of turbulent reacting flows past a bluff body: assessment of accuracy and efficiency[END_REF] , but the strategy has been reported to create divergences in the solution [START_REF] Porumbel | Large eddy simulation of bluff body stabilized premixed flame[END_REF][START_REF] Zou | Large eddy simulation of bluff-body stabilized turbulent premixed flame based on openfoam[END_REF] . To capture the unsteady shedding and vortices breakdown created by the bluff-body, Large Eddy Simulation (LES) presents an attractive alternative, and was adopted in a large number of numerical studies [START_REF] Fureby | Large eddy simulation of reacting flows applied to bluff body stabilized flames[END_REF][START_REF] Zettervall | Large eddy simulation of a premixed bluff body stabilized flame using global and skeletal reaction mechanisms[END_REF][START_REF] Wu | Mvp-workshop contribution: Modeling of volvo bluff body flame experiment[END_REF][START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Cocks | Impact of numerics on the predictive capabilities of reacting flow les[END_REF][START_REF] Chakravarthy | Large eddy simulations of confined bluff body stabilized highly turbulent premixed flames. asme fluids engineering summer meeting paper no[END_REF] . To the authors' knowledge, however, this simulation was never attempted using LBM.

The numerical simulations are carried out with the PROLB a) Electronic mail: pierre.boivin@univ-amu.fr software using a pressure-based Lattice Boltzmann (LB) model [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] coupled with a Finite Differences (FD) solver.

In the LB solver, the probability density function f i (of finding gas particles at x with velocity c i ) is solved using a hybrid regularized collision model [START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[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 streaming and collision process can be expressed under discrete form as

f i (t + δt, x) = f col i (t, x -c i δt), f col i (t, x) = f eq i + 1 - δt τ f neq i + δt 2 F E i ( 1
)
where δt is the time-step, c i is the i th discrete velocity of the D3Q19 lattice set, f eq i and f neq i the equilibrium and nonequilibrium segments of the distribution function and F E i is the forcing term required to correctly recover the stress tensor [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF] . τ is a non-dimensional relaxation time depending on the dynamic viscosity ν as

τ = ν c 2 s + δt 2 (2)
where c s is the lattice sound speed. Full expressions for f eq i , f neq i , F E i and their relation with the macroscopic variables are provided as supplementary material. They correspond [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] to solving the following macroscopic equations

∂ ρ ∂t + ∂ ρu β ∂ x β = 0, , (3a)
∂ ρu α ∂t + ∂ ρu α u β + pδ αβ -T αβ ∂ x β = 0 , ( 3b
)
where ρ is the volume mass, u α is the local velocity vector and p is the pressure, obeying the multi-component perfect gas law. T αβ is the stress tensor,

T αβ = ρν ∂ u α ∂ x β + ∂ u β ∂ x α -δ αβ 2 3 ∂ u γ ∂ x γ . ( 4
)
The FD solver coupled with the above is responsible for solving species and energy conservation. For each species k,

ρ ∂Y k ∂t + ρu α ∂Y k ∂ x α = ∂ ∂ x α (-ρY k V k,α ) + ωk , (5)
where Y k is the mass fraction of species k, ωk is its net chemical production rate , and V k,α the diffusion velocity [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] . In this work, the energy conservation is solved in its enthalpy form 26

h = N ∑ k=1 h k Y k , h k = T T 0 C p,k (T )dT + ∆h 0 f ,k , (6)
where T and h are linked through NASA polynomials leading to the following enthalpy equation

ρ ∂ h ∂t + ρu α ∂ h ∂ x α = Dp Dt - ∂ q α ∂ x α +T αβ ∂ u α ∂ x β , (7)
where

Dp Dt = ∂ p ∂t + u α ∂ p ∂ x α is neglected. The heat flux q α reads q α = -λ ∂ T ∂ x α + ρ N ∑ k=1 h k Y k V k,α , (8)
with λ the thermal conductivity. Diffusion velocities are defined as in [START_REF] Tayyab | Hybrid regularized latticeboltzmann modelling of premixed and non-premixed combustion processes[END_REF] , using constant Schmidt number for each species, and a correction velocity to ensure mass conservation [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] .

The possible spatial discretization being much larger than the smallest flow structures, subgrid modelling is required. First, a turbulent model is implemented for the flow, coming down to modifying the dynamic viscosity in the momentum equation (3b) as ν → ν + ν t , where ν t is the turbulent Vre- mann model [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF] available in PROLB, a model validated for a wide range of applications covering external aerodynamics to urban flows [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann method-based large eddy simulation[END_REF][START_REF] Wilhelm | A new explicit algebraic wall model for les of turbulent flows under adverse pressure gradient[END_REF] . Turbulent Schmidt and Prandtl numbers are introduced [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] , to modify accordingly the transport properties in Eqs. (5, 8).

Lastly, the flame is thickened via the dynamic thickened flame large eddy simulation (TFLES) model [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] . To account for flame/turbulence interactions at the subgrid level, the model consists in modifying the flame structure locally in two steps:

Step 1: Thickening of the flame by a factor F,

Step 2: Applying an efficiency function E , as proposed by Charlette et al. [START_REF] Charlette | A power-law flame wrinkling model for les of premixed turbulent combustion part i: non-dynamic formulation and initial tests[END_REF][START_REF] Charlette | A power-law flame wrinkling model for les of premixed turbulent combustion part ii: dynamic formulation[END_REF] , to enhance the transport and kinetics where the flame is highly wrinkled.

In practice, applying this model corresponds to multiplying all transport terms by E F, and kinetic source terms by E /F; only where a flame sensor S is activated. The corresponding species and energy filtered equations are recalled in the supplementary material. S and E are obtained following Rochette's proposal [START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] . The VOLVO burner [START_REF] Sjunnesson | Lda measurements of velocities and turbulence in a bluff body stabilized flame[END_REF][START_REF] Sjunnesson | Cars measurements and visualization of reacting flows in a bluff body stabilized flame[END_REF] consists of a rectangular box of dimensions (L x , L y , L z ) = 1.55m × 0.12m × 0.24m (see the choice of axes in Fig. 1). In its center lies a bluff-body of equilateral triangular cross section of dimension D = 0.04m. The region upstream the flame holder consists of air and fuel injectors along with a honeycomb structure to obtain a turbulent premixed mixture of propane and air. Following other studies [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF][START_REF] Ghani | Les of longitudinal and transverse self-excited combustion instabilities in a bluff-body stabilized turbulent premixed flame[END_REF] , the numerical setup is simplified to a propaneair turbulent inlet (without injectors or honeycomb). Dimensions of the numerical burner are provided in Fig. 1. In the experiment, top and bottom walls (in the y direction) are water cooled, while side walls (in the z direction) are air cooled. RMS and mean values of velocity and temperatures are reported [START_REF] Sjunnesson | Cars measurements and visualization of reacting flows in a bluff body stabilized flame[END_REF] for various x-isolines on the centerline z = 0, identified in Fig. 1. Velocity profiles are also available along the centerline y = z = 0.

The operating points (equivalence ratio ϕ, bulk velocity U b , temperature T 0 and Reynolds number Re) listed in Tab. I are considered in the following, to allow comparison with available numerical studies [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] . These conditions were selected for no thermo-acoustic instability is expected [START_REF] Sjunnesson | Cars measurements and visualization of reacting flows in a bluff body stabilized flame[END_REF] , while some conditions are reportedly unstable [START_REF] Ghani | Les of longitudinal and transverse self-excited combustion instabilities in a bluff-body stabilized turbulent premixed flame[END_REF] . Numerical boundary conditions are specified as follow.

Turbulent velocity inlet: prescribed velocity U b with a turbulence intensity of 8%, generated using the Random Flow Generation (RFG) technique [START_REF] Smirnov | Random flow generation technique for large eddy simulations and particle-dynamics modeling[END_REF] .

No-slip adiabatic walls: top, bottom and bluff-body walls (along the y-direction).

Periodic conditions: in the z-direction. Note that no-slip adiabatic conditions were also considered for these walls, with results presented in Fig. 2.f.

Non-reflective characteristic boundary condition: outlet, with a target pressure of 1atm.

All thermo-chemical properties were taken from earlier studies [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] : a two/step propane-air mechanism is considered, with constant Prandtl and species Schmidt numbers, with details provided as supplementary material.

For both the hot and cold configurations, two different meshes (Case A, Case B) are considered.

Case A: Uniform mesh, with spatial resolution δ x = 2.0 mm, resulting in 5.5 × 10 6 grid points. Case B: Starting from the case A mesh, the region of interest is refined with δ x = 1.0mm (see the dotted lines in Fig. 1), resulting in 17.5 × 10 6 points.

The time step used is δt = 1.05 × 10 -6 s for the δ x = 2mm region (with the time-step halved in the refined region). For each simulation, three convective times t c = L x /U b are simulated to establish the flow, and statistics (mean and RMS quantities) are obtained over another four convective times. Figure 3 shows qualitatively the Von-Karman vortex shedding structure in the cold case for both Case A and Case B, in the form of Q-criterion. As expected, smaller structures are obtained for Case B. Figure 2 presents mean and RMS values compared quantitatively with the experimental results [START_REF] Sjunnesson | Lda measurements of velocities and turbulence in a bluff body stabilized flame[END_REF] for axial and transverse velocities at three different positions (x 1 , x 3 , x 5 , see Fig. 1).

Results are in good agreement with the experiments and the two mesh levels (Case A, Case B) produce similar results. Furthermore, we compare the axial mean velocity along xaxis in Fig. 2f. It can be viewed, that the refined mesh (Case B) captures the recirculation zone well. However, the results deviate from the experimental data downstream of the bluff body. This result is believed to be due to the periodic boundaries used on the sides of the wall [START_REF] Wu | Mvp-workshop contribution: Modeling of volvo bluff body flame experiment[END_REF] . In order to check this hypothesis, a simulation is carried out using wall boundaries on the sides at coarser grid δ x = 2.0mm and compared with Let us now present the reactive case, with an inlet bulk velocity of 17.3m/s. Following Rochette et al. [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF] , the parameters related to TFLES model, effective filter and averaging filter widths, are set to ∆ = 1.4Fδ L and ∆ avg = 2.5∆, respectively. Here, δ L is laminar flame thickness. Note that results report- for U, V and U ′ V ′ , along with temperature profiles measured at x 4 , x ′ 4 , x ′ 5 . The results are also compared with few available curves in [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF] .

It can be viewed from the figures that an overall agreement with the experimental and reference data is achieved, we also observe discrepancies between Case A and Case B results. This disparity is then highlighted in Fig. 4.f, where we compare the mean axial velocity profiles along x-axis. The results for finer grid have better agreement just after the bluff body. The recirculation zone is well captured, but the results deviate from the experimental data farther from the flame holder. As for the cold flow, this may be partly due to the use of periodic boundaries. Nonetheless, it is also observed that the computed recirculation zone is too short, suggesting that the combustion is too fast. This behavior is known and is reportedly due to the use of a simplified chemistry [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF] (see, comparison between detailed and simplified chemistry in Fig. 4.f), an assumption to be lifted in a future study. Furthermore, the errors observed between reference solution, which uses the same TFLES model, may be associated to the interaction of the LB method with the sub-grid models and its dependency FIG. 6. Temperature profiles at (x 4 , x ′ 4 , x ′ 5 -left to right). Case A (solid line), Case B (dashed line), experimental data (black solid square) and reference data -two/step chemistry 22 (dotted line) on the grid size as reported in [START_REF] Cocks | Impact of numerics on the predictive capabilities of reacting flow les[END_REF] . Hence, in order to characterize such discrepancies a thorough research is required into the phenomenon at play. This has been found to be a challenging task, demonstrated by the different results produced by various solvers [START_REF] Cocks | Impact of numerics on the predictive capabilities of reacting flow les[END_REF] . Present results remain nevertheless comparable to most reported attempts using LES methods [START_REF] Wu | Mvp-workshop contribution: Modeling of volvo bluff body flame experiment[END_REF][START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Ghani | Les of longitudinal and transverse self-excited combustion instabilities in a bluff-body stabilized turbulent premixed flame[END_REF] . Furthermore, Tab. II provides the computational cost for the cold and reacting flow using Reduced Computational Time (RCT) and Reduced Time to Solution (RTTS), computed as:

RCT = TCPU N iter × N p , RTTS = TCPU T sim × N p , ( 9
)
where TCPU is product of number of cores (N c ) and elapsed time to carry out simulation, N iter is number of iterations, N p represents number of points and T sim is the physical simulation time. Results are highly encouraging, with about 1000cpuh required to simulate one burner convective time on the coarse mesh.

In conclusion, the first lattice-Boltzmann large eddy simulation of a turbulent premixed flame stabilized in the wake of a bluff-body has been presented. Results are in line with those obtained with classical Navier-Stokes solvers. The reported CPU costs indicate that LBM is an interesting candidate for future LES of combustion processes. Supplementary material is available for (i) more detail on the pressure-based Lattice-Boltzmann model 10 ; (ii) the TFLES filtered equations [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] and (iii) the transport and thermo-chemical properties retained for the study [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF][START_REF] Rochette | A generic and selfadapting method for flame detection and thickening in the thickened flame model[END_REF] .

B. Rochette is gratefully acknowledged for providing the numerical data associated to his study [START_REF] Rochette | Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on les of premixed turbulent flames[END_REF] . The French Space agency (CNES) is acknowledged for supporting Song Zhao at M2P2. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

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

5 x

 5 FIG. 1. A Mid-plane cut view of the domain in z-direction is presented, indicating all the positions where measurements are available. The axes origin is identified as the red dot. The dimension of the equilateral triangle is D = 0.04m, spanning throughout the width of the burner. The dotted lines indicate the location of the refined region for Case B.

5 FIG. 2 .

 52 FIG. 2. Quantitative comparisons of our results for the cold flow, at two different grid definitions, with experimental data at five different points (x 1 , x 3 , x 5 ) (from left to right of each sub-figure). a) Mean U profiles. b) Mean V profiles. c) RMS U profiles. d) RMS V profiles. e) RMS UV profiles. Case A ( dashed line), Case B ( solid line) and experimental data ( solid square) f) Mean axial velocity profiles starting from the bluff body. Periodic boundaries on the sides (z-direction). Case A (solid line), Case B (dashed line), wall boundaries Case A (dotted line) and experimental data (black solid square).

FIG. 3 .

 3 FIG. 3. Vortex shedding pattern for the cold flow. Top: Q-criterion for Case A at 2.0 × 10 6 . Bottom: Q-criterion for Case B at 2.0 × 10 7

3 FIG. 4 .

 34 FIG. 4. Quantitative comparisons of our results for the reactive flow, at two different grid definitions, with experimental data at five different points (x 1 , x 3 , x 5 ) (from left to right of each sub-figure). a) Mean U profiles. b) Mean V profiles. c) RMS U profiles. d) RMS V profiles. e) RMS UV profiles. f) Mean axial velocity profiles starting from the flame holder. Case A (solid line), Case B (dashed line), experimental data (solid square), reference data -two/step chemistry 22 (dotted line) and reference data -detailed chemistry 22 (dash-dot line)

FIG. 5 .

 5 FIG.5. Iso-surface of the progress variable at 0.5 for the reactive test case (colored by wrinkling factor). Top: Case A. Bottom: Case B

1 )

 1 Aix Marseille Univ, CNRS, Centrale Marseille, M2P2, 13451, Marseille, France 2) CNES launchers directorate, 75612, Paris, France (Dated: 19 January 2021)

TABLE I .

 I VOLVO burner operating conditions[START_REF] Sjunnesson | Lda measurements of velocities and turbulence in a bluff body stabilized flame[END_REF][START_REF] Sjunnesson | Cars measurements and visualization of reacting flows in a bluff body stabilized flame[END_REF] selected for this study, where ν = 1.49 × 10 -5 m 2 /s.

	Case	φ	U b (m/s)	T 0 (K)	Re=U b D/ν
	Cold	0.65	16.6	288.0	45000.0
	Reactive	0.65	17.3	288.0	47000.0

TABLE II .

 II Computational costs for running a given time (≈ 7L x /U b ) of cold and hot ( * ) flows, on both meshes (A and B).

	Case N p	N c N iter T sim δt	TCPU RCT RTTS
	A	5.5M 512 0.6M 0.63s 1.06 µs 6 114 6.6 µs 58.2
	A *		1.06 µs 7 096 7.7 µs 67.5
	B	17.5M 512 1.2M 0.63s 0.52 µs 31 000 5.3 µs 147.6
	B *		0.52 µs 36 270 6.2 µs 172.7