1

An Introduction to Lattice Boltzmann Methods

1. The Study of Fluid Motion
Fluid motion can take a variety of forms ranging from simple flows such as laminar flow in a
pipe, to more complex flows such as vortex shedding behind cylinders, wave motion and
turbulence. It incorporates both liquid and gaseous flows. Many of the different flow situations
have been examined experimentally; however it is advantageous to develop a numerical model
capable of simulating the many flow structures experienced in the motion of different fluids.

2. Numerical Methods in Fluid Study
Computational fluid dynamics (CFD) has developed mainly around using numerical techniques
to solve the Navier-Stokes equation and the continuity equation or an equation derived from
them. Another approach which has proved less popular is the molecular dynamics approach.

2.1 Numerical Solutions of the Navier-Stokes Equation
The most popular method in CFD is the numerical solution of the Navier-Stokes equation. Given
the Navier-Stokes equation and a set of suitable boundary conditions it is possible to solve on a
grid using the standard numerical techniques. This works well for simple flows; however more
complex problems frequently require a more complex approach.

2.2 Molecular Dynamics
One obvious way to simulate a fluid on a computer is to model the individual molecules which
make up the fluid. Then, provided the inter-molecular interactions are modeled correctly, the
system should behave as a fluid. Different situations can be modeled by changing the average
energy of the molecules and their separation.
The main disadvantage with such an approach is that large computer resources are required,
many simulations taking hours to evolve a fraction of a second. The system must be updated in
small time-steps, the new position and velocity of all particles being calculated, at every time-
step, from a knowledge of their previous position and velocity, taking into account any external
forces which are acting on them. Any particles which collided during the previous time-step have
to be identified and their new trajectories calculated. This can be restrictively time consuming
when considering even a very small volume of fluid. Even when a gas is being considered where
there are fewer molecules and a larger time-step can be used, because of the longer mean free
path of the molecules, the number of molecules which can be considered is severely limited.

