See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/290158292

An introduction to Lattice-Boltzmann methods

Article · January 2013

DOI: 10.2174/9781608057160113030004

CITATIONS READS
9 547

3 authors:

Daniel Heubes Andreas Bartel

Bergische Universität Wuppertal Bergische Universität Wuppertal
5 PUBLICATIONS   31 CITATIONS    68 PUBLICATIONS   500 CITATIONS   

SEE PROFILE SEE PROFILE

Matthias Ehrhardt
Bergische Universität Wuppertal
169 PUBLICATIONS   1,597 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Stochastic Lie Group Methods View project

Mathematical Modelling in Energy Markets View project

All content following this page was uploaded by Matthias Ehrhardt on 07 September 2019.

The user has requested enhancement of the downloaded file.

 MA
CM

Bergische Universität Wuppertal

Fachbereich Mathematik und Naturwissenschaften

Institute of Mathematical Modelling, Analysis and Computational

Mathematics (IMACM)

Preprint BUW-IMACM 11/28

Daniel Heubes, Andreas Bartel and Matthias Ehrhardt

An Introduction to the Lattice Boltzmann

Method for Coupled Problems

December 2011

http://www.math.uni-wuppertal.de
1. INTRODUCTION

The circulation of blood in human blood vessels, the river flow at bridge pillars, the air flow passing a car
are examples of complex problems from fluid mechanics. Such problems can only be solved numerically.
Due to the today’s complexity of structures and models, such a fluid simulation must be very efficient. Often
computation is based on continuous models for mass and momentum conservation as the Navier-Stokes
equations. These models operate on the macroscopic level. More recently kinetic models and methods as the
Boltzmann equation and the lattice Boltzmann method compete against the standard concepts. Especially,
the lattice Boltzmann method is quite attractive, since an explicit update scheme is not only very efficient
in computation, but also fits to modern computer architecture as in graphical processing units. Its drawback
and its advantage is the microscopical description on the particle level. Each interaction, boundary condition
etc. have to be formulated on this basis.
This chapter is organized as follows: In the remaining part of this section, we shortly state the Navier-Stokes
equations and we discuss the basics of the Boltzmann equation, the collision integral and an approximation
of it. Our main focus in Section 2 lies on the derivation of the lattice Boltzmann method. Firstly, we study
the historical derivation via lattice gas automata. Next we show how the Boltzmann equation can be derived
from a particle ensemble description. Then we present the more recent derivation via discretization and
Gauß-Hermite quadrature. For the latter part, we discuss in particular the D3Q19 model. Finally, in this
section, we show with a Chapman-Enskøg expansion the relation to the Navier-Stokes equations. Section 3
is devoted to boundary conditions for the lattice Boltzmann method. Then, we discuss enhancements as
additional forces, temperature and an electric circuit coupling in Section 4. Eventually, this chapter ends
with conclusions (Section 5).

1.1. Navier-Stokes Equations

The macroscopic motion of a Newtonian fluid can be described by a balance equation of momenta. To this
end, let u = u (xx,t) denote the fluid velocity at space location x and time t, ρ = ρ(xx,t) the local mass density,
p = p(xx,t) the local pressure and η the dynamic viscosity. For a constant mass density ρ, the Navier-Stokes
equations read [34]:

∇ · u = 0, (1a)
∂u
ρ ∇ p + η∆uu.
+ ρ (uu · ∇ ) u = −∇ (1b)
∂t
The first equation (1a) states the incompressibility of the fluid (conservation of mass). The second equation
(1b) is the balance equation for momenta, where the fluid transport is driven by the pressure p and the stress
term η∆uu. This description gives the fundamental macroscopic fluid model, which we want to simulate. For
further details see [1].
Clearly, if additional body forces are present, this modifies the second balance equation and these additional
forces need to be added as a further addend to this equation.

1.2. Boltzmann Equation

In this section, we state the microscopic view on fluid modeling. Here, the Boltzmann equation without
additional forces is an evolution equation in time (t) for the single particle distribution f (xx, v,t) at space
position x with velocity v ∈ R3 :
∂ f (xx,t)
+ v · ∇ f (xx,t) = Q( f ) (2)
∂t
The value f (xx, v ,t) dxx dvv for all infinitesimal small dxx and dvv represents the probability to find a particle at
time t having position x and velocity v . In fact, this value is scaled with the constant particle mass m, thus
the mass density is given by the integral
Z
ρ(xx,t) = f (xx, v,t) dvv (3)
R3

in three dimensional space. The term Q( f ) in (2) expresses the change in the distribution f due to the
collision of particles with different velocities. It is referred to as collision term. In Section 1.3., we give
some more details of this quantity. Moreover, a rough derivation of the Boltzmann equation (2) is found in
Section 2.2.1. In addition to (3), the other macroscopic quantities as velocity u and temperature T can be
computed by higher moments of the single particle distribution:
Z
ρ(xx,t)uu(xx,t) = v f (xx, v ,t) dvv, (4)
R3
Z
3 kB T (xx,t) |vv − u (xx,t)|2
ρ(xx,t) = f (xx, v ,t) dvv, (5)
2 m 2
R3

with Boltzmann constant kB . Integrating the Boltzmann equation (2) over velocity space, that is, the zeroth
moment, results in the continuity equation
∂
ρ(xx,t) + ∇ · (ρ(xx,t)uu(xx,t)) = 0. (6)
∂t
It was used that the integral of the collision term vanishes [4]. Furthermore, the first moment is computed
by multiplying the Boltzmann equation by v and then integrating over velocity space; this yields
  Z
∂
∇ · (vv − u (xx,t)) (vv − u (xx,t))⊤ f (xx, v ,t) dvv,
(ρ(xx,t)uu(xx,t)) + ∇ · ρ(xx,t)uu(xx,t)uu(xx,t)⊤ = −∇ (7)
∂t
R3

where the integral on the right hand side is the momentum flux tensor, which exhibits the form [16]
Z
(vv − u (xx,t)) (vv − u (xx,t))⊤ f (xx, v ,t) dvv = p(xx,t)I − σ (xx,t),
R3

with identity matrix I, p(xx,t) denoting the pressure and

