P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL

November 24, 1997 15:23 Annual Reviews AR049-12

Annu. Rev. Fluid Mech. 1998. 30:32964
Copyright c 1998 by Annual Reviews Inc. All rights reserved
LATTICE BOLTZMANN METHOD
FOR FLUID FLOWS
Shiyi Chen
1,2
and Gary D. Doolen
2
1
IBM Research Division, T. J. Watson Research Center, P.O. Box 218, Yorktown
Heights, NY 10598;
2
Theoretical Division and Center for Nonlinear Studies, Los
Alamos National Laboratory, Los Alamos, NM 87545; e-mail: syc@cnls.lanl.gov
KEY WORDS: lattice Boltzmann method, mesoscopic approach, uid ow simulation
ABSTRACT
We present an overview of the lattice Boltzmann method (LBM), a parallel and
efcient algorithm for simulating single-phase and multiphase uid ows and
for incorporating additional physical complexities. The LBMis especially useful
for modeling complicated boundary conditions and multiphase interfaces. Recent
extensions of this method are described, including simulations of uid turbulence,
suspension ows, and reaction diffusion systems.
INTRODUCTION
In recent years, the lattice Boltzmann method (LBM) has developed into an
alternative and promising numerical scheme for simulating uid ows and
modeling physics in uids. The scheme is particularly successful in uid ow
applications involving interfacial dynamics and complex boundaries. Unlike
conventional numerical schemes based on discretizations of macroscopic con-
tinuum equations, the lattice Boltzmann method is based on microscopic mod-
els and mesoscopic kinetic equations. The fundamental idea of the LBM is
to construct simplied kinetic models that incorporate the essential physics of
microscopic or mesoscopic processes so that the macroscopic averaged prop-
erties obey the desired macroscopic equations. The basic premise for using
these simplied kinetic-type methods for macroscopic uid ows is that the
macroscopic dynamics of a uid is the result of the collective behavior of many
microscopic particles in the system and that the macroscopic dynamics is not
sensitive to the underlying details in microscopic physics (Kadanoff 1986).
By developing a simplied version of the kinetic equation, one avoids solving
329
0066-4189/98/0115-0329$08.00

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
330 CHEN & DOOLEN
complicated kinetic equations such as the full Boltzmann equation, and one
avoids following each particle as in molecular dynamics simulations.
Even though the LBM is based on a particle picture, its principal focus is
the averaged macroscopic behavior. The kinetic equation provides many of
the advantages of molecular dynamics, including clear physical pictures, easy
implementation of boundary conditions, and fully parallel algorithms. Because
of the availability of very fast and massively parallel machines, there is a current
trend to use codes that can exploit the intrinsic features of parallelism. The LBM
fullls these requirements in a straightforward manner.
The kinetic nature of the LBM introduces three important features that dis-
tinguish it from other numerical methods. First, the convection operator (or
streaming process) of the LBM in phase space (or velocity space) is linear.
This feature is borrowed from kinetic theory and contrasts with the nonlinear
convection terms in other approaches that use a macroscopic representation.
Simple convection combined with a relaxation process (or collision operator)
allows the recovery of the nonlinear macroscopic advection through multi-scale
expansions. Second, the incompressible Navier-Stokes (NS) equations can be
obtained in the nearly incompressible limit of the LBM. The pressure of the
LBM is calculated using an equation of state. In contrast, in the direct nu-
merical simulation of the incompressible NS equations, the pressure satises a
Poisson equation with velocity strains acting as sources. Solving this equation
for the pressure often produces numerical difculties requiring special treat-
ment, such as iteration or relaxation. Third, the LBM utilizes a minimal set of
velocities in phase space. In the traditional kinetic theory with the Maxwell-
Boltzmann equilibrium distribution, the phase space is a complete functional
space. The averaging process involves information from the whole velocity
phase space. Because only one or two speeds and a few moving directions are
used in LBM, the transformation relating the microscopic distribution func-
tion and macroscopic quantities is greatly simplied and consists of simple
arithmetic calculations.
The LBM originated from lattice gas (LG) automata, a discrete particle ki-
netics utilizing a discrete lattice and discrete time. The LBM can also be
viewed as a special nite difference scheme for the kinetic equation of the
discrete-velocity distribution function. The idea of using the simplied kinetic
equation with a single-particle speed to simulate uid ows was employed by
Broadwell (Broadwell 1964) for studyingshockstructures. Infact, one canview
the Broadwell model as a simple one-dimensional lattice Boltzmann equation.
Multispeed discrete particle velocities models have also been used for studying
shock-wave structures (Inamuro & Sturtevant 1990). In all these models, al-
though the particle velocity in the distribution function was discretized, space
and time were continuous. The full discrete particle velocity model, where

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 331
space and time are also discretized on a square lattice, was proposed by Hardy
et al (1976) for studying transport properties of uids. In the seminal work of
the lattice gas automaton method for two-dimensional hydrodynamics, Frisch
et al (1986) recognized the importance of the symmetry of the lattice for the re-
covery of the Navier-Stokes equation; for the rst time they obtained the correct
Navier-Stokes equation starting from the lattice gas automata on a hexagonal
lattice. The central ideas in the papers contemporary with the FHP paper in-
clude the cellular automaton model (Wolfram 1986) and the 3-D model using
the four-dimensional face-centered-hyper-cubic (FCHC) lattice (dHumi` eres
et al 1986).
The lattice gas automaton is constructed as a simplied, ctitious molecular
dynamic in which space, time, and particle velocities are all discrete. From this
perspective, the lattice gas method is often called lattice gas cellular automata.
In general, a lattice gas automaton consists of a regular lattice with particles
residing on the nodes. A set of Boolean variables n
i
(x, t )(i = 1, , M)
describing the particle occupation is dened, where M is the number of direc-
tions of the particle velocities at each node. The evolution equation of the LG
automata is as follows:
n
i
(x +e
i
, t +1) = n
i
(x, t ) +
i
(n(x, t )), (i = 0, 1 , M), (1)
where e
i
are the local particle velocities. Starting from an initial state, the
conguration of particles at each time step evolves in two sequential sub-steps,
(a) streaming, in which each particle moves to the nearest node in the direction
of its velocity, and (b) collision, which occurs when particles arriving at a node
interact and change their velocity directions according to scattering rules. For
simplicity, the exclusion principle (no more than one particle being allowed at
a given time and node with a given velocity) is imposed for memory efciency
and leads to a Fermi-Dirac local equilibrium distribution (Frisch et al 1987).
The main feature of the LBMis to replace the particle occupation variables, n
i
(Boolean variables), in Equation 1 by single-particle distribution functions (real
variables) f
i
= n
i
and neglect individual particle motion and particle-particle
correlations in the kinetic equations (McNamara & Zanetti 1988), where
denotes an ensemble average. This procedure eliminates statistical noise in the
LBM. Inthe LBM, the primitive variables are the averagedparticle distributions,
which are mesoscopic variables. Because the kinetic form is still the same as
the LGautomata, the advantages of locality in the kinetic approach are retained.
The locality is essential to parallelism.
An important simplication of the LBM was made by Higuera & Jim enez
(1989) who linearized the collision operator by assuming that the distribution
is close to the local equilibrium state. An enhanced collision operator approach
which is linearly stable was proposed by Higuera et al (1989). A particular

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
332 CHEN & DOOLEN
simple linearized version of the collision operator makes use of a relaxation time
towards the local equilibriumusinga single time relaxation. The relaxationterm
is known as the Bhatnagar-Gross-Krook (BGK) collision operator (Bhatnagar
et al 1954) and has been independently suggested by several authors (Qian
1990, Chen et al 1991). In this lattice BGK (LBGK) model, the local equilib-
riumdistribution is chosen to recover the Navier-Stokes macroscopic equations
(Qian et al 1992, Chen et al 1992). Use of the lattice BGK model makes the
computations more efcient and allows exibility of the transport coefcients.
LATTICE BOLTZMANN EQUATIONS
LBE: An Extension of LG Automata
There are several ways to obtain the lattice Boltzmann equation (LBE) from ei-
ther discrete velocity models or the Boltzmann kinetic equation. There are also
several ways to derive the macroscopic Navier-Stokes equations from the LBE.
Because the LBM is a derivative of the LG method, we will introduce the LBE
beginning froma discrete kinetic equation for the particle distribution function,
which is similar to the kinetic equation in the LG automata in Equation 1:
f
i
(x +e
i
x, t +t ) = f
i
(x, t ) +
i
( f (x, t )), (i = 0, 1 , M), (2)
where f
i
is the particle velocity distribution function along the i th direction;

i
=
i
( f (x, t )) is the collision operator which represents the rate of change of
f
i
resulting from collision. t and x are time and space increments, respec-
tively. When x/t = |e
i
|, Equations 1 and 2 have the same discretizations.

i
depends only on the local distribution function. In the LBM, space is dis-
cretized in a way that is consistent with the kinetic equation, i.e. the coordinates
of the nearest neighbor points around x are x +e
i
.
The density and momentum density u are dened as particle velocity
moments of the distribution function, f
i
,
=

i
f
i
, u =

i
f
i
e
i
, (3)
where

i


M
i =1,
.
i
is required to satisfy conservation of total mass and
total momentum at each lattice:

i
= 0,

i
e
i
= 0. (4)
If only the physics in the long-wavelength and low-frequency limit are of
interest, the lattice spacing x and the time increment t in Equation 2 can
be regarded as small parameters of the same order . Performing a Taylor

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 333
expansion in time and space, we obtain the following continuum form of the
kinetic equation accurate to second order in :
f
i
t
+e
i
f
i
+
_
1
2
e
i
e
i
: f
i
+e
i

f
i
t
+
1
2

2
f
i
t
2
_
=

i

. (5)
To derive the macroscopic hydrodynamic equation, we employ the Chapman-
Enskog expansion, which is essentially a formal multiscaling expansion (Frisch
et al 1987),
_

t
=

t
1
+
2

t
2
,

x
=

x
1
.
_
The above formula assumes that the diffusion time scale t
2
is much slower than
the convection time scale t
1
. Likewise, the one-particle distribution function
f
i
can be expanded formally about the local equilibrium distribution function
f
eq
i
,
f
i
= f
eq
i
+f
(neq)
i
. (6)
Here f
eq
i
depends on the local macroscopic variables ( and u) and should
satisfy the following constraints:

i
f
eq
i
= ,

i
f
eq
i
e
i
= u. (7)
f
(neq)
i
= f
(1)
i
+ f
(2)
i
+ O(
2
) is the nonequilibrium distribution function,
which has the following constraints:

i
f
(k)
i
= 0,

i
f
(k)
i
e
i
= 0, (8)
for both k = 1 and k = 2.
Inserting f
i
into the collision operator
i
, the Taylor expansion gives:

i
( f ) =
i
( f
eq
) +

i
( f
eq
)
f
j
f
(1)
j
+
2
_

i
( f
eq
)
f
j
f
(2)
j
+

2

i
( f
eq
)
f
j
f
k
f
(1)
j
f
(1)
k
_
+ O(
3
). (9)
From Equation 5, we note that when 0, we have:
i
( f
eq
) = 0. This leads
to a linearized collision operator,

i
( f )

=
M
i j

_
f
j
f
eq
j
_
, (10)
where M
i j


i
( f
eq
)
f
j
is the collision matrix (Higuera & Jim enez 1989), which
determines the scattering rate between directions i and j . For a given lattice,

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
334 CHEN & DOOLEN
M
i j
only depends on the angle between directions i and j and has a limited
set of values. For mass and momentum conservation collision, M
i j
satisfy the
following constraints (Benzi et al 1992):
M

i =1
M
i j
= 0,
M

i =1
e
i
M
i j
= 0. (11)
If we further assume that the local particle distribution relaxes to an equilib-
rium state at a single rate ,
M
i j
=
1

i j
, (12)
we arrive at the lattice BGK collision term (Bhatnagar et al 1954),

=
1

f
neq
i
=
1

_
f
(1)
i
+f
(2)
i
_
, (13)
and the LBGK equation:
f
i
(x +e
i
, t +1) = f
i
(x, t )
f
i
f
eq
i

. (14)
The BGK collision term (Qian 1990, Chen et al 1991, Qian et al 1992, Chen
et al 1992) was used previously in the full Boltzmann simulation for studying
shock formation (Chu 1965) and was used recently for a shock-capturing using
nite-volume methods (Xu & Prendergast). From Equation 5 one obtains the
following equations:
f
eq
i
t
1
+e
i

1
f
eq
i
=
f
(1)
i

, (15)
to order
0
and

t
1
f
(1)
i
+

t
2
f
eq
i
+e
i
f
(1)
i
+
1
2
e
i
e
i
: f
eq
i
+e
i


t
1
f
eq
i
+
1
2

2
t
2
1
f
eq
i
=
1

f
(2)
i
, (16)
to order
1
. Using Equation 15 and some algebra, we can rewrite the rst order
equation as
f
(1)
i
t
2
+
_
1
2

__
f
(1)
i
t
1
+e
i

1
f
(1)
i
_
=
f
(2)
i

. (17)
From Equations 15 and 17 we obtain the following mass and momentum equa-
tions:

t
+ u = 0, (18)

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 335
u
t
+ = 0, (19)
which are accurate to second order in for Equation 2. Here the momentum
ux tensor has the form:

i
(e
i
)

(e
i
)

_
f
eq
i
+
_
1
1
2
_
f
(1)
i
_
, (20)
and (e
i
)

is the component of the velocity vector e

i
in the -coordinate direction.
To specify the detailed form of

, the lattice structure and the correspond-

ing equilibrium distribution have to be specied. For simplicity and without
loss of generality, we consider here the two-dimensional square lattice with
nine velocities: e
i
= (cos(/2(i 1), sin(/2(i 1)) for i = 1, 3, 5, 7,
e
i
=

2(cos(/2(i 1) + /4, sin(/2(i 1) + /4) for i = 2, 4, 6, 8;

e
0
= 0 corresponds to a zero-speed velocity. The requirement for using the
nine-velocity model, instead of the simpler ve-velocity square lattice, comes
from the consideration of lattice symmetry: the LBE cannot recover the cor-
rect Navier-Stokes equations unless sufcient lattice symmetry is guaranteed
(Frisch et al 1986).
Note that the Navier-Stokes equations has a second-order nonlinearity. The
general formof the equilibriumdistribution function can be written up to O(u
2
)
(Chen et al 1992):
f
eq
i
=
_
a +be
i
u +c(e
i
u)
2
+du
2

, (21)
where a, b, c and d are lattice constants. This expansion is valid only for small
velocities, or small Mach numbers u/C
s
, where C
s
is the sound speed. Using
the constraints in Equation 7, the coefcients in Equation 21 can be obtained
analytically (Qian et al 1992):
f
eq
i
= w
i
_
1 +3e
i
u +
9
2
(e
i
u)
2

3
2
u
2
_
, (22)
with w
0
= 4/9, w
1
= w
3
= w
5
= w
7
= 1/9, and w
2
= w
4
= w
6
= w
8
=
1/36. Inserting the above formula into Equation 20 we have,

(0)

i
(e
i
)

(e
i
)

f
eq
i
= p

+u

(1)

=
_
1
1
2
_

i
(e
i
)

(e
i
)

f
1
i
= (

(u

) +

(u

)),
(23)
where p = /3 is the pressure, which gives a constant sound speed, C
s
=
1/

3, and = (2 1)/6 is the kinematic viscosity.

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
336 CHEN & DOOLEN
The resulting momentum equation is

_
u

t
+

_
=

p +

), (24)
which is exactly the same as the Navier-Stokes equations if the density variation
is small enough (Qian & Orszag 1993).
LBE: Approximation to the Continuum Boltzmann Equation
Although the LBE evolved from its Boolean counterpart, the LG automaton
method, it has been shown recently by two groups independently (He & Luo
1997, Abe 1997) that the LBE can be obtained from the continuum Boltzmann
equation for discrete velocities by using a small Mach number expansion. In
these derivations, the starting point is the Boltzmann BGKequation (Bhatnagar
et al 1954, Chu 1965, Xu & Prendergast 1994):
g
t
+ g =
1

(g g
eq
), (25)
where g g(x, , t ) is the single-particle distribution function in continuum
phase space (x, ), and g
eq
is the Maxwell-Boltzmann equilibriumdistribution
function:
g
eq

(2/3)
D/2
exp
_

3
2
( u)
2
_
. (26)
Here D is the spatial dimension, and for simplicity, the particle velocity and
the uid velocity u have been normalized by

3RT, giving a sound speed of

c
s
= 1/

3. T is the temperature. The macroscopic uid variables are the

velocity moments of the distribution function g:
=
_
gd, u =
_
gd, =
1
2
_
( u)
2
gd, (27)
where =
D
2
T is the internal energy. Assuming that the uid velocity in Equa-
tion 26 is a small parameter (compared with the sound speed), the equilibrium
distribution g
eq
up to O(u
2
) has the following form (Koelman 1991):
g
eq
=

(2/3)
D/2
exp
_

3
2

2
__
1 +3( u) +
9
2
( u)
2

3
2
u
2
_
. (28)
For discrete velocity models, only a small set of particle velocities e
i
=
i
(i = 1, . . . , M), and their distribution functions at these velocities, g
i
(x, t ) = g
(x, e
i
, t ), are used; and the kinetic evolution in Equation 25 will only require
the solution of g
i
. The rst two denitions in Equation 27 can be approximated
using the discrete velocities in a Gaussian-type quadrature:
(x, t ) =

i
W
i
g
i
(x, t ), u(x, t ) =

i
W
i
e
i
g
i
(x, t ).
(29)

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 337
This denes an effective distribution function f
i
(x, t ):
f
i
(x, t ) = W
i
g
i
(x, t ), (30)
which satises the simple conservation relations in Equation 3.
Given e
i
, W
i
is constant, and f
i
satises the same equation as g:
f
i
t
+e
i
f
i
=
1

_
f
eq
i
f
i
_
, (31)
with f
eq
i
= w
i
[1 +3(e
i
u) +
9
2
(e
i
u)
2

3
2
u
2
] and w
i
= W
i
/(2/3)
D/2
exp(
3
2
e
2
i
).
To obtain the weight w
i
, He & Luo (1997) use a third-order Hermite formula
to approximate the integrals in Equation 27. Abe (1997) assumes w
i
has a
simple truncated functional form based on e
i
. For the nine-velocity square
lattice, both papers found that w
0
= 4/9, w
i
(i = 1, 3, 5, 7) = 1/9 and w
i
(i =
2, 4, 6, 8) = 1/36, which give the same equilibrium distribution function as
the original lattice Boltzmann model in Equation 22.
If inEquation31the time derivative is replacedbya rst order time difference,
and the rst order upwind discretization for the convective terme
i
f
i
is used
and a downwind collision term(xe
i
, t ) for (x, t ) is used (Cao et al 1997),
then we obtain the nite difference equation for f
i
:
f
i
(x, t +t ) = f
i
(x, t )
_
f
i
(x, t ) f
i
(x xe
i
, t )

_
f
i
(x xe
i
, t ) f
eq
i
(x xe
i
, t )

, (32)
where = t |e
i
|/x, = t /, and t and x are the time step and the
grid step, respectively. Choosing = 1 and = 1, Equation 32 becomes the
standard LBE as shown in Equation 2.
From the discretization process, we note Equation 32 only has rst order
convergence in space and time to Equation 31. However, as shown by Sterling
and Chen (1995), because Equation 32 has a Lagrangian nature in its spatial
discretization, the discretization error has a special form which can be included
in the viscous term, resulting in second-order accuracy both in space and time.
Boundary Conditions in the LBM
Wall boundary conditions in the LBM were originally taken from the LG
method. For example, a particle distribution function bounce-back scheme
(Wolfram 1986, Lavall ee et al 1991) was used at walls to obtain no-slip ve-
locity conditions. By the so-called bounce-back scheme, we mean that when
a particle distribution streams to a wall node, the particle distribution scatters
back to the node it came from. The easy implementation of this no-slip veloc-
ity condition by the bounce-back boundary scheme supports the idea that the

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
338 CHEN & DOOLEN
LBMis ideal for simulating uid ows in complicated geometries, such as ow
through porous media.
For a node near a boundary, some of its neighboring nodes lie outside the
ow domain. Therefore the distribution functions at these no-slip nodes are
not uniquely dened. The bounce-back scheme is a simple way to x these
unknown distributions on the wall node. On the other hand, it was found
that the bounce-back condition is only rst-order in numerical accuracy at the
boundaries (Cornubert et al 1991, Ziegler 1993, Ginzbourg & Adler 1994).
This degrades the LBM, because numerical accuracy of the LBM in Equation 2
for the interior mesh points is second-order. He et al (1997) conrmed this
result by analyzing the slip velocity near the wall node for Poiseuille ow.
To improve the numerical accuracy of the LBM, other boundary treatments
have been proposed. Skordos (1993) suggested including velocity gradients
in the equilibrium distribution function at the wall nodes. Noble et al (1995)
proposed using hydrodynamic boundary conditions on no-slip walls by enforc-
ing a pressure constraint. Inamuro et al (1995) recognized that a slip velocity
near wall nodes could be induced by the bounce-back scheme and proposed
to use a counter slip velocity to cancel that effect. Maier et al (1996) modi-
ed the bounce-back condition to nullify net momentum tangent to the wall
and to preserve momentum normal to the wall. Zou and He (1997) extended
the bounce-back condition for the nonequilibrium portion of the distribution.
Ziegler (1993) noticed that if the boundary was shifted into uid by one half
mesh unit, i.e. placing the nonslip condition between nodes, then the bounce-
back scheme will give second-order accuracy. Simulations demonstrated that
these heuristic models yield good results for uid ows around simple wall
boundaries. The above boundary treatments were also analyzed by Zou et al
(1995) and He et al (1997) by solving the lattice BGK equation (Equation 2).
It appears, however, that the extension of these simple assumptions to arbi-
trary boundary conditions is difcult. Chen et al (1996) viewed the LBM as
a special nite difference scheme of the kinetic equation (Equation 31). They
adopted staggered mesh discretization from traditional nite difference meth-
ods and proposed using a second-order extrapolation scheme of distributions
in the ow to obtain the unknown particle distribution functions. Their extrap-
olation scheme is simple and can be extended to include velocity, temperature,
and pressure (Maier et al 1996, Zou & He 1997) boundary conditions and their
derivatives. In this treatment, boundary conditions can be assigned to grid nodes
or to any locations in space for little computational effort. Numerical simula-
tions including time-dependent Couette ow, a lid-driven square-cavity ow,
and ow over a column of cylinders were carried out, showing good agreement
with analytical solutions and nite difference solutions.

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 339
Numerical Efciency, Stability, and Nonuniform Grid
The discrete velocity equation (Equation 31) is a hyperbolic equation that
approximates the Navier-Stokes equation in the nearly-incompressible limit
(Majda 1984). The numerical accuracy depends on Mach number (Reider &
Sterling 1995). The advective velocity e
i
in Equation 31 is a constant vec-
tor, in contrast to the spatial dependent velocity in compressible hydrodynamic
equations, which prevents the LBMfromsolving nonlinear Riemann problems.
Froma numerical analysis point of view, the LBM, like other kinetic equations,
is a relaxation method (Cao et al 1997) for the macroscopic equations, which
has much in common with the explicit penalty or pseudocompressibility
method (Sterling & Chen 1996). This view was used by Ancona (1994) to
generalize the LBM to include fully Lagrangian methods for the solution of
partial differential equations. The advective velocity e
i
in Equation 31 is con-
stant in contrast to the spatial dependent velocity in compressible hydrodynamic
equations. The kinetic-relaxation method for solving a hyperbolic conservation
system was proposed by Jin and Xin (1995). This approach uses the relaxation
approach to model the nonlinear ux terms in the macroscopic equations and
thus it does not require nonlinear Riemann solvers either. Using this relax-
ation method, Jin and Katsoulakis (1997) calculated curvature-dependent front
propagation. In principle, the LBM and the kinetic-relaxation method are very
much alike. The kinetic-relaxation was developed mainly for shock capture in
Euler systems, whereas the LBM is more focused on viscous complex ows in
the nearly-incompressible limit. Nadiga and Pullin (1994) proposed a simula-
tion scheme for kinetic discrete-velocity gases based on local thermodynamic
equilibrium. Their method seems more general and able to obtain high or-
der numerical accuracy. Their nite volume technique was further developed
(Nadiga 1995) to solve the compressible Euler equation by allowing the discrete
velocities to adapt to the local hydrodynamic state. Elton et al (1995) studied
issues of convergence, consistency, stability, and numerical efciency for lat-
tice Boltzmann models for viscous Burgers equation and advection-diffusion
system.
The numerical efciency of the LBM was studied by several authors. Succi
et al (1991) noted that the number N
f po
of oating point calculations for the
LBM would be r N
D
, while N
f po
(25 log
2
N)N
D
for the pseudospectral
method, where N is the number of lattice points along each direction in D
dimensions. r is 150 for the two-dimensional LBM using the full matrix M
i j
in Equation 10 and 40 for the LBGK model. For 3-D low Reynolds number
simulations, an LBM performance of 2.5 times faster than the pseudo-spectral
method was reported (Chen et al 1992). The scalability of the LBMwas studied
by Noble et al (1996). They found a better-than-linear scalability for the LBM

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
340 CHEN & DOOLEN
when the number of processors increases while simultaneously increasing size
of the computational domain. For a 2-Dproblemwith 2176
2
lattice points, they
achieved an overall 13.9 gigaops using 512 processors of CM-5 machines.
In the continuum Boltzmann equation, the equilibrium distribution function
given by Equation 26 corresponds to the maximum entropy state. Thus, any
initial state will evolve towards a state of higher entropy. This result is known
as Boltzmanns H-theorem, which ensures an increase of entropy and ensures
stability. An H-theorem has been derived also for the LG method (Frisch et al
1987). In the LBM, however, only a small number of discrete velocities is used
and usually the equilibrium distribution is truncated to O(u
2
) (see Equation
21). Consequently, one cannot simultaneously guarantee an H-theorem and
allow the correct form of the macroscopic equations. Therefore, the LBM,
although of a kinetic nature, is subject to numerical instability. Furthermore, we
know also that the traditional LBM equation (Equation 2) is a nite-difference
solution of the discrete-velocity Boltzmann equation in Equation 31 and the
discretization in Equation 32 with = 1 marginally satises the Courant-
Friedrichs-Lewy condition. A von Neumann linearized stability analysis of the
LBMwas carried out by Sterling &Chen (1996). In this analysis, they expanded
f
i
as f
i
(x, t ) = f
(0)
i
+ f

i
(x, t ), where the global equilibrium populations f
(0)
i
were assumed to be constants that depend only on a constant mean density and
constant velocity. Taylor-expanding Equation 2, they arrived at the following
linearized equation:
f

i
(x +e
i
t, t +t ) = J
i j
f

j
(x, t ), (33)
where J
i j
is the Jacobian matrix corresponding to the coefcient of the linear
termin Equation 2 and does not depend on location or time. Spatial dependence
of the stability analysis was carried out by taking the Fourier transformation of
Equation 33 and solving the eigenvalue problem for diag[exp(i k e
i
t )]J
i j
.
For k = 0, they found the linear stability condition: 1/2 which is consistent
with the positivity of viscosity. Detailed stability analysis were carried out
for the hexagonal 2-D LBM model (Chen et al 1992), the square 2-D LBM
model and the 3-D 15-velocity model (Qian et al 1992). They concluded
that a stable LBM requires that the mean ow velocity be below a maximum
velocity that is a function of several parameters, including the sound speed C
s
,
the relaxation time , and the wave number. The numerical stability can also
be improved if one starts from the continuum kinetic equation (Equation 31)
(Ancona 1994, McNamara et al 1995, Cao et al 1997) using different nite
difference discretizations.
Until recently, one of the limitations to the numerical efciency in the LBM,
as compared with other CFD methods, is that the discretization of Equation 31
was constrained on a special class of uniform and regular lattices. To increase

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 341
numerical efciency and accuracy, nonuniform grid methods have been devel-
oped. Koelman (1991) suggested including spatial nonuniformity in the lattice
structure. Nannelli & Succi (1992) and Amati et al (1997) developed the nite
volume LBM (FVLBM) based on Equation 31. They dened a nonuniform
coarse lattice whose cell typically contained several original lattice units. The
evolution equation for the mean value F
i
V
1
C
_
C
f
i
d
3
x in the Cell C re-
quires the evaluation of the ux across the boundaries of C, where V
C
is the
volume. A piece-wise constant or a piece-wise linear interpolation was used to
approximate the ux. This FVLBM has been used to study two-dimensional
owpast a bluff body (Succi &Nannelli 1994) and three-dimensional turbulent
channel ow(Amati et al 1997). The results of the 3-Dchannel owsimulation
reveal that the FVLBM is about one order of magnitude faster than the tradi-
tional nonuniform CFD methods. On the other hand, while most simulation
results agree well with other traditional methods, the comparison of simulation
results is correspondingly less satisfactory, owing partially to the low-order
interpolation schemes. Cao et al (1997) developed a nonuniform nite differ-
ence LBM for Equation 31 and simulated annular ows in polar coordinates.
He et al (1996) adhered to the LBM equation (Equation 2) and developed a
time-dependent interpolation scheme that increases grid density in high-shear
regions. They successfully simulated ows in the 2-D symmetric channel
with sudden expansion and ow around a circular cylinder (He & Doolen
1997).
LATTICE BOLTZMANN SIMULATION
OF FLUID FLOWS
Flows with Simple Boundaries
DRIVEN CAVITY FLOWS For two-dimensional cavity ows, there are many pa-
pers using traditional schemes, including nite-difference (Schreiber & Keller
1983, E & Liu 1996) and multigrid (Vanka 1986). The fundamental charac-
teristics of the 2-D cavity ow are the emergence of a large primary vortex
in the center and two secondary vortices in the lower corners. The values
of the stream function and the locations of the centers of these vortices as a
function of Reynolds numbers have been well studied. The lattice Boltzmann
simulation of the 2-D driven cavity by Hou et al (1995) covered a wide range
of Reynolds numbers from 10 to 10,000 using a 256
2
lattice. They carefully
compared simulation results of the stream function and the locations of the
vortex centers with previous numerical simulations and demonstrated that the
differences of the values of the stream function and the locations of the vortices
between the LBM and other methods were less than 1%. This difference is

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
342 CHEN & DOOLEN
within the numerical uncertainty of the solutions using other numerical meth-
ods. The spatial distributions of the velocity, pressure and vorticity elds were
also carefully studied. An error analysis of the compressibility effect from the
model was carried out. Two-dimensional cavity owwas also studied by Miller
(1995). Instead of using a uniform velocity on the top of the cavity, in Millers
work a shear-ow condition with certain velocity distribution was given which
allowed an analytical solution of the velocity for the whole cavity. Excellent
agreement between the LBM simulations and the analytical solutions for the
velocity and pressure elds was found. The effects of the variation of the cavity
aspect ratio on vortex dynamics were also studied.
Hou (1995) simulated three-dimensional cubic cavity ow using the three-
dimensional 15-velocity LBM (Qian et al 1992, Chen et al 1992). 128
3
lattice
points were used with Re = 3,200. Flows at this Reynolds number had been
extensively studied earlier, including the nite-difference simulation and the
experimental work by Prasad & Koseff (1989). Flow structures, including
the velocity and vorticity elds in different planes were analyzed using the
LBM simulation. They compared well with previous numerical studies. The
correlation of the Taylor-G ortler-like vortices in the yz- and xz-planes was
reported for the rst time. Figure 1 displays the mean velocity proles from
3-D LBM, 2-D LBM simulations, and experimental work (Prasad & Koseff
1989) in the symmetry plane along the vertical and the horizontal centerlines
(left), and the root-mean-square (rms) velocity proles (right). The agreement
Figure 1 Mean velocity proles (left),
u
U
and
v
U
, in the symmetry plane along the vertical and
the horizontal centerlines. Solid line is for the 3-D LBE simulation, dotted line is a 2-D simulation
(Hou 1995), and circles represent experimental results (Prasad & Koseff 1989). rms velocities
(right), U
rms
and V
rms
, in the symmetry plane along vertical and horizontal centerlines. Solid line
is the LBE simulation, upward triangles and downward triangles are experimental results (Prasad
& Koseff 1989).

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 343
between the 3-D LBM results and experimental results shows that the LBM is
capable of simulating complex 3-D unsteady ows.
FLOW OVER A BACKWARD-FACING STEP The two-dimensional symmetric sud-
den expansion channel ow was studied by Luo (1997) using the LBM. The
main interest in Luos research was to study the symmetry-breaking bifurcation
of the ow when Reynolds number increases. In this simulation, an asymmet-
ric initial perturbation was introduced and two different expansion boundaries,
square and sinusoidal, were used. This simulation reproduced the symmetric-
breaking bifurcation for the ow observed previously, and obtained the critical
Reynolds number of 46.19. This critical Reynolds number was compared with
earlier simulation and experimental results of 40.45 and 47.3, respectively.
Qian et al (1996) used the ow over a backward-facing step as a benchmark
for validating the LBM. They studied the recirculation length as a function of
Reynolds number and the channel expansion ratio. The results from the LBM
compared well with previous results using nite-difference methods, spectral-
element methods, and experiments. The appearance of the second vortex and
the transition leading to turbulence were also studied.
FLOW AROUND A CIRCULAR CYLINDER The ow around a two-dimensional
circular cylinder was simulated using the LBM by several groups of people.
The ow around an octagonal cylinder was also studied (Noble et al 1996).
Higuera &Succi (1989) studied owpatterns for Reynolds number up to 80. At
Re = 52.8, theyfoundthat the owbecame periodic after a longinitial transient.
For Re = 77.8, a periodic shedding ow emerged. They compared Strouhal
number, ow-separation angle, and lift and drag coefcients with previous
experimental and simulation results, showing reasonable agreement. Pressure
distributions aroundthe surface of the cylinder as a functionof Reynolds number
were carefully studied by Wagner (1994). He concluded that the difference
between his Strouhal numbers and those from other simulations was less than
3.5%. Considering the compressible nature of the LBM scheme, this result is
signicant. It demonstrates that in the nearly incompressible limit, the LBM
simulates the pressure distribution for incompressible uids well (Martinez
et al 1994). The above two simulations were based on uniform lattices, which
poorly approximate circular geometry. The no-slip boundary was enforced at
the boundary mesh point which is jagged and not necessary on the cylinder.
He & Doolen (1997) revisited the problem of two-dimensional ow around a
circular cylinder basedonthe interpolation-supplement LBM(He et al 1996). In
this simulation, the underlying lattice was square, but a spatial interpolation was
used. A polar coordinate was constructed and the simulation domain consisted
of a circular region. The no-slip wall boundary condition was exactly enforced

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
344 CHEN & DOOLEN
on the cylinder boundary. It is clear that the computational accuracy in this
simulation is better than previous LBM simulations. The streaklines showed
vortex shedding quite similar to experimental observation. The LBMcalculated
Strouhal numbers. Lift and drag coefcients for ows at different Reynolds
numbers compared well with previous simulations using other schemes. The
error was within experimental uncertainty and is the same as errors among other
schemes.
Flows in Complex Geometries
An attractive feature of the LBMis that the no-slip bounce-back LBMboundary
condition costs little in computational time. This makes the LBMveryuseful for
simulating ows in complicated geometries, such as owthrough porous media,
where wall boundaries are extremely complicated and an efcient scheme for
handing wall-uid interaction is essential.
The fundamental problem in the simulation of uid ows through porous
media is to understand and model the macroscopic transport from microscopic
uid dynamics. Astarting point for this problemis Darcys law, which assumes
a linear relation between the pressure gradient, P/L, and the volume ow
rate per unit area q: q = K/(
0
)P/L. Here K is the permeability.
Previous numerical simulations, including nite-difference schemes (Schwarz
et al 1994) and networking models (Koplik & Lasseter), were either limited
to simple physics, small geometry size, or both. Lattice gas automata were
also used for simulating porous ows and verifying Darcys law in simple and
complicated geometries (Rothman 1988). Succi et al (1989) utilized the LBMto
measure the permeability in a 3-D random media. Darcys law was conrmed.
Cancelliere et al (1990) studied the permeability K as a function of solid fraction
in a system of randomly positioned spheres of equal radii. The simulation
covered a wide range of (0.02 0.98). The values of K from the LBM
compared well with the Brinkman approximation: K = K
0
[1 + 3/4(1

8/ 3)] for < 0.2, and agreed well with the semiempirical Kozeny-
Carman equation: K = (1 )
3
/(6s
2
). Here K
0
is the effective permeability
for a single sphere and s is the specic surface area for the system. Later, Heijs
& Lowe (1995) conrmed this result independently. In addition, they studied
the validity of the Kozeny-Carman equation for soil samples where ow occurs
only through some specic continuous connected pore and ows occurring at
smaller scales are negligible. For this condition, the Kozeny-Carman equation
provided a much less successful estimate of the permeability. Flows through
sandstones measured using X-ray microtomography were simulated by Buckles
et al (1994), Soll et al (1994), and Ferr eol & Rothman (1995). They found that
the permeability for these sandstones, although showing large variation in space
and ow directions, in general agreed well with experimental measurements

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 345
within experimental uncertainty. Ferr eol & Rothman (1995) also studied the
effect of grid resolution on permeability. They found that the permeability
strongly depends on the viscosity when the minimum channel in the pores is
several grid points wide.
Spaid & Phelan (1997) investigated the injection process in resin transfer
molding. For this heterogeneous porous media simulation, they used the LBM
simulationof full Navier-Stokes ows tomodel a owarounda circular cylinder.
Inside the cylinder, the velocity is governed by the Brinkman equation (
2
u
u = P/, where is a model parameter). Excellent agreement between the
LBM simulation and lubrication theory for cell permeabilities was reported.
Simulation of Fluid Turbulence
DIRECT NUMERICAL SIMULATION A major difference between the LBM and
the LG method is that the LBM can be used for smaller viscosities. Conse-
quently the LBM can be used for direct numerical simulation (DNS) of high
Reynolds number uid ows. To validate the LBM for simulating turbulent
ows, Martinez et al (1994) studied decaying turbulence of a shear layer using
both the pseudospectral method and the LBM. The initial shear layer consisted
of uniformvelocity reversing sign in a very narrowregion. The initial Reynolds
number was 10,000. The simulation used a 512
2
mesh and ran for 80 large-
eddy turnover times. Martinez et al carefully compared the spatial distribution,
time evolution of the stream functions, and the vorticity elds. Energy spectra
as a function of time, small scale quantities, including enstrophy, palinstrophy,
and 4th-order enstrophy as a function of time were also studied. The correla-
tion between vorticity and stream function was calculated and compared with
theoretical predictions. They concluded that the LBE method provided a so-
lution that was accurate in the sense that time histories of global quantities,
wavenumber spectra, and vorticity contour plots were very similar to those
obtained from the spectral method. In particular, they reproduced details of
the wavenumber spectra at high wavenumbers as well as the detailed structure
of vortex distributions. Their conclusions are signicant for simulations of
turbulence in which small-scale dynamics are important. In that paper they
found discrepancies in the spatial locations of vortical structures at very late
stages.
Figure 2 (top left) compares the energy spectra froma pseudospectral method
(dotted line) and the LBM(solid line) for t = 5, and displays the time evolution
of the enstrophy function (bottom left). Figure 2 also depicts contour plots of
the vorticity function at the same time step (right). The agreement of these
two methods is excellent. Two-dimensional forced isotropic turbulence was
simulated by Benzi & Succi (1990) to study the enstrophy cascade range. Two-
dimensional forced turbulence was also simulated by Qian et al (1995) to study

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
346 CHEN & DOOLEN
Figure 2 Velocity energy spectra versus wavenumber, k (top left). The enstrophy of the system,

2
(bottom left), as a function of time. Right, a comparison of the vorticity distributions from a
pseudospectral method and the LBM (Martinez et al 1994).
the energy inverse cascade range. They reproduced the k
5/3
inertial range
scaling, in good agreement with theoretical predictions (Kraichnan 1967).
Chen et al (1992) made a similar validation of LBM by comparing results
of three-dimensional isotropic turbulence using the LBM and the pseudospec-
tral method. In that work, they simulated three-dimensional Beltrami ows,
the decaying Taylor-Green vortex, and decaying three-dimensional turbulence.
The agreement between the two methods for spatial and time distributions of
velocities and vortices was good. This assessment was echoed by Trevi no &
Higuera (1994), who studied the nonlinear stability of Kolmogorov ows using
the pseudospectral method and the LBM at various Reynolds numbers.
The LBM is a useful direct numerical simulation (DNS) tool for simulating
nonhomogeneous turbulent ows. To validate the generalized extended self-
similarity for anisotropic ows where the simple extended self-similarity is not

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 347
valid, Benzi et al (1996) used the LBMto simulate three-dimensional shear ow
where a sinusoidal force in the x direction varying along the z direction was
applied. This simulation generated useful velocity and vorticity information,
and the scaling exponents of the velocity increment and their dependence on
the shear rate were studied. Based on their analysis of this ow, they concluded
that the generalized extended self-similarity was valid for anisotropic ows.
Succi et al (1991) studied the bifurcation of a two-dimensional Poiseuille ow
and obtained the amplitude of the primary and secondary bifurcated models as
a function of the Reynolds number. Eggels (1996) utilized the LBM on the
FCHC lattice (dHumi` eres et al 1986) as a DNS method for simulating three-
dimensional channel ows with Re

