Lattice Boltzmann methods for combustion applications

The lattice Boltzmann method, after close to thirty years of presence in computational fluid dynamics has turned into a versatile, efficient and quite popular numerical tool for fluid flow simulations. The lattice Boltzmann method owes its popularity in the past decade to its efficiency, low numerical dissipation and simplicity of its algorithm. Progress in recent years has opened the door for yet another very challenging area of application: Combustion simulations. Combustion is known to be a challenge for numerical tools due to, among many others, the large number of variables and scales both in time and space, leading to a stiff multi-scale problem. In the present work we present a comprehensive overview of models and strategies developed in the past years to model combustion with the lattice Boltzmann method and discuss some of the most recent applications, remaining challenges and prospects


Introduction
The lattice Boltzmann (LB) method, proposed in the early 80's has grown popular over the past decades [1,2].The rapid emergence of this numerical method is mainly due to the simplicity and strict locality of the involved time-evolution operators [3,4].The locality of the operators and intrinsic coupling between the pressure and velocity fields through the distribution function (as opposed to pressure-based incompressible or low Mach solvers) allows for better performances on parallel clusters and a much more efficient treatment of flows in complex geometries [4].During the past decade, the LB method originally proposed for computational fluid dynamics (CFD) has been extended to many complex flow configurations ranging from non-Newtonian [5][6][7][8][9], to multi-phase [10][11][12][13][14][15][16], and multi-component flows.Although initially limited to low-Mach isothermal flows with an ideal gas equation of state, the LB approach was later modified to lift many of these restrictions.Releasing the restriction on thermo-compressibility is an essential step to develop LB solvers for many applications such as combustion.
The topic of combustion modeling with LB was first touched upon in 1997 in an article by Succi et al. [17].Since then, and up until very recently, a limited number of publications had appeared on the topic, all limited to simplified 1-D and 2-D test-cases, see for instance [18][19][20][21][22][23][24].The limited progress of the lattice Boltzmann method during that period might be attributed to a number of factors such as the absence of a good compressible realization, persistent issues with stability of solvers, and the absence of multi-species formulations.During the past years a considerable amount of research work has been conducted to extend the lattice Boltzmann method to compressible flows, which has led to a number of stable and efficient realizations, see for instance [25][26][27][28].In parallel, the stability domain of lattice Boltzmann solvers both for incompressible and compressible flows has been considerably expanded through more advanced collision models, see for instance [29][30][31][32][33][34].These two factors along with the development of models for species transport and the idea of hybrid solvers taking advantage of classical numerical methods for the species end energy balance equations led to considerable progress in combustion simulation with the lattice Boltzmann method in recent years.Contrary to the first wave of models, the more recent efforts have been extended and used for many complex configurations involving thermo-acoustics, complex geometries and turbulent flows.
It has to be noted that in parallel with efforts to develop lattice Boltzmann-based models for combustion simulations, a number of attempts at developing discrete velocity Boltzmann-based models with Eulerian discretization in physical space have also been reported, see for instance [35][36][37][38].
In the present contribution we will review developments in the area of lattice Boltzmann simulations of combustion.Different challenges, solutions and models developed in that area in the past years will be presented and discussed.The review starts with a brief overview of basic concepts, i.e. target macroscopic system and basic concepts from the lattice Boltzmann method.In the third section of this review we will discuss topics specific to combustion simulations, i.e. strategies to solve the energy balance equation, models developed for species transport equations, and introduction of compressibility effects into the lattice Boltzmann solver.
The review closes with section four where key points are briefly listed and future prospects and challenges are discussed.

Brief overview of target macroscopic system
Throughout the manuscript, the target set of macroscopic equations is the multi-component system of Navier-Stokes-Fourier equations (see, e.g.[39]) Here u α is the α th component of the fluid velocity, ρ is the mixture density, E is the total energy (sum of internal energy e and kinetic energy u 2 α /2), Y k is the mass fraction of species k, and δ αβ is the Kronecker symbol (1 if α = β, 0 else).The above system is fully closed upon choosing Equation of state: a thermodynamic closure, linking state variables p, ρ, e, T and Y k , e.g.following the perfect gas assumption p = ρ.r.T = ρ R W T .
Transport models: to define the species diffusion velocities V k,β , heat flux term q β and viscous stress tensor τ αβ .
Chemistry model: to define the reaction rates ωk .

Isothermal lattice Boltzmann for incompressible flows
The construction of a discrete kinetic solver like the lattice Boltzmann method has two main ingredients: (a) Reduction of the particles' speed continuous space to a discrete set, and (b) discretization of the resulting system of hyperbolic equations in physical space and time.In this section these two components will be briefly reviewed.

Discrete velocity system and discrete equilibrium state
The rationale behind the construction of the lattice Boltzmann method consists in using a truncated version of the Boltzmann equation with a linear approximation to the collision term to recover the dynamics of the macroscopic equations of interest, here the isothermal Navier-Stokes and continuity equations.
From Boltzmann-BGK to the discrete velocity Boltzmann equations.Consistent with the terminology of the early literature, in the context of the present work we will refer to all methods using a form of the Boltzmann equation with a discrete set of particles' velocities as discrete-velocity models (DVM).In recent years interest in such models has been revived in the form of numerical methods such as the lattice Boltzmann equation.
DVM generally aim at approximating the distribution function with quadrature rules or similar integral approximations and using a discrete set of velocities: changing the Boltzmann-Bhatnagar-Gross-Krook (BGK) equation [40] into a set of coupled hyperbolic partial differential equations: Constraints on the discrete equilibrium function, f eq i , e.g.moments of the equilibrium distribution function to be correctly recovered, are identified via a Chapman-Enskog (CE) multi-scale expansion.For the continuity equation: ∂Π eq 0 ∂t + ∂Π eq α ∂x α = 0, (7) while for the momentum balance equations: where we have made use of the following notation: meaning that for the system of interest one needs to correctly recover moments of orders zero through three of the equilibrium distribution function.
In the specific context of the lattice Boltzmann method, the Gauss-Hermite quadrature is utilized to satisfy most of the above-listed conditions on the discrete distribution functions, i.e.
where P M (v, ρ, u) is a polynomial of order M of v and w(v) is a function of the form: For the quadrature to be applicable the distribution function must be expanded as the product of a polynomial series and w(v) and the integration variable, i.e. v normalized via the reference temperature, i.e. rT 0 .A change of variable as v ′ = (v − u)/ √ rT like Grad's expansion is not possible here as it would lead to changing discrete particle velocities that are also not necessarily space-filling.
Choosing the abscissae, i.e. c i to be the roots of the Hermite polynomial of order Q and the weights as: results in the maximum algebraic degree of precision, i.e. 2Q−1.This means that the quadrature guarantees exact recovery of moments up to order 2Q−1 2 .In the case of the classical lattice Boltzmann stencil, a third-order quadrature is used, i.e.Q = 3 in 1-D, with: and The simplest multi-dimensional extension of this quadrature can be obtained by operating a tensorial product of the 1-D lattice.For instance, in 2-D as illustrated in Fig. 1, this leads to the D2Q9 lattice with: and and A similar procedure is used to obtain the D3Q27 lattice for 3-D simulations.
Discrete equilibrium: polynomial form.A number of different ways for constructing the discrete equilibrium have been proposed over the years.One of the early approaches, first discussed in [41] was to re-write the equilibrium as: and Taylor-expand the last term around Mach Ma= 0, i.e.
ultimately leading to, after discretization of particles velocity space and application of the Gauss-Hermite quadrature, the second-order polynomial discrete equilibrium: An alternative construction based on an expansion of the distribution function with Hermite polynomial was proposed, which led to the following final form: where H n and a eq n are tensors of rank n representing respectively the order n Hermite polynomial and coefficient.
Alternatively the polynomial equilibrium can also be constructed via the product form.The product form of the equilibrium distribution function (EDF) is a special realization of the moments matching approach.
Considering the standard discrete velocity set D3Q27, where D=3 stands for three dimensions and Q=27 is the number of discrete velocities, one first defines a triplet of functions in two variables, ξ α and ζ αα , and considers a product-form associated with the discrete velocities c i (22), All pertinent populations below are determined by specifying the parameters ξ α and ζ αα in the product-form (26).The two-dimensional version of the model on the D2Q9 lattice is obtained by omitting the z-component in all formulas.After matching moments with their continuous counter-parts the parameters are set as, and the local equilibrium populations are represented with the product-form (26), This form of the discrete equilibrium populations, when c 2 s = r0 T 0 /3 is equivalent to third-order quadraturebased scheme with a full expansion of the distribution function.
Alternative to polynomial equilibria: Entropic equilibria.As an alternative to the classical discrete equilibrium construction approach where all degrees of freedom are used to fulfill moments constraints, the entropic approach adds minimization of an entropy functional to the list of constraints, changing the equi-librium construction problem into a constraint minimization problem.While a number of different discrete entropy functions have been proposed in the literature the most commonly used one is: Minimization of this functional under constraints on moments of order zero and one, leads to the following well-known entropic equilibrium: One of the most interesting feature of this equilibrium, contrary to polynomial equilibria is that, as demonstrated in [42,43], it guarantees unconditional linear stability.

