# Multiple-distribution-function lattice Boltzmann method for convection-
diffusion-system based incompressible Navier-Stokes equations

Zhenhua Chai11footnotemark: 122footnotemark: 2<EMAIL_ADDRESS>Baochang
Shi<EMAIL_ADDRESS>Chengjie Zhan33footnotemark: 3 School of Mathematics
and Statistics, Huazhong University of Science and Technology, Wuhan, 430074,
China  Hubei Key Laboratory of Engineering Modeling and Scientific Computing,
Huazhong University of Science and Technology, Wuhan 430074, China

###### Abstract

In this paper, a multiple-distribution-function lattice Boltzmann method (MDF-
LBM) with multiple-relaxation-time model is proposed for incompressible
Navier-Stokes equations (NSEs) which are considered as the coupled convection-
diffusion equations (CDEs). Through direct Taylor expansion analysis, we show
that the Navier-Stokes equations can be recovered correctly from the present
MDF-LBM, and additionally, it is also found that the velocity and pressure can
be directly computed through the zero and first-order moments of distribution
function. Then in the framework of present MDF-LBM, we develop a locally
computational scheme for the velocity gradient where the first-order moment of
the non-equilibrium distribution is used, this scheme is also extended to
calculate the velocity divergence, strain rate tensor, shear stress and
vorticity. Finally, we also conduct some simulations to test the MDF-LBM, and
find that the numerical results not only agree with some available analytical
and numerical solutions, but also have a second-order convergence rate in
space.

###### keywords:

Multiple-distribution-function lattice Boltzmann method, incompressible
Navier-Stokes equations , convection-diffusion equations , direct Taylor
expansion analysis

##  1 Introduction

In the past three decades, the lattice Boltzmann method (LBM), as a discrete
numerical approach to the Boltzmann equation, has attained a great success in
the study of the complex fluid flows (see some review articles [1, 2, 3, 4, 5]
and monographs [6, 7, 8, 9, 10]), including the multicomponent and multiphase
flows [2, 5, 8, 9, 10], thermal flows [7, 8], turbulent flows [1, 2, 6, 7],
particle suspensions [3], microfluidics [3, 8], flows in porous media [2, 4,
7], to name but a few.

Recently, the vectorial lattice Boltzmann method (VLBM) based on vectorial
kinetic model (e.g., the general BGK model [11, 12, 13]) has also been
developed for the shallow water equations [14] and incompressible Navier-
Stokes equations (NSEs) [15, 16, 17]. Actually, for the isothermal and
incompressible fluid flows where the density is assumed to be a constant, the
NSEs can be viewed as a convection-diffusion system consisted of ($d+1$)
equations in $d$-dimensional space [named convection-diffusion-system based
NSEs, see the following Eq. (2)]. The main idea of the VLBM is to construct a
single evolution equation of distribution function for each convection-
diffusion equation (CDE). Like the double-distribution-function LBM for
thermal fluid flows [18, 19, 20], the VLBM can also be considered as a special
multiple-distribution-function lattice Boltzmann method (MDF-LBM), which would
be used in the present work.

Compared to the popular scalar or single LBM for the NSEs, the MDF-LBM (or
VLBM) for convection-diffusion-system based NSEs has some distinct features.
The first is that one can use fewer discrete velocities to construct the MDF-
LBM. For instance, if we consider the two-dimensional problems, the LBM with
D2Q4 or D2Q5 (four or five discrete velocities in two-dimensional space)
lattice model is enough for the CDEs [21], while in the single LBM for the
NSEs, the D2Q9 lattice model is usually adopted since the high-order isotropy
of the discrete velocities is needed [22]. The second is that in the MDF-LBM
for convection-diffusion-system based NSEs, it is more flexible and much
easier to construct the equilibrium distribution function such that the CDEs
can be recovered correctly. And the third is that in the MDF-LBM for
convection-diffusion-system based NSEs, some physical variables, for example,
the velocity gradient, velocity divergence, the strain rate tensor, the shear
stress and the vorticity, can be calculated locally through the first-order
moments of the non-equilibrium distribution function (see Section III for
details), while in the commonly used single LBM for NSEs, the second-order
moments of the non-equilibrium distribution function are needed to compute
strain rate tensor and shear stress [23, 24, 25, 26], and what is more, it
seems more difficult to directly develop the local scheme for velocity
gradient or vorticity [27].

Inspired by the VLBM for the NSEs [15, 16, 17], in this work we developed a
MDF-LBM for the convection-diffusion-system based incompressible NSEs.
Compared to the previous work [15, 16, 17], however, there are four main
differences: (1) we focus on the MDF-LBM for incompressible NSEs, and the
compressible effect in Refs. [15, 16, 17] is neglected; (2) we propose a
special formula with the first-order moment of the distribution function for
the pressure, which is not only different from those considered in the
previous work [15, 16, 17], but also more strict theoretically; (3) the
formula for the pressure is consistent with the continuity equation, and thus
the lattice Boltzmann model for the continuity equation can be omitted; (4) we
develop some local schemes for velocity gradient, velocity divergence, stain
rate tensor, shear stress and vorticity, which have not been presented and
discussed in the previous works [15, 16, 17].

The rest of the paper is organized as follows. In section 2, we develop a MDF-
LBM for convection-diffusion-system based incompressible NSEs, and then the
direct Taylor expansion of the present MDF-LBM is conducted in section 3. In
section 4, we present some numerical results and discussion, and finally, some
conclusions are given in section 5.

##  2 Multiple-distribution-function lattice Boltzmann method for
incompressible Navier-Stokes equations

In this section, we would first write the incompressible NSEs as a coupled
convection-diffusion system, and then present a MDF-LBM for the convection-
diffusion-system based incompressible NSEs.

###  2.1 The convection-diffusion-system based incompressible Navier-Stokes
equations

For incompressible fluid flows where the density $\rho$ is assumed to be a
positive constant $\rho_{0}$, the NSEs can be written as [28, 29]

$\nabla\cdot\mathbf{u}=R,$ (1a) $\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot(\mathbf{u}\mathbf{u})=-\nabla
P+\nabla\cdot(\nu\nabla\mathbf{u})+\mathbf{F},$ (1b)

where $\mathbf{u}=(u_{\alpha})_{\alpha=1-d}$ is the velocity in
$d$-dimensional space, $R$ is a source term, $P$ is the pressure, $\nu$ is the
kinematic viscosity and $\mathbf{F}=(F_{\alpha}$) is the external force. We
would like to point out that above NSEs can also be reformulated as a coupled
convection-diffusion system,

$\frac{\partial\bar{u}_{\alpha}}{\partial
t}+\nabla\cdot(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})=\nabla\cdot(\nu\nabla\bar{u}_{\alpha})+\bar{F}_{\alpha},\
\alpha=0,1,\cdots,d$ (2)

where $\bar{\mathbf{u}}=(\bar{u}_{\alpha})$ with $\bar{u}_{0}=\rho_{0}$ and
$\bar{u}_{\alpha}=u_{\alpha}\ (\alpha=1-d)$,
$\bar{\mathbf{F}}=(\bar{F}_{\alpha})$ with $\bar{F}_{0}=R$ and
$\bar{F}_{\alpha}=F_{\alpha}(\alpha=1-d)$, $\mathbf{E}_{0}=\mathbf{0}$ and
$\mathbf{E}_{\alpha}(\alpha=1-d)$ is the unit vector in $d$-dimensional space.

It is clear that Eq. (2) is composed of ($d+1$) convection-diffusion equations
(CDEs), and in the following, it is considered as the convection-diffusion-
system based NSEs. We note that although the incompressible NSEs (1) are
equivalent to Eq. (2) mathematically, while the latter is just a convection-
diffusion system, and can be solved more efficiently in the LBM. Actually, in
the LBM for the CDE, only the second-order moment isotropy is needed, which
brings more flexibility in the development of the lattice Boltzmann (LB)
models; while in the single LBM for NSEs, we need the fourth-order moment
isotropy, which would give rise to more limitations in the design of the LB
models.

###  2.2 The multiple-distribution-function lattice Boltzmann method for the
convection-diffusion-system based Navier-Stokes-equations

