Derivation of lattice Boltzmann equation via analytical characteristic integral

TL;DR

This paper introduces the analytical characteristic integral (ACI) lattice Boltzmann theory, rigorously deriving LB equations from the BGK Boltzmann equation, clarifying the role of collision number, and analyzing factors affecting LB method accuracy.

Contribution

It provides a rigorous derivation of LB equations from kinetic theory using characteristics, and redefines the characteristic parameter as collision number rather than relaxation time.

Findings

ACI LB theory supports existing LB models including LBGK.

Collision number is the key parameter, not relaxation time.

Accuracy depends on the evolution of equilibrium distribution along characteristics.

Abstract

A lattice Boltzmann (LB) theory, analytical characteristic integral (ACI) LB theory, is proposed in this paper. ACI LB theory takes Bhatnagar-Gross-Krook (BGK) Boltzmann equation as the exact kinetic equation behind Navier-Stokes continuum and momentum equations and constructs LB equation by rigorously integrating BGK-Boltzmann equation along characteristics. It's a general theory, supporting most existed LB equations including the standard lattice BGK (LBGK) equation inherited from lattice-gas automata, whose theoretical foundation had been questioned. ACI LB theory also indicates that the characteristic parameter of LB equation is collision number, depicting the particle-interacting intensity in the time span of LB equation, instead of traditionally assumed relaxation time, and the over relaxation time problem is merely a manifestation of temporal evolution of equilibrium distribution along characteristics under high collision number, irrelevant to particle kinetics. In ACI LB theory, the temporal evolution of equilibrium distribution along characteristics is the determinant of LB method accuracy and we numerically prove it.

Figures

Tables

Table 1. Table 1 : Configurations of unsteady Couette flow demonstrating cases. For all cases, the channel width equals L = 0.1 ​ m 𝐿 0.1 m L=0.1\;\rm{m} , and the driven velocity of upper wall is U = 1.0 ​ m / s 𝑈 1.0 m s U=1.0\;\rm{m/s} .

Case	CN111CN stands for collision number	Grid		Kinematic	Collision
name	type	number222Grid number is the mesh number on cavity side length $L$	Timestep	viscosity	number
	LCN		2.50E$-$4		3.33E$-$1
case1	HCN	201	2.50E$-$6	1.00E$-$3	3.33E$+$1
	LCN		1.00E$-$4		3.33E$-$1
case2	HCN	101	1.00E$-$6	1.00E$-$2	3.33E$+$1
	LCN		2.00E$-$5		3.33E$-$1
case3	HCN	101	2.00E$-$7	5.00E$-$2	3.33E$+$1
	LCN		1.00E$-$5		3.33E$-$1
case4	HCN	101	1.00E$-$7	1.00E$-$1	3.33E$+$1
	LCN		1.00E$-$6		3.33E$-$1
case5	HCN	101	1.00E$-$8	1.00E$+$0	3.33E$+$1

Table 2. Table 2 : Configurations of lid-driven cavity flow. The cavity side length equals 1.0 ​ m 1.0 m 1.0\rm{m} , the driven-lid velocity is 1.0 ​ m / s 1.0 m s 1.0\rm{m/s}

CN333CN stands for collision number	Grid		Kinematic	Collision
type	number444Grid number is the mesh number on cavity side length $L$	Timestep	viscosity	number
LCN		1.00E$-$3		5.09E$-$1
HCN	257	1.00E$-$6	1.00E$-$2	5.09E$+$2

Equations

$$\frac{\partial f}{\partial t}+\vec{\xi }\cdot \frac{\partial f}{\partial \vec{r}}=-\frac{1}{\lambda }\left(f-g\right),$$

$$g\equiv \frac{\rho }{{\left(2\pi RT\right)}^{D/2}}\exp \left(-\frac{{\left(\vec{\xi }-\vec{u}\right)}^2}{2RT}\right),$$

$$\rho =\int fd\vec{\xi }=\int gd\vec{\xi },$$

$$\rho \vec{u}=\int \vec{\xi }fd\vec{\xi }=\int \vec{\xi }gd\vec{\xi },$$

$$\rho RT=\frac{1}{2}\int {\left(\vec{\xi }-\vec{u}\right)}^{2}fd\vec{\xi }=\frac{1}{2}\int {\left(\vec{\xi }-\vec{u}\right)}^{2}gd\vec{\xi }.$$

$$\frac{\partial }{\partial t}=\epsilon \frac{\partial }{\partial t_1}+{\epsilon }^2\frac{\partial }{\partial t_2}+...,$$

$$\frac{\partial }{\partial \vec{r}}=\epsilon \frac{\partial }{\partial {\vec{r}}_1},$$

$$f=g+\sum _{n=1}^{\infty }{\epsilon }^n{f}^{\left(n\right)},$$

$$\int \sum _{n=1}^{\infty }{\epsilon }^n{f}^{\left(n\right)}d\vec{\xi }=0,$$

$$\int \sum _{n=1}^{\infty }{\epsilon }^n{f}^{\left(n\right)}\vec{\xi }d\vec{\xi }=\vec{0},$$

$$\int \sum _{n=1}^{\infty }{\epsilon }^n{f}^{\left(n\right)}{\left(\vec{\xi }-\vec{u}\right)}^{2}d\vec{\xi }=0.$$

$$\int f^{\left(n\right)}\varphi \left(\vec{\xi }\right)d\vec{\xi }=0,$$

$$\varphi \left(\vec{\xi }\right)=A+\vec{B}\cdot \vec{\xi }+C{\xi }^2,$$

$$\frac{\partial g}{\partial t_1}+\vec{\xi }\cdot \frac{\partial g}{\partial {\vec{r}}_1}=-\frac{1}{\lambda }{f}^{\left(1\right)},$$

