Generalized lattice Boltzmann method: Modeling, analysis, and elements

TL;DR

This paper introduces a unified framework for the generalized lattice Boltzmann method (GLBM), compares four analysis techniques, and discusses implementation elements, showing their equivalence and how existing models fit into this framework.

Contribution

It provides a comprehensive, unified approach to modeling and analyzing GLBM, clarifies the relationship between different analysis methods, and offers practical elements for implementation.

Findings

Four analysis methods are mathematically equivalent.

Existing LB models can be derived from the proposed GLBM.

The framework simplifies the understanding and development of LB models.

Abstract

In this paper, we first present a unified framework for the modelling of generalized lattice Boltzmann method (GLBM). We then conduct a comparison of the four popular analysis methods (Chapman-Enskog analysis, Maxwell iteration, direct Taylor expansion and recurrence equations approaches) that have been used to obtain the macroscopic Navier-Stokes equations and nonlinear convection-diffusion equations from the GLBM, and show that from mathematical point of view, these four analysis methods are equivalent to each other. Finally, we give some elements that are needed in the implementation of the GLBM, and also find that some available LB models can be obtained from this GLBM.

Equations

$$f_{j}(\mathbf{x}+\mathbf{c}_{j}\Delta t,t+\Delta t)=f_{j}(\mathbf{x},t)-\mathbf{\Lambda}_{jk}f_{k}^{ne}(\mathbf{x},t)+\Delta t\big{[}G_{j}(\mathbf{x},t)+F_{j}(\mathbf{x},t)+\frac{\Delta t}{2}\bar{D}_{j}F_{j}(\mathbf{x},t)\big{]},$$

$$\displaystyle\textbf{Collison:}\ \ \tilde{f}_{j}(\mathbf{x},t)=f_{j}(\mathbf{x},t)-\mathbf{\Lambda}_{jk}f_{k}^{ne}(\mathbf{x},t)+\Delta t\big{[}G_{j}(\mathbf{x},t)+F_{j}(\mathbf{x},t)+\frac{\Delta t}{2}\bar{D}_{j}F_{j}(\mathbf{x},t)\big{]},$$

$$\displaystyle\textbf{Propagation:}\ \ f_{j}(\mathbf{x}+\mathbf{c}_{j}\Delta t,t+\Delta t)=\tilde{f}_{j}(\mathbf{x},t),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$f_j=f_j^{(0)}+ϵf_j^{(1)}+{ϵ}^2{f}_j^{(2)},$$

$$G_j=ϵG_j^{(1)}+{ϵ}^2{G}_j^{(2)},F_j=ϵF_j^{(1)}+{ϵ}^2{F}_j^{(2)},$$

$${\partial }_t=ϵ{\partial }_{t_1}+{ϵ}^2{\partial }_{t_2},\nabla =ϵ{\nabla }_1,$$

$$D_j{f}_j+\frac{\Delta t}{2}{D}_j^2{f}_j+\cdots =-\frac{1}{\Delta t}{\mathbf{\Lambda }}_{jk}{f}_k^{ne}+G_j+F_j+\frac{\Delta t}{2}{\bar{D}}_j{F}_j,$$

$$\epsilon^{0}:\ \ \mathbf{\Lambda}_{jk}\big{(}f_{k}^{(0)}-f_{k}^{eq}\big{)}=0,$$

$${ϵ}^1:D_{1j}{f}_j^{(0)}=-\frac{1}{\Delta t}{\mathbf{\Lambda }}_{jk}{f}_k^{(1)}+G_j^{(1)}+F_j^{(1)},$$

$${ϵ}^2:{\partial }_{t_2}{f}_j^{(0)}+D_{1j}{f}_j^{(1)}+\frac{\Delta t}{2}{D}_{1j}^2{f}_j^{(0)}=-\frac{1}{\Delta t}{\mathbf{\Lambda }}_{jk}{f}_k^{(2)}+G_j^{(2)}+F_j^{(2)}+\frac{\Delta t}{2}{\bar{D}}_{1j}{F}_j^{(1)}.$$

$$f_j^{(0)}=f_j^{eq},j=0,1,\cdots ,q-1.$$

$$\partial_{t_{2}}f_{j}^{(0)}+D_{1j}\big{(}\delta_{jk}-\frac{\mathbf{\Lambda}_{jk}}{2}\big{)}f_{k}^{(1)}+\frac{\Delta t}{2}D_{1j}\big{(}G_{j}^{(1)}+F_{j}^{(1)}\big{)}=-\frac{1}{\Delta t}\mathbf{\Lambda}_{jk}f_{k}^{(2)}+G_{j}^{(2)}+F_{j}^{(2)}+\frac{\Delta t}{2}\bar{D}_{1j}F_{j}^{(1)}.$$

$$D_j{f}_j^{eq}=-\frac{1}{\Delta t}{\mathbf{\Lambda }}_{jk}{f}_k^{ne}+G_j+F_j+O({ϵ}^2+{ϵ}^2/\Delta t),$$

$$D_{j}f_{j}^{eq}+D_{j}\big{(}\delta_{jk}-\frac{\mathbf{\Lambda}_{jk}}{2}\big{)}f_{k}^{ne}+\frac{\Delta t}{2}D_{j}\big{(}G_{j}+F_{j}\big{)}=-\frac{1}{\Delta t}\mathbf{\Lambda}_{jk}f_{k}^{ne}+G_{j}+F_{j}+\frac{\Delta t}{2}\bar{D}_{j}F_{j}+O(\epsilon^{3}+\epsilon^{3}/\Delta t).$$

$$\mathbf{f}=(f_0,f_1,\cdots ,f_{q-1}{)}^T,{\mathbf{f}}^{eq}=(f_0^{eq},f_1^{eq},\cdots ,f_{q-1}^{eq}{)}^T,{\mathbf{f}}^{ne}=\mathbf{f}-{\mathbf{f}}^{eq}$$

$$\mathbf{G}=(G_0,G_1,\cdots ,G_{q-1}{)}^T,\mathbf{F}=(F_0,F_1,\cdots ,F_{q-1}{)}^T,$$

$$\bar{\mathbf{D}}=\text{diag}({\bar{D}}_0,{\bar{D}}_1,\cdots ,{\bar{D}}_{q-1}),\mathbb{D}=\text{diag}(D_0,D_1,\cdots ,D_{q-1}),$$

$$\mathbf{L}=\mathbf{L}(\Delta t\mathbb{D})=\sum _{k\ge 1}\frac{(\Delta t\mathbb{D}{)}^k}{k!},$$

$$\mathbf{Lf}=-\mathbf{\Lambda }{\mathbf{f}}^{ne}+\Delta t\tilde{\mathbf{F}}=-\mathbf{\Lambda }(\mathbf{f}-{\mathbf{f}}^{eq})+\Delta t\tilde{\mathbf{F}}.$$

$$\mathbf{f}={\mathbf{f}}^{eq}+\Delta t{\mathbf{\Lambda }}^{-1}\tilde{\mathbf{F}}-{\mathbf{\Lambda }}^{-1}\mathbf{Lf}.$$

