Extended Lattice Boltzmann Model

Conventional lattice Boltzmann models for the simulation of fluid dynamics are restricted by an error in the stress tensor that is negligible only for small flow velocity and at a singular value of the temperature. To that end, we propose a unified formulation that restores Galilean invariance and the isotropy of the stress tensor by introducing an extended equilibrium. This modification extends lattice Boltzmann models to simulations with higher values of the flow velocity and can be used at temperatures that are higher than the lattice reference temperature, which enhances computational efficiency by decreasing the number of required time steps. Furthermore, the extended model also remains valid for stretched lattices, which are useful when flow gradients are predominant in one direction. The model is validated by simulations of two- and three-dimensional benchmark problems, including the double shear layer flow, the decay of homogeneous isotropic turbulence, the laminar boundary layer over a flat plate and the turbulent channel flow.


Introduction
The lattice Boltzmann method (LBM) solves a Boltzmann-type kinetic equation on a discrete velocity set, forming the links of a space-filling lattice. Efficiency of the LBM makes it attractive for the simulation of a wide range of problems in fluid dynamics, see, for example, [1,2].
In this paper, we revisit the restrictions of LBM due to the geometry of the discrete velocities. It is well known that standard LBM velocities yield a persistent error in the fluid stress tensor, which breaks Galilean invariance and limits the operation range of LBM to small flow velocities and a singular value of the lattice reference temperature; only under these conditions can the error be ignored. While one can cope with this error in most incompressible flow applications [1,2], it certainly affects high-speed compressible flows [3][4][5][6][7][8][9][10][11][12] and sometimes even low-speed isothermal cases [13]. Moreover, the same error is amplified when using stretched (rectangular) lattices instead of the conventional (cubic) lattice, where in addition to the corrupted Galilean invariance, the stress tensor becomes anisotropic [14,15].
The extension of LBM beyond its classical operation domain has so far been addressed with different techniques, depending on the desired outcome. Introducing lattices with more velocities (multi-speed lattices) [16] is one technique, which comes with a significant increase of computational cost and has been used mostly for high-speed compressible flow applications [16]. In the standard LB setting, on the other hand, one approach is to alter the relaxation rates and use a multi-relaxation time collision operator [9]. Another approach to extend the flow velocity and temperature range of the standard cubic lattices is to add correction terms to the original LBM [4][5][6][7][8][10][11][12][17][18][19][20][21]. The realization varies among different authors but none address the general case of rectangular lattices. On the other hand, rectangular lattices may improve the computational efficiency of the LBM by using a coarser mesh in the direction of smaller gradients in the flow. Unlike other approaches of handling non-uniform grids (e.g., Eulerian [22,23] and semi-Lagrangian offlattice LBM [24][25][26][27] or grid refinement techniques [28,29]), stretched lattices do not require a substantial change in the standard LBM algorithm. Recent work on the stretched LBM restores the isotropy of the stress tensor by using multi-relaxation time LBM models [30][31][32]. However, these approaches do not address the flow velocity and temperature restrictions.
In this paper, we propose a unified view on the three aspects of the problem, the velocity range, the temperature range and grid stretching, which all stem from the same error, induced by constraints of the discrete velocity set. In particular, we propose to use an extended equilibrium, which restores the Galilean invariance and isotropy of the stress tensor, enabling simulations at higher flow velocities, higher temperatures using both cubic and stretched lattices, yielding increased accuracy and efficiency.
The paper is organized as follows-in Section 2, we start by presenting the discrete kinetic equations, following the standard single-relaxation time lattice Bhatnagar-Gross-Krook (LBGK) setting, as well as the equilibrium and extended equilibrium formulation, followed by the derivation of the model's hydrodynamic limit. In addition, the locally corrected LBM of Reference [9] is compared to the extended LBGK in Appendix A.
Subsequently, in Section 3, we assess the validity, accuracy and performance of our model using both two-and three-dimensional benchmark problems. As a first step, we verify Galilean invariance, temperature independence and isotropy of the model on the example of an advected decaying shear wave. It is shown that the theoretical viscosity is recovered for both cubic as well as stretched lattices in a large range of temperatures and advection velocities. This also indicates that the model can readily be extended to high-speed compressible flows, provided that it is augmented with a suitable solver for the total energy. Next, for the example of homogeneous isotropic turbulence, we demonstrate that a speed-up can be achieved by using an operating temperature, which is larger than the lattice reference temperature. The present model can also be viewed as an alternative to preconditioned LBM [33] for accelerating the convergence rate but without the restriction to steady flows. Finally, accuracy and performance are assessed for rectangular lattices using the doubly periodic shear flow, laminar flow over a flat plate and turbulent channel flow as examples. Conclusions are drawn in Section 4.