$$\frac{\partial g}{\partial t_2}+\frac{\partial f^{\left(1\right)}}{\partial t_1}+\vec{\xi }\cdot \frac{\partial f^{\left(1\right)}}{\partial {\vec{r}}_1}=-\frac{1}{\lambda }{f}^{\left(2\right)}.$$

$$\frac{\partial \rho }{\partial t_1}+\frac{\partial \rho \vec{u}}{\partial {\vec{r}}_1}=0,$$

$$\frac{\partial \rho \vec{u}}{\partial t_1}+\frac{\partial \rho \vec{u}\vec{u}}{\partial {\vec{r}}_1}+\frac{\partial \rho RT}{\partial {\vec{r}}_1}=0;$$

$$\frac{\partial \rho }{\partial t_2}=0,$$

$$\frac{\partial \rho \vec{u}}{\partial t_2}+\frac{\partial }{\partial \vec{r}}\cdot \left(-\rho RT\lambda \left(\left[\nabla \vec{u}\right]+{\left[\nabla \vec{u}\right]}^T\right)\right)=0.$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho \vec{u}}{\partial \vec{r}}=0,$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho \vec{u}\vec{u}}{\partial \vec{r}}=-\frac{\partial p}{\partial \vec{r}}+\frac{\partial }{\partial \vec{r}}\cdot \left(\nu \rho \left(\left[\nabla \vec{u}\right]+{\left[\nabla \vec{u}\right]}^T\right)\right),$$

$$\lambda =\frac{\nu }{RT}.$$

$$\frac{df}{dt}=-\frac{1}{\lambda }\left(f-g\right),$$

$$\frac{d}{dt}\equiv \frac{\partial }{\partial t}+\vec{\xi }\cdot \frac{\partial }{\partial \vec{r}}.$$

$$\frac{d}{dt}{e}^{t/\lambda }f=\frac{1}{\lambda }{e}^{t/\lambda }g.$$

$$e^{\Delta t/\lambda }{f}^{n+1}-f^n=\underset{0}{\overset{\Delta t}{\int }}\frac{1}{\lambda }{e}^{t^{\prime }/\lambda }g\left(t^{\prime }\right)d{t}^{\prime }.$$

$$f^{n+1}$$

$$f^n$$

$$g\left(t^{\prime }\right)$$

$$g\left(t^{\prime }\right)=g^n,t^{\prime }\in \left[0,\Delta t\right],$$

$$f^{n+1}-f^n=-\left(1-e^{-\Delta t/\lambda }\right)\left(f^n-g^n\right).$$

$$g\left(t^{\prime }\right)=g^n+\frac{t^{\prime }}{\Delta t}\left(g^{n+1}-g^n\right).$$

$$f^{n+1}-f^n=$$

$$-\left(\frac{-\frac{\Delta t}{\lambda }}{e^{-\frac{\Delta t}{\lambda }}-1}-1-\frac{\Delta t}{2\lambda }\right)\left(f^{n+1}-f^n\right).$$

$$f^{n+1}-f^n=-\frac{\Delta t}{\lambda }\left(f^n-g^n\right)$$

$$h^{n+1}-h^n=-\left(1-e^{-\Delta t/\lambda }\right)\left(h^n-g^n\right),$$

Full text

**Keywords: ** LB equation, analytical characteristic integral; characteristic parameter

**PACS: **47.11.-j, 02.70.-c

1 Introduction

In the past few decades, the lattice Boltzmann (LB) method has been proven as an efficient alternative approach to many computational fluid dynamics (CFD) areas, e.g. hydrodynamic systems\uciteLiuZou-178,BaochangZhaoli-179, multiphase and multicomponent fluids\uciteLaddVerberg-186,SheikholeslamiGorji-Bandpy-163,LiuCheng-181,SemmaElGanaoui-240,MountrakisLorenz-483, porous media flow\uciteMengGuo-2046,LiuHe-2047. With the rapid development of multi-core super computer and parallel computation, LB method is gaining more attention from CFD researchers due to its inherent parallelism\uciteZhangSun-2775,WuLiu-2776,SunJiang-2777. Historically, LB method arose from the pioneering idea of lattice-gas automata\uciteFrischHasslacher-2751, however more recently we tend to link it with the kinetic equation, Bhatnagar-Gross-Krook (BGK) Boltzmann equation\uciteBhatnagarGross-95, which leads to the reconstruction of LB theory\uciteDellar-2074,HeChen-173,Hwang-2023,YongZhao-2019,BoeschKarlin-2075,HeLuo-2016,Shan-2582.

Although the newly constructed LB theory has abandoned the original lattice-gas theory, the interest to reproduce the classical lattice BGK (LBGK) equation inherited from lattice-gas automata, which has been proven numerically stable, is still considerable. To achieve this, numerous approaches have been proposed, e.g. Strang splitting\uciteDellar-2074, Maxwell iteration\uciteYongZhao-2019, Taylor expansion\uciteHeLuo-2016, SterlingChen-252, etc. However, the approximation and truncation employed in these approaches make the reproduction unconvincing. Taking the most popular Taylor expansion scheme as an instance, LB equation is dealt as a characteristic integral of BGK Boltzmann equation with constant collision term. To fix the truncation error introduced by constant collision term assumption, the scheme expands the left-hand-side integral of BGK Boltzmann equation in Taylor series during recovering Navier-Stokes equations with Chapmann-Enskog expansion\uciteSterlingChen-252. Employing the 2nd-degree Taylor expansion, the LBGK equation can be properly recovered. But the truncated Taylor series introduce uncertainty. And a further research\uciteHeChen-173 asserts that employing Taylor expansion scheme to recover Navier-Stokes energy equation, the viscosity in viscous heat dissipation term is inconsistent with momentum equation. Hence He et al\uciteHeChen-173 employ a trapezoidal rule to improve the approximate accuracy of collision term to avoid Taylor expansion. Recently, Boesch and Karlin\uciteBoeschKarlin-2075 use Euler-Maclaurin integral to analytically integrate the collision term, eliminating the truncation error of trapezoidal rule. Compared with trapezoidal rule, the result of Euler-Maclaurin integral (exact LB equation called in the paper), is theoretically more rigorous. The problem is that the elegant but extremely complicated Euler-Maclaurin integral cloaks the important effect of equilibrium distribution. It makes the derivation in Ref. [19] be founded on a quite rough approximation of equilibrium distribution, which leads Boesch and Karlin to the conclusion that classical LBGK equation can not be recovered from analytical characteristic integral of BGK-Boltzman equation. Especially, when the viscosity is minimal, the relaxation time in LBGK equation would be greater than 1 tending to 2, i.e. overrelaxation, meanwhile it is always under 1 in exact LB equation. Thus Boesch and Karlin\uciteBoeschKarlin-2075 assert that the overrelaxation in LBGK equation is beyond the kinetic theory of continuous-time BGK-Boltzmann equation. Then an essential question arises: can LBGK equation be exactly reproduced based on BGK-Boltzmann equation?