2.3 Lattice Gas Method
Over the past fifteen years, a new method for the computer simulation of fluids has been
developed: the lattice gas model. Instead of considering a large number of individual molecules,
the molecular dynamics approach, a much smaller number of fluid particles are considered. A
fluid particles is a large group of molecules which although much larger than a molecule is still
considerably smaller then the smallest length scale of the simulation. This reduces the amount of
data which needs to be stored since large simulations can be performed using less than one
million particles. This is justified on the grounds that the macroscopic properties do not depend
2
directly on the microscopic behavior of the fluid. This can be seen in low Mach flows where,
provided the Reynolds number is the same, experiments carried out in a water tank and a wind
tunnel produce the same results. These two fluids have different microscopic structures, but they
both exhibit the same macroscopic features. In a lattice gas model the `particles' are restricted to
move on the links of a regular underlying grid and the motion evolves in discrete time-steps. The
conservation laws are incorporated into update rules which are applied at each discrete time.
A lattice gas model in which the state of the fluid needs to be known only at the lattice sites and
only at discrete times can run much faster on a computer than a molecular dynamics simulation.
The lattice gas model has another big advantage over molecular simulation since all the
collisions occur at the same time. This is a particular advantage if the simulation is being run on
a parallel computer. These two time saving advantages of the lattice gas model allow simulations
of a significantly large scale to be performed.

2.4 Lattice Boltzmann Method

The lattice Boltzmann method is a powerful technique for the computational modeling of a wide
variety of complex fluid flow problems including single and multiphase flow in complex
geometries. It is a discrete computational method based upon the Boltzmann equation. It
considers a typical volume element of fluid to be composed of a collection of particles that are
represented by a particle velocity distribution function for each fluid component at each grid
point. The time is counted in discrete time steps and the fluid particles can collide with each
other as they move, possibly under applied forces. The rules governing the collisions are
designed such that the time-average motion of the particles is consistent with the Navier-Stokes
equation

3. Why is the Lattice Boltzmann Method Important?

In recent years, the Lattice Boltzmann Method (LBM) has emerged as a promising numerical
method for simulating fluid flows. Unlike conventional methods which solve the discretized
macroscopic Navier-Stokes equations, the LBM is based on microscopic particle models and
mesoscopic kinetic equations. The fundamental concept of the LBM is to construct simplified
kinetic models that incorporate the essential physics of microscopic or mesoscopic processes so
that the macroscopic averaged properties obey the desired macroscopic equations. The LBM is
especially useful for modeling interfacial dynamics, flows over porous media, and multi-phase
flows. In addition, the LBM algorithm tends to be very simple, allowing parallelism in a
straightforward manner.
The LBM is a derivative of the lattice gas automata method which was first proposed about a
dozen years ago by a number of physicists. Nowadays, the method has quickly found its way in
dealing with a number of engineering flow problems. However, the underlying lattice gas (or
mesh) is a rectangle (or a hexagon) in two dimensions or a cube (or a shape with perfect
geometric symmetry) in three dimensions, equivalent to the regular Cartesian grid used by
conventional Navier-Stokes solvers. As a result, solution domains with inclined or curved
boundaries are approximated by staircase-like steps. This restriction severely limits the
applicability of LBM as most of the industrial and practical flows have complex flow
geometries. Looking back on the history of the conventional Navier-Stokes methods, the issue of
mesh flexibility dominates its development. It begins with the Cartesian regular grid, then the use
3
of body-fitted coordinate grids and structured multi-block grids, and finally the widespread
acceptance of unstructured grids, allowing the greatest flexibility in adapting the grid to domain
boundaries. The body-fitted method and the multi-block structured method are merely special
cases of the more general unstructured mesh.

This method naturally accommodates a variety of boundary conditions such as the pressure drop
across the interface between two fluids and wetting effects at a fluid-solid interface. It is an
approach that bridges microscopic phenomena with the continuum macroscopic equations.
Further, it can model the time evolution of systems.

3. Evolution of Lattice Boltzmann Method (LBM)
3.1 The Classical Boltzmann Equation

A statistical description of a fluid system can be made in terms of a distribution function
( ) , , f t x where ( ) , , f t x is defined such that ( ) , , f t d d x x is the number of molecules at
time t positioned between x and d + x x which have velocities in the range d + . Consider
a gas in which an external force per unit mass F acts and assume initially that no collisions take
place between the gas molecules. In time dt the velocity of any molecule will change to
dt + F and its position x will change to dt + x . Thus the number of molecules ( ) , , f t d d x x
is equal to the number of molecules ( ) , , f dt dt t dt d d + + + x F x , that is to say,
( ) ( ) , , , , 0 f dt dt t dt d d f t d d + + + = x F x x x (1)
If, however, collisions do occur between the molecules there will be a net difference between the
number of molecules ( ) , , f dt dt t dt d d + + + x F x and the number of molecules
( ) , , f t d d x x . This can be written as ( ) f d d dt O x where ( ) f O is the collision operator. This
gives the following equation describing the evolution of the distribution function:
( ) ( ) ( ) , , , , f dt dt t dt d d f t d d f d d dt + + + = O x F x x x x (2)
Dividing Eq. 2 by d d dt x and letting 0 dt gives the Boltzmann equation:
( )
( )
f f f
f
t
f
f f f
t

c c c
+ + = O
c c c
c
+ V + V = O
c
F
x
F
(3)
The macroscopic variables, like density ( ), velocity ( u) and internal energy ( 2 DRT c = , with
D is the dimension, R is the universal gas constant and T is the macroscopic temperature) can be
obtained by taking (microscopic velocity) moments of the distribution function f as:
f d =
}
(4)
f d =
}
u (5)
( )
2 1
2
f d c =
}
u (6)
4
Any solution of Boltzmann's equation (Eq. 3) requires that an expression is found for the
collision operator ( ) f O . If the collision is to conserve mass, momentum and energy it is
required that:
( )
2
1
0 f d
(
( O =
(

}

(7)
The form of the collision function can be found by assuming that the gas has a low density so
only binary collisions need be considered. It is also assumed that the molecules are completely
uncorrelated before the collision, this assumption is called molecular chaos. Still it is found that
the collision operator has a very complicated form. It is therefor assumed that ( ) f O can be
replaced by a simplified collision operator which retains only the qualitative and average
properties of the actual collision operator. Any replacement collision function must satisfy the
conservation of mass, momentum and energy expressed by Eq. 7. Such an operator is based on
the idea of a single relaxation time and can be written as:
( )
( ) ( ) , , , ,
eq
f t f t
f
t

O =
x x
(8)
where t is the relaxation time which is of the order of the time between collisions and
( ) , ,
eq
f t x is the Maxwell-Boltzmann equilibrium distribution function. This model is
frequently called the BGK model after Bhatnagar, Gross and Krook who first introduced it.
Although
eq
f is written as an explicit function of t, the time dependence of
eq
f lies solely in
the macroscopic variables , u and T, i.e., ( ) ( ) , , , ; , ,
eq eq
f t f T = x x u , where
( )
( )
2
2
exp
2
2
eq
D
f
RT
RT

t
(

= (
(

u
(9)
With the help of Eq. 8, the continuous Boltzmann equation (Eq. 3) can now be written as:
eq
f f f
f f
t

t
c
+ V + V =
c
F (10)
It is well known that, using the multi-scale Chapman-Enskog expansion the Boltzmann-BGK
equation (Eq. 10) recovers the macroscopic continuity, momentum and energy equations at the
Navier-Stokes level. The Chapman-Enskog expansion parameter is the Knudsen number, K,
defined as the ratio of the mean free path to a typical macroscopic length scale. For small
Knudsen number, K can be introduced into the Boltzmann equation as:
eq
f f f
f f
t K

t
c
+ V + V =
c
F (11)
Setting,
( )
0
n n
n
f K f

=
= , we look for the solution of Eq. 11 such that,
5
( )
( )
( )
( )
2
2
1
for 0
2
1
0 for 1
2
n
n
f d n
f d n

c
(
(
(
= = (
(
(

(

(
(
= >
(

(

}
}
u
u

u
(12)
The zeroth order term,
( ) 0
f , is taken to be the local Maxwell-Boltzmann distribution,
eq
f , and
the corresponding higher order terms are so chosen that they have no contribution to the
moments expressed in Eq. 12. Hence the macroscopic variables, which are obtained by taking
(microscopic velocity) moments of the distribution function f, can be described as:
eq
f d f d = =
} }
(13)
eq
f d f d = =
} }
u (14)
( ) ( )
2 2 1 1
2 2
eq
f d f d c = =
} }
u u (15)
The first-order solution can be found by considering
( )
1
O K

, which gives the Euler equation.

The second-order solution found by considering
( )
0
O K and gives the Navier-Stokes equation,
when the binary collision function is used.
( ) 0
t

c
+V =
c
u (16)
( )
( ) ( ) p
t

v
c
+V = V +V V + V + (

c
u
uu u u F (17)
( )
( ) ( ) ( ) p
t
c
c _ c v
c
+V = V V + V + V : V V (

c
u u u u u (18)
where p RT = is the pressure, RT v t = is the kinematic viscosity and ( ) 2 D RT D _ t = + is
the thermal conductivity.
It needs to be noted here that the Boltzmann equation with single-relaxation-time BGK model
does have one unsatisfactory feature: the energy equation obtained from the second moment of
the distribution function yields a fixed Prandtl number, implying that the thermal diffusivity
cannot be adjusted independent of the kinematic viscosity, which restricts its applicability to a
limited class of problems only. Hence, for simulating general fluid flows with arbitrary Prandtl
numbers, the distribution function f should not be used to calculate the internal energy or
temperature. For obtaining temperature field in case of generalized fluid flow problems, one
evolution equation of internal energy density distribution function needs to be developed, which
is beyond the scope of present discussion. The present model deals only with the isothermal fluid
flow problems.
In summary, the following equations are proposed for simulating a hydrodynamic problem:
eq
f f f
f f
t

t
c
+ V + V =
c
F
6
where
( )
( )
2
2
exp
2
2
eq
D
f
RT
RT

t
(

= (
(

u

and the macroscopic variables are calculated using
( )
2 1
2
f d
f d
f d

c
=
=
=
}
}
}

u
u


3.2 The Lattice Boltzmann Equation (LBE)
The lattice Boltzmann equation (LBE) can be directly derived from the continuous Boltzmann
equation discretized in some special manner in both time and phase space. The analysis shows
that theoretically the lattice Boltzmann equation is independent of the lattice gas automata
(LGA). The lattice Boltzmann equation is a finite difference form of the continuous Boltzmann
equation.
In the following analysis, the external force term is truncated and with this assumption the
continuous Boltzmann equation takes the form:
eq
f f f
f
t t
c
+ V =
c
(19)
Eq. 19 can be formally rewritten in the form of an ordinary differential equation:
1 1
eq
df
f f
dt t t
+ = (20)
where
d
dt t
c
+ V
c

Eq. 20 is a Leibniz linear differential equation. Integrating Eq. 20 over a time step,
t
o to yield:
( ) ( ) ( )
0
1
, , , , , ,
t
t t
t eq
t t
f t e e f t t t dt e f t
o
o t o t t
o o
t
'
' ' ' + + = + + +
}
x x x (21)
Assuming that
t
o is small enough and g is smooth enough locally, the following approximation
can be made:
( ) ( ) ( ) ( )
2
, , 1 , , , ,
eq eq eq
t t t
t t
t t
f t t t f t f t O o o o
o o
| | ' '
' ' + + = + + + +
|
\ .
x x x , 0
t
t o ' s s (22)
The leading terms neglected in the above approximation are of the order of
2
t
o . With this
approximation, Eq. 21 becomes
( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
, , , , 1 , , , ,
1 1 , , , ,
t
t
eq
t t
eq eq
t t
t
f t f t e f t f t
e f t f t
o t
o t
o o
t
o o
o

( + + = +

(
( + + +
(


x x x x
x x
(23)
If we expand
t
e
o t
in its Taylor expansion and, further, neglect the terms of order
( )
2
t
O o or
smaller on the right-hand side of Eq. 23, then Eq. 23 becomes:
7
( ) ( ) ( ) ( ) , , , , , , , ,
eq t
t t
f t f t f t f t
o
o o
t
( + + =

x x x x (24)
Therefore, Eq. 24 is accurate to the first order in
t
o . Equation 24 is the evolution equation of the
distribution function f with discrete time.
To obtain the hydrodynamic moments described by Eqs. 4-6, the velocity space, denoted by ,
must be discretized. With appropriate discretization, integration in momentum space (with
weight function
eq
f ) can be approximated by quadrature up to a certain degree of accuracy, that
is,
( ) ( ) ( ) ( ) , , , ,
eq eq
i i i
i
f x t d w f t =

}
e x e (25)
where ( ) is a polynomial in ,
i
w is the weight coefficient of the quadrature and
i
e is the
discrete velocity set or abscissa of the quadrature. Accordingly, the hydrodynamic moments can
be computed as:

eq
i i
i i
f f = =

(26)

eq
i i i i
i i
f f = =

u e e (27)
( ) ( )
2 2 1 1
2 2
eq
i i i i
i i
f f c = =

e u e u (28)
For calculating the above hydrodynamic moments and to recover the Navier-Stokes equations, a
proper equilibrium distribution function is needed. In the lattice Boltzmann equation, the
equilibrium distribution function is obtained by a truncated small velocity expansion (or low
Mach number approximation). Accordingly,
eq
f is expanded for a low Mach-number flow
(i.e.
2
RT u ), as:
( )
( )
( )
( )
( )
2 2
2
2
2 2
3
2 2
exp exp
2 2
2
exp 1 O
2 2
2 2
eq
D
D
f
RT RT RT
RT
RT RT RT
RT RT

t
| | | |
=
| |
\ . \ .
(
| |
= + + + (
|
( \ .

u u
u
u u
u
(29)
In deriving the Navier-Stokes equations from the Boltzmann equation via the Chapman-Enskog
analysis, the first two order approximations of the distribution function (i.e.
( ) 0
f and
( ) 1
f ) must
be considered. Therefore, given the equilibrium distribution function,
eq
f of Eq. (23), the
quadrature used to evaluate the hydrodynamic moments must be able to compute the following
moments with respect to
eq
f exactly:

: 1, ,
: , ,
i i j
i i j i j k

u
(30)
where
i
is the component of in Cartesian coordinates.
Thus to obtain the Navier-Stokes equation from the isothermal model, the moments that are to be
evaluated are
5
1, ,..., . Calculating the hydrodynamic moments of
eq
f is equivalent to
evaluating the following integral in general:
8
0
1
2
3
4
5 6
7 8
x
y
( )
( )
( )
( )
( )
2
2
2
2
2
exp
2
2
1
2
2
eq
D
I f d
RT
RT
d
RT RT
RT

t
| |
= =
|
\ .
(

+ + (
(

} }


u
u u

(31)
The integral has the structure ( )
2
x
e x dx

}
, which can be calculated numerically with Gaussian-
type quadrature. Applying a third-order Gauss-Hermite quadrature formula for a two-
dimensional case, the above integral can be evaluated to yield a d2q9 LBE model (see
Appendix), with the following set of discrete velocities and weight coefficients (refer to Fig. 1):
( ) ( ) ( )
( ) ( ) ( )
0
cos 1 2 , sin 1 2 1, 2, 3, 4
2 cos 5 2 4 , sin 5 2 4 5, 6, 7, 8
i
i
i i c i
i i c i
t t
t t t t

= = ( (

+ + = ( (

0
e (32a)
where
( )
3 c RT = is the characteristic speed, and
4 9 0
1 9 1, 2, 3, 4
1 36 5, 6, 7, 8
i
i
w i
i
=

= =

(32b)
















Fig. 1 Discretized 2D velocity space (d2q9)

The equilibrium distribution function, described by Eq. (29), can accordingly be simplified to
obtain (see Appendix):
( ) ( )
2
2
2 4 2
3 9
3
1
2 2
i i eq
i i
f w
c c c

(

= + + (
(

e u e u
u
(33)
Finally, the evolution equation of the discrete distribution function
i
f can now be given as:
9
( ) ( ) ( ) ( ) , , , ,
eq t
i i t t i i i
f t f t f t f t
o
o o
t
( + + =

x e x x x (34)
Equation 34 is the final form of the isothermal lattice Boltzmann equation (LBE).
It is to be recognized here that the pressure in the LB model is calculated from an ideal gas law
and is not an independent dynamical variable. Further, the LBE always simulates the
compressible Navier-Stokes equation instead of the incompressible one, because the spatial
density variation is not zero in LBE simulations. Hence an appropriate modification in the LBE
is required for simulating incompressible flows. In an incompressible limit, the density
fluctuation is small (
2
M with M being the Mach number) and accordingly the equilibrium
distribution becomes,
( ) ( )
2
2
0 2 2 4 2
0
3 9
3 3
2 2
i i eq
i i
p
f w
c c c c

(

= + + (
(

e u e u
u
(35)
where
0
is the constant part of the density and the pressure p replaces the density as the
primary variable and given as:
eq
i
p RT f =

.
Finally, to implement the LBM on digital computers, the discretized (temporal, velocity spaces)
Boltzmann equation must be discretized in physical space by a series of grid nodes. The
discretization in the physical space should be such that the calculated distributions must reside on
the grid nodes. To accomplish this, the physical space can be discretized into a regular lattice so
that every
node i t
o + x e is another grid node. The general practice is to choose the lattice constant
as
x t
c o o = , which ensures that the information at all the grid nodes will automatically be
conveyed at the next time step.
4. Chapman-Enskog Expansion
The essential spirits of the Chapman-Enskog expansion are to reveal the macroscopic behavior
of the lattice Boltzmann equation deviating from the equilibrium distribution when disturbed by
a small perturbation. The perturbation is introduced through the Knudsen number ( ) 1 K .
The Chapman-Enskog expansion is first performed over the isothermal LBE to recover the
macroscopic continuity and momentum equations followed by a recovery of the macroscopic
energy equation from the thermal LBE. This is systematically achieved as follows.
First, a Taylor expansion of ( ) ,
i i t t
f t o o + + x e yields,
( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
2 2
3
, ,
1 1
, O
2 2
i i t t i
t t t i t i i t t t i t i t
f t f t
f t
o o o o | | o o
o o
o o o o o
+ +
(
= c + c + c c + c + c c +c +
(

x e x
e e e e x
(36a)
Expanding the distribution function and the time and space derivatives in terms of the Knudsen
number, one obtains:

( )
0
n n
i i
n
f K f

=
=

(36b)

1
n
t nt
n
K

=
c = c

(36c)

1
n
n
n
K

=
c = c
x x
(36d)
10
The isothermal LB equation is rewritten as:
( ) ( ) ( ) ( ) , , , ,
eq t
i i t t i i i
f t f t f t f t
o
o o
t
( + + =

x e x x x (37)
Substituting Eq. 36 into Eq. 37, and in the consecutive order of parameters K, it is possible to
arrive at the following equations:
( )
( ) 0 0
O :
eq
i i
K f f = (38a)
( )
( ) ( ) ( ) 0 0 1 1
1 1
1
O :
t i i i i
K f f f
o o
t
c + c = e (38b)
( )
( ) ( ) ( ) ( ) ( )
( )
( ) ( )
( )
( )
0 1 1 0 0 2
2 1 1 1 1 1
0 0 2
1 1 1
1
O :
2
1 1
2
t i t i i i t t t i i i
t t i i i i i i
K f f f f f
f f f
o o o o
o o | | o
o
o
t
c + c + c + c c + c
+ c c + c =
e e
e e e
(38c)
where the notation ( )
1 1o
o
c = c
x
has been used.
The distribution function f
i
is the normal solution, which is constrained by:

( ) 0
1
i
i i
f

( (
=
( (

e u
(39a)

( )
1
0, 0
n
i
i i
f n
(
= >
(

e
(39b)
For a 9-bit model, the tensor
( )
1 2
0
...
n
i i i in
i
E w
=
=

e e e , where
io
e is the projection of
i
e on o -axis
( or x y o = ), have the following properties:

( ) 2 2
0
1
3
i i i
i
E w c
o | o|
o
=
= =

e e (40a)

( )
( )
4 4
0
1
9
i i i i i
i
E w c
o | o o| o o |o oo |
o o o o o o
=
= = + +

e e e e (40b)
because

4
2
1
2
i i
i
c
o | o|
o
=
=

e e (41a)

8
2
5
4
i i
i
c
o | o|
o
=
=

e e (41b)

4
4
1
2
i i i i
i
c
o | o o|o
o
=
=

e e e e (41c)

( )
8
4 2
5
4 8
i i i i
i
c c
o | o o| o o |o oo | o|o
o o o o o o o
=
= + +

e e e e (41d)
where
o|
o and
o|o
o are the Kronecker delta with two and four indices, respectively.
Also,
( ) 2 1
0 for 0,1,...
n
E n
+
= =
With the above properties of the tensor
( ) n
E , one can have

( ) 0
i
i
f =

(42a)
11

( ) 0
i i
i
f =

e u (42b)

( ) 0 2
1
3
i i i
i
f c u u
o | o| o |
o = +

e e (42c)

( )
( )
0 2
1
3
i i i i
i
f c u u u
o | o| o | | o
o o o = + +

e e e (42d)
Summing Eq. 38b and using Eqs. 42a-b, one obtains:

1 1
0
t
u
o o
c + c = (43)
Multiplying Eq. 38b by
i|
e , one obtains,

( ) ( ) ( ) 0 0 1
1 1
1
t i i i i i i i
f f f
| o | o |
t
c + c = e e e e (44)
Summing Eq. 44 and using Eqs. 42b-c, one obtains:

2
1 1 1
1
3
t
u u u c
| o o | o o|
o
(
c + c = c
(

(45)
Summing Eq. 38c and noting that the second and third terms vanish due to the conservation of
mass and momentum and also fourth and fifth terms are zero from Eqs. 43, 45, one obtains:

2
0
t
c = (46)
Multiplying Eq. 38c by
i
e and summing over I, one obtains:
( ) ( ) ( ) ( ) ( )
( )
( ) ( )
( )
( )
0 1 1 0 0
2 1 1 1 1 1
0 0 2
1 1 1
1
2
1 1
2
t i i t i i i i i t t t i i i i i
i
t t i i i i i i i i i
f f f f f
f f f
o o o o
o o | | o
o
o
t

c + c + c + c c + c

(
+ c c + c =
(

e e e e e e e
e e e e e e
(47)
The second term in left hand side of Eq. 47 is zero by definition of
( ) 1
i
f , and the fourth term is
also zero by Eq. 44. The fifth term is given by Eq. 42c-d, while the third term can be found by
considering Eq. 38b multiplied by
1 i i | |
c e e and summed over i to the order ( ) O u , as:

( )
( )
1 2 2
1 1 1 1 1
1 1
3 3
i i i t
i
f c c u u u
o o o o| o | o| o | | o
t o o o o
(
c = c c + c c + +
(

e e (48)
Using Eq. 48, Eq. 47 to the order ( ) O u becomes,
12
( )
( )
( )
2 2
2 1 1 1 1
2 2
1 1 1 1
2
2 1 1
2
1 1
2 1
1 1
3 3
1 1 1 1
0
2 3 2 3
1 1
3 2
1 1
0
3 2
t t
t t t
t t t
t
t
u c c u u u
c c u u u
u c
c u u u
u
o o| o | o| o | | o
o o| o | o| o | | o
o o|
o | o| o | | o

t o o o o
o o o o o o
t o o
t o o o o
v
(
c c c + c c + +
(

+ c c + c c + + =
(
c c c
(

(
c c + + =
(

c = c
( )
( )
1 1 1
2 1 1 1 1
t
t
u u u
u u u u u
o o| o | o| o | | o
o o o| o | o | | o o|
o v o o o
v o v o o o
c + c c + +
c = c c + c c + +
(49)
where v is the kinematic shear viscosity given as:
( )
2
1 1
0.5
3 2
t t
c RT v t o t o
(
= =
(

(50)
At this point it should be remembered that the viscosity has changed from RT t to
( ) 0.5
t
RT t o . This can be attributed by the fact that the collision operator in the Boltzmann-
BGK equation is assumed as a constant during each time step. This assumption introduces a
second-order truncation error in the LBE which causes the viscosity to change.
Combining the first and second order density and momentum equations and recombining the
derivatives with the expansion parameter (by setting 1 K = ), gives the macroscopic continuity
and momentum equations, as:
( ) 0
t
c +V = u (51a)
( ) ( ) ( )
t
p v c +V = V +V V + V (

u uu u u (51b)
where
2
1
3
p c RT = = , is the pressure.

5. Example: Heat conduction problem

The macroscopic energy conservation equation for a heat conduction problem can be given as:
( ) ( )
t
CT k T c = V V (52a)
2
t
T T o c = V (52b)
The corresponding LBGK model is,
( ) ( ) ( ) ( )
1
, , , ,
eq
t i i i i i
f t f t f t f t
t
( c + V =

x e x x x (53)
Equation (47) can be discretized in time to yield the evolution equation of the particle
distribution function as:
( ) ( ) ( ) ( ) ( ) ( )
, , , ,
eq
i i i i i
f t t t f t t f t f t t + A + A = A x e x x x (54)
The pertinent macroscopic physical quantities, subsequently, can be obtained from the above
particle distribution function information. For instance, for a heat diffusion problem, the
temperature can be obtained as:
13
( ) ( )
0
, ,
b
i
i
T t f t
=
=

x x (55)
where b is the number of lattice connection vectors.
The equilibrium distribution function, appearing in Eq. 54, can be described as:
( ) ( ) , ,
eq
i i
f t wT t = x x (56)
It is found that invoking the Chapman-Enskog expansion, Eq. 54 recovers the macroscopic
energy equation (Eq. 52) for the following modeling of thermal diffusivity in case of d2q9
model:
( )
2
2
6
c
t o t = A (57)

Boundary Conditions
A boundary can be introduced to a lattice Boltzmann model by selecting the grid sites where the
boundary is to be set and evolving the particle distribution function in a different manner at these
sites. Two types of boundary conditions are described for the present context, one is prescribed
temperature and another one is the prescribed heat flux.












Fig. 2 Schematic plot of particle streaming and boundary condition for a d2q9 model
Following Fig. 2, a Dirichlet boundary condition (i.e. prescribed temperature) can be imposed on
the left wall (temperature T
0
), for example. To determine f
1
, f
5
, and f
8
, first Eq. 55 is invoked as:
( )
1 5 8 0 0 2 3 4 6 7
f f f T f f f f f f + + = + + + + + (58)
Now, applying the bounce back rule at the wall,
( )
1 3 3 1
eq eq
f f f f = and
( )
5 7 7 5
eq eq
f f f f = ,
one obtains:
( )
1 3
5 7
8 0 0 2 3 4 6 7
2 2
f f
f f
f T f f f f f f
=
=
= + + + + +
(59)
The boundary conditions for other walls follow the same procedure.
A Neumann boundary condition (i.e. prescribed heat flux) can be implemented by following a
conventional control volume based formulation. Applying an energy balance to the boundary
lattice j (refer to Fig. 2),
1
5
6
3
7 8
2
4
x
y
q
j
q
in
q
out
y A
j
j+1
14
( )
1 2 j t t t t
in out
j t A t
T
c dt dV q q dtdA
t

+ +A +A
c
=
c
} } } }
(60)
where
in j
q q = is the prescribed heat flux and
out
q is the flux leaving the control volume (which
can be described by Fouriers law). Hence, Eq. 60 becomes:
( )
( )
1
1
2
t t
j j
n n
j j j
t
k T T
y
c T T q dt
y

+A
+
+
(

A
= (
A
(

}
(61)
where superscripts n and n+1 represent time levels t andt t + A , respectively. Applying an
explicit scheme, the following linear algebric equation yields form Eq. 61:
( )
1
1 1 1
n n n
j j j j j j j j
a T a a T a T q
+
+ + +
= + + (62)
where 2
j
a c y t = A A and
1 j
a k y
+
= A . Equation 62 can be utilized to convert a prescribed heat
flux boundary condition, as an equivalent pseudo-isothermal boundary condition, which in turn,
can be implemented identical to the generalized formulation for incorporation of Dirichlet type
boundary condition, described as above.

Exercise: Prove that for a d2q9 model the thermal diffusivity assumes the following form,
( )
2
2
6
c
t o t = A


























15



Two Dimensional 9-bit Square Lattice Model

To recover the 9-bit LBE model on a square lattice, the Cartesian coordinate system is used, and
accordingly, ( ) can be set to
( )
,
m n
m n x y
= (A1)
where
x
and
y
are the x and y components of . The integral of the moments, defined by Eq.
(20), can be given as:
( )
( )
( )
2
1 1
,
2 2
2 1 1 2
2
2 1
2 2
2
m n
x m n y m n
eq
m n m n
x m n x y m n y m n
u I I u I I
I f d RT I I
RT RT
u I I u u I I u I I
RT

t
+
+ +
+ + + +
+
| |

= =
|
\ .
+ +

+
`

)
}
u

