Extended Lattice Boltzmann Model for Gas Dynamics

We propose a two-population lattice Boltzmann model on standard lattices for the simulation of compressible flows. The model is fully on-lattice and uses the single relaxation time Bhatnagar-Gross-Krook kinetic equations along with appropriate correction terms to recover the Navier-Stokes-Fourier equations. The accuracy and performance of the model are analyzed through simulations of compressible benchmark cases including Sod shock tube, sound generation in shock-vortex interaction and compressible decaying turbulence in a box with eddy shocklets. It is demonstrated that the present model provides an accurate representation of compressible flows, even in the presence of turbulence and shock waves.


I. INTRODUCTION
The development of accurate and efficient numerical methods for the simulation of compressible fluid flows remains a highly active research field in computational fluid dynamics (CFD), and is of great importance to many natural phenomena and engineering applications. Compressiblilty is usually measured by the Mach number, M a = u/c s , defined as the ratio of the flow velocity to the speed of sound and is mainly characterized by the importance of density and temperature variations and a dilatational velocity component. The presence of shock waves in compressible flows also imposes severe challenges for an accurate numerical simulation. Shock waves are sharp discontinuities of the flow properties across a thin region with the thickness of the order of mean free path. Since in practical simulations, it is impossible to use a grid size fine enough to resolve the physical shock structure defined by the molecular viscosity, most numerical schemes rely on some numerical dissipation to stabilize the simulation and capture the shock over a few grid points 1,2 . The additional numerical dissipation of shock capturing schemes, however, is problematic in smooth turbulent regions of the flow, where a non-dissipative scheme is required to capture the complex physics accurately. Therefore, in recent years, much effort has been devoted to developing numerical schemes capable of treating shocks and turbulence, simultaneously. This has resulted in various improvements of the WENO scheme 3-6 , artificial diffusivity approaches 7 and hybrid schemes 8 , to name a few.
In the past decades, the lattice Boltzmann method (LBM) has received considerable attention for the CFD as a kinetic theory approach based on the discrete Boltzmann equation. LBM has been proved to be a viable and efficient tool for the simulation of complex fluid flows and has been applied to a wide range of fluid dynamics problems including, but not limited to, turbulence 9 , multiphase flows 10 and relativistic hydrodynamics 11 . The attractiveness of the LBM over conventional CFD methods, a) Corresponding author; Electronic mail: ikarlin@ethz.ch lies in the simplicity and locality of its underlying numerical algorithm which can be summarized as "stream populations along the discrete velocities c i and equilibrate at the nodes x". It is, however, well known that LBM faces stiff challenges in dealing with high-speed flows and its success has been mainly limited to low-speed incompressible flow applications.
While LBM on standard lattices recovers the Navier-Stokes (NS) equations in the hydrodynamic limit, there exist Galilean non-invariant error terms in the stress tensor which are negligible only in the limit of vanishing velocities and at a singular temperature, known as the lattice temperature. This prevents LBM from going to higher velocities as well as incorporating temperature dynamics. A natural approach to overcome this limitation is to include more discrete velocities and use the hierarchy of admissible high-order (or multi-speed) lattices 12,13 to ensure the Galilean invariance and temperature independence of the stress tensor. Although models based on high-order lattices 14,15 have been shown to be successful in simulating compressible flows to some extent, they increase significantly the computational cost and suffer from a limited temperature range 16 , as well.
Another approach, which has received considerable attention in recent years, maintains the simplicity and efficiency of the standard lattices and employs correction terms in order to remove the aforementioned spurious terms in the stress tensor 17,18 . Due to intrinsic nonuniqueness of the correction term, different implementations exist in the literature, all recover the same equations in the hydrodynamic limit [19][20][21] . See Hosseini, Darabiha, and Thévenin 22 for a detailed review of different implementations. Besides correction term, to fully recover the Navier-Stokes-Fourier (NSF) equations, one also needs to incorporate the energy equation. For doing that different models have been proposed in the literature which, in general, can be categorized into two main groups: hybrid and two-population methods. Hybrid methods [23][24][25] rely on solving the total energy equation using conventional numerical schemes like finite-difference or finitevolume. However, the majority of hybrid LB schemes suffer from lack of energy conservation, as the energy equation is solved in a non-conservative form 26 . In the two-population approach 19,21,27,28 , however, another population is used for the conservation of total energy. The latter provides a fully conservative and unified kinetic framework for the compressible flows. Previous attempts of the simulation of supersonic flows within the twopopulation framework on standard lattices 21,29,30 have been based on the concept of shifted lattices 31 or adaptive lattices 29 and need some form of interpolation during the streaming step. Therefore, a fully on-lattice conservative scheme capable of capturing the complex physics of compressible flows involving shock waves is still needed.
In this paper, we revisit and propose a two-population realization of the compressible LB model on standard lattices and investigate its accuracy and performance for a range of compressible cases from subsonic to moderately supersonic regime with shock waves and turbulence. The model is fully on-lattice and uses the single relaxation time (SRT) Bhatnagar-Gross-Krook 32 (BGK) collision term along with the product-form formulation 33 of the equilibrium populations with a consistent correction term that restores the correct stress tensor. Through Chapman-Enskog analysis, the model recovers the compressible NSF equations with adjustable Prandtl number and adiabatic exponent in the hydrodynamic limit. It is shown that the model can accurately simulate compressible flows. Moreover, computing the correction terms with a simple upwind scheme provides enough numerical dissipation to avoid the Gibbs oscillations, and effectively capture the shock waves without degrading the accuracy of the scheme and overwhelming the physical dissipation in smooth regions. This is demonstrated through simulation of acoustic waves in the shock-vortex interaction problem. We then investigate a more challenging case of compressible decaying isotropic turbulence at large turbulent Mach numbers and Reynolds number, where interaction of compressibility effects, turbulence and shocks are present in the flow field.
The remainder of the paper is organized as follows: The kinetic equations of the two-population compressible LB model along with the pertinent equilibrium and quasiequilibrium populations are presented in Sec. II. In Sec. III, the model is validated and analyzed through simulation of benchmark test-cases, including Sod shock-tube, shock-vortex interaction and decaying of a compressible isotropic turbulence. Conclusions are drawn in Sec. IV.