= 180 (based on the friction velocity).

Their simulation results, including mean stream-wise velocity proles and rms
velocities as functions of the distance from the wall, were compared with those
reported by Kim et al (1987), showing good agreement.
LBMMODELS FORTURBULENT FLOWS As in other numerical methods for solv-
ing the Navier-Stokes equations, a subgrid-scale (SGS) model is required in the
LBMto simulate ows at very high Reynolds numbers. Direct numerical simu-
lation is impractical due to the time and memory constraints required to resolve
the smallest scales (Orszag & Yakhot 1986). Hou et al (1996) directly ap-
plied the subgrid idea in the Smagorinsky model to the LBM by ltering the
particle distribution function and its equation in Equation 2 using a standard
box lter. To account for the nonlinear term in the equilibrium distribution
(Equation 22), Hou et al (1996) replaced the relaxation time in Equation 13
by (3
0
+ C
2
s

2
|S|
2
) + 1/2, where
0
= (2 1)/6 is the kinetic viscosity,
is the lter width, |S| is the amplitude of the ltered large-scale strain-rate
tensor, and C
s
is the Smagorinsky constant. This simple Smagorinsky-type
LBM subgrid model was used to simulate ows in a two-dimensional cavity
at Reynolds numbers up to 10
6
using a 256
2
mesh (Hou et al 1996). Three-
dimensional turbulent pipe ow was simulated by Somers (1993) using the
LBM SGS model at a Reynolds number of 50,000 for mesh sizes up to 80
3
.
Friction as a function of Reynolds number compared reasonably well with the
Blasius power law for turbulent ows. Succi et al (1995) suggested using the
relaxation idea in the LBM to solve the turbulent K equations. Hayot &
Wagner (1996) used a dependent relaxation (t ) in the lattice BGK equation
(Equation 13), which leads to an effective turbulent viscosity. Eggels (1996)
and Eggels & Somers (1995) included the turbulent stress tensor
t
directly
into the equilibrium distribution, f
eq
i
in Equation 22, by adding an extra term:
f