In this paper, we propose an LB theory to link the LB method and BGK-Boltzmann equation, designated as analytical characteristics integral (ACI) LB theory. ACI theory follows the philosophy of Ref. [19, 16] to construct LB equation, i.e. improving the accuracy of collision term integral to eliminate the truncation error and inconsistence in Taylor expansion scheme. Through analyzing BGK-Boltzmann equation, ACI theory asserts that BGK-Boltzmann equation is accurate enough to describe the particle kinetics behind Navier-Stokes continuum and momentum equations, and LB equation can be constructed via analytically integrating BGK-Boltzmann equation along characteristics with approximated temporal evolution of equilibrium distribution. ACI theory employs differential equation solving skill to achieve the analytical integral. Compared with Euler-Maclaurin integral employed in Ref. [19], differential equation solving skill is mathematically more concise with a clear physical figure. Meanwhile the rigorous characteristic integral avoids the truncation error of trapezoidal rule in Ref. [16]. In ACI theory, the temporal evolution of equilibrium distribution along characteristics is the kernel of constructing LB equation, determining the computational stability and accuracy of LB equations. It extends the generality of BGK-Boltzmann kinetic theory discussed in Ref. [19]. To demonstrate it, we recover several popular LB equations and numerically analyze them. Our derivation also indicates that the characteristic parameter of LB equation is collision number instead of the traditionally assumed relaxation time, which is merely a reflection of the temporal evolution of equilibrium distribution along characteristics, and the over relaxation time problem is a manifestation of the evolution model behind LBGK equation under high collision number, supported by the kinetic theory of continuous-time BGK-Boltzmann equation.

This paper is organized as follows. In Sec. 2   , we propose ACI LB theory and recover popular LB equations. In Sec. 3   , we demonstrate the discretion of LB equation under ACI LB theory. In Sec. 4   , we numerically analyze the derived LB equations with several simple benchmarks. Sec. 5    concludes the paper.

2 Analytical characteristic integral LB theory

BGK-Boltzmann equation is a result of linear approximation on collision term of Boltzmann equation\uciteBhatnagarGross-95,

[TABLE]

where $f\equiv f\left({\vec{r},\vec{\xi},t}\right)$ is the particle distribution, $\vec{\xi}$ is particle microscopic velocity, $\lambda$ is the collision time, and $g$ is the Maxwell-Boltzmann distribution,

[TABLE]

where $R$ is the ideal gas constant, $D$ is the dimension of the space, and $\rho$, $\vec{u}$, and $T$ are the macroscopic density of mass, velocity and temperature, respectively. The macroscopic variables, $\rho$, $\vec{u}$, $T$ are the (microscopic velocity) moments of the distribution function $f$:

[TABLE]

BGK-Boltzmann equation is a kinetic equation describing the transient interaction between particles. In LB theory, we use it to depict the particle kinetics behind Navier-Stokes continuum and momentum equations. It can be mathematically validated with classical Chapman-Enskog (CE) expansion\uciteChapmanCowling-93. CE expansion is a multiple scale technique, proposed to solve the Boltzmann equation at first\uciteChapmanCowling-93. The key concept of CE expansion is formulating the flow phenomena with different time scales. It decomposes the partial derivatives with respect to $t$ and ${\vec{r}}$ into different scales,

[TABLE]

where $\varepsilon$ is a parameter who is small enough to identify the scales, $t_{i}$ and ${\vec{r}}_{1}$ are the relative time and space scales respectively. Specifically, $t_{1}$ and $t_{2}$ represent the time scales of convection and diffusion respectively\uciteHwang-2023. Similarly, we expand the particle distribution $f$ with parameter $\varepsilon$,

[TABLE]

where $g$ is the Maxwell-Boltzmann distribution. The expansion implies that the deviation between particle distribution $f$ and equilibrium distribution $g$ is limited. Combining with the macroscopic variable formulas in Eq. (3), the integrals of the expansion for macroscopic variables yield

[TABLE]

In order to solve the BGK-Boltzmann equation, CE expansion employs an extra assumption that the first three velocity moment integrals of the decomposed distribution equal zero on each order of $\varepsilon$. For sake of simplicity but without losing generality, we use a linear combination of first three moments $\phi\left({\vec{\xi}}\right)$ to describe the assumption,

[TABLE]

where

[TABLE]

where $A$ and $C$ are arbitrary constants, $\vec{B}$ is arbitrary constant vector.

Now, we substitute the expansions in Eq. (4) and Eq. (5) into BGK-Boltzmann equation in Eq. (1) and collect the terms with same order of $\varepsilon$ to generate new equations on difference scales. The first two scaling equations of $\varepsilon$ read:

[TABLE]

Integrating Eq. (8) on the first two velocity moments with the assumption in Eq. (7), we can get

[TABLE]

[TABLE]

Incorporating the formulas in Eq. (9) and Eq. (10) based on their orders of velocity moment respectively, the macroscopic hydrodynamic equations are recovered for BGK-Boltzmann equation with Maxwell-Boltzmann distribution,

[TABLE]

where $p=\rho RT$ is the pressure, $\nu=RT\lambda$ is the kinematic viscosity. Then the collision time $\lambda$ can be solved as

[TABLE]

It should be noted that this recovery of Navier-Stokes equations did not involve Taylor expansion, avoiding the truncation error. The derivation shows that BGK-Boltzmann equation is accurate enough to describe the particle kinetics behind Navier-Stokes continuum and momentum equations. Assuming BGK-Boltzmann equation is the exact kinetic equation for these two equations, all we need to do is constructing LB equation based on BGK-Boltzmann equation. Here we follow the philosophy of Boesh-Karlin approach\uciteBoeschKarlin-2075, i.e. constructing LB equation through analytically integrating BGK-Boltzmann equation along characteristics. Actually, with simple equation transformation\uciteHeLuo-2016, all terms in BGK-Boltzmann equation can be directly integrated except the equilibrium term (Maxwell-Boltzmann distribution). We start with transforming the left side of BGK-Boltzmann equation, treating it as the temporal derivative along characteristic line $\vec{\xi}$,

[TABLE]

where

[TABLE]

The Eq. (13) can be converted into integrable form by introducing an integrating factor $e^{t/\lambda}$,

[TABLE]

Hence LB equation can be generated by directly integrating BGK-Boltzmann equation along characteristics, it reads

[TABLE]

It should be noted that we have introduced short notations to simplify the expression in Eq. (15),

[TABLE]

This notating convention would be kept in the following and extended into other variables.

The integration in Eq. (15) hasn’t been finished, where the integral of $g\left({t^{\prime}}\right)$ still remains unsolved. Since the calculation of $g\left({t^{\prime}}\right)$ involves macroscopic density $\rho$ and velocity $\vec{u}$ (the temperature $T$ would be dealt as a field constant in practice), which evolute nonlinearly and need to be calculated by LB equation, it is impossible to get its analytical distribution along the characteristics. The only solution is using approximation. Since BGK-Boltzmann equation is assumed exact for the particle kinetics behind Navier-Stokes continuum and momentum equation, and LB equation is its analytical integration along characteristics, then the accuracy of LB equation solely depends on the design of $g\left({t^{\prime}}\right)$ evolution along characteristics, i.e. the aforementioned temporal evolution of equilibrium distribution. In other words $g\left({t^{\prime}}\right)$ model forges the final form of LB equation. It greatly extends the kinetic theory of BGK-Boltzmann equation discussed in Boesh-Karlin approach\uciteBoeschKarlin-2075, which we will detailedly demonstrate in the following. In case the evolution of BGK-Boltzmann equation departs the physical process too much, we restrict the value of $g\left({t^{\prime}}\right)$ at $t^{\prime}=0$ to $g^{n}$ , ensuring that the transient equation of particle interaction with respect to $\left({\vec{r},\vec{\xi},t}\right)$ is exact BGK-Boltzmann equation. We call this LB theory as analytical characteristic integral (ACI) LB theory.

We start with typical dealing for an unknown distribution, i.e. assuming $g\left({t^{\prime}}\right)$ constant over $[0,\Delta t]$, namely steady assumption (SA) model,

[TABLE]

then its LB equation reads

[TABLE]

Ref. [20] proposes linear interpolation between $g^{n}$ and $g^{n+1}$,

[TABLE]

Substituting it into Eq. (15), we can get its corresponding LB equation,

[TABLE]

It worths noting that though Ref. [20] got this analytical integration, the authors had expanded it in Taylor series to recover LBGK equation form,

[TABLE]

which makes the analytical integral of collision term degrade into constant approximation. In order to fix the truncation error introduced by constant collision term approximation, Taylor expansion is still inevitable during recovery of Navier-Stokes equations. Therefore, the scheme in Ref. [20] is still a classical Taylor expansion approach regardless of its analytical integrating skill (See Ref. [20] for more detail). The full equation in Eq. (2) equals exact LB equation proposed by Boesh-Karlin approach\uciteBoeschKarlin-2075, denoted as Boesch-Karlin LB equation in this paper. Compared with Euler-Maclaurin integral in Ref. [19], the derivation in ACI theory is mathematically more concise with a clear physical figure. In practice, the implicitness of $g^{n+1}$ would be removed by introducing an artificial distribution $h$, a trick proposed by Ref. [16]. Thus the implemented Boesch-Karlin LB equation is

[TABLE]

where

[TABLE]

In computation, we deal $h$ as $f$, neglecting their difference. Then the practical Boesch-Karlin LB equation equals our aforementioned SA LB equation. As the relaxation time in Eq. (21) is always $\leq 1$, Boesch and Karlin assert that the overrelaxation in LBGK LB equation is beyond the kinetic theory of continuous-time BGK-Boltzmann equation\uciteBoeschKarlin-2075, which is questionable since the $g\left({t^{\prime}}\right)$ model behind Boesch-Karlin LB equation is quite rough.

Now the question is what’s the $g\left({t^{\prime}}\right)$ model behind LBGK equation. As we assume that LB equation is an analytical characteristic integral of BGK-Boltzmann equation and $g\left({t^{\prime}}\right)$ model ensures that the transient equation of particle interaction with respect to $\left({\vec{r},\vec{\xi},t}\right)$ is exact BGK-Boltzmann equation, then the $g\left(t^{\prime}\right)$ model is not difficult to deduce, namely evolutionary correction of deviation (ECD) model,