h i
∇u (xx,t))⊤ .
σ (xx,t) = ρ(xx,t)ν ∇ u (xx,t) + (∇

Here ν is the kinematic viscosity, and the complete term expresses the dynamic portion in the momentum
flux tensor. Hence, equation (7) becomes
∂  
(ρ(xx,t)uu(xx,t)) + ∇ · ρ(xx,t)uu(xx,t)uu(xx,t)⊤ = −∇
∇ p(xx,t) + ∇ · σ (xx,t). (8)
∂t
Now, if ρ is constant, (6) yields (1a) and (8) gives (1b), where the dynamic viscosity is η = ρν. Finally, we
notice that the temperature is derived by (5) from the distribution f and is not a variable of the system.

1.3. Collision Term and BGK Approximation

Next, we derive the Boltzmann integral expression for Q( f ). This needs several assumptions [27]: One
condition is that the particles interact only in two-particle collisions, this means we assume that interactions
involving more than two particles can be neglected. For all two-particle collisions we assume that they
appear locally in the sense that they take place at a single point x . A similar condition holds for time t,
it is assumed that the duration of a collision is negligible. Particles involved in a collision are assumed to
be uncorrelated, and the collision itself is modeled as an elastic collision, meaning that kinetic energy and
especially momentum are conserved. Now for a collision of two particles, we have velocities v , w before
and velocities v ′ , w ′ after the collision, the collision integral due to Boltzmann reads
Z Z  
Q( f ) = σ (Ω)|vv − w | f (xx, v ′ ,t) f (xx, w ′ ,t) − f (xx, v ,t) f (xx, w ,t) dΩ dw
w, (9)
R3 S 2

where σ (Ω) denotes the differential collision cross section and the inner integration is done over all pos-
sible solid angles Ω ∈ S2 . Moreover the post-collision velocities can be computed in dependence of the
pre-collision velocities and the impact angle.
The collision integral (9) possesses a rather complicated form. In the following, we motivate the approx-
imation introduced by Bhatnagar, Gross and Krook (BGK) [3] which simplifies the collision term and is
frequently used nowadays. Besides the BGK term, there exist also different alternatives for the collision
integral. One of them are multi-relaxation time models, e.g. [10]. In BGK schemes the collision integral
is replaced by a single time relaxation term. This term is chosen in such a way that it emulates specific
properties of the original collision integral, especially those which were necessary to derive equations (6)
and (8). It can be shown [4], that there exist five elementary invariants which read

ψ1 = 1, (ψ2 , ψ3 , ψ4 ) = v and ψ5 = |vv|2

and lead to vanishing integrals

Z
ψi Q( f ) dvv = 0, i = 1, . . . , 5.
R3

It even holds the following equivalence

Z 5
Q( f )ϕ(vv) dvv = 0 ⇔ ϕ(vv) = ∑ si ψi = α + β · v + γ|vv|2 ,
i=1
R3

with real scalars α, γ, si ∈ R (i = 1, . . . , 5) and a real vector β ∈ R3 . Since a Maxwellian distribution describes
the equilibrium state for the Boltzmann equation, the BGK approximation is modeled as a relaxation towards
a Maxwellian distribution
 3/2  
m −m
M(vv; ρ, u , T ) := ρ exp |vv − u |2 . (10)
2πkB T 2kB T

We remark the relation of the Boltzmann constant kB to the specific gas constant R = kB /m. Inserting the
Maxwellian distribution for f into (9), the integral over velocity space vanishes. In the BGK approximation
the collision term reads
1h i
Ω( f ) := − f (xx, v ,t) − f (eq) (xx, v ,t) ,
τc

where the equilibrium distribution f (eq) (xx, v ,t) is a Maxwellian one with space and time dependent ρ, u and
T . The values are chosen such that the aforementioned properties are emulated. Therefore they have to be
calculated in terms of f (xx, v ,t) as given in (3), (4) and (5).
The Boltzmann equation with BGK approximation for the collision term then reads

∂ f (xx, v ,t) 1h i
+ v · ∇ f (xx, v ,t) = − f (xx, v ,t) − f (eq) (xx, v ,t) . (11)
∂t τc
Like for the Boltzmann equation (2) the zeroth and first moment of the BGK approximated one (11) lead
also to equations (6) and (8).
2. LATTICE BOLTZMANN METHOD BASIC CONCEPTS

The current section revolves completely around the lattice Boltzmann method (LBM)

∆t h (eq)
i
fi (xx + c i ∆t,t + ∆t) − fi (xx,t) = − fi (xx,t) − fi (xx,t) . (12)
τc

We explain two ways which lead to this equation. On the one hand, in Section 2.1. we demonstrate the
historical development from very simple computational concepts to the LBM. The former could already be
implemented on the first computing devices in the 1960s, and the LBM can be implemented efficiently on
most recent computer architecture, for example a cluster of GPUs. On the other hand, in Section 2.2., we
derive the LBM from the Boltzmann equation. This connection to the Boltzmann equation was not known
directly with the appearance of the first lattice Boltzmann models in literature. Furthermore we also give a
short idea how the Boltzmann equation can be derived from an exact microscopical description of a fluid.
After these sections, which both end with the LBM equation (12), we verify the physical fundament of it
by a technique called Chapman-Enskøg expansion. We complete this section with providing a basis for
beginners who get in touch with the LBM for the first time by briefly explaining how the LBM works and
how it can be implemented.

2.1. Classical Approach

The historical development of the lattice Boltzmann method started with an adaption of lattice gas automata,
which itself can be understood as modified cellular automata. We track this development in the current
section by briefly describing cellular automata and lattice gas automata, afterwards we explain the adaption
of the lattice gas automata which lead to the lattice Boltzmann method. Finally we demonstrate the relation
of the lattice Boltzmann method to the discrete velocity Boltzmann equation.

2.1.1. Cellular Automata

As a universal concept for computation, John von Neumann [35] introduced the cellular automaton (CA) as
a discrete model. Imaging a CA with cells C formally placed on a D−dimensional lattice, we have C ⊂ ZD .
Each cell holds a state from a finite set of possible states S ⊂ Z [22]. At discrete time levels the states are
updated simultaneously by prescribed deterministic rules. These rules operate locally, that is, a new state of
any interior cell C i ∈ C depends on possibly all old states of neighboring cells including the cell itself. Due
to the regular structure of the lattice, we can refer to the neighboring cells by local coordinates N , such that
the set of neighbors is given by C i + N . Boundary cells are updated due to some appropriate update rules.
As an illustrating example, we consider a one-dimensional CA consisting of 400 cells with only two possible
states for each cell: S = {0, 1}. States of cell C1 and C400 are set to 0 and are not updated. The updating
rule for cells in the interior shall only depend on the states of the two adjacent cells and the cell itself, i.e.,
N = {−1, 0, 1}. Hence, any update rule r must assign a new configuration to a local configuration (which
are the corresponding three states), that is,

r(0, 0, 0) = r0 , r(0, 0, 1) = r1 , r(0, 1, 0) = r2 , r(0, 1, 1) = r3 ,

r(1, 0, 0) = r4 , r(1, 0, 1) = r5 , r(1, 1, 0) = r6 , r(1, 1, 1) = r7

with corresponding values ri ∈ S . This amounts to 28 = 256 different possible update rules for this one-
dimensional CA. Then the state si of the ith interior cell is updated to s′i :

s′i = r(si−1 , si , si+1 ), i = 2, . . . , 399.

2.1.2. Lattice Gas Automata

Lattice gas automata (LGA) are derived from classical CA by some modifications, which simplify the con-
struction and application of automata to given physical processes [36]. LGA are capable to simulate fluid
flows successfully, since there exist corresponding automata which lead to the Navier-Stokes equations in
the macroscopic limit [38]. And despite of their relative simple nature, they can be applied in less simple
themes such as for instance in the simulation of flows through porous media [7].
The first introduced LGA, known as HPP due to Hardy, de Pazzis and Pomeau [17], was proposed as a
new technique for the numerical study of the Navier-Stokes equations. Here the direct integration of the
Navier-Stokes equations is replaced by the simulation of a very simple microscopic system: particles of the
same mass are allowed to move on a regular lattice, and local collision rules are introduced on the nodes
which conserve the number of particles and momentum. However the momentum flux tensor in the resulting
equations has not the correct form such as required by the Navier-Stokes equations.
The first LGA which overcame this problem were introduced by Frisch, Hasslacher and Pomeau in 1986
[14], the so called FHP model. The major difference between the HPP and FHP model consists in the fact
that the HPP used a square lattice and in the FHP model a triangular lattice is used. The use of a triangular
lattice leads to a specific isotropic tensor of rank four which is necessary to achieve the correct form of the
momentum flux tensor in sense of the Navier-Stokes equations, see [36] for details. The first step towards
a formal description of LGA is the definition of a regular lattice. Depending on the lattice, a number of k
different lattice vectors d i (i = 1, . . . , k) are introduced which connect nearest neighbors. For the FHP model
(more precisely, the FHP-I model) in two dimensions, with k = 6, they read
 π   π 
d i = cos i , sin i (i = 1, . . . , 6). (13)
3 3
The lattice vectors describe the possible particle velocities (and thus possible the momenta). Now at certain
discrete positions r , lattice nodes are located. Each lattice node consists of k cells, which describe whether
a particle with corresponding velocity is present or not. Thus we have binary states ni (i = 1, . . . , k) with
(
0, cell i is not occupied by a particle,
ni (rr ,t) =
1, cell i is occupied by a particle.

Here, t indicates the time. Due to an exclusion principle each cell is either occupied by one particle or holds
no particle at all. All states are updated at discrete time levels simultaneously by an evaluation of collision
rules Ψi and a streaming:

⊤
ni (rr + d i ,t + ∆t) = ni (rr ,t) + Ψi (nn(rr ,t)) with n (rr ,t) = n1 (rr ,t), n2 (rr ,t), . . . , nk (rr ,t) . (14)

The streaming is encoded in the location of the left-hand side r + d i . The collision rules Ψi (nn(rr ,t)) should
be chosen such that they conserve the number of particles and momentum. We illustrate the collision rules
of the FHP model in Fig. (1). For the two particle head on collision there are two possible outcomes which
conserve the number of particles and momentum and one has to take a choice randomly. Configurations not
listed in Fig. (1) are not affected due to collision, i.e., only two and three particle collisions are considered.

2.1.3. Translation to Continuous Variables

In the following, our main focus is on the historical development of LBM from LGA. Using a Chapman-
Enskøg expansion [5] up to second order terms in velocity u , the original FHP model was shown to yield
the continuous mathematical model [14]:
∂ρ
+ ∇ · (ρuu) = 0,
∂t
∂ ρuα ∂   ∂ ∂
+∑ γ(ρ)ρuα uβ = − p + η∆uα + ξ ∇ · u)
(∇
∂t β
∂ xβ ∂ xα ∂ xα
 (a) Two particle collisions; conservation of mass and momentum
yields two possible configurations after collision. Thus a random
choice is applied (with probability one half).

(b) Three particle collisions.

Fig. (1): Collision rules used in the FHP model. Empty circles indicate a state of 0, whereas solid circles are
occupied cells of state 1.

with γ(ρ) = ρ−3ρ−6 . This result differs from the Navier-Stokes equations mainly by the presence of the term
t
γ(ρ). If this term is a constant, for 0 < ρ < 3, it can be absorbed in a rescaled time t γ(ρ) . Thus in
the incompressible limit, where all density variations are neglected except density variations in the pressure
term, the latter equations become

∇ · u = 0,
∂u (15)
ρ ∇ p + η∆uu.
+ ργ(ρ) (uu · ∇ ) u = −∇
∂t
Now, the generalized substantial derivative
∂u
+ γ · (uu · ∇ )uu
∂t
fulfills a Galilean invariance if and only if γ ≡ 1 (e.g. [36]), which is impossible in the original FHP model.
This lack of Galilean invariance in this FHP model is due to the exclusion principle in combination with
a finite set of allowed velocities. Though, there exist modifications of this model, in which the Galilean
invariance holds [11].
What we have not yet dealt with is the question how the macroscopic quantities in (15) can be derived from
LGA. This is specially demanding, since the quantities appearing in (14) are only binary ones. Both, density
and velocity are computed from local average populations Ni by

ρ(xx,t) = ∑ Ni (xx,t), ρ(xx,t)uu(xx,t) = ∑ d i Ni (xx,t).

i i

The notation Ni indicates an averaging that corresponds to the lattice vector d i . One procedure to compute Ni
is called coarse graining; here one calculates the mean values over large subregions. Another possibility is
ensemble averaging. In this approach one simulates the evolution of LGA many times with different initial
configurations, and computes the average over all simulations. It could also be used in combination with
coarse graining. Although the evolution equation (14) for LGA is quite simple, one needs to spend a high
computational effort to obtain reasonable macroscopic quantities by the averaging procedure.
Two years after the introduction of the FHP McNamara and Zanetti [24] suggested to change the binary
populations of LGA into real values from the interval [0, 1] representing the average values and thus to
avoid the averaging process. Nowadays, this paper is seen as the origin of the lattice Boltzmann method,
even though similar probabilistic views of LGA were common in the literature before [13, 37]. Due to this
suggestions, not only the ni are replaced by continuous variables fi , but also an arithmetic treatment of the
collision term replaces the boolean collision term. The BGK model, see Section 1.3., was independently
introduced for lattices by [26] and [6] in 1992. Using the BGK model for collision, the LGA model changed
into the lattice Boltzmann equation
1h (eq)
i
fi (xx + c i ,t + ∆t) − fi (xx,t) = − fi (xx,t) − fi (xx,t) , (16)
τc
with lattice vectors c i and an equilibrium distribution
 
(eq) 9 3
fi (xx,t) = wi ρ(xx,t) 1 + 3(cci · u (xx,t) + (cci · u (xx,t)2 − |uu|2 ,
2 2
where the weights wi depend on the chosen lattice and the macroscopic quantities are computed by
ρ(xx,t) = ∑ fi (xx,t), ρ(xx,t)uu(xx,t) = ∑ ci fi (xx,t).
i i
The lattice Boltzmann method closes not only the lack of Galilean invariance, see above, but also another
drawback of LGA not considered here which is statistical noise [9]. We will consider different lattices
below.

2.1.4. Relation to the Discrete Velocity Boltzmann Equation and Standard Discretizations

The lattice Boltzmann method (16) uses a discrete velocity space, say:

V = c i ∈ R3 : i = 0, . . . , nv .
Thus it can be understood as a certain discretization of the Boltzmann equation (2). Applying this velocity
discretization directly to the Boltzmann equation gives a discrete velocity Boltzmann equation [32]:
∂ fi (xx,t) 1h (eq)
i
+ c i · ∇ fi (xx,t) = − fi (xx,t) − fi (xx,t) , i = 0, . . . , nv , (17)
∂t τc
(eq)
where the discrete velocity space V replaces R3 , such that f (xx, c i ,t) → fi (xx,t) and f (eq) (xx, c i ,t) → fi (xx,t).
The set V has to be interpreted here as a set of vectors defining the lattice. Frequently used lattices in LBM
are the D2Q7 and D2Q9 in two dimensions and D3Q15, D3Q19, as well as D3Q27 in three dimensions.
In the notation DxQy, which goes back to [26], the value x denotes the spatial dimension and the value
y(= nv + 1) denotes the number of lattice velocities. D2Q7 is the lattice used in the FHP model; the corre-
sponding velocities are given in (13) and depicted as lattice vectors in Fig. (2)(a). For given reference lattice
velocity c, the lattice vectors of the more usual D2Q9 lattice are

! √
!
cos π2 (i − 1) 2 cos π2 ( j − 12 ) i = 1, 2, 3, 4
c0 = 0, ci = c
, cj = c √
,
sin π2 (i − 1) π
2 sin 2 ( j − 2 )1 j = 5, 6, 7, 8

and they are depicted in Fig. (2)(b). For the three dimensional case, we consider a cubic lattice with cubes of
edge length 1. Connecting
√ the nearest neighbors,
√ we achieve six lattice velocities of length 1, twelve velocity
vectors of length 2 and eight of length 3. See Fig. (3) for an illustration. In addition, we have the zero
velocity, which belongs to any discrete velocity set. Different combinations of these lattice velocities yield
the well-known√ three dimensional models: We obtain√ the D3Q15 model by taking the velocity vectors of
length 1 and 3. Extending this model by the twelve 2-velocity vectors, we get√the D3Q27 model. And
the popular D3Q19 model is constructed of lattice velocities with length 1 and 2, see Fig (4), i.e., the
generalized lattice vectors for given reference lattice velocity c are:
c 0 = (0, 0, 0), c 1 = (1, 0, 0)c, c 2 = (0, 1, 0)c c 3 = (−1, 0, 0)c,
c 4 = (0, −1, 0)c, c 5 = (0, 0, 1)c, c 6 = (0, 0, −1)c c 7 = (1, 1, 0)c,
c 8 = (−1, 1, 0)c, c 9 = (−1, −1, 0)c, c 10 = (1, −1, 0)c, c 11 = (1, 0, 1)c, (18)
c 12 = (−1, 0, 1)c, c 13 = (−1, 0, −1)c, c 14 = (1, 0, −1)c, c 15 = (0, 1, 1)c,
c 16 = (0, −1, 1)c, c 17 = (0, −1, −1)c, c 18 = (0, 1, −1)c.
 6 2 5

3 1

7 4 8
(a) Lattice vectors of D2Q7. (b) Lattice vectors of D2Q9.

Fig. (2): Standard two dimensional lattice models.

(a) (b) (c)

Fig. (3): Lattice velocities for cubic lattice in three dimensions.

15

5
12 11

16
2
8 7

3 1

9 10
4
18

z y 13 14
6

x
17

Fig. (4): The lattice vectors of the D3Q19 model.

We come now to the space and time discretization of the discrete velocity Boltzmann equation (17). This
discretization must be suitable for the velocity space, that is, for all velocity vectors c i at an interior grid
point x j , the sum x j + c i ∆t describes a further grid point, where ∆t denotes the time step. Hence, we have a
link between space discretization and time discretization. Then the advection term c i · ∇ fi (xx,t) in (17) can
be discretized by finite differences as
fi (xx + c i ∆t,t + ∆t) − fi (xx,t + ∆t)
c i · ∇ fi (xx,t) = + O(∆t).
∆t
Using the forward Euler to approximate the time derivative, the fully discretized Boltzmann equation with
BGK approximation gives the lattice Boltzmann equation (16). This completes the derivation.

2.2. Theoretical Approach

In the previous section we have outlined the historical development of the lattice Boltzmann method. In
contrast the current section explains roughly the theoretical background. First, we give an idea of how the
Boltzmann equation can be derived from a microscopical description. Then we explain the derivation of
LBM via suitable discretizations.

2.2.1. Sketched Derivation of the Boltzmann Equation

For a microscopical view of a fluid, we consider a system consisting of N particles (e.g. molecules) with
the mass m each. Let the Hamiltonian H = H (t, q , p ) expresses the energy of the system as a function of
generalized position coordinates q = q (t) = (qq1 , . . . , q N ) and generalized momenta p = p (t) = (pp1 , . . . , p N ),
then the evolution of that dynamic system can be described by the Hamilton equations [4]:
∂H ∂H
q̇qk = , ṗpk = − , k = 1, . . . , N. (19)
∂ pk ∂ qk
A direct consideration of all N particles is being ruled out, because in practical all interesting cases the
number of particles is too huge O(10k ), k > 20. Fortunately, in general one is not directly interested in
the movement of each specific particle, instead it is often enough to know how many particles with specific
velocities are in average present. Therefore, it offers to consider a N-particle distribution function PN , where

PN (qq1 , p1 , q2 , p2 , . . . , qN , pN ,t) dqq1 . . . dqqN dpp1 . . . dppN

gives the probability to find this N particle ensemble with ith particle in an infinitesimal small space volume
dqqi around q i with momentum in the infinitesimal small volume dppi around p i at time t. The evolution of
this probability distribution is expressed by the Liouville equation (see e.g. [4])
N N
dPN ∂ PN ∂ PN ∂ PN
= + ∑ q̇qk + ∑ ṗpk = 0, (20)
dt ∂t k=1 ∂ q k k=1 ∂ p k

where the time derivatives q̇q and ṗp are given by the Hamilton equations (19). The Liouville equation can
be transformed into an equivalent system of coupled equations, which successively involves more and more
particles. To this end, we introduce s-particle distribution functions by integration with respect to all particles
from s + 1 till N:
Z
(s)
PN (qq1 , p 1 , q 2 , p 2 , . . . , q s , p s ,t) = PN (qq1 , p 1 , q 2 , p 2 , . . . , q N , p N ,t) dqqs+1 . . . dqqN dpps+1 . . . dppN

s = 1, . . . , N − 1. Moreover, we split the force terms for the kth particle (Hamilton equation)
N
∂H
ṗpk = − = k k,ex + ∑ k k, j
∂ qk j=1, j6=k
into external forces k k,ex and forces coming from interactions with the jth particle k k, j . Now, we plan to
successively remove the dependency on particles by integration, to this end we write the substitution of
these split terms in (20) as follows:
" ! #
N
∂ PN ∂ PN N−1 N−1
∂ PN ∂ PN ∂ PN N−1 ∂ PN
+ ∑ q̇qk +∑ k k,ex + ∑ k k, j + k k,N + k N,ex + ∑ k N, j = 0. (21)
∂t k=1 ∂ q k k=1 j=1, j6=k ∂ p k ∂ p k ∂ p N j=1 ∂ pN

The integration of (21) with respect to the coordinates and momenta of the last N − 1, N − 2, . . ., 1 particles
yields the following equivalent system for s-particle distribution functions:
(1) (1) (1) Z (2)
∂ PN ∂P (1) ∂ P ∂ PN
+ q̇q1 · N + k 1 · N = (1 − N) k 1,2 · dqq2 dpp2
∂t ∂ q1 ∂ p1 ∂ p1
..
.
(s)
s s (s) (s) Z s (s+1)
∂ PN ∂P (s) ∂ P ∂ PN
+ ∑ q̇qk · N + ∑ k k · N = (s − N) ∑ kk,s+1 · dqqs+1 dpps+1 (22)
∂t k=1 ∂ q k k=1 ∂ pk k=1 ∂ pk
..
.
(N−1) N−1 (N−1) N−1 (N−1) Z N−1 (N)
∂ PN ∂ PN (N−1) ∂ PN ∂ PN
∂t
+ ∑ q̇qk ·
∂ qk
+ ∑ kk ·
∂ pk
=− ∑ k k,N ·
∂ pk
dqqN dppN
k=1 k=1 k=1

with abbreviations for the forces

s
k (s)
k = k k,ex + ∑ k k, j .
j=1, j6=k

E.g. for the last equation of (22) the derivatives w.r.t. the Nth particle in (21) balance and therefore these
terms drop out of the equation. System (22) is the so-called BBGKY-hierarchy named after Bogoliubov,
Born, Green, Kirkwood and Yvon (see e.g. [4]). The labeling of the particles was arbitrary, but the system
is given for the specific particle distribution functions. For our purpose, it would be enough to know if
there are particles at given positions, instead of having also the information which specific particles they are.
Therefore we scale the specific particle distribution functions to obtain general particle distribution functions
(s)
P̃N :

(s) N! (s)
P̃N = P .
(N − s)! N
By this, from (22) we get
(s) s s (s) (s) Z s (s+1)
∂ P̃N ∂ P̃ (s) ∂ P̃ ∂ P̃N
+ ∑ q̇qk · N + ∑ k k · N = − ∑ k k,s+1 · dqqs+1 dpps+1 , s = 1, . . . , N − 1. (23)
∂t k=1 ∂ q k k=1 ∂ pk k=1 ∂ pk

The structure of the BBGKY hierarchy is simple, the first equation (i.e., s = 1) is an equation for the single
particle distribution function and it is coupled to the two particle distribution function by the integral on
the right hand side. Then the evolution equation for the two particle distribution (s = 2) couples to the three
particle distribution function, and so on. The full system (23) is equivalent to the Liouville equation, however
we pay special attention to the first equation, since it will yield the Boltzmann equation: its term on the right
hand side covers all the coupling to the other particles and thus describes collisions. An exact modeling of
collision as (23) has to incorporate a correlation of particle velocities due to preceding collisions. Now, the
Boltzmann equation uses this first equation with a certain approximation of the right-hand side term Q( f )
(1) (1)
and the scaled unknown f (xx, v ,t) := mP̃N (qq1 , p 1 ,t). This reads for the external force density K = k k :

∂ f (xx, v ,t) ∂ f (xx, v ,t) ∂ f (xx, v ,t)

+v· +K · = Q( f ).
∂t ∂x ∂v
2.2.2. Time Discretization

For the sake of simplicity we consider the Boltzmann equation without external forces. Later, we incorporate
external forces to the lattice Boltzmann method again in Section 4.2. Also, we substitute the collision integral
by the BGK approximation in the following. Hence, we have
∂ f (xx, v ,t) 1h i
+ v · ∇ f (xx, v ,t) = − f (xx, v ,t) − f (eq) (xx, v ,t) , (24)
∂t τc
with single time relaxation parameter τc . For the discretization, we follow the strategy from He and Luo
[19]. By introducing the advective differential operator dtd = ∂t
∂
+ v · ∇ , equation (24) can formally be written
as a linear first order ordinary differential equation (ODE):
d 1 1
f (xx, v ,t) + f (xx, v ,t) = f (eq) (xx, v ,t).
dt τc τc
Formally, its solution at t + ∆t (with ∆t ≥ 0) reads
    Z∆t  
∆t 1 ∆t s
f (xx + v ∆t, v ,t + ∆t) = exp − f (xx, v ,t) + exp − exp f (eq) (xx + v s, v ,t + s) ds.
τc τc τc τc
0

A Taylor expansion for the integrand (around s = 0)

 
s
exp f (eq) (xx + v s, v ,t + s) = f (eq) (x,t) + O(s),
τc

Taylor expansion for the exponential functions as well as neglecting all terms of order O(∆t 2 ) yields the
following approximation of the Boltzmann equation as a time discrete formulation:
∆t h i
f (xx + v ∆t, v ,t + ∆t) − f (xx, v ,t) = − f (xx, v ,t) − f (eq) (xx, v ,t) . (25)
τc
This equation already resembles the lattice Boltzmann method, but with continuous velocity space. By the
lattice construction, a velocity space discretization will imply a space discretization. This is our next topic.

2.2.3. Finding Appropriate Discrete Velocity Spaces


A discrete velocity space V = c i ∈ R3 : i = 0, . . . , nv is also used to compute the momenta of the distri-
bution f for certain polynomials ψ(vv), i.e.,
Z
I= ψ(vv) f (xx, v ,t) dvv
R3

via a quadrature rule. Now, an appropriate quadrature rule is required to deliver the necessary momenta
sufficiently accurate. In performing a Chapman-Enskøg expansion to recover the Navier-Stokes equations,
one derives that the first two moments of f (0) and f (1) have to be given exactly. The zeroth order term f (0)
is equal to the Maxwellian distribution f (eq) , the first order term is given by
!
(1) ∂ (1) (1)
f = −τc +v ·∇ f (eq)
∂t

and thus also related to the equilibrium distribution. Finally, one derives from the Chapman-Enskøg expan-
sion that the integral I is also given by
Z
I= ψ(vv) f (eq) (xx, v ,t) dvv (26)
R3
for all polynomials ψ(vv) = vm n p
x vy vz of up to third degree (i.e., m + n + p ≤ 3). These integrals are needed for
the computation of the macroscopic variables ρ and u [19].
Now, a Taylor expansion of the Maxwellian distribution f (eq) around u = 0 up to second order yields the
following approximation
 3/2   
1 |vv|2 1 1 2 1 2
f (eq) ≈ ρ exp − 1+ (vv · u ) + (v
v · u ) − |u
u | . (27)
2πRT 2RT RT 2(RT )2 2RT

Consequently, we have to evaluate integrals of type

Z  
|vv|2
Ib = exp − Ψ(vv) dvv
2RT
R3

to recover the Navier-Stokes equations, where Ψ(vv) = vm n p

x vy vz is a polynomial of degree five at most, i.e.,
m + n + p ≤ 5. In order to evaluate the integral numerically but without further approximation error, we
have to use a quadrature formula which is exact up to this degree. Using a quadrature formula leads to
a discrete equilibrium distribution, which contains the weights and lattice velocities. He and Luo [19]
separated the integral into products of one dimensional integrals, where they used 3-point Gauß-Hermite
quadrature formula, thus they obtained the D3Q27 lattice Boltzmann model. Our purpose here is to introduce
a quadrature formula which yields the D3Q19 model. Using the approximation (27) for the integral (26), we
obtain
" 
ρ |uu|2 1
I≈ 1 − Im,n,p + (ux Im+1,n,p + uy Im,n+1,p + uz Im,n,p+1 )
(2πRT )3/2 2RT RT
1
+ u2 Im+2,n,p + u2y Im,n+2,p + u2z Im,n,p+2 (28)
2(RT )2 x
#
+2ux uy Im+1,n+1,p + 2ux uz Im+1,n,p+1 + 2uy uz Im,n+1,p+1 )

with
Z  
|vv|2
Im,n,p := exp − vm n p
x vy vz dv
v.
2RT
R3

Approximating the integrals Im,n,p by a Gauß-Hermite type formula on the velocity set of D3Q19, we have
Z

m+n+p+3 m+n+p+3
Im,n,p = (2RT ) 2 exp −|ζζ |2 ζxm ζyn ζzp dζζ ≈ (2RT ) 2 ∑ wi, j,k ζim ζ jn ζkp =: Qm,n,p , (29)
(i, j,k)∈J
R3
p
with nodes ζ1/3 = ∓ 3/2, ζ2 = 0, index set

J = {(1, 1, 2); (1, 2, 1); (1, 2, 2); (1, 2, 3); (1, 3, 2); (2, 1, 1); (2, 1, 2); (2, 1, 3); (2, 2, 1); (2, 2, 2);
(2, 2, 3); (2, 3, 1); (2, 3, 2); (2, 3, 3); (3, 1, 2); (3, 2, 1); (3, 2, 2); (3, 2, 3); (3, 3, 2)} ,

and weights

1

3 for (i, j, k) = (2, 2, 2)
1 1

wi, j,k = 18 for (i, j, k) ∈ (1, 2, 2) ; (2, 1, 2) ; (2, 2, 1) ; (2, 2, 3) ; (2, 3, 2) ; (3, 2, 2)
π 3/2 
1
36 otherwise.
 α 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
i 2 3 2 1 2 2 2 3 1 1 3 3 1 1 3 2 2 2 2
j 2 2 3 2 1 2 2 3 3 1 1 2 2 2 2 3 1 1 3
k 2 2 2 2 2 3 1 2 2 2 2 3 3 1 1 3 3 1 1

Table 1: Index mapping for α.

The quadrature Qm,n,p yields an exact integration for Im,n,p for m + n + p ≤ 5. Thus, the substitution of these
Qm,n,p in (28) gives

18  
1 1 2 |uu|2 3


I= ∑ ρwα ξim ξ jn ξkp 1+ u ξ
(u · α ) + (u
u · ξ α ) − + O |u
u |
α=0 RT 2(RT )2 2RT
√ 1
with ξi = 2RT ζi , wα = π 3/2 wi, j,k and ξ α = (ξi , ξ j , ξk )⊤ corresponding to the index mapping defined by
Table 1. For the discrete equilibrium distribution it follows from (27):
( )
(eq) 3 (ccα · u ) 9 (ccα · u )2 3|uu|2
fα = ρwα 1 + + − , (30)
c2 2c4 2c2
√ √
with c = k 2RT ζ1 k = 3RT and corresponding velocities c α = ξ α , which are exactly the velocities given
in (18). In summary, we have
Z 18 18  
(eq) (eq) 3/2 |ccα |2
ψ(vv) f (xx, v ,t) dvv = ∑ ψ(ccα ) fα (xx,t) = ∑ ψ(ccα )(2πRT ) exp wα f (eq) (xx, c α ,t),
α=0 α=0 2RT
R3

and we conclude that moments of f (xx, v ,t) can be computed by

Z 18  
|ccα |2
ψ(vv) f (xx, v ,t) dvv = ∑ ψ(ccα ) fα (xx,t) with fα (xx,t) := (2πRT )3/2 exp wα f (xx, c α ,t).
α=0 2RT
R3

Since the discrete velocity set is determined, equation (25) can be fully discretized, and it follows the system

∆t h (eq)
i
fα (xx + c α ∆t,t + ∆t) − fα (xx,t) = − fα (xx,t) − fα (xx,t) , α = 0, . . . , nv . (31)
τc
In particular, it follows that the macroscopic quantities mass density ρ and velocity u are computed as
18 18
ρ(xx,t) = ∑ fα (xx,t), ρ(xx,t)uu(xx,t) = ∑ c α fα (xx,t). (32)
α=0 α=0

This completes the theoretical derivation.

2.3. Verification of physical fundament: Chapman-Enskøg-Expansion

In this section, we draw the relation of the lattice Boltzmann equation (12) to the Navier-Stokes equations
for incompressible fluids. To this end, a Taylor expansion in ∆t of the lattice Boltzmann equation (31) up to
second order terms is performed. This gives
   2
∂ 1 2 ∂
fα (xx + c α ∆t,t + ∆t) − fα (xx,t) = ∆t + c α · ∇ fα (xx,t) + ∆t + c α · ∇ fα (xx,t) + O(∆t 3 ).
∂t 2 ∂t
Neglecting the third order terms in this expansion and inserting the result into (31) gives:
   2
∂ 1 2 ∂ ∆t h (eq)
i
∆t + α·
c ∇ fα (x ,t) + ∆t
x + α·
c ∇ fα (xx,t) = − fα (xx,t) − fα (xx,t) . (33)
∂t 2 ∂t τc

Next, we expand the distribution fα (xx,t) around the equilibrium distribution:

(0) (1) (2)
fα (xx,t) = fα (xx,t) + ε fα (xx,t) + ε 2 fα (xx,t) + O(ε 3 ), (34)
(0) (eq)
where the zeroth order fα (xx,t) represents the equilibrium distribution fα (xx,t) and ε is a small parameter.
The macroscopic density and velocity compute as follows
nv nv
(0)
ρ(xx,t) = ∑ fα (xx,t) = ∑ fα (xx,t),
α=0 α=0
nv nv
(35)
(0)
ρ(xx,t)uu(xx,t) = ∑ c α fα (xx,t) = ∑ c α fα (xx,t).
α=0 α=0

Hence, we can conclude

nv nv
(k) (k)
∑ fα (xx,t) = 0, ∑ cα fα (xx,t) = 0. (36)
α=0 α=0

for any k ≥ 1. The differential operator with respect to time is expanded as well

∂ ∂ (1) ∂ (2)
=ε + ε2 + O(ε 3 ), (37)
∂t ∂t ∂t
whereas the spatial differential operator obeys

∇(1) .
∇ = ε∇ (38)

Substituting the expansions (34), (37) and (38) in (33) and splitting terms with respect to the order of ε yield
the first order equation
!
∂ (1) (1) (0) 1 (1)
ε + cα · ∇ fα (xx,t) = − ε fα (xx,t), (39)
∂t τc

and the second order equation

"  ! #
2 ∆t ∂ (1) (1) ∂ (2) (0) 1 (2)
ε 1− + cα · ∇ (1)
fα (xx,t) + fα (xx,t) = − ε 2 fα (xx,t) (40)
2τc ∂t ∂t τc

which was simplified by using (39). Summing up (39) for α = 0, . . . , nv yields

∂ (1)
ρ(xx,t) + ∇(1) · (ρ(xx,t)uu(xx,t)) = 0 (41)
∂t
by (35) and (36). This corresponds to a zeroth moment. Correspondingly, the summation of (39) after
multiplying with cα results in

∂ (1)
(ρ(xx,t)uu(xx,t)) + ∇ (1) · P(0) (xx,t) = 0, (42)
∂t
with
nv
(0)
P(0) (xx,t) = ∑ c α c α⊤ fα (xx,t).
α=0
This corresponds to a first moment. Due to (34), also the momentum flux tensor
nv
P(xx,t) = ∑ c α c α⊤ fα (xx,t)
α=0

exhibits an expansion in terms of ε

P(xx,t) = P(0) (xx,t) + εP(1) (xx,t) + ε 2 P(2) (xx,t) + O(ε 3 ).

The known form of f (0) (xx,t) implies

P(0) (xx,t) = ρ(xx,t)uu(xx,t)uu⊤ (xx,t) + p(xx,t)I

with identity matrix I and the pressure in lattice Boltzmann methods

c2
p(xx,t) = ρ(xx,t). (43)
3
Hence, equation (42) becomes

∂ (1)  
(ρ(xx,t)uu(xx,t)) + ∇ (1) · ρ(xx,t)uu(xx,t)uu(xx,t)⊤ = −∇
∇(1) p(xx,t). (44)
∂t

Accordingly, the zeroth and first moment of (40) have to be computed. Like above, the zeroth moment is
readily given as

∂ (2)
ρ(xx,t) = 0. (45)
∂t
Summation of (40) after multiplying with c α gives the terms of the momentum equation of second order. It
follows
  
∆t ∂ (2)
∇ (1) · 1 − P(1) (xx,t) + (ρ(xx,t)uu(xx,t)) = 0. (46)
2τc ∂t

The crucial part is now the computation of the first order momentum flux tensor
nv
(1)
P(1) (xx,t) = ∑ c α c α⊤ fα (xx,t).
α=0

Using (39) makes the computation straightforward, a detailed calculation is given for instance in [36]. By
using the previous results of lower order, equation (43) and neglecting terms of O(|uu|2 ) one obtains
 h i⊤ 
(1) c2 (1) (1)
P (xx,t) = −τc ρ(xx,t) ∇ u (xx,t) + ∇ u (xx,t) .
3

Eventually, by substituting the latter equation into (46), the first moment of (40) is given by
    ⊤ 
∂ (2) 2τc − ∆t 2 (1) (1) (1)
(ρ(xx,t)uu(xx,t)) = ∇ u (xx,t) + ρ(xx,t) ∇ u (xx,t)
c ∇ · ρ(xx,t)∇ . (47)
∂t 6

Adding the equations for the zeroth moments, (41) and (45), leads to an approximation up to second order
of the continuity equation

∂
ρ(xx,t) + ∇ · (ρ(xx,t)uu(xx,t)) + O(ε 3 ) = 0.
∂t
Similarly, by addition of the first moments, both first and second order in ε, i.e., (44) and (47), yields
∂  
(ρ(xx,t)uu(xx,t)) + ∇ · ρ(xx,t)uu(xx,t)uu(xx,t)⊤ + O(ε 3 )
∂t  
2τc − ∆t 2 h i
= −∇ ∇ p(xx,t) + c ∇ · ρ(xx,t)∇ ∇u (xx,t))⊤ .
∇u (xx,t) + ρ(xx,t)(∇
6
Under the assumption that density fluctuations are negligible these equations coincide with the Navier-Stokes
equations for incompressible fluids (1). We derive that the kinematic viscosity ν in the lattice Boltzmann
method is tuned by varying the relaxation parameter τc as:
η 2τc − ∆t 2
ν= = c .
ρ 6
This completes the link of LBM to the Navier-Stokes equations.

2.4. LBM: How It Works

Basically, the lattice Boltzmann method can be divided in two separate steps, one called the streaming or
propagation step and the other called the collision step. These terms are frequently used when dealing with
the lattice Boltzmann method, and we will explain them here now.
For simplification, we call each of the fi ’s in (12) a population of the single particle distribution. Note that
each population is linked with a direction according to the lattice vectors. Fig. (5) shows the 8 populations
for one point in a D2Q9 lattice, where the different lengths of the arrows symbolize the different magnitudes
of the populations. The magnitude of the population corresponding to the rest velocity is omitted in this
picture. Writing the lattice Boltzmann equation (12) as
∆t h (eq)
i
fi (xx + c i ∆t,t + ∆t) = fi (xx,t) − fi (xx, ,t) − fi (xx,t)
τc
separates the two steps. An evaluation of the right hand side, which is done at a fixed time point t yields new
values for the populations at any x:
∆t h (eq)
i
f˜i (xx,t) = fi (xx,t) − fi (xx,t) − fi (xx,t) .
τc
This process is called collision. It is illustrated in the first line of Fig. (5), the dashed arrows correspond
to the populations of f˜i (xx,t). The so-called streaming or propagation step is described by the following
equation:

fi (xx + c i ∆t,t + ∆t) = f˜i (xx,t).

It is noteworthy that the time evolution is exclusively specified by this step. In the streaming step, the post-
collision population at position x corresponding to direction c i becomes the post-streaming population at
position x + c i ∆t, see also the second line of Fig. (5).
We finish this section with a pseudocode algorithm for the lattice Boltzmann method on D3Q19 with BGK
collision term in Algorithm 1. For sake of simplicity the constant lattice velocity c is set to 1. The algorithm
illustrates the principal idea. For an implementation it is necessary to restrict the fluid to a bounded domain,
which is the topic of the following section.

3. STANDARD BOUNDARY CONDITIONS

For an implementation of the LBM, a restriction to a bounded domain is necessary. An interior lattice point
has the full number of neighbors. For these points the main concept of the lattice Boltzmann method holds.
Initialization:;
Set time horizon Tend for simulation ;
Define lattice velocities:;
c0 = [ 0; 0; 0], c1 = [ 1; 0; 0],
c2 = [ 0; 1; 0], c3 = [−1; 0; 0],
c4 = [ 0;−1; 0], c5 = [ 0; 0; 1],
c6 = [ 0; 0;−1], c7 = [ 1; 1; 0],
c8 = [−1; 1; 0], c9 = [−1;−1; 0],
;
c10 = [ 1;−1; 0], c11 = [ 1; 0; 1],
c12 = [−1; 0; 1], c13 = [−1; 0;−1],
c14 = [ 1; 0;−1], c15 = [ 0; 1; 1],
c16 = [ 0;−1; 1], c17 = [ 0;−1;−1],
c18 = [ 0; 1;−1]
Set t = 0;
Choose relaxation parameter τc ;
Determine time step ∆t;
Set initial values:;
foreach lattice point x do
for i=0 to 18 do
Set initial values fi (xx, 0);
end
end
Time Evolution:;
while t < Tend do
Collision:;
foreach lattice point x do
Compute density: ρ = ∑18 i=0 f i (x
x,t);
Compute velocity: u = ρ1 ∑18 i=0 ci f i (x
x,t);
Compute equilibrium distribution:;
(eq)
f0 = 13 ρ [1 − 32 |u|2 ];
for i=1 to 6 do
(eq) 1
 
fi = 18 ρ 1 + 3(ci · u) + 92 (ci · u)2 − 32 |u|2 ;
end
for i=7 to 18 do
(eq) 1
 
fi = 36 ρ 1 + 3(ci · u) + 92 (ci · u)2 − 32 |u|2 ;
end
for i=0 to 18 do h i
(eq)
fi (xx,t) = fi (xx,t) − ∆t
τc f i (x ,t) − f i
x ;
end
end
Streaming:;
foreach lattice point x do
for i=0 to 18 do
fi (xx + ci ∆t,t + ∆t) = fi (xx,t);
end
end
Update:;
t = t + ∆t;
end
Algorithm 1: Pseudocode for LBM on the D3Q19 lattice.
 before after

Collision

before after

Streaming

Fig. (5): Illustration of streaming and collision.

Nevertheless, there are lattice points which do not have the full number of neighbors, we call these points
boundary lattice points.
For interior points, after the streaming step all populations of the density distribution are determined by pre-
streaming values of adjacent points with the only exception of the value for the rest particle. For boundary
lattice points, there are unknown populations after the streaming step due to missing neighbors. Now,
boundary conditions in lattice Boltzmann models have to define the missing populations, in a way which
leads to the desired macroscopic behavior at the boundary. There are various possibilities. In general, we
find local and nonlocal conditions. A condition is nonlocal if information from neighboring lattice points is
used. This can be combined with conditions which solely compute the unknown populations at the boundary
or with conditions which additionally modify the known populations at boundary lattice points. An example
of a nonlocal boundary condition which recomputes not only unknown populations is presented in [31].
If a local boundary condition only determines the unknown populations, then the action of a boundary
condition can be viewed as an additional streaming step, where the unknown populations are streamed from
imaginary lattice points located outside of the actual numerical domain.
In the following subsections we will give a short description of the most popular boundary conditions when
dealing with the lattice Boltzmann method. These are: periodic, bounceback, Inamuro, and Zou/He condi-
tions. We focus on, how they are applied and their relation to physical conditions as no slip, pressure and
velocity conditions. For this discussion, we restrict the description of boundary conditions to straight bound-
aries which are aligned with the main lattice directions. For boundaries not aligned with lattice directions,
schemes based on extrapolations like in [8] can be used. A consideration of curved boundaries can be found
e.g. in [25].

3.1. Periodic Boundary Conditions

The boundary conditions which are simplest implemented are the so-called periodic boundary conditions.
Here the domain is wrapped around, i.e., corresponding boundary counterpart are identified. The missing
 Bounceback

(a) Initial state at the begin of a time (b) Bounceback condition applied
step. (the streaming step is to follow).

Fig. (6): Illustration of Bounceback Condition for one lattice point (the shaded area represents a boundary
or an obstacle).

neighboring lattice points are replaced with boundary lattice points of the corresponding boundary counter-
part, and vice versa. Thus, populations leaving the domain on one side in the streaming step, re-enter the
domain at the other side.
This kind of boundary condition is usually applied if the domain of interest and the physical process has a
regular symmetry. They are also typically applied in order to get rid of actual physical boundaries, e.g. if
one is interested in the behavior of a multi-component mixture for a given initial state where surface effects
play a negligible role [33].

3.2. Bounceback Boundary Conditions

This condition is often applied for solid walls and obstacles. Its implementation is also not difficult, because
one does not need to consider the orientation of the boundary. For those lattice points where a bounceback
condition is assigned, the collision step, i.e., the relaxation towards equilibrium in the BGK model, is omit-
ted. Instead of the collision step all unknown populations are adopted from those of opposite directions.
The lattice point with bounceback condition normally does not belong to the fluid domain, but the effective
physical boundary lies half way between these bounceback lattice points and the neighboring fluid lattice
points, see Fig. (6). Certainly, the macroscopic quantities (32) cannot be computed in these points, but since
they do not belong to the fluid domain, this is not a problem.
This kind of boundary condition is often implemented if no-slip boundaries are present, however they do not
model a no-slip condition exactly [20]. In fact, in the time average, bounceback yields the no-slip condition.

3.3. Inamuro No-Slip Boundary Conditions

Inamuro et al. [21] presented an approach such that no-slip boundaries can be realized exactly. To this end,
we consider a wall which may move; we denote its velocity by u wall and the outer normal (of the wall
as a boundary) is denoted by n . For a stationary wall, we simply have to set u wall = 0 . Moreover, the
effective boundary of this wall shall pass through lattice points, such that n is aligned with the main lattice
directions. Let us consider a lattice point on this wall, then the macroscopic velocity in the normal direction
un := u wall · n is fixed in that lattice point (by the wall velocity).
In the following, we aim at assigning only the unknown populations. Before we introduce the needed
additional degrees of freedom, we show that the known value un already fixes the macroscopic density at the
wall. Now, let I− = {i ∈ {0, . . . , nv } | fi unknown after streaming step} be the set of indices corresponding
to unknown populations. Then we abbreviate the corresponding sum by

ρ− = ∑ fi .
i∈I−
Moreover, for a given lattice velocity index i let i− be the index of the opposite lattice velocity, i.e., c i− = −cci .
Analogously, we have ρ+ as

ρ+ = ∑ fi with I+ = i ∈ {0, . . . , nv } | i− ∈ I− .
i∈I+

Denoting the remaining populations by ρ0 , the macroscopic quantities (32) are given by

ρ = ρ− + ρ+ + ρ0 , ρun = ρ+ − ρ− . (48)

The equations (48) can be combined to

2ρ+ + ρ0
ρwall := ρ = , (49)
1 + un
which shows that the density is fixed exclusively by the normal velocity un and the known populations.
When modeling the collision as a relaxation towards equilibrium, like in the BGK model, then the discrete
equilibrium distribution can be computed already by (30). By simply replacing all populations, also the
given ones, by their corresponding equilibrium values would result in a less accurate method, cf. [23]. The
idea of the so-called Inamuro boundary conditions is to replace the unknown populations by equilibrium
values corresponding to (30) with ρ → ρ ′ and velocity u → u ′ , such that, we satisfy the velocity condition
from (32) at the wall:

! 1
u wall =
ρwall ∑ c i fi (xx,t). (50)
i

The density ρ ′ is a free parameter and the velocity u ′ is the given wall velocity plus a so-called counterslip
velocity w ′ , i.e., u ′ = u wall + w ′ . For w ′ only the tangential wall components are non-zero, that is, w ′ · n = 0.
Hence the number of unknown components in u ′ is  the space dimension
 minus one. For example in a 2D
u + w′
x,wall
simulation with a bottom wall, it would read u′ = x with the unknown w′ . Thus, the number
x
uy,wall
′ ′
of free velocity component(s) in w plus one unknown density ρ matches the number of constraints in (50).
We note that the wall density is then also automatically fulfilled:

ρwall = ∑ fi (xx,t).
i

This also justifies the need of the free density parameter ρ ′ . For the D2Q9 model, the solution of (50) can
be stated explicitly. For instance at a bottom wall we have
ρwall uy,wall + ρ+ 6 ρwall ux,wall − f1 + f3 + f7 − f8
ρ′ = 6 , w′x = · − ux,wall
1 + uy,wall + (uy,wall )2 ρ′ 1 + 3uy,wall

with populations labeled as indicated in Fig. (2)(b). In the D3Q19 model, this procedure yields a nonlinear
system of equations, which can be solved for instance by the Newton-Raphson method.
Finally we note that Inamuro boundary conditions can also be used to assign velocity boundary conditions
(e.g. for an inlet or outlet) by setting artificial wall velocities.

3.4. Zou/He Boundary Conditions

Zou and He [39] present a possibility to implement both a pressure boundary condition and a velocity
boundary condition in lattice Boltzmann simulations. Next, we address both conditions.
a) Velocity Condition. Again, the aim is to find the unknown populations, such that the macroscopic
velocity matches a prescribed velocity u wall , that is, (50) shall hold. Consider a simulation in d space
dimensions and suppose a velocity is prescribed for a boundary lattice point. Then this will result in d
constraints for the missing populations. For instance using the D2Q9 model, we have three unknowns
(inward directed populations), but only two constraints. This similarly holds for the D3Q19 model. There
we have five missing populations but only three relations to fulfill. Hence additional relations are required.
Now, Zou and He suggest that the non-equilibrium part of the population normal to the boundary is bounced
back. Following this suggestion, we get
(eq) (eq)
fn = fn− − fn− + fn ,

where the index n corresponds to an inward normal direction and n− corresponds to the opposite direction
(outer normal). Note that the discrete equilibrium distribution can be computed when a velocity is prescribed,
since the consideration leading to (49) is also valid here. Thereby, we now have for the D2Q9 model three
equations, which determine the three unknown populations. Unfortunately, the three-dimensional case of
D3Q19 is not covered by this consideration, because there are still more unknowns than equations. One way
to overcome this lack of constraints is the following. First one computes proposals fi⋆ for the still unknown
populations by extending the idea of bouncing back the non-equilibrium part to all unknown populations.
Thus for D3Q19 and a boundary with outer normal vector c 1 = (1, 0, 0)⊤ (i.e., n = 3, see Fig. (4)), we obtain
(eq) (eq) (eq) (eq)
f8⋆ = f8 = f10 − f10 + f8 , f9⋆ = f9 = f7 − f7 + f9 ,
⋆ (eq) (eq) ⋆ (eq) (eq)
f12 = f12 = f14 − f14 + f12 , f13 = f13 = f11 − f11 + f13 .

By this, the velocity normal to the boundary computed by (32) with the proposals matches the prescribed
one. However, generally there is a mismatch for the velocity components tangential to the boundary, i.e., the
term

∆ j := ∑ c i fi − ρwall u wall
i

is in general not equal to the zero vector. Therefore the proposals are adapted to match these velocities as
well. This is done by redistributing ∆ j over the proposed populations. In the example above this reads
1
fi = fi⋆ − c i · ∆ j , (51)
2
with i = 8, 9, 12, 13. Note that in the current example (∆ j)1 = 0 and therefore does not contribute to the
equation above. Also note the coefficient being here 12 , because in each case only two populations contribute
to a direction, hence the momentum excess ∆ j is split in portions of 12 . In D3Q19, the formula (51) can
be used for all walls aligned with the main lattice direction. When using the D3Q15 model, the coefficient
would read 14 instead.
We conclude this paragraph with some remarks: (i) As we have seen above, a prescribed velocity fixes the
density in the boundary lattice points by equation (49). (ii) In principle the statement of part a) can also
be used to model a no-slip boundary by setting the velocity appropriately. (iii) Unlike Inamuro boundary
conditions, the Zou/He boundary conditions are explicit also in three dimensional cases.
b) Pressure Condition. In lattice Boltzmann methods the pressure is proportional to the density, see (43).
Hence, prescribing a pressure in boundary lattice points is equivalent to prescribing a density ρwall . Writing
equation (49) as
2ρ+ + ρ0
un = −1
ρwall
shows that a prescribed pressure determines the normal velocity. Now, the easiest way to continue is to
enforce additionally the magnitude of the tangential velocity. This leads to the same situation as described
for velocity boundary conditions and the missing populations can be found as stated above. A typical choice
for the tangential velocity is setting it equal to zero, which gives a no-slip condition.
This completes the standard boundary conditions. We consider now a collection of enhancements of the
basic LBM concept.
4. STANDARD ENHANCEMENTS OF THE BASIC CONCEPT