Lattice Boltzmann equations
From discrete velocity Boltzmann to the lattice Boltzmann method.To go to the final form of the lattice Boltzmann equations, two main ingredients are to be used: (a) integration of the discrete velocity Boltzmann equation along their constant characteristics and (b) a re-definition of the discrete distribution functions.
The former step results in: where the term on the right-hand side, representing collision, has to be approximated via the trapezoidal rule, in order to keep the scheme second-order accurate: However, as observed here, application of the trapezoidal rule would make the scheme implicit and therefore not attractive regarding efficiency.The second ingredient, i.e. redefinition of the discrete distribution changes the system of equations into a fully explicit scheme: where τ is now defined as: The lattice Boltzmann method for incompressible flows.A multi-scale analysis of the above-listed system of equations shows that it recovers the isothermal continuity plus Navier-Stokes system of equations for an ideal equation of state at reference temperature: where the coefficient 1/3 is specific to the third-order quadrature.Note that the recovered system of macroscopic equations further admits a defect in the viscous stress tensor; For the classical second-order polynomial expansion defects are present both in shear and bulk viscosity scaling in both cases as ∝ U 2 α δt 2 /δx 2 , where U is the characteristic velocity along the α-axis.For the full polynomials expansion or product-form equilibria, this defect is only present in the bulk viscosity.This means that under acoustic scaling, i.e. δx/δt = const the solver converges to the compressible isothermal Navier-Stokes equations, as for a fixed characteristic velocity U the Mach number remains constant in the limit of δt → 0. Furthermore, under acoustic scaling the defects in the effective viscosities do not vanish.
Under diffusive scaling on the other hand, i.e. δx 2 /δt = constant, the solver converges to the incompressible Navier-Stokes equations, as U/c s → 0 in the limit of δt → 0 and the defect in the effective viscosities also goes to zero.

Lattice Boltzmann models for compressible reacting flows
We have now introduced the main ingredients of classical Lattice-Boltzmann methods, and shown that they allow to recover the continuity (1) and momentum (2) equations of a weakly compressible gas at constant temperature T 0 .This is indeed not sufficient for combustion applications, where energy (3) and species (4) are absolutely required.This Section is divided into 3 main subsections.First, we shall explore different alternatives for the resolution of the additional equations (energy and species) in Section 3.1.Second, we will detail the required changes to the lattice Boltzmann formulation in Section 3.2.Finally, we will list the reactive flow configurations successfully simulated by LBM solvers and discuss their performance in Section 3.3.

