Low Dissipative Entropic Lattice Boltzmann Method

: In the entropic lattice Boltzmann approach, the stability properties are governed by the parameter α , which in turn affects the viscosity of a ﬂow. The variation of this parameter allows one to guarantee the fulﬁllment of the discrete H -theorem for all spatial nodes. In the ideal case, the alteration of α from its normal value in the conventional lattice Boltzmann method ( α = 2) should be as small as possible. In the present work, the problem of the evaluation of α securing the H -theorem and having an average value close to α = 2 is addressed. The main idea is to approximate the H function by a quadratic function on the parameter α around α = 2. The entropy balance requirement leads to a closed form expression for α depending on the values of the H -function and its derivatives. To validate the proposed method, several benchmark problems are considered: the Sod shock tube, the propagation of shear, acoustic waves, and doubly shear layer. It is demonstrated that the obtained formula for α yields solutions that show very small excessive dissipation. The simulation results are also compared with the essentially entropic and Zhao–Yong lattice Boltzmann approaches.


Introduction
The lattice Boltzmann (LB) method has established itself as an alternative approach to the Navier-Stokes equations for the modeling of fluid motion. Due to its simplicity and high potential for parallel computing, the LB approach can be efficiently implemented in numerical solvers aimed at addressing academic and application-oriented problems [1][2][3][4]. Except for hydrodynamic problems, the LB approach has found applications in the modeling of porous media [5,6], microflows [7], multiphase, and other complex flows [8][9][10][11].
The standard realization of LB models is based on local collision and linear advection to adjacent lattice nodes in the direction of the corresponding discrete velocities. Such an approach is rather appealing due to its simplicity but is prone to instabilities when Mach (Ma) and Reynolds (Re) numbers are increased. Hence, additional stabilization techniques are usually adopted. These can be divided into several classes.
The instabilities usually manifest in oscillations, unphysical (negative) values of the hydrodynamic variables (density), and eventually in the breakdown of the solutions. It is possible to enhance the dissipative properties (and therefore stabilize the computations) of LB models by tuning the bulk viscosity [12]. This idea is generalized in the multiplerelaxation-time (MRT) LB method [13][14][15]. The drawback of this method is that acoustic waves show over-dissipative behavior due to the increased bulk viscosity. To maintain correct dissipation rates of acoustic waves, one can adopt selective viscosity filtering [16], which suppresses only high wave-number instabilities. In this approach, viscosity becomes dependent on the values of the wave-number. A simpler way of stabilizing the simulations is enforcement of the positivity of the solutions by adapting the relaxation time in the LB models [17]. This stabilization technique is rather straightforward but is less accurate than more sophisticated methods, such as the entropic LB (ELB) method [18]. In the ELB method, the stabilization is achieved by tuning the relaxation rates in order to satisfy the H-theorem. The ELB approach requires two components [19,20]: the equilibrium state, which minimizes the H-function (entropy taken with minus sign), and the adjustable (in a special way) relaxation rate (the parameter α), which is evaluated from the significantly nonlinear entropy balance equation. Note that the standard equilibrium in the form of the polynomial function on the flow velocity is not suitable, since for this case the H-theorem is not valid [21,22]. The problem of construction of the entropic equilibrium states is non-trivial [23][24][25][26]. The solution α to the entropy balance equation can be obtained by applying the Newton-Raphson method, the bisection methods [27][28][29], or a combination of the Newton-Raphson method and analytical estimates [30]. In the essentially entropic approach, the solution to the entropy balance equation, which is reformulated as an inequality, has been obtained in an exact form by exploiting Padé approximations for the logarithms entering the entropy function [31]. This method yields solutions for α that have little deviation from its normal value (α = 2). The entropic approach of Zhao and Yong relies on geometric considerations [32] (secant algorithm) and gives concise closed-form expressions for α. The comparison of several entropic methods has been performed recently [33]. Another class of entropic methods is based on the entropic limiters. In this approach, the non-equilibrium part of the distribution function is multiplied by a coefficient that depends on local non-equilibrium entropy [34][35][36]. This type of limiter has been successfully used for the modeling of supersonic flows with the LB method [37,38]. The entropic method can be extended for multiple relaxation time LB models. This approach uses decomposition of the distribution function into equilibrium, non-equilibrium hydrodynamic, and non-equilibrium non-hydrodynamic components. Non-equilibrium parts of the distribution function have different relaxation rates [39][40][41][42]. In order to satisfy the H-theorem, only the relaxation times for high-order moments are varied; hence, the viscosity remains unchanged.
In the regularized LB approach, the distribution function is projected onto some finite basis during each time step. The removal of the moments that cannot be described by the model leads to noticeable improvements in stability [43][44][45]. Furthermore, the concept of recursive regularization has been introduced, including thermal flows [46][47][48][49][50][51][52]. In the recursively regularized LB models, several non-equilibrium high-order moments are kept. Their values can be computed analytically as a special combination of low-order moments. In addition, variants of the regularized LB models in the relative reference frame have been proposed [53][54][55]. Note that the entropic LB and regularization methods share common features: regularized LB models maximize quadratic entropy [56]. In practice, however, the regularized LB models are usually more stable than un-regularized. Linear stability analysis shows that the regularization significantly changes the dissipation and dispersion properties of LB models [57][58][59][60]. Even more so, the regularized LB models can be unstable under some specific conditions for which un-regularized models are stable.
Another stabilization approach is focused on modification of the advection step. Employment of the standard streaming step means that the Courant number is equal to unity, and this can cause instabilities. Then, in order to reduce the Courant number, one can apply a generalized streaming step, for which the values of the distribution function in additional nodes are needed [61][62][63][64][65]. Numerical experiments indicate that the maximal stable Reynolds number for this method depends exponentially on the reduction of the Courant number. On the other hand, this method is more complicated than the standard one-for instance, the formulation of the boundary conditions requires additional efforts [66][67][68] The aim of the present paper is to develop a variant of the entropic LB method in which the averaged deviations of the parameter α from its normal value, which equals 2, are minimized. This in turn means that the variations of the viscosity are also small and the dissipation properties of the modeled system are not distorted. This paper is organized as follows. In Section 2, the main features of the entropic LB approach are considered. The formulation of the variational problem leading to the entropic equilibrium compatible with the hydrodynamic equations is performed by applying geometrical treatments; in addition, some pitfalls of the entropic method are discussed. In Section 3, the derivation of a novel formula for α is presented. Since the entropy balance equation is a nonlinear transcendental equation, one cannot find an exact analytical solution α, but it is possible to obtain an approximate solution. This is performed in several steps. First, one finds an interval that contains the solution α; the lower (dissipative) bound and the upper bound of this interval are evaluated by applying Padé approximations for the logarithms. Then, the quadratic approximation of the H-function for α in the interval between the lower and upper bounds is constructed, and the entropy balance equation reduces to a quadratic algebraic equation for the variable α, which is solved analytically. In Section 4, the results of the numerical simulations are addressed: the shock tube problem, the propagation of shear and acoustic waves, and the doubly shear layer problem are solved by employing the essentially entropic, Zhao-Yong, and present method. In Section 5, the major results of the paper are outlined.