LBGK
We consider the LBGK equation for the populations f i , associated to the discrete velocities v i for i = 0, . . . , Q − 1, where x denotes the location in space and δt is the time step. The extended equilibrium f ex i , which will be specified below, satisfies the local conservation laws for the density ρ and momentum ρu, As we will show below, the relaxation parameter ω is related to the kinematic viscosity ν, where T is the temperature and R is the gas constant. We now proceed with identifying the extended equilibrium.

Discrete Velocities and Factorization
We use the D3Q27 lattice, where D = 3 denotes the spatial dimension and Q = 27 is the number of discrete speeds, which are given by With (5), we define the particles' velocities v i in a stretched cell as where λ α is the stretching factor in the direction α. To maintain on-lattice propagation, the cell size is changed accordingly to δx α = λ α δt. The (normalized, M 000 = 1) moments M lmn are defined using the convention l → x, m → y, n → z; l, m, n = 0, 1, 2, . . . , and thus For convenience, we use a more specific notation for the first-order and the diagonal second-order moments, M 200 = P xx , M 020 = P yy , M 002 = P zz .
We essentially follow [34] and consider a class of factorized populations. To that end, we define a triplet of functions in the three variables, u α , P αα and λ α , For the vectors u, P, and λ, P = (P xx , P yy , P zz ), we consider a product-form, associated with the discrete velocities v i (6), The normalized moments of the product-form (17), are readily computed thanks to the factorization, where For any stretching (16), the six-parametric family of normalized populations (17) is identified by the flow velocity (14) and the diagonal of the pressure tensor at unit density (15), and was termed the unidirectional quasi-equilibrium in Ref. [34]. We make use of the product-form (17) to construct all pertinent populations, the equilibrium and the extended equilibrium.

Equilibrium and Extended Equilibrium
The equilibrium f eq i is defined by setting P αα (10) equal to the equilibrium diagonal element of the pressure tensor at unit density, Substituting (24) into (17), we get With (18), we find the pressure tensor and the third-order moment tensor at the equilibrium (25) as follows, The isotropic parts, P MB and Q MB , are the Maxwell-Boltzmann (MB) 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 (25) manifests with the deviationQ = Q eq − Q MB in (29), where only the diagonal elements are non-vanishing, The origin of the diagonal anomaly (30) is the geometric constraint, v 3 iα = λ 2 α v iα , which is imposed by the choice of the discrete speeds (5), and is well known in the case of the standard (unstretched) lattice with λ α = 1. A remedy in the latter case is to minimize spurious effects of anisotropy by fixing the temperature T = T L , In order to eliminate the linear term ∼ u α in (30). Thus, using the equilibrium (25) in the LBGK Equation (1) imposes a two-fold restriction on the operation domain: the temperature cannot be chosen differently from (31) and the flow velocity has to be maintained asymptotically vanishing. Moreover, for stretched lattices, the anisotropy becomes even more pronounced since it is not possible to eliminate the linear deviation in all directions simultaneously by fixing any temperature. Alternatively, the spurious anisotropy effects can be canceled out by extending the equilibrium such that the third-order moment anomaly is compensated in the hydrodynamic limit. Because the anomaly only concerns the diagonal (unidirectional) elements of the third-order moments, the cancellation can be achieved by redefining the diagonal elements of the second-order moments. As demonstrated below, in order to cancel the errors, the diagonal elements P ex αα for the extended equilibrium must be chosen as where spatial derivative is evaluated using a second-order central difference scheme. Hence, the extended equilibrium f ex i is specified by using the product-form (17), We shall now proceed with the derivation of the Navier-Stokes equations in the hydrodynamic limit of the proposed extended LBGK model.