$$\bar{\mathbf{L}}={\mathbf{\Lambda }}^{-1}\mathbf{L},\bar{\mathbf{G}}={\mathbf{\Lambda }}^{-1}\tilde{\mathbf{F}},{\bar{\mathbf{f}}}^{eq}={\mathbf{f}}^{eq}+\Delta t{\mathbf{\Lambda }}^{-1}\tilde{\mathbf{F}}={\mathbf{f}}^{eq}+\Delta t\bar{\mathbf{G}},$$

$$\mathbf{f}={\bar{\mathbf{f}}}^{eq}-\bar{\mathbf{L}}\mathbf{f}.$$

$$\mathbf{f}=\sum _{i=0}^{+\infty }(-\bar{\mathbf{L}}{)}^i{\bar{\mathbf{f}}}^{eq}=\sum _{i=0}^k(-\bar{\mathbf{L}}{)}^i{\bar{\mathbf{f}}}^{eq}+O(\Delta t^{k+1})=\sum _{i=0}^k(-\bar{\mathbf{L}}{)}^i{\mathbf{f}}^{eq}+\Delta t\sum _{i=0}^{k-1}(-\bar{\mathbf{L}}{)}^i\bar{\mathbf{G}}+O(\Delta t^{k+1}),$$

$$\mathbf{f}=(\mathbf{I}-\bar{\mathbf{L}}){\mathbf{f}}^{eq}+\Delta t\bar{\mathbf{G}}=(\mathbf{I}-\bar{\mathbf{L}}){\mathbf{f}}^{eq}+\Delta t{\mathbf{\Lambda }}^{-1}(\mathbf{G}+\mathbf{F})+O(\Delta t^2),$$

$$\mathbf{f}=(\mathbf{I}-\bar{\mathbf{L}}+{\bar{\mathbf{L}}}^2){\mathbf{f}}^{eq}+\Delta t(\mathbf{I}-\bar{\mathbf{L}}){\mathbf{\Lambda }}^{-1}(\mathbf{G}+\mathbf{F})+\frac{\Delta t^2}{2}{\mathbf{\Lambda }}^{-1}\bar{\mathbf{D}}\mathbf{F}+O(\Delta t^3),$$

$$\bar{\mathbf{L}}={\mathbf{\Lambda }}^{-1}\mathbf{L}$$

$${\bar{\mathbf{L}}}^2=\Delta t^2({\mathbf{\Lambda }}^{-1}\mathbb{D}{)}^2+O(\Delta t^3).$$

$${\mathbf{f}}^{ne}=\mathbf{f}-{\mathbf{f}}^{eq}=-\Delta t{\mathbf{\Lambda }}^{-1}\mathbb{D}{\mathbf{f}}^{eq}+\Delta t{\mathbf{\Lambda }}^{-1}(\mathbf{G}+\mathbf{F})+O(\Delta t^2),$$

$$\mathbf{f}^{ne}=\mathbf{f}-\mathbf{f}^{eq}=-\Delta t\mathbf{\Lambda}^{-1}\mathbb{D}\mathbf{f}^{eq}+\Delta t\bar{\mathbf{G}}+\Delta t^{2}\mathbf{\Lambda}^{-1}\mathbb{D}\big{[}\big{(}\mathbf{I}-\frac{\mathbf{\Lambda}}{2}\big{)}\Lambda^{-1}\mathbb{D}\mathbf{f}^{eq}-\Lambda^{-1}(\mathbf{G}+\mathbf{F})\big{]}+O(\Delta t^{3}).$$

$${\mathbf{f}}^{ne}=O(\Delta t),$$

$$\mathbb{D}{\mathbf{f}}^{eq}=-\frac{\mathbf{\Lambda }}{\Delta t}{\mathbf{f}}^{ne}+\mathbf{G}+\mathbf{F}+O(\Delta t),$$

$$\mathbb{D}\mathbf{f}^{eq}-\Delta t\mathbb{D}\big{[}\big{(}\mathbf{I}-\frac{\mathbf{\Lambda}}{2}\big{)}\mathbf{\Lambda}^{-1}\mathbb{D}\mathbf{f}^{eq}-\mathbf{\Lambda}^{-1}(\mathbf{G}+\mathbf{F})\big{]}=-\frac{\mathbf{\Lambda}}{\Delta t}\mathbf{f}^{ne}+\mathbf{G}+\mathbf{F}+\frac{\Delta t}{2}\bar{\mathbf{D}}\mathbf{F}+O(\Delta t^{2}).$$

$$\mathbb{D}\mathbf{f}^{eq}+\mathbb{D}\big{(}\mathbf{I}-\frac{\mathbf{\Lambda}}{2}\big{)}\mathbf{f}^{ne}+\frac{\Delta t}{2}\mathbb{D}(\mathbf{G}+\mathbf{F})=-\frac{\mathbf{\Lambda}}{\Delta t}\mathbf{f}^{ne}+\mathbf{G}+\mathbf{F}+\frac{\Delta t}{2}\bar{\mathbf{D}}\mathbf{F}+O(\Delta t^{2}).$$

$$\sum _{l=1}^N\frac{\Delta t^l}{l!}{D}_j^l{f}_j+O(\Delta t^{N+1})=-{\mathbf{\Lambda }}_{jk}{f}_k^{ne}+\Delta t{\tilde{F}}_j,$$

$$f_j^{ne}=O(\Delta t),$$

$$\sum _{l=1}^{N-1}\frac{\Delta t^l}{l!}{D}_j^l(f_j^{eq}+f_j^{ne})+\frac{\Delta t^N}{N!}{D}_j^N{f}_j^{eq}=-{\mathbf{\Lambda }}_{jk}{f}_k^{ne}+\Delta t{\tilde{F}}_j+O(\Delta t^{N+1}).$$

$$\Delta t{D}_j{f}_j^{eq}=-{\mathbf{\Lambda }}_{jk}{f}_k^{ne}+\Delta t(G_j+F_j)+O(\Delta t^2),$$

Full text

1 Introduction

Fluid flows, heat and mass transfer are the most common basic phenomena in nature, science and engineering, and usually, they can be described by the classical Navier-Stokes equations (NSEs) and nonlinear convection-diffusion equations (CDEs) [1]. However, it is difficult to obtain their analytical solutions due to the coupling between the nonlinear NSEs and CDEs. With the development of computing technique, the numerical simulation has been becoming more popular in the study of the complex fluid flows, heat and mass transfer. The lattice Boltzmann (LB) method, as one of mesoscopic numerical approaches, has attained increasing attention, and also gained a great success in the simulation of the complex physical systems for its distinct advantages in kinetic background, treatment of complex boundary conditions, and computational efficiency in parallel systems [2, 3, 4, 5].

