Journal of Computational Physics 165, 288–306 (2000)

doi:10.1006/jcph.2000.6616, available online at http://www.idealibrary.com on

Lattice BGK Model for Incompressible

Navier–Stokes Equation
Zhaoli Guo,∗ Baochang Shi,† and Nengchao Wang†
∗ National Laboratory of Coal Combustion, and Department of Computer Science, Huazhong University of
Science and Technology, Wuhan 430074, People’s Republic of China; and †Department of Mathematics,
Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China
E-mail: sbchust@public.wuhan.cngb.com

Received May 10, 1999; revised February 28, 2000

Most of the existing lattice Boltzmann BGK models (LBGK) can be viewed as
compressible schemes to simulate incompressible fluid flows. The compressible ef-
fect might lead to some undesirable errors in numerical simulations. In this paper a
LBGK model without compressible effect is designed for simulating incompressible
flows. The incompressible Navier–Stokes equations are exactly recovered from this
incompressible LBGK model. Numerical simulations of the plane Poiseuille flow,
the unsteady 2-D shear decaying flow, the driven cavity flow, and the flow around
a circular cylinder are performed. The results agree well with the analytic solutions
and the results of previous studies. °c 2000 Academic Press
Key Words: Lattice BGK method; Incompressible Navier–Stokes equation.

1. INTRODUCTION

The Lattice Boltzmann BGK (LBGK) method is a new numerical scheme for simulating
viscous compressible flows in the subsonic regime [2]. In recent years, LBGK has achieved
great success in simulations of fluid flows and in modeling physics in fluids. Through mul-
tiscaling expansion [7], the compressible Navier–Stokes equations can be recovered from
the lattice Boltzmann BGK equation on the assumptions that (i) the Mach number is small,
and (ii) the density varies slowly. Therefore, theoretically the LBGK model can only be
used to simulate compressible flows in the incompressible limit. When used for incompress-
ible flows, it must be viewed as an artifical compressible method for the incompressible
Navier–Stokes equations. In such circumstances, the LBGK solutions might depart from
the direct solutions of the incompressible Navier–Stokes equations [14]; at least part of the
departures might be attributable to the effects of the compressibility of the LBGK model.
Some efforts have been made to reduce or eliminate such errors [7, 9, 13, 22]. However,
most of these existing incompressible LBGK models can be used only to simulate steady
288
0021-9991/00 $35.00
Copyright ° c 2000 by Academic Press
All rights of reproduction in any form reserved.
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 289

flows. By neglecting terms of higher order Mach number in the equilibrium density dis-
tribution function, He and Luo [9] proposed an incompressible LBGK model in which
the distribution function is of pressure representation. From the model, the incompressible
Navier–Stokes equations in artificial compressible form can be derived. He’s model can be
used for both incompressible steady and unsteady flows. However, before the model is used,
the average pressure of the flow must be specified in advance. In some cases, especially in
practical problems, the average pressure is not known or cannot be prescribed precisely.
Furthermore, when used to simulate unsteady incompressible flows, the model requires an
additional condition, T À L/cs (T and L are characteristic time and length, respectively),
to neglect the artifical compressible effect.
Considering the significance of the incompressible Navier–Stokes equations in theory
and applications, it is necessary to establish a LBGK model which can exactly model the
incompressible Navier–Stokes equations in general. It is well known that the small Mach
number limit is equivalent to the incompressible limit, so it is possible to set up a LBGK
model which can properly model the incompressible Navier–Stokes equations only with the
small Mach number limit. In this paper such an incompressible LBGK model is proposed.
The rest of the paper is organized as follows. In Section 2, the 2-D incompressible LBGK
model is designed. In Section 3, numerical results are performed to test the incompressible
LBGK model. The simulations include the steady Poiseuille flow, the unsteady decaying
shear flow, the driven cavity flow, and flow around a circular cylinder. Section 4 summarizes
the results and concludes the paper.

2. 2-D INCOMPRESSIBLE LBGK MODEL

In this section, we will propose a 9-bit LBGK model in two-dimensional space. The
approach can also be used to develop other incompressible LBGK models in either two- or
three-dimensional space.

2.1. d2q9 LBGK Model

The velocity directions of the d2q9 LBGK model [17] are defined as

(0, 0)
 i =0
ei = (cos[(i − 1)π/2], sin[(i − 1)π/2]) i = 1, 2, 3, 4 (2.1)

√
2(cos[(i − 5)π/2 + π/4], sin[(i − 5)π/2 + π/4]) i = 5, 6, 7, 8.
The evolution equation of the density distribution function reads
1£ ¤
f i (x + cei 1t, t + 1t) − f i (x, t) = − f i (x, t) − f i(0) (x, t) , (2.2)
τ
where c = 1x/1t, 1x, and 1t are the lattice grid spacing and the time step, respectively;
τ is the dimensionless relaxation time; and f i(0) (x, t) is the equilibrium density distribution
function, which is determined by
f i(0) (x, t) = ωi ρ + ρsi (u(x, t)), (2.3)