i
= [e
i
e
i

t

M
D
tr(
t
)]. Using this model, Eggels (1996) carried out
a large-eddy simulation of the turbulent ow in a bafed stirred tank reactor.
The impact of a mechanical impeller was modeled via a spatial and temporal

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
348 CHEN & DOOLEN
dependent force on the momentum equation. It was reported that the subgrid
LBE turbulent model coupled with the impeller model produced satisfactory
results. Both mean ow and turbulent intensities compared well with exper-
imental data. This simulation demonstrated the potential of the LBM SGS
model as a useful tool for investigating turbulent ows in industrial applica-
tions of practical importance.
LBM SIMULATIONS OF MULTIPHASE
AND MULTICOMPONENT FLOWS
The numerical simulation of multiphase and multicomponent uid ows is
an interesting and challenging problem because of difculties in modeling
interface dynamics and the importance of related engineering applications, in-
cluding ow through porous media, boiling dynamics, and dendrite formation.
Traditional numerical schemes have been successfully used for simple inter-
facial boundaries (Glimm et al 1981, Brackbill et al 1992, Chang et al 1996).
The LBM provides an alternative for simulating complicated multiphase and
multicomponent uid ows, in particular for three-dimensional ows.
Method of Gunstensen et al
Gunstensen et al (1991) were the rst to develop the multicomponent LBM
method. It was based on the two-component LG model proposed by Rothman
& Keller (1988). Later, Grunau et al (1993) extended this model to allow
variations of density and viscosity. In these models, red and blue particle
distribution functions f
(r)
i
(x, t ) and f
(b)
i
(x, t ) were introduced to represent
two different uids. The total particle distribution function (or the color-blind
particle distribution function) is dened as: f
i
= f
(r)
i
+ f
(b)
i
. The LBM
equation is written for each phase:
f
k
i
(x +e
i
, t +1) = f
k
i
(x, t ) +
k
i
(x, t ), (34)
were k denotes either the red or blue uid, and