Based on the collision term, the LB models developed in the past years can be classified into three different kinds. The first one is the lattice BGK model in which only one relaxation time is introduced into the collision term, hence it is also usually called the single-relaxation-time (SRT) model. Historically, this model, as one of the most popular LB models for its simplicity and efficiency, has been proposed independently by two different research groups [6, 7], and later, it is also found that the SRT model can be considered as a special discretized form of the continuous BGK-Boltzmann equation [8]. Almost at the same time, the second one named multiple-relaxation-time (MRT) or generalized LB model is developed by d Humières [9]. Compared to the SRT model, the collision of MRT model is depicted with a matrix that allows for a decoupled relaxation of different moments, and includes all possible degrees of freedom that can be used to optimize the LB method in terms of stability and accuracy [10, 11, 12, 13, 14, 15]. The last one is two-relaxation-time (TRT) model, which is first developed by Ginzburg [16, 17]. This model, through introducing two different relaxation parameters in the collision term, possesses the flexibility of the MRT model and the simplicity of SRT model, and also retains the computational efficiency [13, 18]. Actually, besides the Ginzburg’s TRT model [16, 17], the lattice kinetic (LK) scheme [19], the modified lattice kinetic (MLK) schemes [20, 21], and the regularized LB (RLB) models [22, 23, 24, 25, 26] can also be viewed as some special TRT models, as pointed out in Refs. [18, 27]. Meng et al. [28] also developed another version of MRT model which combines the commonly used MRT and MLK schemes [9, 10, 19, 20, 21], and found that compared to the popular MRT model [9, 10], their model is more stable. Recently, Zhao et al. [29] proposed a block triple-relaxation-time LB model for nonlinear anisotropic CDEs in which a block diagonal relaxation matrix composed of three blocks is introduced. They also found that that above mentioned RLB models and the MLK schemes are the particular forms of the block triple-relaxation-time LB model. However, one can show the model developed by Zhao et al. [29] is also a particular MRT version.

As discussed above, although many different LB models have been developed in the past decades, all of them can be viewed as the special forms of MRT model. It should also be noted that as a more general version, however, the MRT model still has some issues needed to be addressed. For instance, it is unclear how to choose the eigenvectors of the collision matrix in the MRT model besides the ones corresponding to the conserved variables (e.g., the density and momentum appeared in the NSEs), this is mainly caused by their nonuniqueness. In addition, the analysis of MRT model is usually limited to the specified dimensional space and lattice structure (e.g., D2Q9 lattice [9, 4, 14]), and compared to the SRT model, it seems more difficult to present a general analysis to derive the macroscopic governing equations. Here we would like to point out that Adhikari and Succi have proposed a guideline to choose the eigenvectors through introducing the notion of duality [30]. Furthermore, based on the Taylor expansion, Kaehler and Wagner have performed a general analysis (not restricted to a specific lattice structure) to the MRT model and also gave some necessary requirements on the equilibrium distribution function and the collision matrix, while in their analysis, some assumptions on the space and time derivatives are made [31]. On the other hand, we also note that some different analysis methods (e.g., the Chapman-Enskog analysis [4, 32], Maxwell iteration [33, 34, 35, 36], direct Taylor expansion or the Taylor expansion in Refs. [31, 37, 38], and recurrence equations [39, 40] approaches) used to derive the macroscopic governing equations from the LB models are developed from different physical and mathematical points of view, while it is still unclear whether they are equivalent to each other in obtaining these macroscopic governing equations. In this paper, we will first present a unified framework for the modeling of GLBM. Then we also carry out a comparison among four different analysis methods, and following the previous work [31], a general analysis on the GLBM with some different methods is also performed. To implement the GLBM, some elements are also provided.

2 Generalized lattice Boltzmann method

The semi-discrete evolution equation of the GLBM with the D$d$Q$q$ ($q$ discrete velocities in $d$-dimensional space, $d=1\sim 3$) lattice model is written as (see Appendix A for details)

[TABLE]

where $f_{j}(\mathbf{x},t)$ is the distribution function at position $\mathbf{x}$ in $d$-dimensional space and time $t$ along the velocity $\mathbf{c}_{j}$, $\mathbf{\Lambda}=(\mathbf{\Lambda}_{jk})$ is a $q\times q$ invertible collision matrix, $f_{j}^{ne}(\mathbf{x},t)=f_{j}(\mathbf{x},t)-f_{j}^{eq}(\mathbf{x},t)$ is the non-equilibrium distribution function, and $f_{j}^{eq}(\mathbf{x},t)$ is the equilibrium distribution function. $F_{j}(\mathbf{x},t)$ is the distribution function of a source or forcing term, and $G_{j}(\mathbf{x},t)$ is the auxiliary source distribution function for removing additional terms. $\Delta t$ is the time step, $\bar{D}_{j}=\partial_{t}+\gamma\mathbf{c}_{j}\cdot\nabla$ with $\gamma\in\{0,1\}$ being a parameter to be determined later.

Generally, the evolution process based on Eq. (1) can be divided into two sub-steps, i.e., collision and propagation,

[TABLE]

where $\tilde{f}_{j}(\mathbf{x},t)$ is the post-collision distribution function. It should be noted that the GLBM (1) has the form of popular MRT model with collision matrix $\mathbf{\Lambda}$ [9], and the commonly used LB models for the NSEs and nonlinear convection-diffusion equations (NCDEs) can be viewed as its special cases (see section 5 for details). Here we would also like to point out that in almost all of the MRT models, the collision process in Eq. (1) is carried out in the moment space rather than the velocity (distribution function) space, and the analysis method (e.g., the Chapman-Enskog analysis) is usually more complicated and often depends on the specified lattice structure or discrete velocity set [9, 4, 41, 42, 14]. In this work, we aim to provide a unified framework for the modeling and theoretical analysis of GLBM (1) in the velocity space.

In the GLBM, besides the evolution equation (1), some other key elements including the equilibrium distribution function $f_{j}^{eq}$, distribution function of the source or forcing term $F_{j}$, collision matrix $\mathbf{\Lambda}$ and auxiliary source distribution function $G_{j}$, must also be given properly. In addition, in the implementation of the GLBM, there are two popular schemes that can be used to discrete the last term on the right hand side of Eq. (1), i.e., $\Delta t\bar{D}_{j}F_{j}(\mathbf{x},t)/2$. Actually, if $\gamma=0$, the first-order explicit difference scheme $\partial_{t}F_{j}(\mathbf{x},t)=[F_{j}(\mathbf{x},t+\Delta t)-F_{j}(\mathbf{x},t)]/\Delta t$ is adopted for NCDEs [45, 14]. While for the case of $\gamma=1$, we can use the first-order implicit difference scheme $(\partial_{t}+\mathbf{c}_{j}\cdot\nabla)F_{j}(\mathbf{x},t)=[F_{j}(\mathbf{x}+\mathbf{c}_{j}\Delta t,t+\Delta t)-F_{j}(\mathbf{x},t)]/\Delta t$, and then the linear transform considered in Refs. [43, 44, 14] is applied to avoid the implicitness.

3 The analysis approaches for GLBM: Chapman-Enskog analysis, Maxwell iteration, direct Taylor expansion and recurrence equations

In the literature, there are mainly four basic analysis approaches for LB models, i.e.,

Chapman-Enskog (CE) analysis. 2. 2.

Maxwell iteration (MI) method. 3. 3.

Direct Taylor expansion (DTE) method. 4. 4.

Recurrence equations (RE) method.