In the LBM, the LB models can be classified into three basic kinds, i.e., the
single-relaxation-time LB (SRT-LB) model or lattice BGK model [22, 30], the
two-relaxation-time LB (TRT-LB) model [31, 32] and the multiple-relaxation-
time LB (MRT-LB) model [33, 34]. In this work, we would consider the MRT-LB
model for its generality and accuracy [35, 36, 37, 38]. In the past years,
some different MRT-LB models were developed for the isotropic and anisotropic
CDEs [37, 39, 40, 41, 42]. Recently, Chai and Shi proposed a unified framework
for the modeling of the MRT-LB models for the NSEs and nonlinear CDEs [43].
Following this work and inspired by the VLBM [14, 15, 16, 17], The evolution
equation of the MDF-LBM for the convection-diffusion-system based NSEs (2) can
be written as

$f_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta t,t+\delta
t)=f_{i,\alpha}(\mathbf{x},t)-\mathbf{\Lambda}_{ik}\big{[}f_{k,\alpha}(\mathbf{x},t)-f_{k,\alpha}^{eq}(\mathbf{x},t)\big{]}+\delta
t\big{[}G_{i,\alpha}(\mathbf{x},t)+F_{i,\alpha}(\mathbf{x},t)+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)\big{]},$ (3)

where $f_{i,\alpha}(\mathbf{x},t)$ is the distribution function corresponding
to the variable $\bar{u}_{\alpha}$ at position $\mathbf{x}$ and time $t$ along
the discrete velocity $\mathbf{c}_{i}$. $\delta t$ is the time step,
$\mathbf{\Lambda}=(\mathbf{\Lambda}_{ik})$ is a $q\times q$ invertible
collision matrix with $q$ representing the number of discrete velocities.
$\bar{D}_{i}=\partial_{t}+\gamma\mathbf{c}_{i}\cdot\nabla$ with the parameter
$\gamma\in\\{0,1\\}$, which can be discretized by some different first-order
difference schemes [43]. $f_{i,\alpha}^{eq}(\mathbf{x},t)$ is the equilibrium
distribution function, $G_{i,\alpha}(\mathbf{x},t)$ is the auxiliary
distribution function, $F_{i,\alpha}(\mathbf{x},t)$ is the distribution
function of the source term $\bar{F}_{\alpha}$, to derive the correct
macroscopic equation (2), they should be defined by [43]

$f_{i,\alpha}^{eq}=\omega_{i}\big{[}\bar{u}_{\alpha}+\frac{\mathbf{c}_{i}\cdot(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})}{c_{s}^{2}}\big{]},$
(4)
$G_{i,\alpha}=\big{(}1-\frac{s_{1}}{2}\big{)}\omega_{i}\frac{\mathbf{c}_{i}\cdot\partial_{t}(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})}{c_{s}^{2}},$
(5) $F_{i,\alpha}=\omega_{i}\bar{F}_{\alpha},$ (6)

where the simple linear equilibrium distribution function (4) with respect to
the discrete velocity is considered. $\omega_{i}$ is the weight coefficient,
$s_{1}$ is the relaxation parameter corresponding to the kinematic viscosity
$\nu$. $c_{s}$ is the lattice sound speed related to lattice speed $c=\delta
x/\delta t$ with $\delta x$ being the lattice spacing. In the DdQq (q discrete
velocities in $d$-dimensional space) lattice model, the discrete velocity
$\mathbf{c}_{i}$, the weight coefficient $\omega_{i}$ and the lattice sound
speed $c_{s}$ should satisfy the following relations,

$\sum_{i}\omega_{i}=1,$ (7a) $\sum_{i}\omega_{i}\mathbf{c}_{i}=\mathbf{0},$
(7b) $\sum_{i}\omega_{i}\mathbf{c}_{i}\mathbf{c}_{i}=c_{s}^{2}\mathbf{I},$
(7c)

where $\mathbf{I}$ is the unit matrix. Here we list some special cases that
have been widely used in the LBM.

D1Q2 lattice model:

$\mathbf{c}_{i}=(1,-1)c,$ (8a) $\omega_{1}=\omega_{2}=\frac{1}{2},$ (8b)
$c_{s}=c.$ (8c)

D1Q3 lattice model:

$\mathbf{c}_{i}=(0,1,-1)c,$ (9a) $\omega_{0}=\frac{2}{3},\
\omega_{1}=\omega_{2}=\frac{1}{6},$ (9b) $c_{s}=\frac{c}{\sqrt{3}}.$ (9c)

D2Q4 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{cccc}1&0&-1&0\\\
0&1&0&-1\end{array}\right)c$ (10a)
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{4}=\frac{1}{4},$ (10b)
$c_{s}=\frac{c}{\sqrt{2}}.$ (10c)

D2Q5 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{ccccc}0&1&0&-1&0\\\
0&0&1&0&-1\end{array}\right)c$ (11a) $\omega_{0}=\frac{1}{3},\
\omega_{1}=\omega_{2}=\omega_{3}=\omega_{4}=\frac{1}{6},$ (11b)
$c_{s}=\frac{c}{\sqrt{3}}.$ (11c)

D2Q9 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{ccccccccc}0&1&0&-1&0&1&-1&-1&1\\\
0&0&1&0&-1&1&1&-1&-1\end{array}\right)c$ (12a) $\omega_{0}=\frac{4}{9},\
\omega_{1-4}=\frac{1}{9},\ \omega_{5-8}=\frac{1}{36},$ (12b)
$c_{s}=\frac{c}{\sqrt{3}}.$ (12c)

D3Q6 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{cccccc}1&-1&0&0&0&0\\\ 0&0&1&-1&0&0\\\
0&0&0&0&1&-1\end{array}\right)c$ (13a) $\omega_{1-6}=\frac{1}{6},$ (13b)
$c_{s}=\frac{c}{\sqrt{3}}.$ (13c)

D3Q7 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{ccccccc}0&1&-1&0&0&0&0\\\
0&0&0&1&-1&0&0\\\ 0&0&0&0&0&1&-1\end{array}\right)c$ (14a)
$\omega_{0}=\frac{1}{4},\ \omega_{1-6}=\frac{1}{8},$ (14b)
$c_{s}=\frac{c}{2}.$ (14c)

D3Q15 lattice model:

$\mathbf{c}_{i}=\left(\begin{array}[]{ccccccccccccccc}0&1&-1&0&0&0&0&1&1&1&-1&-1&-1&1&-1\\\
0&0&0&1&-1&0&0&1&1&-1&1&-1&1&-1&-1\\\
0&0&0&0&0&1&-1&1&-1&1&1&1&-1&-1&-1\end{array}\right)c$ (15a)
$\omega_{0}=\frac{2}{9},\ \omega_{1-6}=\frac{1}{9},\
\omega_{7-14}=\frac{1}{72}$ (15b) $c_{s}=\frac{c}{\sqrt{3}}.$ (15c)

We would like to point out that different lattice models can be applied for
different evolution equations represented by different $\alpha$. In the MDF-
LBM for the convection-diffusion system-based NSEs, the macroscopic variable
$\bar{u}_{\alpha}$ is computed by

$\bar{u}_{\alpha}=\sum_{i}f_{i,\alpha}.$ (16)

##  3 The direct Taylor expansion of the multiple-distribution-function
lattice Boltzmann method

Historically, there are some asymptotic analysis methods that can be used to
derive the macroscopic governing equation (2) from the MDF-LBM [43], including
the Chapman-Enskog analysis [2, 44], Maxwell iteration method [45], direct
Taylor expansion method [46], recurrence equations method [47]. However, it
has been shown that at the second order of expansion parameters, these four
analysis methods can give the same macroscopic equations [43]. For this
reason, we only consider the direct Taylor expansion (DTE) method for its
simplicity, and additionally, compared to the commonly used Chapman-Enskog
analysis, this method only includes a single expansion parameter $\delta t$.

###  3.1 The direct Taylor expansion

Based on the previous works [43, 46], when the Taylor expansion is applied to
Eq. (3), we have

$\sum_{j=1}^{N}\frac{\delta t^{j}}{j!}D_{i}^{j}f_{i,\alpha}+O(\delta
t^{N+1})=-\mathbf{\Lambda}_{ik}(f_{k,\alpha}-f_{k,\alpha}^{eq})+\delta
t\big{[}G_{i,\alpha}(\mathbf{x},t)+F_{i,\alpha}(\mathbf{x},t)+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)\big{]},$ (17)

where $D_{i}=\partial_{t}+\mathbf{c}_{i}\cdot\nabla$. Introduce
$f_{i,\alpha}^{ne}=f_{i,\alpha}-f_{i,\alpha}^{eq}$ and substitute it into the
collision term of Eq. (17), we can derive the following equations,