So far, we have introduced the basic version of LBM. Many applications e.g. from computational physics,
mechanics and life sciences need enhancements and further couplings to other models. We discuss here
multiphase and multicomponent flows and the incorporation of additional forces (as gravitation). Further-
more, we address couplings to other phenomena as temperature and we roughly present an application and
a coupling for chip design. There are other important couplings: for example applications in fluid struc-
ture interaction (from civil engineering), coupling to chemical reaction kinetics (for e.g. the behavior of the
atmosphere), which we just name here.

4.1. Multiphase/Multicomponent Flows

Up to now we have considered fluids consisting of a single component only. Also our consideration was
limited to a single phase of the fluid. However, in practical applications one often has a mixture of several
components. The most common way to incorporate these extensions to the lattice Boltzmann method is the
multiphase/multicomponent model of Shan and Chen [29]. In this model nonlocal interactions among the
different species are incorporated.
Next, we shortly summarize how a fluid consisting of C components is modeled in a lattice Boltzmann
method. Each component evolves, correspondingly to (12), as follows
∆t h m m(eq)
i
fim (xx + c i ∆t,t + ∆t) − fim (xx,t) = − fi (xx,t) − fi (xx,t) , m = 1, . . . ,C
τm
where fim is the single-particle distribution for the mth component. Note that there are C independent re-
m(eq)
laxation parameter τm , which control the viscosity of each component. As usual the fi , m = 1, . . . ,C,
are identical to the equilibrium distribution of a single component simulation, i.e. (30), with accordingly
computed macroscopic quantities ρm and u ′m [30]:

ρm (xx,t) = ∑ fim (xx,t),

i
! !−1
C C
ρ j (xx,t)uu j (xx,t) ρ j (xx,t) τm
u′m (xx,t) = ∑ ∑ τj + Fm (xx,t)
j=1 τj j=1 ρm (xx,t)

where Fm (xx,t) models the momentum change of the mth component due to interactions. The velocities
u m (xx,t) for m = 1, . . . ,C are given by

ρm (xx,t)uum (xx,t) = ∑ c i fim (xx,t).

i

The interaction force Fm (xx,t) can be expressed by

C
Fm (xx,t) = −ρm (xx,t) ∑ Gm, j ∑ ρ j (xx + c i ,t)cci ,
j=1 i

with a symmetric matrix G ∈ RC×C having diagonal elements equal to zero. The entry Gm, j controls the
strength of the interaction between the components m and j.

4.2. Additional Forces

In absence of external forces driving a fluid, the motion of a fluid decreases due to internal friction. The
amount of friction is given by the fluid’s viscosity. Typically in real experiments there is at least one force
acting on the fluid coming from gravity. There are many other circumstances which make it desirable to be
able to incorporate internal or external forces in numerical simulations.
On the macroscopic scale, where the Navier-Stokes equations describe the motion of the fluid, forces can be
incorporated by simply adding an acceleration term corresponding to the desired force. Thus, in the lattice
Boltzmann sense it is the task to find a modification which leads to
∇ · u = 0,
∂u (52)
∇ p̂ + ν∆uu + k ,
+ (uu · ∇ ) u = −∇
∂t
in the multiscale expansion where k is the acceleration due to the force density K = ρkk (and p̂ = p/ρ).
One possibility to realize it is given by a modification proposed by Guo et al. [15]. They state that the fluid
velocity has to be redefined and the contribution of the force to the momentum flux must be canceled. In their
article a modified lattice Boltzmann method is proposed fulfilling these requirements. We briefly explain the
necessary modifications. Let K (xx,t) be a given force density, then a corresponding discrete quantity reads
   
