In the last three decades, the lattice Boltzmann model (LBM) has provided a reliable mesoscopic computational method for fluid dynamics, with applications ranging from high Reynolds number flows to flows in porous media and relativistic hydrodynamics.^[1–9] We refer to Refs. [3]–[5] for reviews of the method and its applications. As a discrete space-time kinetic theory for hydrodynamics, the LBM employs discretized particle distribution functions associated with discrete velocities to describe the flow field.^[1,2] By fitting the discrete velocities into a regular lattice, the LBM realizes the propagation and collision of the distribution functions efficiently.

Though simple and efficient, the standard LBM is limited to moderate Reynolds number flows due to the lack of numerical stability.^[1,10,11] To alleviate this obstacle, the entropic LBM (ELBM) restoring the discrete H-theorem has been proposed in Refs. [12]–[21] as a paradigm shift for computational fluid dynamics. The introduction of the discrete H-theorem in the ELBM significantly extends the application range of the discrete kinetic theory to turbulent flows,^[22,23] binary droplet collisions,^[24] and multiphase fluid–solid interface problems.^[25] The solid physical background and successful applications make the ELBM a powerful approach for the study of complex flows.^[26–29]

A crucial problem in the ELBM is to determine a relaxation parameter α introduced in 1999^[13] for ensuring the discrete H-theorem. It needs to solve a complicated nonlinear algebraic equation at each grid point, which greatly affects the efficiency of the ELBM. Considerable efforts have been made in the past^{[14,15,18,30]} to improve the efficiency and it was only recognized in Ref. [31] that the value of α should not lead to a numerical entropy increase. This requirement essentially guarantees the discrete H-theorem at the numerical level and was further emphasized in Ref. [21] recently. Therein some analytical approximate expressions of α, which also guarantee the discrete H-theorem numerically, were derived by relaxing the entropy equality^[21] and by making a near-equilibrium assumption. Nevertheless, the first-order approximation is too dissipative while higher-order approximations seem difficult to be explicitly obtained. Therefore, a satisfactory determination of the relaxation parameter is not available yet and a critical breakthrough is much needed for using the ELBM.

In this work, we solve the long-standing problem by proposing a simple formula for the relaxation parameter α of the ELBM. This formula is based on a novel combination of the consistency with hydrodynamics and the constraint that the entropy must not increase within a discrete time step. It is free of any near-equilibrium assumption. Besides compliance of the discrete H-theorem at the numerical level, an excellent property of the formula is that it is applicable to arbitrary convex entropy functions. Additionally, numerical simulations demonstrate that the value of α given with the formula effectively characterizes the drastic changes of flow fields.

The entropic lattice Boltzmann model reads as^[16]

With the equal entropy auxiliary distribution f^* : = f + α(f^(eq) − f), the ELBM (1) can be rewritten as

In spite of this H-theorem, the entropy balance equation (4) requires an additional step of searching for the parameter α. The efficiency of this search is crucial for the ELBM.

In this section, we derive a simple formula for α. Observe that the discrete H-theorem (6) holds true for all f^*(α): = f + α (f^(eq) − f) satisfying

First, we follow Ref. [15] and define for non-negative f_i,

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Function H(f*(α)) − H(f) which is convex and increasing on [1, αmax]. is the solution to Eq. (4).

Recall from Ref. [15] that the viscosity relation (2) is derived with α = 2. Then α should be close to 2 as much as possible in order to maintain the consistency. To this end, we introduce

About this formula, three remarks are in order. (I) Formula (10) with Eq. (8) not only guarantees the nonnegativity of the distributions but also maintains the discrete H-theorem (6) The former is due to α ∈ [0, α_max] and the latter follows from H(f^*(α_*)) ≤ H(f) or H(f^*(α⁽¹⁾)) ≤ H(f). (II) The introduction of α_* (≤ 2) is a key, it reduces the computational cost drastically. Indeed, α_* is often much smaller than α_max used in Refs. [14], [15], [30], and [31]. It also extracts lots of grid points where the entropy balance (4) is irrelevant. Moreover, our numerical example shows that the number of grid points where H(f^*(2)) ≤ H(f) is about half of the total grid point number, see Fig. 5. (III) Formula (10) is much simpler than those given in Ref. [21]. It relies neither on any near-equilibrium assumption nor on the specific form of the entropy function H = H(f), while those in Ref. [21] do.