$f_{i,\alpha}^{ne}=O(\delta t),$ (18a) $\sum_{j=1}^{N-1}\frac{\delta
t^{j}}{j!}D_{i}^{j}(f_{i,\alpha}^{eq}+f_{i,\alpha}^{ne})+\frac{\delta
t^{N}}{N!}D_{i}^{N}f_{i,\alpha}^{eq}=-\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}+\delta
t\big{[}G_{i,\alpha}(\mathbf{x},t)+F_{i,\alpha}(\mathbf{x},t)+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)\big{]}+O(\delta t^{N+1}).$ (18b)

Then the equations at first and second orders of $\delta t$ can be obtained,

$D_{i}f_{i,\alpha}^{eq}=-\frac{1}{\delta
t}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}+(G_{i,\alpha}+F_{i,\alpha})+O(\delta
t),$ (19a) $D_{i}(f_{i,\alpha}^{eq}+f_{i,\alpha}^{ne})+\frac{\delta
t}{2}D_{i}^{2}f_{i,\alpha}^{eq}=-\frac{1}{\delta
t}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}+\big{(}G_{i,\alpha}+F_{i,\alpha}+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}\big{)}+O(\delta t^{2}).$ (19b)

From Eq. (19a), we have

$\frac{\delta
t}{2}D_{i}^{2}f_{i,\alpha}^{eq}=-\frac{1}{2}D_{i}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}+\frac{\delta
t}{2}D_{i}(G_{i,\alpha}+F_{i,\alpha})+O(\delta t^{2}).$ (20)

Substituting Eq. (20) into Eq. (19b) yields

$D_{i}f_{i,\alpha}^{eq}+D_{i}\big{(}\delta_{ik}-\frac{1}{2}\mathbf{\Lambda}_{ik}\big{)}f_{k,\alpha}^{ne}+\frac{\delta
t}{2}D_{i}(G_{i,\alpha}+F_{i,\alpha})\\\ =-\frac{1}{\delta
t}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}+G_{i,\alpha}+F_{i,\alpha}+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}+O(\delta t^{2}).$

As pointed out in the previous work [43], to recover Eq. (2) from the MDF-LBM
(3), some appropriate requirements on the collision matrix $\mathbf{\Lambda}$
are needed,

$\displaystyle\sum_{i}\mathbf{e}_{i}\mathbf{\Lambda}_{ik}=s_{0}\mathbf{e}_{k},\
\ \sum_{i}\mathbf{c}_{i}\mathbf{\Lambda}_{ik}=s_{1}\mathbf{c}_{k},$ (21)

where $\mathbf{e}=(1,1,\cdots,1)\in R^{q}$, $s_{0}$ and $s_{1}$ are
eigenvalues of the collision matrix $\mathbf{\Lambda}$ or the relaxation
parameters corresponding to the zero- and first-order moments of the
distribution function. Here it should be noted that a more general case of Eq.
(21) shown in Ref. [43] can also be considered.

According to Eqs. (4), (5) and (6), one can also determine the moments of
$f_{i,\alpha}^{eq}(\mathbf{x},t)$, $G_{i,\alpha,}(\mathbf{x},t)$ and
$F_{i,\alpha}(\mathbf{x},t)$,

$\sum_{i}f_{i,\alpha}^{eq}=\bar{u}_{\alpha},\ \
\sum_{i}\mathbf{c}_{i}f_{i,\alpha}^{eq}=\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha},\
\
\sum_{i}\mathbf{c}_{i}\mathbf{c}_{i}f_{i,\alpha}^{eq}=\bar{u}_{\alpha}c_{s}^{2}\mathbf{I},$
(22a) $\sum_{i}G_{i,\alpha}=0,\ \
\sum_{i}\mathbf{c}_{i}G_{i,\alpha}=\big{(}1-\frac{s_{1}}{2}\big{)}\partial_{t}(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha}),$
(22b) $\sum_{i}F_{i,\alpha}=\bar{F}_{\alpha},\ \
\sum_{i}\mathbf{c}_{i}F_{i,\alpha}=\mathbf{0}.$ (22c)

From Eqs. (16) and (22a), we can first obtain

$\sum_{i}f_{i,\alpha}^{ne}=0.$ (23)

Then with the help of Eq. (21), one can also derive the following equations,

$\sum_{i}\mathbf{e}_{i}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}=s_{0}\sum_{k}\mathbf{e}_{k}f_{k,\alpha}^{ne}=\mathbf{0},$
(24)
$\sum_{i}\mathbf{c}_{i}\mathbf{\Lambda}_{ik}f_{k,\alpha}^{ne}=s_{1}\sum_{k}\mathbf{c}_{k}f_{k,\alpha}^{ne}.$
(25)

###  3.2 Derivation of the convection-diffusion-system based Navier-Stokes-
equations

We now present some details on how to derive the convection-diffusion-system
based NSEs (2) from MDF-LBM (3). To this end, we conduct a summation of Eq.
(3.1), and obtain the following equation,

$\frac{\partial\bar{u}_{\alpha}}{\partial
t}+\nabla\cdot(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})+\nabla\cdot\big{[}\big{(}1-\frac{s_{1}}{2}\big{)}\sum_{i}\mathbf{c}_{i}f_{i,\alpha}^{ne}\big{]}+\frac{\delta
t}{2}\nabla\cdot\big{[}\big{(}1-\frac{s_{1}}{2}\big{)}\partial_{t}(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})\big{]}\\\
=\bar{F}_{\alpha}+O(\delta t^{2}).$

Now let us give an evaluation to the term
$\sum_{i}\mathbf{c}_{i}f_{i,\alpha}^{ne}$. Actually, from Eq. (19a) we have

$\displaystyle\sum_{i}\mathbf{c}_{i}f_{i,\alpha}^{ne}$ $\displaystyle=$
$\displaystyle-\delta
t\sum_{i}\mathbf{c}_{i}\mathbf{\Lambda}_{ik}^{-1}\big{(}D_{k}f_{k,\alpha}^{eq}-G_{k,\alpha}-F_{k,\alpha}\big{)}+O(\delta
t^{2})$ (26) $\displaystyle=$ $\displaystyle-\frac{\delta
t}{s_{1}}\sum_{i}\mathbf{c}_{i}\big{(}D_{i}f_{i,\alpha}^{eq}-G_{i,\alpha}-F_{i,\alpha}\big{)}+O(\delta
t^{2})$ $\displaystyle=$ $\displaystyle-\frac{\delta
t}{2}\partial_{t}(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})-\frac{\delta
t}{s_{1}}c_{s}^{2}\nabla\bar{u}_{\alpha}+O(\delta t^{2}),$

where Eqs. (21) and (22) have been used. Substituting Eq. (26) into Eq. (3.2)
yields

$\frac{\partial\bar{u}_{\alpha}}{\partial
t}+\nabla\cdot(\bar{u}_{\alpha}\mathbf{u}+P\mathbf{E}_{\alpha})=\nabla\cdot\big{[}\big{(}\frac{1}{s_{1}}-\frac{1}{2}\big{)}c_{s}^{2}\delta
t\nabla\bar{u}_{\alpha}\big{]}+\bar{F}_{\alpha}+O(\delta t^{2}).$ (27)

If we neglect the truncation error $O(\delta t^{2})$, one can obtain the
macroscopic convection-diffusion-system based NSEs (2) with the following
viscosity,

$\nu=\big{(}\frac{1}{s_{1}}-\frac{1}{2}\big{)}c_{s}^{2}\delta t.$ (28)

Here we would also like to give a special discussion on how to calculate the
pressure $P$. From Eq. (26) we can get

$\displaystyle\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}^{ne}$ $\displaystyle=$
$\displaystyle\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}-\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}^{eq}$
(29) $\displaystyle=$
$\displaystyle\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}-(\bar{u}_{\alpha}u_{\alpha}+P)$
$\displaystyle=$ $\displaystyle-\frac{\delta
t}{s_{1}}c_{s}^{2}\nabla_{\alpha}\bar{u}_{\alpha}+O(\delta t^{2}+\delta
t\mbox{Ma}^{2}),\ (\alpha\neq 0),$

where Ma is the Mach number. Neglecting the truncation error term $O(\delta
t^{2}+\delta t\mbox{Ma}^{2})$ and summing Eq. (29) over $\alpha$, one can
derive

$\displaystyle\sum_{\alpha=1}^{d}\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}-(|\mathbf{u}|^{2}+dP)=-\frac{\delta
t}{s_{1}}c_{s}^{2}R,$ (30)

where Eq. (1a) has been used. From Eq. (30) we can give an expression to
compute the pressure,