∆t (cci − u ⋆ (xx,t)) (cci · u ⋆ (xx,t))
Ki (xx,t) = 1 − wi 3 +9 c i · K (xx,t), (53)
2τc c2 c4

where the adapted velocity u ⋆ is now computed by

!
nv
⋆ 1 K (xx,t)
u (xx,t) =
ρ(xx,t) ∑ c i fi (xx,t) + 2 (54)
i=0

instead of using (32). The discrete force density (53) is used in a modified collision step. The complete
modified lattice Boltzmann method satisfying (52) reads

∆t h (eq)
i
fi (xx + c i ∆t,t + ∆t) − fi (xx,t) = − fi (xx,t) − fi (xx,t) + ∆t Ki (xx,t).
τc

4.3. Thermal Coupling

The goal of this subsection is to develop a numerical method to simulate the temperature evolution. There-
fore we propose a method which extends the lattice Boltzmann method (with BGK approximation) such
that the temperature is advected by the fluid’s flow and diffuses in it, meaning the temperature obeys an
advection-diffusion equation

∂
(ρ(xx,t)E(xx,t)) + ∇ ·(ρ(xx,t)uu(xx,t)E(xx,t)) = χ∆ (ρ(xx,t)E(xx,t)) , (55)
∂t
with thermal diffusivity χ. The temperature T is related to the internal energy by