k
i
=
_

k
i
_
1
+
_

k
i
_
2
(35)
is the collision operator. The rst term in the collision operator, (
k
i
)
1
, rep-
resents the process of relaxation to the local equilibrium similar to the LBGK
model in Equation 13:
_

k
i
_
1
=
1

k
_
f
k
i
f
k(eq)
i
_
. (36)
Here, f
k(eq)
i
is the local equilibrium distribution depending on the local macro-
scopic variables
(k)
and u
(k)
.
k
is the characteristic relaxation time for species

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 349
k. The viscosity of each uid can be selected by choosing the desired
k
. Con-
servation of mass for each phase and total momentumconservation are enforced
at each node during the collision process:

r
=

i
f
r
i
=

i
f
r(eq)
i
,
b
=

i
f
b
i
=

i
f
b(eq)
i
,
u =

i,k
f
k
i
e
i
=

i,k
f
k(eq)
i
e
i
,
(37)
where =
r
+
b
is the total density and u is the local total momentum.
The form of f
k(eq)
i
can be chosen to be similar to Equation 22.
The additional collision operator (
k
)
2
contributes to the dynamics in the
interfaces and generates a surface tension:
_

k
i
_
2
=
A
k
2
|F|
_
(e
i
F)
2
/|F|
2
1/2
_
, (38)
where F is the local color gradient, dened as:
F(x) =