Hydrodynamic Limit with Extended Equilibrium
Taylor expansion of the shift operator in (1) to second order gives, where D i is the derivative along the characteristics, Introducing the multi-scale expansion, substituting into (34) and using the notation, we obtain, from the zeroth to second order in the time step δt, With (40), the mass and the momentum conservation (2) and (3) imply the solvability conditions, With the equilibrium (25), taking into account the solvability conditions (43) and (44), and also making use of the equilibrium pressure tensor (26) and (28), the first-order Equation (41) implies the following relations for the density and the momentum, Moreover, the first-order constitutive relation for the nonequilibrium pressure tensor P (1) is found from (41) as follows, where With the help of Equations (28) and (29), the first-order constitutive relation (47) is transformed to make explicit the contribution of the anomalous term (30), The last term is evaluated using (45) and (46) to give, where (·) † denotes transposition. Combining (51) and (50), the first-order constitutive relation becomes, Note that, if we had used the equilibrium f eq i instead of the extended equilibrium f ex i in (1), at this stage of the derivation, we get instead of (52), The anomalous term ∇ ·Q cannot be canceled in the latter expression, rather, by choosing T = T L (31), its effect can be ignored but only under the assumption of an asymptotically small flow velocity. Moreover, for a quasi-incompressible flow (Ma → 0, density variation ∇ρ ∼ Ma 2 , where Ma is a characteristic Mach number), it is possible to further reduce the effect of the anomaly by rescaling the relaxation parameter [9], see a discussion in Appendix A.
In a contrast, using the present formulation, the cancellation is possible by finding the corresponding expression for the correction term P ex (1) , to which end we need to consider the second-order contribution to the momentum equation. Applying the solvability condition (43) and (44) to the second-order Equation (42), we obtain, The latter is transformed by virtue of (52), The last (anomalous) term is canceled out by choosing, Combining the result (57) with the zeroth-order (equilibrium) value, we arrive at the extended pressure tensor Since the anomalous contribution is a diagonal tensor, cf. Equation (30), the result (58) is implemented with the extended equilibrium in the product-form by choosing the normalized (at unit density) diagonal elements of the pressure tensor as follows, which is equivalent to (32). Finally, combining the first-and second-order contributions to the density and the momentum equation, (45), (46), (54) and (56), using a notation, t , and also taking into account the cancellation of the anomalous term in (56), we arrive at the continuity and the flow equations as follows, where p is the pressure of ideal gas at constant temperature T, Π is the viscous pressure tensor, with S the rate of strain, and µ the dynamic viscosity, The above considerations can be summarized as follows-because of the third-order moment anomaly (30), the LBGK Equation (1) with the product-form equilibrium (25) is restricted in several ways, namely: (i) The temperature is restricted to a single value, the lattice reference temperature T L (31); (ii) The flow velocity has to be asymptotically vanishing; (iii) Stretched velocities amplify these restrictions by making it impossible to cancel even the linear (in velocity) anomaly in all the directions simultaneously.
Note that, in addition to all of the restrictions above, when using the conventional second-order equilibrium obtained by retaining the terms up to the order of ∼ u α u β in (25), the anomaly becomes not only confined to the diagonal elements Q eq ααα but also contaminates the off-diagonal elements Q eq αββ . While the diagonal anomaly (30) is genuine, that is, it is caused by the geometry of the discrete velocities, this additional off-diagonal deviation is due to an unsolicited second-order truncation of the product-form equilibrium (25).
The proposed revision of the LBGK model is based on extending the product-form equilibrium such that the anomaly of the diagonal third-order moment is compensated in the hydrodynamic limit by counter terms, which are added to the diagonal of the equilibrium pressure tensor. With this, all three restrictions mentioned above are addressed at once, without making a special distinction between the temperature, flow velocity or stretching as separate causes for the anisotropy.