In above four approaches, the CE analysis is the most widely used method in the LB community, and combines the physical principle with mathematical technique (Taylor expansion). Different from the CE analysis, however, in the MI, DTE and RE methods, only the Taylor expansion is adopted. In what follows, we will present an introduction to these four analysis methods, and also conduct a comparison among them.

3.1 The Chapman-Enskog analysis

To derive the macroscopic equation from Eq. (1), one can consider the popular CE analysis in which the distribution functions $f_{j}$, $F_{j}$, $G_{j}$, the time and space derivatives can be expanded as [4, 14]

[TABLE]

where $\epsilon$ is a small expansion parameter.

Applying the Taylor expansion to Eq. (1), we have

[TABLE]

where $D_{j}=\partial_{t}+\mathbf{c}_{j}\cdot\nabla$. Denoting $D_{1j}=\partial_{t_{1}}+\mathbf{c}_{j}\cdot\nabla_{1},\bar{D}_{1j}=\partial_{t_{1}}+\gamma\mathbf{c}_{j}\cdot\nabla_{1}$, and substituting Eq. (3) into Eq. (4) yields the following equations at different orders of $\epsilon$,

[TABLE]

Based on the fact that $\mathbf{\Lambda}$ is invertible, Eq. (5a) gives

[TABLE]

With the help of Eq. (5b), we can rewrite Eq. (5c) as

[TABLE]

In addition, multiplying $\epsilon$ on both sides of Eq. (5b) and including some terms at the order of $\epsilon^{2}$, we have

[TABLE]

Through a combination of Eqs. (5b) and (7), i.e., $\epsilon^{1}\times$Eq. (5b)+$\epsilon^{2}\times$Eq. (7), we can also obtain

[TABLE]

3.2 The Maxwell iteration method

To simplify the following discussion on the MI method, we first introduce the following notations,

[TABLE]

then with the Taylor expansion, Eq. (1) can be written as an operator equation,

[TABLE]

where $\tilde{\mathbf{F}}=\mathbf{G}+\mathbf{F}+(\Delta t/2)\bar{\mathbf{D}}\mathbf{F}$.

From Eq. (11), one can obtian

[TABLE]

Denote

[TABLE]

we can write Eq. (12) in a simple form,

[TABLE]

In the MI method [33, 34, 35, 36] where the distribution function $\mathbf{f}$ is substituted into the right hand side of Eq. (14), one can obtain the following formula from Eq. (14),

[TABLE]

which also leads to

[TABLE]

where $\bar{\mathbf{L}}=O(\Delta t)$ has been applied.

In addition, from Eq. (10) one can also derive

[TABLE]

Substituting Eq. (17) into Eq. (16) yields

[TABLE]

Based on Eq. (18) and using $\bar{\mathbf{G}}=\mathbf{\Lambda}^{-1}\tilde{\mathbf{F}}=\mathbf{\Lambda}^{-1}(\mathbf{G}+\mathbf{F})+O(\Delta t)$, we have

[TABLE]

Inserting Eq. (19b) into Eq. (19c) gives rise to the following equation,

[TABLE]

We note that Eqs. (19b) and (20) are the same as Eqs. (8) and (9) in the CE analysis except the truncation errors neglected. Actually, when $\epsilon=O(\Delta t)$, the truncation errors of the CE analysis and MI method are of the same order in $\Delta t$.

3.3 The direct Taylor expansion method

From above discussion, it can be seen clearly that Eqs. (20) and (9) have the similar form. In addition, one can also find that the aim of the MI method is to give an approximate expression of distribution function in terms of $\bar{\mathbf{f}}^{eq}$ or $\mathbf{f}^{eq}$ [see Eq. (15)], which usually makes the MI method a little more complicated than the CE expansion. However, above analysis [see Eqs. (19) and (20)] based on the MI method also implies that it may be possible to get the same results by directly using the Taylor expansion. In the following, we present this simple approach called direct Taylor expansion method, which is similar to those adopted in the previous works [31, 37, 38].

The basic idea of the DTE method is to give a series of equations on $\mathbf{f}^{eq}$ and $\mathbf{f}^{ne}$ with different orders of truncation errors from Eq. (1), which is different from that in the MI method. If we apply the Taylor expansion to Eq. (1), one can get

[TABLE]

Based on the relation $f_{j}=f_{j}^{eq}+f_{j}^{ne}$ and Eq. (21), the following equations are obtained,

[TABLE]

Then from Eq. (22), one can derive the equations at different orders of $\Delta t$,

[TABLE]

or equivalently,

[TABLE]

According to Eq. (23a), we have

[TABLE]

Substituting Eq. (25) into Eq. (24b), one can obtain the following two equivalent equations,

[TABLE]

[TABLE]

It is worth noting that Eqs. (24a) and (3.3) are the same as Eqs. (19b) and (20) which means that mathematically, the DTE method is equivalent to MI method, while the DTE method seems more straightforward and much simpler. In addition, we would also like to point out that different from the DTE method used in Ref. [31, 38], we do not make any assumptions in above analysis.

3.4 The recurrence equations method

Like MI method, the RE method is also an approach based on Taylor expansion. In the RE method [39, 40], a kind of difference equations about $f_{j}^{eq}$ and $f_{j}^{ne}$ named recurrence equations is first constructed from the evolution equation of LB model, then one can apply Taylor expansion to these recurrence equations to obtain the serious equations at different orders of truncation errors.

In order to derive REs, we first rewrite Eq. (1) as

[TABLE]

Let

[TABLE]

then

[TABLE]

which can be used to derive the following equations,

[TABLE]

Substituting Eq. (30) into Eq. (27), one can obtain

[TABLE]

The sum and the difference of Eqs. (31a) and (31b) become [39, 40]

[TABLE]

where the following central difference schemes are adopted for the time and space derivatives of the variable $\phi=\{\bar{f}_{j}^{eq},\bar{g}_{j}\}$,

[TABLE]

Here It should be noted that unlike the previous works limited to cases with the constant coefficients [39, 40], the results of the present RE method can be used for a more general case with variable coefficients.

In addition, we would like to point out that Eqs. (31) and (32) are the so-called REs obtained from Eq. (1). Through using Taylor expansion to the REs, we can also determine a series of equations on $f_{j}^{eq}$ and $g_{j}$ (or $f_{j}^{ne}$) with different orders of time step. Actually, applying the Taylor expansion to Eq. (32b) gives

[TABLE]

Based on the relation $\bar{g}_{j}=O(\Delta t)$ [see Eq. (28)], we have

[TABLE]

Then we can also rewrite Eqs. (35) and (34) as

[TABLE]

where Eqs. (13), (28) and (29) have been used. Keeping in mind that $\tilde{F}_{j}=G_{j}+F_{j}+(\Delta t/2)\bar{D}_{j}F_{j}$, the following equation can be obtained from Eq. (36b),

[TABLE]

It should be noted that Eqs. (36a) and (37) are the same as Eqs. (24a) and (3.3), which also means that the MI, DTE and RE methods are equivalent to each other.

4 Derivation of the macroscopic equations