$P=\frac{1}{d}\big{(}\sum_{\alpha=1}^{d}\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}+\frac{\delta
t}{s_{1}}c_{s}^{2}R-|\mathbf{u}|^{2}\big{)}.$ (31)

From above procedure, it can be found that the continuity equation (1a) has
been used to derive Eq. (31), this means that the formula for the pressure
[Eq. (31)] is consistent with the continuity equation. In other words, we do
not need to consider the continuity equation (1a) in the present MDF-LBM, and
the evolution equation (3) with $\alpha=0$ can be omitted. We also note that
the term related to $G_{i,\alpha}$ in Eq. (3) can also be neglected since it
is order of $O(\delta t\mbox{Ma}^{2})$. In this case, the evolution equation
(3) can be simplified by

$f_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta t,t+\delta
t)=f_{i,\alpha}(\mathbf{x},t)-\mathbf{\Lambda}_{ik}\big{[}f_{k,\alpha}(\mathbf{x},t)-f_{k,\alpha}^{eq}(\mathbf{x},t)\big{]}+\delta
t\big{[}F_{i,\alpha}(\mathbf{x},t)+\frac{\delta
t}{2}\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)\big{]}.$ (32)

Then one can derive two special schemes [37] from Eq. (32).
Scheme I: $\gamma=0$. Under this condition,
$\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)=\partial_{t}F_{i,\alpha}(\mathbf{x},t)$.
Considering the Euler method for the time derivative, we have

$f_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta t,t+\delta
t)=f_{i,\alpha}(\mathbf{x},t)-\mathbf{\Lambda}_{ik}\big{[}f_{k,\alpha}(\mathbf{x},t)-f_{k,\alpha}^{eq}(\mathbf{x},t)\big{]}+\frac{\delta
t}{2}\big{[}3F_{i,\alpha}(\mathbf{x},t)-F_{i,\alpha}(\mathbf{x},t-\delta
t)\big{]}.$ (33)

Scheme II: $\gamma=1$. With this choice, we have
$\bar{D}_{i}F_{i,\alpha}(\mathbf{x},t)=D_{i}F_{i,\alpha}(\mathbf{x},t)$, which
can be discretized by the following implicit finite-difference scheme,

$D_{i}F_{i,\alpha}(\mathbf{x},t)=\frac{F_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta
t,t+\delta t)-F_{i,\alpha}(\mathbf{x},t)}{\delta t}.$ (34)

Substituting above equation into Eq. (32) and introducing a new variable
$\bar{f}_{i,\alpha}(\mathbf{x},t)=f_{i,\alpha}(\mathbf{x},t)-\delta
tF_{i,\alpha}(\mathbf{x},t)/2$, one can obtain

$\bar{f}_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta t,t+\delta
t)=\bar{f}_{i,\alpha}(\mathbf{x},t)-\mathbf{\Lambda}_{ik}\big{[}\bar{f}_{k,\alpha}(\mathbf{x},t)-f_{k,\alpha}^{eq}(\mathbf{x},t)\big{]}+\delta
t\big{(}\delta_{ik}-\frac{1}{2}\mathbf{\Lambda}_{ik}\big{)}F_{k,\alpha}(\mathbf{x},t).$
(35)

In this scheme, the macroscopic variable $\bar{u}_{\alpha}$ is computed by

$\bar{u}_{\alpha}=\sum_{i}f_{i,\alpha}=\sum_{i}\bar{f}_{i,\alpha}+\frac{\delta
t}{2}\bar{F}_{\alpha},$ (36)

while the pressure is still given by Eq. (31) in which $f_{i,\alpha}$ is
replaced by $\bar{f}_{i,\alpha}$.

###  3.3 The computational schemes for the velocity gradient, velocity
divergence, strain rate tensor, shear stress and vorticity

Besides the fluid velocity and pressure mentioned above, usually we also need
to consider some other physical variables. For instance, the strain rate
tensor (or the symmetric velocity gradient tensor), shear stress and vorticity
(it is related to the antisymmetric velocity gradient tensor) are also
important in the study of the non-Newtonian fluid flows and turbulence [48,
49], and in the framework of LBM, they also received increasing attention in
the past years [23, 24, 25, 26, 27, 50]. In the SRT-LB model, Artoli et al.
[23] first developed a local scheme for the shear stress where the second-
order moment of non-equilibrium distribution function is adopted, and applied
the scheme to study the blood flows in a symmetric bifurcation. Then Krüger et
al. [24] conducted a theoretical analysis, and found that the the local scheme
for the shear stress has a second-order convergence rate. Chai and Zhao [25]
further considered the forcing term effect on the computation of the strain
rate tensor and shear stress, and proposed two more general local schemes for
the strain rate tensor and shear stress. They also performed the theoretical
analysis and numerical simulations, and demonstrated that the local schemes
have a second-order accuracy in space. Yong and Luo [26] carried out an
asymptotic analysis on the SRT-LB model coupled with the simple bounce-back
boundary condition, and illustrated that the shear stress computed with the
second-order moment of the non-equilibrium distribution function has a second-
order accuracy in space. In addition, the velocity divergence can also be
calculated locally by the second-order moment of non-equilibrium distribution
function [25].

In the LB method, it seems difficult to construct the local schemes for the
velocity gradient and vorticity, and for this reason, some non-local finite-
difference schemes are usually adopted (e.g., Ref. [36]). Recently, through
the careful design of the high-order moments of equilibrium distribution
function, Peng et al. [27] developed a local scheme for the vorticity in the
MRT-LB model with D3Q27 lattice structure. However, this scheme cannot be
extended to other lattice models, e.g., the commonly used D2Q9 lattice model
in two-dimensional space. To overcome this limitation, Hajabdollahi and
Premnath [50] proposed another scheme for the computation of vorticity in the
double-distrbution-function LBM. However, besides the LB model for NSEs, they
also need to introduce another LB model for an additional CDE, which would
bring more computational cost. In this work, we will develop a local scheme
for the velocity gradient, which can be further extended to calculate the
velocity divergence, the strain rate tensor, shear stress and vorcitity
without introducing any additional requirements. Compared to the previous
works [23, 24, 25, 26, 27, 50] where second-order moments of the non-
equilibrium distribution function are adopted, the present local schemes for
velocity gradient, velocity divergence, the strain rate tensor, shear stress
and vorticity only include the first-order moments of the non-equilibrium
distribution function, and what is more, they are not restricted to the
special lattice models.

Now let focus on how to calculate the the velocity gradient, velocity
divergence, strain rate tensor shear stress and vorticity in the MDF-LBM.
Actually, from Eq. (26) we can first obtain the velocity gradient,

$\displaystyle\nabla_{\beta}u_{\alpha}=-\frac{s_{1}}{c_{s}^{2}\delta
t}\sum_{i}\mathbf{c}_{i,\beta}f_{i,\alpha}^{ne},\ (\alpha\neq 0),$ (37)

where the truncation error in Eq. (29) has been neglected. We note that Eq.
(37) is similar to the schemes reported in the previous works [51, 41], and
can be used to derive the schemes for the velocity divergence
$\nabla\cdot\mathbf{u}$, strain rate tensor $\mathbf{S}$, shear stress
$\sigma$ and the antisymmetric velocity gradient tensor $\mathbf{\Omega}$,

$\displaystyle\nabla\cdot\mathbf{u}=\sum_{\alpha=1}^{d}\nabla_{\alpha}u_{\alpha}=-\frac{s_{1}}{c_{s}^{2}\delta
t}\sum_{\alpha=1}^{d}\sum_{i}\mathbf{c}_{i,\alpha}f_{i,\alpha}^{ne},$ (38)
$\displaystyle
S_{\alpha\beta}=\frac{1}{2}(\nabla_{\beta}u_{\alpha}+\nabla_{\alpha}u_{\beta})=-\frac{s_{1}}{2c_{s}^{2}\delta
t}\sum_{i}(\mathbf{c}_{i,\beta}f_{i,\alpha}^{ne}+\mathbf{c}_{i,\alpha}f_{i,\beta}^{ne}),\
(\alpha\neq 0,\beta\neq 0),$ (39)
$\displaystyle\sigma_{\alpha\beta}=2\rho_{0}\nu
S_{\alpha\beta}=-\rho_{0}\big{(}1-\frac{s_{1}}{2}\big{)}\sum_{i}(\mathbf{c}_{i,\beta}f_{i,\alpha}^{ne}+\mathbf{c}_{i,\alpha}f_{i,\beta}^{ne}),\
(\alpha\neq 0,\beta\neq 0),$ (40)
$\displaystyle\Omega_{\alpha\beta}=\frac{1}{2}(\nabla_{\beta}u_{\alpha}-\nabla_{\alpha}u_{\beta})=-\frac{s_{1}}{2c_{s}^{2}\delta
t}\sum_{i}(\mathbf{c}_{i,\beta}f_{i,\alpha}^{ne}-\mathbf{c}_{i,\alpha}f_{i,\beta}^{ne}),\
(\alpha\neq 0,\beta\neq 0).$ (41)