where
· ¸
(ei · u) (ei · u)2 |u|2
si (u) = ωi 3 + 4.5 − 1.5 (2.4)
c c2 c2
290 GUO, SHI, AND WANG

with the weight coefficient

4

 i =0
9
ωi = 1
i = 1, 2, 3, 4


9
1
36
i = 5, 6, 7, 8.

The velocity in the above equilibrium distribution function is required to be small; i.e.,
|u|/c ≈ M ¿ 1, where M is the Mach number. The fluid density ρ and velocity u are
obtained from the density distribution function f i (x, t):
X
ρ= fi (2.5)
i
X
ρu = cei f i . (2.6)
i

The mass and momentum equations can be derived from the model via multiscaling expan-
sion as
∂ρ
+ ∇ · (ρu) = 0 (2.7)
∂t
∂(ρu)
+ ∇ · (ρuu) = −∇ p + ν[∇ 2 (ρu) + ∇(∇ · (ρu))], (2.8)
∂t
√
where p = cs2 ρ is the pressure, cs = c/ 3 is the sound speed, and ν = (2τ − 1)c2 1t/6 is
the kinematic viscosity. Clearly, the mass and momentum equations are exactly the same
as the compressible Navier–Stokes equations if the density variation is small enough.

2.2. Incompressible d2q9 LBGK Model

The d2q9 LBGK model described in the above section requires two conditions: (i) small
Mach number and (ii) small density variation. That is, the model is only appropriate for
compressible flows in incompressible limit; the incompressible Navier–Stokes equations
cannot be recovered properly. In this section, we will give a novel incompressible LBGK
model from which the incompressible Navier–Stokes equations can be exactly recovered
with the small Mach number limit.
The discrete velocity directions of the model are the same as those of the d2q9 model, and
a new type of distribution function gi (x, t) is introduced with the equilibrium distribution
function gi(0) (x, t) defined by


 −4σ cp2 + s0 (u) i = 0

gi(0) = λ cp2 + si (u) i = 1, 2, 3, 4 (2.9)


 p
γ c2 + si (u) i = 5, 6, 7, 8,

where σ , λ, and γ are parameters satisfying

λ+γ = σ
1
λ + 2γ = .
2
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 291

The distribution function satisfies the following conservation laws:

X
8 X
8
gi = gi(0) (2.10)
i=0 i=0

X
8 X
8
cei gi = cei gi(0) . (2.11)
i=0 i=0

The evolution equation of the system is

1¡ ¢
gi (x + cei 1t, t + 1t) − gi (x, t) = − gi (x, t) − gi(0) (x, t) . (2.12)
τ
The velocity and pressure of flow are given by

X
8
u= cei gi (2.13)
i=1
" 8 #
c2 X
p= gi + s0 (u) . (2.14)
4σ i=1

Through multiscaling expansion, the incompressible Navier–Stokes equations can be

derived from this incompressible LBGK model as (see the Appendix for details)

∇ ·u = 0 (2.15)
∂u
+ ∇ · (uu) = −∇ p + ν∇ 2 u, (2.16)
∂t
where the kinematic viscosity is determined by

(2τ − 1) (1x)2
ν= . (2.17)
6 1t

3. NUMERICAL RESULTS

To test the incompressible LBGK model proposed in the above section, numerical simula-
tions of the steady plane Poiseuille flow with pressure boundary condition, the unsteady 2-D
shear decaying flow, the driven cavity flow with Re = 400, 1000, 2000, and 5000 and the
flow around a circular cylinder with Re = 10, 20, 40, and 82.67 are performed. In the simula-
tions, the parameters of the incompressible LBGK model are taken as σ = 5/12, λ = 1/3,
and γ = 1/12. The distribution function gi is initialized by setting to equal gi(0) for all
nodes at t = 0. For boundary conditions, the extrapolation scheme proposed by Chen et al.
[3] is used. The original scheme uses second-order extrapolation, which is consistent with
LBGK methods. However, we found that the second-order extrapolation scheme has poor
stability for higher Reynolds numbers. Therefore, we used the second-order extrapola-
tion scheme for boundary conditions in simulating the Poiseuille flow, the shear decay-
ing flow, and the flow around a circular cylinder with Re = 10, 20, and 40, but used the
first-order extrapolation scheme for the flow around a circular cylinder with Re = 82.67
and the driven cavity flow with different Reynolds numbers considered for the sake of
stability.
292 GUO, SHI, AND WANG

3.1. Steady Plane Poiseuille Flow

A steady plane Poiseuille flow with pressure boundary conditions is defined in the region
{(x, y) | 0 ≤ x ≤ 2, 0 ≤ y ≤ 1} with the initial and boundary conditions

u(x, y, 0) = v(x, y, 0) = 0; p(x, y, 0) = p0 ;