i
e
i
(
r
(x +e
i
)
b
(x +e
i
)). (39)
Note that in a single-phase region of the incompressible uid model, Fvanishes.
Therefore, the second term of the collision operator (
k
i
)
2
only contributes to
interfaces and mixing regions. The parameter A
k
is a free parameter, which
determines the surface tension. The additional collision term in Equation 38
does not cause the phase segregation. To maintain interfaces or to separate
the different phases, the LBM by Gunstensen et al follows the LG method of
Rothman & Keller (1988) to force the local color momentum, j =

i
( f
r
i

f
b
i
)e
i
, to align with the direction of the local color gradient after collision. In
other words, the colored distribution functions at interfaces were redistributed
to maximize j F. Intuitively, this step will force colored uids to move
toward uids with the same colors.
The multicomponent LBM by Gunstensen et al (1991) has two drawbacks.
First, the procedure of redistribution of the colored density at each node requires
time-consuming calculations of local maxima. Second, the perturbation step
with the redistribution of colored distribution functions causes an anisotropic
surface tension that induces unphysical vortices near interfaces. Dortona et al
(1994) modied the model of Gunstensen et al. In their model, the recoloring
step is replaced by an evolution equation for f
k
i
that increases computational
efciency.
Method of Shan & Chen
Shan & Chen (1993) and Shan & Doolen (1995) used microscopic interactions
to modify the surface-tensionrelated collision operator for which the surface

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
350 CHEN & DOOLEN
interface can be maintained automatically. In these models, (
k
i
)
2
in Equation
38 was replaced by the following formula:
_

k
i
_
2
= e
i
F
k
. (40)
Here F
k
is an effective force on the kth phase owing to a pairwise interaction
between the different phases:
F
k
(x) =

i
V
kk
(x, x +e
i
)e
i
. (41)
Here V
kk
is an interaction pseudopotential between different phases (or com-
ponents):
V
kk
(x, x

) = G
kk
(x, x

)
k
(x)
k

(x

). (42)
G
kk
(x) is the strength of the interaction; and
k
(x) is a function of density
for the k phase at x which has taken the following empirical form: =

0
[1 exp(/
0
)], where
0
is a constant free parameter. It was shown that
this form of the effective density gives a non-ideal equation of state, which
separates phases or two component uids. Qian et al (1995) have recently
demonstrated that other pseudopotential forms, including fractional potential
and a van der Waals potential, will produce similar results.
It should be noted that the collision operator (
k
i
)
2
in the Shan-Chen model
does not satisfy local momentum conservation. This treatment is physically
plausible because the distant pairwise interactions between phases change the
point-wise local momenta at the positions involved in the interactions. This
feature differs from the model of Gunstensen et al where the total local mo-
mentum is conserved (see the third operation in Equation 37). It is arguable
that this spurious conservation might be one reason why the model exhibits
unphysical features near interfaces.
For simplicity, only the nearest-neighbor interactions were involved in the
Shan-Chen model, in which G
kk
(x, x

) is constant when x x

= e
i
, and zero
otherwise. G
kk
acts like a temperature; when G is smaller than the critical value
G
c
(depending on the lattice structure and initial density), the uids separate.
Theoretical calculations indicate (Shan & Chen 1994) that the surface tension
MG/[2D(D +1)], which was also veried by numerical simulation.
In the Shan-Chen model, the separation of uid phases or components is
automatic (Chen 1993). This is an important improvement in numerical ef-
ciency compared with the original LBM multiphase models. The Shan-Chen
model also improves the isotropy of the surface tension.
Free Energy Approach
The above multiphase and multicomponent lattice Boltzmann models are based
on phenomenological models of interface dynamics and are probably most

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 351
suitable for isothermal multicomponent ows. One important improvement
in models using the free-energy approach (Swift et al 1995, 1996) is that the
equilibrium distribution can be dened consistently based on thermodynamics.
Consequently, the conservation of the total energy, including the surface energy,
kinetic energy, and internal energy can be properly satised (Nadiga & Zaleski
1996).
The van der Waals formulation of quasilocal thermodynamics for a two-
component uid in thermodynamic equilibrium at a xed temperature has the
following free-energy functional:
(r) =
_
dr[(T, ) + W()]. (43)
The rst term in the integral is the bulk free-energy density, which depends
on the equation of state. The second term is the free-energy contribution from
density gradients and is related to the surface tension. For simple multiphase
uid ows, W =

2
()
2
, where is related to the surface tension. The full
pressure tensor in a nonuniform uid has the following form:
P

= p

. (44)
The pressure is obtained from the free energy: p(r) =

(r) = p
0

2


2
()
2
, where p
0
=

() () is the equation of state. For a

van der Waals equation of state, = TIn /(1 b) a
2
, where a and b
are free parameters.
To incorporate the above formulation into the lattice Boltzmann equation,
Swift et al (1995) added a term to the original equilibrium distribution in Equa-
tion 21: f
eq
i

= f
eq
i
+ G

e
i
e
i
. They vary the coefcients in f
eq
i
(see Equa-
tion 21) and G

to guarantee the conservation of mass and momentum and to

obtain the following stress tensor condition:

i
f
eq
i

e
i
e
i
= P

+u

.
Numerical Verication and Applications
Two fundamental numerical tests associated with interfacial phenomena have
been carried out using the multiphase and multicomponent lattice Boltzmann
models. In the rst test, the lattice Boltzmann models were used to verify
Laplaces formula by measuring the pressure difference between the inside and
the outside of a droplet: P
inner
P
outer
=

R
, where P
inner
and P
outer
denote
the pressures inside and outside of the droplet, respectively. R is the radius
of the droplet and is the surface tension. The simulated value of has
been compared with theoretical predictions, and good agreement was reported
(Gunstensen et al 1991, Shan and Chen 1993, Swift et al 1995). Alternatively,

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
352 CHEN & DOOLEN
the surface tension can also be measured using the mechanical denition
=
_

(P
N
P
T
)dz, (45)
where P
N
and P
T
are the normal and the tangential pressure along a at inter-
face, respectively, and the integration is evaluated in a direction perpendicular
to the interface. It is shown from LBM simulations (Gunstensen et al 1991,
Grunau et al 1993, Swift et al 1995) that the surface tension obtained from the
mechanical denition agrees with the Laplace formula denition.
In the second test of LBM interfacial models, the oscillation of a capillary
wave was simulated (Gunstensen et al 1991, Shan & Chen 1994, Swift et al
1995). A sine wave displacement of a given wave vector was imposed on an
interface that had reached equilibrium. The resulting dispersion relation was
measured and compared with the theoretical prediction (Laudau & Lifshitz
1959):
2
= /k
3
. Good agreement was observed, validating the LBM
surface tension models.
Because the lattice Boltzmann multiphase and multicomponent models do
not track interfaces, it is easy to simulate uid ows with very complicated
interfaces, such as the domain growth in the spinodal decomposition process for
binary uids (Alexander et al 1993a, Rybka et al 1995). In these simulations,
hydrodynamics plays an important role and multiphase or multicomponent
uids evolved from an initial mixed state owing to surface tension, viscosity,
and inertial force.
A convenient way to characterize the growth kinetics is to use the correlation
function of the order parameter, G(r, t ) (r)(0), where is the order
parameter, dened as (
r

b
)/(
r
+
b
). The Fourier transform of G(r, t )
is then the structure function S(k,t). As time evolves, the structure function
becomes more sharply peaked, and its maximum value moves toward small
k. In a wide variety of phase segregating systems, S(k,t) has a simple form at
late times: S(k, t ) = R
D
(t )F(x). Here R(t ) is the average domain size, D
is the dimension, and F(x) is a universal function depending on x = k R(t ).
Numerical calculations of R(t ) and F(x) have been carried out using the LBM
(Alexander et al 1993a). It was found that R(t ) scales t
2/3
in two dimen-
sions and t in three dimensions for two-component uids. This agrees with
theoretical predictions. F(x) was found to agree well with Porods law for
large x, x
D+1
, but small scales showed strong dependence of ow properties,
which was later conrmed by Wu et al (1995) using the Langevin equation.
Osborn et al (1995) studied similar scaling problems for liquid-gas two-phase
domain growth. They found that R(t ) t
1/2
and t
2/3
at high and low vis-
cosities, respectively. A similar crossover for two-component uids from t
1/3
to t
2/3
was also reported as time increased. The hydrodynamic effects on the

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 353
spinodal decomposition for three and four components in two and three di-
mensions in critical and off-critical quenching were also studied by Chen &
Lookman (1995). Cieplak (1995) carried out a study of two-dimensional rup-
ture dynamics under initial perturbation and compared results with molecular
dynamics simulations. The coalescence process was also studied for a droplet
in a background uid falling in a gravitational eld to the bottom of a con-
tainer in which the bottom is a bare wall, a shallow liquid with the same
kind as the droplet, or a deeper liquid. Interesting interfacial structures were
observed.
In Figure 3 we showcondensation and subsequent coalescence of liquid drops
insupersaturatedvapor at twodifferent times. The LBEmodel withinterparticle
interaction (Shan & Chen) was used in this simulation on a 128
3
lattice. The
uid obeys the equation of state: p =
1
3
+
3
2
G(1 e

)
2
. G (chosen to
be 0.55) is a parameter that gives the relative strength of the interaction.
This equation of state gives a non-monotonic isotherm with a critical point at
G =
4
9
. The initial density is constant except for a randomly distributed small
perturbation.
One feature of the free-energy approach is that one can parameterize the free
energy to simulate complicated uid ows, including lamellar uids. In their
study, Gonnella et al (1997) used the following Ginzburg-Landau potential
for the function in Equation 43 based on the order parameter, : =

2
+

4

4
. They included high-order gradients in the W function: W =

2
()
2
+

2
(
2
)
2
where and are positive. < 0 and k < 0 lead to
lamellar phases. The magnitude of is directly connected to the elasticity of
Figure 3 Condensation and subsequent coalescence of liquid drops in supersaturated vapor. Early,
left, late, right (Shan 1997).

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
354 CHEN & DOOLEN
Figure 4 Conguration of a lamellar uid when a shear ow is applied. The shear velocities are
v = 0 (a), v = 0.1 (b), and v = 0.3 (c) (Gonnella et al 1997).
the uid. In Figure 4, we present ow patterns for a lamellar uid under a
shear induced by forcing a constant velocity v in the left and v in the right
simulation domains, respectively. As the shear was applied, ordered lamellae
quickly formed, aligned at an angle to the direction of shear. A further increase
of the shear rate resulted in a decrease in the angle between the direction of
the lamellae and that of the shear. This pattern directly reects the balance
between the shear force and elastic force of the lamellar uids. It would be
of great interest to use this model for simulating bio-uids where elasticity is
important.
Lattice Boltzmann multiphase uid models have been extensively used to
simulate multicomponent owthrough porous media in order to understand the
fundamental physics associated with enhanced oil recovery, including relative
permeabilities (Vangenabeek et al 1996). The LBM is particularly useful for
this problem because of its capability of handling complex geometrical bound-
ary conditions and varying physical parameters, including viscosities (Grunau
et al 1993) and wettabilities (Blake et al 1995). In the work by Buckles et al
(1994), Martys & Chen (1995), and Ferreol & Rothman (1995), realistic sand-
stone geometries from oil elds were used. Very complicated ow patterns
were observed. The numerical values of the relative permeability as a function
of percent saturation of wetting uid agree qualitatively with experimental data.

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 355
Figure 5 A sandstone sample used in the multiphase LBM simulation (top left). Other panels
display two-phase uid ows through the sandstone at different times (Buckles et al 1994).
Gunstensen & Rothman (1993) studied the linear and nonlinear multicompo-
nent ow regimes corresponding to large and small ow rates, respectively.
They found, for the rst time, that the traditional Darcys law must be modi-
ed in the nonlinear regime because of the capillary effect. Figure 5 shows a
typical pore geometry (top left). It is a portion of a sandstone sample obtained
from a Mobil offshore oil reservoir. Two-phase uid ows through the sand-
stone are displayed at successive times (top right, bottom left, bottom right).
The sandstone is transparent. The dark color indicates invading water and the
grey color indicates oil. As seen in the plot, the lattice Boltzmann simulation
preserves the fundamental phenomena observed in experiments that the water
phase forms long ngers through the porous medium because of the wettability
properties of the water. The LBM is becoming an increasingly popular means
of modeling multiphase uid ows in porous media because of its ability to
simulate the exact Navier-Stokes equation in a parallel fashion, to handle com-
plicated geometry, and to simulate surface dynamics and wettability (or contact
angles).

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
356 CHEN & DOOLEN
LATTICE BOLTZMANN SIMULATION
OF PARTICLES IN FLUIDS
To simulate particles suspended in uids, an approximate treatment of uid-
particle interaction must be incorporated. Much of the pioneering work and
some interesting applications in this area were carried out by Ladd (1993,
1994a,b, 1997). Important variations and applications were mostly associated
with Behrend (1995) and Aidun & Lu (1995).
In the model of Ladd, a solid boundary was mapped onto the lattice and a set
of boundary nodes, r
b
, was dened in the middle of links, whose interior points
represent a particle. A no-slip boundary condition on the moving particle
requires the uid velocity to have the same speed at the boundary nodes as
the particle velocity u
b
which has the translational part U, and the rotational
part
b
. Assuming that the center position of the particle is R, then, u
b
=
U +
b
(r
b
R). The distribution function f
i
is dened for grid points
inside and outside the particle. To account for the momentum change when u
b
was not zero, Ladd proposed adding a term to the distribution function for both
sides of the boundary nodes:
f