In this section, we will adopt above analysis methods to derive the macroscopic equations including the nonlinear anisotropic convection-diffusion equation (NACDE) and classical NSEs, and mainly focus on the CE analysis and some of other methods. This is because from the discussion in the previous section, we can see that the CE method has incorporated some physical principles, and compared to some other mathematical methods based Taylor expansion, it also presents some different results at the higher order equations. On the other hand, it is clear that the MI, DTE and RE methods are equivalent to each other since they can give the same equations at the first and second order of time step, thus for brevity, only some of them are considered.

4.1 Derivation of NACDE

The $d$-dimensional NACDE with a source term can be expressed as [14]

[TABLE]

where $\phi$ is a unknown scalar function of position $\mathbf{x}$ and time $t$, $S$ is a scalar source term. $\mathbf{B}=(\mathbf{B}_{\alpha})$ is a vector function, $\mathbf{K}=(\mathbf{K}_{\alpha\beta})$ and $\mathbf{D}=(\mathbf{D}_{\alpha\beta})$ are symmetric tensors (matrices), and they can be functions of $\phi,\mathbf{x}$, and $t$. We note that Eq. (38) can be considered as a general form of CDEs, and many different kinds of CDEs considered in some previous works [16, 45] are its special cases. When $\mathbf{K}_{\alpha\beta}=\kappa\delta_{\alpha\beta}$ and $\mathbf{D}_{\alpha\beta}=h\delta_{\alpha\beta}$ with $\kappa$ and $h$ being two scalar functions, Eq. (38) becomes an isotropic CDE.

To recover Eq. (38) from the GLBM (1), we only need to give some appropriate constraints on the collision matrix $\mathbf{\Lambda}$, the moments of $f_{j}^{eq}(\mathbf{x},t)$, $G_{j}(\mathbf{x},t)$ and $F_{j}(\mathbf{x},t)$. Based on these constrains, one can determine the expressions of $\mathbf{\Lambda}$ and these distribution functions. Actually, to derive Eq. (38), only the zeroth- to second-order moments of $f_{j}^{eq}$, the zeroth- and first-order moments of $G_{j}$ and $F_{j}$, and the constraints on $\mathbf{\Lambda}$ corresponding to these moments are needed. In the GLBM, the unknown conserved scalar $\phi$ can be calculated by $\phi=\sum_{j}f_{j}$, the distribution functions $f_{j}$, $f_{j}^{eq}$, $G_{j}$ and $F_{j}$ should satisfy the following relations,

[TABLE]

where $c_{s}$ is the lattice sound speed related to lattice speed $c=\Delta x/\Delta t$ with $\Delta x$ being the lattice spacing. $\beta$ is a parameter that can be used to adjust the relaxation matrix [see following Eq. (47)], $\mathbf{C}$ is an auxiliary moment to be determined later, $\mathbf{M}_{1G}$ and $\mathbf{M}_{1F}$ are the first-order moments of $G_{j}$ and $F_{j}$.

In addition, for the collision matrix $\mathbf{\Lambda}$, the following requirements are needed,

[TABLE]

where $\mathbf{e}=(1,1,\cdots,1)\in R^{q}$, $\mathbf{S}_{1}$ is an invertible $d\times d$ relaxation matrix corresponding to the diffusion matrix $\mathbf{K}$. Based on Eqs. (39a) and (40), one can find that the relaxation parameter $s_{0}$ does not appear in the recovered macroscopic equation, and its value can be chosen arbitrarily.

4.1.1 Derivation of the NACDE through CE analysis

In this part, we will present some details on how to derive NACDE (38) from GLBM (1). In the CE analysis, we can first express the source term as $S=\epsilon S^{(1)}+\epsilon^{2}S^{(2)}$, then from Eqs. (3), (6) and (39) we have

[TABLE]

where $k\geq 1$. Summing Eqs. (5b) and (7) over $j$ and using Eqs. (39), (40) and (41), one can obtain

[TABLE]

With the aid of Eqs. (5b), (39), (40) and (41), we get

[TABLE]

Substituting above equation into Eq. (42b), we have

[TABLE]

where

[TABLE]

If we take $RH=0$, Eq. (44) would reduce to

[TABLE]

with

[TABLE]

Combining Eqs. (42a) with (46), we obtain the macroscopic NACDE (38).

Here we would like to point out that in the expression of diffusion matrix $\mathbf{K}$ [see Eq. (47)], the part $(-\Delta t\beta c_{s}^{2}\mathbf{I}/2)$ is the numerical diffusion caused by the discrete effect. Actually, when we use the CE analysis to obtain the NACDE from the discrete-velocity Boltzmann equation (see Appendix B for details), the term $(-\Delta t\beta c_{s}^{2}\mathbf{I}/2)$ is not included in the matrix $\mathbf{K}$. In addition, from above CE analysis, we can also see that the functions $\mathbf{C},\mathbf{M}_{1G}$ and $\mathbf{M}_{1F}$ must be chosen properly to make $RH=0$, and in this case, the NACDE (38) can be recovered correctly without additional assumptions. In the following, we present some special cases.

Case 1: $\mathbf{C}=0,\ \mathbf{M}_{1F}=0$. Under the condition of $RH=0$, we have

[TABLE]

which leads to the following $\mathbf{M}_{1G}$,

[TABLE]

Case 2: $\mathbf{B}=\mathbf{B}(\phi),\ \mathbf{M}_{1F}=0$. If $\mathbf{B}$ is a differentiable function of $\phi$, $\mathbf{C}$ and $\mathbf{M}_{1G}$ should satisfy the following equation,

[TABLE]

which leads to

[TABLE]

Furthermore, if we take $\mathbf{C}$ as

[TABLE]

which satisfies

[TABLE]

then it follows from Eq. (42a) that

[TABLE]

Substituting above equation into Eq. (50), one can determine $\mathbf{M}_{1G}$ by

[TABLE]

4.1.2 Derivation of NACDE through MI method

Based on the similarity between CE analysis and MI method, we can also recover NACDE (38) through MI method. To do this, we first introduce the following notations,

[TABLE]

then Eqs. (39) and (40) can be written in the matrix forms,

[TABLE]

According to above equations, we have

[TABLE]

Multiplying $\mathbf{e}$ on both sides of Eqs. (19b) and (20), and using Eq. (58), we have

[TABLE]

Based on Eqs. (57) and (58), Eq. (59) can be written as

[TABLE]

On the other hand, Eq. (19b) gives

[TABLE]

Substituting Eq. (61) into Eq. (60b), one can obtain

[TABLE]

where the diffusion matrix $\mathbf{K}$ is the same as that appeared in Eq. (47), and $RH_{1}$ is defined by

[TABLE]

Under the requirement of $RH_{1}=0$, we can derive the following condition,

[TABLE]

which is also consistent with Eq. (51). Under the condition of Eq. (64), we can correctly recover the NACDE (38) at the order of $\Delta t^{2}$,

[TABLE]

We now present a special discussion on how to calculate the diffusion flux [$-\mathbf{K}\cdot(\nabla\cdot\mathbf{D})$], which is similar to the previous works [46, 47]. Actually, from Eqs. (61) and (64) we have

[TABLE]