Then based on the following relation between the antisymmetric velocity
gradient tensor $\mathbf{\Omega}$ and vorticity
$\mathbf{\omega}$=$\nabla\times\mathbf{u}$,

$\Omega_{\alpha\beta}=-\frac{1}{2}\epsilon_{\alpha\beta\gamma}\omega_{\gamma},$
(42)

one can determine the component of vorticity $\omega_{\gamma}$,
$\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor.

Finally, we also give some remarks on the present MDF-LBM for the convection-
diffusion-system based incompressible NSEs.
Remark I: In the above analysis, if the diffusive scaling ($\delta
t\propto\delta x^{2}$) is considered [42, 52], the evolution equation (3) can
be simply written as

$f_{i,\alpha}(\mathbf{x}+\mathbf{c}_{i}\delta t,t+\delta
t)=f_{i,\alpha}(\mathbf{x},t)-\mathbf{\Lambda}_{ik}\big{[}f_{k,\alpha}(\mathbf{x},t)-f_{k,\alpha}^{eq}(\mathbf{x},t)\big{]}+\delta
tF_{i,\alpha}(\mathbf{x},t).$ (43)

Then we can also derive the macroscopic equation (27) from Eq. (43), while the
truncation error is $O(\delta x^{2})$ rather than $O(\delta t^{2})$. On the
other hand, for the specified kinematic viscosity $\nu$ and the relaxation
parameter $s_{1}$, one can also obtain the diffusive scaling $\delta
t\propto\delta x^{2}$ from Eq. (28), which means that the diffusive scaling
used in the LBM is reasonable. Additionally, the computational schemes for the
velocity, pressure, velocity gradient, velocity divergence, strain rate
tensor, shear stress and vorticity are the same as Eqs. (16), (31), (37),
(38), (39), (40) and (41). In the following, for simplicity we would consider
the MDF-LBM (43) for the convection-diffusion-system based incompressible
NSEs.
Remark II: In the single LBM for the NSEs [23, 24, 25, 26], the velocity
divergence, strain rate tensor and shear stress can be computed locally with
the second-order moments of the non-equilibrium distribution function, while
in the present MDF-LBM for the convection-diffusion-system based
incompressible NSEs, only the first-order moments of the non-equilibrium
distribution function are needed, as seen from Eqs. (38), (39) and (40).
Remark III: In the commonly used LBM for NSEs, it is difficult to compute the
velocity gradient locally. However, in the present MDF-LBM for the convection-
diffusion-system based incompressible NSEs, the velocity gradient can be
calculated locally from Eq. (37), which can also be used to determine some
other physical variables, including the velocity divergence
$\nabla\cdot\mathbf{u}$, strain rate tensor $\mathbf{S}$, shear stress
$\sigma$ and vorticity $\omega$. Additionally, we also note that although some
local schemes have been developed for vorticity, there are still some
limitations in these available works [27, 50]. For example, the local scheme
proposed by Peng et al. [27] is only suitable for the D3Q27 lattice model, and
cannot be extended to some other lattice models. Compared to the scheme in
Ref. [27], the scheme developed by Hajabdollahi and Premnath [50] is more
general, while an additional CDE must be introduced and solved by another LB
model, which would make the computational cost more expensive.
Remark IV: In the single LBM for the NSEs [22], the pressure is related to
density through the relation $P=\rho c_{s}^{2}$ where the density is
calculated with the zeroth-order of distribution function, while in the LBM
for incompressible NSEs [53], the pressure is determined by the zeroth-order
moment of distribution function without considering the one at zeroth
direction. In the present MDF-LBM, the pressure is computed by the first-order
moment of distribution function, but the one with $\alpha=0$ is not included.
In addition, it should be noted that the formula for the pressure [see Eq.
(31)] is consistent with the continuity equation (1a), and thus we do not need
to consider the evolution equation (3) with $\alpha=0$.

##  4 Numerical results and discussion

For simplicity but without loss of generality, we only considered two-
dimensional problems in this section, and adopted some two-dimensional
benchmark problems, including the Poiseuille flow, the simplified four-roll
mill problem and lid-driven cavity flow, to test the developed MDF-LBM with
D2Q5 lattice model. In this lattice model, the collision matrix
$\mathbf{\Lambda}$ is taken as

$\mathbf{\Lambda}=\mathbf{M}^{-1}\mathbf{S}_{d}\mathbf{M}.$ (44)

The transformation matrix $\mathbf{M}$ and relaxation matrix $\mathbf{S}_{d}$
appeared in above equation are defined as [21]

$\mathbf{M}=\left({\begin{array}[]{*{20}{c}}1&1&1&1&1\\\ 0&c&0&{-c}&0\\\
0&0&c&0&{-c}\\\ 0&c^{2}&{-c^{2}}&c^{2}&{-c^{2}}\\\
{-4c^{2}}&c^{2}&c^{2}&c^{2}&c^{2}\\\ \end{array}}\right),$ (45a)
$\mathbf{S}_{d}=\mathbf{diag}(s_{0},s_{1},s_{1},s_{2},s_{2}),$ (45b)

where $s_{i}$ is the relaxation parameter corresponding to the $i$th-order
moment of the distribution function. Additionally, according to Eq. (7), the
relations among the weigh coefficients and lattice speed in the D2Q5 lattice
model can be generally expressed as

$\omega_{i}=\frac{1-\omega_{0}}{4},\ 2\omega_{i}c^{2}=c_{s}^{2},\ i=1-4,$ (46)

where $\omega_{0}$ is considered as a free parameter within the range $(0,1)$.
Actually, if $\omega_{0}=1/5$, one can obtain $\omega_{i}=1/5\ (i=1-4)$ and
$c_{s}^{2}=2c^{2}/5$. However, if $\omega_{0}=1/3$, we have $\omega_{i}=1/6\
(i=1-4)$ and $c_{s}^{2}=c^{2}/3$, which is the same as Eq. (11b) and would be
used in the following simulations. Unless otherwise stated, the initialization
of the distribution functions is realized by their equilibrium distribution
functions where the pressure and velocity are given by $P=1$ and
$u_{1}=u_{2}=0$. The anti-halfway bounce-back scheme for the MDF-LBM [17] is
used to treat Dirichlet boundary conditions of velocity.

In our simulations, to test the accuracies of the MDF-LBM and the local
schemes for the physical variables, the following $L^{2}$ norm of relative
error ($E_{L^{2}}$) is adopted,

$E_{L^{2}}(\phi)=\frac{\|\phi_{a}(x,y)-\phi_{n}(x,y)\|_{2}}{\|\phi_{a}(x,y)\|_{2}},$
(47)

where $\phi$ denotes velocity, velocity gradient, velocity divergence, strain
rate tensor or vorticity, the subscripts $a$ and $n$ represent the analytical
and numerical solutions of $\phi$.

###  4.1 The two-dimensional Poiseuille flow

The first problem we considered is the simple Poiseuille flow, which is driven
by an external constant force in $x$ direction ($F_{1}=1.0\times 10^{-6}$), as
shown in Fig. 1. For this problem, one can obtain its analytical solutions of
velocity $\mathbf{u}=(u_{1},u_{2})^{\top}$, velocity gradient
$\nabla\mathbf{u}$, velocity divergence $\nabla\cdot\mathbf{u}$, strain rate
tensor $\mathbf{S}$ and vorticity $\omega$,

$u_{1}=\frac{F_{1}}{2\nu}H^{2}\big{[}\frac{y}{H}-\big{(}\frac{y}{H}\big{)}^{2}\big{]},\
u_{2}=0,$ (48a) $\frac{\partial u_{1}}{\partial x}=\frac{\partial
u_{2}}{\partial x}=\frac{\partial u_{2}}{\partial y}=0,\ \frac{\partial
u_{1}}{\partial y}=\frac{F_{1}}{2\nu}H\big{(}1-2\frac{y}{H}\big{)},$ (48b)
$\nabla\cdot\mathbf{u}=\frac{\partial u_{1}}{\partial x}+\frac{\partial
u_{2}}{\partial y}=0,$ (48c) $S_{xx}=S_{yy}=0,\
S_{xy}=S_{yx}=\frac{F_{1}}{4\nu}H\big{(}1-2\frac{y}{H}\big{)},$ (48d)
$\omega=\frac{\partial u_{2}}{\partial x}-\frac{\partial u_{1}}{\partial
y}=-\frac{F_{1}}{2\nu}H\big{(}1-2\frac{y}{H}\big{)}.$ (48e)