[TABLE]

It is emphasized that this equation is mathematical development for LBGK equation under ACI theory. Substituting it into the characteristic integral of BGK-Boltzmann equation, we reproduce LBGK equation,

[TABLE]

It worths noting that after substituting ECD model, BGK-Boltzmann equation is integrable directly without requiring the introduce of integrating factor $e^{t/\lambda}$. Similarly, the constant collision model which assumes that the collision changing rate along the characteristics is constant, namely direct correction of deviation (DCD) model, reads

[TABLE]

and its corresponding LB equation yields

[TABLE]

We can also deduce the actual $g\left({t^{\prime}}\right)$ model behind the LB equation proposed by He et al\uciteHeChen-173, designated as He-Chen-Dong LB equation in this paper,

[TABLE]

its LB equation reads

[TABLE]

Removing the implicitness of $g^{n+1}$, the equation turns into

[TABLE]

where

[TABLE]

Ignoring the difference between $h$ and $f$, He-Chen-Dong LB equation is equivalent with the standard LBGK equation (named ECD LB equation in this paper).

Summarizing all LB equations derived in this section, SA, Boesh-Karlin, ECD, DCD and He-Chen-Dong LB equations, they all are analytical characteristic integrals of BGK-Boltzmann equation. The difference is their designs of $g\left(t^{\prime}\right)$ along the characteristics. In practice, Boesh-Karlin and He-Chen-Dong LB equation equal SA and ECD LB equation respectively. All these LB equations are implemented with a unified form,

[TABLE]

with

[TABLE]

It indicates that $1/\tau$, the relaxation time of LB method in classical LB theories\uciteBoeschKarlin-2075,HeChen-173,Hwang-2023, is actually a reflection of $g\left(t^{\prime}\right)$ model along the characteristics. It also worths noting that LB equation, the characteristic integral of BGK-Boltzmann equation, is by no means limited to the form in Eq. (31).

We close this section with a few remarks:

1. $\bullet$

In this section, we propose ACI LB theory. In ACI LB theory, we take BGK-Boltzmann equation as the exact kinetic equation behind Navier-Stokes continuum and momentum equations and construct LB equation through rigorously integrating BGK-Boltzmann equation along characteristics. BGK-Boltzmann equation can be analytically integrated except Maxwell-Boltzmann distribution, which evolutes nonlinearly and needs to be calculated by LB equation. Thus the design of $g\left(t^{\prime}\right)$, temporal evolution of equilibrium distribution along characteristics, decides the accuracy of LB equation. In case the whole evolution of BGK-Boltzmann equation departs the physical process too much, we restrict the value of $g\left(t^{\prime}\right)$ at $t^{\prime}=0$ to $g^{n}$, ensuring that the transient equation of particle interaction with respect to $\left({\vec{r},\vec{\xi},t}\right)$ is exact BGK-Boltzmann equation. 2. $\bullet$

In ACI LB theory, there are two important parameters, the collision time $\lambda$, depicting the time required for the distribution to reach equilibrium through collision, and the step time $\Delta t$, denoting the time interval of the LB equation. Their ratio, the collision number $\Delta t/\lambda$, describes the times of the distribution reaching equilibrium state. It reflects the intensity of particle interaction in the time span of LB equation, independent from detailed $g\left(t^{\prime}\right)$ models. By contrast, the traditionally assumed characteristic parameter of LB equation, relaxation time, is merely a reflection of $g\left(t^{\prime}\right)$ model. Thus the actual characteristic parameter of LB equation should be collision number $\Delta t/\lambda$, which is supported by Eq. (32). 3. $\bullet$

ACI LB theory is a physically general theory. By tuning $g\left(t^{\prime}\right)$ model, ACI LB theory can easily reproduce SA, Boesh-Karlin, ECD, DCD and He-Chen-Dong LB equations which are developed under different theories. It makes the comparison of LB equations avoid the ambiguous theoretical analysis of their derivations, leads to a simple investigation of their designs of $g\left(t^{\prime}\right)$ model. Another direct implication of this generality is that the form of LB equation is by no means limited to Eq. (31) since the flexibility of $g\left(t^{\prime}\right)$ model is far more great than the presented ones in this section. 4. $\bullet$

ACI LB theory is a mathematically rigorous but concise theory. It avoids the viscosity inconsistence\uciteHeChen-173 introduced by truncating Taylor series in Taylor expansion schemes\uciteHeLuo-2016, SterlingChen-252. As an analytical characteristic integral approach like Boesch-Karlin approach, ACI theory eliminates the uncertainty of trapezoidal rule in He-Chen-Dong approach\uciteHeChen-173 (See Ref. [19] for more detailed discussions). Compared with Boesch-Karlin approach\uciteBoeschKarlin-2075, the integrating skill of ACI theory is mathematically more concise with a clear physical figure. And the asserted concept that $g\left(t^{\prime}\right)$ model forges LB equation significantly extends the understanding of BGK-Boltzmann kinetic theory in Boesch-Karlin approach. ACI theory also avoids the uncertainty of collision and streaming splitting operation in Strang splitting approach\uciteDellar-2074. 5. $\bullet$

ACI LB theory is sensitive. The difference among LB equations will clearly manifest in their corresponding $g\left(t^{\prime}\right)$ models. For example, the $g\left(t^{\prime}\right)$ models corresponding to Boesh-Karlin and He-Chen-Dong LB equation are different from their implemented equations, i.e. SA and ECD LB equation.

3 Discretization of LB equation