Numerical Results
In this section, we shall access the accuracy and performance of the proposed LB model in a variety of scenarios of activating spurious anisotropy. First, we test Galilean invariance, isotropy and temperature independence of the model with both regular and rectangular lattices in the simulation of a decaying shear wave. Second, we validate the model for the more complex case of decaying homogeneous isotropic turbulence and show the effectiveness of using higher temperatures in saving compute time. Third, we investigate the applicability of the proposed model with stretched lattices in a periodic double shear layer flow, in a laminar flow over a flat plate, and finally in the case of the turbulent channel flow. In the simulations below, the gas constant was set to R = 1, the time step is δt = 1 and Grad's approximation, as proposed in [35], was used for wall boundary conditions.

Galilean Invariance, Isotropy and Temperature Independence Test
To probe the Galilean invariance and temperature independence of the model, the kinematic viscosity ν = µ/ρ (4) is measured for the decay of a plane shear wave. First, we consider the axis-aligned setup, with the initial condition, where Ma = u 0 / √ T is the advection Mach number, a 0 = 0.001 is the amplitude, L y = 200 is number of grid nodes in the y direction, ρ 0 = 1. The nominal kinematic viscosity is set to ν = 0.01. Periodic boundary conditions are imposed in both xand ydirections. The numerical viscosity (ν num ) is measured by fitting an exponential to the time decay of maximum flow velocity u x,max ∼ e −νt(2π/L y ) 2 . In this special case, the diagonal anomaly (30) is dormant and does not trigger any spurious effects because the derivatives ∂ xQxxx and ∂ yQyyy both vanish. Consequently, the extended equilibrium (33) becomes equivalent to the product-form equilibrium (25) in this case.
In order to compare with the standard LBGK, the standard velocities λ α = 1 were used in this simulation. Figures 1 and 2 show the importance of using the product-form equilibrium (25) as opposed to the conventional LBGK model with the second-order equi-librium. A strong dependence of the viscosity on the reference frame for the second-order equilibrium can be seen in Figure 1, where the viscosity drops with increasing advection Mach number. This well-known artifact of the second-order equilibrium is due to the nonvanishing anomaly in the off-diagonal moments Q eq αββ and, unlike the diagonal anomaly, is caused only by the approximate treatment of the product-form equilibrium. Moreover, as shown in Figure 2, even at a small enough velocity this spurious feature improves only at the lattice reference temperature T L . In contrast, as is shown in Figures 1 and 2, the product-form equilibrium of the present model is able to accurately predict the viscosity in this setup for a wide range of temperatures and reference frame velocities. Next, in order to trigger the anisotropy of the deviation terms (30) and to show the necessity of using the extended equilibrium, the shear wave is rotated by π/4. The anisotropy is further increased by also conducting simulations on a stretched grid with λ x = 2. The temperature is kept at T = 1/3. The viscosity measurement is shown in Figure 3 for different advection Mach numbers and stretching factors. It can be observed that the model lacks Galilean invariance for larger velocities when using the product-form equilibrium without correction (25). Furthermore, the stretching factor λ x = 2 results in a significant hyper-viscosity since the deviation (30) in this case amounts to a large positive number. However, once the correction term is included and the extended equilibrium (33) is used, the present model recovers the imposed viscosity, independent of the frame velocity and stretching factor.