Figure 1: The schematic of the Poiseuille flow.

We performed some simulations with the present MDF-LBM, and the lattice size
$32\times 32$ is adopted for the computational domain $[0,1]\times[0,1]$ with
the periodic boundary condition in $x$ direction. We presented some numerical
results in Figs. 2-6 where the kinematic viscosity is $\nu=0.001$, the
relaxation parameters are set to be $s_{0}=1$, $s_{1}=1.2$ and
$s_{2}=8(2-s_{1})/(8-s_{1})$ [21, 54]. As shown in these figures, the
numerical results of the velocity $\mathbf{u}$, the components of velocity
gradient ($\partial{u_{1}}/\partial x$ and $\partial{u_{1}}/\partial y$),
velocity divergence $\nabla\cdot\mathbf{u}$, the components of strain rate
tensor ($S_{xx}$ and $S_{xy}$), and the vorticity $\omega$ are in good
agreement with analytical solutions. In addition, we also measured the
relative errors of velocity $u_{1}$, the component of velocity gradient
$\partial{u_{1}}/\partial y$, the component of strain rate tensor $S_{xy}$ and
the vorticity $\omega$, and found that their values are about
$E_{L^{2}}(u_{1})=5.9427\times 10^{-4}$, $E_{L^{2}}(\partial{u_{1}}/\partial
y)=9.7541\times 10^{-17}$, $E_{L^{2}}(S_{xy})=9.7541\times 10^{-17}$ and
$E_{L^{2}}(\omega)=9.7541\times 10^{-17}$. It should be noted that the
component of velocity gradient $\partial{u_{1}}/\partial y$, the component of
strain rate tensor $S_{xy}$ and the vorticity $\omega$ can achieve the machine
precision. This is because the local schemes (37), (39) and (41) are of
second-order accuracy, while the distributions of these physical variables are
only linear. Finally, it is also observed from Fig. 4 that the continuity
equation $\nabla\cdot\mathbf{u}=0$ is preserved automatically.

Figure 2: The numerical and analytical solutions of the velocity [(a):
$u_{1}$, (b): $u_{2}$].

Figure 3: The numerical and analytical solutions of the velocity gradient
[(a): $\partial{u_{1}}/\partial x$, (b): $\partial{u_{1}}/\partial y$]. Figure
4: The numerical and analytical solutions of the velocity divergence.

Figure 5: The numerical and analytical solutions of the strain rate tensor
[(a): $S_{xx}$, (b): $S_{xy}$]. Figure 6: The numerical and analytical
solutions of the vorticity $\omega$.

###  4.2 The simplified two-dimensional four-roll mill problem

The second example we used to test the present MDF-LBM is the simplified two-
dimensional four-roll mil problem. The schematic of the problem is given in
Fig. 7 where the four rollers are replaced by a body force to drive fluid
flow, the physical domain of the problem is $[0,2\pi]\times[0,2\pi]$ with the
periodic boundary conditions in both $x$ and $y$ directions. As a benchmark
problem, it has also been used to test the accuracy of the LBM [25, 55] for
the following two reasons. The first is that the problem has analytical
solutions of velocity, velocity gradient, velocity divergence, strain rate
tensor and vorticity. The second is that the boundary condition of the problem
is periodic such that its effect on the numerical results can be excluded
[25]. For the incompressible flow driven by the following force,

$F_{1}=U_{0}^{2}\sin(x)\cos(x)+2\nu U_{0}\sin(x)\cos(y),$ (49a)
$F_{2}=U_{0}^{2}\sin(y)\cos(y)-2\nu U_{0}\sin(y)\cos(x),$ (49b)

one can obtain the analytical solutions of velocity, velocity gradient,
velocity divergence, stain rate tensor and vorticity,

$u_{1}=U_{0}\sin(x)\cos(y),\ u_{2}=-U_{0}\cos(x)\sin(y),$ (50a)
$\frac{\partial u_{1}}{\partial x}=U_{0}\cos(x)\cos(y)=-\frac{\partial
u_{2}}{\partial y},\ \frac{\partial u_{1}}{\partial
y}=-U_{0}\sin(x)\sin(y)=-\frac{\partial u_{2}}{\partial x},$ (50b)
$\nabla\cdot\mathbf{u}=\frac{\partial u_{1}}{\partial x}+\frac{\partial
u_{2}}{\partial y}=0,$ (50c) $S_{xx}=U_{0}\cos(x)\cos(y)=-S_{yy},\
S_{xy}=S_{yx}=0,$ (50d) $\omega=\frac{\partial u_{2}}{\partial
x}-\frac{\partial u_{1}}{\partial y}=2U_{0}\sin(x)\sin(y).$ (50e)

Figure 7: The schematic of the two-dimensional four-roll mill problem.

We carried out some numerical experiments with a lattice size $64\times 64$,
and the results are shown in Figs. 8-12 where $U_{0}=0.0001$, $\nu=0.01$, the
relaxation parameters are the same as those used in the first problem. As seen
from these figures, the numerical results of velocity, velocity gradient,
velocity divergence, strain rate tensor and vorticity are very close to the
corresponding analytical solutions, and specially, from the contour lines of
vorticity shown in Fig. 12(b), one can clearly observe that there are four
vortices formed at the locations of the rollers.

Figure 8: The numerical and analytical solutions of velocity at different
positions [(a): $u_{1}$, (b): $u_{2}$; solid line: analytical solution,
symbol: numerical solution].

Figure 9: The numerical and analytical solutions of velocity gradient at
different positions [(a): $\partial u_{1}/\partial x$, (b): $\partial
u_{1}/\partial y$; solid line: analytical solution, symbol: numerical
solution]. Figure 10: The numerical and analytical solutions of the velocity
divergence at different positions (solid line: analytical solution, symbol:
numerical solution).

Figure 11: The numerical and analytical solutions of the components of strain
rate tensor at different positions [(a): $S_{xx}$, (b) $S_{yy}$; solid line:
analytical solution, symbol: numerical solution].

Figure 12: The numerical and analytical solutions of vorticity at different
positions (a) (solid line: analytical solution, symbol: numerical solution),
and the contour lines of vorticity (b) (black line: analytical solution,
colored line: numerical solution).

This problem is also used to test the convergence rates of the present MDF-LBM
for velocity and the local schemes for the strain rate tensor and vorticity.
To this end, we calculated the relative errors of velocity, the strain rate
tensor and the vorticity at different lattice sizes ($N\times N=16\times 16$,
$32\times 32$, $48\times 48$ and $64\times 64$), and plotted them in Fig. 13.
From this figure, one can find that the present MDF-LBM and the local schemes
are of second-order accuracy in space.

In addition, it is well known that the eigenvalue (or relaxation parameter)
$s_{1}$ of the collision matrix $\mathbf{\Lambda}$ is a key parameter in the
LBM, and may also affect the numerical results. Here we also conducted some
simulations to test the effect of the relaxation parameter $s_{1}$, and
presented the relative errors of the velocity, the strain rate tensor and the
vorticity at three different values of the relaxation parameter $s_{1}$ in
Table 1 where the viscosity and lattice size are fixed as $\nu=0.01$ and
$64\times 64$. As shown in this table, the relaxation parameter $s_{1}$ has no
apparent influences on the numerical results, especially on the stain rate
tensor and vorticity.

Figure 13: The convergence rates of the MDF-LBM and the local schemes for the strain rate tensor and vorticity. Table 1: The relative errors of velocity, strain rate tensor and vorticity under different relaxation parameter $s_{1}$. | $s_{1}=0.7$ | $s_{1}=1.2$ | $s_{1}=1.7$
---|---|---|---
$E_{L^{2}}(u_{1})$ | $3.9869\times 10^{-3}$ | $2.0745\times 10^{-3}$ | $1.2871\times 10^{-3}$
$E_{L^{2}}(u_{2})$ | $3.9869\times 10^{-3}$ | $2.0745\times 10^{-3}$ | $1.2871\times 10^{-3}$
$E_{L^{2}}(S_{xx})$ | $8.0330\times 10^{-4}$ | $8.0328\times 10^{-4}$ | $8.0327\times 10^{-4}$
$E_{L^{2}}(S_{yy})$ | $8.0330\times 10^{-4}$ | $8.0328\times 10^{-4}$ | $8.0327\times 10^{-4}$
$E_{L^{2}}(\omega)$ | $8.0333\times 10^{-4}$ | $8.0330\times 10^{-4}$ | $8.0330\times 10^{-4}$