In the former section, we have demonstrated the derivation of LB equation from BGK-Boltzmann equation. But the whole argument is based on velocity-continuous equation, which should be discretized before implementation. In ACI LB theory, the discretization of LB equation in velocity space is independent from the construction of LB equation, it is directly performed on the BGK-Boltzmann equation. The discretion includes discretizing the velocity space of particle distribution and constructing the equilibrium distribution based on the discrete-velocity model. Technically speaking, the discrete equilibrium distribution model determines the final recovery form of hydrodynamic equations, such as compressible or incompressible equations, conduct equation and so on, but it is out of this paper’s discussion. Here we tend to offer a detail demonstration and illustrate the framework of ACI theory, i.e. the discretization of velocity space on BGK-Boltzmann equation, determining the practical LB algorithm, is independent from detailed LB equation. We take a classical discretization as the detail demonstration, small Mach-number approximation (SMA) approach, proposed by Ref. [20]. For the sake of simplicity but without losing generality, we employ 2-dimensional Maxwell-Boltzmann distribution to construct the classical D2Q9 model. The derivation can be easily extended to 3 dimensions and constructing new models.

The key idea of velocity discretization in ACI LB theory is to ensure that the assumed macroscopic hydrodynamic equations still can be recovered by the discretized BGK-Boltzmann equation. To achieve this, SMA approach employs keeping the relative velocity moment integrals of discretized equilibrium distribution consistent with Maxwell-Boltzmann distribution. For an instance, the following velocity moment integrals in Eq. (33) have been used in Sec. 2    to recover the Navier-Stokes continuum and momentum equations,

[TABLE]

where $\xi_{\alpha}$ is the component of $\vec{\xi}$ in Cartesian coordinates. So the discrete model should keep equivalent with Maxwell-Boltzmann distribution on moments $1,\vec{\xi},\dots,\vec{\xi}^{3}$.

To construct the discrete equilibrium distribution, SMA approach firstly decomposes the exponent part in the Maxwell-Boltzmann distribution: one for the weight function, the other for the distribution with respect to macro velocity,

[TABLE]

For the convenience of calculating weights, the truncated small velocity expansion (or small-Mach-number approximation) is employed\uciteHeLuo-2016, where the terms above 2nd macro-velocity order are neglected,

[TABLE]

It worths noting that the order of small-Mach-number approximation decides the velocity moment integrals’ accuracy of $g^{\prime}$, i.e. the highest order of velocity moment integral that $g^{\prime}$ can restore. Taking Eq. (3) for example, the 2nd order truncation could only restore the velocity moment integral up to 2nd order. The 3rd order moment integral of the truncated Maxwell-Boltzmann distribution in Eq. (37) would deviate from the original in Eq. (3),

[TABLE]

[TABLE]

To keep equivalent on 3rd velocity moment integral, $g^{\prime}$ needs to maintain macro-velocity terms up to 3rd order. Here are 3rd order approximation (seeing Eq. (38)) and its 3rd velocity moment integral (seeing Eq. (3)),

[TABLE]

[TABLE]

Unfortunately, the implementation of 3rd order approximation requires a more complicated discrete velocity set. In this paper, we still employ the 2nd-order approximation.

After small-Mach-number approximation, the target continuous equilibrium distribution has been $g^{\prime}$ instead of original Maxwell-Boltzmann distribution $g$. Right now, $g^{\prime}$ in Eq. (3) bears a strong resemblance to the classical D2Q9 equilibrium distribution, all remains to be accomplished is calculating the weight of discrete velocity. As we discussed before, the discrete equilibrium model $f^{eq}$ should keep the velocity moment integrals required in Eq. (33) consistent with $g^{\prime}$,

[TABLE]

where $\psi\left({\vec{\xi}}\right)$ is a polynomial of $\vec{\xi}$ up to order 3. For the sake of simplicity but without losing generality, assuming

[TABLE]

then the velocity moment integrals of $g^{\prime}$ can read as

[TABLE]

where

[TABLE]

which can be calculated numerically with Gaussian-type quadrature. Our object is using proper discretization of velocity space to evaluate the integral $I_{m}$ with $m$ up to 5. Naturally, the third-order Hermite formula\uciteHeLuo-2016 is the optimal choice for the purpose of deriving the D2Q9 model,

[TABLE]

The integrals need to be evaluated are, $I_{0}$, $I_{2}$, $I_{4}$. Other integrals with odd orders of $x$ equal 0 due to the symmetry of Eq. (44) which can be easily solved by using symmetrical abscissas of the quadrature. Setting the abscissas of quadrature as $x_{0}=0$,$x_{1}=-\zeta$, $x_{2}=\zeta$, we have the following three equations:

[TABLE]

with symmetrical weights

[TABLE]

The three abscissas of the quadrature are

[TABLE]

and the corresponding weight coefficients are

[TABLE]

Then the evaluation Eq. (40) can read as

[TABLE]

where ${{\vec{\xi}}_{\alpha}}={{\vec{\xi}}_{i,j}}=\sqrt{2RT}\left({{x_{i}},{x_{j}}}\right)$. It’s straightforward to identify the discrete equilibrium distribution with

[TABLE]

Combining with the notations of discrete velocity and corresponding weights, employing the relation $RT=c_{s}^{2}=c^{2}/3$, a complete D2Q9 model is constructed,

[TABLE]

with discrete velocity,

[TABLE]

where

[TABLE]

and weights,

[TABLE]

where $\rho=\sum{{f_{\alpha}}}$, $\vec{u}=\sum{{f_{\alpha}}{{\vec{\xi}}_{\alpha}}}/\rho$.

The derived D2Q9 model is directly constructed on BGK-Boltzmann equation, independent from detail LB equation. Thus it is applicable to all $g\left(t^{\prime}\right)$ models.

4 Numerical analysis