Decaying Homogeneous Isotropic Turbulence
In order to further validate the model as a reliable method for the simulation of complex flows and to show the application of using higher temperatures, decaying homogeneous isotropic turbulence was considered. The initial condition, in a box of the size L × L × L, was set at 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 A is the parameter that controls the initial kinetic energy [36]. The initial velocity field is generated using a kinematic simulation as proposed in [37]. The turbulent Mach number is defined as where c s = √ T is the speed of sound. The Reynolds number is based on the Taylor microscale, and is given by where u rms = √ u · u/3 is the root mean square (rms) of the velocity and overbar denotes the volume average over the entire computational domain.
Simulations were performed at Ma t = 0.1, Re Λ = 72, κ 0 = 16π/L, at two different temperatures, T = 1/3 and T = 0.55, and with L = 256 grid points. Figure 4 shows a snapshot of the velocity magnitude √ u · u at time t * = t/τ = 1.0, where τ = L I /u rms,0 is the eddy turnover time, which is defined based on the initial rms of the velocity and the integral length scale L I = √ 2π/κ 0 . To quantitatively assess the accuracy of the model at different temperatures, the time evolution of the turbulent kinetic energy, normalized with its initial value (K 0 ), and of the Taylor microscale Reynolds number are compared in Figures 5 and 6 with results from direct numerical simulations (DNS) [36]. It is apparent that the two working temperatures yields almost identical results that agree well with the DNS simulation. This indicates that the correction terms do not degrade the accuracy of the model at higher temperatures, even though the magnitude of error term (30) is higher due to amplification of the linear term. The immediate advantage of using the present model at a temperature higher than the lattice temperature T L = 1/3 is that it effectively increases the characteristic velocity (here u rms,0 ) and therefore the time step by a factor of √ T/T L . A larger time step is equivalent to fewer number of time steps. The present model, therefore, speeds up the simulation by a factor of √ T/T L compared to the conventional LBM, which can operate only at the lattice temperature T L . Furthermore, this speedup strategy can be used for both steady and unsteady flows. This is in contrast to the preconditioned LBM [33], which works by altering the effective Mach number and therefore reduces the disparity between the speeds of the acoustic wave propagation and the waves propagating with the fluid velocity, cf. [33]. This makes preconditioned LBM restricted to steady state applications. In contrast, the present model enables us to increase the speed of sound without changing the Mach number. This increases the effective time step of the solver. Therefore, the present model increases the computational efficiency by decreasing the number of required time steps. Note that the theoretical temperature range of the model (like any other models based on the D1Q3 lattice) is 0 ≤ T ≤ 1, beyond that the populations become negative and the model is unstable. Therefore, while small temperature is possible but not beneficial, large temperature greater than 1 is out of the stability domain

Periodic Double Shear Layer
The next validation case to test the accuracy of the proposed model with the stretched lattice is the periodic double shear layer flow with the initial condition, where L is the domain length in both x and y directions, u 0 = 0.1 is characteristic velocity, δ = 0.05 is a perturbation of the y-velocity and α = 80 controls the width of the shear layer. The Reynolds number is set to Re = u 0 L/ν = 10 4 and the temperature is T = 1/3. Figure 7 shows the vorticity field at non-dimensional time t * = tu 0 /L = 1 using the conventional square lattice λ x = λ y = 1 and the rectangular lattice with λ x = 2, λ y = 1. Both lattice models perform qualitatively the same.

Laminar Boundary Layer over a Flat Plate
The next test case validates our model for wall-bounded flows. We consider the laminar flow over a flat plate with an incoming Mach number Ma ∞ = u ∞ / √ T ∞ = 0.1, temperature T ∞ = 1/3 and Reynolds number Re = ρ ∞ u ∞ L/µ = 4000, where L is the length of flat plate. Since the flow gradients in the transverse y-direction are much larger compared to the gradients in the streamwise x-direction, the mesh can be stretched in x-direction without significantly affecting the accuracy of the results. The computational domain was set to [L x × L y ] = [200 × 200] and a rectangular lattice with λ x = 2 was used. The flat plate starts at a distance of L x /4 from the inlet and symmetry boundary conditions were imposed at 0 ≤ x ≤ L x /4. In Figure 9, the horizontal velocity profile at the end of the plate is compared with the results of a regular lattice and with the Blasius similarity solution, where η is the dimensionless coordinate [38], It can be seen that results for the regular and the rectangular lattice nearly coincide and agree well with the Blasius solution. Thus, the model achieves accurate results with half of grid points compared to the regular lattice. Furthermore, the distribution of skin friction coefficient over the plate, with the wall shear stress τ wall = µ( ∂u ∂y ) y=0 , is shown in Figure 10 in comparison with the analytical solution C f = 0.664/ √ Re x , where Re x = u ∞ x/ν [38]. Also here, the results of the model with the regular and the stretched velocities are almost identical and in good agreement with the analytical solution.