From Eq. (68) we can obtain a local scheme for the diffusion flux with a second-order accuracy in time,

[TABLE]

where Eq. (47) has been used. Additionally, if $\mathbf{D}$ is only a function of $\phi$, we can also obtain a local scheme for $\nabla\phi$ from Eq. (68) or (71).

Remark: We noted that following the idea in Ref. [45], Zhang et al. also developed a MRT model for a general isotropic CDE in a recent work [36], and found that under the diffusive scaling ($\Delta t=\eta\Delta x^{2}$, $\eta$ is a parameter) [42], the general isotropic CDE can be recovered correctly from the MRT model. However, it should be noted that the diffusive scaling is an assumption, and also brings an additional constrain on the time step and lattice spacing. Besides, we would also like to point out that under the diffusive scaling, the present model and the one in Ref. [14] would reduce to the MRT model in Ref. [36]. Actually, under the diffusive scaling, the last term on the right hand side of Eq. (1) can be neglected since it is of fourth order of $\Delta x$, and also the terms $\partial_{t}\phi$ and $\partial_{t}\mathbf{B}$ would be absent in Eqs. (60a) and (61). In this case, the auxiliary moment $\mathbf{C}$ and first-order moment $\mathbf{M}_{1G}$ are not needed, and can be set to be zero, which also means that the auxiliary distribution function $G_{i}$ can be ignored in the evolution equation (1). Thus under the diffusive scaling, the evolution equation (1) can reduce to that in Ref. [36].

4.1.3 The equilibrium, auxiliary and source distribution functions of GLBM for NACDE

As discussed above, it can be found that to recover the NACDE (38) correctly, some proper requirements on the equilibrium, auxiliary and source distribution functions should be satisfied. Based on Eq. (39) for a general D$d$Q$q$ lattice model, we can give the following expressions of $f_{j}^{eq}$, $G_{j}$ and $F_{j}$,

[TABLE]

where $\mathbf{M}_{1G}$ is given by Eq. (64).

We note that when $\mathbf{D}=\phi\mathbf{I}$, the NACDE (38) would become

[TABLE]

If we further take $\beta=1$ and $\mathbf{C}=0$, the equilibrium and auxiliary distribution functions $f_{j}^{eq}$ and $G_{j}$ defined by Eqs. (72a) and (72b) can be simplified as

[TABLE]

where the term $\partial_{t}\mathbf{B}$ can be computed by the first-order explicit difference scheme, i.e., $\partial_{t}\mathbf{B}=[\mathbf{B}(\mathbf{x},t+\Delta t)-\mathbf{B}(\mathbf{x},t)]/\Delta t$.

In addition, we would also like to point out that at the diffusive scaling, Eq. (72) can be simplified by

[TABLE]

When $\beta=1$ and $\mathbf{D}=D(\phi)\mathbf{I}$, the equilibrium distribution function (75a) would reduce to the one in Ref. [36].

4.2 Derivation of NSEs

We now consider the following $d$-dimensional NSEs with the source and forcing terms,

[TABLE]

where the shear stress $\mathbf{\sigma}$ is defined by

[TABLE]

where $\mu$ is the dynamic viscosity, $\lambda=\mu_{b}-2\mu/d$ with $\mu_{b}$ being the bulk viscosity [48, 49].

Similar to the derivation of NACDE, to recover NSEs from GLBM (1) we also need to give some constraints on $\mathbf{\Lambda}$, $f_{j}$, $f_{j}^{eq}$, $G_{j}$, and $F_{j}$. In addition, compared to Eqs. (39), (40) and (41), there is another requirement on the higher-order moments of the distribution functions for NSEs. Here the following conditions should be satisfied,

[TABLE]

where $\mathbf{M}_{2G}$ is a second-order tensor to be determined below, $\Delta_{\alpha\beta\gamma\theta}=\delta_{\alpha\beta}\delta_{\gamma\theta}+\delta_{\alpha\gamma}\delta_{\beta\theta}+\delta_{\alpha\theta}\delta_{\beta\gamma}$, $\mathbf{S}_{1}$ is an invertible $d\times d$ relaxation matrix, $S_{2s}$ and $S_{2b}$ are the relaxation parameters corresponding to the second-order moments.

4.2.1 Derivation of NSEs through DTE method

Due to the similarity between the CE analysis and and some other methods, here we only use the DTE method, i.e., Eqs. (24a) and (3.3), to recover NSEs. With the help of Eqs. (28) and (29), we can rewrite Eqs. (24a) and (3.3) as

[TABLE]

On the other hand, from Eqs. (28) and (78) we have

[TABLE]

and

[TABLE]

Summing Eq. (79b) and adopting Eqs. (78) and (80), we can obtain

[TABLE]

If we take $\gamma=1$, then the continuity equation (76a) is recovered correctly to the order of $O(\Delta t^{2})$.

Multiplying $\mathbf{c}_{j}$ on both sides of Eq. (79), and through a summation over $j$, we have

[TABLE]

where Eqs. (78), (80) and (81) have been used.

Additionally, from Eqs. (78) and (79) we get

[TABLE]

Based on the following equations,

[TABLE]

we can rewrite Eq. (84) as

[TABLE]

where the term $\nabla_{\gamma}\cdot(\rho\mathbf{u}_{\alpha}\mathbf{u}_{\beta}\mathbf{u}_{\gamma})$ have been neglected, this is because for incompressible fluid flows, it is the order of $O(Ma^{3})$ with $Ma$ being the Mach number. Substituting Eq. (86) into Eq. (83b) and using Eq. (82), we can obtain

[TABLE]

where

[TABLE]

To obtain the correct NSEs, $(RH_{2})_{\alpha\beta}=0$, which gives rise to following equations,

[TABLE]

Introduce the matrix $\mathbf{W}=\bar{\mathbf{F}}\mathbf{u}+\mathbf{u}\bar{\mathbf{F}}+(c_{s}^{2}\mathbf{I}-\mathbf{u}\mathbf{u})\bar{S}$, then Eq. (89b) can be written as

[TABLE]

where $\textbf{tr}(\mathbf{W})$ is the trace of the matrix $\mathbf{W}$. We note that when $\bar{S}=0$ and $S_{2s}=S_{2b}$, Eq. (90) would be the same as the second-order moment of the forcing term distribution function in the previous work [50].

Based on $\mathbf{M}_{2G}$ defined by Eq. (90), we can write Eq. (87) in the following form,

[TABLE]

When $p=\rho c_{s}^{2}$ and the truncation error $O(\Delta t^{2})$ is neglected, Eq. (91) would reduce to the momentum equation (76b) with the following dynamic and bulk viscosities,

[TABLE]

Finally, we also present a local scheme to calculate shear stress or strain rate tensor in the framework of LBM [41, 51, 52]. To this end, from Eq. (86a) we can first obtain the following expression of strain rate tensor with a first-order accuracy in time,

[TABLE]

Then based on Eq. (93), one can also give a local scheme for shear stress with a second-order accuracy in time,

[TABLE]

We note that Eqs. (93) and (94) are the same as the results in Ref. [41] once $\bar{S}=0$.

4.2.2 The equilibrium, auxiliary and source distribution functions of GLBM for NSEs