###  4.3 The two-dimensional lid-driven cavity flow

The last problem we considered is the two-dimensional lid-driven cavity flow,
and the schematic of the problem is shown in Fig. 8 where the length of the
square cavity is $L=1$. Compared to the previous problems, it is more
complicated since there is no exact solution available. Although the geometry
of the problem is very simple, the lid-driven cavity flow is of great
scientific interest because it displays rich fluid mechanical phenomena,
especially the complex vortex dynamics [56]. The flow in the square cavity is
driven by the top moving wall with a constant velocity
$(U_{1},U_{2})^{\top}=(1,0)^{\top}$, and a primary vortex in the center and
some secondary vortices at the corners would be formed with the increase of
the Reynolds number (Re$=LU_{1}/\nu$). Actually, the lid-driven cavity flow,
as a classic benchmark problem, has also been widely used to test the capacity
of the numerical methods [57, 58, 59, 60, 36].

In this part, we would conduct some numerical simulations of lid-driven cavity
flows at different Reynolds numbers, and to ensure the incompressible
condition to be valid, the discrete velocity $c=10$ is adopted to give a small
Mach number (Ma$=U_{1}/c_{s}$). In our simulations, to obtain accurate
results, the lattice size is set to be $512\times 512$. We first presented the
velocity profiles along the vertical and horizontal lines through geometric
center of the cavity in Fig. 9 where Re$=100,400$ and $1000$. As seen from
this figure, the numerical results are in good agreement with the available
data [57]. In addition, we also plotted the distributions of the pressure in
Fig. 10, and found that these results qualitatively agree with the previous
works [58, 36]. To show the complex dynamics of lid-driven cavity flows, we
also plotted the streamlines and vorticity contours in Figs. 11 and 12. From
these two figures, one can observe that a primary vortex in the center of the
cavity and two secondary vortices at two bottom corners are formed.
Simultaneously, it is also found that with the increase of the Reynolds
number, the primary vortex would move towards the center of the cavity. To
quantify these results, we measured the vorticities and locations of the
primary and secondary vortices, and listed them in Table 2. As shown in this
table, the present results are very close to those reported in some previous
studies [57, 58, 36].

Figure 14: The schematic of the two-dimensional lid-driven cavity flow.

Figure 15: The profiles of velocity along vertical and horizontal lines
through geometric center of the cavity [(a): $u_{1}$ along vertical line; (b):
$u_{2}$ along horizontal line].

Figure 16: The contours of pressure at Re=100 (a), 400 (b) and 1000 (c).

Figure 17: The streamlines of lid-driven cavity flows at Re=100 (a), 400 (b)
and 1000 (c).

Figure 18: The contours of the vorticity at Re=100 (a), 400 (b) and 1000 (c). Table 2: The vorticity ($\omega$) and location ($x$, $y$) of the vortex. Vortex | Reference |  | Re=100 | Re=400 | Re=1000
---|---|---|---|---|---
Primary vortex | Present | $\omega$ | 3.1655 | 2.2797 | 1.9794
$x$ | 0.6167 | 0.5540 | 0.5303
$y$ | 0.7373 | 0.6053 | 0.5651
Ref. [57] | $\omega$ | 3.1665 | 2.2947 | 2.0497
$x$ | 0.6172 | 0.5547 | 0.5313
$y$ | 0.7344 | 0.6055 | 0.5626
Ref. [58] | $\omega$ | 3.1348 | 2.2910 | 2.0760
$x$ | 0.6196 | 0.5608 | 0.5333
$y$ | 0.7373 | 0.6078 | 0.5647
| Ref. [36] | $\omega$ | 3.1629 | 2.2950 | 2.0678
| $x$ | 0.6150 | 0.5546 | 0.5312
| $y$ | 0.7378 | 0.6053 | 0.5663
Bottom-left vortex | Present | $\omega$ | -0.0149 | -0.0583 | -0.3514
$x$ | 0.0343 | 0.0511 | 0.0833
$y$ | 0.0345 | 0.0472 | 0.0778
Ref. [57] | $\omega$ | -0.0155 | -0.0570 | -0.3618
$x$ | 0.0313 | 0.0508 | 0.0859
$y$ | 0.0391 | 0.0469 | 0.0781
Ref. [58] | $\omega$ | - | - | -
$x$ | 0.0392 | 0.0549 | 0.0902
$y$ | 0.0353 | 0.0510 | 0.0784
| Ref. [36] | $\omega$ | 0.0145 | 0.0602 | 0.3557
| $x$ | 0.0341 | 0.0517 | 0.0828
| $y$ | 0.0341 | 0.0478 | 0.0789
Bottom-right vortex | Present | $\omega$ | -0.0352 | -0.4486 | -1.0690
$x$ | 0.9425 | 0.8855 | 0.8650
$y$ | 0.0619 | 0.1222 | 0.1120
Ref. [57] | $\omega$ | -0.0331 | -0.4335 | -1.1547
$x$ | 0.9453 | 0.8906 | 0.8594
$y$ | 0.0625 | 0.1250 | 0.1094
Ref. [58] | $\omega$ | - | - | -
$x$ | 0.9451 | 0.8902 | 0.8667
$y$ | 0.0627 | 0.1255 | 0.1137
| Ref. [36] | $\omega$ | -0.0348 | -0.4451 | -1.1039
| $x$ | 0.9425 | 0.8560 | 0.8645
| $y$ | 0.0614 | 0.1218 | 0.1121

##  5 Conclusions

In this work, a MDF-LBM coupled with the MRT model is developed for
incompressible NSEs. To do this, the NSEs are first reformulated into a
convection-diffusion system, and then the MDF-LBM is proposed for the
convection-diffusion system. What is more, besides the macroscopic pressure
and velocity, we also proposed some local schemes for the velocity gradient,
velocity divergence, strain rate tensor, shear stress and vorticity in the
framework of MDF-LBM. To test the capacity of the MDF-LBM and local schemes,
three benchmark problems, including the two-dimensional Poiseuille flow, the
simplified four-roll mill problem and the lid-driven cavity flow, are
considered. The numerical results show that the present MDF-LBM and local
schemes are efficient, and also have a second-order convergence rate in space.

It should be noted that compared to the classic LBM for two-dimensional
incompressible NSEs where the D2Q9 lattice model should be used, the MDF-LBM
is more flexible since the D2Q4, D2Q5 or D2Q9 lattice model can be adopted to
give correct incompressible NSEs. Additionally, in the MDF-LBM for convection-
diffusion-system based NSEs, some physical variables (e.g., the velocity
gradient, velocity divergence, the strain rate tensor, the shear stress and
the vorticity) can be computed locally through the first-order moments of the
non-equilibrium distribution function, while in the commonly used single LBM
for incompressible NSEs, usually only the velocity divergence, strain rate
tensor and the shear stress can be determined locally by the the second-order
moments of the non-equilibrium distribution function.

Finally, we would also like to point out that the present MDF-LBM can be also
extended to study the thermal flows and multiphase fluid system governed by
the incompressible NSEs and CDE, which would be considered in another work.

## Acknowledgements

This work was financially supported by the National Natural Science Foundation
of China (Grants No. 12072127 and No. 51836003).