Turbulent Channel Flow
In the final test case, we assess the accuracy and performance of the extended LBM for the turbulent flow in a rectangular channel, for which many numerical [39][40][41] and experimental [42,43] results are available. The channel geometry was chosen as [5.6H × 2H × 2H], where H is the channel half-width. The friction Reynolds number, based on the friction velocity u τ = τ w /ρ, was set to Re τ = 180. The initial friction velocity was estimated by where K = 0.41 is the von Kármán constant and u 0 = 0.1 is the mean center-line velocity. Periodic boundary conditions were imposed in the streamwise x-direction and the spanwise z-direction. The flow was driven by a constant body force in the x-direction, In order to accelerate the transition to turbulence, a non-uniform divergence-free forcing field as proposed in [44] was added to the flow for some period of time, until t * = tH/u τ = 5. Similar to the previous test case, grid stretching in x-direction with λ x = 1.4 was used in order to reduce the number of grid points in that direction while the temperature was set to T = 0.55, same as in Section 3.2 A snapshot of the velocity magnitude √ u · u is shown in Figure 11. Quantitatively, we compare the mean velocity profile with the DNS results of [40] in Figure 12. In wall units, the mean velocity is given by u + =ū/u τ and the spatial coordinate is y + = yu τ /ν. The statistics are collected after 30 eddy turnover times, i.e., after t * = 30. It is apparent that the viscous sublayer (y + < 5), the buffer layer (5 < y + < 30) and the log-law region (y + > 30) are captured well with our model and the mean velocity profile agrees well with that of the DNS. For a more thorough analysis, we compare the root mean square of the velocity fluctuations with the DNS data in Figure 13. Here, u x,rms = u x u x and u y,rms and u z,rms are defined in a similar way. It can be seen that the results are in excellent agreement with the DNS results [40]. This demonstrates that the LBGK model, also in the presence of a severe anisotropy triggered by stretched velocities, can be used for the simulation of high Reynolds number wall-bounded flows once the corrections are incorporated with the extended equilibrium.

Conclusions
While even with the standard discrete speeds (5) it is possible to develop an error-free, fully Galilean invariant kinetic model in the co-moving reference frame, it does require off-lattice particles' velocities [45,46]. Sticking with the fixed, lattice-conform velocities (6), one is faced with an inevitable and persistent error, which spoils the hydrodynamic equations whenever the flow velocity is increased or the temperature deviates from the lattice reference value, or the discrete speeds are stretched differently in different directions. We proposed an upgrade of the LBGK model to enlarge its operation domain in terms of velocity, temperature and grid stretching by suggesting an extended equilibrium. The extended equilibrium is realized through a product-form, which allows us to compensate the diagonal third-order moment anomaly in the hydrodynamic limit by adding consistent correction terms to the diagonal elements of the second-order moment. As a result, the extended LBGK model restores Galilean invariance and temperature independence in a sufficiently wide range, and can also be used with rectangular lattices. Similar to previous proposals [4][5][6]10,12], the relaxation term of the present model remains almost local as it uses only nearest-neighbor information for computation of the first-order derivatives in the extended equilibrium populations. The extended LBGK model was validated in a range of benchmark problems, probing different aspects of anomaly triggered either by increased velocity or temperature deviation from the lattice reference temperature, or by grid stretching. In all cases, the extended LBGK model featured excellent performance and accuracy in both two and three dimensions. In particular, the simulation of homogeneous isotropic turbulence demonstrated the expected speed-up when a higher temperature was used, while simulations of the laminar boundary layer and of the turbulent channel flow using stretched grids demonstrated good accuracy with a reduced number of grid points.
Furthermore, the present model can be extended to other applications including but not limited to high-speed compressible flows, which can be achieved by incorporating another solver for the total energy (see, e.g., the models proposed in [10,12]). Advanced collision models, such as multiple relaxation times (MRT) schemes, can also readily be employed in the present approach, which can be beneficial when running under-resolved simulations. These two avenues shall be subject of further development and application of the extended LBM. Acknowledgments: We thank anonymous referee who suggested comparison with the local correction method of Ref. [9]. Authors are grateful to Seyed Ali Hosseini at ETHZ for computing spectral dissipation presented in Figure A1.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: LBM Lattice Boltzmann method LBGK Lattice Bhatnagar-Gross-Krook