From above discussion, one can clearly observe that to recover the macroscopic NSEs (76) from GLBM (1), the equilibrium, auxiliary and source distribution functions should satisfy some necessary requirements, as depicted by Eq. (78). For a general D$d$Q$q$ lattice model, the explicit expressions of $f_{j}^{eq}$, $G_{j}$ and $F_{j}$ can be given by

[TABLE]

where $\mathbf{M}_{2G}$ is determined by Eq. (90).

5 Collision matrix and some special cases of GLBM

When the discrete velocity set $V_{q}=\{\mathbf{c}_{j},0\leq j\leq q-1\}$ is given, we can construct different forms of collision matrix $\mathbf{\Lambda}$ in GLBM. In general, the collision matrix $\mathbf{\Lambda}$ in the commonly used MRT model has the following form [4, 5],

[TABLE]

where the rows of transformation matrix $\mathbf{M}$ are made up of discrete velocities in $V_{q}$, $\mathbf{S}$ is a block diagonal matrix, $\mathbf{S}_{k}\in R^{n_{k}\times n_{k}}$ is a relaxation matrix corresponding to the $k$-th ($0\leq k\leq m$) order moment of discrete velocity space. Denote $\mathbf{M}$ with the following form,

[TABLE]

then we have

[TABLE]

with

[TABLE]

For the GLBM (1) discussed above, we only need the first two left-eigenvectors of the collision matrix for NACDE,

[TABLE]

and the other eigenvectors can be chosen arbitrarily. However, for the Navier-Stokes equations, only the first three left-eigenvectors are needed,

[TABLE]

where $\mathbf{M}_{2}^{(1)}=(c_{0\alpha}c_{0\alpha},c_{1\alpha}c_{1\alpha},\cdots,c_{q-1,\alpha}c_{q-1,\alpha})$ and $\mathbf{M}_{2}^{(2)}=(c_{0\alpha}c_{0\beta},c_{1\alpha}c_{1\beta},\cdots,c_{q-1,\alpha}c_{q-1,\beta})_{\alpha<\beta}$. Based on Eq. (78f), we can determine the relaxation matrix $\mathbf{S}_{2}$ as

[TABLE]

where $\mathbb{E}_{d}$ is a $d\times d$ matrix with all elements $(\mathbb{E}_{d})_{ij}=1$. Here it should also be noted that for the orthogonal eigenvectors, we can obtain the similar results as those in the previous section, and refer the reader to Ref. [31] for some details.

We now point out that some commonly used LB models in the literature can be obtained from our GLBM when the matrices $\mathbf{\Lambda}$, $\mathbf{M}$ and $\mathbf{S}$ are specified with some special forms.

Lattice BGK or SRT model [6, 7]:

[TABLE] 2. 2.

TRT model [16, 17]:

[TABLE]

Let

[TABLE]

the collision matrix $\mathbf{\Lambda}$ can be written as

[TABLE]

where $s^{+}+s^{-}=2\omega,s^{+}-s^{-}=2\bar{\omega}$, or equivalently, $s^{+}=\omega+\bar{\omega},s^{-}=\omega-\bar{\omega}$ with $s^{+}$ and $s^{-}$ being the relaxation parameters corresponding to the symmetric and anti-symmetric parts of distribution functions.

In addition, we also note that if $\mathbf{c}_{0}=\mathbf{0}$ is not contained in the discrete velocity set $\mathbf{V}_{q}$, the matrix $\mathbf{\Lambda}$ in Eq. (113) would be given by

[TABLE] 3. 3.

RLB models [22, 23, 24, 26]: For simplicity, we first introduce the following tensors,

[TABLE]

In the RLB model, the collision matrix $\mathbf{\Lambda}$ is taken as $\mathbf{\Lambda}=\mathbf{I}-(1-1/\tau)\mathbf{R}$ for NACDE, while for NSEs, $\mathbf{\Lambda}=\mathbf{I}-(1-1/\tau)\mathbf{P}$. It is clear that as a general form, we can take $\mathbf{\Lambda}=\mathbf{I}-(1-1/\tau)(\mathbf{R+P})$ for both NACDE and NSEs. 4. 4.

MLK schemes [20, 21]: For the NACDE (38), the collision matrix $\mathbf{\Lambda}$ in the MLK scheme is given by

[TABLE]

while for the NSEs (76), it should be taken as

[TABLE]

As a general case, the collision matrix $\mathbf{\Lambda}$ can be specified by

[TABLE]

where $\mathbf{R}$ and $\mathbf{P}$ are the same as those in Eq. (113). 5. 5.

Classical MRT models [9, 10, 11]: In the classical MRT models, the collision matrix $\mathbf{\Lambda}$ can be determined by

[TABLE]

where $\mathbf{S}=\textbf{diag}(S_{0},\ S_{1},\ S_{2},\ \cdots,\ S_{q-1})$ is a standard diagonal relaxation matrix, $\mathbf{M}$ is the transformation matrix and is composed of the orthogonal eigenvectors [4, 5]. 6. 6.

Block triple-relaxation-time LB model [29]: For brevity we first introduce the following $\bar{\mathbf{R}}$ and $\bar{\mathbf{P}}$ in the block triple-relaxation-time LB model,

[TABLE]

where $\hat{\mathbf{I}}$ is a matrix with $\hat{\mathbf{I}}_{ij}=1$, the symbol $\circ$ represents the Hadamard product, $\mathbf{K}_{2}$ is a tensor related to $\mathbf{S}_{2}$ [29]. Then the collision matrix $\mathbf{\Lambda}$ can be given by

[TABLE]

6 Conclusions

In this work, we developed a general framework for the modeling of the GLBM for the NACDE and NSEs, and also presented a detailed analysis on the C-E analysis, MI, DTE and RE approaches that have been adopted to derive the macroscopic equations from the LB models. The results show that mathematically, these four analysis methods are equivalent to each other, and can obtain the same macroscopic governing equations from present GLBM. In addition, we provided some details on the equilibrium, auxiliary and source distribution functions, and some special discussion on how to develop the local schemes for the diffusion flux (and/or $\nabla\phi$) and shear stress (and/or stain rate tensor). Finally, we also pointed out that the existing LB models, including the lattice BGK or SRT model, TRT model, RLB model, MLK Schemes, the classical MRT model and the block triple-relaxation-time LB model, are some special cases of the present GLBM.

Acknowledgements

This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 51836003 and 51576079) and the National Key Research and Development Program of China (Grant No. 2017YFE0100100).

7 Appendix

7.1 Appendix A: From discrete-velocity Boltzmann equation to lattice Boltzmann equation

The GLBM [Eq. (1)] can be obtained from the discrete-velocity Boltzmann equation (DVBE) with invertible collision matrix $\tilde{\mathbf{\Lambda}}$,

[TABLE]

which can be also considered as a basic equation to develop some other mesoscopic methods (e.g., the discrete-unified gas-kinetic scheme [53, 54]).

Compared to Eq. (1), the source distribution function in DVBE (115) is the sum of $F_{j}$ and $\tilde{G}_{j}$, $F_{j}$ is used as the source or forcing term, $\tilde{G}_{j}$ is adopted to remove some additional terms.