Entropic Lattice Boltzmann Method
In the present paragraph, the main ideas of the entropic lattice Boltzmann (ELB) approach are briefly recalled. The LB equation reads as follows [20]: where f i are the distribution functions corresponding to the particle velocities c i , i = 1 . . . N, f eq i are the local equilibrium distribution functions (their exact form is discussed below), t, x are time and spatial coordinates, and ∆t is the time step. The variable 0 < β < 1 is computed as where τ is the relaxation time. In the standard LB method, the variable α = 2 (α is also called the path length), τ is evaluated as ν = c 2 s τ, c s is sound speed, and ν is viscosity. In the ELB method, the variable α is not constant and may change from node to node, and ν = (2/α)c 2 s τ, while α is the solution of the following equation: where H[ f ] is a smooth convex function, f = { f i , i = 1 . . . N}, and f eq = { f eq i , i = 1 . . . N}. The form of the condition (2) is motivated by the following inequality: where the convexity property for the function H was used, and the condition (2) was applied. Hence, the condition (2) guarantees the fulfillment of the H-theorem, and the entropy is increased by the collisions (the H-function is decreasing). Note that the explicit form of the function H[ f ] is not needed in the derivation of the H-theorem. The usual form of the H-function in the LB theory is Boltzmann-like, that is, where w i > 0 are the lattice weights (discrete analogs of the absolute Maxwell distribution). One should emphasize that the ELB approach is based on two components: the equilibrium f eq i (entropic equilibrium), which is a minimum of the H function, and a parameter α, which is tuned in order to guarantee a discrete H-theorem. The presence of both components is crucial. The entropic equilibrium is a non-polynomial function of the flow velocity. For the one-dimensional three-velocity model (D1Q3), the following entropic equilibrium can be used [69]: where f ±1 , f 0 correspond to the particles with the velocities ±c, 0; ρ, u are the density and the flow velocity; and c s = c/ √ 3 = 1/ √ 3, since in the present study c = 1. The extension of the equilibrium (3)-(4) for the two-dimensional nine-velocity model (D2Q9) is based on the product of two equilibrium functions for the D1Q3 model.
The variation of the path length variable α controls the dissipation properties of the flow via increase or decrease of the viscosity. Equation (2) can have solutions α < 2; this means that the stabilization is needed and the viscosity of the flow is increased. In the case α > 2, the viscosity is decreased. In the perfect case, the alteration of α from 2 should be as small as possible since this guarantees that the distortion of the viscosity of the modeled flow is minimal.
Following the previous studies, it would be convenient to introduce a scalar product as follows: Then, the macroscopic variables-density ρ, momentum ρu-can be computed as where γ denotes the spatial coordinates x, y, z.
Note that in terms of the introduced scalar products, the H-function can rewritten as follows: The formulation of the problem for finding f eq i has some subtleties. It would be instructive to consider this problem in detail. For athermal flows, the hydrodynamic fields ρ and u are conserved during the collisions, but the energy ρc 2 s + ρu 2 = ∑ i f i c 2 i is not a collision invariant. Then, D + 1 quantities remain unchanged due to collisions, where D is the number of spatial dimensions. For a given ρ and u, the discrete distribution functions f i , i = 1 . . . N define a manifold denoted by S in R N , an the dimension of S equals N − (D + 1). The entropic local equilibrium state f eq i is then a point in S for which a minimum of H[ f ] is attained. One has the following: It is important to mention that it is not correct to replace f by f eq in the constraints in Equation (5). The geometrical interpretation of the problem is presented in Figure 1. Consider a point f in S, and let f ( ) be a curve in S parameterized by ≥ 0 such that f (0) = f . By applying the Taylor expansion one obtains the following: where δ f is a vector in the tangent space defined as T f (S) to the manifold S at the point f .
On the other hand, Remember that the manifold S is specified by the conditions G j = 0, j = 1 . . . N; then, ∇G j are normal vectors to S, and they can be considered as a basis in T f (S) ⊥ . Therefore, where λ j are expansion coefficients, which can be recognized as Lagrange multipliers. Equation (6) is solved jointly with the constraints G j = 0. Not all equilibrium states are allowed from the physical point of view. Obviously, the considered model should lead to correct hydrodynamics. In order to recover Navier-Stokes equations (at least with O(u 3 ) error), it is sufficient to include the following additional condition [70]: It is worth emphasizing that one cannot replace f eq by f in Equation (7), since the condition (7) holds only for the equilibrium part of the distribution function. The equilibrium distribution function, which also respects the condition (7), is termed perfect.
The form of the solution for the problem (5), (7) depends on the particular choice of H[ f ]. As has been mentioned before, the Boltzmann H[ f ] is the most popular vari-ant. However, one may consider a simpler variant-quadratic entropy, which takes the following form: This function plays a role in the regularization of LB models, namely, regularized LB models minimize the quadratic H-function [56]. It is tempting to apply quadratic entropy in the ELB approach. Firstly, consider the Equation (2) in the form of equality. The solution is as follows [69]: Hence, one has a closed-form expression for α that exactly solves the entropy balance equation. Next, one needs to find the equilibrium states from the problem (5). From (5) and (6), one gets the following: The obtained equilibrium depends linearly on the flow velocity u; hence, it cannot be perfect. In addition, as was mentioned in the previous studies [69], this equilibrium admits negative values, and this can lead to negative density. As a final remark, it can be argued that some caution is necessary in the selection of the H-function and the local equilibrium in the ELB method. In the present study, the equilibrium (3)-(4) is used, which is entropic (minimizes Boltzmann H-function) and perfect.