We conduct several numerical experiments in this section to demonstrate our relaxation-rate formula. We firstly compare formula (10) and the essentially ELBM (EELBM) with its first-order approximation^[21] by simulating the one-dimensional shock tube problem in Ref. [10]. The exact density profile is displayed in Fig. 2 and the two formulae both produce solutions oscillating near the shock. We see that in a narrow region of the shock front, α obtained with the formulae is not larger than 2 and our formula gives larger α than the EELBM. At the point of maximum departure, the deviation of α from 2 is 10.71% for the EELBM while that for the present formula is only 3.96%. This clearly shows that our formula gives a better approximation to the solution of Eq. (4) than the first-order EELBM.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (a) The exact density profile for the shock tube problem at time t = 500. The initial state is ρ(0 ≤ x ≤ 400) = 1.5, ρ(400 < x ≤ 800) = 0.75, and u(0 ≤ x ≤ 800) = 0. (b) Distributions of α with viscosity ν = 10−5 at the compressive shock front from formula (10) (red solid line) and the first-order EELBM (black dashed line).

Furthermore, we simulate the double shear flow with periodic boundary conditions. For this problem, the spatial domain is [0,1] × [0,1], and the initial state is^[20,33]

To show the effectiveness and efficiency of formula (10), we solve this problem with three approaches: formula (10), the first-order EELBM,^[21] and the bisection method used in Ref. [31] for Eq. (4). The stop criterion for the bisection method is −10⁻¹³ ≤ H(f^*(α)) − H(f) ≤ 0. Here we use the D2Q9 lattice, for which the equilibrium minimizing entropy (3) has an analytical expression,^[16] and the mesh size is taken as M² = 128 × 128. With this set, running 3200 steps of computations with the Intel(R) Core(TM) i5-6600 CPU takes 55.237 s, 21.564 s, and 607.411 s for the present implementation, the first-order EELBM, and the bisection method, respectively. Clearly, the present implementation is much faster than the bisection method and it is only about 2.5 times slower than the first-order EELBM with explicit expression of α. The latter is because several log functions need to be computed for the entropy function in formula (10). Note that though the first-order EELBM is fast, it is too dissipative to produce incorrect results. To see this, we display the results produced with the three approaches in Fig. 3. They are the contour lines of vorticity at t = 1 (t = TU/M, where T = 3200 is the number of time steps and M = 128 is the mesh size). It can be seen that, except for the EELBM, the other two approaches yield almost the same shape of vortex that is consistent with that in Refs. [20] and [33]. Furthermore, we also plot the distributions of α in Fig. 4 for the approaches above. It clearly shows that our α is smaller than 2 near the vortexes and close to 2 elsewhere. These demonstrate the efficiency of formula (10).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Vorticity contours of the double shear flow at t = 1 with ν = 10−5 obtained by (a) formula (10), (b) the EELBM, and (c) the bisection method.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Distributions of α at t = 1 with ν = 10−5: (a) formula (10), (b) the EELBM, and (c) the bisection method.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The proportion P[H(f*(2)) ≤ H(f)]/M2 as a function of time t, where P[H(f*(2)) ≤ H(f)] denotes the number of grid points where H(f*(2)) ≤ H(f) and M2 is the total grid point number.

To have a closer look at formula (10), we also count the number, by P[H(f^*(2)) ≤ H(f)], of the grid points where H(f^*(2)) ≤ H(f) at each time step. The proportion P[H(f^*(2)) ≤ H(f)]/M² as a function of time t is plotted in Fig. 5. We see that the proportion is around 0.5 for most of the time. Namely, at each time step, equation (4) needs to be solved only for about half of the total grid points. Therefore, the introduction of α = 2 in Eq. (8) enhances the efficiency significantly.

We have presented a simple relaxation-rate formula (10) for the ELBM. By doing this, we solve a critical outstanding problem in the discrete kinetic theory for hydrodynamics. The formula not only guarantees the nonnegativity of the distributions and the discrete H-theorem, but also gives full consideration to the consistency with hydrodynamics. We demonstrate that, with our new implementation based on the relaxation-rate formula, the entropy balance equation (4) only needs to be solved for about half of the grid points, the value of α given with the formula effectively characterizes the drastic changes of flow fields, and the computational overhead of the ELBM is greatly reduced. Consequently, the ELBM is substantially advanced and now the model can be efficiently used for a broad range of hydrodynamics applications including turbulent flows. The formula can also be adapted to the high-order ELBM recently proposed in Ref. [29].

Furthermore, we show that formula (9) can be improved with the following iteration (a modified secant algorithm)

Finally, we would like to provide some remarks on the force terms in the ELBM. It is a very important problem to design an ELBM for flows with external forces. If the designed LBM can guarantee non-negativeness of the distributions functions and admit an H theorem, our formula can easily be adapted to the new model. However, such an ELBM does not seem available in the literature, while a treatment of the force terms was proposed for the ELBM in Ref. [34]. Thus, how to effectively deal with the source term in the framework of ELBM remains a critical problem and is beyond the scope of this paper.