where σ(Ω) is the differential collision cross section for the two particle collision
which transform the velocities: {~v1, ~v2} (before collision) into {~v′1, ~v′2} (after colli-
sion) [1, 3].
The Boltzmann equation:
(
∂
∂t
+ ~v ∂
∂~x
+ ~Fext
∂
∂~v
)
f(~x,~v, t) = Q(f1, f2) (1.9)
proved its utility in various areas of modern physics, especially those related to
transport phenomena: fluid mechanics (hydrodynamics and aerodynamics), neu-
tron transport, plasma physics, electron transport through semiconductors, etc.
If we wish to get the hydrodynamics equations from the Boltzmann equation,
we should first define the local thermodynamical equilibrium. The equilibrium
state is described by the equilibrium distribution function f eq, which satisfies the
Boltzmann equation. At local equilibrium, we have [3, 5]:
Q(f eq, f eq) = 0 (1.10)
This condition leads to:
f ′1f
′
2 = f1f2 (1.11)
which does not mean that the particles stand still; in the equilibrium state the
effects induced by direct and inverse collisions are balanced.
In the logarithmic description, the condition (1.11) becomes:
lnf ′1 + lnf ′2 = lnf1 + lnf2 (1.12)
The quantity lnf is a collision invariant. An immediate consequence is that, at
thermodynamic equilibrium, lnf depends on the collision invariants [3, 5]:
• mass density
ρ = m
∫
f d3v (1.13)
• momentum density
ρuα = m
∫
fvα d
3v (1.14)
• energy density
ρe =
∫
1
2
mv2αf d3v (1.15)
and can be written as:
lnf = A+Bαvα +
1
2
Cv2 (1.16)
where A,Bα, C (α = 1, . . . , 3) are Lagrangian multipliers.
3
The fluid temperature T is defined from the kinetic energy density DnkBT/2,
where D ∈ {1, 2, 3} is the space dimension, kB is the Boltzmann constant and u
is the local speed:
D
2
nkBT =
1
2
m
∫
f(v − u)2 dDv (1.17)
Solving the equation (1.16) with the Lagrange multiplier method leads to the
Maxwell-Boltzmann distribution function [1, 2, 3, 5]:
f eq = ρ
[ m
2pikBT
]D/2
exp
[
− m
2kBT
(v − u)2
]
(1.18)
The thermal speed of the fluid particles is given by [2, 3]:
vT =
√
kBT
m (1.19)
The equilibrium distribution function (1.18) satisfies the following relations
[1, 2, 5, 7, 8, 9, 10]:
m
∫
f eq dDv = ρ (1.20)
m
∫
f eqvα d
Dv = ρuα (1.21)
m
∫
f eqvαvβ d
Dv = nkBTδαβ + ρuαuβ (1.22)
m
∫
f eqvαvβvγd
Dv = nkBT (uαδβγ + uβδγα + uγδαβ)
+ ρuαuβuγ (1.23)
m
∫
f eqvαvβvγvδδγδ d
Dv = nkBT
[
(D + 2) kBTm + u
2
]
δαβ
+ ρuαuβ
[
(D + 4) kBTm + u
2
]
(1.24)
Using the equations (1.22) - (1.23) and the relation δαβδαβ = D (here we apply
the summation convention over repeating indices) we can easily get:
1
2
mδαβ
∫
f eqvαvβ d
Dv =
[ D
2
kB T
m +
1
2
u2
]
ρ (1.25)
1
2
mδβγ
∫
f eqvαvβvγ d
Dv =
[ D + 2
2
kBT
m +
1
2
u2
]
ρuα (1.26)
A major problem that appears when solving the integro-diferential Boltzmann
equation (1.9) is the complicated form of the collision term. A possible simplifica-
tion of this term was proposed almost simultaneously in 1954 by Bhatnagar, Gross,
4
Krook [11] and independently, by Welander [12]. This simplification, widely known
as BGK aproximation [3, 4, 5], is based on the operator L(f), which replaces the
Boltzmann collision operator Q(f, f) and satisfies the following conditions [4, 5]:
1. L(f) conserves the collision invariants ψ ∈ {ρ, ρuα, 12ρu
2}, so we have:
∫
L(f)ψ(v)d3v = 0
2. The distribution function f(x,v, t) relaxes to the Maxwell - Boltzmann equi-
librium distribution function (1.18):
L(f) = −f(x,v, t)− f
eq(x,v, t)
τ
(1.27)
The relaxation time τ in (1.27) is correlated with the mean time between two
consecutive collisions of the particle [3, 5, 11, 12].
If we replace the collision operator in the Boltzmann equation (1.9) with the
linearized collision term (1.27), we get the following form of the evolution equation
of the distribution function f :
∂tf + vβ∂βf + aβ∂vβf = −
1
τ
[ f − f eq ] (1.28)
(we used the notation ∂t = ∂/∂t, ∂β = ∂/∂xβ , ∂vβ = ∂/∂vβ , as well as the
summation rule over the Greek indices representing the space coordinates from 1
to D, where D is the space dimension).
The particle acceleration can be written as a sum of two terms:
aβ =
Fβ
m = a
extern
β + along rangeβ (1.29)
The first term in the right hand side of equation (1.29) is due to the action of
external forces (including gravity). The second one, due to long range interparticle
forces, is acting when the fluid is not an ideal gas. We may assume that both aexternβ
and along rangeβ , as well as the total acceleration aβ, are independent of the particle
velocities v.
Since the distribution function f(x,v, t) goes to zero when the velocity v = |v|
goes to infinity [1], the following relation is valid for an arbitrary function Φ =
Φ(x,v, t): ∫
Φaβ∂vβfdDv = − aβ
∫
f∂vβΦdDv (1.30)
This equation will be used in the next chapter, where we present the procedure
to recover the conservation equations for mass, momentum and energy from the
Boltzmann equation (1.28). For this purpose we will use the Grad’s moments
method [13, 14], where Φ ∈ { 1, mvα, mvαvβδαβ/2 }.
5
2Lattice Boltzmann models for
isothermal fluids
2.1 Chapman - Enskog procedure and the con-
servation equations
The distribution function f can be expanded to a Chapman - Enskog series with
respect to the Knudsen number ε [3, 5]
f ' f (0) + εf (1) + ε2f (2) + . . . (2.1)
We can use two time scales t1, t2, and a length scale x1 [3]. The time and space
derivatives in the evolution equation (1.28) are replaced with:
∂t = ε∂t1 + ε2∂t2 , ∂β = ε∂β1 (β = 1, . . . , D) (2.2)
Using (2.1) and (2.2) in the Boltzmann equation (1.28) and sepparating the
corresponding powers of ε, we get the evolution equations to zero, first and second
order:
1
τ
[
f (0) − f eq
]
= 0 (2.3)
∂t1f
(0) + vβ∂β1f (0) = −
1
τ
f (1) − aβ∂vβf (0) (2.4)
∂t2f
(0) + ∂t1f (1) + vβ∂β1f (1) = −
1
τ
f (2) − aβ∂vβf (1) (2.5)
The first and second order conservation equations for mass, momentum and
energy are obtained from (2.4) and (2.5) after muliplication by
Φ ∈ {m, mvα, mvαvβδαβ/2 } (2.6)
6
followed by integration over the velocity space.
From Eqs. (1.13) - (1.15), (1.20) - (1.22) and (2.3) we get:
f (0) = f eq (2.7)
as well as some useful relations concerning the momenta of the functions f (l),
l > 0:
∫
f (l) dDv = 0 ∀l > 0 (2.8)
∫
f (l)vα d
Dv = 0 ∀l > 0 (2.9)
∫
f (l)vαvβ d
Dv = 0 ∀l > 0 (2.10)
These relations will be used further to get the conservation equations.
2.2 Conservation equations
2.2.1 First order equations
By multiplying the Eq. (2.4) with Φ = 1, Φ = m and Φ = mvα, respectively,
followed by integration over the whole velocity space we get the first order con-
servation equations for particle density n, mass density ρ and momentum density
ρuα:
∂t1n + ∂β1 (nuβ) = 0 (2.11)
∂t1ρ + ∂β1 (ρuβ) = 0 (2.12)
∂t1(ρuα) + ∂β1 Π(0)αβ = ρaα (2.13)
where:
Π(0)αβ = m
∫
f (0)vαvβd
Dv = nkBTδαβ + ρuαuβ (2.14)
The first order energy equation is obtained in a similar way, using Φ =
mvαvβδαβ/2, as well as equations (1.17) and (1.26):
∂t1
[ D
2
nkBT +
1
2
ρu2
]
+ ∂β1
[( D + 2
2
nkBT +
1
2
ρu2
)
uβ
]
= ρaβuβ
(2.15)
7
The first order equations above lead to:
∂t1(ρuαuβ) = −
[
uα∂γ1Π
(0)
βγ + uβ∂γ1Π(0)αγ
]
+ uαuβ∂γ1(ρuγ)
+ ρ(aαuβ + aβuα) (2.16)
∂t1
(
1
2
ρu2
)
= −uβ [ ∂β1(nkBT ) + ∂γ1(ρuβuγ) − ρaβ ]
+ 1
2
u2∂γ1(ρuγ) (2.17)
∂t1(kBT ) = −uβ∂β1(kBT ) −
2
D kBT∂β1uβ (2.18)
2.2.2 Mass conservation equation
By multiplying Eq. (2.5) with the term Φ = m and taking into account Eqs. (2.8)
and (2.9), we get:
∂t2ρ = 0 (2.19)
The final form of the mass equation (up to second order with respect to Knudsen
number ε) is obtained from (2.12) and (2.19) taking into account Eq. (2.2):
∂tρ + ∂β (ρuβ) = 0 (2.20)
2.2.3 Momentum conservation equation
From Eq. (2.5) multiplied by Φ = mvα, we get:
∂t2(ρuα) + ∂β1 Π(1)αβ = 0 (2.21)
where
Π(1)αβ = m
∫
f (1)vαvβd
Dv (2.22)
We have replaced f (1) in the definition of Π(1)αβ . Then, using the properties
(1.30), we get:
Π(1)αβ = − τ
[
∂t1Π
(0)
αβ + ∂γ1Π
(0)
αβγ
]
+ τρ [ aαuβ + aβuα ] (2.23)
In the above relation, we have:
Π(0)αβγ = m
∫
f (0)vαvβvγd
Dv = uαΠ(0)βγ + uβΠ(0)γα + uγΠ
(0)
αβ − 2ρuαuβuγ (2.24)
We replace Π(1)αβ in Eq. (2.21) and get, using Eqs. (2.16), (2.18) and (2.24):
∂t2(ρuα) + ∂β1
{
τnkBT
[
∂αuβ + ∂βuα −
2
D δαβ(∂γuγ)
] }
= 0 (2.25)
8
When adding the above equation to the first order equation (2.13), we get the
final form [8, 9, 10] of the momentum conservation equation (up to second order
with respect to the parameter ε)
∂t(ρuα) + ∂β (ρuαuβ ) =
− ∂αp + ∂β
{
η
[
(∂αuβ + ∂βuα) −
2
D δαβ(∂γuγ)
] }
+ ρaα (2.26)
where p is the ideal gas pressure, defined by the equation of state:
p = nkBT (2.27)
and
η = τnkBT (2.28)
is the dynamic viscosity of the fluid [8, 15].
2.2.4 Energy conservation equation
After multiplication of Eq. (2.5) with Φ = mvαvβδαβ/2 and use of the tempera-
ture definition (1.17), we get:
∂t2
[ D
2
nkBT +
1
2
ρu2
]
+ ∂γ1
[
1
2
δαβΠ(1)αβγ
]
= 0 (2.29)
where
Π(1)αβγ = m
∫
f (1)vαvβvγd
Dv = (2.30)
= − τ
[
∂t1Π
(0)
αβγ + ∂δ1Π
(0)
αβγδ
]
+ τ aδ
[
δγδΠ(0)αβ + δαδΠ
(0)
βγ + δβδΠ(0)αγ
]
We now use the following form of the first order equations (2.12), (2.13) and (2.18):
∂t1ρ = − ρ∂δ1uδ − uδ∂δ1ρ (2.31)
∂t1uα = aα − ∂α1
(
kBT
m
)
−
(
kBT
mρ
)
∂α1ρ − uδ∂δauα (2.32)
∂t1
(
kBT
m
)
= −uδ∂δ1
(
kBT
m
)
− D
2
τ
kBT
m ∂δ1uδ (2.33)
as well as Eqs. (1.24) - (1.26) to get
1
2
δαβΠ(1)αβγ = −
[ D + 2
2
τn
kBT
m
]
∂γ1(kBT )
− τnkBTuδ
[
∂γ1uδ + ∂δ1uγ −
2
D δγδ(∂χ1uχ)
]
(2.34)
9
If we add Eq. (2.29) to the first order equation (2.15), we get the energy
equation up to second order with respect to ε [8, 9, 10]
∂t
[ D
2
nkBT +
1
2
ρu2
]
+ ∂β
[( D
2
nkBT +
1
2
ρu2
)
uβ
]
−
− ∂γ
[( D + 2
2
τn
kBT
m
)
∂γ(kBT )
]
= ρaγuγ (2.35)
2.3 Phase space discretization
The starting point of the Lattice Boltzmann models is the Boltzmann transport
equation, where the collision term has been linearized using the BGK approxima-
tion [1, 3, 4, 5, 11]
(
∂
∂t
+ v · ∇ + Fm · ∇v
)
f(x,v, t) = − 1
τ
[ f(x,v, t) − f eq(x,v, t) ] (2.36)
As mentioned in the Introduction, the equilibrium distribution function f eq(x,v, t)
for an ideal gas is the well-known Maxwell - Boltzmann distribution (1.18). If the
system is not too far from the equilibrium state, we may assume
∇vf(x,v, t) ' ∇vf eq(x,v, t) = − m
kBT
[v − u(x, t) ] f eq(x,v, t) (2.37)
After introducing the expression (2.37) in (2.36), we get the following form of the
Boltzmann equation:
∂tf(x,v, t) + v · ∇f(x,v, t) = 1
kBT
F · [v − u(x, t) ] f eq(x,v, t)
− 1
τ
[ f(x,v, t) − f eq(x,v, t) ] (2.38)
The independent variables x and v generate the phase space and belong to
the continuum spectrum. The basic ideea of Lattice Boltzmann models is the
discretization of the phase space [5, 16, 17, 18]. According to the discretization
procedure, the distribution functions are now defined only in the nodes of a lattice
L, which belongs to the D-dimensional space (D ∈ { 1, 2, 3 }). The velocity space
is reduced to a finite set of vectors { ei }, i = 0, 1, . . . N . After discretization of
the phase space, the distribution function f(x,v, t) is replaced by a finite set of N
distribution functions fi(x, t) ≡ f(x, ei, t). Thus, the Boltzmann equation (2.38)
is replaced by a system of N evolution equations:
∂tfi(x, t) + ei · ∇fi(x, t) =
1
kBT
F · [ ei − u(x, t) ] f eqi (x, t)
10
− 1
τ
[ fi(x, t) − f eqi (x, t) ] (2.39)
(i = 0, 1, . . . N )
As a result of the discretization of the phase space, the integrals over the
velocity space are replaced by sums computed over the discrete set { ei } [3, 5, 19,
20]. The expression of the mass density (1.13) becomes:
ρ ≡ ρ(x, t) =
N∑
i=0
fi(x, t) (2.40)
and the local speed (1.14) is replaced with
u ≡ u(x, t) = 1
ρ(x, t)
N∑
i=0
ei fi(x, t) (2.41)
In the sequel we will consider only the isothermal Lattice Boltzmann models,
where temperature is constant. In these models the equilibrium distribution func-
tions are expressed as series expansion with respect to the local fluid velocity u,
up to the second order:
f eqi = A+Beiαuα + Cu2 +Duαuβeiαeiβ
i = 0, 1, . . . , N (2.42)
where eiα ≡ (ei)α and uα (i = 1, 2, . . . N , α = 1, 2, 3) are the Cartesian
components of the vectors ei and u. The coefficients A, B, . . . of this series are
determined from the condition that the zero and first order momentum of the
Boltzmann equation (2.39) are equivalent with the continuity equation and Navier-
Stokes equation, respectively [16, 17, 18].
The analytical expression of the equilibrium distribution functions f eqi was
obtained using the Gauss - Hermite quadrature to compute the hydrodynamic
moments up to fourth order [21, 22]:
f eqi = f eqi (x, t) = wi ρ
[
1 + ei · u
χ c2
+ (ei · u)
2
2χ2 c4
− u · u
2χ c2
]
(i = 0, 1, . . . N ) (2.43)
The parameters wi and χ are characteristics of the isothermal Lattice Boltzmann
model. The speed c in Eq. (2.43) is correlated to the thermal speed (1.19) of the
fluid particles [23]:
c =
√
kB T
χm (2.44)
11
To get the above expression of the speed c, we impose the condition that the
discretization procedure conserves the hydrodynamic moments of the Boltzmann
equation [14, 21, 22] up to second order.
In the isothermal Lattice Boltzmann models, the Cartesian components of the
vectors {ei} satisfy the following relations [16, 18, 20, 24]:∑
i
wieiα = 0 (2.45)
∑
i
wieiαeiβ = χ c2 δαβ (2.46)
∑
i
wieiαeiβeiγ = 0 (2.47)
∑
i
wieiαeiβeiγeiδ = χ2 c4 (δαβδγδ + δγβδαδ + δδβδγα) (2.48)
Using the above relations, we can easily compute the sums:∑
i
fi = ρ (2.49)
∑
i
eiα fi = ρ uα (2.50)
∑
i
eiα eiβ fi = ρ
[
χ c2 δαβ + uα uβ
]
(2.51)
∑
i
eiα eiβ eiγ fi = χ c2 ρ [ δαβuγ + δαγuβ + δβγuα ] (2.52)
2.4 The Lattice Boltzmann algorithm
The evolution of the distribution functions in Lattice Boltzmann models is gov-
erned by the discretized equation (2.39). The value of the distribution function
fi(x, t+ δt) at time t+ δt can be evaluated from the Taylor expansion series:
fi(x, t+ δt) = fi(x, t) +
δt
1!
·
∂fi(x, t)
∂t
+ (δt)
2
2!
·
∂(fi(x, t))2
∂t2
+ . . . (2.53)
If we neglect the high order terms, we get:
∂fi(x, t)
∂t
= fi(x, t+ δt)− fi(x, t)
δt
(2.54)
In the relation (2.39), we replace Eq. (2.54) to get:
fi(x, t+ δt) = fi(x, t)− δteiα∂αfi(x, t) −
δt
τ
[ fi(x, t) − f eqi (x, t) ]
+ δt
kBT
[ eiα − uα(x, t) ] Fα f eqi (x, t) (2.55)
i = 0, 1, . . . N
12