Low-Dissipative ELB Method
For the Boltzmann entropy, the Equation (2) is rewritten as [31,71] where the variable x measures the departure from the local equilibrium It would be convenient to introduce the following scalar products: where a, b are two vectors with the components a i , b i , i = 1 . . . N. The general idea for solving Equation (2) or (8) is as follows. At the first step, the bounds for the minimal values of α denoted as α (α ≤ 2) are evaluated from (8). At the second step, the upper bounds for α denoted as α (α ≥ 2) are computed. Next, the function H[ f + α( f eq − f )] is approximated by a quadratic function of the variable α − 2 in the interval [α, 2], and a similar approximation is performed in the interval [2, α]. Finally, the estimates for α ∈ [α, α] (termed low-dissipative) are evaluated by solving a quadratic equation. The geometric interpretation of the proposed method is presented in Figure 2.

Bounds for the Logarithm
To obtain the values of α, α, one needs to find lower and upper bounds of the logarithm. The most natural way of doing this is to apply the Padé approximations, which are based on rational functions. Such an approach has been used in the essentially entropic LB method [31,71]. Following Topsøe [72], one has a chain of inequalities for z > 0 as follows: where the first lower approximants φ n (z) in (9) are given by and the first upper approximants ψ n (z) are In order to derive the bounds for z ∈ (−1, 0], the duality between functions in the intervals z ∈ (−1, 0] and z ∈ [0, ∞) is introduced. For the function g(z), the dual function g * (z) is defined by the rule . Note that the dual function for log(1 + z) is log(1 + z) (i.e., this function is self-dual). Therefore, the following inequalities are obtained for −1 < z ≤ 0 by applying the duality properties to (9) In addition, the following inequality is valid for any z > −1: where In the present study, it would be convenient to apply the following upper bound for z ∈ (−1, 0] instead of log(1 + z) ≤ φ * n (z):