u(x, 0, t) = u(x, 1, t) = v(x, 0, t) = v(x, 0, t) = 0;
p(0, y, t) = pin ; p(1, y, t) = pout ,

where p0 = 0.5( pin + pout ), and pin and pout are the pressure maintained at the entrance
and the exit, respectively. The Poiseuille flow has an analytic solution,

1p
u(x, y, t) = y(1 − y)
4ν
v(x, y, t) = 0 (3.1)
1p
p(x, y, t) = pin − x,
2
where 1p = pin − pout .
The parameters used in the simulations are: 1x = 1/16, ν = 1.0, pin = 1.1, pout = 1.0.
The lattice size is 16 × 32. To reach the steady state, a number of iterations are performed.
The criterion of steady state is
P ¯ (n+1) ¯
¯ − u i(n) ¯
i j ui j
≤ 5.0 × 10−9 ,
j
P ¯¯ (n+1) ¯¯
i j ui j

where u inj = u(xi , y j , n1t). Several different values of τ are used in the simulations. With
the fixed kinematic viscosity ν, the corresponding time steps are determined by Eq. (2.17).
It should be pointed out that the profiles of u along different horizontal lines are almost the
same, and the profiles of pressure p along different vertical lines are also almost the same.
Figure 1 shows profiles of the horizontal component of velocity, u, at the mid-width of the
channel, x = 1.0, and the pressure, p, at the mid-height of the channel, y = 0.5. In all the
simulations, the vertical component of the velocity, v, are found to be smaller than 10−10 .
We can see that the simulation results of the velocity profiles are of parabolic shape along
the channel, and the pressure distribution is linear along the channel.

3.2. Unsteady Decaying Shear Flow

The unsteady decaying shear flow has an analytic solution,

u(x, y, t) = A
v(x, y, t) = B cos(kx − k At)e−k νt
2
(3.2)
p(x, y, t) = p(t),

where the constants A, B, k, ν, and p(t) are chosen as A = B = k = ν = 1.0, p(t) =

p0 + 0.01 sin(20πt), and p0 = 3.0 is the average pressure. The domain is defined in the
region {(x, y) | 0 ≤ x ≤ 2π, 0 ≤ y ≤ 0.16π }. The initial and boundary conditions are im-
plemented according to the analytic solution given by Eqs. (3.2).
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 293

FIG. 1. Profiles of velocity u (left), at the mid-width of the channel (x = 1.0), and pressure p (right), at the
mid-height of the channel (y = 0.5), in steady plane Poiseuille flow with different relaxation times. Solid lines:
analytic solutions. ×: Simulation results. (a): 1/τ = 0.8; (b): 1/τ = 1.0; (c): 1/τ = 1.2.

In the simulation, the lattice size is fixed at 8 × 100. Data at several different times are
measured with fixed relaxation time τ such that 1/τ = 0.95. We observed that in the simu-
lation, the relative global errors of the horizontal component of the velocity, u, are of order
10−4 , and the profiles of the vertical component of the velocity, v, are similar for different
values of y. Figure 2 shows the simulation results of the vertical component of velocity, v,
at y = 0.08π. It shows that the simulation results agree well with the analytic solution.
To compare with the existing incompressible LBGK model, the relative global error of u
and v at different times in the simulation are also measured, as listed in Table I. The relative
global error is given by
P ¡ (n) (n) ¢2
i, j u i j − ū i j
ke(u)k =2
P ¡ (n) ¢2 , (3.3)
i, j ū i j
294 GUO, SHI, AND WANG

FIG. 2. Velocity profiles, v(x, 0.08π, t), of unsteady decaying shear flow. Lattice size: 8 × 100. 1/τ = 0.95.
Solid lines: analytic solutions. ×: simulation results. a: t = 0.3684211; b: t = 0.7368421; c: t = 1.105263; d: t =
1.473684; e: t = 1.842105.

where the summation is over the entire system, and ū is the analytic solution given by
Eq. (3.2). In Table I, u and v are simulation results obtained using the incompressible LBGK
model presented in this paper, and u 1 and v1 are results by using He’s incompressible LBGK
model [9]. In the simulation, the average pressure p0 of He’s model is set to be 3.0. Table I
shows that the relative global errors produced by the model proposed in this paper vary
slowly at different times, and they are much smaller than the errors produced by He’s model.

3.3. Driven Cavity Flow

Numerical simulations for the driven cavity flow have been studied by many authors
using traditional schemes such as finite difference [5, 18] and multigrid [8, 19] methods. The
lattice Boltzmann simulation was performed by Hou and Zou [12] with detailed analysis.
The configuration of driven cavity flow considered here consists of a two-dimensional
square cavity whose top plate moves from left to right with constant velocity, while the

TABLE I
Relative Global Error of Velocity Field in 2-D Decaying Shear Flow

0.3684211 0.7368421 1.105263 1.473684 1.842105