i
(x) = f
i
(x) B(e
i
u
b
), (46)
where B is a coefcient proportional to the mass density of the uid and depends
on the detailed lattice structure; the +sign applies to boundary nodes at which
the particle is moving toward the uid and for moving away from the uid.
Equation 46 includes a mass exchange between the uid inside and the uid
outside the particle. Aidun & Lu (1995) modied Equation 46 by adding
an extra term that forces mass conservation for the uid inside the particle.
The details of the boundary rule and variations were studied extensively by
Behrend (1995), who improved the efciency of the algorithm used by Ladd
while retainingsimilar accuracyfor translational motion, althoughthe rotational
motion is less accurate.
The approximations used so far to simulate moving boundaries are compu-
tationally convenient and efcient. In fact, the work for simulating N particles
in Ladds scheme scales linearly with N, in contrast to nite-element schemes,
which scale as N
2
(Hu 1996) when the exact incompressible condition is being
enforced. The accuracy of the scheme was carefully and extensively studied for
creeping ows and ows at nite Reynolds numbers (Ladd 1994b). It compared
well with nite-difference and nite-element methods. The drag coefcient of
a circular particle in a 2-D channel with gravitational eld and the settling tra-
jectory compared well with results from traditional numerical methods and the
solution of the Stokes equation (Aidun & Lu 1996).

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 357
In the LBM, uctuations in the distribution functions are ignored. These
uctuations are important for the study of Brownian motion and other effects.
To include the uctuations in the uid necessary to drive Brownian motion,
Dufty &Ernst (1996) and Ladd (1993) added stochastic terms to the distribution
function. These distribution function uctuations are constrained to conserve
mass and momentum, but they contribute a uctuating part to the stress tensor
of the uid. The variance of the stress-tensor uctuations is related to the
effective temperature andthe viscosityof the uidvia the uctuation-dissipation
theorem. This method allowsapparently for the rst timetreatment of
the Brownian short-time regime and the pre-Brownian time regime. Close
quantitative agreement is found between experiment and the uctuating LBM
in the short-time regime (Segre et al 1995). Ladd, using up to 32,000 suspended
particles in the uctuation LBM, showed that there is no evidence that the long-
range hydrodynamic interactions are screened by changes in the pair correlation
function at large distances. His results leave open the explanation of the
experimentally reported absence of the theoretically expected divergence of
velocity uctuations (Ladd 1997).
SIMULATION OF HEAT TRANSFER
AND REACTION-DIFFUSION
The lattice Bhatnagar-Gross-Krook (LBGK) models for thermal uids have
been developed by several groups. To include a thermal variable, such as tem-
perature, Alexander et al (1993b) used a two-dimensional 13-velocity model
on the hexagonal lattice. In this work, the internal energy per unit mass E
was dened through the second-order moment of the distribution function,
E =

,i
f
,i
(e
,i
u)
2
/2. Here the index = 0, 1, and 2, indicates
velocity magnitudes. To have a consistent compressible hydrodynamic equa-
tion, the collision operator
,i
was chosen to satisfy local energy conserva-
tion:

i