Energy and species balance equations
3.1.1.Double distribution function lattice Boltzmann approach for thermal flows Kinetic models.Historically, the starting point of double distribution function (DDF) approaches is rooted in the need for simulations with variable Prandtl numbers and specific heat capacities, as alternatives to Holway's ellipsoidal statistics [44] or Shakhov's model [45].This is usually achieved by introducing a second distribution function g to carry a form of energy following [46], which is not uniquely defined.The earliest occurrence of a double distribution function approach is documented in [47] where the authors introduced: Multiplying the Boltzmann equation, i.e. the balance law for f by the coefficients in the definition of the g-distribution function one obtains the balance law for the latter: where the additional non-homogeneous contribution q is: In this model, the total energy E is computed as: Some comments on this approach are necessary: • This model, through the choice of parameter τ g allows for a variable Prandtl number Pr.
• The model assumes a mono-atomic molecule as no degrees of freedom in addition to translational are taken into account.
• The model involves space and time derivatives of macroscopic fields.
To alleviate the last issue, Guo et al. proposed to carry total energy with g instead [48]: While this choice of a second distribution function leads to a much simpler balance law for g it also comes with a limitation of the Prandtl number.Contrary to the previous choice of g carrying internal energy where one could easily vary Pr by changing τ g , here the relaxation time in the collision operator controls both relaxation of internal and kinetic energy, therefore also affecting viscous heating.To allow for variable Pr, the authors proposed to decompose the collision term into kinetic and internal contributions, leading to the following balance law: where and Since g carries the total energy it is computed solely as its zeroth-order moment: In the same contribution the authors proposed a more generalized framework allowing to incorporate additional non-translational degrees of freedom into the model by defining g as: where η is a vector with δ component, with δ the number of additional degrees of freedom, and summation over both α and β is assumed.In this model the equilibrium distribution function is: and the total energy is computed as: While Guo et al. originally proposed this decoupling for the low Mach limit, it was extended to compressible regimes in [49] where the authors used a thirteen-velocity lattice for the f distribution function to eliminate the deviations in the third-order moments of the equilibrium distribution function.A realization of this model on the standard third-order quadrature-based lattice was proposed in [50].The approach originally proposed in [48] has been routinely used since then in a wide number of publications for both compressible and incompressible flows, see for instance [25,[51][52][53].Another realization of the double distribution function method relying on internal energy was also proposed in [54].As noted by the authors, re-writing the balance equation of Eq. ( 43) as: in the case of the model in [48]: while for [54]: Note that both realizations lead to the same hydrodynamic equation and in the case of the third-order quadrature-based lattices, even to the same discrete equations [54].This realization has also been used for a variety of compressible flow simulations, see for instance [27,55].
Lattice Boltzmann equations.Discretization in the space of particles velocities v proceeds very similarly to that of the probability distribution function f , either through projection onto the space of Hermite polynomials or via the product form construction.In the product form approach discussed in [27,54] similar to Eq. ( 29): where the operator O α acts on any smooth function A(rT, u α ) as: The discrete-in-particles speed-space can then be integrated along characteristics to obtain the corresponding collision-streaming equations: where the new distribution function is: and The total energy can then be obtained by summing discrete distribution functions: A multi-scale analysis shows that the models above, at the Euler level, recover: where H is the enthalpy.At the Navier-Stokes level: i ) Π α (g where the first term is the Fourier diffusive flux while the second term is viscous heating.The double distribution function approach in combination with a proper lattice Boltzmann solver for density and momentum balance, detailed in next sections, has been used to model trans-and supersonic flows, for instance in [27].
A few results are illustrated in Fig. 2.
Multi-species flows.For the extension of this model to the case of multi-species flows in the context of a mixture-averaged formulation a few points must be noted: • In all models discussed here, at the ϵ 2 order a Fourier flux of this form is recovered, see Eq. ( 62): which holds if enthalpy is only function of temperature.For a mixture-averaged formulation with multiple species, H = H(T, Y k ), which would lead to: where H k is the enthalpy of the k th species.
• In multi-species flows, diffusive mass flux leads to a net transport of enthalpy which is absent in single-component flows.
A solution to both these shortcomings was proposed in [56] to recover consistent hydrodynamics, where the pseudo-equilibrium g * i is amended with two additional terms, one neutralizing the error in the Fourier diffusion and one introducing enthalpy flux through mass diffusion:

Kinetic models for species balance equations
Over the past decades and starting in the early 2000's [57,58] various attempt at developing lattice Boltzmann-based models for mixtures have been documented, see for instance [59][60][61][62][63].Some of these models are reviewed in this section.
Thermal mixture-averaged model of Kang et al.In [64], the authors proposed a multi-component thermal model for catalytic systems.The model is an extension on previous work documented in [65][66][67][68].It consists of N sp sets of lattice Boltzmann solvers, i.e. one per species: The first point to note in this model is that post-streaming discrete distribution functions migrate to x+c ki δt meaning each species' lattice have different discrete velocity sizes.As discussed in previous sections, in lattice Boltzmann time-step, grid-size and reference temperature are tied through: which in the context of this model where W = W k is different for each species, and assuming that the time-step size is the same for all solvers, would mean: i.e. not all species will propagate on lattice.To overcome this issue, and following [69] the authors proposed to set the grid-size to that needed for the lightest species in the system, and for other species to use interpolation in order to reconstruct distribution functions on the grid.The equilibrium, g eq ki (ρ k , u) follows the product-form equilibrium of Eq. ( 29) with a few differences, namely: while for the pseudo-equilibrium g * ki (ρ k , u k ): In this model the second relaxation time, τk2 sets the diffusivity to: where D k is the mixture-average diffusion coefficient.The viscosity is set through the first relaxation time and using Wilke's formula: with The term ψ ki in Eq. ( 66) is a correction term accounting for: (a) a correction velocity ensuring that the global diffusive mass flux is null, and (b) corrections for equilibrium moments of order three and four not recovered by the first-neighbour lattice.The latter terms allow the scheme to recover the proper viscous stress tensor and non-equilibrium heat flux.In this model the macroscopic properties are computed as: A multi-scale expansion of the model shows that it recovers the mixture-averaged multi-species equations and the Hirschfelder-Curtiss approximation with the mass corrector.
Force-based approach of Vienne et al.In [70], following the kinetic model of [71], Vienne et al. proposed a lattice Boltzmann model for isothermal multi-species mixtures recovering the Maxwell-Stefan system of equations.Considering a mixture made up of N sp individual species, they proposed a coupled system of N sp lattice Boltzmann solvers: where S i is here to introduce external body forces, realized using Guo's approach [72]: where F k represents the body force.In this model and The interaction between different species driving diffusion is introduced via a body force defined as: where D kk ′ represents the binary diffusion coefficients.As noted by the authors, the circular inter-dependence between the force of Eq. ( 83) and the momenta of individual species of Eq. ( 82) make the scheme implicit.
A multi-scale analysis shows that this model recovers the following multi-component isothermal mass Figure 3: Illustration of application of multi-species model of [70].Evolution of viscous fingering instability for a system with two species.Image reproduced from [73].
and momentum balance equations where the bulk viscosities µ k are defined as: The model has been successfully used to study miscible multi-component flow behaviors such as viscous fingering, see Fig. 3.An extension of this model to thermal and reacting cases is yet to done.

Model of Sawant et al.
In [56,74], the authors proposed a kinetic model to recover the Stefan-Maxwell diffusion model.Each component is described by a set of populations g ki .The discrete-velocity time evolution equation is, The species densities are computed as zeroth-order moment of the discrete distribution functions: The symmetric set of relaxation times θ kk ′ = θ k ′ k is related to the binary diffusion coefficients.The firstorder moments of the distribution functions are, The quasi-equilibrium populations g * ki satisfy the following constraints on moments, The momenta of the individual species sum up to the mixture momentum, The equilibrium populations g eq ki are subject to the following constraints: In Eq. ( 95), the partial pressure p k depends on the mixture temperature T which is obtained from the energy balance lattice Boltzmann solver.Noting that and using the equation of state the kinetic model can be re-written as: This equation, for the sake of convenience, is recast in the form of a relaxation equation: where and

Passive-scalar lattice Boltzmann models
So-called passive-scalar lattice Boltzmann solvers, are models where only conservation of the zerothorder moment of the distribution function is ensured by the collision operator.In such models, to solve an advection-diffusion-reaction partial differential equation for a field Ψ, a distribution function g i is defined such that: where S is the source term.A classical non-homogeneous collision-streaming equation of the form: with and leads to a macroscopic equation of the form: Note that in the literature, Eq. ( 104) has been used both with first and second-order polynomials expansions.
Depending on the choice of the order of expansion the diffusion term will admit errors of different forms.
For instance, for a linear equilibrium, On that note, let us now discuss the passive scalar approach in the specific context of the species mass and energy balance equations.Although a number of different works have presented modified passive scalar approaches for non-linear dependence of the diffusion driving force on the zeroth-order moment even in the context of multi-species flows (see for instance [76,77]), here we will limit our discussion to models that have been used for combustion simulation.
Energy balance equation.The energy balance equation can be written in a variety of different ways, see [39].
Here we will only discuss the form involving temperature: The classical approach to recover this balance equation is to set Ψ = T which would lead to the following shortcomings: • An error of the form T ∂u α /∂x α in the convection term.
• An error of the form: in the Fourier diffusion term.
• The enthalpy flux due to species mass diffusion is missing.
While one can overcome these issues via alternative forms of the equilibrium distribution function (see for instance [78]), the simplest way to circumvent these issues is to define the source term S in Eq. ( 105) as: Such an approach, among others, has been used in [79,80].Note that in [79,80] the last term was not considered.A similar approach can be undertaken for the case were the enthalpy or energy balance equation is targeted.
Species mass balance equations.For the sake of simplicity let us assume that the species mass fraction balance equation is targeted.Taking the zeroth-order moment of the distribution function to be mass fraction, Y k and neglecting for the time being leading-order errors in the multi-scale analysis one would recover a diffusion term of the form: which using τ /δt = D k /c 2 s + 1/2 recovers the generalized Fick approximation.With that the passive scalar approach is confronted with a number of issues: • This form of diffusion is only valid in the limit of vanishing density changes as in the non-conservative form of the balance equation there is a factor 1/ρ in front of the diffusion term.
• The form of the convection term as recovered in Eq. ( 105) admits an error of the form Y k ∂u α /∂x α , which only vanishes for incompressible flows.
• It is well known that the generalized Fick approximation does not conserve overall mass, unless either To deal with that issue there are two approaches; If mass fraction of one particular species is dominant everywhere (e.g. that of N 2 in combustion with air) the balance equation for that species is not explicitly solved and one sets Y N2 = Nsp−1 k=0,k̸ =N2 .A more general approach, valid also for non-dilute mixture is to introduce a mass corrector.
• In cases where the driving force of the diffusive flux is a linear function of the variable for which the balance equation is solved, for instance the Fick approximation, the passive scalar model can be used as is.However, for models where the driving force of diffusive flux is a non-linear function that depends on variables other than the zeroth-order moment, for instance the Hirschfelder-Curtiss approximation, the passive scalar approach would lead to errors of the same order as the diffusive flux itself.
A number of different approaches have been proposed in the literature to account for these shortcomings, see for instance [76,77].One of the most straightforward approaches, as used in [79,80], is to put all corrections into a source term.For instance, assuming one targets the mass fraction balance equation with the generalized Fick approximation, the source term S would be: Note that this approach as used in [79,80] still comes short with respect to the mass corrector and the more appropriate diffusion velocity closure.A number of solutions to account for the mass corrector have been proposed in the literature, taking advantage of the non-equilibrium part of the distribution function, see for instance [76].Let us now list the main advantages of the hybrid method: 1.The memory footprint is reduced, as only 1 additional scalar needs to be stored for each energy or species equation, vs. 27 for a D3Q27 distribution.
2. They are by construction free of Prandtl, Schmidt or γ numbers limitations, since energy/species resolution is tackled separately.
3. Since they use the same formalism as classical reactive flow solvers for energy and species equations, it is straightforward to take into account combustion-specific terms (turbulent combustion closures, advanced transport models, Soret effect or multi-component diffusion etc.), based on the experience accumulated over many decades using FD/FV solvers.
In turn, hybrid methods suffer from the following drawbacks: 1. Ensuring consistency between the LBM scheme (for continuity and momentum equations) and FD/FV schemes is not straightforward (at the opposite of, e.g.multi-speed or multiple distribution approaches).This can typically lead to disastrous spurious currents, as illustrated later in Fig. 11.
2. FD/FV schemes based only on nearest-neighbor stencils (as used in most LBM solvers) are typically much more dissipative than LBM schemes [81].
The first point is crucial in designing hybrid LBM schemes, and is therefore discussed at length hereafter.
The impact of the second point is limited for most applications as long as the vortical and acoustic modes are left within the LBM part of the solver.
Which form to use for species/energy equations in a hybrid LBM scheme?.Energy and species equations may be written under a large variety of forms (based on total energy, internal energy, temperature,...).While these forms are indeed equivalent for a continuous formulation, their coupling under discrete form with the LBM scheme may be very different.Let us recall that for small perturbations and neglecting all dissipation terms (reducing to the multicomponent Euler equations), the system (1-4) may be linearized, and each perturbation can be decomposed into the so-called Kovasznay modes [82,83] (acoustic mode, 3 components of the vorticity modes, entropy mode, and 1 per species).
For instance, the entropy mode of the Euler system follows the equation and is only weakly coupled with the rest of the system.
For this reason, hybrid methods using an entropy equation were shown to provide reasonable results for moderately compressible flows [84-86] using several classical convective numerical schemes a priori unrelated with LBM: • Second-order central difference schemes, potentially blended with upwind, [87][88][89][90] • Lax-Wendroff scheme [91,92] • MUSCL schemes [84][85][86] • Heun scheme [93], • . . .Indeed, for reactive flows, the entropy equation is complex to derive in its general case.However, the enthalpy equation under non-conservative form is also a characteristic mode of the system -provided the pressure work dP dt is neglected, a very common assumption for low-Mach reactive flows.
Species also directly follow characteristic equations, provided they are written under non-conservative form There is an alternative but not equivalent way of understanding how crucial this choice is.Consider the species equation in conservative form It is clear that this equation is the sum of the non-conservative form (a system characteristic) and the continuity equation.Therefore, any numerical error between the continuity equation solved by LBM and the one hidden into the conservative form leads to an inconsistency.
To summarize, provided the equations to be solved using FD/FV are only weakly coupled with the rest of the system, the resulting hybrid LBM solver has been shown to provide reasonable results for a wide number of cases.
Restoring the conservativity of hybrid LBM.Equations (110,111) are equivalent to the initial total energy and species equations (3,4), but the discrete formulation is not.This has two disadvantages: • Global energy conservation is not numerically enforced (while the LBM scheme is numerically mass and momentum preserving).
• Rankine-Hugoniot relationships are not satisfied across discontinuities.
The latter is clearly visible in Fig. 5, which presents a reference 2-D Riemann problem and the solution as obtained with hybrid LBM using Muscl-Hancock scheme to solve the entropy equation.Wissocq et al. [92] recently presented a method to construct a linearly equivalent total energy scheme from any FD/FV scheme that is linearly stable for hybrid LBM (e.g., for the entropy equation).

Compressible continuity and momentum balance equations
All strategies presented above for the resolution of the additional energy and species equations coupled to a LBM solver require modifications to the LBM core.The major options are presented hereafter.

Lattices with higher-order quadratures
Standard approach with polynomial equilibria.As discussed in the isothermal models sections the Maxwell-Boltzmann phase-space continuous equilibrium can be matched at the discrete level via a number of methods following the same general principle, matching the different moments of the continuous equilibrium with the discretized version.As such the classical moment-matching method routinely used for Eulerian discrete velocity Boltzmann models and truncated Hermite expansion approach both fall into that category.In the case of the former one, once the number of constraints on moments of the equilibrium and degrees of freedom in the form of number of discrete velocities have been set, construction of the discrete equilibria boils down to solving the following linear system: where Π M B is a vector of size 1 × Q, Q being the number of constraints on moments, with moments of the Maxwell-Boltzmann continuous distribution corresponding to the targeted constraints with: with n = p + q + r.The quantity f eq is the vector of sizes 1 × Q containing discrete equilibria and M the transformation matrix from discrete equilibria to moments.For instance, in 1-D, for a solver targeting the Navier-Stokes-Fourier dynamics a minimum of five discrete velocities are needed as one must correctly recover moments of order zero to four.This approach to construct a discrete solver, while being quite flexible has a number of shortcomings, namely: The matrix M is not necessarily invertible for any choice of moments and discrete velocities as illustrated by the introduction of so-called anti-symmetric discrete velocities in some higher-order discrete velocity Boltzmann models, see for instance [95]; while the number of velocities is set by the constraints the sizes and size-ratios of these velocities have no a priori closures and are usually tuned via trial and error.A possible closure for the size of the discrete velocities would be to use Hermite polynomials roots of the corresponding order.The only issue with that choice is that above order three Hermite roots do not guarantee space-filling lattices and therefore on-lattice propagation.
Nevertheless, a large number of publications using larger discrete velocity sets are documented in the literature: • A group of these publications do not rely on Lagrangian approaches to discretize physical space and time and use Eulerian approaches such as finite differences or finite volumes to discretize the coupled system of hyperbolic equations of the discrete velocity Boltzmann model.In doing so the sizes of discrete velocities can be freely tuned to stabilize the simulation for a specific test-case.
• Another group of publications use Lagrangian approach to discretize physical space and time and overcome the issue of non-space-filling lattices by supplementing the collision-propagation step with an interpolation step to bring back post-streaming discrete velocities on lattice.These approaches are sometimes referred to as semi-Lagrangian, see for instance [96,97].
• Another category of publications, relying on the classical on-lattice method, proposes to stabilize multi-speed lattice Boltzmann solvers for compressible flows through different collision models, such as multiple-relaxation time or regularized, see for instance [98,99].
All of the previously-listed models have had limited success in modeling generic high-speed compressible flows with large temperature variations in the domain.A number of alternatives have been proposed since then to considerably widen the stability domain of multi-speed lattice Boltzmann solvers.They will be discussed next.
Extension of stability domain: Entropic equilibria.The entropic construction of the discrete equilibrium state introduced for isothermal models, can be reformulated in a more general form as a minimization problem subject to M constraints: The formal solution of this constrained minimization leads to a function of the following form: Note that other form of the minimizer without the weights w i have also been proposed and used in the literature [100], most notably for entropic Grad moments methods [100,101].For instance, a model imposing only constraints on collisional invariants, i.e.
would lead to the following discrete equilibrium [102]: It is interesting to note that while, for the most part, entropic equilibria construction has been done by enforcing constraints on collisional invariants, one may reduce higher-order moments error by adding corresponding constraint in Eq. ( 115).This is sometimes referred to as guiding the equilibrium and corresponding discrete equilibria are referred to as guided equilibria [103,104].In the context of the lattice Boltzmann method, this extension of constraints was discussed for the first time in [105] through the concept of auxiliary and target equilibria.There, auxiliary equilibria were constructed by enforcing constraints on collisional invariants and target equilibria, a combination of auxiliary equilibria and additional degrees of freedom, by enforcing constraints on higher order moments.
Once the form of the equilibrium distribution function has been determined, its construction consists of finding the expression of the different Lagrangian multiplicators.This is done by introducing back the discrete equilibrium into the set of constraints which would lead to a system of M equations with M unknowns, i.e.
the Lagrange multipliers, to be determined.While an analytical expression was derived for the isothermal case with D + 1 constraints, for larger systems no such solutions exist.In the absence of a closed form solution one can use numerical methods such as Newton iterations to find the Lagrange multipliers at every grid-point and every time-step [106].
As shown in previous sections, one systematic approach to choose an optimal set of discrete velocities is to rely on the Gauss-Hermite quadrature and roots of Hermite polynomials.However, apart from the third-order quadrature leading to the DdQ3 d lattices, all other higher order quadratures result in off-lattice propagation of some of the discrete distribution functions.In [107], starting from a set of discrete velocities the authors proposed an approach to find a reference temperature and corresponding weights.This is achieved through the closure relation and matching conditions.For a set of discrete velocities V with Q vectors c i , the Q th power of c i can be written as a linear combination of lower order odd-powers from Q − 2 to 1, i.e.
For instance, in the case of the D1Q3 lattice one has c 3 i = c i .This essentially means that the moment of order Q is not an independent moment and can not be set at one's will.The only possibility is to set the linear in u term of the Q th order to its Maxwell-Boltzmann counter-part and in doing so determine the reference temperature, which is referred to as the matching condition.Consider for instance the D1Q3 lattice again.The third-order moment is going to be u x while the Maxwell-Boltzmann distribution leads to u 3 x + 3T 0 u x .To match the linear term one must have 3T 0 = 1.Note that not any choice of lattice admits a reference temperature.For example the velocity set V = {−2, −1, 0, +1, +2} will lead to a closure relation of the form c 5 i = 5c 3 i − 4c i and a matching condition, 15T 2 0 − 15T 0 + 4 = 0, which does not admit any solutions.This explains why the shortest admissible five-velocity lattice is V = {−3, −1, 0, +1, +3} with T 0 = 1 ± 2/5.Once the reference temperature is determined, the weights are readily found by matching the moments of the discrete equilibrium at ρ = 1 and u x = 0 to their Maxwell-Boltzmann counter-parts.
Considering the condition of positivity of the weights one also find the range of temperature that can be covered by the chosen system of discrete velocities.The closure relations and reference temperatures of a number of 1-D lattice are summarized in Table 1.
One successful example of such lattices is the D2Q49 shown in Fig. 6.The closure relation for the 1-D set is  The 1-D weights read: which lead to T min = 1 − 2/5 and T max = 1 + 2/5.Note that the range of accessible temperatures can be further extended by changing the ratio of the largest and shortest discrete velocities, here ±3 and ±1.In [106] the author also proposed pruning strategies to reduce the number of discrete velocities in 2-D and 3-D, leading to the D3Q39 lattice, which reduces the discrete velocities by one order of magnitude compared to the tensor product of the D1Q7, i.e.D3Q343.
Adaptive reference frame models.As observed for both isothermal and compressible models, errors in higherorder moments scale with the deviations of local temperature and velocity from the lattice reference temperature and velocity.For all symmetric lattices considered up to that point the lattice reference velocity is U = 0.In [108] the authors proposed to challenge the idea of a reference frame at rest by introducing a non-zero shift U .It was noted that the discrete entropy functional is uniquely defined by the weights w i .The weights of a lattice with Q discrete velocities, as shown in the previous section, are determined by  matching the first Q moments of the Maxwell-Boltzmann equilibrium distribution function at temperature T and u x = 0:  demand (PonD) method [28].In this approach the collision streaming operation is performed on a reference frame corresponding to the local velocity and temperature.This allows to minimize higher-order moments deviations of the discrete equilibrium from the Maxwell-Boltzmann distribution function and in doing so allows for arbitrarily large variations in speed and temperature.The particles on demand method has been used to model high Mach number cases in recent years [110].It is also currently used in combination with a Lee-Tarver reaction model to simulate detonation at high Mach numbers, see [111].

Standard lattice density-based solvers
Coupling to temperature field.As discussed in the first section, the original lattice Boltzmann method was targeting the incompressible limit as the asymptotic behavior of the incompressible Navier-Stokes equations.
To that end the temperature appearing in the equilibrium distribution function was that of a reference state guaranteeing validity of the low Mach assumption.The compressible Navier-Stokes equations can be recovered by replacing the reference temperature with the local fluid temperature obtained from the second distribution function or the FD/FV solver used for energy balance; Considering for instance Eq. 28 it changes into: where θ = rT /r 0 T 0 .Introducing this term allows for a correct recovery of Euler level pressure while setting the relaxation time to: allows for correct recovery of the Navier-Stokes level viscous stress coefficient.With the temperature now an independent space-and time-varying parameter, it will inevitably deviate from the reference temperature, i.e. θ = 1 which is the optimal operation temperature of the third-order quadrature-based lattice Boltzmann model.Deviation from the reference temperature comes with a number of difficulties.The first one is the reduced domain of stability, illustrated best by the linear stability domain shown in Fig. 10.The second is that to properly recover the full Navier-Stokes viscous stress tensor a number of additional considerations have to be taken into account.These are discussed in the next paragraph.
Galilean-invariance of third-order moments and corrections.As discussed for the isothermal lattice Boltzmann method in previous sections, a simple CE analysis shows that at the NS level, moments of orders two and three of the EDF must be correctly recovered.Diagonal components of the third-order moments tensor, i.e. moments of the form Π ααα , can not be correctly recovered due to the c 3 iα = c iα bias of the third-order quadrature-based lattice.While the continuous Maxwell-Boltzmann equilibrium distribution leads to: any of the discrete equilibrium distribution functions discussed here recovers: which for θ = 1 introduces a cubic-in-velocity error and for θ ̸ = 1 a linear one.As such the issue of Galileanvariance of the third-order moments becomes quite critical in the case of compressible flows where θ ̸ = const.
To account for this error, corrections in the form of source terms in the kinetic equation are introduced: The form of the source term, Ψ i , can be derived through the order-two-in-ϵ (NS level) momentum balance equation: leading to: where For the stress tensor to be correctly recovered at this scale one must have: Note that to get this expression the correction term was assumed to involve first-order derivatives via the i .A different form of the correction term can be obtained with a different expansion, i.e.
i .Such an expansion would lead to the following NS-level equation: and a correction term of the form: The above-listed corrections were derived for the discrete kinetic equations.The classical lattice Boltzmann approach to space/time discretization would lead to the following redefined discrete distribution function: which in turn would lead to the following final algebraic system: This consistent derivation of the extended LBM holds for any realization of the correction term, whether it is introduced simply as a Hermite-expanded term [52] or by using the extended equilibrium approach [27].

Pressure-based solvers
While the density-based model detailed in the previous section was successfully used for a number of applications, see [88,113], it was observed that it led to spurious currents near curved flame interfaces, see Fig. 11 for a circular flame.A detailed study of numerical properties of the model showed that this is a result of a non-physical coupling between entropic and vorticity mode, see [114].A pressure-based formulation was then proposed [84] to cure the problem, the detailed reason only being understood later [85].The pressure-based algorithm is presented hereafter.
Note the differences compared to the density-based model presented in the previous section: • The zeroth moment of f i is now p = ρθ instead of ρ.
• The macroscopic reconstruction of Step 4 now includes a correction ρ(t + δt, x) = This second point was presented by the authors as a predictor-corrector procedure, close to early artificial compressibility methods [119,120].It is important to note, however, that despite being pressure-based, the algorithm is mass conserving globally, as density-based methods.
Link between pressure-based and density based formulations.Since many mesh transition algorithms, boundary conditions, etc. were initially developed for density-based algorithms [121], there is an interest in establishing a rigorous link between pressure-based and density-based algorithm.
This can be obtained noting that the correction ρ(t, x)[1 − θ(t, x)] in the macroscopic reconstruction can be equivalently embedded directly in the f 0 term corresponding to the stationary discrete velocity by introducing the density-based function: where δ 0i is the Kronecker symbol.By projecting δ 0i on a Hermite polynomial basis, it was further shown [85] that this change is equivalent to adding to the classical f ρ i a fourth order contribution, leading to an equilibrium function of the generic form with additional information projected onto fourth order polynomials A i , B i and C i .In the model, κ and ζ are free parameters.For instance, (κ, ζ) = (1, 1 − θ) corresponds to the density-based model of [117], while (κ, ζ) = (0, 0) yields the pressure-based model of (136).
Successful applications of the resulting generic model include the modelling of • Hele-Shaw cell [87]; • Turbulent premixed combustion burners; [90,122] • Turbulent lifted H 2 -air jet flame [123]; • Thermo-acoustic instabilities [122,124]; • Cellular detonation structure [93]; with the last two points being illustrated in Fig. 12.  and more generally of dilatable flows [125].A low Mach reduction of the fully compressible models of the previous section was also proposed in [126].The scheme is categorized as low Mach in the sense that it follows the philosophy of Majda's zero-Mach model [127] where after a Helmholtz decomposition of the velocity field, the divergence-free part is obtained via Poisson's equation and the curl-free part from the species and energy fluxes.At the difference of Majda's zero-Mach model, here the divergence-free component solver allows for a certain level of compressibility, i.e. spurious acoustic waves.To recover this modified set of macroscopic equations the model makes use of the following modified kinetic system [125]: where and the source term Ξ i is defined as: In this model, hydrodynamic pressure p h and velocity u are coupled via the distribution function via: while density ρ is now a variable computed locally through the ideal equation of state: The source of divergence appearing in Eq. ( 146) is computed via the continuity equation combined with the energy and species balance equations: where summation of α is assumed.A multi-scale analysis of this kinetic model shows that the following balance equation is effectively applied to the hydrodynamic pressure [34]: while for momentum, as for the classical lattice Boltzmann with second-order equilibrium, the Navier-Stokes equation is recovered with a deviation of order ∝ ∥u∥ 3 δt 3 /δx 3 in both diagonal and deviatoric components of the viscous stress tensor.
Note that after integration along characteristics, the discrete time evolution equations for the re-defined distribution function ḡ′ i are obtained as: while moments are computed as: While only the most basic single relaxation time form of this model is introduced here, interested readers can readily get access to more advanced realizations, for instance via the Cumulants-based multiple relaxation time collision operator in [32,128].
Over the past couple of years this low Mach lattice Boltzmann solver, in combination with shock-capturing finite differences schemes like the weighted essentially non-oscillatory (WENO) have been used to model a variety of complex reacting flow configurations, including • Turbulent combustion [129]; • Combustion in swirl burner [128]; • Combustion in porous media [130].
Some of these applications are illustrated in Fig. 13.A simpler form of this model was also used in [131] to  model droplet combustion.In terms of limitations, the model is -as obvious by its name -limited to the low Mach regime.Furthermore, as mentioned previously, the Navier-Stokes level viscous stress tensor admits deviation, for which corrections are to be published in an upcoming article.Finally, exact conservation is a topic that needs improvement here as the form of the energy and species balance equations along with the finite-difference discretizations and curved boundary treatment all lead to loss of exact conservation.

Performance of multi-physics LBM solvers
Let us now provide details regarding the advantages and limitations of multi-physics LBM solvers, and how they compare to classical LBM solvers targeting mainly low-Mach aerodynamic and aero-acoustic applications.
Classical LBM solvers have found their success due to three main reasons : (i) ability to tackle complex geometries in a simple way, (ii) low dissipation owing to the streaming step, and (iii) a reduced computational cost due to the stencil compactness and octree structure.These advantages come at the cost of a heavier memory load (with more degrees of freedom), a complex treatment of boundary conditions (owing to the non body-fitted mesh), and filtering problems in grid refinement areas where both spatial and time discretizations are halved/doubled.
Let us now consider the three aforementioned advantages and see if they apply to LBM multiphysics solvers.(i) Ability to tackle complex geometries is preserved, as the discretization remains identical.(ii) The low dissipation property is more model dependent.For multiple distribution functions, each distribution is solved using the same algorithm so that dissipation is supposed to remain identical.For hybrid approaches, however, a separate numerical scheme is used for energy and species equations, which may lead to additional dissipation [84] on the entropy/enthalpy Kovasznay modes.Comparing the computational cost (iii) between LBM and Navier-Stokes solvers is a long standing issue.While a thorough comparison is still lacking, some LBM studies report reduced computational times (RCT), consistently below 5µs per time step and grid point for relevant combustion problems [89,90,122,132].
We list in Tab. 2 different references tackling combustion problems using multi-physics LB solvers (regardless of the strategy).This fast-expanding list is a clear indication that multi-physics LBM solvers have now reached sufficient maturity for combustion applications.

Conclusion and discussion
While extension of the lattice Boltzmann method to the simulation of combustion happened slower than other areas of application such as incompressible and multi-phase flows, the progress documented in recent years has laid the ground for widespread use and application of lattice Boltzmann to combustion.The complex simulations reported in the literature, see for instance [93,122,128,130], show that the numerical approach has reached a level of maturity allowing it to be applied to realistic configurations.In this contribution we focused on the development of lattice Boltzmann-based approaches to model combustion and discussed some of the most pressing challenges and solutions proposed in the literature.
The evolution of the literature shows that one of the major challenges preventing progress in that area was the development of stable and efficient solutions to extend the lattice Boltzmann solvers to compressible regimes.Stability has been one of the most restricting issues in that area.While use of higher order lattices is one straightforward approach to move up to compressible flows, it has been observed that apart from the additional memory and computation costs stemming from the larger number of discrete velocities, higherorder quadrature are subject to more restrictive stability conditions, especially on the temperature.This has led the community to opt for approaches relying on low-order quadratures, i.e. third-order, which had shortcomings regarding the Navier-Stokes-level viscous stress tensor due to insufficient degrees of freedom in the model.Introduction of correction terms for the viscous stress tensor along with more robust collision operators has now paved the way for simulations involving large temperature variations.
Other issues specific to the simulation of combustion in the context of the lattice Boltzmann method are tied to additional balance equations for species and energy.A number of different strategies have been devised for these additional fields.Some rely on developing kinetic models and, therefore, lattice Boltzmann solvers for multi-species fluids, while others prefer either lattice-Boltzmann-based passive scalar solvers or classical FV/FD solvers for the balance equations, leading to a hybrid formulation.
While state-of-the-art lattice Boltzmann solvers are now routinely used for complex combustion configurations involving complex geometries and turbulent flows, a number of technical challenges still persist: • One of the remaining challenges is to get exactly conservative curved boundary conditions for the lattice Boltzmann solver.While the bare half-way bounce-back method resulting in a stair-case approximation to the geometry [134,135] ensures mass conservation, all curved treatment presented in the literature, see for instance [136][137][138], result in loss of conservativity of the boundary condition.For a more detailed discussed on conservation issues for curved boundary treatment interested readers are referred to [29,139,140].A number of routes can be taken to overcome this issue such as the use of immersed boundaries which would come at the cost of diffuse interfaces, or the use of volumetric/fluxed based boundary treatments, see for instance [141,142].
• Development and implementation of conservative and efficient dynamic grid-refinement strategies is also another topic to be further developed in the future.Although grid-refinement in the context of the lattice Boltzmann method has been developed and used since the early 2000's, see for instance [143][144][145][146][147], mass-conservation, spurious currents at refinement interfaces and dynamic refinement are still topics of discussion in the literature.
At the end of this detailed review regarding past achievements obtained for combustion thanks to LB simulations, it is now time to look briefly to the future.This work is obviously not the first review concerning lattice Boltzmann simulations, and previous studies have sometimes included long-term perspectives, most prominently in the work by Succi [148] -recently updated in [149].It is obviously necessary to evaluate again recent evolutions in the light of those predictions.One must keep in mind that applications are in the focus of the present review, as stated in the title, so that it is does not appear meaningful to include exceedingly exotic concepts here.
In the same manner, the future of high-fidelity combustion simulations has been discussed in previous reviews, for instance [150][151][152].Here also, reflecting on the corresponding statements will be useful.Of course, the aspects already discussed at length in the core of the present article will not be repeated here in the interest of space.Since the focus has been placed on important methodological aspects of LB in this review, a bird's view seems more appropriate to finish.As such, the main points emerging for the foreseeable future would read as follows.
• Coupling of different model approaches: even after having reviewed the many advantages of LB in this article, it remains clear that alternative model formulations may be very attractive from a numerical or an algorithmic point of view for specific applications.The coupling of LB to finite differences or finite volumes has been extensively described in the section concerning hybrid models.In a similar manner, coupling LB to a variety of other approaches is feasible and attractive for corresponding problems, also when involving reacting flows.Integrating the Immersed Boundary Method (IBM) within a LB framework is particularly relevant for reactive particles [153].To take into account even more complex physical interactions within particulate flows, a coupling of LB with the Discrete Element Method (DEM) is the next logical step [154].For configurations involving clearly separated fluid regions, adding a Volume of Fluid (VOF) approach -or related methods like Volume of Pixel -appears as a promising solution [155].
• Multi-physics and multiscale applications, adaptivity: the growing performance of existing computational platforms coming with the ever-increasing complexity of the target applications, problems involving a variety of physical and chemical processes -mostly coupled in a strongly non-linear manner -and taking place over a broad diversity of scales in time and space progressively become the rule.While concentrating here on combustion applications, it must still be recognized that LB has now virtually been used successfully for almost any kind of flows (and even beyond fluid mechanics).
The next step -certainly soon to be reached -will for example involve accurate LB simulations of turbulent multiphase reacting flows including realistic chemistry, thermodynamics, and transport models for species and heat, up to radiative heat transfer.
Such multi-physics configurations typically involve a broad range of relevant scales in time and space [156], from micrometers to centimeters for living beings [157], or even "from inside protons to the outer Universe", citing the optimistic statement of [149].In that case, an optimal combination of different numerical approaches will become necessary to allow for a numerical solution of the full configuration, for instance by coupling LB to Molecular Dynamics simulations [158][159][160].Multiscale issues can be partly mitigated by using adaptivity -in particular in space for LB (local grid refinement/coarsening), a solution that has already been discussed in this review [122,161].
• Integration of Machine Learning: the incredibly fast development of Machine Learning (ML) approaches will obviously not leave LB untouched, since it is increasingly supporting all approaches used for flow simulations.Regarding more specifically LB for combustion applications -usually taking place in the turbulent regime -two lines of research will certainly be followed in the near future.One will concentrate on the thermoreactive part of the problem, for instance using Deep Neural Networks to describe kinetics, potentially leading to very large savings in computing time and memory [162].The other line will concentrate on ML at subgrid scale (SGS) when using LB for Large-Eddy Simulations (LES) of turbulent reacting flows [122,128].In that case, either the behaviour of the pure turbulent flow will be described by a SGS model based on ML [163] -an approach now well established in conventional CFD [164]; or ML could be used to represent directly turbulent/combustion coupling at subgrid scale, an even more challenging but potentially more rewarding solution, for which much remains to be done [165].
• New computational architectures: Though this might not be immediately clear for young scientists, the performance of a numerical approach is not constrained only by considerations from applied mathematics (stability, convergence, dissipation, dispersion), but is also directly impacted by the details of the computational architecture on which large-scale simulations are carried out.In that sense, the comparison of different methods in terms of computational performance completely depends on the employed system.While method 1 might be order-of-magnitude faster than method 2 on a conventional, single-core system, the comparison might be completely reversed on a large, fine-grain parallel computer, or when computing on Graphical Processing Units (GPU).In that sense, an essential advantage of LB (in addition to the error-free streaming step and linearity of the basic operator) is its locality.In the standard LB formulation, only first-neighbour communications are requested, making it perfectly suited for the currently employed computer architectures; this is one essential explanation to understand the growing success of LB for many applications.Nobody knows which computer architecture will dominate in 20 years.In his review from 2015, Succi [148] already mentioned the suitability of LB for Quantum Computing (QC).Indeed, porting high-order CFD methods involving unstructured grids on QC systems sounds like a nightmare, and LB could here again profit from its apparent simplicity and locality.Lattice Boltzmann on Quantum Computers is a subject of current research [166,167].Still, QC systems being currently barely available for researchers, the impact of Quantum Computing on future LB simulations cannot be reliably estimated.The same applies to even more exotic architectures like biological computers, that have not even entered a preliminary test-phase.

Figure 1 :
Figure 1: Illustration of the tensorial product process to build D2Q9 lattice from D1Q3.

Figure 2 :
Figure 2: Illustration of applications of the model of Eq. (55).(Left) Sound pressure field for shock-vortex interaction with advection Mach number of 1.2 and vortex Mach number set to 0.25 at t * = 6.(Right) Iso-surface of velocity divergence colored by local Mach number for compressible decaying turbulence at Ma= 0.5 and t * = 0.4.Images are reproduced from [27].

3. 1 . 4 .
Hybrid models: Finite difference and finite volume solvers for energy and species Hybrid models are closest to multiple distribution functions approaches.In multiple distribution functions, let us remind that ρE (or an alternative energy form) and ρY k correspond to the zeroth order of separate distribution functions, see Eq.(46).When the number of species to be considered becomes large -typically O(10−100), the memory required to solve all scalars increases very quickly, which may become prohibitive for detailed chemistry descriptions.Hybrid models reduce the memory load by introducing a single scalar for (ρE, ρY k ) instead of the number of discrete velocities.Each additional conserved scalar ρϕ (where ϕ may represent E or Y k ) is solved by classical finite-difference or finite-volume (FD/FV) schemes, while continuity (1) and momentum (2) equations are still solved via their associated distribution function f i .

Figure 8 :
Figure 8: Drag coefficient c d as a function of the free stream Mach number for the Busemann biplane simulations.Inset: snapshots of the pressure distribution around the biplane for three different Mach numbers: Ma = 1.5, top; Ma = 1.7, bottom left; Ma = 2.0, bottom right.Figure reproduced from [102].

)
It was shown through the Galilean-invariance of the moments of the Maxwell-Boltzmann distribution function and the binomial theorem that the weights are also Galilean-invariant and therefore untouched by the change of reference frame.The immediate consequences of that observation are: (a) construction of 3-D lattices via tensorial product of the 1-D lattice remains as before, (b) assuming U = kδx/δt with k ∈ Z the propagation remains on-lattice and (c) the discrete entropy functional is Galilean invariant and therefore equilibrium populations are form invariant under the shift of reference frame.This point along with the effect of the shift on operation range has also been discussed for standard isothermal lattices[109].The process of changing the reference frame and resulting discrete lattice is illustrated in Fig.7through the D2Q49 lattice.The use of the shifted lattice along with the entropic equilibrium has been successfully used to model a wide variety of high Mach number flows as illustrated in Fig.8.The idea of shifted reference frames was later generalized to local adaptive reference velocity and temperature through the particles on