(A2)
where
2
, 2
m
m
I e d RT
,
, , ,
+

= =
}
is the m-th order moment of the weight function
2
e
,
on the real axis. Consequently, the third-order Gauss-Hermite quadrature formula is the
optimal choice to evaluate
m
I for the purpose of deriving a 9-bit LBE model:

3
1
m
m j j
j
I e ,
=
=

(A3)
The three abscissas of the quadrature are:

1 2 3
3 2, 0, 3 2 , , , = = = (A4)
and the corresponding weight coefficients are:

1 2 3
6, 2 3, 6 e t e t e t = = = (A5)
Then, the integral of the moment in Eq. (E2) becomes
( )
( ) ( )
( )
2
2 3
, ,
, 2
, 1
1
2
2
i j i j
i j i j
i j
I
RT RT
RT

ee
t
=



= + +
`

)

u u
u
(A6)
where
( ) ( )
,
, 2 ,
i j i j i j
RT , , = = , and the equilibrium distribution function can be identified
as:

( ) ( )
( )
2
2
, ,
, 2
1
2
2
i j i j i j eq
i j
f
RT RT
RT
ee

t



= + +
`

)
u u
u
(A7)
Employing the notations
( )
( ) ( ) ( )
( ) ( ) ( )
0, 0 0
cos 1 2 , sin 1 2 1, 2, 3, 4
2 cos 5 2 4 , sin 5 2 4 5, 6, 7, 8
c
c
o
o
o t o t o
o t t o t t o

= = ( (

+ + = ( (

e (A8)
16
and

4 9 2, 0
1 9 1, 2,... 1, 2, 3, 4
1 36 1,... 5, 6, 7, 8
i j
i j
w i j
i j
o
o
ee
o
t
o
= = =

= = = = =

= = =

(A9)
and substituting
( )
2 2
1
3 or 2 3
s
RT c c RT RT c , = = = = , the equilibrium distribution function
can be obtained as:

( ) ( )
2
2
, 2 4 2
3 9
3
1
2 2
eq
i j
f w
c c c
o o
o

= + +
`

)
e u e u
u
(A10)