Appendix A. Comparison of Extended LBGK to Locally Corrected LBM
Below, we compare the locally corrected lattice Boltzmann model (LC LBM) [9] with both the standard and the present extended LBGK. To that end, it suffices to consider the one-dimensional D1Q3 lattice. In order to introduce the LC LBM, we begin with the standard LBGK (δt = 1, R = 1), The equilibrium populations in (A1) are given by (24) and (25), Thanks to the diagonal anomaly, the second-order asymptotic analysis of Section 2.4 results in the following viscous stress in the one-dimensional version of the Navier-Stokes Equation (61), For the LBGK model, the anomalous (second) term in (A3) reads, Upon realizing that the first term of the anomalous contribution (A4) is similar in its structure to the relevant (first) term in the LBGK stress (A3), the locally corrected (LC) LBGK groups these two terms together and replaces the relaxation parameter ω in (A1) with a new relaxation ω LC , which depends on the flow velocity. While the original work [9] addressed the case of the lattice temperature, T = T L = 1/3 (31), we first consider a slightly more general formulation for a flexible temperature parameter. Consequently, the locally corrected relaxation parameter ω LC in the LBGK Equation (A1) is defined as, where the renormalization factor X reads, The LBGK model with the locally corrected relaxation parameter ω LC (A5) results in the viscous stress of the form (A3), with the remaining error term, For the sake of a discussion, let us introduce the local Mach number, Ma x = u x / √ T. For a quasi-incompressible (slow) flow, the density variation scales as ∂ x ρ ∼ Ma 2 x . Thus, for T = T L , the error (A7) can be estimated as,Π LC xx ∼ Ma 3 x . This is two orders of magnitude lower than the error of the original LBGK at T = T L , cf. Equation (A4), at small Mach number. Moreover, by setting the temperature T = T L = 1/3, it was first realized in Ref. [9] that the error (A7) reduces to, In this case, the scaling at Ma x → 0 becomes,Π LC xx ∼ Ma 5 x . In other words, the local correction at T = T L provides a gain of two orders of magnitude in accuracy with respect to the standard LBGK under the quasi-incompressible flow conditions [9]. This consideration extends straightforwardly to the D2Q9 and D3Q27 lattices by constructing a multiple relaxation time LBM that corrects the relaxation of each diagonal component of the pressure tensor [9].
However, for a generic isothermal flow, the error (A7) becomes amplified through the renormalization factor (A6) as the velocity increases and eventually diverges when u 2 x → (1 − T)/3. This error persists also for the special case T = T L (A8). On the other hand, the second-order analysis of Section 2.4 reveals that the present LBGK with the extended equilibrium (33) removes the entire anomalous term,Π ex xx = 0. Thus, the difference between the extended LBGK and the LC LBM [9] is expected beyond the asymptotic Ma x → 0.
In order to demonstrate this point, a spectral analysis was performed for the twodimensional D2Q9 lattice (see [21] for details of the spectral analysis in the LBM context). The normalized spectral dissipation of acoustic modes (ω κ )/νκ 2 x , is shown in Figure A1, for T = T L and the background flow velocity (u x , u y ) = (0.3, 0), for the three models, the standard LBGK, the present extended LBGK and the LC LBM of Ref. [9].
It can be seen that, the extended LBGK recovers the correct dissipation rate in the continuum limit (vanishing wave number κ x ), confirming its Galilean invariance. However, both the standard LBGK and the LC LBM show deviations in the form of under-dissipation at low wave numbers, while the deviation for the LC LBM is indeed smaller. This nonvanishing deviation is amplified in cases with different working temperature and/or non-unit stretching factor for both the standard LBGK and the LC LBM, which makes their applications limited to the quasi-incompressible flow regime at the lattice temperature. Finally, it is interesting to note that, in the case of shock capturing, all the three models are expected to behave similarly, given that their respective dissipation rates at the wave number κ x = π are close in value, see Figure A1. This observation is confirmed by the simulation of a shock tube with the following initial condition, (ρ, u x , T) = (ρ l , 0, 1/3), x ≤ L/2, (ρ r , 0, 1/3), x > L/2, with L = 800 grid points and viscosity ν = 0.04. Results are presented in Figures A2 and A3 for ρ l = 1.5, ρ r = 0.5, corresponding to the initial density ratio ρ l /ρ r = 3. Figures A2 and A3 demonstrate that all models produce almost indistinguishable results, with a similar oscillation pattern at the shock front.