Compressible lattice Boltzmann methods with adaptive velocity stencils: An interpolation-free formulation
CCompressible lattice Boltzmann methods with adaptive velocity stencils:An interpolation-free formulation
C. Coreixas a) and J. Latt b) Department of Computer Science, University of Geneva, 1204 Geneva, Switzerland (Dated: 13 October 2020)
Adaptive lattice Boltzmann methods (LBMs) are based on velocity discretizations that self-adjust to local macroscopicconditions such as velocity and temperature. While this feature improves the accuracy and the stability of LBMsfor large velocity and temperature fluctuations, it also strongly impacts the efficiency of the algorithm due to spaceinterpolations that are required to get populations at grid nodes. To avoid this defect, the present work proposesnew formulations of adaptive LBMs which do not rely anymore on space interpolations, hence, drastically improvingtheir parallel efficiency for the simulation of high-speed compressible flows. To reach this goal, the adaptive phasediscretization is restricted to particular states that are compliant with the efficient “collide and stream” algorithm,and as a consequence, it does not require additional interpolation steps. The development of proper state-adaptivesolvers with on-grid propagation imposes new restrictions and challenges on the discrete stencils, namely the need foran extended operability range allowing for the transition between two phase discretizations. Achieving the minimumoperability range for discrete polynomial equilibria requires rather large stencils (e.g. D2Q81, D2Q121) and is thereforenot competitive for compressible flow simulations. However, as shown in the article, the use of numerical equilibria canprovide for overlaps in the operability ranges of neighboring discrete shifts at acceptable cost using the D2Q21 lattice.Through several numerical validations, the present approach is shown to allow for an efficient realization of discretestate-adaptive LBMs for high Mach number flows even in the low viscosity regime.
I. INTRODUCTION
The Boltzmann equation (BE) describes the space and timeevolution of the velocity distribution function f ( x , ξ , t ) . Thelatter accounts for the number of fictitious particles at the lo-cation x , time t , that are propagating with a given velocity ξ : ∂ t f + ξ α ∂ α f = Ω ( f ) , (1)where Greek letters stand for space coordinates, and Ein-stein’s summation rule is implied. Roughly speaking, the BE(1) illustrates the balance between transport and collision ofparticles. Assuming that f is close to its equilibrium state f eq ,the collision term can be approximated by a relaxation mech-anism Ω ( f ) ≈ − f − f eq τ , (2)which makes f tend towards f eq under a relaxation time τ ,as originally proposed by Bhatnager, Gross and Krook, andmore commonly known as BGK collision model . Even if thelatter operator leads to a Prandtl number of unity and does notaccount for extra internal (rotational and vibrational) degreesof freedom, simple extensions exist to allow for the modelingof polyatomic gases with flexible Prandtl number .Since the 1950s, several types of deterministic solvers ofthe BGK-BE were proposed . For efficiency reasons, mostof them rely on a (physical) discretization of the phase space,and can then be encompassed in the framework of discrete a) Corresponding author: [email protected] b) Electronic mail: [email protected] velocity methods (DVMs). The latter focus on the resolutionof the discrete velocity Boltzmann equation (DVBE): ∂ t f i + ξ i α ∂ α f i = Ω ( f i ) . (3)Contrary to the BE (1), the DVBE (3) is a set of partial dif-ferential equations of finite size V , where V is the number ofdiscrete velocities ξ i that compose the phase space of interest( i ∈ (cid:74) , V (cid:75) ).For DVMs, the phase space discretization is chosen in sucha way that all macroscopic velocity and temperature fluctu-ations of interest can be accounted for even in extreme con-ditions –which usually leads to large velocity stencils com-pared to the standards of LBMs. By further relying on Eu-lerian discretizations of the DVBE (finite-difference, finite-volume, etc), DVMs are able to simulate –in a stable manner–flows with strong departures from equilibrium that are en-countered, e.g., in rarefied conditions . LBMs on the otherhand also solve the DVBE (3), but contrary to DVMs, theyrely on the collide-and-stream algorithm which is a very effi-cient Lagrangian-based numerical scheme . The correspond-ing BGK-DVBE reads as f i ( x + ξ i , t + ) = f i ( x , t ) − τ f (cid:0) f i − f eqi (cid:1) ( x , t ) , (4)where lattice units with unitary time and space step are im-plied, and τ f = τ + /
2. In addition, the phase space dis-cretization is chosen according to quadrature or moment-matching rules. This aims at recovering the physics of inter-est in terms of moments , instead of velocity/temperaturefluctuations. Hence, the size of the discrete phase space di-rectly depends on the targeted physics. This is the reasonwhy standard LBMs –that are based on 19 (or 27) discretevelocities– are usually used as efficient alternatives to Navier-Stokes-Fourier (NSF) solvers for the simulation of isothermal a r X i v : . [ phy s i c s . c o m p - ph ] O c t
1. The redarrow/dot corresponds to the reference population used for the shift: (0,0), (-1,0) and (1,1) from left to right.weakly compressible flows , and also for more complexconfigurations . To conclude the overview of method-ologies, it is worth noting that discrete Boltzmann methods(DBMs) and discrete unified gas kinetic schemes (DUGKS),that share common features with both DVMs and LBMs, canalso be used for the simulation of out-of-equilibrium phenom-ena that are encountered, e.g., with multi-component mix-tures, reactive flows as well as multiscale rarefied gasflows . Like DVMs, DBMs and DUGKS are basedon Eulerian numerical schemes. Nonetheless, the moment-matching approach and quadrature rules are adopted (simi-larly to LBMs) for the derivation of the velocity stencil, whichallows for a drastic reduction of its size.All the aforementioned methodologies are usually based onstatic velocity discretizations. This means that if either strongout-of-equilibrium effects or important fluctuations of macro-scopic quantities are encountered in a simulation, the accuracyand the stability of the numerical solver is at risk, because thevelocity discretization does not match anymore the simulatedphysics. To tackle this issue, one can rely on adaptive velocitydiscretizations of the BE, which consist of sets of discrete ve-locities that self-adjust to local macroscopic quantities (suchas velocity u and temperature T ): c i α = ξ i α √ θ + u α , (5)with θ = T / T and T being the reduced and lattice (or ref-erence) temperatures, respectively. This adaptive methodol-ogy was introduced in the 1990s to increase the accuracy andthe stability of DVMs , as well as LBMs , in the con-text of high-speed compressible flow simulations. Interest-ingly, adaptive DVMs are interpolation-free solvers for theevolution of discrete populations due to their Eulerian nature.On the contrary, adaptive LBMs rely on a Lagrangian-basedscheme ( collide-and-stream ) meaning that populations are notnecessarily streamed from one grid node to another one, e.g.,if the discrete velocity components are non-integer. Hence,the algorithm must be supplemented with a space interpo-lation step to recover populations at grid nodes. This addi-tional step is computationally expensive, potentially prone toanisotropy issues (depending on the considered stencil of in- terpolation points), and deteriorates the natural parallel effi-ciency of LBMs. This reduces the interest one might have inmoving from NSF solvers to adaptive LBMs for the simula-tion of high-speed compressible flows.To get rid of the interpolation step, one must ensure that alldiscrete velocity components are integer valued in any flowcondition ( ξ i α , c i α ∈ Z ), which guarantees an on-grid stream-ing step. These conditions on the velocity/temperature shiftsof discrete velocities are illustrated in Fig. 1 for the D2Q21 lat-tice –whose characteristics can be found in the work by Zhanget al. . Yet, jumping from a lattice to another one in a dis-crete manner imposes new restrictions and challenges on thediscrete stencils to guarantee that the transition is achievedin an accurate and stable manner. Such an investigation isthe starting point of our work, which is organized as follows.In Section II, the viability of the transition between velocitydiscretizations is investigated by looking at operability rangeoverlaps for LBMs based on polynomial and numerical equi-libria. After finding a good compromise between stability,accuracy and efficiency, shifting criteria and interpolation-free strategies are discussed in Section III. The ability ofinterpolation-free adaptive LBMs to handle both linear (prop-agation of shear, entropy and acoustic waves) and non-linearregimes (1D and 2D Riemann problems) is then investigatedin the low-viscosity limit (Section IV). Conclusions are drawnin Section V. Eventually, a number of appendices are providedto help the reader further understanding concepts used in thiswork: (A) approximation of the Maxwellian in a given refer-ence frame, (B) linear stability analyses, (C) detailed descrip-tion of numerical equilibria, (D) preliminary study regardingthe choice of lattice and set of constraints, and (E) extensionto variable Prandtl number. II. REALIZABILITY OF PHASE SPACE TRANSITIONSA. Motivation
In the LB framework, the velocity space discretization ismost commonly based on stencils that are static (in space andtime), symmetric, and consequently, centered around a pop-ulation at rest. This implies that these lattices are dedicatedto the simulation of quasi-quiescent flows with little temper-ature fluctuations, i.e., ( u , T ) ≈ ( , T ) . This is further sup-ported by the fact that corresponding equilibrium distributionfunctions (EDFs) are usually based on polynomial expansionsof the Maxwellian about ( u , T ) = ( , T ) . When deviatingfrom the reference state ( , T ) , both the accuracy and the sta-bility of these LBMs are impacted, and this can be quanti-fied through macroscopic error evaluations and linear stabilityanalyses . In fact, one can already identify for a continu-ous phase space issues induced by: (1) the approximation ofthe Maxwellian and (2) the associated reference frame. Theinterested reader may refer to Appendix A for an in-depth dis-cussion regarding the latter point.To make sure the correct physics can be simulated in a sta-ble and accurate manner, one should adapt the direction andnorm of discrete velocities to the mean flow conditions , asshown in Fig. 1. In addition, one must ensure that the tran-sition between two states of the same lattice can be done ina stable and accurate manner. Henceforth, corresponding re-quirements for stable and accurate phase space transitions areinvestigated. They directly depend on how the LBM is de-rived, and more precisely, on the type of EDF used to recoverthe macroscopic behavior of interest. For LBMs based onpolynomial EDFs, it is proposed to perform linear stabilityanalyses to find the minimal lattice size allowing for stabil-ity domain overlaps between two shifted versions of the samelattice. Regarding numerical EDFs, the latter methodologycannot be used anymore, and instead, the convergence of theroot-finding solver is considered as an alternative criterion forthe evaluation of stability domain overlaps. B. Realization based on polynomial discrete equilibria
A variety of systematic approaches to approximate theMaxwell-Boltzmann equilibrium in a discretized phase spacehave been developed. Of these methods, two have becomerather popular : (a) truncated Hermite expansion and(b) moment-matching. While they both converge to the samediscrete equilibrium state for tensor-product-based stencils(e.g. DdQ3 d ), the latter is more flexible and allows to de-rive a discrete equilibrium for almost any stencil structure. Asclearly indicated by its name, it mainly consists in matchingthe moments of the discrete equilibrium state to those of theMaxwell-Boltzmann distribution up to the highest order of in-terest.
1. Polynomial EDF construction and corresponding lattices
For the scheme to correctly recover the targeted physics, thenumber of discrete (shifted) velocities in the stencil c i mustbe higher than the minimum number of constraints, i.e., mo-ments of the EDF needed to match their continuous counter-parts. As such the discrete equilibrium construction processcomes down to solving a system of algebraic equations of the following form: Gf eq = M MB , (6)where M MB is the constraints vector containing the con-tinuous moments of the Maxwellian, with the n th moment( n = p + q ∈ N ) defined as: M MB pq = (cid:90) c px c qy f MB d ξ . (7)and f MB = ρ ( πθ ) d / exp (cid:20) − ( ξ − u ) θ (cid:21) , (8)where d is the number of physical dimensions. Assuming G is invertible, the system of equations (6) can be solved, andthe discrete (shifted) equilibria are then obtained as: f eq = G − M MB . (9)For a one-dimensional physical space, the DVBE correctly re-covers the NSF dynamics if and only if the first five momentsof the EDF are exactly matched. If we restrict the study totensorial products of 1D velocity stencils, then the smalleststencils that satisfies these conditions are of the form DdQ5 d .A well-known illustration of this moment-matching require-ment is the third-order family of stencils (DdQ3 d ) widely usedin the LB community for low-Mach isothermal flows. The lat-ter restriction results from the number of degrees of freedomin the system that prevents these lattices to correctly recovermoments higher than two.
2. Assessment of operability range: linear stability domain
The operability range of a given DVBE solver, targetingthe Navier-Stokes-Fourier (NSF) system of equations, can beassessed following the Lax-Richtmyer equivalence theoremconditions, namely stability and consistency.The latter can be readily evaluated by looking at the de-viations of moments (appearing at the NSF level) of the dis-cretized equilibrium function from their continuous counter-parts. This tool can be especially useful for low order sten-cils that fail to impose all constraints tied to the NSF levelmoments. However, for larger stencils that explicitly enforcephysical constraints on all moments intervening at the NSFlevel, deviations reduce down to zero and are therefore use-less in assessing the operation range of the solver. While anecessary condition for the equivalence, consistency is notsufficient and must be completed by a proof of stability ofthe numerical scheme.The linear stability analysis (also called von Neumann orFourier analysis) is a rather popular and widely accepted ap-proach to assess the stability and accuracy of a given linearsystem of discrete equations . In the LB context, thistool was used for various comparative studies, and to quantifythe impact of several parameters: numerical discretization ofthe DVBE , spatial filtering , accuracy and efficiency com-parison with NSF solvers , optimization of multi-relaxation V − − u x T opt T opt V T U x = 0 U x = 1 U x = 0 U x = 1 FIG. 2: Minimal stability condition which allows adaptive shifting with LBMs based on polynomial equilibria. Stabilitydomain obtained for different 1D lattices of size V : (solid) U x =
0, and (dashed) U x =
1. Configurations T opt and T corresponds to the maximum achievable velocities at optimal reference temperature, and with at least 50 percent variation intemperature, respectively. The dotted lines represent the minimal stencil required for an overlap of ∆ u x ≥ . T opt and T respectively.time collision models , stabilization mechanism of colli-sion models , shifted stencils , etc.While intended for linear equations, it is also widely usedto assess the stability properties of non-linear systems, by re-placing the target equations with a linearized version . ForLBMs, this boils down to the linearization of the collisionterm since the convective part is already expressed in a lin-ear form. While -strictly speaking- the linearization limits thevalidity of the analysis to the linear regime, one can argue thatin practice this just changes the value of the outcome from sufficient to merely a necessary condition. Based on the latterassertion, the linear stability analysis is still relevant as it es-tablishes upper bounds for the operability range of the solver,making it especially interesting for higher-order stencils withpolynomial equilibria.As such, we establish a trend for the linear stability do-main as a function of the order of the stencil in the context ofpolynomial equilibria. Following the approach described in and briefly recalled in Appendix B, the analysis is performedin the { u , θ , τ } parametric space for 1D symmetric stencilsof different orders, both with and without shifts. The largestvelocity resulting in a stable system in the limit of vanish-ing viscosity is reported as the stability domain of the sten-cil. However it is common knowledge that while allowing forlarger Mach numbers, larger stencils result in narrower tem-perature ranges. As such a second stability domain, based onthe maximum velocity allowing for fifty percent variation intemperature is also reported. The obtained results for shiftedand non-shifted stencils are shown in Fig. 2. It can be read-ily observed that the smallest stencil even allowing for a sta-ble unit shift is the D1Q9 lattice ( ξ i α ∈ { , ± , ± , ± , ± } ).For the stencil to allow for a stable shift and -more or less-pronounced temperature variations the stencil size goes upto ξ i α ∈ { , ± , ± , ± , ± , ± , ± } , making it rather inef- ficient in terms of memory consumption, processing powerand communication overhead in parallel computations.Therefore, while our study of shifted stencils with polyno-mial equilibria remains of theoretical interest, practical appli-cation are better served with the approach of numerical equi-libria described below. C. Realization based on numerical discrete equilibria
While polynomial EDFs impose strong constraints onthe velocity discretization through the moment-matchingapproach, numerical EDFs allow for the derivation ofquadrature-free LBMs, hence providing more freedom regard-ing the size of the velocity stencil. During the past threedecades, the latter EDFs have been extensively used in thecontext of rarefied gas flow simulations for both static andadaptive phase space discretizations . Hereafter,they are introduced as interesting alternatives to polynomialEDFs for the transition between two velocity space discretiza-tions of smaller size.
1. Construction of quadrature free discrete equilibria
This kind of EDF results from the minimization of the H -functional H = ∑ i f i [ ln ( f i / a )] (10)under the constraints G pq = ∑ i f eqi ξ pix ξ qiy − M MB pq = . (11)and reads as f eqi = a exp [ − ( + ∑ p , q λ M MB pq ξ pix ξ qiy )] , (12)where λ M MB pq are the Lagrange multipliers corresponding to theconstraints (11). Following our previous work , the prefactor a = ρ is adopted hereafter.Assuming the set of constraints (11) corresponds to theconservation of mass, momentum and total energy (the ze-roth, first and trace of second-order Maxwellian moments),the Chapman-Enskog (CE) expansion gives : ∂ t ρ + ∂ χ ( ρ u χ ) = , ∂ t ( ρ u α ) + ∂ β ( ρ u α u β + p δ αβ ) + ∆ = ∂ β ( Π αβ ) + ∆ , ∂ t ( ρ E ) + ∂ α (( ρ E + p ) u α ) + ∆ tr = ∂ α ( Φ α ) + ∆ tr . (13) ∆ n and ∆ trn are errors that emerge because their correspond-ing constraints (moments of order n or their trace) are notaccounted for in the computation of the exponential equilib-rium (12). While diffusive errors ( ∆ and ∆ tr ) can usuallybe neglected at moderate or high Reynolds numbers , this isnot the case for deviations related to convective terms unlessone adopts large velocity discretizations. In order to deriveefficient and accurate LBMs based on the numerical equilib-rium (12), one of the following sets of constraints should beadopted :• 8-moment: ( , ξ i α , ξ i α ξ i β , ξ i χ ξ i χ ξ i α ) (14)• 9-moment: ( , ξ i α , ξ i α ξ i β , ξ i χ ξ i χ ξ i α , ξ i χ ξ i χ ξ i η ξ i η ) (15)• 10-moment: ( , ξ i α , ξ i α ξ i β , ξ i α ξ i β ξ i γ ) (16)• 11-moment: ( , ξ i α , ξ i α ξ i β , ξ i α ξ i β ξ i γ , ξ i χ ξ i χ ξ i η ξ i η ) (17)• 13-moment: ( , ξ i α , ξ i α ξ i β , ξ i α ξ i β ξ i γ , ξ i χ ξ i χ ξ i α ξ i β ) (18)The set of eight constraints (18) is the minimal configura-tion allowing to accurately simulate convective phenomena( ∆ = ∆ tr = M MB αβγ and M MB αβ χχ , one ends up with the 13-moment ap-proach (18) which exactly leads to the Navier-Stokes-Fourierlevel of physics in the continuum regime.
2. Operability range based on the convergence of theroot-finding algorithm and non-linear kinetic stabilization
Increasing the number of constraints automatically leads tobetter macroscopic properties. Nevertheless, this generallycomes at the cost of a reduced operating range in terms ofvelocity and temperature fluctuations. This is explained bythe fact that increasing the number of constraints puts moreeffort on the root-finding algorithm used for the computationof the numerical equilibrium. Depending on the type of al-gorithm considered (bisection, secant, Newton, etc), the ro-bustness of the resulting LBM can be strongly impacted. Al-gorithms based on exact formulations of the Jacobian matrixlead to wider stability ranges, and we noticed that the largerthe lattice size, the wider the stability domain of the LBM.Yet, a trade-off has to be made between lattice size and num-ber of constraints in order to obtain an accurate, robust andefficient LBM.In the early stage of this work, several lattices were stud-ied (D2Q9, D2Q13, D2Q17, D2Q21, D2Q49, D2Q81) incombination with a large number of constraints, i.e., M ∈{ , , , , , } . In the context of supersonic flow simula-tions, and following the methodology proposed in our previ-ous work , a good trade-off between accuracy, stability andefficiency was found with the D2Q21 lattice based on theEDF computed via M = ε = ε = (cid:12)(cid:12) M MB pq − M eqpq (cid:12)(cid:12) M MB pq , (19)as originally proposed by Kornreich and Scalo . To illus-trate this point, deviations obtained with the present LBM(D2Q21, M = T = .
7) are com-piled in Fig. 3 for two moments related to diffusive phenom-ena, namely, the third-order moment M MB xxx = M MB30 = ρ u x ( u x + T ) , (20)and the trace of the fourth-order moment M MB xx αα = M MB40 + M MB22 = ρ [( E + T ) u x + T ( E + T )] . (21)One can observe that error levels remain under the 10%threshold for a large range of Mach numbers. In addition, byadapting the lattice in a discrete manner, it is confirmed thatthe root-finding solver can converge for even higher values ofthe Mach number while leading to acceptable deviations withrespect to the targeted physics. More importantly, large over-lapping zones are now present for several shifts of the D2Q21lattice which would solve the problem encountered with poly-nomial equilibria, i.e., the need for large lattices to allow thetransition between reference states. u x − − − − ε M xxx u x M xxαα U x = 0 U x = 1 U x = 2 FIG. 3: Velocity and shift impact on macroscopic errors ε for the D2Q21 lattice with 8-moment approach and T = .
7. Theflow propagates along the x-axis, and the grey zone starts at ε = – thatcombines the fast convergence of Newton’s approach with agradient descent (see Chaps 6 and 7 in for more details re-garding this algorithm). Eventually, to further improve thestability in case of strong departure from equilibrium and/ordue to under-resolved conditions, we propose to also includea non-linear stabilization technique that shares some similar-ities with flux limiters . Its main characteristics are brieflyrecalled in Section IV A 2. III. COUPLING STRATEGIES FOR SHIFTED LBMS
In the previous section, the possibility to switch betweendifferent versions of the same lattice in a stable manner hasbeen detailed. We now focus our discussion on the two maintechnicalities related to the coupling of different LB solvers:(1) when to switch between the two solvers, and (2) how toexchange the information at the interface between them. Theformer will first be tackled by discussing local criteria for theshift of velocity space discretization. The second issue will bedealt with investigating the different strategies for the recon-struction of missing populations at the interface between twostates of the (same) velocity discretization. This discussion isconducted in the context of polyatomic gas flow simulationsin which rotational (and vibrational) degrees of freedom areaccounted for through a second set of populations g i ,as discussed in Section IV A 1. A. Criteria for lattice shifting
Hereafter, we focus on criteria that can be used to decidewhich version of the lattice better fits the physics of interest ata given node x and time t .
1. Breakdown criteria
In the literature dedicated to the coupling of LBMs withother kinetic based solvers, e.g., for rarefied gas flow simula-tions (DSMC , DVMs or high-order LBMs ),valuable information is found regarding static and dynamicdomain decomposition relying on particular criteria such asthe breakdown of the continuum assumption.In the static case, more advanced solvers are used close towalls in order to capture out-of-equilibrium effects, such asvelocity and temperature jumps, which cannot be simulated bycontinuum based solvers (NSF or low-order LBMs), at least,for Knudsen numbers close to unity. The domain decomposi-tion is then designed manually, as similarly done for mesh re-finement. . The latter approach assumes an a priori knowl-edge of domain portions related to out-of-equilibrium phe-nomena. When this information is unavailable, a dynamic do-main decomposition criterion is necessary. Most breakdowncriteria are based on the (local) Knudsen number, gradientsof macroscopic quantities, or deviations with respect to (non-)equilibrium moments/populations .A criterion for the local shift of velocity stencils, as targetedin this paper, seems however more straightforward to formu-late, as it is sufficient to formulate bounds on the macroscopicquantities of interest.
2. Criteria based on macroscopic quantities
In the context of numerical EDFs, it is possible to evaluateat the same time deviations from the macroscopic behavior ofinterest, as well as stability limits, through studies similar tothose conducted in Fig. 3. This is the starting point for oursearch for adequate phase space transition criteria.The first criterion one might think of is a dimensionless one,namely, the Mach number Ma = u / √ γ r T , with u the norm ofthe velocity vector u , and γ r the specific heat ratio. Althoughthis criterion is rather convenient, as it is based on a dimen-sionless quantity, it fails to reveal the true nature of instabil-ities. For example, a simulation based on the D2Q21-LBMwith numerical equilibirum (M=8) shows that numerical sta-bility is limited by a range of lattice-unit velocities rather thanthe Mach number. Indeed, assuming a uniform flow propagat-ing along the x-axis at Ma = .
85 with T = . θ =
1) and γ r = /
3, the root-finding solver becomes unstable for u x ≥ T to 0 . θ = = T is limited by the stabilityrange of the root-finding solver. This observation could leadto the desire to adjust the reference temperature from a portionof the domain to another, just as it is done for the reference ve-locity of the lattice. However, temperature fluctuations remainsufficiently low in the conducted experiments ( T max / T min < n − . < u α ≤ n + . ⇐⇒ U α = n (22)where n ∈ Z . For the D2Q21-LBM based on the numericalEDF with M =
8, the above intervals ensure that the previousand the updated lattices will share similar macroscopic devi-ations. This is justified, once again, through analyses that areconducted in Fig. 3 where error levels are almost identical for U x = u x = . U x = u x = . ( U x , U y ) =( , ) , ( − , ) and ( , ) for the following macroscopic con-ditions ( − . , − . ) < ( u x , u y ) ≤ ( . , . ) , ( − . , − . ) < ( u x , u y ) ≤ ( − . , . ) and ( . , . ) < ( u x , u y ) ≤ ( . , . ) , re-spectively.Eventually, to reduce the number of interface cells, or thefrequency of change of reference velocity, it can be usefulto use velocity ranges of a width larger than 1 (e.g., a range [ − , ] for U α = U α is dictated either by its history (to avoid changes of refer-ence) or by the cell’s neighborhood (to reduce the number of interface cells). Such enlarged ranges can only be achievedwith LBM models that are stable over a sufficiently large ve-locity range. B. Population reconstruction strategies at shift interfaces
Now that criteria have been defined to allocate differentvelocity stencils to each part of the simulation domain, aproper algorithm for the transfer of information across the in-terface, between two areas with different stencils, must be de-fined. Hereafter, all reconstruction strategies take place duringthe streaming step, which is considered hereafter in a ‘pull’manner (i.e., from the perspective of the cell receiving thestreamed population) to ease our discussion without any lossof generality. In the conventional LB algorithm, a populationcan only be pulled from a location x + c i α to the current node x if and only if it belongs to the same velocity discretization.Otherwize, the population is considered to be missing and areconstruction step is then required, as already done, for ex-ample, in the context of boundary conditions.
1. Existing strategies
In first works on adaptive LBMs by Sun et al., missingpopulations are approximated by their equilibrium value .While it is computationally cheap and general, this method islimited by the fact that a higher number of discrete velocitiesis required to recover the correct physics, which is in contra-diction with our goal to propose efficient adaptive LBMs.The above reconstruction strategy was recently improvedby computing post-collision populations through their mo-ment space : M (cid:48) h ∗ , (cid:48) = M h ∗ . (23) h ∗ , (cid:48) is the vector of (missing) post-collision populations in thepost-streaming velocity discretization, and h ∗ = h − M − S M ( h − h eq ) (24)is the vector of (known) post-collision populations in the pre-streaming velocity discretization. M and S are the mo-ment and relaxation matrices in the pre-streaming velocitydiscretization, whereas M (cid:48) is the moment matrix in the post-streaming one. This approach is the most expensive strategyto reconstruct missing populations since it requires the com-putation of all moments of h i = f i or g i . Even if the size ofthe moment space remains small for the most basic velocitydiscretizations, such as D2Q9, the number of floating pointoperations required by this reconstruction grows as V , where V is the size of the lattice. In addition, the full set of momentsis only defined in an unique manner for tensor-product-basedlattices (D2Q9, D2Q25, D2Q49, D2Q81, etc), whereas addi-tional steps are required to derive full moment sets for morecompact lattices (D2Q13, D2Q17, D2Q21, etc).One possible way to improve the efficiency would be tocompute missing populations via their Hermite polynomialexpansion, as already done, e.g., for hybrid solvers (DVM-LBM) in the context of rarefied gas flow simulations .This was very recently proposed through the adoption ofa standard (non-recursive) regularized approach in theplace of the reconstruction in the moment space . This al-lowed Zipunova et al. to reduce the computational overheadby one order of magnitude . This latter result is explained bythe fact that standard regularization steps discard high-ordermoments to reconstruct populations, hence reducing both thenumber of floating point operations and memory consumptionat the same time – as originally highlighted by Ladd and Ver-berg for standard LBMs . Even if the latter approach canbe applied to more compact lattices, it is based on formulasderived in the context of quadrature-based LBMs, and conse-quently, it cannot be used, as it is, in the present quadrature-free formalism.
2. From Chapman-Enskog to Grad-type reconstructions
In order not to restrict ourselves to tensor-product based lat-tices, we compute post-collision populations through the de-composition h ∗ i = h eqi + ( − / τ h ) h neqi , (25)where the BGK operator is adopted for the sake of simplicity,and τ h is the relaxation time corresponding to populations h i .The interest reader may refer to Appendix E for its extensionto variable Prandtl numbers.For (quadrature-free) LBMs based on numerical EDFs,while h eqi can be computed through the root-finding algorithmwith respect to the lattice of interest, h neqi is not directly avail-able. One possible way to compute it is to rely on the CEexpansion truncated at the first-order in Knudsen number : h neqi ≈ h ( ) , CE i = − τ h [ ∂ t h eqi + ξ i α ∂ α h eqi ] , (26)which is the most general form of regularizationsteps , in the sense that it does not rely onany assumption regarding high-order moments or theirrelaxation frequency. Instead, it requires the computationof the space and time evolution of EDFs. This is achievedby evaluating space and time derivatives through standardapproaches such as finite-difference (FD) approximations.This may strongly impact the accuracy of the transition strat-egy because errors at the level of populations might inducedrastic discrepancies at the macroscopic level. Regarding theefficiency of this reconstruction, the computation of the timederivative requires the storage of equilibrium populations atprevious time steps. Even in the most optimistic scenario(first-order Euler FD approximation), this strategy roughlydoubles the memory consumption of the resulting LBM,leading to a decrease of parallel efficiency, especially in thecontext of GPU acceleration.The last strategy considered in this work originates from thekinetic theory of gases, and is based on Grad’s description ofpopulations h i ≈ h eqi + h ( ) , Grad i = h eqi ( + φ h ) (27) where φ h is a deviation that accounts for (small) non-equilibrium contributions of populations h i . These deviationscan either be computed in terms of Hermite coefficients or through the maximum entropy principle . In both cases,the same form is obtained for monatomic gases: φ f = σ αβ c i α c i β ρ T + q α c i α ρ C p T (cid:32) c i χ T − C p (cid:33) , (28)where the traceless viscous stress tensor reads as − σ αβ = µ [ ∂ α u β + ∂ β u α − ( / D ) ∂ χ u χ δ αβ ] (29)and Fourier’s heat flux is q α = − κ∂ α T , (30)with c i α = c i α − u α shifted peculiar discrete velocities. µ = ρν , µ b = ( / D − / C v ) µ , ν and κ are the dynamic viscosity,bulk viscosity, kinematic viscosity and thermal conductivitycoefficients. Regarding the polyatomic non-equilibrium cor-rection imposed through population g i , one possible formula-tion is φ g = q α c i α ρ C p . (31)The latter only impacts the definition of the heat flux, andit naturally results from Rykov’s model for polyatomic gaseswhen non-elastic contributions are neglected .Similarly to the reconstruction based on the CE expan-sion (26), this methodology is lattice-independent, and con-sequently, more flexible than the reconstruction through themoment space (23). In addition, while the (numerical) equi-librium part h eqi is computed via macroscopic quantities ex-pressed in the lattice of interest, gradients are evaluatedthrough standard finite-difference (FD) approximations. Thereconstruction of missing information can then be done in anefficient and rather local manner, that is compliant with HPCarchitecture.All of this makes Grad’s approximation (27) a very appeal-ing reconstruction technique from the theoretical and practicalviewpoints. This strategy is therefore applied in the numericalsection of the article. IV. NUMERICAL APPLICATIONS AND VALIDATIONA. Implementation details
1. Underlying LB scheme
In the context of LBMs, the DVBE (3) is numericallydiscretized using the very efficient collide-and-stream algo-rithm (4), which can be divided into two consecutive steps:- Collision (here with the BGK approximation): h ∗ i ( x , t ) = h i ( x , t ) − τ h (cid:0) h i − h eqi (cid:1) ( x , t ) , (32a)- Streaming (pull): h i ( x , t + ) = h ∗ i ( x − c i , t ) , (32b)where the latter step has been formulated, as usual for adap-tive LBMS, from a "pull" perspective. . Here, adouble distribution function (DDF) formulation of LBMs isused ( h i = f i or g i ) to account for internal degrees of freedomof molecules via g i , hence, leading to a flexible specific heatratio γ r .To allow transitions between shifted versions of the D2Q21lattice, 8-constraint numerical EDFs (C29) are considered inthe rest of the paper (see Appendix C for more details). Inter-estingly, the equilibrium of the second set of populations ( g eqi )is not computed via the root-finding solver, but instead, it isobtained from its monatomic counterpart ( f eqi ) through g eqi = ( C v − D ) T f eqi , (33)with the heat capacity at constant volume being C v = / ( γ r − ) . Consequently, the polyatomic behavior is obtained at alow cost in terms of floating point operations, but at the costof doubled memory needs.Due to the BGK approximation of the collision model, theabove DDF-LBM is restricted to the simulation of flows witha unity Prandtl number, i.e., Pr = ρ C p ν / κ =
1. This can becorrected adopting a more advanced collision model, such asShakhov’s or Ryjkov’s . Both can be deduced from Grad’sformulation (27) by discarding the term related to stresses inthe definition of φ f (28). More details about this extension canbe found in Appendix E.Eventually, Grad’s reconstruction technique (27) will alsobe used as an extension to the compressible case of (regu-larized) initial and boundary conditions , using second-order FD approximations to approximate the viscous stresstensor σ αβ and the heat flux q α . For boundary conditions, itis worth noting that one could also rely on the bounce-back orthe extrapolation of non-equilibrium populations for the com-putation of σ αβ and q α . Nevertheless, an extensive compari-son of these approaches is out of the scope of this work.
2. Non-linear stabilization technique
By relying on numerical EDFs, the above collision mod-els have been shown to be more stable than their polynomialcounterparts for compressible flows . Yet, this is not suf-ficient to obtain stable simulations (1) in the low-viscosityregime, (2) for severely under-resolved conditions, and (3)when strong compressibility effects (shock waves) are en-countered. Following best practices in the CFD community,one should rely on stabilization techniques that are able toaccurately capture sharp gradients and discontinuities with-out suffering from Gibbs oscillations while keeping smoothregions of the flow intact . To that end, the kineticstabilization methodology proposed in our previous work isadopted, in order to tackle both points in a fairly good (but notperfect) manner . The latter consists in locally evaluating the departure from equilibrium through an approximation of theKnudsen number ε Kn = V V − ∑ i = | f i − f eqi | f eqi , (34)which allows the distinction of all the features encounteredduring simulations. When ε Kn > .
01, the recovery of theproper macroscopic behavior is at risk, and dissipation mustbe locally added to damp phenomena related to departuresfrom equilibrium (e.g., shockwaves) or Gibb’s oscillations in-duced by under-resolved conditions. This simple yet pow-erful kinetic sensor is coupled with the common BGK colli-sion operator through the dynamic relaxation time τ f ( ε Kn ) = τ f α ( ε Kn ) . Interestingly, this stabilization mechanism can alsobe seen as a limiter for changes induced by the collision term.The interested reader may refer to the work by Gorban andPackwood (and therein references) for other types of non-equilibrium limiters .In addition to its compliance with parallelism paradigms(local and fast evaluation), this stabilization technique wasshown to lead to stable and accurate simulations of the invis-cid Sod shock tube, and viscous flows past a 2D airfoil and 3Dsphere in the supersonic regime . Eventually, this stabiliza-tion technique barely depends on the chosen phase discretiza-tion, and is independent of the considered way of computingthe equilibrium (polynomial or numerical). It is therefore agood candidate to further increase the stability of the proposedadaptive LBMs when needed be.
3. Algorithm overview
In summary, a time iteration of the proposed algorithmpresent itself as follows:- Domain decomposition for the application of a local ve-locity space discretization (22)- Computation of monatomic numerical EDFs f eqi (C29)based on the constraints (14)- Inclusion of polyatomic contributions to numericalEDFs through g eqi (33)- If needed, computation of Knudsen-dependent relax-ation times based on the sensor (34)- Computation of post-collision populations (32a)- If the same lattice is used at x and x − c i then the nor-mal streaming step is applied (32b). Otherwise, miss-ing post-collision populations h ∗ i at x are reconstructedusing Eq. (27), where macroscopic quantities and theirgradients are evaluated at x − c i .It is interesting to point out that, the present adaptive ap-proach barely modifies the standard collide-and-stream algo-rithm and, most importantly, does not require any space in-terpolation. Extra steps consist in (1) the evaluation of thedomain decomposition for the velocity discretization, and (2)0reconstruction of missing post-collision populations when lat-tices at x and x − c i differ.The remainder of the paper is dedicated to the validationof the present approach, in terms of accuracy and robustness,for the simulation of linear and non-linear phenomena in thelow-viscosity regime (Sections IV B and IV C respectively). B. Preliminary investigations: Shifting impact on linearproperties
It is first proposed to investigate the linear behavior of thepresent approach for a wide range of flow velocities. Usu-ally, this is done via linear stability analyses ,which canhowever not been carried out presently due to the numericalnature of the EDF. Alternatively, one can directly simulatethe space-time evolution of small perturbations (shear, ther-mal and acoustic waves) that are superimposed to a uniformflow.Hereafter, our adaptive formulation will be validated us-ing the D2Q21 lattice based on the numerical EDF computedthrough M = γ r = . T = . ρ =
1. Allof this is done in the inviscid context ( ν =
0) to further high-light the accuracy and stability of the present approach. If nototherwise stated, reference data are obtained using the D2Q49lattice based on the numerical EDF computed through M = ( γ r , T ) = ( . , . ) .Due to the linear nature of the following flow configura-tions, the lattice self adjusts at the beginning of the simulationand does not change afterwards. It then restricts the recon-struction (27)-(31) to the initialization step, which highlightsits interesting properties as an extended initial condition. Allboundary conditions are periodic.
1. Shear wave
We first investigate the transport of a shear wave by an in-viscid ( ν =
0) uniform mean flow. This wave corresponds to asmall sinusoidal perturbation (in the transverse velocity field)that is superimposed to a uniform mean flow : ρ = ρ , u x = u , u y = A sin ( π x / L x ) , T = T , (35)where the perturbations’ amplitude is A = . ( ρ , T ) = ( , . ) , and the flow velocity u is a parameter that is varied to evaluate the accuracy and ro-bustness of the present approach in high-speed flow condi-tions. The simulation domain is quasi-one-dimensional, i.e. [ L x × L y ] = [ × ] , where L α is the number of points indirection α = x or y .In order to quantify the numerical dissipation and disper-sion of the proposed approach, it is mandatory to exactly prop-agate the shear wave over the same distance for any value ofthe mean flow velocity u . This is done by carefully choosing u and the time ( t = Nt c ) at which results are compared, where . . . . . . x/L x − . − . − . − . . . . . . u y / A u = 2 . u = 2 . u = 3 . FIG. 4: Convection of a shear wave in an inviscid flow with agrid mesh composed of L x =
100 points. Normalizedtransverse velocity fluctuations u y / A are plotted after t = t c for three different flow velocities u =
2, 2 . u = L x =
500 points.the characteristic time is defined as t c = L x / u , and N ∈ N ∗ .By taking u ∈ { . , . , . } , it is possible to compare all re-sults at t = t c , which corresponds to an integer number ofiterations for all values of u (i.e., t = L x = ν num = | max ( u ref y ) − max ( u y ) | max ( u ref y ) , (36)is 0 . .
9% and 0 .
7% for u =
2, 2 .
2. Entropy spot
The linear properties of our approach are further investi-gated through the transport of an (inviscid) entropy spot. Thelatter is initialized as a (Gaussian) hot spot which is superim-posed to a uniform velocity mean flow at constant pressure.Such an initial state was proposed by Fabre et al. , and re-cently investigated in the context of hybrid LBMs : ρ = ρ [ − A exp ( − r )] , T = T [ + A exp ( − r )] , (37)with u x = u , u y = ρ = T = . A = . r =[( x − x c ) + ( y − y c ) ] / R where ( x c , y c ) are the coordinates1 . x/L . y / L Reference u = 2 . u = 2 . u = 3 . . . . . . . . . . . . . x/L . . . . . . ∆ T / ( T A ) u = 2 . u = 2 . u = 3 . FIG. 5: Convection of an entropy spot in an inviscid flow with a grid mesh composed of L =
100 points in each direction:global view and slice along the x -axis (from left to right). Normalized temperature fluctuations ∆ T / ( AT ) are plotted at t = t c for three different flow velocities u =
2, 2 . L =
100 and 500 points for the global view and the slice respectively. The numerical dissipation with respect to thetheoretical temperature peak is 1 . .
9% and 1 .
2% for u =
2, 2 . R = L /
10 is related to the spot width,and L is the characteristic length of the simulation domain. Inaddition to the evaluation of the spectral properties (dispersionand dissipation) related to thermal waves, this benchmark testallows us to further quantify possible isotropy issues inducedby our approach.Here, the entropy spot is convected using u =
2, 2 . t = t c , as for the previous testcase. Results arecompiled in Fig. 5 for a relatively coarse mesh: L x = L y = L =
100 points, leading to a full width at half height of 2 R =
20 points. All configurations prove that the present approachcan transport small thermal fluctuations over long distancesat the correct speed (negligible dispersion error), with onlylittle loss of information (about 1% of the peak amplitude forall configurations), while keeping the isotropic nature of thetemperature field.
3. Acoustic pulse
The last of these linear testcases deals with the generationand propagation of (isentropic) acoustic waves. The latterare induced by a (small) pressure disturbance that is super-imposed to a mean velocity flow field: P = P [ + A exp ( − r )] , u x = u , u y = ρ = ρ ( P / P ) / γ r , T = T ( ρ / ρ ) γ r − , (38b)with P = ρ T . In addition to being convected by the meanflow velocity, this perturbation evolves over time and space. Hence, the comparison of the pressure field at the same phys-ical time for various values of u is not as straightforward asbefore. We then propose to recenter the pressure field around( x c , y c ) after a fixed number of iterations, so that, we will stillbe able to quantify numerical dissipation and dispersion of theproposed approach.Corresponding results are plotted in Fig. 6, for three meanvelocities ( u = .
3, 1 . . t / t c ≈ .
297 where the acousticcharacteristic time is defined as t c = L / √ γ r T and L x = L y = L =
100 points. The different pressure flow fields confirm thevery good isotropic dispersion properties of our approach, andthe pressure profiles further highlight its low numerical dissi-pation property.
C. Further validation: Shifting impact on non-linearphenomena
In the above section, it was shown that, in the linear regime,it is possible to extend the stability domain of our compress-ible DDF-LBM by adjusting the lattice to the local velocityfield. To further validate the adaptive reconstruction of miss-ing post-collision populations, we then propose to move to-wards more complex testcases, that include non-linearities,such as shock waves.In the following, the kinetic sensor (34) is used to reducethe generation of spurious Gibbs oscillations induced by dis-continuities, and which are more prominent in under-resolvedconditions and/or in the inviscid regime. All initial and bound-ary conditions are based on Grad’s reconstruction (27)-(31),where second-order FD approximations are used for the com-putation of σ αβ and q α .2 . x/L . y / L Reference u = 0 . u = 1 . u = 2 . − . − . − . . . . . .
20 0 . . . . . . x/L − . − . . . . . . ∆ P / ( P A ) u = 0 . u = 1 . u = 2 . FIG. 6: Convection of an acoustic pulse in an inviscid flow with a grid mesh composed of L =
100 points in each direction:global view (left) and slice along the x -axis (right). Normalized pressure fluctuations ∆ P / ( AP ) are plotted after 30 timeiterations for four different flow velocities u =
0, 1 .
3, 2 . .
9. Data corresponding to the latter three simulation has beenrecentered to ease to comparison with the reference state. Reference results correspond to u =
0, and are based on a meshcomposed of L =
100 and 500 points for the global view and the slice respectively.
1. Sod shock tube
To assess the ability of the proposed solver to deal withhighly compressible phenomena, we first consider a 1-D Rie-mann problem commonly referred to as Sod’s shock tube .In its most popular form, it consists of a 1-D simulation do-main initially divided into two sub-domains with differentdensities and temperatures, i.e., ( ρ L / ρ R , T L / T R ) = ( , . ) (39)with ( ρ R , T R ) = ( , T ) and u x = u y = . The rarefactionwave propagates towards the high-density region of the simu-lation domain and induces smooth variations of density, tem-perature and velocity fields. On the contrary, both the contactdiscontinuity and the shock wave propagate towards the low-density region and lead to discontinuous macroscopic fields,with the exception of the (normal) velocity and pressure fieldswhich remain constant at the contact discontinuity.Hereafter, we investigate the ability of our adaptive LBMto accurately reproduce the above features in the vanishingviscosity limit ( ν = [ L x × L y ] = [ × ] . More precisely, we focus our attentionon the impact of static and dynamic phase space transitions ondiscontinuities. Corresponding results are complied in Fig. 7for three different configurations: (1) no shift U x =
0, staticshift U x = x ∈ [ L / , L / ] , and (3) dynamic shift basedon our criterion (22).Interestingly, all waves are properly generated and correctlyevolve over time and space. The static configuration proves that discontinuities can travel through phase space transitionswithout generating high-amplitude spurious waves. This val-idates our reconstruction strategy (27)-(31). In addition, thedynamic configuration shows how the transition follows thevelocity field, as well as, its (almost) negligible impact on themacroscopic flow fields. It is also worth noting that bound-ary conditions based on Grad’s reconstruction (27)-(31) donot generate spurious oscillations and lead to accurate results.Eventually, only the contact discontinuity is overdissipated.This can be easily corrected by either increasing the meshresolution, or by fine tuning our sensor –as usually done forLBMs based on shock-capturing techniques .
2. Riemann 2D configuration
While the ability of our adaptive approach to cope withnon-linearities has been proven in the 1D case, hereafter, wefurther intend to validate it against a 2D Riemann problem.The latter corresponds to four 1D Riemann problems (or Sodshock tubes) that are initialized by dividing the simulation do-main into four quadrants. Depending on the density, pressureand velocity conditions imposed for each quadrant, a largepanel of compressible phenomena arise at quadrant interfaces,and especially in the vicinity of the simulation domain center.These phenomena were extensively studied in several sem-inal works (Schulz-Rinne et al. , Lax and Lui , Kurganovand Tadmor , etc). Based on the latter works, we focus ourattention on the configuration depicted in Table. I. This bench-mark test is of particular interest to validate adaptive transi-tions of phase space due to the complex interplay occurringbetween all waves.Generally speaking, several features/phenomena are ex-3 . . . . . . ρ / ρ R . . . . . . . . . . . . T / T R . . . . . . x/L x P / P R No shiftStaticDynamicReference 0 . . . . . . x/L x . . . . . . u x . . . . . . U x FIG. 7: Comparison of different fields (density ρ , temperature T , pressure P , velocity u x ) against exact solution of the Sodshock tube configuration. The spatial evolution of the shift ( U x ) is superimposed to the velocity field u x The two configurations‘No Shift’ and ‘Static’ corresponds to a fixed shift U x = ≤ x ≤ L , and U x = L / ≤ x ≤ L / n = . ρ [ ] = ρ [ ] = . u [ ] x = . u [ ] x = u [ ] y = u [ ] y = P [ ] = P [ ] = . ρ [ ] = . ρ [ ] = u [ ] x = u [ ] x = u [ ] y = u [ ] y = . P [ ] = P [ ] = TABLE I: Initialisation of each quadrant [ q ] of the simulationdomain ( q ∈ (cid:74) , (cid:75) ). This setup was proposed in previousworks, and it corresponds to configurations F , and12 respectively.pected when simulating this particular 2D Riemann problem.First, the initial state is symmetric with respect to the diago- nal axis ( x = y ), hence, it can be supposed that macroscopicfields will keep this symmetry property. Second, quadrants[2] and [4] show velocities ( u x and u y respectively) higherthan 0 .
5, whereas other quadrants ([1] and [3]) are initializedwith a flow at rest. Consequently, it is guaranteed that at leastthree versions of the D2Q21 lattice will coexist in the follow-ing simulation: ( U x , U y ) = ( , ) , ( , ) , and ( , ) . Third, allquadrants but [1] share the same pressure field, and have iden-tical perpendicular velocities at the quadrant interfaces, whichimplies that contact discontinuities are to be expected at theseinterfaces. Eventually, more complex calculations show thattwo shock waves form at interfaces with quadrant [1], and theypropagate towards the upper right corner of the simulation do-main. Doing so, they interact with each other, and generate acomplex pattern in this part of the simulation domain. Thelatter pattern further forces contact discontinuities to roll upinto a pair of vortices inside the third quadrant .Hereafter, the above Riemann 2D problem is simulated ina domain [ L x × L y ] using three grid mesh resolution: L x = L y = L = . y / L L = 200 L = 500 L = 1000 . . . . . .
60 0 . x/L . y / L . x/L . x/L (0 , , , , FIG. 8: 2D Riemann configuration: density ρ and velocity shift fields U x , U y (from top to bottom). Black plain lines are densityiso-contours: 30 equally-spaced levels in the density interval ρ ∈ [ .
54 1 . ] following the instructions provided bySchulz-Rinne et al. . The most common features of the flow match results obtained by previous works . For the sake ofclarity, only positive velocity shifts are plotted.in under- and well-resolved conditions, whereas the interme-diate case is typical of mesh resolution encountered in the lit-erature . The post-processing proposed in the latter workis also adopted hereafter, i.e., only half of the simulation do-main is plotted. Density ρ and velocity shifts ( U x , U y ) fieldsare compiled in Fig. 8 for the three grid mesh resolutions.Interestingly, the above-mentioned features are well recov-ered (symmetry, complex flow patterns), even if contact dis-continuities and the related ‘mushroom’-shaped dipole aregenerally over-dissipated –as it was already the case for theSod chock tube. Nevertheless, by increasing the mesh reso-lution, finer structures appear close to ( x , y ) = ( L / , L / ) andinside the complex pattern induced by the interplay betweenthe two shock waves. All of this favorably compares withprevious studies , even if high-frequency oscillationsare visible in the first quadrant. Since the present approachshows low numerical dissipation levels (as proved by resultspresented in Section IV B), this suggests that our stabilizationstrategy is too naive. It would then benefit from fine tuning(i) to better damp these high-frequency waves, as well as, (ii)to decrease the unnecessary amount of artificial viscosity thatimpacts the growth of the ‘mushroom’-shaped dipole.Regarding phase space transitions, it is interesting to note that they accurately follow the shock waves –and the complexpattern generated by the latter– in their propagation towardsthe upper right side of the simulation domain. Velocity flowconditions also generate a new phase space, which is basedon the shift ( , ) , upstream the intersection between the twoshock waves. Even though they are not plotted for the sakeof clarity, it is worth noting that three negative shifts also ap-pear inside the dipole when the mesh resolution is increased: ( U x , U y ) = ( − , ) , ( , − ) , and ( − , − ) .These simulation results confirm the viability of our ap-proach for interpolation-free adaptive LBMs in the context ofhigh-speed compressible flow simulations. V. CONCLUSION
Since the introduction of adaptive stencils in the LB con-text, no interpolation-free formulation has been proposed inthe literature. Knowing that space interpolations are com-putationally expensive, potentially prone to anisotropy issues(depending on the considered stencil for interpolation points),and deteriorate the parallel efficiency of the solver, it seemsappropriate to propose a solution to this issue. In this context,5the present work aims at proving the viability of interpolation-free adaptive LBMs for high-speed compressible flow simula-tions.The latter formulation leads to a series of new restrictionsand challenges to realize the transition between two phasespaces in a stable and accurate manner. In the context ofLBMs based on polynomial equilibrium distribution functions(EDFs), we investigated both points through linear stabilityanalyses of tensor-product based lattices. As a striking result,the overlap of stability domains is only possible for DdQq d lattices with q ≥ ξ i α ∈ { , ± , ± , ± , ± } ). Obvi-ously, this type of adaptive LBMs has only little interest inpractice. Alternatively, LBMs based on numerical EDFs aremore flexible due to their quadrature-free nature, and theyshow interesting stability properties for more compact lattices,and notably, the D2Q21 lattice.To couple the two shifted versions of the solver, one mustanswer the following questions: (1) when to switch betweenthe two solvers, and (2) how to exchange the information atthe interface between them. To tackle the former issue, weexplored a number of solutions based on either already ex-isting, or new switching criteria. In the end, it is proposedto locally adjust the lattice depending on velocity fluctua-tions, in accordance with the original formulation of adaptiveLBMs, and stability ranges of numerical EDFs. Regardingthe transfer of information at phase space transitions, miss-ing populations are recontructed through their equilibrium andnon-equilibrium contributions. While the EDF can be easilycomputed for both polynomial and numerical approaches, thenon-equilibrium part cannot be obtained in a straightforwardmanner for the latter case. Two (regularized) approaches areproposed to compute it, i.e., Chapman-Enskog’s and Grad’sformulations. The former relies on the space-time evolutionof EDFs, and requires their storage at previous time step(s),which is not suitable from the point of view of memory con-sumption. On the contrary, Grad’s approach is a generaliza-tion of already existing regularized approaches which hereassumes that non-equilibrium contributions are proportionalto EDFs and diffusive fluxes (viscous stress tensor and heatflux). Ultimately, this provides an efficient, on-the-fly recon-struction strategy that can be applied to both polynomial andnumerical EDFs.To account for internal degrees of freedom –mandatoryfor the simulation of polyatomic gas flows– the aboveinterpolation-free strategy is build on top of a double-distribution-function formulation (DDF-LBM). Following themethodology provided in , the D2Q21 lattice is coupled witha numerical EDF based on all convective constraints (8 mo-ments in 2D) in order to obtain a purely LB solver that offersa good trade-off between accuracy, efficiency and stability. Incase better stability is required, we further propose to locallyadjust the kinematic viscosity depending on the departure ofpopulations from their equilibrium state, which shares simi-larities with shock capturing techniques and flux limiters.The above DDF-LBM is validated in the low viscosityregime through several benchmark tests of increasing com-plexity. The first three tests aim at assessing the accuracyand stability of the present approach in the linear regime through the transport of shear, thermal and acoustic distur-bances. They confirm the excellent numerical properties (lowdispersion and dissipation) of our approach even for very highflow velocities ( || u || ≤ ACKNOWLEDGMENTS
The authors received no financial support for the research,authorship, and/or publication of this article.
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
Appendix A: Approximation of the Maxwellian in a givenreference frame
In the seminal work by Grad , velocity distribution func-tions are quantities that propagate at a temperature-normalizedpeculiar velocity ξ α = ( c α − u α ) / √ θ ⇐⇒ c α = ξ α √ θ + u α , (A1)which is the continuous counterpart of Eq. (5). Both formu-las correspond to a change of reference frame that leads toa Galilean and temperature invariant description of velocitydistribution functions (resp. populations). In other words, byworking with populations that propagate at a speed c α (resp. c i α ), macroscopic errors emerging from the approximation ofthe Maxwellian (e.g., Hermite polynomial expansion, numer-ical equilibrium, etc) boil down to zero, even for strong devi-ations from the reference state ( u , θ ) = ( , ) .In the discrete case, one can prove it by comparing the val-ues of moments computed via discrete EDFs with those of6 − − ξ x − . . . . . f MB ρ u x = 1 . θ = 1 . − − ξ x u x = 0 . θ = 2 . − − ξ x u x = 0 . θ = 2 . Maxwellian Hermite (Static) Hermite (Adapt) Hermite (Adapt Norm)
FIG. 9: Comparison of the Maxwellian with three Hermite polynomial expansion truncated at order 4. The static, adaptive andnormalized adaptive expansions are based on Hermite, central Hermite, and θ -normalized central Hermite momentsrespectively. Only the latter set of moments lead to an error-free approximation of the Maxwellian, whereas the adaptiveconfiguration shows θ -dependent errors. The static configuration suffers from both u - and θ -dependent errors.the Maxwellian, for different values of the macroscopic ve-locity u α and θ (see Section II C 2, Appendix D as well asRefs. ). Hereafter, we will focus on an elegant alter-native way to achieve the same goal in the continuous caseby simply comparing the shape of the Maxwellian with itsHermite polynomial expansion. The latter relies on the factthat moments are quantitative measures related to the shape ofa probability distribution function. Hence, distribution func-tions with a similar shape will have similar moments and viceversa.With this aim in mind, let us first recall some basics aboutthe polynomial expansion of the Maxwellian for a given ref-erence frame. By definition, the Maxwell-Boltzmann EDFreads as f MB ( x , ξ , t ) = ρ ( πθ ) d / exp (cid:20) − ( ξ − u ) θ (cid:21) , (A2)and its Hermite polynomial expansion (truncated at order N )is expressed as f MB , N ( x , ξ , t ) = ρθ d / w ( η ) N ∑ n = n ! a ( n ) MB ( x , t ) : H ( n ) ( η ) . (A3) H ( n ) are Hermite polynomials at order n computed via η α ( α = x or y in 2D) which is a quantity yet to be defined. a ( n ) MB are Hermite coefficients normalized by the density ρ , i.e., a ( n ) MB ( x , t ) = ρ (cid:90) f MB ( x , ξ , t ) H ( n ) ( η ) d η . (A4)In Eq. (A3), “:” further stands for Frobenius inner product alsoknown as full tensor contraction, and w is the weight function w ( η ) = ( π ) d / exp (cid:0) − η / (cid:1) (A5) used for the Gauss-Hermite quadrature. Several definitionsof η α exist depending on the reference frame considered .For a reference frame at rest, η α = ξ α and a ( n ) MB is a Her-mite moment of the Maxwellian, whereas in the co-movingreference frame, η α = ξ α − u α and a ( n ) MB is now a central Her-mite moment . In the most general case, the weightfunction further accounts for temperature variations through η α = ( ξ α − u α ) / √ θ , hence leading to w = θ d / ρ f MB . (A6)Here, polynomial coefficients are central Hermite momentsnormalized by the temperature, and they satisfy the ele-gant mathematical property of all being null but the zeroth-order term – central Hermite moments only share thesame property in the isothermal case. Consequently, a ( n ) MB are velocity- and temperature-independent, whereas for η α = ξ α − u α , they will depend on temperature variations. For areference frame at rest, these moments will finally depend onboth velocity and temperature variations.This behavior is illustrated in Fig. 9, where the Maxwellianand its three Hermite polynomial expansions are plotted inthe 1D case. These polynomial expansions are done based onHermite (static), central Hermite (adapt), and θ -normalizedcentral Hermite moments (adapt norm). All of them includeterms up to N =
4, which is the minimal configuration to re-cover the macroscopic behavior of the compressible Navier-Stokes-Fourier equations . As expected, the static con-figuration shows errors when the operating point is far fromthe reference state ( u , θ ) = ( , ) . In addition, describing theMaxwellian in terms of central Hermite moments prevents anyvelocity-dependent mismatch, but it cannot correct errors in-duced by temperature variation. In the end, accounting for7velocity and temperature variations in the propagation speedis the only way to obtain an error-free approximation of theMaxwellian, and consequently, the correct macroscopic be-havior. Appendix B: Linear stability analysis in a nutshell
For non-linear systems of equations, this approach consistsof introducing a perturbation into the linearized version of theconsidered system of equations in a fully periodic domain,then following its time evolution . In the context of LBsolvers, under the assumption of a linear regime, one can ap-proximate the field though a first-order Taylor-McLaurin ex-pansion: f i ≈ f i + f (cid:48) i . (B1)Introducing this expansion into the target discrete system ofequations one recovers the discrete linearized equations forthe perturbation.Given that the topic has been thoroughly treated in the liter-ature , details and derivation of the finalequations will be omitted. As an example, the linearization ofthe BGK-LBE leads to f (cid:48) i ( x + ξ i , t + ) = (cid:20) δ i j (cid:18) − τ (cid:19) + τ J eqi j (cid:21) f (cid:48) j ( x , t ) . (B2)where J eqi j = ∂ f eqi / ∂ f j is the Jacobian matrix of the EDF thatis evaluated at f j = f j . For a detailed derivation of these Jaco-bians, interested readers are referred to .Introducing the standing waves (wave number k ∈ R and timefrequency ω ∈ C ) into the linearized discrete time evolutionequation f (cid:48) j = F (cid:48) j exp ( i ω t − k · x ) , (B3)one obtains a system of equations, through which the time-amplification factor of the perturbation [exp [ i ω ] ) can be ob-tained for different wave numbers k . In order for the sys-tem to be linearly stable for the chosen set of parameters, i.e.velocity, temperature and non-dimensional viscosity, the realcomponent of exp [ i ω ] must remain negative for all possiblevalues of k .For all stability domains presented in Section II B, the 1Dwave number space k ∈ [ , π ] is discretized using a resolu-tion of 100 points , i.e. ∆ k x = π / . Appendix C: Detailed description of numerical equilibria1. Exact formulation for polyatomic gases
For LBMs to fully recover the macroscopic behavior ofNSF equations, the discrete EDF must mimick up to the (traceof the) fourth-order moment of the Maxwellian. In the bidi-mensional case, it should then satisfy the following 13 con-straints: G = ∑ i f eqi − ρ , (C1) G , x = ∑ i f eqi ξ ix − ρ u x , (C2) G , y = ∑ i f eqi ξ iy − ρ u y , (C3) G , xx = ∑ i f eqi ξ ix − ρ ( u x + T ) , (C4) G , xy = ∑ i f eqi ξ ix ξ iy − ρ u x u y , (C5) G , yy = ∑ i f eqi ξ iy − ρ ( u y + T ) , (C6) G , xxx = ∑ i f eqi ξ ix − ρ u x ( u x + T ) , (C7) G , xxy = ∑ i f eqi ξ ix ξ iy − ρ ( u x + T ) u y , (C8) G , xyy = ∑ i f eqi ξ ix ξ iy − ρ u x ( u y + T ) , (C9) G , yyy = ∑ i f eqi ξ iy − ρ u y ( u y + T ) , (C10) G , xx αα = ∑ i f eqi ξ ix ξ i α − ρ [( E + T ) u x + ( E + T ) T ] , (C11) G , xy αα = ∑ i f eqi ξ ix ξ iy ξ i α − ρ [( E + T ) u x u y ] , (C12) G , yy αα = ∑ i f eqi ξ iy ξ i α − ρ [( E + T ) u y + ( E + T ) T ] , (C13)where the repetition of (Greek) indices stands for Einstein’ssummation rule on geometrical coordinates, e.g., ξ i α = ξ ix + ξ iy .The above constraints (C1)-(C13) are then used for thecomputation of the numerical equilibrium f eqi = ρ exp (cid:2) − (cid:0) + λ + (cid:8) λ , x ξ ix + λ , y ξ iy (cid:9) + (cid:8) λ , xx ξ ix + λ , xy ξ ix ξ iy + λ , yy ξ iy (cid:9) + (cid:8) λ , xxy ξ ix ξ iy + λ , xyy ξ ix ξ iy + λ , yyy ξ iy (cid:9) + (cid:8) λ , xx ξ ix ξ i α + λ , xy ξ ix ξ iy ξ i α + λ , yy ξ iy ξ i α (cid:9)(cid:1)(cid:3) , (C14)where λ n are Lagrange multipliers that are obtained from a root-finding algorithm as zeros of the constraints G n .8In addition, it is important to understand that populations f i do not account for extra internal degrees of freedom (ro-tational and vibrational). Hence, the total energy used in theabove constraints is defined as 2 E = u x + u y + DT , and D = g i is used to impose the correct specific heat ratio γ r . While populations f i transfer the monatomic informationto g i through g eqi (33), the feedback from g i to f i is implicit(i.e., no forcing term is used in the collision of f i ), and it oc-curs in the computation of the temperature T = ρ C v (cid:20) ∑ i ( ξ i α f i + g i ) − ρ u α (cid:21) (C15)with C v = / ( γ r − ) being the polyatomic heat capacity atconstant volume. This “polyatomic” temperature is then in-jected in the constraints (C1)-(C13) that are used for the com-putation of f eqi (C14), hence, closing the loop.The latter methodology allows the user to compute f eqi ina numerical manner, i.e., even if the system (C1)-(C13) doesnot have an analytical solution. The resulting LBM is thenfreed from the very constraining quadrature rules that imposea minimal lattice size to recover a macroscopic behavior of in-terest. Nevertheless, the user must cautiously choose the con-vergence criterion of the root-finding algorithm. If not, con-servation issues will appear because the constraints will not becorrectly imposed. In the end, it is sufficient to impose a con-vergence criterion of 10 − , which leads to constraint errorsthat oscillates around 10 − , i.e., close to machine precision.
2. Reduced models
For both DVMs and LBMs, numerical equilibria were firstbased on constraints corresponding to the conservation ofmass, momentum, and energy. In the bidimensional case, thiscorresponds to the following set of four constraints G = ∑ i f eqi − ρ , (C16) G , x = ∑ i f eqi ξ ix − ρ u x , (C17) G , y = ∑ i f eqi ξ iy − ρ u y , (C18) G , αα = ∑ i f eqi ξ i α − ρ E . (C19)The latter are then used to compute the numerical equilibrium f eqi = ρ exp (cid:2) − (cid:0) + λ + (cid:8) λ , x ξ ix + λ , y ξ iy (cid:9) + λ , αα ξ i α (cid:1)(cid:3) . (C20)Even if this methodology reduces the CPU cost –induced bythe iterative computation of the equilibrium– as compared tothe more demanding 13-moment approach (18), it also comesat the expense of part of physics since both convective and dif-fusive terms are not properly recovered (See Eq. 13). Usually,one can only recover the correct physics by increasing the sizeof the lattice, as pointed out in the parametric study conductedin Section 3(c) of our previous work .In the context of high-speed and high-Reynolds numberflows, (numerical) errors related to diffusive phenomena havea lower impact on the accuracy of the solver than those re-lated to convective phenomena. This is the reason why, NSsolvers are usually based on second-order (centered) numeri-cal schemes for diffusive fluxes whereas convective terms arediscretized using higher-order schemes . Applying thelatter reasoning to numerical EDFs, one can improve the 4-moment methodology by accounting for all constraints relatedto convective fluxes –so that ∆ = ∆ tr = . The corresponding set of eight constraints reads as G = ∑ i f eqi − ρ , (C21) G , x = ∑ i f eqi ξ ix − ρ u x , (C22) G , y = ∑ i f eqi ξ iy − ρ u y , (C23) G , xx = ∑ i f eqi ξ ix − ρ ( u x + T ) , (C24) G , xy = ∑ i f eqi ξ ix ξ iy − ρ u x u y , (C25) G , yy = ∑ i f eqi ξ iy − ρ ( u y + T ) , (C26) G , x αα = ∑ i f eqi ξ ix ξ i α − ρ u x ( E + T ) , (C27) G , y αα = ∑ i f eqi ξ iy ξ i α − ρ u y ( E + T ) , (C28)and it is used to compute the equilibrium f eqi = ρ exp (cid:2) − (cid:0) + λ + (cid:8) λ , x ξ ix + λ , y ξ iy (cid:9) + (cid:8) λ , xx ξ ix + λ , xy ξ ix ξ iy + λ , yy ξ iy (cid:9) + (cid:8) λ , x αα ξ ix ξ i α + λ , y αα ξ iy ξ i α (cid:9)(cid:1)(cid:3) . (C29)Several alternative formulations can be derived depending on the targeted physics, i.e., which error terms ∆ n and/or ∆ trn shouldbe cancelled in Eq. (13). As an example, if diffusive fluxes are negligible in the energy equation but not in the momentumequation, one could rely on the set of ten constraints (C1)-(C10) so that only ∆ tr would be non-zero. One could further wantto better fit the definition of “admissible” spaces which requires that the maximal order of constraints should always be even .This is done adding the constraint G , ααββ = G , xx ββ + G , yy ββ = ∑ i f eqi ξ i α ξ i β − ρ [( E + T )( u x + u y ) + ( E + T ) T ] , (C30)to the set of either eight or ten constraints leading to 9- and 11-moment approaches [see Eqs. (15) and (17) respectively]. Appendix D: Preliminary study on the choice of lattice and number of constraints
To find a good trade off between accuracy, stability and efficiency, we followed the methodology introduced in our previouspaper . The latter consists in evaluating the stability domain and accuracy of numerical EDFs for several lattices. This is doneby computing macroscopic deviations of the numerical EDF moments with respect to the Mach number, and at a given referencetemperature.This aspect was investigated in a preliminary study for lattices composed of 9, 13, 17, 21, 25 (zero-one-three or “ZOT” formula-tion) and 49 discrete velocities respectively (see Fig. 10 and Tab. II for their reference temperature T ). Deviations with respectto the Maxwellian moments are reported for a numerical EDF based on M = ε =
10% threshold, especially with the8-moment approach. In summary, the D2Q21 lattice is less accurate than the D2Q49 lattice, but it is far more efficient both interms of simulation time and memory storage. In addition, it shares similar accuracy properties with the D2Q13 lattice whilebeing more stable (the root-finding solver converges up to Ma ≈ .
85 for the D2Q21 lattice). For all these reasons, coupling theD2Q21 lattice with the 8-moment approach seems to be a good trade off in terms of accuracy, stability and efficiency.
FIG. 10: Velocity discretizations considered in the preliminary study. From top left to bottom right: D2Q9, D2Q13, D2Q17,D2Q21, D2Q25ZOT and D2Q49. Data are compiled from the works and therein references.
Lattice Q9 Q13 Q17 Q21 Q25ZOT Q49 T / . . . . . TABLE II: Reference temperatures T obtained following instructions provided in our previous work .0 Ma10 − − − − − ε Q9 Q13 M xxx M xx M xyy M yy M xxαα M xαα Q17 − − − − − ε Q21
Q25ZOT
Q49
Ma10 − − − − − ε Q9 Q13 Q17 . . . . . Ma − − − − − ε Q21 . . . . . Ma Q25ZOT . . . . . Ma Q49
FIG. 11: Mach number impact on macroscopic errors ε for the 4-moment (top two rows) and 8-moment (bottom two rows)approaches. The flow propagates along the x -axis, and the grey zone starts at ε = − (i.e., machine precision). For all configurations, T = T with the reference temperature T defined in Tab. II for each lattice.1 Appendix E: Flexible Prandtl number extension
By relying on a single relaxation time approximation (32a), the present DDF-LBM is restricted to Pr =
1. This can be correctedby adopting more sophisticated collisions models, e.g., the ES-BGK or Shakov’s collision model. Both were introduced inthe monatomic case, and extended to polyatomic gases by Andries et al. and Rykov respectively. In the following, we adoptRykov’s approach where we further suppose that collisions are elastic. This basically modifies the collision step (32a) as follows: f ∗ i = f eqi + (cid:18) − τ f (cid:19) ( f i − f eqi ) + (cid:18) τ f − τ Pr (cid:19) q α ξ i α ρ C p T (cid:18) ξ i T − C p (cid:19) f eqi , (E1) g ∗ i = g eqi + (cid:18) − τ g (cid:19) ( g i − g eqi ) + (cid:18) τ g − τ Pr (cid:19) q α ξ i α ρ C p g eqi , (E2)with τ f = τ g = . + ν / T , and τ Pr = . + ( ν / Pr ) / T .In order to extend the reconstruction of missing information at transition interfaces, one simply needs to replace the non-equilibrium contributions h neqi = h i − h eqi by h ( ) , CE i or h ( ) , Grad i with h i = f i , g i . Interestingly, when the reconstruction step isapplied to all cells of the simulation domain (at each time step), one ends up with an extended regularized collision model, whoseproperties will be investigated in a future paper. P. Bhatnagar, E. Gross, and M. Krook, “A model for collision processes ingases. I. Small amplitude processes in charged and neutral one-componentsystems,” Phys. Rev. , 511–525 (1954). L. H. Holway, “New statistical models for kinetic theory: Methods of con-struction,” Phys. Fluids , 1658–1673 (1966). E. M. Shakhov, “Generalization of the Krook kinetic relaxation equation,”Fluid Dyn. , 95–96 (1968). P. Andries, P. Le Tallec, J.-P. Perlat, and B. Perthame, “The Gaussian-BGK model of Boltzmann equation with small Prandtl number,” Eur. J.Mech. - B/Fluids , 813 – 830 (2000). V. A. Rykov, “A model kinetic equation for a gas with rotational degreesof freedom,” Fluid Dyn. , 959–966 (1975). L. Mieussens, “A survey of deterministic solvers for rarefied flows (In-vited),” AIP Conf. Proc. , 943–951 (2014). L. Mieussens, “Discrete velocity model and implicit scheme for the BGKequation equation of rarefied gas dynamics,” Math. Models Methods Appl.Sci. , 1121–1149 (2000). P. Le Tallec and J.-P. Perlat, “Numerical Analysis of Levermore’s MomentSystem,” Research Report RR-3124 (INRIA, 1997) projet M3N. P. Charrier, B. Dubroca, J.-L. Feugeas, and L. Mieussens, “Modèles àvitesses discrètes pour le calcul d’écoulements hors équilibre cinétique,”C. R. Acad. Sci. Paris, Ser. I , 1347 – 1352 (1998). V. A. Titarev, “Efficient deterministic modelling of three-dimensional rar-efied gas flows,” Commun. Comput. Phys. , 162–192 (2012). C. Baranger, J. Claudel, N. Hérouard, and L. Mieussens, “Locally refineddiscrete velocity grids for stationary rarefied flow simulations,” J. Comput.Phys. , 572 – 593 (2014). V. Aristov, O. Ilyin, and O. Rogozin, “Kinetic multiscale scheme basedon the discrete-velocity and lattice-Boltzmann methods,” J. Comput. Sci. , 101064 (2020). F. Schornbaum and U. Rüde, “Massively parallel algorithms for the latticeBoltzmann method on nonuniform grids,” SIAM J. Sci. Comput , C96–C126 (2016). X. Shan and X. He, “Discretization of the velocity space in the solution ofthe Boltzmann equation,” Phys. Rev. Lett. , 65–68 (1998). X. Shan, X.-F. Yuan, and H. Chen, “Kinetic theory representation of hy-drodynamics: A way beyond the Navier–Stokes equation,” J. Fluid Mech. , 413–441 (2006). Z. Guo and C. Shu,
Lattice Boltzmann Method and Its Applications inEngineering (World Scientific, 2013). E. Manoha and B. Caruelle, “Summary of the LAGOON solutions fromthe benchmark problems for airframe noise computations-III workshop,”in (2015) p. 2846. M. F. Barad, J. G. Kocheemoolayil, and C. C. Kiris, “Lattice Boltzmannand Navier-Stokes Cartesian CFD approaches for airframe noise predic-tions,” in . C. L. Rumsey, J. P. Slotnick, and A. J. Sclafani, “Overview and summaryof the third aiaa high lift prediction workshop,” J. Aircraft , 621–644(2018). Y. Hou, D. Angland, A. Sengissen, and A. Scotto, “Lattice-Boltzmann andnavier-stokes simulations of the partially dressed, cavity-closed nose land-ing gear benchmark case,” in (2019) p. 2555. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M.Viggen,
The Lattice Boltzmann Method: Principles and Practice (SpringerInternational Publishing, 2017). S. Succi,
The Lattice Boltzmann Equation: For Complex States of FlowingMatter (Oxford University Press, 2018). M. Bauer, S. Eibl, C. Godenschwager, N. Kohl, M. Kuron, C. Ret-tinger, F. Schornbaum, C. Schwarzmeier, D. Thönnes, H. Köstler,and U. Rüde, “waLBerla: A block-structured high-performance frame-work for multiphysics simulations,” Comput. Math. Appl. (2020),https://doi.org/10.1016/j.camwa.2020.01.007. M. J. Krause, A. Kummerländer, S. J. Avis, H. Kusumaatmaja, D. Dapelo,F. Klemens, M. Gaedtke, N. Hafen, A. Mink, R. Trunk, J. E. Mar-quardt, M.-L. Maier, M. Haussmann, and S. Simonis, “OpenLB -Open source lattice Boltzmann code,” Comput. Math. Appl. (2020),https://doi.org/10.1016/j.camwa.2020.04.033. J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. Lagrava, F. Brogi,M. B. Belgacem, Y. Thorimbert, S. Leclaire, S. Li, F. Marson, J. Lemus,C. Kotsalos, R. Conradin, C. Coreixas, R. Petkantchin, F. Raynaud,J. Beny, and B. Chopard, “Palabos: Parallel lattice Boltzmann solver,”Comput. Math. Appl. (2020), 10.1016/j.camwa.2020.03.022. A. De Rosis, R. Huang, and C. Coreixas, “Universal formulation ofcentral-moments-based lattice Boltzmann method with external forcingfor the simulation of multiphysics phenomena,” Phys. Fluids , 117102(2019). A. De Rosis and C. Coreixas, “Multiphysics flow simulations usingD3Q19 lattice Boltzmann methods based on central moments,” arXivpreprint arXiv:2010.01628 (2020). C. Lin, A. Xu, G. Zhang, and Y. Li, “Double-distribution-function dis-crete Boltzmann model for combustion,” Combust. Flame , 137 – 151(2016). C. Lin and K. H. Luo, “Discrete Boltzmann modeling of unsteady reactiveflows with nonequilibrium effects,” Phys. Rev. E , 012142 (2019). C. Lin, K. H. Luo, A. Xu, Y. Gan, and H. Lai, “MRT discrete Boltzmannmodeling of multicomponent mixture with nonequilibrium effects,” arXivpreprint arXiv:2002.02668 (2020), arXiv:2002.02668. Z. Guo, K. Xu, and R. Wang, “Discrete unified gas kinetic scheme forall Knudsen number flows: Low-speed isothermal case,” Phys. Rev. E ,033305 (2013). Z. Guo, R. Wang, and K. Xu, “Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case,” Phys. Rev. E ,033313 (2015). B. Nadiga and D. Pullin, “A method for near-equilibrium discrete-velocitygas flows,” J. Comput. Phys. , 162 – 172 (1994). B. Nadiga, “An adaptive discrete-velocity model for the shallow waterequations,” J. Comput. Phys. , 271 – 280 (1995). B. Nadiga, “An euler solver based on locally adaptive discrete velocities,”J. Stat. Phys. , 129–146 (1995). J. Huang, F. Xu, M. Vallières, D. H. Feng, Y.-H. Qian, B. Fryxell, andM. R. Strayer, “A thermal LBGK model for large density and temperaturedifferences,” Int. J. Mod. Phys. C , 827–841 (1997). C. Sun, “Lattice-Boltzmann models for high speed flows,” Phys. Rev. E , 7283–7287 (1998). C. Sun, “Adaptive lattice Boltzmann model for compressible flows: Vis-cous and conductive properties,” Phys. Rev. E , 2645–2653 (2000). C. Sun, “Simulations of compressible flows with strong shocks by an adap-tive lattice Boltzmann model,” J. Comput. Phys. , 70 – 84 (2000). C. Sun and A. T. Hsu, “Three-dimensional lattice Boltzmann model forcompressible flows,” Phys. Rev. E , 016303 (2003). C. Sun and A. T. Hsu, “Multi-level lattice Boltzmann model on squarelattice for compressible flows,” Comput. Fluids , 1363 – 1385 (2004). R. Zhang, X. Shan, and H. Chen, “Efficient kinetic method for fluidsimulation beyond the Navier-Stokes equation,” Phys. Rev. E , 046703(2006). S. A. Hosseini, C. Coreixas, N. Darabiha, and D. Thévenin, “Extensiveanalysis of the lattice Boltzmann method on shifted stencils,” Phys. Rev. E , 063301 (2019). I. Karlin and P. Asinari, “Factorization symmetry in the lattice Boltzmannmethod,” Physica A , 1530 – 1548 (2010). J. Von Neumann and A. W. Burks,
Theory of self-reproducing automata (University of Illinois Press Urbana, 1996). D. A. Wolf-Gladrow,
Lattice-gas cellular automata and lattice Boltzmannmodels: an introduction (Springer, 2004). C. Hirsch,
Numerical computation of internal and external flows: The fun-damentals of Computational Fluid Dynamics (Elsevier, 2007). D. Wilde, A. Krämer, K. Küllmer, H. Foysi, and D. Reith, “Multisteplattice Boltzmann methods: Theory and applications,” Int. J. Numer. Meth.Fluids , 156–169 (2019). D. Ricot, S. Marié, P. Sagaut, and C. Bailly, “Lattice Boltzmann methodwith selective viscosity filter,” J. Comput. Phys. , 4478 – 4490 (2009). S. Marié, D. Ricot, and P. Sagaut, “Comparison between lattice Boltz-mann method and Navier–Stokes high order schemes for computationalaeroacoustics,” J. Comput. Phys. , 1056 – 1070 (2009). P. Lallemand and L.-S. Luo, “Theory of the lattice Boltzmann method:Dispersion, dissipation, isotropy, galilean invariance, and stability,” Phys.Rev. E , 6546–6562 (2000). H. Xu and P. Sagaut, “Optimal low-dispersion low-dissipation LBMschemes for computational aeroacoustics,” J. Comput. Phys. , 5353 –5382 (2011). H. Xu, O. Malaspinas, and P. Sagaut, “Sensitivity analysis and determina-tion of free relaxation parameters for the weakly-compressible MRT-LBMschemes,” J. Comput. Phys. , 7335 – 7367 (2012). M. Chávez-Modena, E. Ferrer, and G. Rubio, “Improving the stabilityof multiple-relaxation lattice Boltzmann methods with central moments,”Comput. Fluids , 397 – 409 (2018). S. A. Hosseini, C. Coreixas, N. Darabiha, and D. Thévenin, “Stabilityof the lattice kinetic scheme and choice of the free relaxation parameter,”Phys. Rev. E , 063305 (2019). F. Dubois, T. Février, and B. Graille, “On the stability of a relative velocitylattice Boltzmann scheme for compressible Navier–Stokes equations,” C.R. Mécanique , 599 – 610 (2015). C. Coreixas, G. Wissocq, B. Chopard, and J. Latt, “Impact of collisionmodels on the physical properties and the stability of lattice boltzmannmethods,” Phil. Trans. R. Soc. A , 20190397 (2020). G. Wissocq, C. Coreixas, and J. Boussuge, “Linear stability of athermalregularized lattice Boltzmann methods,” arXiv preprint arXiv:2006.07353(2020), arXiv:2006.07353 [physics.comp-ph]. S. A. Hosseini, N. Darabiha, and D. Thévenin, “Compressibility in lat-tice boltzmann on standard stencils: effects of deviation from referencetemperature,” Phil. Trans. R. Soc. A , 20190399 (2020). B. Dubroca and J.-L. Feugeas, “Etude théorique et numérique d’une hiérar-chie de modèles aux moments pour le transfert radiatif,” C. R. Acad. Sci.Paris, Ser. I , 915 – 920 (1999). L. Mieussens, “Discrete-velocity models and numerical schemes for theboltzmann-bgk equation in plane and axisymmetric geometries,” J. Com-put. Phys. , 429 – 466 (2000). L. Mieussens and H. Struchtrup, “Numerical comparison of Bhatna-gar–Gross–Krook models with proper Prandtl number,” Phys. Fluids ,2797–2813 (2004). H. C. Öttinger, H. Struchtrup, and M. Torrilhon, “Formulation of momentequations for rarefied gases within two frameworks of non-equilibriumthermodynamics: RET and GENERIC,” Philos. Trans. R. Soc. A ,20190174 (2020). J. Latt, C. Coreixas, J. Beny, and A. Parmigiani, “Efficient supersonic flowsimulations using lattice Boltzmann methods based on numerical equilib-ria,” Phil. Trans. R. Soc. A (2020), 10.1098/rsta.2019.0559. S. Chapman and T. Cowling,
The Mathematical Theory of Non-uniformGases: An Account of the Kinetic Theory of Viscosity, Thermal Conductionand Diffusion in Gases (Cambridge University Press, 1970). C. D. Levermore, “Moment closure hierarchies for kinetic theories,” J.Stat. Phys. , 1021–1065 (1996). P. J. Kornreich and J. Scalo, “Supersonic lattice gases: Restoration ofGalilean invariance by nonlinear resonance effects,” Physica D , 333– 344 (1993). GSL Library 2.6, “Multidimensional root-finding,” (2020). P. Rabinowitz,
Numerical methods for nonlinear algebraic equations , 1sted. (Routledge, 1970). B. Dubroca and L. Mieussens, “A conservative and entropic discrete-velocity model for rarefied polyatomic gases,” ESAIM: Proceedings ,127–139 (2001). V. Titarev, “Conservative numerical methods for advanced model kineticequations,” in
ECCOMAS CFD 2006: Proceedings of the European Con-ference on Computational Fluid Dynamics, Egmond aan Zee, The Nether-lands, September 5-8, 2006 (2006). X. Nie, X. Shan, and H. Chen, “Thermal lattice Boltzmann model forgases with internal degrees of freedom,” Phys. Rev. E , 035701 (2008). N. Frapolli,
Entropic lattice Boltzmann models for thermal and compress-ible flows , Ph.D. thesis, ETH-Zürich (2017). G. Di Staso, H. Clercx, S. Succi, and F. Toschi, “Dsmc–lbm mappingscheme for rarefied and non-rarefied gas flows,” J. Comput. Sci. , 357 –369 (2016), discrete Simulation of Fluid Dynamics 2015. G. Di Staso, S. Srivastava, E. Arlemark, H. Clercx, and F. Toschi, “Hy-brid lattice Boltzmann-direct simulation monte carlo approach for flows inthree-dimensional geometries,” Comput. Fluids , 492 – 509 (2018). V. V. Aristov, O. V. Ilyin, and O. A. Rogozin, “Fluid-kinetic coupling ofthe BGK and lattice Boltzmann equations,” in (2017). V. V. Aristov, O. V. Ilyin, and O. A. Rogozin, “A hybrid numerical schemebased on coupling discrete-velocities models for the BGK and LBGKequations,” AIP Conf. Proc. , 060007 (2019). J. Meng, Y. Zhang, and X. Shan, “Multiscale lattice Boltzmann approachto modeling gas flows,” Phys. Rev. E , 046701 (2011). A. Dupuis and B. Chopard, “Theory and applications of an alternativelattice Boltzmann grid refinement algorithm,” Phys. Rev. E , 066707(2003). D. Lagrava, O. Malaspinas, J. Latt, and B. Chopard, “Advances in multi-domain lattice Boltzmann grid refinement,” J. Comput. Phys. , 4808 –4822 (2012). F. Gendre, D. Ricot, G. Fritz, and P. Sagaut, “Grid refinement for aeroa-coustics in the lattice Boltzmann method: A directional splitting ap-proach,” Phys. Rev. E , 023311 (2017). T. Astoul, G. Wissocq, J. françois Boussuge, A. Sengissen, andP. Sagaut, “Lattice boltzmann method for computational aeroacoustics onnon-uniform meshes: a direct grid coupling approach,” arXiv preprintarXiv:2004.14887 (2020), arXiv:2004.14887 [physics.comp-ph]. I. D. Boyd, G. Chen, and G. V. Candler, “Predicting failure of the con-tinuum fluid equations in transitional hypersonic flows,” Phys. Fluids ,210–219 (1995). S. Tiwari, “Coupling of the Boltzmann and Euler equations with automaticdomain decomposition,” J. Comput. Phys. , 710 – 726 (1998). W.-L. Wang and I. D. Boyd, “Predicting continuum breakdown in hyper-sonic viscous flows,” Physics of Fluids , 91–100 (2003). D. A. Lockerby, J. M. Reese, and H. Struchtrup, “Switching criteria forhybrid rarefied gas flow solvers,” Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences , 1581–1598 (2009). J. Meng, N. Dongari, J. M. Reese, and Y. Zhang, “Breakdown parameterfor kinetic modeling of multiscale gas flows,” Phys. Rev. E , 063305(2014). B. Dorschner, F. Bösch, and I. V. Karlin, “Particles on demand for kinetictheory,” Phys. Rev. Lett. , 130602 (2018). A. V. Zakirov, B. A. Korneev, V. D. Levchenko, and A. Y. Perepelkina,“On the conservativity of the particles-on-demand method for solution ofthe discrete Boltzmann equation,” Keldysh Inst. Prepr. , 1–19 (2019). J. Latt and B. Chopard, “Lattice Boltzmann method with regularizedpre-collision distribution functions,” Math. Comput. Simul. , 165–168(2006). H. Chen, R. Zhang, I. Staroselsky, and M. Jhon, “Recovery of full rota-tional invariance in lattice Boltzmann formulations for high Knudsen num-ber flows,” Physica A , 125–131 (2006). A. Z. E. Zipunova, A. Perepelkina, “Regularization and the Particles-on-Demand method for the solution of the discrete Boltzmann equation,”(2020), presented at DSFD 29. A. J. C. Ladd and R. Verberg, “Lattice-Boltzmann simulations of particle-fluid suspensions,” J. Stat. Phys. , 1191–1251 (2001). O. Malaspinas, “Increasing stability and accuracy of the lattice Boltzmannscheme: Recursivity and regularization,” arXiv preprint arXiv:1505.06900(2015). C. Coreixas, G. Wissocq, G. Puigt, J.-F. Boussuge, and P. Sagaut, “Recur-sive regularization step for high-order lattice Boltzmann methods,” Phys.Rev. E , 033306 (2017). K. K. Mattila, P. C. Philippi, and L. A. Hegele Jr., “High-order regulariza-tion in lattice-Boltzmann equations,” Phys. Fluids , 046103 (2017). H. Grad, “On the kinetic theory of rarefied gases,” Commun. Pure Appl.Math. , 331–407 (1949). H. Grad, “Statistical mechanics, thermodynamics, and fluid dynamics ofsystems with an arbitrary number of integrals,” Commun. Pure Appl.Math. , 455–494 (1952). A. Kogan, “Derivation of Grad’s type equations and study of their relax-ation properties by the method of maximization of entropy,” J. Appl. Math.Mech. , 130 – 142 (1965). P. A. Skordos, “Initial and boundary conditions for the lattice Boltzmannmethod,” Phys. Rev. E , 4823–4842 (1993). J. Latt, B. Chopard, O. Malaspinas, M. Deville, and A. Michler, “Straightvelocity boundaries in the lattice Boltzmann method,” Phys. Rev. E ,056703 (2008). B. Dorschner, S. Chikatamarla, F. Bösch, and I. Karlin, “Grad’s approxi-mation for moving and stationary walls in entropic lattice Boltzmann sim-ulations,” J. Comput. Phys. , 340 – 354 (2015).
E. Garnier, N. Adams, and P. Sagaut,
Large eddy simulation for compress-ible flows (Springer Science & Business Media, 2009).
S. Pirozzoli, “Numerical methods for high-speed flows,” Annu. Rev. FluidMech , 163–194 (2011). A. Gorban and D. Packwood, “Enhancement of the stability of latticeBoltzmann methods by dissipation control,” Physica A , 285 – 299(2014).
X. Shan, “Central-moment-based galilean-invariant multiple-relaxation- time collision model,” Phys. Rev. E , 043308 (2019).
F. Renard, Y. Feng, J.-F. Boussuge, and P. Sagaut, “Improved com-pressible hybrid lattice Boltzmann method on standard lattice for sub-sonic and supersonic flows,” arXiv preprint arXiv:2002.03644 (2020),arXiv:2002.03644 [physics.flu-dyn].
D. Fabre, L. Jacquin, and J. Sesterhenn, “Linear interaction of a cylindricalentropy spot with a shock,” Phys. Fluids , 2403–2422 (2001). G. Farag, S. Zhao, T. Coratger, P. Boivin, G. Chiavassa, and P. Sagaut, “Apressure-based regularized lattice-boltzmann method for the simulation ofcompressible flows,” Phys. Fluids , 066106 (2020). C. W. Schulz-Rinne, J. P. Collins, and H. M. Glaz, “Numerical solution ofthe Riemann problem for two-dimensional gas dynamics,” SIAM Journalon Scientific Computing , 1394–1414 (1993). P. D. Lax and X.-D. Liu, “Solution of two-dimensional Riemann problemsof gas dynamics by positive schemes,” SIAM J. Sci. Comput , 319–340(1998). A. Kurganov and E. Tadmor, “Solution of two-dimensional Riemann prob-lems for gas dynamics without Riemann problem solvers,” Numer. Meth-ods Partial Differ. Equ. , 584–608 (2002). G. A. Sod, “A survey of several finite difference methods for systemsof nonlinear hyperbolic conservation laws,” Journal of ComputationalPhysics , 1–31 (1978). C. Coreixas,
High-order extension of the recursive regularized latticeBoltzmann method , Ph.D. thesis, INP Toulouse (2018).
D. Wilde, A. Krämer, D. Reith, and H. Foysi, “Semi-Lagrangian latticeBoltzmann method for compressible flows,” Phys. Rev. E , 053306(2020).
C. Coreixas, B. Chopard, and J. Latt, “Comprehensive comparison of col-lision models in the lattice Boltzmann framework: Theoretical investiga-tions,” Phys. Rev. E , 033305 (2019).
X. Li, Y. Shi, and X. Shan, “Temperature-scaled collision process for thehigh-order lattice Boltzmann model,” Phys. Rev. E , 013301 (2019).
J. D. Sterling and S. Chen, “Stability analysis of lattice Boltzmann meth-ods,” J. Comput. Phys. , 196 – 206 (1996).
R. A. Worthing, J. Mozer, and G. Seeley, “Stability of lattice Boltzmannmethods in hydrodynamic regimes,” Phys. Rev. E , 2243–2253 (1997). P. J. Dellar, “Nonhydrodynamic modes and a priori construction of shal-low water lattice Boltzmann equations,” Phys. Rev. E , 036309 (2002). R. Adhikari and S. Succi, “Duality in matrix lattice Boltzmann models,”Phys. Rev. E , 066701 (2008). S. A. Hosseini, N. Darabiha, D. Thévenin, and A. Eshghinejadfard, “Sta-bility limits of the single relaxation-time advection–diffusion lattice Boltz-mann scheme,” Int. J. Mod. Phys. C , 1750141 (2017). G. Wissocq, P. Sagaut, and J.-F. Boussuge, “An extended spectral analysisof the lattice Boltzmann method: modal interactions and stability issues,”J. Comput. Phys. , 311 – 333 (2019).
P.-A. Masset and G. Wissocq, “Linear hydrodynamics and stability of thediscrete velocity boltzmann equations,” J. Fluid Mech. , A29 (2020).
F. Renard, G. Wissocq, J. Boussuge, and P. Sagaut, “A linear stability anal-ysis of compressible hybrid lattice Boltzmann methods,” arXiv preprintarXiv:2006.08477 (2020), arXiv:2006.08477 [physics.comp-ph].
S. Le Bras, H. Deniau, C. Bogey, and G. Daviller, “Development of com-pressible large-eddy simulations combining high-order schemes and wallmodeling,” AIAA J. , 1152–1163 (2017). P. C. Philippi, L. A. Hegele, L. O. E. dos Santos, and R. Surmas, “From thecontinuous to the lattice Boltzmann equation: The discretization problemand thermal models,” Phys. Rev. E73