Evaluation of α
The value of α is obtained as a particular solution of (8). For the present method, it serves as a lower bound. Note that α can be used as a stabilizer for LB models, but in practice, it leads to dissipative behavior.
Consider the Equation (8). Since α is not inside the logarithmic function in the right hand side of (8), approximation of only the left-hand side is needed. Assume that one finds a function A(α) such that for α ≥ 1; then, the fulfillment of A(α) ≤ α( f , x log(1 + x)) means that the inequality (8) is satisfied.
One has where the first term in this sum is approximated with use of the inequality from (9).
For the second term, the inequality (12) for M = 3 is applied.
Hence, A(α) is a sum of (13) and (14). Using ( f , x) = 0, the inequality A(α) ≤ α( f , x log(1 + x)) is rewritten as Now, one can prove the following: Proof. Note that from (11), it follows that ( f , x log(1 + x)) ≤ ( f , x 2 ). In addition, since Applying Lemma 1, one has for the left-hand side of (15) Finally, redefining the function A(α) as α 1 2 ( f , As a remark, one can notice that it is possible to take M > 3 in the inequality (12) applied in (14), but the numerical simulations show that the difference will be relatively small in comparison to Equation (16).

Evaluation of α
From the condition f i + α( f eq i − f i ) ≥ 0, i = 1 . . . N, one obtains the upper bound α * [28] for the variable α as follows: In the present paragraph, an additional upper bound defined as α is evaluated. Assume that one finds a function B(α) such that for all α ≥ 1. Then, the solution to the equation B(α) = α( f , x log(1 + x)) defines the upper bound. One has Note that ψ * 1 = ψ 1 has been used. Then,

Proposition 1. The inequality (8), that is,
has the following solution.

Numerical Experiments
In the present paragraph, the numerical solutions of several test problems are considered: the Sod shock tube, the propagation of shear waves and acoustic waves for different angles between a flow direction vector and a perturbation wave-vector, and the double shear layer. The results are obtained with use of the following entropic LB methods: essentially entropic (EELB) of the first order (used only for Sod shock tube) and the third order [31], the method of Zhao and Yong (ZY) [32], the lower bound α evaluated from (16) (used only for the Sod shock tube), and the method described in Proposition 1, defined hereafter as the low-dissipative (LD). During the simulations, the following variables are recorded: the mean, minimal value of α measured over the modeled spatial domain at some moment of time t and L 1 , L 2 deviations of α from 2 in the form where x denotes spatial nodes, n is the number of the spatial nodes, and, for all modeled problems, ∆t = 1 is used. For convenience, the list of main notations is also given in Appendix A, Table A1.

Shock Tube
The first test problem is the athermal one-dimensional Sod shock tube. The initial condition is a step density profile: ρ = 1.5, x ≤ H/2, ρ = 0.5, x > H/2 (i.e., the initial density ratio is 3:1), and H is the length of the domain. For this problem, a grid consisting of 500 spatial nodes is used, and the viscosity is ν = 10 −5 . Since the problem is onedimensional, the D1Q3 LB model is applied. The equilibrium is taken in the form (3)-(4).
The simulation results for the different entropic LB methods are presented in Figure 3 at the time step t = 250. As expected, the density profiles show oscillations [73], and the amplitudes of the oscillations are larger for the less dissipative methods. The first-order EELB method leads to the smallest average values of α (and potentially to the largest dissipation), followed by the method based on α from (16), and the method of Zhao and Yong. LD and the third-order EELB are the least dissipative. Moreover, the LD method has the closest to 2 average value of α, that is, it produces the smallest amount of additional dissipation compared to the other methods, but has the largest variations of D 1 , D 2 , since relatively large values of α > 2 are observed ( Table 1). Remember that the condition α > 2 means that the viscosity is decreased. Note that the dispersive oscillations, which are typical for the ELB method, can be removed by applying LB models with additional artificial viscosity and dispersion terms [74]. In addition, it would be interesting to compare the accuracy of the LB schemes with Navier-Stokes numerical solvers without artificial viscosity [75].
Note that the methods in which α cannot exceed 2 tend to have smaller values of D 1 , D 2 . Then, in order to test the models under the same conditions, the following experiment has been performed: for the third-order EELB and the present LD method, the variable α has been bounded by 2 in the simulations. In this case, these models show significantly smaller values of D 1 , D 2 variances ( Table 1). The methods based on α from (16) and the first-order EELB give the least accurate solutions to the entropy balance equations and serve as predictors for the LD and third-order EELB methods. In subsequent problems, only the latter approaches are addressed.   (16); LD is the low-dissipative entropic method described in the Proposition 1; EELB 1 is the first-order EELB method; EELB 3 is the third-order EELB method; ZY is the method of Zhao and Yong; LD, max(α) = 2 EELB 3, max(α) are the low-dissipative entropic and the third-order EELB methods for which α is bounded by 2.

Shear Waves
The simulation of shear waves serves as an excellent test for the assessment of the dissipative properties of the LB models and validation of the linear stability analysis predictions [58]. For this problem, a periodic two-dimensional spatial domain is considered, the initial conditions for the hydrodynamic fields are as follows: u y | t=0 = εc s Ma cos(ϕ(k)) cos(k x x + k y y), where ε = 10 −4 , Ma = 0.2, 0.4 is the flow Mach number, ν = 10 −5 , k = (k x , k y ) is the perturbation wave-number, and ϕ(k) = arctan(k y /k x ). In addition, the following variable is introduced: u 0,y = εc s Ma cos(ϕ(k)). It would be convenient to introduce the time variable in the form Fo = k 2 νt (the Fourier number). The first considered case for the problem (26)-(28) is as follows: ϕ = 0 and k y = 0, u x is constant, and u y depends only on x. The size of the spatial domain is 32 × 2 for this case. In addition, k x = 2π/8, this means that the perturbation wavelength equals 8 spacings between lattice nodes. The Navier-Stokes equations have the following solution for the velocity u y : u y = u 0,y exp(−νk 2 x t). The simulation results are shown in Figure 4 for Ma = 0.2 and Ma = 0.4 (left slides). The second case is k x = 2π/16, k y = 2π/12. The simulation spatial domain contains 48 × 36 nodes. Interestingly, this problem is unstable for the third-order recursively regularized LB model [58]. In the present study, all the entropic models are stable but have different dissipation properties. The results are depicted in Figure 4 (right slides).
The ratios ν e /ν (termed effective viscosity, where ν e is the viscosity measured during the numerical simulations, and ν is the theoretical viscosity) are shown in Table 2 for both considered cases. Clearly, from the results in Figure 4 and Table 2 one can see that the present LD method has the smallest dissipation among all considered entropic approaches and gives the values of ν e /ν closest to the ones for the conventional D2Q9 model (α = 2). It is important to mention that for the D2Q9 model, the viscosity decreases if the flow velocity is increased [76]. Such a behavior is stipulated by the cubic defects ( ∇ · (ρuuu)), which enter the Navier-Stokes equations derived from low-order LB models like D2Q9 [70,76,77]. This effect is seen in Tables 2 and 3. The magnitude of the effective viscosity tends to be smaller than unity for the D2Q9, α = 2 model. : k x = 2π/8, k y = 0; and the right slides: k x = 2π/16, k y = 2π/12. LD is the low-dissipative entropic method described in Proposition 1; EELB 3 is the third-order EELB method; ZY is the method of Zhao and Yong. The Navier-Stokes solutions are denoted as: (−). Table 2. Shear waves, the ratios ν e /ν (ν e is the viscosity measured during the numerical simulations, ν is the theoretical viscosity) are presented. α = 2 is the standard LB method; LD is the low-dissipative entropic method described in Proposition 1; EELB 3 is the third-order EELB method; ZY is the method of Zhao and Yong.

Acoustic Waves
This problem involves the variations of density. The initial conditions are given by the following expressions [58]: In the simulations, the following values of the parameters are adopted: ε = 10 −4 , Similarly to the previous case, the problem (29)-(31) is solved for two combinations of k x , k y . The spatial domain is considered as periodic, and the size of the domain is taken as 48 × 36. Initially, the case ϕ = 0 and k x = 2π/8, k y = 0 is considered. The results are plotted in Figure 5 (left slides). Next, the results of the simulations for the inclined acoustic waves k x = 2π/16, k y = 2π/12 are presented in Figure 5 (right slides). Again, the application of the LD method only slightly varies the dissipation properties compared to the standard LB model results (Table 3). Acoustic waves, upper slides: Ma = 0.2, lower slides: Ma = 0.4; left slides: k x = 2π/8, k y = 0 and the right slides: k x = 2π/16, k y = 2π/12. LD is the low-dissipative entropic method described in Proposition 1; EELB 3 is the third-order EELB method; ZY is the method of Zhao and Yong. The Navier-Stokes solutions are denoted as: (−).
In conclusion, one can notice that the considered entropic methods are significantly better suited for the modeling of shear and acoustic waves than the regularized LB models [58,60] due to the relatively small deviations of the measured efficient viscosities from the theoretical values.

Double Shear Layer
In this problem, a flow in a doubly periodic two-dimensional square domain of a size L × L is considered. The initial conditions for the flow velocities are defined as follows: where the parameters k, δ are as follows: The initial density equals 1, and the Reynolds number is Re = U 0 L/ν, where ν is the viscosity.
Note that the form of the initial distribution function noticeably affects the dissipation properties of the flow, so in order to avoid initial spurious oscillations, the distribution function at t = 0 is initialized in the regularized form f (eq) + f (1) , where the elements of the stress tensor entering the non-equilibrium component of the distribution function f (1) are computed by applying finite-difference approximations [41] of ∇u .
The numerical simulations of the problem (32)- (34) are performed for Re = 3 × 10 4 and U 0 = 0.04. This means that the Mach number equals Ma = U 0 /c s = √ 3 × 0.04 ≈ 0.0693. It is convenient to introduce the characteristic time variable (convection time) as t c = L/U 0 . The flow is assumed to be stable if the solutions remain finite in the time interval t ∈ [0, 2t c ].
The sensitivity of the entropic methods to the grid resolution can be tested for the present problem. For coarse meshes 64 × 64 and 128 × 128, the third-order EELB, Zhao-Yong, and LD methods show qualitatively correct results but produce spurious vortices. These vortices may lead to blow-up instabilities, but in the present case, the solutions are finite because the fulfillment of the entropy balance equation suppresses unstable oscillations. Note that the standard D2Q9 model (α = 2) is unstable for grid resolutions of 64 × 64 and 128 × 128. For grid resolutions larger than 128 × 128, spurious vortices are not observed, and all approaches produce accurate vortical profiles. In the case of the spatial resolution 256 × 256, the normalized vorticity profiles are presented for the all considered methods at t = t c in Figure 6. One can observe that all approaches give very similar results. Next, the time histories of the variables of α − 2 and D 1 are presented in Figure 7. Clearly, the LD method gives the smallest average values of |α − 2| among the all considered methods. In addition, the LD method shows smaller D 1 variances than the third-order EELB method.
Finally, the evolution of the following variables is considered: the averaged on the spatial variables kinetic energy u 2 = u 2 x + u 2 y normalized on U 2 0 and the averaged on the spatial variables squared vorticity ω 2 (enstrophy) normalized on (U 0 /L) 2 . As a benchmark, the simulation results based on the standard D2Q9 LB model for a grid with the resolution 1024 × 1024 are taken (these simulations are stable). From Figure 8, one can observe that all entropic methods show similar results: they are close to the benchmark solution but more dissipative to some extent than the benchmark solution, which uses well-resolved grid. In addition, the LD method is slightly less dissipative than the other entropic treatments. Figure 6. Double shear layer. The grid resolution is 256 × 256, Ma = 0.0693, Re = 3 × 10 4 . The normalized vorticity plots are presented at the moment of time t = t c ; LD is the low-dissipative entropic method described in Proposition 1; EELB 3 is the third-order EELB method; ZY is the method of Zhao and Yong.

Performance
The performance of the considered entropic stabilization approaches is affected by the implementation details and the modeled problems. For instance, the method of Zhao and Yong is applied only for the nodes in which the H-theorem is violated. Therefore, this method can be faster than the others if the proportion of the spatial nodes prone to instabilities is small. On the other hand, the Zhao and Yong and the present low-dissipative methods exploit the values of the H-function, which require the computation of logarithms; the latter is relatively slow operation. On average, the LB D2Q9 model stabilized by the considered entropic techniques is approximately two to three times slower than the LB D2Q9 model without additional stabilization. For example, for the shear waves, the Zhao and Yong method is 2.2 times slower than the conventional D2Q9 model followed by the low-disspative and essentially entropic methods, which are 2.38, 2.77 times slower than the conventional D2Q9 model. For acoustic waves, the methods of Zhao-Yong, low-dissipative, and essentially entropic are 1.81, 2.0, and 2.38 times slower than the conventional D2Q9 model. Note that, in the present study, the simulations are performed with use of selfstanding C++ code (and the figures are generated with use of MATLAB). In order to increase the speed of code execution, the dynamic arrays, such as vector class from the standard C++ library, are avoided in the parts of the code related to the entropic stabilizers.

Conclusions
In the present study, a novel approach for finding the solutions to the entropy balance Equation (2) is proposed: by employing the bounds for the logarithms, the H-function is approximated by quadratic function on the variable α around α = 2; then, the solution for α is obtained analytically from the quadratic algebraic equation. The proposed stabilization technique (termed as low-dissipative) is very "mild" (i.e., the departure of α from its normal value 2 is small). Therefore, the injection of the additional dissipation in the modeled systems is very moderate.
In order to validate the developed method, the following test problems have been considered: the shock tube problem, the propagation of shear waves, acoustic waves, and doubly shear layer problem. In addition, for comparison purposes, these problems have been also solved employing the approaches in which α is given in the form of a closed analytical expression: the essentially entropic LB method [31] and the method of Zhao and Yong [32]. For all problems, the presented method gives the solutions for α that are closest to its normal value 2. For the Sod shock tube problem, one can notice that all entropic methods lead to dispersive oscillations near shock. This behavior is typical for these models [73]. The entropic methods for which alpha is close to 2 (such as essentially entropic and the present approach) show larger oscillations than more dissipative methods. In addition, it has been shown that, compared to the regularized LB models, the present entropic method is able to reproduce correct attenuation rates for shear and acoustic waves for relatively coarse grids (8-16 nodes per wavelength).
Funding: This work is supported by the Ministry of Science and Higher Education of the Russian Federation, Project No. 075-15-2020-799.

Data Availability Statement:
The simulation code that supports the findings of this study is available from the author upon reasonable request.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

ELB
Entropic lattice Boltzmann EELB Essentially entropic lattice Boltzmann LB Lattice Boltzmann LD Low dissipative ZY Zhao Yong