Relaxation-rate formula for the entropic lattice Boltzmann method

In the last three decades, the lattice Boltzmann method (LBM) has become a popular mesoscopic numerical method for fluid dynamics, with applications ranging from high Reynolds number flows to flows at porous media and relativistic hydrodynamics Su ; Gu . We refer to BSV ; CD ; AC 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 Su ; Gu . By fitting the discrete velocities into a regular lattice, the LBM realizes the propagation and collision of distribution functions efficiently.

Though simple and efficient, the standard LBM is limited to moderate Reynolds number flows due to the lack of numerical stability Su ; AK ; BYCW . To alleviate this obstacle, the entropic LBM (ELBM) restoring the discrete H-theorem has been proposed in KGSB ; KFO ; AK2002 ; AK2002jsp ; AKO2003 ; BLY2004 ; CAK2006 ; KVYSV2007 ; KBC2014 ; AKTA2017 as a paradigm shift for computational fluid dynamics. The introduction of the discrete H-theorem in the ELBM significantly extends the operation range of the discrete kinetic theory to turbulent flows CFKTB2010 ; DBCBK2016 , binary droplet collisions MCK2015 ; MCK2016 , multiphase fluid-solid interface problems MCK2015pre and compressible flows FCK2015 . The solid physical background and successful applications make the ELBM a powerful approach for the study of complex flows.

A key step in the ELBM is to determine the relaxation-rate involving a parameter $\alpha$ introduced in KFO for ensuring the discrete H-theorem. It needs to solve a complicated nonlinear algebraic equation, which greatly affects the efficiency of the ELBM. Though considerable efforts have been made in the past AK2002 ; AK2002jsp ; CAK2006 ; TUSK2007 to improve the efficiency, it was only recognized by Brownlee et al. BGL2007 that the computational value of $\alpha$ should not lead to a numerical entropy increase. This requirement essentially guarantees the discrete H-theorem at the numerical level and is further emphasised in AKTA2017 recently. Therein some analytical approximate expressions of $\alpha$, which also guarantee the discrete H-theorem numerically, are derived by relaxing the entropy equality AKTA2017 and by making a near-equilibrium assumption. Nevertheless, the first-order approximation is too dissipative while higher-order approximations are difficult to be explicitly obtained. Therefore, a critical breakthrough is much needed for the efficient determination of $\alpha$.

In this Letter, we solve this long-standing problem by proposing an elegant and uniform formula for the parameter $\alpha$, or equivalently the relaxation-rate of the ELBM. This formula is based on a novel combination of the consistency of the ELBM and the constraint that the entropy must not increase within a discrete time step. Besides compliance of the discrete H-theorem at numerical level, an excellent property of the formula is that it is applicable to arbitrary convex entropy functions. Additionally, numerical simulations demonstrate that the grid points where flow fields change drastically are effectively marked by the formula.

Before presenting our formula, we recall from AKO2003 that the entropic lattice Boltzmann method (ELBM) reads as

with $i=1,2,\cdots,N$. Here $f_{i}=f_{i}(\bm{x},t)$ is the $i$-th distribution function for particles with velocity $\bm{c}_{i}$ at position $\bm{x}$ and time $t$, $\delta_{t}$ is the time step, $\alpha\beta$ is the relaxation-rate with $\beta\in(0,1)$ related to the fluid viscosity $\nu$ via

$c_{s}$ is the sound speed and $f_{i}^{(eq)}=f_{i}^{(eq)}(\bm{x},\,t)>0$ is the equilibrium minimizing the convex entropy function

subject to the conservation laws of mass and momentum (for the isothermal case):

In the above equations, $f$ stands for the vector $(f_{1},f_{2},\cdots,f_{N})$, $W_{i}>0$ is the $i$-th weight and $\rho$ and $\bm{u}$ are the macroscopic fluid density and velocity, respectively. A key point of the ELBM is the parameter $\alpha$ in (1) that maintains the entropy balance

If $\alpha=2$ and $f_{i}^{(eq)}$ is taken as polynomials, (1) degenerates to the standard lattice Bhatnagar-Gross-Krook (LBGK) model QDL .

With the equal entropy auxiliary distribution $f^{*}:=f+\alpha(f^{(eq)}-f)$, the ELBM (1) can be rewritten as

By the convexity of the function $H$ defined in (3), we see from (4) and (5a) that

Furthermore, for a periodic domain $\Omega$ we have

This is a discrete H-theorem for the ELBM.

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

Here we derive a simple approximate solution formula for $\alpha$. Observe that the discrete H-theorem (6) holds true for all $f^{*}(\alpha):=f+\alpha(f^{(eq)}-f)$ satisfying

instead of the entropy balance (4). Then we can relax Eq. (4) and replace it with the inequality (7) as in AKTA2017 . Based on this observation, we propose an efficient implementation different from that in AKTA2017 . In what follows, we assume that $f^{(eq)}\neq f$. Otherwise, $\alpha$ can be any number.

First, we follow AK2002jsp and define for non-negative $f_{i}$:

which partly measures the departure of distribution $f$ from equilibrium. Because $\sum_{i}f_{i}=\sum_{i}f_{i}^{(eq)}$, the above set is non-empty. Notice that $\alpha\in[0,\alpha_{\max}]$ ensures the nonnegativity of the distribution $f^{*}(\alpha)$ and thereby the entropy function $H(f^{*}(\alpha))$ is well defined AK2002jsp . Moreover, from the convexity of $H(f)$ it follows that $H(f^{*}(\alpha))$ is convex and increasing on $[1,\alpha_{\max}]$ (see Fig. 1 and the Supplementary Material).

Recall from AK2002jsp that the viscosity relation (2) is derived with $\alpha=2$. Then $\alpha$ should be close to 2 as much as possible in order to maintain the consistency. To this end, we introduce

Notice that $\alpha_{\max}\gg 1$ near equilibrium. If $H(f^{*}(\alpha_{*}))\leq H(f)$, then $\alpha$ is taken as $\alpha_{*}$, which maintains the entropy inequality (7). Otherwise, we refer to Fig.  1 and take

From Fig. 1 we can see that this $\alpha^{(1)}$ will be very close to the solution of the entropy balance equation (4) if $[H(f^{*}(\alpha_{*}))-H(f)]$ is small, which is of frequent occurrence. Thanks to the convexity of function $H$, it is proved in the Supplementary Material that

About the above implementation, three remarks are in order. (a). The implementation not only guarantees the nonnegativity of distributions but also maintains the discrete H-theorem (6). The former is due to $\alpha\in[0,\alpha_{\max}]$ and the latter follows from $H(f^{*}(\alpha_{*}))\leq H(f)$ or $H(f^{*}(\alpha^{(1)}))\leq H(f)$. (b). The introduction of $\alpha_{*}(\leq 2)$ is a key in our implementation, it reduces the computational cost drastically. Indeed, $\alpha_{*}$ is often much smaller than $\alpha_{\max}$ used in AK2002jsp ; AK2002 ; TUSK2007 ; BGL2007 . 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))\leq H(f)$ is about half of the total grid point number, see Fig. 5. (c). Formula (9) is much simpler than those given in AKTA2017 . It relies neither on any near-equilibrium assumption nor on the specific form of the entropy function $H=H(f)$, while those in AKTA2017 do.

To compare Formula (9) and the essentially ELBM (EELBM) with its first-order approximation AKTA2017 , we simulate the one-dimensional shock tube problem in AK . The exact density profile is displayed in Fig. 2 and the two implementations both produce solutions oscillating near the shock. We see that in a narrow region of the shock front, $\alpha$ obtained with the implementations is not larger than 2 and our implementation gives larger $\alpha$ than the EELBM. At the point of maximum departure, the deviation of $\alpha$ from 2 is 10.71% for the EELBM while that for the present implementation is only 3.96%. This clearly shows that Formula (9) gives a better approximation to the solution of the equation (4) than the first-order EELBM.

Furthermore, we simulate the double shear flow with periodic boundary conditions. For this problem, the spatial domain is $[0,1]\times[0,1]$, and the initial state is MB1997 ; KBC2014

and $\rho(x,y,0)=1$. Here $u_{x}$ and $u_{y}$ are the $x-$ and $y-$component of the fluid velocity, respectively, $U=0.04$ and $\lambda=80$ determines the slop of the shear layer.

To show the effectiveness and efficiency of our implementation, we solve this problem with three approaches: our implementation, the first-order EELBM AKTA2017 , and the bisection method used in BGL2007 for the equation (4). The stop criterion for the bisection method is $-10^{-13}\leq H(f^{*}(\alpha))-H(f)\leq 0$. Here we use the D2Q9 lattice, for which the equilibrium minimizing the entropy (3) has an analytical expression AKO2003 , and the mesh size is taken as $M^{2}=128\times 128$. For this setup, the results produced with the three approaches are displayed 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 methods yield almost the same shape of vortex that is consistent with those in MB1997 ; KBC2014 . Furthermore, we also plot the distributions of $\alpha$ in Fig. 4 for the three implementations above. It clearly shows that our $\alpha$ is smaller than 2 near the vortexes and close to 2 elsewhere.

These demonstrate the efficiency of our implementation.

To have a closer look at our implementation, we also count the number, by $P[H(f^{*}(2))\leq H(f)]$, of the grid points where $H(f^{*}(2))\leq H(f)$ at each time step. The proportion $P[H(f^{*}(2))\leq H(f)]/M^{2}$ 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 Eq. (4) needs to be solved only for about half of the total grid point number. Therefore, the introduction of $\alpha=2$ in (8) enhances the efficiency significantly.

In conclusion, we have presented a simple and uniform relaxation-rate formula for the ELBM. The formula solves a long-standing critical problem in kinetic theory approach: it 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 algebraic equation only needs to be solved for about half of the grid points, where the flow fields change drastically is effectively marked, and the computational overhead of the ELBM is greatly reduced. Thus the novel relaxation-rate formula makes a significant step for the development of the ELBM and now the method can be efficiently used for a broad range of hydrodynamics applications including turbulence flows.

Finally, we point out that Formula (9) can be improved with the following iteration (a modified secant algorithm)

for $k=1,2,3,\cdots$. In the Supplementary Material it is proved that this iteration converges unconditionally to the solution of Eq. (4) and $H(f^{*}(\alpha^{(k)}))\leq H(f)$ for all $k$.