e(u) × 10 4
3.5252809 2.3628879 1.5602352 1.0394268 0.71491523
e(v) × 103 2.5844411 2.5898510 2.5948778 2.5968753 2.5945725
e(u 1 ) × 103 0.80553308 1.1336378 1.6208783 0.2162959 1.4642498
e(v1 ) × 103 3.8908607 7.2841784 14.286716 3.3058078 27.141689
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 295

other three boundaries are fixed. The fundamental characteristics of the driven cavity flow
are the emergence of a large primary vortex in the center and of two secondary vortices in
the lower two corners. The values of the stream function and the locations of the centers of
the vortices are functions of the Reynolds number, which is determined by Re = LU/ν,
where L is the height or width of the cavity, U is the uniform velocity of the top plate, and
ν is the kinematic viscosity.
Numerical simulations are carried out using the present incompressible LBGK model for
the driven cavity flow with Re = 400, 1000, 2000, and 5000 on a 257 × 257 lattice. The re-
laxation parameter ω = 1/τ is set to be 1.5, 1.75, 1.85, and 1.85, respectively. The flow with
Re = 400 is first simulated, where the initial condition is set as u = v = 0, p = 1. The sim-
ulations for Re = 1000, 2000, and 5000 start from the steady states for Re = 400, 1000, and
2000, respectively. In the simulations, steady state is reached when the difference between
the velocities at the center of the cavity for successive 1,000 steps is less than 5 × 10−6 .
Figure 3 shows the contours of the stream function of the R flow for the Reynolds numbers
considered. The stream function, ψ, is defined as ψ = v d x − u dy, and is calculated
from the discrete velocity field obtained from the LBGK simulation. In this paper, the
integral is calculated in the y direction using the Simpson’s rule with zero value on the
bottom boundary. These plots show clearly the effect of the Reynolds number on the flow
pattern. For low Re(≤1000), only three vortices appear in the cavity, a primary one near
the center and a pair of secondary ones in the lower corners of the cavity. At Re = 2000, a
third secondary vortex is seen in the upper left corner. When Re reaches to 5000, a tertiary

FIG. 3. Streamlines of the driven cavity flow for different Reynolds numbers. (a) Re = 400; (b) Re = 1000;
(c) Re = 2000; (d) Re = 5000.
296 GUO, SHI, AND WANG

FIG. 4. The velocity components, u and v, along the vertical and horizontal lines through the cavity center for
the driven cavity flow at different Reynolds numbers. (a) Re = 400; (b) Re = 1000; (c) Re = 2000; (d) Re = 5000.

vortex appears in the lower right corner. We can also see that the center of the primary
vortex moves toward the center of the cavity as Re increases. The velocity components, u
and v, along the vertical and horizontal center lines for different Re are shown in Fig. 4.
The profiles are found to become near linear in the center core of the cavity as Re becomes
large. These observations show that the present LBGK simulation is in agreement with the
previous studies [8, 12, 18, 19].
To quantify the results, the strength and locations of the primary center vortex and the two
secondary ones are listed in Table II. From the table, we can see that the strength and locations
of the vortices predicted by the LBGK model agree well with those of previous work for all
the Reynolds numbers considered. The velocity components along the vertical and horizon-
tal lines through the cavity center are also compared with Ghia’s benchmark solutions [8]; the
maximum, average, and relative global error are listed in Table III. Here the relative global

TABLE II
Strength and Locations of Vortex of the Driven Cavity Flow: (·)c Primary Vortex; (·)l Lower
Left Vortex; (·)r Lower Right Vortex

Re ψc xc yc ψl xl yl ψr xr yr

400 a 0.1136 0.5563 0.6000 −1.46e-5 0.0500 0.0500 −6.45e-4 0.8875 0.1188
b 0.1139 0.5547 0.6055 −1.42e-5 0.0508 0.0469 −6.42e-4 0.8906 0.1250
c 0.1121 0.5608 0.6078 −1.30e-5 0.0549 0.0510 −6.19e-4 0.8902 0.1255
d 0.1126 0.5547 0.6094 −1.36e-5 0.0508 0.0469 −6.23e-4 0.8867 0.1250
1000 a 0.1173 0.5438 0.5625 −2.24e-4 0.0750 0.0813 −1.74e-3 0.8625 0.1063
b 0.1179 0.5313 0.5625 −2.31e-4 0.0859 0.0781 −1.75e-3 0.8594 0.1094
c 0.1178 0.5333 0.5647 −2.22e-4 0.0902 0.0784 −1.69e-3 0.8667 0.1137
d 0.1170 0.5313 0.5625 −2.21e-4 0.0859 0.0781 −1.68e-3 0.8672 0.1172
2000 a 0.1116 0.5250 0.5500 −6.90e-4 0.0875 0.1063 −2.60e-3 0.8375 0.0938
c 0.1204 0.5255 0.5490 −7.26e-4 0.0902 0.1059 −2.44e-3 0.8471 0.0980
d 0.1186 0.5234 0.5469 −7.13e-4 0.0898 0.1016 −2.40e-3 0.8438 0.1016
5000 a 0.0921 0.5125 0.5313 −1.67e-3 0.0625 0.1563 −5.49e-3 0.8500 0.0813
b 0.1190 0.5117 0.5352 −1.36e-3 0.0703 0.1367 −3.08e-3 0.8086 0.0742
c 0.1214 0.5176 0.5373 −1.35e-3 0.0784 0.1373 −3.03e-3 0.8078 0.0745
d 0.1120 0.5159 0.5391 −1.29e-3 0.0781 0.1328 −2.85e-3 0.8086 0.0781