3 kB T (xx,t)
E(xx,t) =
2 m
in R3 . Like viscosity in the lattice Boltzmann method, a parameter which can be chosen to tune the thermal
diffusivity is available in the numerical scheme.
Basically, the advection-diffusion equation (55) is coupled to the Navier-Stokes equations only in one way.
The fluid is not affected due to the temperature, this means it behaves like an isothermal fluid, which es-
pecially means that the viscosity is temperature independent. Consequently, this limits the application of
the proposed method to those problems where the influence of the temperature on the viscosity is small.
However, an external force dependent on the temperature or more precisely temperature variations can be
incorporated to the basic lattice Boltzmann method. By this we indirectly achieve a two way coupling. Fur-
thermore, the compression work done by the pressure and the viscous heat dissipation are neglected in the
method we state here. When simulating the motion of real incompressible fluids in applications these terms
are (very often) negligible. A method which incorporates also the compression work done by the pressure
and the viscous heat dissipation is presented in [18].
We define a new distribution, calling it the energy distribution, by

|vv − u (xx,t)|2
g(xx, v ,t) = f (xx, v ,t).
2
The derivatives of g can be computed in terms of the single particle distribution f . The latter is a solution of
the Boltzmann equation, therefore we get by straightforward calculations

∂ g(xx, v ,t) |vv − u (xx,t)|2