A. Kinetic equations
In the two-population approach, conservation laws are split between the two sets. A set of f -populations f i represents mass and momentum while another set of gpopulations g i is earmarked for the energy conservation. Following Karlin, Sichau, and Chikatamarla 28 , we consider 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 − 1, The extended equilibrium f ex i , the equilibrium g eq i and the quasi-equilibrium g * i satisfy the local conservation laws for the density ρ, momentum ρu and energy ρE, We consider a general caloric equation of state of ideal gas. Without loss of generality, the reference temperature is set at T = 0 and the internal energy at unit density U is written as, where T is the temperature and C v (T ) is the mass-based specific heat at constant volume. The energy at unit density E is, The relaxation parameters ω and ω 1 are related to viscosity and thermal conductivity, as it will be shown below. We now proceed with specifying the equilibria and quasiequilibria for the standard lattice.

B. Discrete velocities and factorization
We consider the D3Q27 set of three-dimensional discrete velocities c i , where D = 3 is the space dimension and Q = 27 is the number of discrete speeds, Below, we make use of a product-form to represent all pertinent populations, the extended f -equilibrium, and the g-equilibrium and g-quasi-equilibrium, featured in the relaxation terms of (1) and (2). We follow Karlin and Asinari 33 and consider a triplet of functions in two variables ξ and P, For vector-parameters (ξ x , ξ y , ξ z ) and (P xx , P yy , P zz ), we consider a product associated with the speeds c i (8), The moments of the product-form (12), are readily computed thanks to the factorization, where M 000 = 1, and where M 00n = ξ z , n odd P zz , n even .
With the product-form (12), we proceed to specifying the extended equilibrium f -populations f ex i in (1), and the equilibrium 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 LBGK equation (1) has been already introduced by Saadat, Dorschner, and Karlin 34 for the fixed temperature case. We shall summarize the construction for the purpose of the present compressible flow situation. At first, we define the equilibrium f eq i by specifying, Substituting (18) and (19) into (12), we obtain, The factorization (14) implies that equilibrium (20) verifies the maximal number Q = 27 of the moment relations established by the Maxwell-Boltzmann (MB) distribution, where Furthermore, with (14), we find the pressure tensor and the third-order moment tensor at the equilibrium (20), Here, the isotropic parts, P MB and Q MB , are the Maxwell-Boltzmann pressure tensor and the third-order moment tensor, respectively, where P = ρRT is the pressure, sym(. . . ) denotes symmetrization and I is the unit tensor. The anisotropy of the equilibrium (20) manifests with the deviationQ = Q eq − Q MB in (24), where only the diagonal elements are non-vanishing, The origin of the diagonal anomaly (27) is the geometric constraint featured by the discrete speeds (8), c 3 iα = c iα , for any i = 0, . . . , 26. Put differently, the equilibrium pressure tensor (23) and the off-diagonal elements of the equilibrium third-order moments (24) are included in the set of independent moments (21), hence they verify the Maxwell-Boltzmann moment relations by the productform. 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 LBM for incompressible flow simulations, is to minimize the spurious effects of the said anisotropy by fixing the lattice reference temperature, RT L = 1/3 in order to eliminate the linear term O(u α ) in (27). Thus, the use of the equilibrium (20) in the LBGK equation (1) imposes a two-fold restriction: the temperature cannot be chosen differently from T L while at the same time the flow velocity has to be maintained asymptotically vanishing. While the equilibrium (20) can still be used for the thermal LBM in the Bussinesq approximation 28 , they make (20) insufficient for a compressible flow setting.
Instead, as was proposed by Saadat, Dorschner, and Karlin 34 , the equilibrium (20) needs to be extended in such a way that the third-order moment anomaly (27) is compensated in the hydrodynamic limit. Because the anomaly only concerns the diagonal elements of the thirdorder moments, the cancellation is achieved by redefining the diagonal elements of the second-order moments P αα . As was demonstrated in 34 , in order the achieve cancellation of the errors, the diagonal elements P ex αα must be extended as where ∂ α = ∂/∂x α andQ ααα is the diagonal element of the anomaly (27), With (28) instead of (19), the extended equilibrium f ex i is defined using the product form as before, The pressure tensor of the extended equilibrium is thus As it has been shown in 34 , when the extended equilibrium (30) is used in the LBGK equation (1) at a fixed temperature T , the Navier-Stokes equation for the flow momentum is recovered in a range of flow velocities and temperatures. However, in the problem under consideration, the temperature is input from the g-population dynamics, specifically, by solving the integral equation (7). We thus turn our attention to specifying the equilibrium and 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, Let us introduce operators O α acting on any smooth function A(u, T ) as follows 28 , The Maxwell-Boltzmann energy moments (32) can be written as the result of repeated application of operators (33) on the generating function, where the generating function E MB is the energy of the ideal monatomic gas at unit density (three translational degrees of freedom, Next, we extend the Maxwell-Boltzmann energy moments (34) to a general caloric ideal gas equation of state (7). This amounts to replacing the generating function (35) with the energy (7), Among the higher-order moments (36), we recognize those pertinent to the hydrodynamic limit of the energy equation to be analyzed below. These are the equilibrium energy flux q eq and the flux of the energy flux tensor R eq , Here H is the specific enthalpy, while C p is the specific heat at constant pressure, satisfying Mayer's relation, The equilibrium populations g eq i are specified with the operator version of the product-form (12). To that end, we consider parameters ξ α and P αα as operator symbols, With the operators (40) and (41) substituted into the product form (12), the equilibrium populations g eq i are written using the generating function (7), With (14), it is straightforward to see that the equilibrium (42) verifies a subset of the equilibrium energy moments (36), Thus, by construction, the g-equilibrium (42) recovers the maximal number Q = 27 of the energy moments (36), including the energy flux (37) and the flux of the energy flux (38). Finally, similarly to 28 , the quasi-equilibrium populations g * i are needed for adjusting the Prandtl number of the model. To that end, the quasi-equilibrium g * i differs from g eq i by the non-equilibrium energy flux only, Here q * is a specified quasi-equilibrium energy flux. Indeed, (44) and (43) imply for l, m, n ∈ {0, 1, 2}, While the above construction holds for any specified q * , the quasi-equilibrium flux required for the consistent realization of the adjustable Prandtl number by the LBM system (1) and (2) reads, where P is the pressure tensor, Note that unlike in the original incompressible thermal model 28 , the quasi-equilibrium flux (46) now includes an extension due to the diagonal anomaly. With all the elements of the LBM system (1) and (2) specified, we now proceed with working out its hydrodynamic limit.

E. Hydrodynamic limit
Taylor expansion of the shift operator in (1) and (2) to second order gives, where D i is the derivative along the characteristics, Introducing a multi-scale expansion, substituting into (48) and (49), and using the notation, we obtain, from zeroth through second order in the time step δt, for the f -populations, and similarly for the g-populations, With (57) and (60), the mass, momentum and energy conservation (3), (4) and (5) imply the solvability conditions, 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 equilibrium pressure tensor (23) and (25), and the equilibrium energy flux (37), the first-order kinetic equations (58) and (61) imply the following first-order balance equations for the density, momentum and energy, The first-order energy equation (68) can be recast into the temperature equation by virtue of (66) and (67), Thus, to first order, the LBM recovers the compressible Euler equations for a generic ideal gas. Moreover, the first-order constitutive relation for the nonequilibrium pressure tensor P (1) is found from (58) as follows, using (25), (24), (26) and (27), where Using (66), (67) and (69), we find in (70), where we have introduced a short-hand notation for the total stress, including both the shear and the bulk contributions, and where (·) † denotes transposition. With (73) and (74), the nonequilibrium pressure tensor (70) becomes, A comment is in order. In (75), the first term is the conventional contribution from both the shear and the bulk stress. The second term is anomalous due to the diagonal anisotropy (27) while the third is the counter-term required to annihilate the spurious contribution in the next, second-order approximation. According to (31), Similarly, the first-order constitutive relation for the nonequilibrium energy flux q (1) is found from (61), Evaluating the right hand side of (77) with the help of the first-order relations (66), (67) and (69), we obtain, With (78), the nonequilibrium energy flux (77) becomes, The quasi-equilibrium energy flux q * (1) is evaluated according to (46) and by taking into account the first-order constitutive relation for the pressure tensor (75), We comment that the first term in the nonequilibrium energy flux (79) is a precursor of the Fourier law of thermal conductivity while the second and the third terms combine to the viscous heating contribution, as we shall see it below. The quasi-equilibrium flux (80) is required for consistency of the viscous heating with the prescribed Prandtl number 28 .
With the first-order constitutive relations for the nonequilibrium fluxes (75) and (79) in place, we proceed to the second-order approximation. Applying the solvability condition (63) and (64) to the second-order f -equation (59), we obtain, The second-order momentum equation (82) is transformed by virtue of (75) and (76) to give, Note that, the anomalous terms cancel out and the result (83) is manifestly isotropic. Finally, applying solvability condition (65) to the second-order g-equation (62), we find Taking into account the first-order energy flux (79) and the quasi-equilibrium energy flux (80), we obtain in (84), While the first term leads to the Fourier law, it is important to note that the second term represents viscous heating consistent with the momentum equation (83). The latter consistency is implied by the construction of the quasi-equilibrium energy flux (46) and (80). This concludes the second-order accurate analysis of the hydrodynamic limit of the LBM system (1) and (2), and we proceed with a summary of the gas dynamics equations thereby recovered.

F. Equations of gas dynamics
Combining the first-and second-order contributions to the 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 at the continuity, the flow and the energy equations of gas dynamics as follows, Here, π is the pressure tensor with P the pressure of ideal gas, with the strain rate tensor and the dynamic viscosity µ and the bulk viscosity ς, The heat flux q in the energy equation (88) reads with the thermal conductivity coefficient κ, The Prandtl number due to (92) and (95) is, while the adiabatic exponent, is defined by the choice of the caloric equations of state (6) and Mayer's relation, C p − C v = R. The mass, momentum and energy equations, (86), (87) and (88) are the standard equations of the macroscopic gas dynamics. We shall conclude the model development with a summary of 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 D3Q27 discrete velocity set introduced by Karlin, Sichau, and Chikatamarla 28 is extended to the compressible flow simulation following the three key modifications: • The product-form extended equilibrium for the momentum lattice, Eqs. (12), (30), (28): • The operator product-form equilibrium for the energy lattice, Eqs. (7), (33), (12), (42): • The quasi-equilibrium for the energy lattice is made consistent with both of the above, Eqs. (44), (46): We shall proceed with the implementation of the compressible lattice Boltzmann model and numerical validation.

A. General implementation issues
The spatial discretization of the deviationQ in Eqs. (31) and (46) has important effect on stability of the model, especially in the case of supersonic flows where discontinuities emerge in the flow field. It has been shown through linear stability analysis 22 that, while secondorder central difference scheme provides good stability domain in the subsonic regime, the first-order upwind scheme is necessary for maintaining the stability in the supersonic regime and capturing shock wave. We, therefore, employ the first-order upwind scheme in order to have a wider stability domain.
For example, the x-derivative of the deviationQ xxx at grid point x i,j,k can be written as where (omitting xxx and j, k indices)Q i+1/2 andQ i−1/2 are upwind reconstruction ofQ at the interface x i±1/2,j,k , The performance and accuracy of the proposed LBM for compressible flow is tested numerically through the simulation of benchmark cases. First, the Sod shock tube problem is considered. Second, we show the ability of the model in capturing moderately supersonic shock waves, through simulation of shock-vortex interaction. Finally, the model is tested with a compressible turbulence problem, i.e, decaying of compressible homogeneous isotropic turbulence at different turbulent Mach numbers. All simulations are performed assuming constant specific heats with gas constant R = 1, adiabatic exponent γ = 1.4 and the D3Q27 lattice.

B. Sod's shock tube
Sod's shock tube benchmark 35 is a classical Riemann problem, which is often used to test capability of a compressible flow solver in capturing shock waves, contact discontinuities and expansion fans. The initial flow field is given by

C. Shock-vortex interaction
Sound generation by a vortex passing through a shock wave 36 is studied to assess the performance and accuracy of the developed model for supersonic flows involving shock. This problem consists of an isentropic vortex, with vortex Mach number Ma v , initially in the upstream shock region, which is passed through a station- ary shock wave at advection Mach number Ma a = 1.2 with the left state (ρ, T, u x , u y ) l = (1, 0.05, Ma a √ γT l , 0) and Rankine-Hugoniot right state. The initial field with standing shock (ρ ∞ , P ∞ , u x,∞ , u y,∞ ) is perturbed with an isentropic vortex with radius r v centered at (x v , y v ) 36 dius and the shock is initially located at x s = 8r v . We perform a simulation with Ma a = 1.2, Ma v = 0.25, where the Reynolds number is set to Re = ρ L c s,l rv µ = 800, c s,l is the speed of sound upstream of the shock, and the Prandtl number is Pr = 0.75. The computational domain size is L x × L y = 1680 × 1440, the vortex radius is r v = L x /28 and the vortex center is at (x v , y v ) = (6r v , L y /2). Fig. 3 shows the sound pressure contours at time t * = 6, where the sound pressure is defined as, ∆P = (P − P s )/P s , and P s is the pressure behind the shock wave. The shock wave deformation caused by the interaction with the vortex is observed. To quantify the accuracy of the computations, the radial sound pressure distribution is plotted in Fig. 4  Note that the sound pressure is typically a small perturbation (around 1%) on top of the hydrodynamic pressure. This shows that the present model with the LBGK collision term can accurately capture moderately supersonic shock waves.

D. Decaying of compressible isotropic turbulence
To demonstrate that the present compressible model is a reliable method for the simulation of complex flows involving both turbulence and shocks, decaying compressible homogeneous isotropic turbulence in a periodic box is considered as the final test-case. This problem has been studied extensively [37][38][39][40][41][42][43] and is a challenging test-case, as it contains both compressibility effects and shocks, as well as turbulent structures in the flow field 40 .
The initial condition in a cubic domain with the side L is set to be unit density and constant temperature along with a divergence-free velocity field which follows the specified energy spectrum, where κ is the wave number, κ 0 is the wave number at which the spectrum peaks and the amplitude A controls the initial kinetic energy 39 . The method of kinematic simulation 44 is used to generate the velocity field. Control parameters for this problem are the turbulent Mach number, and the Reynolds number based on the Taylor microscale, where u rms = u · u /3 is the root mean square (rms) of the velocity and notation . . . stands for the volume averaging over the entire computational domain while λ is the Taylor microscale, The dynamic viscosity is following a power law dependence on temperature, with T 0 being the initial temperature. The Prandtl number for all the simulations is Pr = 0.7 in accordance with the DNS 39 .

Low turbulent Mach number
The simulation is first performed at a relatively low turbulent Mach number Ma t = 0.3 with Re λ = 72, κ 0 = 8(2π/L) and initial temperature T 0 = 0.15. Fig.  5 illustrates the instantaneous iso-surface of the velocity divergence ∇ · u colored with the local Mach number at the non-dimensional time t * = t/τ = 0.4, where τ = L I /u rms,0 is the large eddy turnover time defined based on the initial rms of the velocity and the integral length scale, It is observed that in this case, the flow is in a moderately compressible to a high-subsonic regime, with the maximum local Mach number Ma max ∼ 0.8. In order to quantify the validity of the model, a grid convergence study is performed by using three domain sizes, 64 3 , 128 3 and 256 3 . The decay of the turbulent Mach number and of the turbulent kinetic energy K = 1/2 ρu 2 are shown in Fig. 6 and Fig. 7, where the convergence to the DNS results 39 can be observed.
To assess the effect of compressibilty, time evolution of the rms of dilatation, is compared in Fig. 8 with the DNS, where dilatation is normalized with the initial rms of vorticity, ω rms,0 = |ω 0 | 2 , and ω = ∇ × u. Strong compressibility effects can be seen at the initial stage, where dilatation rapidly increases from its initial value θ rms,0 = 0, followed by a monotonic decay. Furthermore, the rms of the density ρ rms = ρ 2 − ρ 2 normalized by Ma 2 t,0 is shown in Fig. 9. Also here the agreement with the DNS is quite satisfactory with 256 3 grid points.
The enstrophy defined as, is a sensitive variable to analyze the performance of a numerical scheme for turbulent flows, as it is closely related to small-scale turbulence motions 45,46 . The temporal evolution of the enstrophy normalized with its initial value is compared in Fig. 10 with the DNS results of the spectral method reported in Fang et al. 45 . It can be seen that in all cases the enstrophy increases in the beginning due to vortex stretching, which generates small-scale tur- bulence structures. This makes the viscous dissipation stronger, which leads to a decrease of enstrophy 46 . Furthermore, coarse simulations result in under-prediction of peak value and also fast decay rate, due to strong suppression of small-scale fluctuations. Here, contrary to the previous cases, 256 3 grid size is not enough to accurately capture the statistics. By increasing the resolution to 512 3 , the peak value and decay rate of enstrophy can be captured with good accuracy. order of the model is evaluated based on the L ∞ error of enstrophy with respect to the DNS results. As shown in Fig. 11, the overall accuracy in space is slightly below second-order.

Effect of deviation discretization on the accuracy
As pointed out earlier, first-order upwind discretization of the deviation term is necessary for preventing the Gibbs phenomenon and maintaining the stability of the model in supersonic regime.
Here, we investigate the effect of the discretization scheme on the accuracy in subsonic turbulent regime, by comparing the results to the case with second-order central evaluation of derivatives of deviation term. It can be seen from Fig. 12 that, the time history of enstrophy is insensitive to the evaluation of deviation term. All other turbulence statistics showed similar behaviour, but are not presented here for the sake of brevity. This indicates that the use of first-order scheme does not degrade the formal accuracy of the solver (shown in Fig. 11), although it provides sufficient dissipation for stabilizing the solver and capturing the shock.

High turbulent Mach number
We now move on to a higher turbulent Mach number. It is well known that at sufficient high turbulent Mach numbers, random shock waves commonly known as eddyshocklets appear in the flow 37,39,40 , due to compressiblity and turbulent motions. This scenario can, therefore, be considered as a rigorous test case for the validity of the present model. We increase the turbulent Mach number to Ma t = 0.5 and perform the simulation with 256 3 and 512 3 grid points and the same Reynolds number Re λ = 72. The iso-surface of the velocity divergence colored by local Mach number is shown in Fig. 13, which confirms the presence of local supersonic regions during the decay. 14, 16 and ??, respectively. Here also the results agree well with the reference DNS results. As a final validation case, we investigate the performance of the model at a relatively high Reynolds number of Re λ = 175, while keeping the turbulent Mach number high enough Ma t = 0.488. The initial spectrum peaks at κ 0 = 4 in this case. History of turbulent kinetic energy, solenoidal dissipation rate = µω 2 and Taylor miroscale Reynolds number (102) are plotted in Figs. 17, 18 and 19, using 768 3 grid points. The results agree well with the reference DNS solution 39 . The energy spectrum at various times is shown in Fig.20. It is observed that initially, large scales contain most of the energy and as time evolves the energy is transferred to small scales. Moreover, since the κ η E(κ) Reynolds number is high enough, the spectrum shows the inertial range with slope of κ −5/3 which further confirms the accuracy of the results and shows the ability of the model in capturing broadband turbulent motions in the presence of shocks.

IV. CONCLUSION
In this work, we proposed a lattice Boltzmann framework for the simulation of compressible flows on standard lattices. The product-form factorization was used to represent all pertinent equilibrium and quasi-equilibrium populations. The well-known anomaly of the standard lattices was eliminated by redefining the diagonal components of the equilibrium pressure tensor through adding an appropriate correction term. The analysis of the model was conducted through simulation of the Sod's shock tube, sound generation in shock-vortex interaction and compressible decaying homogeneous isotropic turbulence.
It was demonstrated that the present fully on-lattice model with the single relaxation time LBGK collision term is able to properly predict the relevant features of the compressible flows. In particular, the model can capture moderately supersonic shock waves up to Ma ∼ 1.5. Furthermore, the simulation of compressible decaying turbulence demonstrated that the model can accurately capture compressibilty effects, turbulence fluctuations and shocks. It was shown that the model performs well even at high turbulent Mach number, where eddy-shocklets exist in the flow field and interact with turbulence. The results of the model were found to be in good agreement with DNS results.
Overall, the promising results of the proposed model on standard lattices open interesting prospects towards the numerical simulation of more complex applications such as compressible jet flow 47 or shock boundary-layer interaction 48 . Moreover, the model could be augmented with more sophisticated collision terms in order to enhance the stability of under-resolved simulations. This will be the focus of our future research.