Note. a, Vanka [19]; b, Ghia et al. [8]; c, Hou and Zou [12]; d, present work.
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 297

TABLE III
Errors of the Velocity Components through the Cavity Center: emax ,
Maximum Error; eave , Average Error; erg , Relative Global Error

Re

400 1000 5000

u(L/2, y) emax 1.346 × 10−2 2.266 × 10−2 1.884 × 10−2

eave 4.280 × 10−3 5.942 × 10−3 9.877 × 10−3
erg 1.450 × 10−2 2.219 × 10−2 2.818 × 10−2
v(x, L/2) emax 3.837 × 10−3 4.976 × 10−3 2.405 × 10−2
eave 2.268 × 10−3 1.413 × 10−3 1.047 × 10−2
erg 1.027 × 10−2 6.075 × 10−3 3.258 × 10−2

error is calculated as in Eq. (3.3) with ū replaced by Ghia’s benchmark solutions, and the
summation is over the corresponding lines. The relative global errors about these velocity
components are found to be less than 3.3% for all values of Re (data for Re = 2000 were
not given in [8]) in comparision with Ghia’s benchmark solutions. This indicates that the
agreement between the present LBGK solutions and Ghia’s benchmark solutions [8] is
satisfactory.

3.4. Flow Around a Circular Cylinder

Although the flow in the square cavity is complex, the geometry is nevertheless simple
because only flat boundaries are involved. To demonstrate the capability of the present
incompressible LBGK model, we apply the model to the two-dimensional flow past a
circular cylinder. The flow has been studied using lattice Boltzmann methods based on
uniform meshes [3, 11, 20, 21] or nonuniform meshes [6, 10, 15] by several groups.
The cylinder is set into a channel with uniform velocity U at the inlet in the streamwise
direction x, with origin at the center of the cylinder. The computation domain considered
here is −2.5D ≤ x ≤ 11.5D and −3.5D ≤ y ≤ 3.5D, where D is the diameter of the
cylinder. The boundary conditions are as follows:
inlet
u = U, v = 0

outlet
∂u ∂v
= =0
∂x ∂x
top and bottom
∂u
= 0, v = 0
∂y
solid surface of the cylinder

u = v = 0.

The Reynolds number(Re = DU/ν) of the flow is based on the uniform inlet velocity
U and the cylinder diameter D. Numerical simulations are carried out for steady flow with
298 GUO, SHI, AND WANG

FIG. 5. Streamlines of the steady flow around a circular cylinder for different Reynolds numbers.

small Reynolds numbers (Re = 10, 20, 40) and unsteady flow with Re = 82.67 based on
an N y × N x = 211 × 421 grid. In the simulations, the surface of the circular cylinder is
approximated by the lattice nodes closest to it. The relaxation parameter ω = 1/τ is set to
be 0.7 for Re = 10, 1.2 for Re = 20 and 40, and 1.6 for Re = 82.67, respectively.
Figure 5 shows the streamlines of the flow when it reaches its final steady state. In all
cases, a pair of stationary recirculating eddies appear behind the cylinder. The wake length,
L, the distance from the rearmost point of the cylinder to the end of the wake, and the
separation angle, θ, are measured. The quantitative geometrical parameters are listed in
Table IV as well as related previous computational and experimental data. Both the wake
length and the separation angle agree well with the results of previous studies for all the
three Reynolds numbers considered.
For Re = 82.67, a lattice with larger size should be used. However, we still use the
211 × 421 lattice because of limited computer resources. In the simulation, a periodic
vortex shedding flow is obtained after a sufficient number of iterations. Figure 6 shows the
instantaneous streamlines at different times within a cycle of the periodic vortex shedding.
For comparison, this flow is also simulated using He’s model [9] with the same lattice
parameters and boundary conditions. A final periodic vortex shedding flow is also obtained,
which is quite similar to the present incompressible LBGK solution. Figures 7 and 8 show the
instantaneous pressure fields and the velocity components along the horizontal and vertical
lines through the cylinder center. Close agreement between the two LBGK solutions can
be observed from these figures. Table V lists the maximum, average, and relative global
errors of the velocity components corresponding to Fig. 8. The relative global errors for the
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 299

TABLE IV
Geometrical Parameter of Flow around Circular Cylinder

Re

10 20 40