,i
e
2
,i
/2 = 0, and the third-order velocity dependent terms, such
as (e
,i
u)
3
, (e
,i
u)u
2
and u
3
, had been included in the equilibrium distri-
bution (Equation 21). The transport coefcients, including viscosity and heat
conductivity derived from the Chapman-Enskog expansion, agreed well with
numerical simulation. ACouette owwith a temperature gradient between two
parallel planes was also simulated and the results agreed well with theoretical
predictions when the temperature difference was small. Vahala et al (1995)
examined this model by studying the effect of 2-D shear velocity turbulence on
a steep temperature gradient prole. Qian (1993) developed 3-D thermal LBM
models based on 21 and 25 velocities. Chen et al (1994) extended the above
thermal models by including more speeds based on general square lattices in

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
358 CHEN & DOOLEN
D dimensions. They utilized the freedom in the equilibrium distribution to
eliminate nonphysical high-order effects (Qian & Orszag 1993). Because the
multispeed BGK collision operator constrained the momentum and the thermal
modes to relax at the same speed, the Prandtl number was xed. The inclusion
of additional properties to provide variation of the Prandtl number was studied
by Chen et al (1995) and McNamara et al (1995).
Two limitations in multispeed LBM thermal models severely restrict their
application. First, because only a small set of velocities is used, the varia-
tion of temperature is small. Second, all existing LBM models suffer from
numerical instability (McNamara et al 1995), owing to the absence of an
H-theorem. It seems that the numerical instability is more severe in the thermal
LBM models than in the isothermal models. Both problems can be partially
resolved if one treats the temperature equation as an active scalar and solves
it using the relaxation method by adding an independent distribution func-
tion (Bartoloni et al 1993). On the other hand, it is difcult for the active
scalar approach to incorporate the correct and full dissipation function. Two-
dimensional Rayleigh-B enard (RB) convection was simulated using this active
scalar scheme for studying scaling laws (Bartoloni et al 1993) and probabil-
ity density functions (Massaioli et al 1993) at high Prandtl numbers. Two-
dimensional free-convective cavity ow was also simulated (Eggels & Somers
1995), and the results compared well with benchmark data. Two-D and 3-D
Rayleigh-B enard convections were carefully studied by Shan (1997) using a
passive scalar temperature equation and a Boussinesq approximation. This
scalar equation was derived based on the two-component model of Shan &
Chen (1993). The calculated critical Rayleigh number for the RB convection
agreed well with theoretical predictions. The Nusselt number as a function of
Rayleigh number for the 2-D simulation was in good agreement with previous
numerical simulation using other methods.
The LBM was extended by Dawson et al (1993) to describe a set of reaction-
diffusion equations advected by velocities governed by the Navier-Stokes equa-
tion. They studied hexagonal patterns and stripe patterns caused by the Turing
instability in the Selkov model. The effect of uid ows on the chemical re-
action at solid surfaces was simulated (Chen et al 1995) to study geochemical
processes, including dissolution and precipitation on rock surfaces and chimney
structures at seaoor hydrothermal vents. Asimilar study on the effect of nutri-
ent diffusion and ow on coral morphology was carried out by Kaandorp et al
(1996). Alvarez-Ramirez et al (1996) used the model of Dawson et al (1993) to
study the effective diffusivity of a heterogeneous medium when the inclusion
was impermeable or permeable with a different diffusivity. They found that
the LBM simulation results compared well with Monte Carlo simulations and

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 359
agreed well with the predictions using Maxwells equation. Holme &Rothman
(1992) carried out a study of viscous ngering in two-dimensional Hele-Shaw
patterns with P eclet number = 1700. Creeping ow in a Hele-Shaw cell was
also simulated by Flekky et al (1996) to investigate the inertial effect at very
small Reynolds number, which is sensitive to effects of hydrodynamic irre-
versibility. Anomalous diffusion of LBM uids in fractal media and Taylor
hydrodynamic dispersion in a two-dimensional channel were simulated by Cali
et al (1992). They found that the diffusion scaling exponents agreed well with
the exact solution for a Sierpinski gasket. The measured longitudinal diffusion
coefcient as a function of ow speed in 2-D Taylor hydrodynamic dispersion
agreed well with an asymptotic prediction. This research demonstrates that the
LBMis capable of simulating the diffusion process in irregular geometries. The
effects of turbulent mixing behind a wake on Turing patterns were studied by
Weimar & Boon (1996). They found that the mixing caused by turbulence in-
creased eddy diffusivity and destroyed the pattern formation. Qian et al (1995)
simulated front dynamics and scaling properties for the irreversible reaction:
A + B C.
CONCLUDING REMARKS
This reviewis intended to be an overviewof the subject of the lattice Boltzmann
method and its applications in uid mechanics. The LBM is so diverse and
interdisciplinary that it is not possible to include all interesting topics. We refer
readers to related review articles written by Benzi et al (1992), Rothman &
Zaleski (1994), Chen et al (1995), and Qian et al (1995).
We have introduced the fundamentals of the lattice Boltzmann method, in-
cluding the lattice Boltzmann equation and its relation to the macroscopic
Navier-Stokes equation. Some LBM methods, including multiphase LBM
models and uid-particle interaction models, have been discussed in detail.
We have demonstrated that simulation results fromthe lattice Boltzmann meth-
ods are in good quantitative agreement with experimental results and results
from other numerical methods.
We wish to emphasize the kinetic nature of the LBM. Because the LBM is
a mesoscopic and dynamic description of the physics of uids, it can model
problems wherein both macroscopic hydrodynamics and microscopic statistics
are important. The LBMcan be considered to be an efcient numerical method
for computational uid dynamics. It is also a powerful tool for modeling new
physical phenomena that are not yet easily described by macroscopic equations.
Froma computational point of view, the lattice Boltzmann equation is hyper-
bolic and can be solved locally, explicitly, and efciently on parallel computers.

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
360 CHEN & DOOLEN
Not only is the scheme computationally comparable with traditional numerical
methods, but it is also easy to program and to include new physics because of
the simplicity of the form of the LBM equations.
The lattice Boltzmann method is still undergoing development. Many mod-
els, including simulation of granular ows (Flekky & Herrmann 1993, Tan
et al 1995), viscoelastic ows (Aharonov & Rothman 1993), magnetohydro-
dynamics (Martinez et al 1994), and microemulsions (Boghosian et al 1996)
were recently proposed. Although innovative and promising, these existing
LBM methods, including multicomponent LBM models, require additional
benchmarking and verication. The current lattice models for multiphase and
reacting systems are most suitable for isothermal problems. The development
of a reliable LBM for thermal systems will allow the simulation of heat transfer
and surface phenomena simultaneously (Kato 1997). This would open many
new areas of application.
ACKNOWLEDGMENTS
We are grateful to N Cao, H Chen, B Hasslacher, X He, S Hou, S Jin, AJC
Ladd, T Lookman, L Luo, D Martinez, W Matthaeus, B Nadiga, Y Qian, X
Shan, S Succi, and MR Yeomans for useful discussions.
Visit the Annual Reviews home page at
http://www.AnnualReviews.org.
Literature Cited
Abe T. 1997. Derivation of the lattice Boltz-
mann method by means of the discrete ordi-
nate method for the Boltzmann equation. J
Comp. Phys. 131:24146
Aharonov E, Rothman DH. 1993. Non-
Newtonian ow (through porous media): a
lattice Boltzmann method. Geophys. Res.
Lett. 20:67982
Aidun CK, Lu Y-N. 1995. Lattice Boltzmann
simulation of solid particles suspended in
uid. J. Stat. Phys. 81:4961
Alexander FJ, Chen S, Grunau DW. 1993a. Hy-
drodynamic spinodal decomposition: growth
kinetics and scaling. Phys. Rev. B 48:634
37
Alexander FJ, Chen S, Sterling JD. 1993b. Lat-
tice Boltzmann thermohydrodynamics. Phys.
Rev. E 47:R224952
Alvarez-Ramirez J, Nieves-Mendoza S, Gon-
zalez-Trejo J. 1996. Calculation of the effec-
tive diffusivity of heterogeneous media using
the lattice Boltzmann method. Phys. Rev. E
53:22982303
Amati G, Succi S, Benzi R. 1997. Turbu-
lent channel ow simulation using a coarse-
grained extension of the lattice Boltzmann
method. Fluid Dyn. Res. 19:289302
Ancona MG. 1994. Fully-Lagrangian and
lattice-Boltzmann methods for solving sys-
tems of conservation equations. J. Comp.
Phys. 115:10720
Bartoloni A, Battista C, Cabasino S, Paolucci
PS, Pech J, et al. 1993. LBE simulations of
Rayleigh-Bernard convection on the APE100
parallel processor. Int. J. Mod. Phys. C4:993
1006
Behrend O. 1995. Solid-uid boundaries in par-
ticle suspension simulations via the lattice
Boltzmann method. Phys. Rev. E 52:1164
75
Benzi R, Succi S. 1990. Two-dimensional tur-
bulence with the lattice Boltzmann equation.
J. Phys. A 23:L15
Benzi R, Succi S, Vergassola M. 1992. The lat-
tice Boltzmann-equationtheory and appli-
cations. Phys. Rep. 222:14597
Benzi R, Struglia MV, Tripiccione R. 1996. Ex-
tended self-similarity in numerical simula-

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 361
tions of three-dimensional anisotropic turbu-
lence. Phys. Rev. E 53:R556568
Bhatnagar PL, Gross EP, Krook M. 1954. A
model for collision processes in gases. I:
small amplitude processes in charged and
neutral one-component system. Phys. Rev.
94:511525
Blake TD, Deconinck J, Dortona U. 1995. Mod-
els of wetting: immiscible lattice Boltzmann
automata versus molecular kinetic theory.
Langmuir. 11:458892
Boghosian BM, Coveney PV, Emerton AN.
1996. A lattice-gas model for microemul-
sions. Proc. Roy. Soc. London. Ser. A
452:122150
Brackbill JU, Kothe DB, Zemach C. 1992. A
continuum method for modeling surface ten-
sion. J. Comp. Phys. 100:33554
Broadwell JE. 1964. Study of rareed shear
owby the discrete velocity method. J. Fluid
Mech. 19:40114
Buckles J, Hazlett R, Chen S, Eggert KG, Gru-
nau DW, Soll WE. 1994. Flow through por-
ous media using lattice Boltzmann method.
Los Alamos Sci. 22:11221
Cali A, Succi S, Cancelliere A, Benzi R,
Gramignani M. 1992. Diffusionandhydrody-
namic dispersion with the lattice Boltzmann
method. Phys. Rev. A 45:577174
Cao N, Chen S, Jin S, Martinez D. 1997. Physi-
cal symmetry and lattice symmetry in lattice
Boltzmann method. Phys. Rev. E 55:R21
24
Chang YC, Hou TY, Merriman B, Osher S.
1996. A level set formulation of Eulerian in-
terface capturing methods for incompressible
uid ows. J. Comp. Phys. 124:44964
Chen H. 1993. Discrete Boltzmann systems and
uid ows. Comp. Phys. 7:63237
Chen H, Chen S, Matthaeus WH. 1992. Recov-
ery of the Navier-Stokes equations using a
lattice-gas Boltzmann method. Phys. Rev. A.
45:R533942
Chen S, Chen HD, Martinez D, Matthaeus W.
1991. Lattice Boltzmann model for simula-
tion of magnetohydrodynamics. Phys Rev.
Lett. 67:377679
Chen S, Dawson SP, Doolen GD, Janecky DR,
Lawniczak A. 1995. Lattice methods and
their applications to reacting systems. Comp.
Chem. Eng. 19:61746
Chen S, Lookman T. 1995. Growth kinetics in
multicomponent uids. J. Stat. Phys. 81:223
35
Chen S, Martinez D, Mei R. 1996. On bound-
ary conditions in lattice Boltzmann methods.
Phys. Fluids 8:252736
Chen S, Wang Z, Shan XW, Doolen GD. 1992.
Lattice Boltzmann computational uid dy-
namics in three dimensions. J. Stat. Phys.
68:379400
Chen Y, Ohashi H, Akiyama M. 1994. Thermal
lattice Bhatnagar-Cross-Krook model with-
out nonlinear deviations in macrodynamic
equations. Phys. Rev. E 50:277683
Chen Y, Ohashi H, Akiyama M. 1995. Prandtl
number of lattice Bhatnagar-Gross-Krook
uid. Phys. Fluids 7:228082
Chu CK. 1965. Kinetic-Theoretic description of
the formation of a shock wave. Phys. Fluids
8:1221
Cieplak M. 1995. Rupture and coalescence in 2-
dimensional cellular-automata uids. Phys.
Rev. E 51:435361
Cornubert R, dHumi` eres D, Levermore D.
1991. A Knudsen layer theory for lattice
gases. Physica D 47:24159
Dawson SP, Chen S, Doolen GD. 1993. Lat-
tice Boltzmann computations for reaction-
diffusion equations. J. Chem. Phys. 98:1514
23
dHumi` eres D, Lallemand P, Frisch U. 1986.
Lattice gas model for 3Dhydrodynamics. Eu-
rophys. Lett. 2:29197
Dortona U, Salin D, Cieplak M, Banavar JR.
1994. Interfacial phenomena in Boltzmann
cellular automata. Europhys. Lett. 28:31722
Dufty JW, Ernst MH. 1996. Lattice Boltzmann-
Langevin equation. Fields Inst. Comm. 6:99
107
E W-N, Liu J-G. 1996. Essentially compact
scheme for unsteady viscous incompressible
ows. J. Comp. Phys. 126:12238
Eggels JGM. 1996. Direct and large-eddy simu-
lation of turbulent uid owusing the lattice-
Boltzmann scheme. Int. J. Heat Fluid Flow
17:30723
Eggels JGM, Somers JA. 1995. Numerical
simulation of free convective ow using
the lattice-Boltzmann scheme. J. Heat Fluid
Flow 16:35764
Elton BH, Levermore CD, Rodrigue H. 1995.
Covergence of convective-diffusive lattice
Boltzmann methods. SIAM J. Sci. Comp.
32:132754
Ferreol B, Rothman DH. 1995. Lattice-
Boltzmann simulations of ow-through
Fontainebleau sandstone. Transp. Por. Med.
20:320
Flekky EG, Herrmann HJ. 1993. Lattice Boltz-
mann models for complex uids. Physica A.
199:111
Flekky EG, Rage T, Oxaal U, Feder J.
1996. Hydrodynamic irreversibility in creep-
ing ow. Phys. Rev. Lett. 77:417073
Frisch U, Hasslacher B, Pomeau Y. 1986.
Lattice-gas automata for the Navier-Stokes
equations. Phys. Rev. Lett. 56:15058
Frisch U, dHumi eres D, Hasslacher B, Lalle-
mand P, Pomeau Y, Rivet J-P. 1987. Lattice
gas hydrodynamics in two and three dimen-
sions. Complex Syst. 1:649707

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
362 CHEN & DOOLEN
Ginzbourg I, Adler PM. 1994. Boundary ow
condition analysis for the three-dimensional
lattice Boltzmann model. J. Phys. II 4:191
214
Glimm JE, Issacson E, Marchesin D, McBryan
O. 1981. Front tracking for hyperbolic sys-
tems. Adv. Appl. Math. 2:91119
Gonnella G, Orlandini E, Yeomans JM. 1997.
Spinodal decomposition to a lamellar phase:
effects of hydrodynamic ow. Phys. Rev. Lett.
78:169598
Grunau D, Chen S, Eggert K. 1993. A lattice
Boltzmann model for multiphase uid ows.
Phys. Fluids A 5:255762
Gunstensen AK, Rothman DH. 1993. Lattice
Boltzmann studies of immiscible two phase
ow through porous media J. Geophys. Res.
98:643141
Gunstensen AK, Rothman DH, Zaleski S,
Zanetti G. 1991. Lattice Boltzmann model
of immiscible uids. Phys. Rev. A. 43:4320
27
Hardy J, de Pazzis O, Pomeau Y. 1976. Molec-
ular dynamics of a classical lattice gas: trans-
port properties and time correlation func-
tions. Phys. Rev. A 13:194961
Hayot F, Wagner L. 1996. Anon-local modica-
tion of a lattice Boltzmann model. Europhys.
Lett. 33:43540
He X, Doolen GD. 1997. Lattice Boltzmann
method on curvilinear coordinates system:
vortex shedding behind a circular cylinder.
Phys. Rev. E. In press
He X, Luo L-S. 1997. A priori derivation of
the lattice Boltzmann equation. Phys. Rev. E.
55:633336
He X, Luo L-S, Dembo M. 1996. Some progress
in lattice Boltzmann method. Part I nonuni-
form mesh grids. J. Comp. Phys. 129:357
63
He X, Zou Q, Luo L-S, Dembo M. 1997. An-
alytic solutions of simple ow and analysis
of non-slip boundary conditions for the lat-
tice Boltzmann BGK model. J. Stat. Phys.
87:11536
Heijs AWJ, Lowe CP. 1995. Numerical evalua-
tion of the permeability and the Kozeny con-
stant for two types of porous media. Phys.
Rev. E 51:434652
Higuera FJ, Jim enez J. 1989. Boltzmann ap-
proach to lattice gas simulations. Europhys.
Lett. 9:66368
Higuera FJ, Succi S. 1989. Simulating the ow
around a circular cylinder with a lattice Boltz-
mann equation. Europhys. Lett. 8:51721
Higuera FJ, Succi S, Benzi R. 1989. Lattice
gas dynamics with enhanced collisions. Eu-
rophys. Lett. 9:34549
Holme R, Rothman DH. 1992. Lattice-gas and
lattice-Boltzmann models of miscible uids.
J. Stat. Phys. 68:40929
Hou S. 1995. Lattice Boltzmann method for
incompressible viscous ow. PhD thesis.
Kansas State Univ., Manhattan, Kansas
Hou S, Zou Q, Chen S, Doolen G, Cogley AC.
1995. Simulation of cavity ow by the lattice
Boltzmann method. J. Comp. Phys. 118:329
47
Hou S, Sterling J, Chen S, Doolen GD. 1996.
A lattice Boltzmann subgrid model for high
Reynolds number ows. Fields Inst. Comm.
6:15166
Hu HH. 1996. Direct simulation of ows of
solid-liquid mixtures. Int. J. Multiphase ow.
22:33552
Inamuro T, Sturtevant B. 1990. Numerical
study of discrete-velocity gases. Phys. Flu-
ids 2:21962203
Inamuro T, Yoshino M, Ogino F. 1995. A non-
slip boundary condition for lattice Boltzmann
simulations. Phys. Fluids 7:292830
Jin S, Katsoulakis M. Relaxation approxima-
tions to front propagation. J. Diff. Eq. In press
Jin S, Xin ZP. 1995. The relaxation schemes
for systems of conservation laws in arbitrary
space dimensions. Comm. Pure Appl. Math.
48:23576.
Kaandorp JA, Lowe CP, Frenkel D, Sloot PMA.
1996. Effect of nutrient diffusion and owon
coral morphology. Phys. Rev. Lett. 77:2328
31
Kadanoff L. 1986. On two levels. Phys. Today
39:79
Kato T, Kono K, Seta K, Martinez D, Chen S.
1997. Boiling two-phase ow by the lattice
Boltzmann method Int. J. Mod. Phys. In press
Kim J, Moin P, Moser R. 1987. Turbulence
statistics in fully-developed channel ow
at low Reynolds number. J. Fluid Mech.
177:13366
Koelman JMVA. 1991. A simple lattice Boltz-
mann scheme for Navier-Stokes uid ow.
Europhys. Lett. 15:6037
Koplik J, Lasseter T. 1985. One- and two-phase
ow in network models of porous media.
Chem. Eng. Comm. 26:28595
Kraichnan RH. 1967. Inertial ranges in two-
dimensional turbulence. Phys. Fluids 10:
141723
Ladd AJC. 1993. Short-time motion of colloidal
particles: numerical simulation via a uctu-
ating lattice-Boltzmann equation Phys. Rev.
Lett. 70:133942
Ladd AJC. 1994a. Numerical simulations of
particulate suspensions via a discretized
Boltzmannequation. Part. 1. theoretical foun-
dation. J. Fluid Mech. 271:285309
Ladd AJC. 1994b. Numerical simulations of
particulate suspensions via a discretized
Boltzmann equation. Part 2: Numerical re-
sults. J. Fluid Mech. 271:31139
Ladd AJC. 1997. Sedimentation of homoge-

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
LATTICE BOLTZMANN METHOD 363
neous suspensions of non-Brownian spheres
Phys. Fluids 9:49199
Laudau LD, Lifshitz EM. 1959. Fluid Dyn. New
York: Pergamon
Lavall ee P, Boon JP, Noullez A. 1991. Bound-
aries in lattice gas ows. Physica D 47:233
40
Luo L. 1997. Symmetry breaking of ow in 2-
Dsymmetric channels: simulations by lattice
Boltzmann method. Int. J. Mod. Phys. C. In
press
Maier RS, Bernard RS, Grunau DW. 1996.
Boundary conditions for the lattice Boltz-
mann method. Phys. Fluids 8:17881801
Majda A. 1984. Compressible Fluid and Sys-
tems of Conservation Laws in Several Space
Variables. New York: Springer-Verlag
Martys NS, Chen H. 1996. Simulation of
multicomponent uids in complex three-
dimensional geometries by the lattice Boltz-
mann method. Phys. Rev. E 53:74350
Martinez DO, Chen S, Matthaeus WH. 1994.
Lattice Boltzmann magnetohydrodynamics.
Phys. Plasmas 1:185067
Martinez DO, Matthaeus WH, Chen S, Mont-
gomery DC. 1994. Comparison of spectral
method and lattice Boltzmann simulations of
two-dimensional hydrodynamics. Phys. Flu-
ids 6:128598
Massaioli F, Benzi R, Succi S. 1993. Expo-
nential tails in two-dimensional Rayleigh-
B enard convection. Europhys. Lett. 21:305
10
McNamara GR, Garcia AL, Alder BJ. 1995.
Stabilization of thermal lattice Boltzmann
models. J. Stat. Phys. 81:395408
McNamara GR, Zanetti G. 1988. Use of the
Boltzmann equation to simulate lattice-gas
automata. Phys. Rev. Lett. 61:233235
Miller W. 1995. Flowin the driven cavity calcu-
lated by the lattice Boltzmann method. Phys.
Rev. E 51:365969
Nadiga BT. 1995. An Euler solver based on
locally adaptive discrete velocities. J. Stat.
Phys. 81:12946
Nadiga BT, Pullin DI. 1994. A method for
near-equilibrium discrete-velocity gas ows.
J. Comp. Phys. 112:16272
Nadiga BT, Zaleski S. 1996. Investigations of a
two-phase uidmodel. Eur. J. Mech. B/Fluids
15:88596
Nannelli F, Succi S. 1992. The lattice
Boltzmann-equation on irregular lattices. J.
Stat. Phys. 68:4017
Noble DR, Chen S, Georgiadis JG, Buckius RO.
1995. A consistent hydrodynamic boundary
condition for the lattice Boltzmann method.
Phys. Fluids 7:2039
Noble DR, Georgiadis JG, Buckius RO. 1996.
Comparison of accuracy and performance for
lattice Boltzmann and nite difference simu-
lations of steady viscous ow. Int. J. Numer.
Meth. Fluids 23:118
Orszag SA, Yakhot V. 1986. Reynolds number
scaling of cellular automaton hydrodynam-
ics. Phys. Rev. Lett. 56:169193
Osborn WR, Orlandini E, Swift MR, Yeomans
JM, Banavar JR. 1995. Lattice Boltzmann
study of hydrodynamic spinodal decompo-
sition. Phys. Rev. Lett. 75:403134
Prasad AK, Koseff JR. 1989. Reynolds-number
and end-wall effects on a lid-driven cavity
ow. Phys. Fluids A 1:20818
Qian YH. 1990. Lattice gas and lattice kinetic
theory applied to the Navier-Stokes equa-
tions. PhD thesis. Universit e Pierre et Marie
Curie, Paris
Qian YH. 1993. Simulating thermohydrody-
namics with lattice BGK models. J. Sci.
Comp. 8:23141
Qian YH, dHumi` eres D, Lallemand P. 1992.
Lattice BGKmodels for Navier-Stokes equa-
tion. Europhys. Lett. 17:47984
Qian YH, Orszag SA. 1993. Lattice BGK mod-
els for the Navier-Stokes equation: nonlinear
deviationincompressible regimes. Europhys.
Lett. 21:25559
Qian YH, Succi S, Massaioli F, Orszag SA.
1996. A benchmark for lattice BGK model:
owover a backward-facing step. Fields Inst.
Comm. 6:20715
Qian YH, Succi S, Orszag SA. 1995. Recent
advances in lattice Boltzmann computing.
Annu. Rev. Comp. Phys. 3:195242
Rothman DH. 1988. Cellular-automaton uids:
a model for ow in porous media. Geophys.
53:50918
Rothman DH, Keller JM. 1988. Immisci-
ble cellular-automaton uids. J. Stat. Phys.
52:111927
RothmanDH, Zaleski S. 1994. Lattice-gas mod-
els of phase separation: interfaces, phase
transitions and multiphase ow. Rev. Mod.
Phys. 66:141779
Rybka RB, Cieplak M, Salin D. 1995. Boltz-
mann cellular-automata studies of the spin-
odal decomposition. Physica A 222:10518
Reider MB, Sterling JD. 1995. Accuracy of
discrete-velocity BGK models for the sim-
ulation of the incompressible Navier-Stokes
equations. Comput. Fluids 24:45967
Schreiber R, Keller HB. 1983. Driven cavity
ows by efcient numerical techniques. J.
Comp. Phys. 49:31033
Schwarz LM, Martys N, Bentz DP, Garboczi
EJ, Torquato S. 1993. Cross-property rela-
tions and permeability estimation in model
porous media. Phys. Rev. E 48:458491
Segre PN, Behrend OP, Pusey PN. 1995. Short-
time Brownian motion in colloidal suspen-
sions: experiment and simulation. Phys. Rev.
E 52:507083

P1: ARK/ary P2: MBL/vks QC: MBL/bsa T1: MBL
November 24, 1997 15:23 Annual Reviews AR049-12
364 CHEN & DOOLEN
Shan X. 1997. Simulation of Rayleigh-B enard
convection using a lattice Boltzmann method.
Phys. Rev. E 55:278088
Shan X, Chen H. 1993. Lattice Boltzmann
model for simulating ows with multi-
ple phases and components. Phys. Rev. E
47:181519
Shan X, Chen H. 1994. Simulation of non-
ideal gases and liquid-gas phase transitions
by the lattice Boltzmann equation. Phys. Rev.
E 49:294148
Shan X, Doolen G. 1995. Multicomponent
lattice-Boltzmann model with interparticle
interaction. J. Stat. Phys. 81:37993
Skordos PA. 1993. Initial and boundary condi-
tions for the lattice Boltzmann method. Phys.
Rev. E 48:482342
Soll W, Chen S, Eggert K, Grunau D, Ja-
necky D. 1994. Applications of the lattice-
Boltzmann/lattice gas techniques to multi-
uid owin porous media. In Computational
Methods in Water Resources X, ed. A Peters,
pp. 99199. Dordrecht: Academic
Somers JA. 1993. Direct simulation of uid
ow with cellular automata and the lattice-
Boltzmann equation. Applied Sci. Res.
51:12733
Spaid MAA, Phelan FR Jr. 1997. Lattice
Boltzmann methods for modeling microscale
ow in brous porous media. Phys. Fluids.
9:246874
Sterling JD, Chen S. 1996. Stability analysis of
lattice Boltzmann methods. J. Comp. Phys.
123:196206
Swift MR, Osborn WR, Yeomans JM. 1995.
Lattice Boltzmannsimulationof nonideal u-
ids. Phys. Rev. Lett. 75:83033
Swift MR, Orlandini SE, Osborn WR, Yeo-
mans JM. 1996. Lattice Boltzmann simula-
tions of liquid-gas and binary-uid systems.
Phys. Rev. E 54:504152
Succi S, Amati G, Benzi R. 1995. Challenges
in lattice Boltzmann computing. J. Stat. Phys.
81:516
Succi S, Benzi R, Higuera F. 1991. The lattice-
Boltzmann equationa new tool for compu-
tational uid dynamics. Physica D 47:219
30
Succi S, Foti E, Higuera F. 1989. Three-
dimensional ows in complex geometries
with the lattice Boltzmann method. Euro-
phys. Lett. 10:43338
Succi S, Nannelli F. 1994. The nite volume for-
mulation of the lattice Boltzmann equation.
Transp. Theory Stat. Phys. 23:16371
Tan M-L, Qian YH, Goldhirsch I, Orszag SA.
1995. Lattice-BGK approach to simulating
granular ows. J. Stat. Phys. 81:87103
Trevi no C, Higuera F. 1994. Lattice Boltzmann
and spectral simulations of non-linear sta-
bility of Kolmogorov ows. Rev. Mex. Fis.
40:87890
Vahala G, Pavlo P, Vahala L, Chen H. 1995.
Effect of velocity shear on a strong tempera-
ture gradienta lattice Boltzmann approach.
Phys. Lett. A 202:37682
Vangenabeek O, Rothman DH. 1996. Macro-
scopic manifestations of microscopic ows
through porous-media phenomenology from
simulation. Annu. Rev. Earth Planet. Sci.
24:6387
Vanka SP. 1986. Block-implicit multigrid solu-
tion of Navier-Stokes equations in primitive
variables. J. Comp. Phys. 65:13858
Wagner L. 1994. Pressure in lattice Boltzmann
simulations of ow around a cylinder. Phys.
Fluids 6:351618
Weimar JR, Boon JP. 1996. Nonlinear reac-
tions advected by a ow. Physica A 224:207
15
Wolfram S. 1986. Cellular automaton uids. 1:
Basic theory. J. Stat. Phys. 45:471526
Wu Y, Alexander FJ, Lookman T, Chen S. 1995.
Effects of hydrodynamics on phase transi-
tion kinetics in two-dimensional binary u-
ids. Phys. Rev. Lett. 74:385255
Xu K, Prendergast KH. 1994. Numerical Na-
vier-Stokes solutions fromgas kinetic theory.
J. Comp. Phys. 114:917
Ziegler DP. 1993. Boundary conditions for lat-
tice Boltzmann simulations. J. Stat. Phys.
71:117177
Zou Q, He X. 1997. On pressure and veloc-
ity boundary conditions for the lattice Boltz-
mann BGK model. Phys. Fluids. 9:159198
Zou Q, Hou S, Doolen GD. 1995. Analytical so-
lutions of the lattice Boltzmann BGK model.
J. Stat. Phys. 81:31934