Within ACI LB theory, all popular LB equations are analytical characteristic integrals of BGK-Boltzmann equation, identified by their designs of $g\left(t^{\prime}\right)$ model, then all we need to analyze is their designs of $g\left(t^{\prime}\right)$ model. ACI LB theory also proposes a new characteristic parameter of LB equation, collision number $\Delta t/\lambda$, depicting the particle-interacting intensity in the time span of LB equation. Unlike the traditionally assumed characteristic parameter, relaxation time $1/\tau$, which is merely a reflection of $g\left(t^{\prime}\right)$ model, the collision number is independent from detail $g\left(t^{\prime}\right)$ models. In this section, we would analyze the numerical performance of derived LB equations under different collision numbers. In numerical computation, Boesh-Karlin and He-Chen-Dong LB equation would be replaced by SA and ECD LB equation, ignoring the implicitness as their proposers suggested\uciteBoeschKarlin-2075,HeChen-173. We start with the comparison of relaxation time curves along collision number as all derived LB equations share a unified equation form in Eq. (31), only distinguished in relaxation time formulas. Then two benchmarks, unsteady Couette flow and lid-driven cavity flow, will be calculated by SA and ECD LB equations, employing the derived D2Q9 model. The results of SA and ECD will be compared to analyze their accuracy under different collision numbers and validate the effect of $g\left(t^{\prime}\right)$ model. In the calculation of benchmarks, the DCD approximation will be abandoned due to its instability.

4.1 Relaxation time

As Eq. (32) shows, for the same calculation, i.e. the same collision number $\Delta t/\lambda$, the proposed $g\left(t^{\prime}\right)$ models, SA, DCD and ECD, only differ in relaxation time $1/\tau$. Fig. 1 illustrates the relaxation time profiles of SA, DCD and ECD with respect to collision number. For $\Delta t/\lambda\to 0$, SA, DCD and ECD approaches each other. As the collision number increases, all three relaxation times increase with different rate, which leads to $1/{\tau_{DCD}}>1/{\tau_{ECD}}>1/{\tau_{SA}}$. When the collision number goes into infinite, SA and ECD relaxation time would converge to 1 and 2 respectively, while DCD increase unlimitedly, which causes the instability of DCD. Fig. 1 indicates that SA, DCD and ECD LB equations would have same numerical performance under low collision number meanwhile distinguish from each other under high collision number. This conclusion is supported by the design of $g\left(t^{\prime}\right)$. Once collision number is small enough, which indicates small collision change, the particle distribution $f\left(t^{\prime}\right)$ along characteristics would not depart $f^{n}$ too much. Then the modification parts in ECD model Eq. (23) and DCD model Eq. (25) would both approach 0,

[TABLE]

Hence the $g\left(t^{\prime}\right)$ models of ECD and DCD can be approximated as

[TABLE]

which is exactly SA model. When the collision change increases as the collision number increases, the modification parts in Eq. (54) are not neglectable. Thus SA, ECD and DCD model would show different results. It worths noting that under large collision change, SA and DCD model are theoretically unsuitable as SA takes the equilibrium distribution constant and DCD takes the collision change rate constant, which contradict the physical process. What’s more, Fig. 1 also indicates the over relaxation time problem is merely a manifestation of $g\left(t^{\prime}\right)$ model under high collision number, irrelevant to particle kinetics.

4.2 Unsteady Couette flow

Unsteady Couette flow, also as known as transient plane Couette flow, is a typical numerical benchmark. It evolves in a straight channel, which infinitely extends in the x direction. The walls of the channel are parallel to the x axis and defined by the equations $y=0$ for the lower wall and $y=L$ for the upper wall. The flow is stationary at the beginning and driven by the upper wall with a constant velocity after that. Its equation can be described as

[TABLE]

with the following initial and boundary conditions in the range of $0\leq y\leq L$,

[TABLE]

where $\nu$ is kinematic viscosity, $U$ is the upper wall’s driven velocity, $L$ is the width of channel. Unsteady Couette flow can be solved analytically:

[TABLE]

The periodic boundary is implemented on x direction to simulate the infinite extension. For the upper and lower wall, the regularized boundary Ref. [24] is applied. Actually, in unsteady Couette flow, Zou-He boundary\uciteZouHe-251 and Inamuro boundary\uciteInamuroYoshino-238 would share the same macro result with regulariezed boundary, though the particle distributions $f$ may differ. 5 cases are designed to discuss the SA and ECD models (See Table 1). Each case has two states: low collision number(LCN) and high collision number(HCN). Fig. 2 shows the SA and ECD velocity profiles of unsteady Couette flow at $t=1.00E-3\;\rm{s}$. Fig. 3 illustrates the error evolution on time. The error formula is defined as

[TABLE]

where $u$ is the computing velocity profile, $u_{ana}$ is the analytical velocity profile. For the sake of space saving, only the error evolutions of case1,case3 and case5 are plotted.

The velocity profiles of SA LB equation in Fig. 2 show that at LCN, SA results are close to the analytical solutions meanwhile at HCN, they evidently deviate, manifesting magnified viscosities. By contrast, ECD results always keep consistent with analytical solutions regardless of collision number. The deviation of SA LB equation under HCN is explainable, due to the intensive collision change, the equilibrium distribution would be greatly changed, which leads to the failure of the constant equilibrium distribution assumption in SA model. And the results of ECD LB equation indicate that the modification part in ECD model significantly improves its numerical performance under HCN. Comparing the results of SA and ECD LB equation, it clearly shows that $g\left(t^{\prime}\right)$ model determines the accuracy of LB equation. The consistent numerical performance of SA and ECD LB equation under different collision numbers also validates that the collision number is a good characteristic parameter to depict LB equation. It reflects the collision intensity in the time span of LB equation, independent from detail $g\left(t^{\prime}\right)$ model.

4.3 Lid-driven cavity flow