2L/D θ 2L/D θ 2L/D θ

a 0.434 27.96 1.786 43.37 4.357 53.34

b 0.68 32.5 1.86 44.8 4.26 53.5
c 0.474 26.89 1.842 42.9 4.490 52.84
d 0.498 30.0 1.804 42.1 4.38 50.12
e 0.533 31.61 1.867 42.27 4.400 53.13

Note. a, Nieuwstadt and Keller [16]; b, Coutanceau and Bouard [4]; c, He and
Doolen [10]; d, Mei and Shyy [15]; e, present work.

FIG. 6. Instantaneous streamlines of the flow around a circular cylinder for Re = 82.67.
300 GUO, SHI, AND WANG

FIG. 7. Comparison of the pressure fields of the flow around a circular cylinder. (a) Present work; (b) He’s
[9] model.

FIG. 8. Comparison of the velocity components of the flow around a circular cylinder. (a) u(x, 0); (b) v(x, 0);
(c) u(0, y); (d) v(0, y). Solid lines: He’s model; ×: Present work.
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 301

TABLE V
Errors of the Velocity Components through the Cylinder Center: emax ,
Maximum Error; eave , Average Error; erg , Relative Global Error

u(x, 0) v(x, 0) u(0, y) v(0, y)

emax 2.5614 × 10−3 1.6127 × 10−3 3.0158 × 10−4 1.0618 × 10−4

eave 8.2109 × 10−4 5.8597 × 10−4 1.0825 × 10−4 4.2629 × 10−5
erg 1.5685 × 10−2 2.1428 × 10−2 1.2034 × 10−3 3.9605 × 10−3

velocity components are found to be less than 2.2% between the simulation results of the
present LBGK model and He’s model [9].

4. SUMMARY AND CONCLUSION

In the above sections, a LBGK model for incompressible flows has been developed.
Through multiscaling expansion, the Navier–Stokes equations are precisely recovered from
the model accurate to the order of O(² 2 ) in the continuity equation and O(² 2 + ² M 2 ) in the
momentum equation. By choosing the distribution function appropriately, the compressible
effect is eliminated in the incompressible LBGK model.
Numerical simulations are performed to test the performance of the incompressible
LBGK model. The test problems include the plane Poiseuille flow condition, the unsteady
shear decaying flow, the driven cavity flow, and the flow around a circular cylinder. For
the steady plane Poiseuille flow with constant pressure gradient boundary condition, the
simulation results agree well with the analytic solution of the incompressible Navier–Stokes
equations.
As to the unsteady decaying shear flow, simulation is measured at several different times,
and all of the results are in excellent agreement with the analytic solutions of the problem.
The relative global errors of the results are also measured to compare the incompressible
LBGK model presented in this paper with He’s model [9]. It shows that our model has
better performance than He’s in dealing with this problem. The reason might be that our
model completely eliminates the compressible effect, whereas He’s model only reduces it.
Another advantage of our model is that the average pressure does not need to be specified
in advance; this enables us to use the model for general incompressible flows.
For the driven cavity flow, simulation is carried out with four different Reynolds num-
bers.The LBGK model correctly predicts the emergence of the third secondary vortex near
the upper left corner for Re = 2000 and the tertiary vortex in the lower right corner for
Re = 5000. The strength and locations of the primary and two secondary vortices in the
cavity are measured. Compared with the results of previous studies, the LBGK simulation
for the driven cavity flow is satisfactory.
For flow around a circular cylinder, steady-state solutions are obtained for Re = 10, 20,
and 40. The wake length and separation angle of the steady flow agree well with previous
experimental and computational data. For Re = 82.67, the present model predicts the peri-
odic vortex shedding flow, and the simulation results are in excellent agreement with those
obtained from He’s model for the same problem.
The results of these numerical tests indicates the capability of the present incompressible
LBGK model in handling steady and unsteady flows. The close agreement with the ana-
lytic solutions or experimental or computational data of previous studies shows the good
302 GUO, SHI, AND WANG

performance of the model. In conclusion, we have developed a LBGK model for the incom-
pressible Navier–Stokes equations. The model completely eliminates the compressibility
effect that lies in other existing LBGK models and can be used for both steady and unsteady
flows.

APPENDIX: RECOVERING THE INCOMPRESSIBLE NAVIER–STOKES EQUATIONS

It is well known that in incompressible flows,

O(δp) = O(δρ) = O(M 2 ) (A1a)

O(u) = O(M) (A1b)

in the limit M → 0, where δp and δρ are the pressure and density fluctuations, respectively.
To derive the incompressible Navier–Stokes equations from the LBGK model presented
in this paper, we first expand the distribution function as

gi = gi(0) + ²gi(1) + ² 2 gi(2) + · · · , (A2)

where ² = 1t. From Eqs. (2.10) and (2.11), we know that

X
8
gi(m) = 0 (A3a)
i=0

X
8
cei gi(m) = 0 (A3b)
i=0