+ v · ∇ g(xx, v ,t) = Q( f ) − f (xx, v ,t)q(xx, v ,t), (56)
∂t 2
with
 
∂ u (xx,t)
q(xx, v,t) = (vv − u(xx,t))· v ∇
+ (v · )u (x ,t) .
u x
∂t

The term f (xx, v ,t)q(xx, v ,t) only contributes to the compression work and the heat dissipation. As explained
above these ingredients are neglected here, therefore we can neglect this term in (56) as well. Furthermore
we approximate

|vv − u (xx,t)|2 1 h i
Q( f ) ≈ − g(xx, v ,t) − g(eq) (xx, v ,t)
2 τd
in a BGK consistent way, with

|vv − u (xx,t)|2 (eq)

g(eq) (xx, v ,t) := f (xx, v ,t). (57)
2
Thus, we end up with an evolution equation for the energy distribution

∂ g(xx, v ,t) 1 h i
+ v · ∇ g(xx, v ,t) = − g(xx, v ,t) − g(eq) (xx, v ,t) . (58)
∂t τd
Due to (5) it holds
Z
g(xx, v,t) dvv = ρ(xx,t)E(xx,t),
R3

and the first two moments of the energy equilibrium distribution can be computed explicitly
Z
g(eq) (xx, v ,t) dvv = ρ(xx,t)E(xx,t),
R3
Z (59)
v g(eq) (xx, v ,t) dvv = ρ(xx,t)uu(xx,t)E(xx,t).
R3

When performing a Chapman-Enskøg expansion, analogue to the procedure of Section 2.3., these results are
used. The Chapman-Enskøg expansion verifies that the outcome matches the advection-diffusion equation
(55).
A comparison of the structure of equations (58) and (11) shows that they differ only in a renaming of the
distribution f → g. Hence a scheme like (12) can be used, it follows
∆t h (eq)
i
gα (xx + c α ∆t,t + ∆t) − gα (xx,t) = − gα (xx,t) − gα (xx,t) . (60)
τd
The parameter τd in this equation can be used to tune the thermal diffusivity [28]. The equation (60) is not
well defined, since the discrete energy equilibrium distribution is not determined. However, using the defi-
nitions of the energy equilibrium distribution (57) and the Maxwellian distribution (10), a Taylor expansion
like in the derivation of (27) yields
 3/2   2  2  
(eq) 1 |vv|2 |vv| |vv| 2 (vv · u ) (vv · u )2 |uu|2
g = ρE exp − + − + −
2πRT 2RT 3RT 3RT 3 RT 2(RT )2 2RT
 3/2    2   2  
1 |vv|2 |vv| 7 (vv · u )2 |vv| 5 |uu|2
+ ρE exp − − − − + O(|uu|3 ).
2πRT 2RT 3RT 3 2(RT )2 3RT 3 2RT

The second line has no contribution to the integrals (59), hence this addend can be neglected. Thus the
quadrature formula (29) can be used for the discrete computation of the integrals (59):
Z 18
(eq)
ψ(vv)g(eq) (xx, v ,t) dvv ≈ ∑ ψ(ccα )gα ,
α=0
R3

with the discrete energy equilibrium distribution

   
(eq) |ccα |2 |ccα |2 2 3(ccα · u ) 9(ccα · u )2 3|uu|2
gα = ρEwα + − + − . (61)
c2 c2 3 c2 2c4 2c2

4.4. Circuit Coupling

As one specific example, we want to model the coupling of integrated electric circuits using a fluid (air)
current produced for instance by some fan. Let Ω ⊂ R3 denote the fluid region, where air can flow. Due to
(11) and (58) the fluid and temperature inside Ω is governed by

∂ f (xx, v ,t) 1h i
+ v · ∇ f (xx, v ,t) = − f (xx, v ,t) − f (eq) (xx, v ,t) ,
∂t τc
∂ g(xx, v ,t) 1 h i
+ v · ∇ g(xx, v ,t) = − g(xx, v ,t) − g(eq) (xx, v ,t) .
∂t τd
The relation to the temperature is as follows (cf. (5)):
Z
2
T (xx,t) = g(xx, v ,t) dvv.
3ρ(x ,t)R
x
R3

The electric network model shall be described by the modified nodal analysis (MNA), see for instance [12]:

d
(Aqq(ee, j )) = h (ee, j ,t)
dt
with unknowns e and j being the node voltages and currents through voltage defining elements (voltage
sources, inductors), respectively. The term q models charges and fluxes, h is the so-called static part includ-
ing independent sources.
We notice that the term h includes thermally relevant resistors with incidence matrix AT . We separate h for
the following using the corresponding conduction matrix GT :

h (ee, j ,t) = h̃h(ee, j ,t) + AT GT A⊤

Te

We remark that within MNA no geometry information is present, which we need here to supply such that
the physical location of thermally relevant elements is defined.
The coupling of the models for fluid/temperature and the circuit consists in two ways. In the circuit model
devices depend on an environment temperature, and on the other hand, the dissipated power of the electric
circuit is an energy input for the temperature model. Let the ith thermal active element show a common
interface Γi with the fluid, the environment temperature Tenv,i of this ith thermal element can be obtained
directly out of the fluid/temperature model, e.g. by averaging:
Z Z Z
1 2
Tenv,i = T (xx,t) dx = g(xx, v ,t) d v d x.
kΓi k 3ρRkΓi k
Γi Γi R 3

To set up a coupling in the other direction, we need to compute the dissipated power of the electric circuit.
Using all thermally relevant resistors for heat dissipation we obtain the dissipated power for the ith element
as [2]
1 2
Pi = ei · ji = e
Ri i
(with resistance Ri ). This is a production term for the energy balance of the ith thermal element. Further-
more, we model the energy transfer to the fluid via Newton cooling on the active boundary Γi , such that we
obtain
d
cv Ti = Pi − γkΓi k(Ti − Tenv,i )
dt
with the lumped heat capacity cv and the parameter γ denoting the transfer coefficient. Then the boundary
condition can be formulated for T (for the heat flux):
 
∂ T (xx,t)
λ n· = −γ(Ti − Tenv,i ) on Γi , (62)
∂x

where n is the unit outer normal vector of the boundary Γi and λ denotes the heat conductivity (plus further
boundary conditions). Assuming the boundary conditions for the fluid, i.e., for f , and the conditions for
temperature on remaining boundaries are determined, the continuous problem is fully defined. Next we
have to consider how the continuous description can be transferred to the discrete model. We again restrict
this consideration to boundaries aligned with the main lattice direction. We focus only on the transfer of
the boundary condition (62) to the lattice Boltzmann method (60). Due to the assumption of alignment of
the boundaries with the main lattice direction, except for sign the normal vector is a unit vector, hence the
left hand side of (62) is simply a partial derivative. Especially, there is a lattice velocity with cβ = cnn. The
non-central finite difference quotient of second order can be used to compute a temperature on the boundary,
which embodies the von Neumann condition (62):
 
1 1
T (xx,t) = 4T (xx + c β ∆t,t) − T (xx + 2ccβ ∆t,t) − 2c∆tγ (Ti − Tenv,i ) . (63)
3 λ

Therefore we have transferred the von Neumann boundary condition to a Dirichlet boundary condition.
With the temperature (63) an equilibrium distribution according to (61) can be computed. Similar to the
Zou/He velocity boundary conditions, the missing populations of g can be obtained by setting them to the
equilibrium values whereas the non-equilibrium part is bounced back.

5. CONCLUSIONS

We have recapitulated the derivation of LBM from two perspectives: one was the historical development
from lattice gas automata and the other was a more recent derivation using certain discretizations and Gauß-
Hermite quadrature. The second approach was here transferred to the D3Q19 lattice. A pseudocode for LBM
using D3Q19 concluded this section. Then the basic and standard boundary conditions were explained with
 remarks for implementation and applicability. After this part, we have discussed a couple of enhancements
for LBM and further trends. We have also added a rough mathematical modeling for a coupled system from
chip design, where the fluid and temperature could be simulated via LBM. This thermal electric coupling
becomes more and more important due to the increasing power challenges in chip design.

REFERENCES

[1] D. J. Acheson. Elementary fluid dynamics. Clarendon Press, Oxford, 1998.

[2] A. Bartel. Partial Differential-Algebraic Models in Chip Design-Thermal and Semiconductor Problems. VDI Verlag, Düsseldorf, Germany, 2004.
[3] P. L. Bhatnagar, E. P. Gross, and M. Krook. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral
One-Component Systems. Physical Review, 94(3):511–525, May 1954.
[4] C. Cercignani. The Boltzmann equation and its applications. Springer Verlag, New York, USA, 1988.
[5] S. Chapman and T. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, London, UK, 1970.
[6] H. Chen, S. Chen, and W. H. Matthaeus. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Physical Review A,
45(8):R5339–R5342, Apr 1992.
[7] S. Chen, K. Diemer, G. D. Doolen, K. Eggert, C. Fu, S. Gutman, and B. J. Travis. Lattice gas automata for flow through porous media. Physica
D: Nonlinear Phenomena, 47(1-2):72–84, 1991.
[8] S. Chen, D. Martinez, and R. Mei. On boundary conditions in lattice Boltzmann methods. Physics of Fluids, 8(9):2527–2536, 1996.
[9] J. P. Dahlburg, D. Montgomery, and G. D. Doolen. Noise and compressibility in lattice-gas fluids. Physical Review A, 36(5):2471–2474, Sep 1987.
[10] D. d’Humieres, I. Ginzburg, M. Krafczyk, P. Lallemand, and L.-S. Luo. Multiple-relaxation-time lattice Boltzmann models in three dimensions.
Philosophical Transactions of the Royal Society A, 360:437–451, 2002.
[11] D. d’HumiÃĺres, P. Lallemand, and G. Searby. Numerical Experiments on Lattice Gases: Mixtures and Galilean Invariance. Complex Systems,
1:633–647, 1987.
[12] U. Feldmann and M. Günther. CAD-based electric-circuit modeling in industry I: mathematical structure and index of network equations. Surveys
on Mathematics for Industry, 8(2):97–129, 1999.
[13] U. Frisch, D. d’Humieres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.-P. Rivet. Lattice gas hydrodynamics in two and three dimensions.
Complex Systems, 1:649–707, 1987.
[14] U. Frisch, B. Hasslacher, and Y. Pomeau. Lattice-Gas Automata for the Navier-Stokes Equation. Physical Review Letters, 56(14):1505–1508, Apr
1986.
[15] Z. Guo, C. Zheng, and B. Shi. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Physical Review E, 65(4):046308, Apr
2002.
[16] D. Hänel. Molekulare Gasdynamik. Springer Verlag, 2004.
[17] J. Hardy, Y. Pomeau, and O. de Pazzis. Time evolution of a two-dimensional model system. I. Invariant states and time correlation functions.
Journal of Mathematical Physics, 14(12):1746–1759, 1973.
[18] X. He, S. Chen, and G. D. Doolen. A novel thermal model for the lattice Boltzmann method in incompressible limit. Journal of Computational
Physics, 146(1):282–300, 1998.
[19] X. He and L.-S. Luo. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Physical Review
E, 56(6):6811–6817, Dec 1997.
[20] X. He, Q. Zou, L.-S. Luo, and M. Dembo. Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice Boltzmann
BGK model. Journal of Statistical Physics, 87(1-2):115–136, 1997.
[21] T. Inamuro, M. Yoshino, and F. Ogino. A non-slip boundary condition for lattice Boltzmann simulations. Physics of Fluids, 7:2928–2930, 1995.
[22] M. Kutrib, R. Vollmar, and T. Worsch. Introduction to the special issue on cellular automata. Parallel Computing, 23(11):1567–1576, 1997.
[23] J. Latt, B. Chopard, O. Malaspinas, M. Deville, and A. Michler. Straight velocity boundaries in the lattice Boltzmann method. Physical Review E,
77(5):056703, May 2008.
[24] G. R. McNamara and G. Zanetti. Use of the Boltzmann Equation to Simulate Lattice-Gas Automata. Physical Review Letters, 61(20):2332–2335,
Nov 1988.
[25] R. Mei, W. Shyy, and D. Yu. Lattice Boltzmann Method for 3-D Flows with Curved Boundary. Journal of Computational Physics, 161:680–699,
2000.
[26] Y. Qian, D. d’Humieres, and P. Lallemand. Lattice BGK Models for Navier-Stokes Equation. Europhysics Letters, 17(6):479–484, 1992.
[27] L. Saint-Raymond. Hydrodynamic Limits of the Boltzmann Equation. Springer Verlag, Berlin, Germany, 2009.
[28] Y. Peng, C. Shu, and Y. T. Chew. A 3D incompressible thermal lattice Boltzmann model and its application to simulate natural convection in a
cubic cavity. Journal of Computational Physics, 193(1):260–274, 2004.
[29] X. Shan and H. Chen. Lattice Boltzmann model for simulating flows with multiple phases and components. Physical Review E, 47(3):1815–1819,
Mar 1993.
[30] X. Shan and G. Doolen. Multicomponent lattice-Boltzmann model with interparticle interaction. Journal of Statistical Physics, 81:379–393, 1995.
10.1007/BF02179985.
[31] P. A. Skordos. Initial and boundary conditions for the lattice Boltzmann method. Physical Review E, 48(6):4823–4842, Dec 1993.
[32] J. D. Sterling and S. Chen. Stability Analysis of Lattice Boltzmann Methods. Journal of Computational Physics, 123(1):196 – 206, 1996.
[33] S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford University Press, Oxford, UK, 2001.
[34] R. Temam. Navier-Stokes equations: theory and numerical analysis. AMS Chelsea Publishing, Providence, RI, USA, 2001.
[35] J. von Neumann. Theory of Self-Reproducing Automata. University of Illinois Press, Champaign, IL, USA, 1966.
[36] D. Wolf-Gladrow. Lattice-Gas Cellular Automata and Lattice Boltzmann Models. Springer Verlag, Berlin, Germany, 2000.
[37] S. Wolfram. Cellular automaton fluids 1: Basic theory. Journal of Statistical Physics, 45:471–526, 1986. 10.1007/BF01021083.
[38] G. Zanetti. Hydrodynamics of lattice-gas automata. Physical Review A, 40(3):1539–1548, Aug 1989.
[39] Q. Zou and X. He. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Physics of Fluids, 9(6):1591–1598, 1997.

View publication stats