Extended Lattice Boltzmann Model for Gas Dynamics
CCompressible LBM
Extended Lattice Boltzmann Model for Gas Dynamics
M. H. Saadat, S. A. Hosseini, B. Dorschner, and I. V. Karlin a) Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich,Switzerland (Dated: 19 February 2021)
We propose a two-population lattice Boltzmann model on standard lattices for the simulation of compressibleflows. The model is fully on-lattice and uses the single relaxation time Bhatnagar–Gross–Krook kinetic equa-tions along with appropriate correction terms to recover the Navier-Stokes-Fourier equations. The accuracyand performance of the model are analyzed through simulations of compressible benchmark cases includingSod shock tube, sound generation in shock-vortex interaction and compressible decaying turbulence in abox with eddy shocklets. It is demonstrated that the present model provides an accurate representation ofcompressible flows, even in the presence of turbulence and shock waves.
I. INTRODUCTION
The development of accurate and efficient numericalmethods for the simulation of compressible fluid flows re-mains a highly active research field in computational fluiddynamics (CFD), and is of great importance to manynatural phenomena and engineering applications. Com-pressiblilty is usually measured by the Mach number,
M a = u/c s , defined as the ratio of the flow velocity to thespeed of sound and is mainly characterized by the impor-tance of density and temperature variations and a dilata-tional velocity component. The presence of shock wavesin compressible flows also imposes severe challenges foran accurate numerical simulation. Shock waves are sharpdiscontinuities of the flow properties across a thin regionwith the thickness of the order of mean free path. Sincein practical simulations, it is impossible to use a grid sizefine enough to resolve the physical shock structure de-fined by the molecular viscosity, most numerical schemesrely on some numerical dissipation to stabilize the sim-ulation and capture the shock over a few grid points .The additional numerical dissipation of shock capturingschemes, however, is problematic in smooth turbulent re-gions of the flow, where a non-dissipative scheme is re-quired to capture the complex physics accurately. There-fore, in recent years, much effort has been devoted todeveloping numerical schemes capable of treating shocksand turbulence, simultaneously. This has resulted in var-ious improvements of the WENO scheme , artificial dif-fusivity approaches and hybrid schemes , to name a few.In the past decades, the lattice Boltzmann method(LBM) has received considerable attention for the CFDas a kinetic theory approach based on the discrete Boltz-mann equation. LBM has been proved to be a viable andefficient tool for the simulation of complex fluid flows andhas been applied to a wide range of fluid dynamics prob-lems including, but not limited to, turbulence , multi-phase flows and relativistic hydrodynamics . The at-tractiveness of the LBM over conventional CFD methods, a) Corresponding author; Electronic mail: [email protected] lies in the simplicity and locality of its underlying numer-ical algorithm which can be summarized as ”stream pop-ulations along the discrete velocities c i and equilibrateat the nodes x ”. It is, however, well known that LBMfaces stiff challenges in dealing with high-speed flows andits success has been mainly limited to low-speed incom-pressible flow applications.While LBM on standard lattices recovers the Navier–Stokes (NS) equations in the hydrodynamic limit, thereexist Galilean non-invariant error terms in the stress ten-sor which are negligible only in the limit of vanishingvelocities and at a singular temperature, known as thelattice temperature. This prevents LBM from going tohigher velocities as well as incorporating temperature dy-namics. A natural approach to overcome this limitationis to include more discrete velocities and use the hierar-chy of admissible high-order (or multi-speed) lattices to ensure the Galilean invariance and temperature inde-pendence of the stress tensor. Although models based onhigh-order lattices have been shown to be success-ful in simulating compressible flows to some extent, theyincrease significantly the computational cost and sufferfrom a limited temperature range , as well.Another approach, which has received considerable at-tention in recent years, maintains the simplicity and ef-ficiency of the standard lattices and employs correctionterms in order to remove the aforementioned spuriousterms in the stress tensor . Due to intrinsic non-uniqueness of the correction term, different implementa-tions exist in the literature, all recover the same equationsin the hydrodynamic limit . See Hosseini, Darabiha,and Th´evenin for a detailed review of different imple-mentations. Besides correction term, to fully recover theNavier-Stokes-Fourier (NSF) equations, one also needs toincorporate the energy equation. For doing that differ-ent models have been proposed in the literature which,in general, can be categorized into two main groups: hy-brid and two-population methods. Hybrid methods rely on solving the total energy equation using conven-tional numerical schemes like finite-difference or finite-volume. However, the majority of hybrid LB schemes suf-fer from lack of energy conservation, as the energy equa-tion is solved in a non-conservative form . In the two- a r X i v : . [ phy s i c s . f l u - dyn ] F e b ompressible LBM 2population approach , however, another popula-tion is used for the conservation of total energy. Thelatter provides a fully conservative and unified kineticframework for the compressible flows. Previous attemptsof the simulation of supersonic flows within the two-population framework on standard lattices havebeen based on the concept of shifted lattices or adap-tive lattices and need some form of interpolation duringthe streaming step. Therefore, a fully on-lattice conserva-tive scheme capable of capturing the complex physics ofcompressible flows involving shock waves is still needed.In this paper, we revisit and propose a two-populationrealization of the compressible LB model on standardlattices and investigate its accuracy and performance fora range of compressible cases from subsonic to moder-ately supersonic regime with shock waves and turbulence.The model is fully on-lattice and uses the single relax-ation time (SRT) Bhatnagar–Gross–Krook (BGK) col-lision term along with the product-form formulation ofthe equilibrium populations with a consistent correctionterm that restores the correct stress tensor. ThroughChapman-Enskog analysis, the model recovers the com-pressible NSF equations with adjustable Prandtl numberand adiabatic exponent in the hydrodynamic limit. It isshown that the model can accurately simulate compress-ible flows. Moreover, computing the correction termswith a simple upwind scheme provides enough numericaldissipation to avoid the Gibbs oscillations, and effectivelycapture the shock waves without degrading the accuracyof the scheme and overwhelming the physical dissipationin smooth regions. This is demonstrated through simu-lation of acoustic waves in the shock-vortex interactionproblem. We then investigate a more challenging case ofcompressible decaying isotropic turbulence at large tur-bulent Mach numbers and Reynolds number, where in-teraction of compressibility effects, turbulence and shocksare present in the flow field.The remainder of the paper is organized as follows:The kinetic equations of the two-population compressibleLB model along with the pertinent equilibrium and quasi-equilibrium populations are presented in Sec. II. In Sec.III, the model is validated and analyzed through simula-tion of benchmark test-cases, including Sod shock-tube,shock-vortex interaction and decaying of a compressibleisotropic turbulence. Conclusions are drawn in Sec. IV. II. MODEL DESCRIPTIONA. Kinetic equations
In the two-population approach, conservation laws aresplit between the two sets. A set of f -populations f i represents mass and momentum while another set of g -populations g i is earmarked for the energy conservation.Following Karlin, Sichau, and Chikatamarla , we con-sider a single relaxation time lattice Bhatnagar–Gross–Krook (LBGK) equations for the f -populations and a quasi-equilibrium LBM equation for the g -populations,corresponding to discrete velocities c i , where i =0 , . . . , Q − f i ( x + c i δt, t + δt ) − f i ( x , t ) = ω ( f ex i − f i ) , (1) g i ( x + c i δt, t + δt ) − g i ( x , t ) = ω ( g eq i − g i )+ ( ω − ω )( g eq i − g ∗ i ) . (2)The extended equilibrium f ex i , the equilibrium g eq i andthe quasi-equilibrium g ∗ i satisfy the local conservationlaws for the density ρ , momentum ρ u and energy ρE , ρ = Q − (cid:88) i =0 f ex i = Q − (cid:88) i =0 f i , (3) ρ u = Q − (cid:88) i =0 c i f ex i = Q − (cid:88) i =0 c i f i , (4) ρE = Q − (cid:88) i =0 g eq i = Q − (cid:88) i =0 g ∗ i = Q − (cid:88) i =0 g i . (5)We consider a general caloric equation of state of idealgas. Without loss of generality, the reference temperatureis set at T = 0 and the internal energy at unit density U is written as, U = (cid:90) T C v ( T ) dT, (6)where T is the temperature and C v ( T ) is the mass-basedspecific heat at constant volume. The energy at unitdensity E is, E = U + u . (7)The relaxation parameters ω and ω are related to viscos-ity and thermal conductivity, as it will be shown below.We now proceed with specifying the equilibria and quasi-equilibria for the standard lattice. B. Discrete velocities and factorization
We consider the D Q
27 set of three-dimensional dis-crete velocities c i , where D = 3 is the space dimensionand Q = 27 is the number of discrete speeds, c i = ( c ix , c iy , c iz ) , c iα ∈ {− , , } , i = 0 , . . . , . (8)Below, we make use of a product-form to represent allpertinent populations, the extended f -equilibrium, andthe g -equilibrium and g -quasi-equilibrium, featured inthe relaxation terms of (1) and (2). We follow Karlinand Asinari and consider a triplet of functions in twovariables ξ and P ,Ψ ( ξ, P ) = 1 − P , (9)Ψ ( ξ, P ) = 12 ( ξ + P ) , (10)Ψ − ( ξ, P ) = 12 ( − ξ + P ) . (11)ompressible LBM 3For vector-parameters ( ξ x , ξ y , ξ z ) and ( P xx , P yy , P zz ), weconsider a product associated with the speeds c i (8),Ψ i = Ψ c ix ( ξ x , P xx )Ψ c iy ( ξ y , P yy )Ψ c iz ( ξ z , P zz ) . (12)The moments of the product-form (12), M lmn = (cid:88) i =0 c lix c miy c niz Ψ i , (13)are readily computed thanks to the factorization, M lmn = M l M m M n , (14)where M = 1, and where M l = (cid:26) ξ x , l odd P xx , l even , (15) M m = (cid:26) ξ y , m odd P yy , m even , (16) M n = (cid:26) ξ z , n odd P zz , n even . (17)With the product-form (12), we proceed to specifying theextended equilibrium f -populations f ex i in (1), and theequilibrium g -populations g eq i and the quasi-equilibrium g -populations g ∗ i in (2). C. Extended f -equilibrium The extended equilibrium featured in the LBGKequation (1) has been already introduced by Saadat,Dorschner, and Karlin for the fixed temperature case.We shall summarize the construction for the purpose ofthe present compressible flow situation. At first, we de-fine the equilibrium f eq i by specifying, ξ α = u α , (18) P eq αα = RT + u α . (19)Substituting (18) and (19) into (12), we obtain, f eq i = ρ Ψ c ix ( u x , P eq xx )Ψ c iy ( u y , P eq yy )Ψ c iz ( u z , P eq zz ) . (20)The factorization (14) implies that equilibrium (20) veri-fies the maximal number Q = 27 of the moment relationsestablished by the Maxwell–Boltzmann (MB) distribu-tion, (cid:88) i =0 c lix c miy c niz f eq i = F MB lmn , l, m, n ∈ { , , } , (21)where F MB lmn = ρ (2 πRT ) − (cid:90) c lx c my c nz e − ( c − u )22 RT d c . (22) Furthermore, with (14), we find the pressure tensor andthe third-order moment tensor at the equilibrium (20), P eq = (cid:88) i =0 c i ⊗ c i f eq i = P MB , (23) Q eq = (cid:88) i =0 c i ⊗ c i ⊗ c i f eq i = Q MB + ˜ Q . (24)Here, the isotropic parts, P MB and Q MB , are theMaxwell–Boltzmann pressure tensor and the third-ordermoment tensor, respectively, P MB = P I + ρ u ⊗ u , (25) Q MB = sym( P I ⊗ u ) + ρ u ⊗ u ⊗ u , (26)where P = ρRT is the pressure, sym( . . . ) denotes sym-metrization and I is the unit tensor.The anisotropy of the equilibrium (20) manifests withthe deviation ˜ Q = Q eq − Q MB in (24), where only thediagonal elements are non-vanishing,˜ Q αβγ = (cid:26) ρu α (1 − RT ) − ρu α , if α = β = γ, , otherwise . (27)The origin of the diagonal anomaly (27) is the geometricconstraint featured by the discrete speeds (8), c iα = c iα ,for any i = 0 , . . . ,
26. Put differently, the equilibriumpressure tensor (23) and the off-diagonal elements of theequilibrium third-order moments (24) are included in theset of independent moments (21), hence they verify theMaxwell–Boltzmann moment relations by the product-form. Contrary to that, the diagonal components Q eq ααα are not among the set of moments (8), hence the anomaly.A remedy, commonly employed in the conventional LBMfor incompressible flow simulations, is to minimize thespurious effects of the said anisotropy by fixing the lat-tice reference temperature, RT L = 1 / O ( u α ) in (27). Thus, the use ofthe equilibrium (20) in the LBGK equation (1) imposesa two-fold restriction: the temperature cannot be cho-sen differently from T L while at the same time the flowvelocity has to be maintained asymptotically vanishing.While the equilibrium (20) can still be used for the ther-mal LBM in the Bussinesq approximation , they make(20) insufficient for a compressible flow setting.Instead, as was proposed by Saadat, Dorschner, andKarlin , the equilibrium (20) needs to be extended insuch a way that the third-order moment anomaly (27)is compensated in the hydrodynamic limit. Because theanomaly only concerns the diagonal elements of the third-order moments, the cancellation is achieved by redefiningthe diagonal elements of the second-order moments P αα .As was demonstrated in , in order the achieve cancel-lation of the errors, the diagonal elements P ex αα must beextended as P ex αα = P eq αα + δt (cid:18) − ω ρω (cid:19) ∂ α ˜ Q ααα , (28)ompressible LBM 4where ∂ α = ∂/∂x α and ˜ Q ααα is the diagonal element ofthe anomaly (27),˜ Q ααα = ρu α (1 − RT ) − ρu α . (29)With (28) instead of (19), the extended equilibrium f ex i is defined using the product form as before, f ex i = ρ Ψ c ix ( u x , P ex xx )Ψ c iy ( u y , P ex yy )Ψ c iz ( u z , P ex zz ) . (30)The pressure tensor of the extended equilibrium is thus P ex = P eq + δt (cid:18) − ω ω (cid:19) ∇ · ˜ Q . (31)As it has been shown in , when the extended equilib-rium (30) is used in the LBGK equation (1) at a fixed temperature T , the Navier–Stokes equation for the flowmomentum is recovered in a range of flow velocities andtemperatures. However, in the problem under considera-tion, the temperature is input from the g -population dy-namics, specifically, by solving the integral equation (7).We thus turn our attention to specifying the equilibriumand the quasi-equilibrium in the g -kinetic equation (2). D. g -equilibrium and g -quasi-equilibrium We first consider the moments of the Maxwell–Boltzmann energy distribution function, G MB lmn = ρ (2 πRT ) − (cid:90) c lx c my c nz (cid:18) c e − ( c − u )22 RT (cid:19) d c . (32)Let us introduce operators O α acting on any smoothfunction A ( u , T ) as follows , O α A = RT ∂A∂u α + u α A. (33)The Maxwell–Boltzmann energy moments (32) can bewritten as the result of repeated application of operators(33) on the generating function, G MB lmn = ρ O lx O my O nz E MB , (34)where the generating function E MB is the energy of theideal monatomic gas at unit density (three translationaldegrees of freedom, C v = (3 / R ), E MB = 32 RT + u . (35)Next, we extend the Maxwell–Boltzmann energy mo-ments (34) to a general caloric ideal gas equation of state(7). This amounts to replacing the generating function(35) with the energy (7), G eq lmn = ρ O lx O my O nz E. (36) Among the higher-order moments (36), we recognizethose pertinent to the hydrodynamic limit of the energyequation to be analyzed below. These are the equilibriumenergy flux q eq and the flux of the energy flux tensor R eq , q eq α = ρ O α E = (cid:18) H + u (cid:19) ρu α , (37) R eq αβ = ρ O α O β E = (cid:18) H + u (cid:19) P eq αβ + P u α u β . (38)Here H is the specific enthalpy, H = (cid:90) T C p ( T ) dT, (39)while C p is the specific heat at constant pressure, satis-fying Mayer’s relation, C p − C v = R .The equilibrium populations g eq i are specified with theoperator version of the product-form (12). To that end,we consider parameters ξ α and P αα as operator symbols, ξ α = O α , (40) P αα = O α . (41)With the operators (40) and (41) substituted into theproduct form (12), the equilibrium populations g eq i arewritten using the generating function (7), g eq i = ρ Ψ c ix ( O x , O x )Ψ c iy ( O y , O y )Ψ c iz ( O z , O z ) E. (42)With (14), it is straightforward to see that the equilib-rium (42) verifies a subset of the equilibrium energy mo-ments (36), (cid:88) i =0 c lix c miy c niz g eq i = G eq lmn , l, m, n ∈ { , , } . (43)Thus, by construction, the g -equilibrium (42) recoversthe maximal number Q = 27 of the energy moments(36), including the energy flux (37) and the flux of theenergy flux (38).Finally, similarly to , the quasi-equilibrium popula-tions g ∗ i are needed for adjusting the Prandtl number ofthe model. To that end, the quasi-equilibrium g ∗ i differsfrom g eq i by the non-equilibrium energy flux only, g ∗ i = g eq i + 12 c i · ( q ∗ − q eq ) , if c i = 1 ,g eq i , otherwise . (44)Here q ∗ is a specified quasi-equilibrium energy flux. In-deed, (44) and (43) imply for l, m, n ∈ { , , } , (cid:88) i =0 c lix c miy c niz g ∗ i = q ∗ x , if l = 1 , m = 0 , n = 0 q ∗ y , if l = 0 , m = 1 , n = 0 q ∗ z , if l = 0 , m = 0 , n = 1 G eq lmn , otherwise . (45)ompressible LBM 5While the above construction holds for any specified q ∗ ,the quasi-equilibrium flux required for the consistent re-alization of the adjustable Prandtl number by the LBMsystem (1) and (2) reads, q ∗ = q eq + u · (cid:18) P − P eq + δt ∇ · ˜ Q (cid:19) , (46)where P is the pressure tensor, P = (cid:88) i =0 c i ⊗ c i f i . (47)Note that unlike in the original incompressible thermalmodel , the quasi-equilibrium flux (46) now includes anextension due to the diagonal anomaly. With all theelements of the LBM system (1) and (2) specified, wenow proceed with working out its hydrodynamic limit. E. Hydrodynamic limit
Taylor expansion of the shift operator in (1) and (2)to second order gives, (cid:20) δtD i + δt D i D i (cid:21) f i = ω ( f ex i − f i ) , (48) (cid:20) δtD i + δt D i D i (cid:21) g i = ω ( g eq i − g i )+ ( ω − ω )( g eq i − g ∗ i ) , (49)where D i is the derivative along the characteristics, D i = ∂ t + c i · ∇ . (50)Introducing a multi-scale expansion, f i = f (0) i + δtf (1) i + δt f (2) i + O ( δt ) , (51) f ex i = f ex(0) i + δtf ex(1) i + δt f ex(2) i + O ( δt ) , (52) g i = f (0) i + δtf (1) i + δt f (2) i + O ( δt ) , (53) g ∗ i = g ∗ (0) i + δtg ∗ (1) i + δt g ∗ (2) i + O ( δt ) , (54) ∂ t = ∂ (1) t + δt∂ (2) t + O ( δt ) , (55)substituting into (48) and (49), and using the notation, D (1) i = ∂ (1) t + c i · ∇ , (56)we obtain, from zeroth through second order in the timestep δt , for the f -populations, f (0) i = f ex(0) i = f eq i , (57) D (1) i f (0) i = − ω (cid:16) f (1) i − f ex(1) i (cid:17) , (58) ∂ (2) t f (0) i + c i · ∇ f (1) i − ω D (1) i (cid:16) f (1) i − f ex(1) i (cid:17) = − ωf (2) i + ωf ex(2) i , (59) and similarly for the g -populations, g (0) i = g ∗ (0) i = g eq i , (60) D (1) i g (0) i = − ω g (1) i − ( ω − ω ) g ∗ (1) i , (61) ∂ (2) t g (0) i + c i · ∇ g (1) i − ω D (1) i g (1) i − ω − ω D (1) i g ∗ (1) i = − ω g (2) i − ( ω − ω ) g ∗ (2) i . (62)With (57) and (60), the mass, momentum and energyconservation (3), (4) and (5) imply the solvability condi-tions, (cid:88) i =0 f ex( k ) i = (cid:88) i =0 f ( k ) i = 0 , k = 1 , . . . ; (63) (cid:88) i =0 c i f ex( k ) i = (cid:88) i =0 c i f ( k ) i = 0 , k = 1 , , . . . ; (64) (cid:88) i =0 g ∗ ( k ) i = (cid:88) i =0 g ( k ) i = 0 , k = 1 , , . . . . (65)With the f -equilibrium (20) and the g -equilibrium (42),while taking into account the solvability conditions (63),(64) and (65), and also making use of the equilibriumpressure tensor (23) and (25), and the equilibrium energyflux (37), the first-order kinetic equations (58) and (61)imply the following first-order balance equations for thedensity, momentum and energy, ∂ (1) t ρ = − ∇ · ( ρ u ) , (66) ∂ (1) t ( ρ u ) = − ∇ · ( P I + ρ u ⊗ u ) . (67) ∂ (1) t ( ρE ) = − ∇ · q eq . (68)The first-order energy equation (68) can be recast intothe temperature equation by virtue of (66) and (67), ρC v ∂ (1) t T = − ρC v u · ∇ T − P ( ∇ · u ) . (69)Thus, to first order, the LBM recovers the compressibleEuler equations for a generic ideal gas.Moreover, the first-order constitutive relation for thenonequilibrium pressure tensor P (1) is found from (58)as follows, using (25), (24), (26) and (27), − ω P (1) + ω P ex(1) = ∂ (1) t P MB + ∇ · Q MB + ∇ · ˜ Q , (70)where P (1) = Q − (cid:88) i =0 c i ⊗ c i f (1) i , (71) P ex(1) = Q − (cid:88) i =0 c i ⊗ c i f ex(1) i . (72)Using (66), (67) and (69), we find in (70), ∂ (1) t P MB + ∇ · Q MB = Z , (73)ompressible LBM 6where we have introduced a short-hand notation for thetotal stress, including both the shear and the bulk con-tributions, Z = P (cid:18) ∇ u + ∇ u † −
23 ( ∇ · u ) I (cid:19) + P (cid:18) − RC v (cid:19) ( ∇ · u ) I , (74)and where ( · ) † denotes transposition. With (73) and (74),the nonequilibrium pressure tensor (70) becomes, P (1) = − ω Z − ω ∇ · ˜ Q + P ex(1) . (75)A comment is in order. In (75), the first term is the con-ventional contribution from both the shear and the bulkstress. The second term is anomalous due to the diag-onal anisotropy (27) while the third is the counter-termrequired to annihilate the spurious contribution in thenext, second-order approximation. According to (31), P ex(1) = (cid:18) − ω ω (cid:19) ∇ · ˜ Q . (76)Similarly, the first-order constitutive relation for thenonequilibrium energy flux q (1) is found from (61), − ω q (1) − ( ω − ω ) q ∗ (1) = ∂ (1) t q eq + ∇ · R eq . (77)Evaluating the right hand side of (77) with the help ofthe first-order relations (66), (67) and (69), we obtain, ∂ (1) t q eq + ∇ · R eq = P C p ∇ T + ( u · Z ) . (78)With (78), the nonequilibrium energy flux (77) becomes, q (1) = − ω P C p ∇ T − ω ( u · Z ) − ω − ω ω q ∗ (1) . (79)The quasi-equilibrium energy flux q ∗ (1) is evaluated ac-cording to (46) and by taking into account the first-orderconstitutive relation for the pressure tensor (75), q ∗ (1) = u · (cid:18) P (1) + 12 ∇ · ˜ Q (cid:19) = − ω ( u · Z ) . (80)We comment that the first term in the nonequilibriumenergy flux (79) is a precursor of the Fourier law of ther-mal conductivity while the second and the third termscombine to the viscous heating contribution, as we shallsee it below. The quasi-equilibrium flux (80) is requiredfor consistency of the viscous heating with the prescribedPrandtl number .With the first-order constitutive relations for thenonequilibrium fluxes (75) and (79) in place, we pro-ceed to the second-order approximation. Applying thesolvability condition (63) and (64) to the second-order f -equation (59), we obtain, ∂ (2) t ρ = 0 , (81) ∂ (2) t ( ρ u ) = − ∇ · (cid:104)(cid:16) − ω (cid:17) P (1) + ω P ex(1) (cid:105) . (82) The second-order momentum equation (82) is trans-formed by virtue of (75) and (76) to give, ∂ (2) t ( ρ u ) = − ∇ · (cid:20) − (cid:18) ω − (cid:19) Z (cid:21) . (83)Note that, the anomalous terms cancel out and the result(83) is manifestly isotropic.Finally, applying solvability condition (65) to thesecond-order g -equation (62), we find ∂ (2) t ( ρE ) = − ∇ · (cid:20)(cid:16) − ω (cid:17) q (1) − ω − ω q ∗ (1) (cid:21) . (84)Taking into account the first-order energy flux (79) andthe quasi-equilibrium energy flux (80), we obtain in (84), ∂ (2) t ( ρE ) = − ∇ · (cid:20) − (cid:18) ω − (cid:19) C p P ∇ T (cid:21) − ∇ · (cid:20) − (cid:18) ω − (cid:19) ( u · Z ) (cid:21) . (85)While the first term leads to the Fourier law, it is impor-tant to note that the second term represents viscous heat-ing consistent with the momentum equation (83). Thelatter consistency is implied by the construction of thequasi-equilibrium energy flux (46) and (80). This con-cludes the second-order accurate analysis of the hydro-dynamic limit of the LBM system (1) and (2), and weproceed with a summary of the gas dynamics equationsthereby recovered. F. Equations of gas dynamics
Combining the first- and second-order contributions tothe density, the momentum and the energy equation, (66)and (81), (67) and (83), and (68) and (85), respectively,and using a notation, ∂ t = ∂ (1) t + δt∂ (2) t , we arrive atthe continuity, the flow and the energy equations of gasdynamics as follows, ∂ t ρ + ∇ · ( ρ u ) = 0 , (86) ∂ t ( ρ u ) + ∇ · ( ρ u ⊗ u ) + ∇ · π = 0 , (87) ∂ t ( ρE ) + ∇ · ( ρE u ) + ∇ · q + ∇ · ( π · u ) = 0 . (88)Here, π is the pressure tensor π = P I − µ (cid:18) S −
23 ( ∇ · u ) I (cid:19) − ς ( ∇ · u ) I , (89)with P the pressure of ideal gas, P = ρRT, (90)with the strain rate tensor S = ∇ u + ∇ u † , (91)ompressible LBM 7and the dynamic viscosity µ and the bulk viscosity ς , µ = (cid:18) ω − (cid:19) P δt, (92) ς = (cid:18) − RC v (cid:19) µ. (93)The heat flux q in the energy equation (88) reads q = − κ ∇ T, (94)with the thermal conductivity coefficient κ , κ = (cid:18) ω − (cid:19) C p P δt. (95)The Prandtl number due to (92) and (95) is,Pr = C p µκ = ω (2 − ω ) ω (2 − ω ) , (96)while the adiabatic exponent, γ = C p C v , (97)is defined by the choice of the caloric equations of state(6) and Mayer’s relation, C p − C v = R . The mass, mo-mentum and energy equations, (86), (87) and (88) are thestandard equations of the macroscopic gas dynamics. Weshall conclude the model development with a summaryof the key elements of the LBM system (1) and (2). G. Summary of the lattice Boltzmann model
The two-population lattice Boltzmann model (1) and(2) on the standard D Q
27 discrete velocity set intro-duced by Karlin, Sichau, and Chikatamarla is extendedto the compressible flow simulation following the threekey modifications: • The product-form extended equilibrium for the mo-mentum lattice, Eqs. (12), (30), (28): f ex i = ρ Ψ c ix ( u x , P ex xx )Ψ c iy ( u y , P ex yy )Ψ c iz ( u z , P ex zz ); • The operator product-form equilibrium for the en-ergy lattice, Eqs. (7), (33), (12), (42): g eq i = ρ Ψ c ix ( O x , O x )Ψ c iy ( O y , O y )Ψ c iz ( O z , O z ) E ; • The quasi-equilibrium for the energy lattice is madeconsistent with both of the above, Eqs. (44), (46): g ∗ i = g eq i + 12 c i · ( q ∗ − q eq ) if c i = 1 ,g eq i , otherwise . We shall proceed with the implementation of the com-pressible lattice Boltzmann model and numerical valida-tion.
III. NUMERICAL RESULTSA. General implementation issues
The spatial discretization of the deviation ˜ Q inEqs.(31) and (46) has important effect on stability of themodel, especially in the case of supersonic flows wherediscontinuities emerge in the flow field. It has been shownthrough linear stability analysis that, while second-order central difference scheme provides good stabilitydomain in the subsonic regime, the first-order upwindscheme is necessary for maintaining the stability in thesupersonic regime and capturing shock wave. We, there-fore, employ the first-order upwind scheme in order tohave a wider stability domain.For example, the x -derivative of the deviation ˜ Q xxx atgrid point x i,j,k can be written as ∂ x ˜ Q xxx, ( i,j,k ) = ˜ Q xxx, ( i +1 / ,j,k ) − ˜ Q xxx, ( i − / ,j,k ) ∆ x , (98)where (omitting xxx and j, k indices) ˜ Q i +1 / and ˜ Q i − / are upwind reconstruction of ˜ Q at the interface x i ± / ,j,k ,˜ Q i +1 / = (cid:40) ˜ Q i , if u x > , ˜ Q i +1 , otherwise , (99)The performance and accuracy of the proposed LBMfor compressible flow is tested numerically through thesimulation of benchmark cases. First, the Sod shock tubeproblem is considered. Second, we show the ability of themodel in capturing moderately supersonic shock waves,through simulation of shock-vortex interaction. Finally,the model is tested with a compressible turbulence prob-lem, i.e, decaying of compressible homogeneous isotropicturbulence at different turbulent Mach numbers. All sim-ulations are performed assuming constant specific heatswith gas constant R = 1, adiabatic exponent γ = 1 . D Q
27 lattice.
B. Sod’s shock tube
Sod’s shock tube benchmark is a classical Riemannproblem, which is often used to test capability of a com-pressible flow solver in capturing shock waves, contactdiscontinuities and expansion fans. The initial flow fieldis given by( ρ, u x , P ) = (cid:26) (1 . , , . , x/L x ≤ . , (0 . , , . , x/L x > . , where L x = 600 is the number of grid points. Simulationresults with the viscosity µ = 0 .
015 for the density andreduced velocity u ∗ = u/ √ T l , where T l is temperatureon the left half of tube, at non-dimensional time t ∗ = t √ T l /L x = 0 .
2, are shown in Figs. 1 and 2. It can be seenthat, apart from a small oscillation, the results match thenon-viscous exact solution well.ompressible LBM 8 x/L x ρ PresentExact
FIG. 1. Density profile for Sod’s shock tube simulation atnon-dimensional time t ∗ = 0 .
2. Symbols: present model; line:exact solution. x/L x u * PresentExact
FIG. 2. Reduced velocity profile for Sod’s shock tube simu-lation at non-dimensional time t ∗ = 0 .
2. Symbols: presentmodel; line: exact solution.
C. Shock-vortex interaction
Sound generation by a vortex passing through a shockwave is studied to assess the performance and accu-racy of the developed model for supersonic flows involv-ing shock. This problem consists of an isentropic vor-tex, with vortex Mach number Ma v , initially in the up-stream shock region, which is passed through a station- FIG. 3. The sound pressure field ∆ P for the shock-vortexinteraction with Ma a = 1 .
2, Ma v = 0 .
25 and Re = 800 at t ∗ =6. The contour levels are from ∆ P min = − .
48 to ∆ P max =0 .
16 with an increment of 0 . ary shock wave at advection Mach number Ma a = 1 . ρ, T, u x , u y ) l = (1 , . , Ma a √ γT l , ρ ∞ , P ∞ , u x, ∞ , u y, ∞ ) is perturbed withan isentropic vortex with radius r v centered at ( x v , y v ) ρ = ρ ∞ (cid:20) − γ −
12 Ma v e (1 − r ) (cid:21) / ( γ − ,P = P ∞ (cid:20) − γ −
12 Ma v e (1 − r ) (cid:21) γ/ ( γ − ,u x = u x, ∞ + (cid:112) γT l Ma v ( y − y v ) r v e (1 − r ) / ,u y = u y, ∞ − (cid:112) γT l Ma v ( x − x v ) r v e (1 − r ) / , where r = (cid:113) ( x − x v ) + ( y − y v ) /r v is the reduced ra-dius and the shock is initially located at x s = 8 r v .We perform a simulation with Ma a = 1 .
2, Ma v = 0 . ρ L c s,l r v µ = 800, c s,l is the speed of sound upstream of the shock, and thePrandtl number is Pr = 0 .
75. The computational domainsize is L x × L y = 1680 × r v = L x /
28 and the vortex center is at ( x v , y v ) = (6 r v , L y / t ∗ = 6, where the sound pressure is defined as, ∆ P =( P − P s ) /P s , and P s is the pressure behind the shockwave. The shock wave deformation caused by the inter-action with the vortex is observed. To quantify the accu-racy of the computations, the radial sound pressure dis-tribution is plotted in Fig. 4 in comparison with the DNSompressible LBM 9 r / r v ∆ P t * = 10t * = 6 t * = 8 FIG. 4. Comparison of radial sound pressure distribution ∆ P for Ma a = 1 .
2, Ma v = 0 .
25 and Re = 800 with the DNSresults at three different times t ∗ = 6 , ,
10. Lines: presentmodel; symbol: DNS . results . The sound pressure is measured in the radialdirection with the origin at the vortex center, at an angle θ = − ◦ and at three different non-dimensional times t ∗ = 6 , ,
10, where t ∗ = tc s,l /r v . Excellent agreementis observed between the present model and the DNS .Note that the sound pressure is typically a small pertur-bation (around 1%) on top of the hydrodynamic pressure.This shows that the present model with the LBGK colli-sion term can accurately capture moderately supersonicshock waves. D. Decaying of compressible isotropic turbulence
To demonstrate that the present compressible model isa reliable method for the simulation of complex flows in-volving both turbulence and shocks, decaying compress-ible homogeneous isotropic turbulence in a periodic box isconsidered as the final test-case. This problem has beenstudied extensively and is a challenging test-case,as it contains both compressibility effects and shocks, aswell as turbulent structures in the flow field .The initial condition in a cubic domain with the side L is set to be unit density and constant temperaturealong with a divergence-free velocity field which followsthe specified energy spectrum, E ( κ ) = Aκ exp (cid:0) − κ/κ ) (cid:1) , (100)where κ is the wave number, κ is the wave number atwhich the spectrum peaks and the amplitude A controlsthe initial kinetic energy . The method of kinematicsimulation is used to generate the velocity field. FIG. 5. Iso-surface of velocity divergence ∇ · u = 0 . t = 0 .
3, Re λ = 72 and t ∗ = 0 . Control parameters for this problem are the turbulentMach number, Ma t = (cid:112) (cid:104) u · u (cid:105)(cid:104) c s (cid:105) , (101)and the Reynolds number based on the Taylor microscale,Re λ = (cid:104) ρ (cid:105) u rms λµ , (102)where u rms = (cid:112) (cid:104) u · u (cid:105) / (cid:104) . . . (cid:105) stands for the volumeaveraging over the entire computational domain while λ is the Taylor microscale, λ = (cid:104) u x (cid:105)(cid:104) ( ∂ x u x ) (cid:105) . (103)The dynamic viscosity is following a power law depen-dence on temperature, µ = µ (cid:18) TT (cid:19) . , (104)with T being the initial temperature. The Prandtl num-ber for all the simulations is Pr = 0 . .
1. Low turbulent Mach number
The simulation is first performed at a relatively lowturbulent Mach number Ma t = 0 . λ = 72, κ = 8(2 π/L ) and initial temperature T = 0 .
15. Fig.5 illustrates the instantaneous iso-surface of the veloc-ity divergence ∇ · u colored with the local Mach num-ber at the non-dimensional time t ∗ = t/τ = 0 .
4, whereompressible LBM 10 τ = L I /u rms, is the large eddy turnover time definedbased on the initial rms of the velocity and the integrallength scale, L I = 3 π (cid:82) ∞ [ E ( κ ) /κ ] dκ (cid:82) ∞ E ( κ ) dκ = √ πκ . (105)It is observed that in this case, the flow is in a moder-ately compressible to a high-subsonic regime, with themaximum local Mach number Ma max ∼ . t * M a t / M a t , DNS
FIG. 6. Decay of the turbulent Mach number for compressibledecaying turbulence at Ma t = 0 . λ = 72. Lines:present model; symbol: DNS . In order to quantify the validity of the model, a gridconvergence study is performed by using three domainsizes, 64 , 128 and 256 . The decay of the turbu-lent Mach number and of the turbulent kinetic energy K = 1 / (cid:104) ρu (cid:105) are shown in Fig. 6 and Fig. 7, where theconvergence to the DNS results can be observed.To assess the effect of compressibilty, time evolution ofthe rms of dilatation, θ rms = (cid:112) (cid:104) ( ∇ · u ) (cid:105) , (106)is compared in Fig. 8 with the DNS, where dilatationis normalized with the initial rms of vorticity, ω rms, = (cid:112) (cid:104)| ω | (cid:105) , and ω = ∇ × u . Strong compressibility effectscan be seen at the initial stage, where dilatation rapidlyincreases from its initial value θ rms, = 0, followed by amonotonic decay. Furthermore, the rms of the density ρ rms = (cid:112) (cid:104) ρ (cid:105) − (cid:104) ρ (cid:105) normalized by Ma t, is shown inFig. 9. Also here the agreement with the DNS is quitesatisfactory with 256 grid points.The enstrophy defined as,Ω = 12 (cid:104) ω (cid:105) , (107) t * K / K DNS
FIG. 7. Decay of the turbulent kinetic energy for compressibledecaying turbulence at Ma t = 0 . λ = 72. Lines:present model; symbol: DNS . t * θ r m s / ω r m s , DNS
FIG. 8. Time history of root mean square of dilitation forcompressible decaying turbulence at Ma t = 0 . λ = 72.Lines: present model; symbol: DNS . is a sensitive variable to analyze the performance of anumerical scheme for turbulent flows, as it is closely re-lated to small-scale turbulence motions . The tempo-ral evolution of the enstrophy normalized with its initialvalue is compared in Fig. 10 with the DNS results of thespectral method reported in Fang et al. . It can be seenthat in all cases the enstrophy increases in the beginningdue to vortex stretching, which generates small-scale tur-ompressible LBM 11 t * ρ r m s / M a t , DNS
FIG. 9. Time history of root mean square of density for com-pressible decaying turbulence at Ma t = 0 . λ = 72.Lines: present model; symbol: DNS . t * Ω / Ω DNS
FIG. 10. Time history of enstrophy for compressible decayingturbulence at Ma t = 0 . λ = 72. Lines: present model;symbol: DNS . bulence structures. This makes the viscous dissipationstronger, which leads to a decrease of enstrophy . Fur-thermore, coarse simulations result in under-predictionof peak value and also fast decay rate, due to strong sup-pression of small-scale fluctuations. Here, contrary to theprevious cases, 256 grid size is not enough to accuratelycapture the statistics. By increasing the resolution to512 , the peak value and decay rate of enstrophy can be captured with good accuracy. This further confirms theaccuracy of the present model in capturing the physicsof compressible turbulence. Moreover, the convergence Mesh spacing L ∞ E rr o r Second order slope
FIG. 11. Convergence rate of enstrophy for grid resolutions64 to 512 . Symbols: L ∞ error of enstrophy with respect tothe DNS results; dashed line: second-order slope. order of the model is evaluated based on the L ∞ errorof enstrophy with respect to the DNS results. As shownin Fig. 11, the overall accuracy in space is slightly belowsecond-order.
2. Effect of deviation discretization on the accuracy
As pointed out earlier, first-order upwind discretiza-tion of the deviation term is necessary for preventing theGibbs phenomenon and maintaining the stability of themodel in supersonic regime.Here, we investigate the effect of the discretizationscheme on the accuracy in subsonic turbulent regime, bycomparing the results to the case with second-order cen-tral evaluation of derivatives of deviation term. It can beseen from Fig. 12 that, the time history of enstrophy isinsensitive to the evaluation of deviation term. All otherturbulence statistics showed similar behaviour, but arenot presented here for the sake of brevity. This indicatesthat the use of first-order scheme does not degrade theformal accuracy of the solver (shown in Fig. 11), althoughit provides sufficient dissipation for stabilizing the solverand capturing the shock.
3. High turbulent Mach number
We now move on to a higher turbulent Mach number.It is well known that at sufficient high turbulent Machompressible LBM 12 t * Ω / Ω FIG. 12. Effect of deviation discretization on the enstrophyevolution for compressible decaying turbulence at Ma t = 0 . λ = 72. Lines: present model with first-order up-wind discretization of deviation term; symbols: present modelwith second-order central difference discretization of deviationterm. numbers, random shock waves commonly known as eddy-shocklets appear in the flow , due to compressiblityand turbulent motions. This scenario can, therefore, beconsidered as a rigorous test case for the validity of thepresent model.We increase the turbulent Mach number to Ma t = 0 . and 512 gridpoints and the same Reynolds number Re λ = 72. The FIG. 13. Iso-surface of velocity divergence ∇ · u = 0 . t = 0 .
4, Re λ = 72 and t ∗ = 0 . iso-surface of the velocity divergence colored by localMach number is shown in Fig. 13, which confirms thepresence of local supersonic regions during the decay. t * K / K Present 256 Present 512 DNS
FIG. 14. Decay of the turbulent kinetic energy for compress-ible decaying turbulence at Ma t = 0 . λ = 72. Lines:present model; symbol: DNS . t * ρ r m s / M a t , Present 256 Present 512 DNS
FIG. 15. Time history of root mean square of density forcompressible decaying turbulence at Ma t = 0 . λ = 72.Lines: present model; symbol: DNS . Moreover, to show that the model can accurately predictturbulent statistics in the presence of shocks, time evo-lution of the turbulent kinetic energy, rms of density andTaylor microscale Reynolds number are plotted in Figs.ompressible LBM 13 t * R e λ Present 256 Present 512 DNS
FIG. 16. Time history of Taylor microscale Reynolds numberfor compressible decaying turbulence at Ma t = 0 . λ =72. Lines: present model; symbol: DNS .
14, 16 and ?? , respectively. Here also the results agreewell with the reference DNS results.As a final validation case, we investigate the perfor-mance of the model at a relatively high Reynolds numberof Re λ = 175, while keeping the turbulent Mach numberhigh enough Ma t = 0 . κ = 4 in this case. t * K / K Present 768 DNS
FIG. 17. Decay of the turbulent kinetic energy for compress-ible decaying turbulence at Ma t = 0 .
488 and Re λ = 175.Line: present model; symbol: DNS . t * ε Present 768 DNS
FIG. 18. Time history of the dissipation rate for compressibledecaying turbulence at Ma t = 0 .
488 and Re λ = 175. Line:present model; symbol: DNS . t * R e λ Present 768 DNS
FIG. 19. Time history of Taylor microscale Reynolds num-ber for compressible decaying turbulence at Ma t = 0 .
488 andRe λ = 175. Line: present model; symbol: DNS . History of turbulent kinetic energy, solenoidal dissipa-tion rate (cid:15) = (cid:104) µω (cid:105) and Taylor miroscale Reynolds num-ber (102) are plotted in Figs. 17, 18 and 19, using 768 grid points. The results agree well with the referenceDNS solution . The energy spectrum at various timesis shown in Fig.20. It is observed that initially, largescales contain most of the energy and as time evolves theenergy is transferred to small scales. Moreover, since theompressible LBM 14 κ η E ( κ ) κ / t * = 3t * = 0 t * = 1 FIG. 20. Energy spectrum at various times ( t ∗ = 0 , Ma t = 0 .
488 and Re λ = 175. Dashed line is the initial spectrum ( t ∗ = 0).Here, η is the Kolmogorov length scale. Reynolds number is high enough, the spectrum shows theinertial range with slope of κ − / which further confirmsthe accuracy of the results and shows the ability of themodel in capturing broadband turbulent motions in thepresence of shocks. IV. CONCLUSION
In this work, we proposed a lattice Boltzmann frame-work for the simulation of compressible flows on stan-dard lattices. The product-form factorization was used torepresent all pertinent equilibrium and quasi-equilibriumpopulations. The well-known anomaly of the standardlattices was eliminated by redefining the diagonal compo-nents of the equilibrium pressure tensor through addingan appropriate correction term. The analysis of themodel was conducted through simulation of the Sod’sshock tube, sound generation in shock-vortex interactionand compressible decaying homogeneous isotropic turbu-lence.It was demonstrated that the present fully on-latticemodel with the single relaxation time LBGK collisionterm is able to properly predict the relevant featuresof the compressible flows. In particular, the modelcan capture moderately supersonic shock waves up toMa ∼ .
5. Furthermore, the simulation of compressibledecaying turbulence demonstrated that the model can ac-curately capture compressibilty effects, turbulence fluc-tuations and shocks. It was shown that the model per-forms well even at high turbulent Mach number, whereeddy-shocklets exist in the flow field and interact with turbulence. The results of the model were found to be ingood agreement with DNS results.Overall, the promising results of the proposed modelon standard lattices open interesting prospects towardsthe numerical simulation of more complex applicationssuch as compressible jet flow or shock boundary-layerinteraction . Moreover, the model could be augmentedwith more sophisticated collision terms in order to en-hance the stability of under-resolved simulations. Thiswill be the focus of our future research. ACKNOWLEDGMENTS
This work was supported by the ETH research grantETH-13 17-1 and the European Research Council (ERC)Advanced Grant No. 834763-PonD. The computationalresources at the Swiss National Super Computing CenterCSCS were provided under the grant s897. D. A. Caughey,
Computational aerodynamics (Elsevier, 2003). J. von Neumann and R. D. Richtmyer, “A method for the nu-merical calculation of hydrodynamic shocks,” Journal of appliedphysics , 232–237 (1950). C.-W. Shu, “High order ENO and WENO schemes for computa-tional fluid dynamics,” in
High-order methods for computationalphysics (Springer, 1999) pp. 439–582. A. Subramaniam, M. L. Wong, and S. K. Lele, “A high-orderweighted compact high resolution scheme with boundary closuresfor compressible turbulent flows with shocks,” Journal of Com-putational Physics , 108822 (2019). L. Fu, X. Y. Hu, and N. A. Adams, “A new class of adaptive high-order targeted ENO schemes for hyperbolic conservation laws,”Journal of Computational Physics , 724–751 (2018). L. Fu, “A very-high-order TENO scheme for all-speed gas dynam-ics and turbulence,” Computer Physics Communications ,117–131 (2019). T. Haga and S. Kawai, “On a robust and accurate localized ar-tificial diffusivity scheme for the high-order flux-reconstructionmethod,” Journal of Computational Physics , 534–563(2019). M. Visbal and D. Gaitonde, “Shock capturing using compact-differencing-based methods,” in (2005) p. 1265. B. Dorschner, F. B¨osch, S. S. Chikatamarla, K. Boulouchos, andI. V. Karlin, “Entropic multi-relaxation time lattice boltzmannmodel for complex flows,” Journal of Fluid Mechanics , 623(2016). A. Wagner, “Thermodynamic consistency of liquid-gas latticeBoltzmann simulations,” Physical Review E , 056703 (2006). M. Mendoza, B. Boghosian, H. J. Herrmann, and S. Succi, “Fastlattice Boltzmann solver for relativistic hydrodynamics,” Physi-cal Review Letters , 014502 (2010). X. Shan, X.-F. Yuan, and H. Chen, “Kinetic theory representa-tion of hydrodynamics: a way beyond the Navier–Stokes equa-tion,” Journal of Fluid Mechanics , 413–441 (2006). S. S. Chikatamarla and I. V. Karlin, “Lattices for the latticeBoltzmann method,” Physical Review E , 046701 (2009). N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, “Entropiclattice Boltzmann model for gas dynamics: Theory, boundaryconditions, and implementation,” Physical Review E , 063302(2016). D. Wilde, A. Kr¨amer, D. Reith, and H. Foysi, “Semi-Lagrangianlattice Boltzmann method for compressible flows,” Physical Re-view E , 053306 (2020). ompressible LBM 15 N. Frapolli,
Entropic lattice Boltzmann models for thermal andcompressible flows , Ph.D. thesis, ETH Zurich (2017). N. I. Prasianakis and I. V. Karlin, “Lattice Boltzmann methodfor simulation of compressible flows on standard lattices,” Phys-ical Review E , 016704 (2008). N. I. Prasianakis, I. V. Karlin, J. Mantzaras, and K. B. Boulou-chos, “Lattice Boltzmann method with restored Galilean invari-ance,” Physical Review E , 066702 (2009). Z. Guo, C. Zheng, B. Shi, and T. Zhao, “Thermal lattice Boltz-mann equation for low Mach number flows: decoupling model,”Physical Review E , 036704 (2007). Y. Feng, P. Sagaut, and W. Tao, “A three dimensional latticemodel for thermal compressible flow on standard lattices,” Jour-nal of Computational Physics , 514–529 (2015). M. H. Saadat, F. B¨osch, and I. V. Karlin, “Lattice Boltzmannmodel for compressible flows on standard lattices: VariablePrandtl number and adiabatic exponent,” Physical Review E ,013306 (2019). S. A. Hosseini, N. Darabiha, and D. Th´evenin, “Compressibilityin lattice Boltzmann on standard stencils: effects of deviationfrom reference temperature,” Philosophical Transactions of theRoyal Society A , 20190399 (2020). Y. Feng, P. Sagaut, and W.-Q. Tao, “A compressible latticeBoltzmann finite volume model for high subsonic and transonicflows on regular lattices,” Computers & Fluids , 45–55 (2016). Y. Feng, P. Boivin, J. Jacob, and P. Sagaut, “Hybrid recursiveregularized thermal lattice Boltzmann model for high subsoniccompressible flows,” Journal of Computational Physics , 82–99 (2019). S. Guo, Y. Feng, and P. Sagaut, “Improved standard thermallattice Boltzmann model with hybrid recursive regularization forcompressible laminar and turbulent flows,” Physics of Fluids ,126108 (2020). S. Zhao, G. Farag, P. Boivin, and P. Sagaut, “Toward fullyconservative hybrid lattice Boltzmann methods for compressibleflows,” Physics of Fluids , 126118 (2020). Q. Li, K. Luo, Y. He, Y. Gao, and W. Tao, “Coupling latticeBoltzmann model for simulation of thermal flows on standardlattices,” Physical Review E , 016710 (2012). I. Karlin, D. Sichau, and S. Chikatamarla, “Consistent two-population lattice Boltzmann model for thermal flows,” PhysicalReview E , 063310 (2013). B. Dorschner, F. B¨osch, and I. V. Karlin, “Particles on demandfor kinetic theory,” Physical review letters , 130602 (2018). M. H. Saadat, F. B¨osch, and I. V. Karlin, “Semi-Lagrangianlattice Boltzmann model for compressible flows on unstructuredmeshes,” Physical Review E , 023311 (2020). N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, “Lattice kinetictheory in a comoving galilean reference frame,” Physical reviewletters , 010604 (2016). P. L. Bhatnagar, E. P. Gross, and M. Krook, “A model for colli-sion processes in gases. I. Small amplitude processes in chargedand neutral one-component systems,” Physical Review , 511(1954). I. Karlin and P. Asinari, “Factorization symmetry in the latticeBoltzmann method,” Physica A: Statistical Mechanics and itsApplications , 1530–1548 (2010). M. H. Saadat, B. Dorschner, and I. V. Karlin, “Extended LatticeBoltzmann Model,” (2021), arXiv:2101.04550 [physics.flu-dyn]. G. A. Sod, “A survey of several finite difference methods forsystems of nonlinear hyperbolic conservation laws,” Journal ofComputational Physics , 1–31 (1978). O. Inoue and Y. Hattori, “Sound generation by shock-vortex in-teractions,” Journal of Fluid Mechanics , 81–116 (1999). S. Lee, S. K. Lele, and P. Moin, “Eddy shocklets in decayingcompressible turbulence,” Physics of Fluids A: Fluid Dynamics , 657–664 (1991). N. Mansour and A. Wray, “Decay of isotropic turbulence at lowReynolds number,” Physics of Fluids , 808–814 (1994). R. Samtaney, D. I. Pullin, and B. Kosovi´c, “Direct numerical sim-ulation of decaying compressible turbulence and shocklet statis-tics,” Physics of Fluids , 1415–1430 (2001). E. Johnsen, J. Larsson, A. V. Bhagatwala, W. H. Cabot, P. Moin,B. J. Olson, P. S. Rawat, S. K. Shankar, B. Sj¨ogreen, H. C. Yee, et al. , “Assessment of high-resolution methods for numerical sim-ulations of compressible turbulence with shock waves,” Journalof Computational Physics , 1213–1237 (2010). G. Kumar, S. S. Girimaji, and J. Kerimo, “WENO-enhanced gas-kinetic scheme for direct simulations of compressible transitionand turbulence,” Journal of Computational Physics , 499–523(2013). G. Cao, L. Pan, and K. Xu, “Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence I: criterion fordirect numerical simulation,” Computers & Fluids , 104273(2019). T. Chen, X. Wen, L.-P. Wang, Z. Guo, J. Wang, andS. Chen, “Simulation of three-dimensional compressible decay-ing isotropic turbulence using a redesigned discrete unified gaskinetic scheme,” Physics of Fluids , 125104 (2020). D. W. Meyer and M. L. Eggersdorfer, “Simulating particle colli-sions in homogeneous turbulence with kinematic simulation – Avalidation study,” Colloids and Surfaces A: Physicochemical andEngineering Aspects , 57–64 (2014). J. Fang, Y. Yao, Z. Li, and L. Lu, “Investigation of low-dissipation monotonicity-preserving scheme for direct numericalsimulation of compressible turbulent flows,” Computers & Fluids , 55–72 (2014). E. Garnier, M. Mossi, P. Sagaut, P. Comte, and M. Deville, “Onthe use of shock-capturing schemes for large-eddy simulation,”Journal of Computational Physics , 273–311 (1999). V. Vuorinen, J. Yu, S. Tirunagari, O. Kaario, M. Larmi,C. Duwig, and B. Boersma, “Large-eddy simulation of highlyunderexpanded transient gas jets,” Physics of Fluids , 016101(2013). S. Pirozzoli, M. Bernardini, and F. Grasso, “Direct numericalsimulation of transonic shock/boundary layer interaction underconditions of incipient separation,” Journal of Fluid Mechanics657