for m = 1, 2, . . ..
With the multiscaling expansion (A2), the evolution equation (2.12) of the incompressible
LBGK model can be written in continuous form as
µ ¶
1 2 (0) 1
Di gi(0) + ² Di gi + Di gi(1) + O(² 2 ) = − gi(1) , (A4)
2 τ

where Di = ( ∂t∂ + cei · ∇). And therefore,

gi(1) = −τ Di gi(0) + O(²). (A5)

Applying Eq. (A5) to the left hand of Eq. (A4), we can rewrite Eq. (A4) as
µ ¶
(0) 1 (1) 1
Di gi = − gi + ² τ − Di2 gi(0) + O(² 2 ). (A6)
τ 2

Note that E(2n+1) = 0 for n = 0, 1, . . . , where E(n) are the tensors defined as E(n) =
P
α eα1 eα2 · · · eαn , and

X
4
eiα eiβ = 2δαβ (A7a)
i=1

X
8
eiα eiβ = 4δαβ (A7b)
i=5
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 303

X
4
eiα eiβ eiζ eiθ = 2δαβζ θ (A7c)
i=1

X
8
eiα eiβ eiζ eiθ = 41αβζ θ − 8δαβζ θ , (A7d)
i=5

where δαβ and δαβζ θ are the Kronecker tensors, and

1αβζ θ = δαβ δζ θ + δαζ δβθ + δαθ δβζ . (A8)

With these identical equations, we obtain the moments of gi(0) :

X
gi(0) = 0 (A9a)
i
X
cei gi(0) = u (A9b)
i
X
5(0) =: c2 ei ei gi(0) = uu + p I (A9c)
i

(0)
X c2
0αβζ =: c3 eiα eiβ eiζ gi(0) = (δαβ uζ + δαζ uβ + δβζ uα ). (A9d)
i
3

With the aids of Eqs. (A9), the moments of Eq. (A6),

X µ ¶X
1
Di gi(0) =² τ− Di2 gi(0) + O(² 2 ) (A10a)
i
2 i

and
X µ ¶X
1
cei Di gi(0) =² τ− cei Di2 gi(0) + O(² 2 ), (A10b)
i
2 i

lead to the Euler equations,

∇ · u = 0 + O(²) (A11a)
∂u
+ ∇ · 5(0) = 0 + O(²), (A11b)
∂t

to the O(²) order, and the following equations,

¶ µ
1
∇ ·u=² τ − P + O(² 2 ) (A12a)
2
µ ¶
∂u (0) 1
+∇ ·5 =² τ − Q + O(² 2 ), (A12b)
∂t 2

to the O(² 2 ) order, where

X µ ¶
∂ ∂u
P= Di2 gi(0) = (∇ · u) + ∇ · + ∇ · 5(0) = O(²)
i
∂t ∂t
304 GUO, SHI, AND WANG

and
X µ ¶ µ ¶
∂ ∂u ∂5(0)
Q= cei Di2 gi(0) = + ∇ · 50 +∇ · + ∇ · 0 (0)
i
∂t ∂t ∂t
µ ¶
∂5(0)
c2 2
= O(²) + ∇ · +
(∇ u + 2∇(∇ · u))
3∂t
µ ¶
∂ p ∂(uu) c2
= O(²) + ∇ · + + (∇ 2 u)
∂t ∂t 3
c2 2
= O(²) + O(M 2 ) + (∇ u).
3

The term ∂∂tp + ∂(uu)

∂t
is of order O(M 2 ) due to Eqs. (A1). By applying the above results
about P and Q to Eqs. (A12), the incompressible Navier–Stokes equations are derived
accurate to the order of O(² 2 ) in the continuity equation and O(² 2 + ² M 2 ) in the momentum
equation

∇ · u = 0 + O(² 2 ) (A13a)
∂u
+ ∇ · (uu) = −∇ p + ν∇ 2 u + O(² 2 + ² M 2 ), (A13b)
∂t

where the kinematic viscosity is determined by

(2τ − 1) (1x)2
ν= . (A14)
6 1t

Now we discuss how to compute the pressure p from the distribution function. According
to the expression of g0(0) , we have

4σ
p(x, t + ²) = −g0(0) (x, t + ²) + s0 (u(x, t + ²))
c2
= −g0 (x, t + ²) + ²g0(1) (x, t + ²) + s0 (u(x, t + ²))
∂g0(1)
= −g0 (x, t + ²) + ²g0(1) (x, t) + s0 (u(x, t + ²)) + ² 2
∂t
= −g0 (x, t + ²) + ²g0(1) (x, t) + s0 (u(x, t + ²)) + O(² 2 ),

while the term ²g0(1) (x, t) can be neglected because from the evolution equation we have