Figure 10 :
Figure 10: Linear stability domain of lattice Boltzmann at different non-dimensional temperatures and viscosities as obtained from von Neumann analysis.Reproduced from [112].

Figure 11 :
Figure 11: Streamlines of the 2-D circular flame simulation colored by velocity magnitude (in m/s).Left column: densitybased model [88], right column: pressure-based model [87].Note the very different ranges of velocity magnitude from the two methods.The yellow contour is the heat release rate peak indicating the flame front.

Figure 12 :
Figure 12: Illustration of successful applications of the generic HRR model (143): (left) a cycle of a thermo-acoustic instability in the PRECCINSTA burner[122], and (right) 2-D detonation cellular structure with varying activation energy[93].

3. 2 . 4 .
Low Mach thermo-compressible pressure-based solver In 2019, Hosseini et al. proposed a low Mach lattice Boltzmann scheme for simulations of combustion,

Figure 13 :
Figure 13: Illustration of some of the recent applications of the low Mach model of [125].(a) Simulation of PRECCINSTA swirl injector at equivalence ratio 0.83 [128], the red surface illustrating the flame structure.(b) Simulation of deflagrating flame in a chamber with obstacles [128].(c) Simulation of flame propagation in randomly-generated porous media [130].

Table 2 :
A list of canonical combustion problems treated using lattice Boltzmann methods in the recent literature