## References

  * [1] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: Theory and applications, Phys. Reports 222 (1992) 145-197.
  * [2] S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329-364.
  * [3] C. K. Aidun J. R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439-472.
  * [4] A. Xu, W. Shyy, T. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. 33 (2017) 555-574.
  * [5] H. Wang, X. Yuan, H. Liang, Z. Chai, B. Shi, A brief review of the phase-field-based lattice Boltzmann method for multiphase flows, Capillary 2 (2019) 33-52.
  * [6] D. A. Wolf-Gladrow, Lattice-gas cellular automata and lattice Boltzmann models: An introduction, Springer, Berlin, 2000.
  * [7] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, 2001.
  * [8] Z. Guo, C. Shu, Lattice Boltzmann Method and Its Applications in Engineering, World Scientific Publishing Co. Pte. Ltd., Singapore, 2013.
  * [9] H. Huang, M. C. Sukop, X.-Y. Lu, Multiphase lattice Boltzmann methods: Theory and applications, John Wiley & Sons, UK, 2016.
  * [10] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, Switzerland, 2017.
  * [11] F. Bouchut, Construction of BGK models with a family of kinetic entropies for a given system of conservation laws, J. Stat. Phys. 95 (1999) 113-170.
  * [12] M. F. Carfora, R. Natalini, A discrete kinetic approximation for the incompressibe Navier-Stokes equations, ESAIM: M2AN 42 (2008) 93-112.
  * [13] R. Bianchini, R. Natalini, Covergence of a vector-BGK approximation for the incompressible Navier-Stokes equations, Kinet. Relat. Mod. 12 (2019) 133-158.
  * [14] F. Dubois, Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes, Int. J. Modern Phys. C 25 (2014) 1441014.
  * [15] J. Zhao, Discrete-velocity vector-BGK models based numerical methods for the incompressible Navier-Stokes equations, Commun. Comput. Phys. 29 (2021) 420-444.
  * [16] J. Zhao, Z. Zhang, Seven-velocity three-dimensional vectorial lattice Boltzmann method including various types of approximations to the pressure and two-parameterized second-order boundary treatments, Comput. Math. Appl. 80 (2020) 2764-2779.
  * [17] J. Zhao, Z. Zhang, W.-A. Yong, Vector-type boundary schemes for the lattice Boltzmann method based on vector-BGK models, SIAM J. Sci. Comput. 42 (2020) B1250-B1270.
  * [18] X. Shan, Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method, Phys. Rev. E 55 (1997) 2780-2788.
  * [19] X. He, S. Chen, G. D. Doolen, A novel thermal model for the lattice Boltzmann method in incompressible limit, J. Comput. Phys. 146 (1998) 282-300.
  * [20] Z. Guo, B. Shi, C. Zheng, A coupled lattice BGK model for the Boussinesq equations, Int. J. Numer. Meth. Fluids 39 (2002) 325-342
  * [21] S. Cui, N. Hong, B. Shi, Z. Chai, Discrete effect on the halfway bounce-back boundary condition of multiple-relaxation-time lattice Boltzmann model for convection-diffusion equations, Phys. Rev. E 93 (2016) 043311.
  * [22] Y. H. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479-484.
  * [23] A. M. Artoli, D. Kandhai, H. C. J. Hoefsloot, A. G. Hoekstra, P. M. A. Sloop, Lattice BGK simulations of flow in a symmetric bifurcation Future Gener. Comput. Syst. 20 (2004) 909-916.
  * [24] T. Krüger, F. Varnik, D. Raabe, Shear stress in lattice Boltzmann simulations, Phys. Rev. E 79 (2009) 046704.
  * [25] Z. Chai, T. S. Zhao, Effect of the forcing term in the multiple-relaxation-time lattice Boltzmann equation on the shear stress or the strain rate tensor, Phys. Rev. E 86 (2012) 016705.
  * [26] W.-A. Yong, L.-S. Luo, Accuracy of the viscous stress in the lattice Boltzmann equation with simple boundary conditions, Phys. Rev. E 86 (2012) 065701(R).
  * [27] C. Peng, Z. Guo, L.-P. Wang, Lattice Boltzmann model capable of mesoscopic vorticity computation, Phys. Rev. E 96 (2017) 053304.
  * [28] L. D. Landau, E. M. Lifshitz, Fluid Mechanics, second Edition, Pergamon Press, Oxford, 1987.
  * [29] P. K, Kundu, I. M. Cohen, D. Dowling, Fluid Mechanics, Elsevier, California, 2016.
  * [30] H. Chen, S. Chen, W. H. Matthaeus, Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method, Phys. Rev. A 45 (1992) R5339-R5342.
  * [31] I. Ginzburg, Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation, Adv. Water Resour. 28 (2005) 1171-1195.
  * [32] I. Ginzburg, F. Verhaeghe, D. d’Humières, Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions, Commun. Comput. Phys. 3 (2008) 427-478.
  * [33] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weave (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., Vol. 159, AIAA, Washington, DC, 1992, pp. 450-458.
  * [34] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546-6562.
  * [35] C. Pan, L.-S. Luo, C. T. Miller, An evaluation of lattice Boltzmann schemes for porous medium flow simulation, Comput. Fluids 35 (2006) 898-909.
  * [36] L.-S. Luo, W. Liao, X. Chen, Y. Peng, W. Zhang, Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations, Phys. Rev. E 83 (2011) 056710.
  * [37] Z. Chai, B. Shi, Z. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection Cdiffusion equations, J. Sci. Comput. 69 (2016) 355-390.
  * [38] Z. Chai, C. Huang, B. Shi, Z. Guo, A comparative study on the lattice Boltzmann models for predicting effective diffusivity of porous media, Int. J. Heat Mass Transfer 98 (2016) 687-696.
  * [39] I. Rasin, S. Succi, W. Miller, A multi-relaxation lattice kinetic method for passive scalar diffusion, J. Comput. Phys. 206 (2005) 453-462.
  * [40] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774-7795.
  * [41] Z. Chai, T. S. Zhao, Nonequilibrium scheme for computing the flux of the convection-diffusion equation in the framework of the lattice Boltzmann method, Phys. Rev. E 90 (2014) 013305.
  * [42] M. Zhang, W. Zhao, P. Lin, Lattice Boltzmann method for general convection-diffusion equations: MRT model and boundary schemes, J. Comput. Phys. 289 (2019) 147-163.
  * [43] Z. Chai, B. Shi, Multiple-relaxation-time lattice Boltzmann method for the Navier-Stokes and nonlinear convection-diffusion equations: Modeling, analysis, and elements, Phys. Rev. E 102 (2020) 023306.
  * [44] S. Chapman, T. G. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, Cambridge, 1970.
  * [45] W.-A. Yong, W. Zhao, L.-S. Luo, Theory of the lattice Boltzmann method: Derivation of macroscopic equations via the Maxwell iteration, Phys. Rev. E 93 (2016) 033310.
  * [46] A. Wagner, Thermodynamic consistency of liquid-gas lattice Boltzmann simulations, Phys. Rev. E 74 (2006) 056703.
  * [47] D. d’Humières, I. Ginzburg, Viscosity independent numerical errors for Lattice Boltzmann models: From recurrence equations to ”magic” collision numbers, Comput. Math. Appl. 58 (2009) 823-840.
  * [48] R. B. Bird, R. C. Armstrong, O. Hassager, Dynamics of Polymeric liquids (Volume 1: Fluid Mechanics), John Wiley & Sons, York, USA, 1987.
  * [49] W. D. McComb, The Physics of Fluid Turbulence, Clarendon Press, Oxford, 1990.
  * [50] F. Hajabdollahi, K. N. Premnath, Local vorticity computation approach in double distribution functions based lattice Boltzmann methods for flow and scalar transport, Int. J. Heat Fluid Flow 83 (2020) 108577.
  * [51] Z. Chai, T. S. Zhao, Lattice Boltzmann model for the convection-diffusion equation, Phys. Rev. E 87 (2013) 063309.
  * [52] M. Junk, A. Klar, L.-S. Luo, Asymptotic analysis of the lattice Boltzmann equation, J. Comput. Phys. 210 (2005) 676-704.
  * [53] Z. Guo, B. Shi, N. Wang, Lattice BGK model for incompressible Navier-Stokes equation, J. Comput. Phys. 165 (2000) 288-306.
  * [54] Y. Wu, Y. Zhao, Z. Chai, B. Shi, Discrete effects on some boundary schemes of multiple-relaxation-time lattice Boltzmann model for convection-diffusion equations, Comput. Math. Appl. 80 (2020) 531-551.
  * [55] O. Malaspinas, N. Fiétier, M. Deville, Lattice Boltzmann method for the simulation of viscoelastic fluid flows, J. Non-Newtonian Fluid Mech. 165 (2010) 1637-1653.
  * [56] P. N. Shankar, M. D. Deshpande, Fluid Mechanics in the driven cavity, Annu. Rev. Fluid Mech. 32 (2000) 93-136.
  * [57] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387-411.
  * [58] S. Hou, Q. Zou, S. Chen, G. Doolen, Simulation of cavity flow by the lattice Boltzmann method, J. Comput. Phys. 118 (1995) 329-347.
  * [59] O. Botella, R. Peyret, Benchmark spectral results on the lid-riven cavity flow, Comput. Fluids 27 (1998) 421-433.
  * [60] E. Erturk, C. Gökçöl, Fourth-order compact formulation of Navier-Stokes equations and deriven cativy flow at high Reynolds number, Int. J. Numer. Meth. Fluids 50 (2006) 421-436.