1
g0(1) (x, t) = − [g0 (x, t + ²) − g0 (x, t)]
τ²
· ¸
∂g0 (x, t)
= −τ + O(²)
∂t
" #
∂g0(0) (x, t)
= −τ + O(²)
∂t
· ¸
∂ p(x, t) ∂s0 (u(x, t))
= τ 4σ − + O(²).
∂t ∂t
 LBGK FOR INCOMPRESSIBLE NAVIER–STOKES EQUATION 305

Note that ∂p
∂t
= O(M 2 ) and s0 (u) = O(|u|)2 = O(M 2 ), so g0(1) (x, t) = O(² + M 2 ).
Therefore,

4σ
p(x, t + ²) = −g0 (x, t + ²) + s0 (u(x, t + ²)) + O(² 2 + ² M 2 ) (A15a)
c2
or

4σ X 8

2
p(x, t + ²) = gi (x, t + ²) + s0 (u(x, t + ²)) + O(² 2 + ² M 2 ) (A15b)
c i=1

due to Eqs. (2.10) and (A9a). As a result, the pressure p can be computed as
" 8 #
c2 X
p= gi + s0 (u) (A16)
4σ i=1

with accuracy of order O(² 2 + ² M 2 ), which is consistent with the order of the incompress-
ible Navier–Stokes equations (A13).

ACKNOWLEDGMENTS

The author Zhaoli Guo thanks Professor Haiping Fang for his helpful discussions. We also thank the reviewers
for comments which greatly improved the work. This work is supported in part by the National Key Project of Fun-
damental Research, the Ministry of Science and Technology, People’s Republic of China (Grant G1999022207),
and the National Science Foundation (Grant 59676030).

REFERENCES

1. Reference removed in proofs.

2. S. Chen and G. Doolen, Lattice Boltzmann method for fluid flows, Ann. Rev. Fluid Mech. 30, 329 (1998).
3. S. Chen, D. Martinez, and R. Mei, On boundary conditions in lattice Boltzmann methods, Phys. Fluids. 8,
2527 (1996).
4. M. Coutanceau and R. Bouard, Experimental determination of the main features of the viscous flow in the
wake of a circular cylinder in uniform translation. 1. Steady flow, J. Fluid Mech. 79, 231 (1977).
5. W.-N. E, and J. Liu, Essential compact scheme for unsteady viscous incompressible flows, J. Comput. Phys.
126, 122 (1996).
6. O. Filippova and D. Hänel, Grid refinement for lattice-BGK models, J. Comput. Phys. 147, 219 (1998).
7. U. Frisch, D. d’Humiéres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.-P. Rivet, Lattice gas hydrodynamics
in two and three dimensions. Complex Syst. 1, 649 (1987).
8. U. Ghia, K. N. Ghia, and C. T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes
equations and a multigrid method, J. Comput. Phys. 48, 387 (1982).
9. X. He and L.-S. Luo, Lattice Boltzmann model for the incompressible Navier–Stokes equation, J. Stat. Phys.
88, 927 (1997).
10. X. He and G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: Flow around a circular
cylinder, J. Comput. Phys. 134, 306 (1997).
11. F. J. Higuera and S. Succi, Simulating the flow around a circular cylinder with a lattice Boltzmann equation,
Europhys. Lett. 8, 517 (1989).
12. S. Hou and Q. Zou, Simulation of cavity flow by the lattice Boltzmann method, J. Comput. Phys. 118, 329
(1995).
306 GUO, SHI, AND WANG

13. Z. Lin, H. Fang, and R. Tao, Improved lattice Boltzmann model for incompressible two-dimensional steady
flows, Phys. Rev. E 54, 6323 (1997).
14. D. O. Martinez, W. H. Matthaeus, S. Chen, and D. C. Montgomery, Comparison of spectral method and lattice
Boltzmann simulations of two-dimensional hydrodynamics, Phys. Fluids. 6, 1285 (1994).
15. R. Mei and Q. Shyy, On the finite difference-based lattice Boltzmann method in curvilinear coordinates,
J. Comput. Phys. 143, 426 (1998).
16. F. Nieuwstadt and H. B. Keller, Viscous flow past circular cylinders, Comput. Fluids. 1, 59 (1973).
17. Y. Qian, D. d’Humiéres, and P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett.
17, 479 (1992).
18. R. Schreiber and H. Keller, Driven cavity flow by efficient numerical techniques, J. Comput. Phys. 49, 310
(1983).
19. S. P. Vanka, Block-implicit multigrid solution of Navier–Stokes equations in primitive variables, J. Comput.
Phys. 65, 138 (1986).
20. L. Wagner, Pressure in lattice Boltzmann simulations of flow around a cylinder, Phys. Fluids. 6, 3516 (1994).
21. L. Wagner and F. Hayot, Lattice Boltzmann simulations of flow past a cylindrical obstacle, J. Stat. Phys. 81,
63 (1995).
22. Q. Zou, S. Hou, S. Chen, and G. Doolen, An improved incompressible lattice Boltzmann model for time-
independent flows, J. Stat. Phys. 81, 35 (1995).