Integrating Eq. (115) along the characteristic line $\mathbf{x}^{\prime}=\mathbf{x}+\mathbf{c}_{j}t^{\prime}$ with $t^{\prime}\in[0,\Delta t]$, we obtain

[TABLE]

Based on the following results,

[TABLE]

we can write Eq. (7.1) as

[TABLE]

where

[TABLE]

Actually, Eq. (118) can also be written as

[TABLE]

From above equation, we have

[TABLE]

Let

[TABLE]

then

[TABLE]

According to the following equation,

[TABLE]

where Eq. (121) has been used, we can write Eq. (122) as

[TABLE]

Substituting above equation into Eq. (117) and removing the term $O(\Delta t^{3})$, we obtain

[TABLE]

where \mathbf{G}=\big{(}\mathbf{I}-\mathbf{\Lambda}/2\big{)}\tilde{\mathbf{G}}, the relaxation matrix $\mathbf{\Lambda}$ is given by

[TABLE]

We note that Eq. (125) has the same form as Eq. (1) except the last term on the right side. Actually if $\gamma=1$, Eq. (1) would reduce to Eq. (125), while $\gamma=0$ can also be used derive correct NACDE [14].

7.2 Appendix B: From discrete-velocity Boltzmann equation to the macroscopic equations through CEM

It is known from above discussion that two classes of macroscopic equations, i.e., NACDE and NSEs, can be recovered correctly from GLBM (1) through some analysis methods (CE analysis, MI, DTE and RE methods) when the moments of the distribution functions $f_{j}^{eq}$, $G_{j}$ and $F_{j}$ are given properly. Here only the CE analysis is considered to obtain the NACDE and NSEs.

Based on the CE analysis, we have

[TABLE]

Substituting Eq. (127) into Eq. (115), we can obtain the following equations at the order of $\epsilon^{0}$, $\epsilon^{1}$ and $\epsilon^{2}$,

[TABLE]

It can be seen that Eq. (128) derived from DVBE (115) is simpler than Eq. (5) obtained from GLBM (1) in which some discrete terms [e.g., $\Delta t^{2}D_{1j}f_{j}^{(0)}/2$ in Eq. (5c)] caused by the numerical scheme are included.

7.2.1 From DVBE to NACDE through CE analysis

To correctly recover NACDE (38) from the DVBE (115), the unknown conserved scalar $\phi$, $f_{j}$ and $f_{j}^{eq}$ satisfy Eq. (39), while moments of $\tilde{G}_{j}$, $F_{j}$ and $\tilde{\Lambda}$ are given by

[TABLE]

where the moment $\mathbf{M}_{1\tilde{G}}$ to be determined below. $\tilde{\mathbf{S}}_{1}$ is an invertible $d\times d$ relaxation matrix corresponding to the diffusion matrix $\mathbf{K}$.

Similar to Eq. (41), it follows from Eqs. (39), (127), (128a) and (129) that

[TABLE]

where $k\geq 1$. Summing Eq. (128b) and Eq. (128c) over $j$, one can obtain

[TABLE]

where Eqs. (39), (129) and (130) have been used.

Multiplying $\mathbf{c}_{j}$ on both sides of Eqs. (128b) and (128c), and summing them over $j$, we have

[TABLE]

where Eqs. (129) and (130) have been adopted to derive above equation.

Substituting Eq. (132a) into Eq. (131b) gives

[TABLE]

where $\mathbf{K}=\beta c_{s}^{2}\tilde{\mathbf{S}}_{1}^{-1}$ with $\tilde{\mathbf{S}}_{1}^{-1}=(\mathbf{S}_{1}^{-1}-\frac{1}{2}\mathbf{I})\Delta t$, RH_{1}=\tilde{\mathbf{S}}_{1}^{-1}\big{(}\partial_{t_{1}}\mathbf{B}+\nabla_{1}\cdot\mathbf{C}-\mathbf{M}_{1\tilde{G}}^{(1)}\big{)}. Taking $\mathbf{M}_{1\tilde{G}}$ such that

[TABLE]

we have $RH_{1}=0$, which leads to

[TABLE]

Combining Eq. (132a) with Eq. (135) yields the NACDE (38),

[TABLE]

In addition, it should be noted that from Eq. (134) and the relation \mathbf{G}=\big{(}\mathbf{I}-\mathbf{\Lambda}/2\big{)}\tilde{\mathbf{G}}, we can determine $\mathbf{M}_{1G}$ as

[TABLE]

which is the same as Eq. (51).

7.2.2 From DVBE to NSEs through CEM

Similar to above discussion, to recover NSEs (76) from the DVBE (115), the macroscopic variables $\rho$ and $\mathbf{u}$, and the distribution functions $f_{j}$ and $f_{j}^{eq}$ satisfy Eqs. (78a) and (78b), while the moments of $\tilde{G}_{j}$, $F_{j}$ and $\tilde{\mathbf{\Lambda}}$ should satisfy the following relations,

[TABLE]

where the moment $\mathbf{M}_{2\tilde{G}}$ will be determined later. $\tilde{S}_{2s}$ and $\tilde{S}_{2b}$ are two relaxation parameters related to dynamic viscosity $\mu$ and bulk viscosity $\mu_{b}$.

Similarly, from Eqs. (78), (127), (128a) and (138) we have

[TABLE]

where $k\geq 1$. Summing Eqs. (128b) and (128c) over $j$, and using Eqs. (78a), (138) and (139), we can obtain

[TABLE]

Multiplying $\mathbf{c}_{j}$ and $\mathbf{c}_{j}\mathbf{c}_{j}$ on both sides of Eqs. (128b) and (128c), and summing these equations over $j$, one can obtain

[TABLE]

[TABLE]

where Eqs. (138) and (139) have been used.

From Eq. (140), the continuity equation (76a) is recovered correctly to the order of $O(\epsilon^{2})$,

[TABLE]

On the other hand, with the help of Eq. (142), we can obtain the momentum equation (76b) from Eq. (141). According to Eqs. (138a) and (141a), we can rewrite Eq. (142a) as

[TABLE]

To derive correct momentum equation, the following conditions should be satisfied,

[TABLE]

Then from Eq. (145), we can obtain

[TABLE]

Substituting Eqs. (146) into Eq. (141b) yields the following equation,

[TABLE]

Combining Eq. (141a) with Eq. (147) we have

[TABLE]

where the pressure $P$, the dynamic and bulk viscosities are given by

[TABLE]

From Eqs. (143) and (148), it is clear that through the CE analysis, the NSEs (76) can be correctly recovered from the DVBE (115) at the order of $O(\epsilon^{2})$. In addition, we would also like to point out that the relaxation parameters $\tilde{S}_{2s}$ and $\tilde{S}_{2b}$ in the DVBE are related to $S_{2s}$ and $S_{2b}$ in the GLBM through $\tilde{S}_{2s}^{-1}=(S_{2s}^{-1}-1/2)\Delta t$ and $\tilde{S}_{2b}^{-1}=(S_{2b}^{-1}-1/2)\Delta t$.