Lid-driven cavity flow is a classical constant-property benchmark for fluid computation, for its simple geometry and complicated flow behaviors, especially the corner flow phenomena. Lid-driven cavity flow is a 2-D case characterized by Renold (Re) number, which develops in a square cavity with the side length $L$. The fluid would keep stationary at the beginning and be driven by the top lid with constant velocity after that. Except the top lid, all other walls keep rest. This case has no analytic solution. As a convention, the results in Ref. [27] are taken as a calibration on velocity profiles. As the definition of collision number shows, with small kinematic viscosity, the collision number is inherently large. It’s very difficult to simulate a high Re-number cavity flow under low collision number without mesh explosion. Fortunately, our mission is to evaluate $g\left(t^{\prime}\right)$ models under collision number instead of Re number. A small Re-number case with different collision numbers will fulfill the task. Hence the case with smallest Re number in Ref. [27] , $\rm{Re}=100$, is calculated for analysis.

The simulation is configured as Table 2. The cavity sides and corners employ the regularized boundary\uciteLattChopard-124. Fig. 4 illustrates the velocity profiles along the vertical and horizontal lines through cavity center. The models keep consistent with their performance in Sec. 4.2   , under LCN , SA and ECD share a close profile; under HCN, ECD is stick with the LCN profile while SA fails. To give a full picture of the calculations, the vorticity contours are plotted in Fig. 5, which can be quantified by post-processing software compared to streamline. The values of vorticity along contours are $\rm{0}$, $\pm\rm{0.5}$, $\pm\rm{1.0}$, $\pm\rm{2.0}$, $\pm\rm{3.0}$, $\rm{4.0}$, $\rm{5.0}$. As Fig. 5 illustrates, under LCN, SA (seeing Fig. 5(b)) and ECD (seeing Fig. 5(a)) have the same vorticity contour; under HCN, though the contour of ECD in Fig. 5(c) is a bit dispersed, the shape agrees well with LCN result while SA (seeing Fig. 5(d)) is quite different.

The results of SA model illuminate that, under LCN, as the collision change is minimal during a steptime, the constant equilibrium distribution assumption in SA LB equation can properly approximate its temporal evolution, meanwhile under HCN, due to the great collision change, it seriously deviates from the physical temporal evolution. The accurate result of ECD model under HCN indicates that the modification part in ECD model could accurately approximate the change of equilibrium distribution caused by intensive collision. The different numerical performances of SA and ECD LB equation under HCN confirm the effect of $g\left(t^{\prime}\right)$ model on accuracy of LB equation. Compared with unsteady Couette flow in Sec. 4.2   , the governing equations of lid-driven cavity flow, i.e. full incompressible-flow Navier-Stokes equations, is more complicated and close to a real CFD problem. Thus the consistence of models’ numerical performances between unsteady Couette and lid-driven cavity flow validates the universality of our assertion on the accuracy of SA and ECD LB equation under different collision numbers. Hence, we can conclude that, for a fluid calculation, SA and ECD are both available under LCN meanwhile ECD is a better choice under HCN.

5 Conclusion

In this paper we propose a mathematically rigorous but physically general LB theory, ACI LB theory. In ACI LB theory, we assume that BGK-Boltzmann equation is exact kinetic equation behind Navier-Stokes continuum and momentum equations and construct LB equation by analytically integrating BGK-Boltzmann equation along characteristics. Since the integral of $g\left(t^{\prime}\right)$ along characteristics can not be analytically solved, the approximation is employed. To ensure that the evolution of LB equation does not depart the physical process too much, we restrict the value of $g\left(t^{\prime}\right)$ at $t^{\prime}=0$ to $g^{n}$ so that the transient equation of particle interaction with respect to $\left({\vec{r},\vec{\xi},t}\right)$ is exact BGK-Boltzmann equation. In ACI LB theory, the $g\left(t^{\prime}\right)$ model is the determinant of LB equation accuracy and collision number is the characteristic parameter of equation, depicting the particle-interacting intensity in the time span. As a demonstration, we recover the approximating models behind popular LB equations, i.g. SA, Boesh-Karlin, ECD, DCD and He-Chen-Dong LB equation, and numerically analyze them, where Boesh-Karlin and He-Chen-Dong equation are replaced by SA and ECD since they are equivalent respectively after removing the implicitness of $g^{n+1}$.

To conclude, we highlight the following conclusions drawn from our work:

1. $\bullet$

LB equation can be constructed through analytically integrating BGK-Boltzmann equation along characteristics with approximated temporal evolution of equilibrium distribution $g\left(t^{\prime}\right)$, avoiding the approximation and truncation in classical schemes\uciteDellar-2074, HeLuo-2016, SterlingChen-252, HeChen-173. 2. $\bullet$

The kinetic theory depicted by BGK-Boltzmann equation is quite general, far beyond the discussion in Boesh-Karlin scheme\uciteBoeschKarlin-2075. By tuning $g\left(t^{\prime}\right)$ model, most popular LB equations can be recovered as analytical characteristic integrals of BGK-Boltzmann equation, including the questioned LBGK equation. 3. $\bullet$

LB equation is identified by its $g\left(t^{\prime}\right)$ model. The approximation accuracy of $g\left(t^{\prime}\right)$ model determines LB equation accuracy. Among presented $g\left(t^{\prime}\right)$ models, ECD is a better choice for numerical simulations. 4. $\bullet$

The characteristic parameter of LB equation is collision number, $\Delta t/\lambda$, instead of the traditional relaxation time, $1/\tau$, which is merely a reflection of $g\left(t^{\prime}\right)$ model. It depicts the particle-interacting intensity in the time span of LB equation, independent from detail $g\left(t^{\prime}\right)$ model. And the over relaxation time problem is merely a manifestation of ECD model under high collision number.

Acknowledgment

Huanfeng Ye would like to express his gratitude to Dr. Zecheng Gan for helpful discussion, and to Dr. Mathias J. Krause for generous support on OpenLB software. Huanfeng Ye is also grateful to the OpenLB team for their open source sharing (www.openlb.net).