faculteit Wiskunde en Natuurwetenschappen

## Deriving Navier Stokes with particle methods.

### Bacheloronderzoek Technische Wiskunde

Juni 2014

Student: M. S. Redeman

Eerste Begeleider: Prof. Dr. Ir. R. W. C. P. Verstappen Tweede Begeleider: Prof. Dr. E. C. Wit.

ii

### Contents

Abstract v

1 Introduction and preliminaries 1

1.1 Fluid Dynamics . . . 2

1.1.1 The Euler model . . . 2

1.1.2 Navier Stokes model . . . 2

1.1.3 Galilean Invariance . . . 3

1.2 Index notation . . . 4

2 The Lattice Boltzmann Method 9 2.1 Deriving the Navier Stokes equations . . . 10

2.1.1 Moments . . . 10

2.1.2 The Boltzmann Equation and the Maxwell Boltzmann distribution . . . 11

2.1.3 Conservation laws . . . 15

2.1.4 Deriving the Navier Stokes equations . . . 19

2.1.5 Galilean invariance and Boltzmann’s H-theorem . . . 25

2.2 The Lattice Boltzmann method . . . 26

2.2.1 Discretizing the Boltzmann equation . . . 27

2.2.2 Accuracy of the Lattice Boltzmann method . . . 29

2.2.3 Boundary conditions . . . 36

2.3 Summary . . . 40

3 Multiparticle Collision Dynamics 41 3.1 The Multiparticle Collision Dynamis method . . . 41

3.2 Deriving the Navier Stokes equations . . . 44

3.3 Boundary Conditions . . . 44

4 Concluding remarks 45

iii

iv CONTENTS

### Abstract

Simulating fluid flow can be very complex. We can model fluids using the Navier Stokes equations, which can be solved numerically using for example finite-difference or finite-elements. But these methods are generally hard to implement on a computer. Particle methods like the Lattice Boltzmann method and multiparticle collision dynamics can easily be implemented as a numerical model, but have the disadvantage that it hard to prove that they model the Navier Stokes equations.

This thesis uses kinetic theory and perturbation theory to derive the Navier Stokes equations from the continuous and discrete Boltzmann equation. A numerical implementation of the Lattice Boltzmann method is discussed along with basic boundary conditions. Lastly we give an outline on how to do a similar derivation for the multiparticle collision method.

v

vi CONTENTS

### Chapter 1

### Introduction and preliminaries

Fluid flow problems are hard to solve analytically. We therefore want to find fast and accurate numerical methods that can solve these problems. This bachelor thesis discusses two such methods: the Lattice Boltzmann Method and the multiparticle collision cynamics method. Most numerical methods such as finite-difference and finite-elements work on a macroscopic scale and solve these problems by discretizing the Euler or Navier Stokes equations.

The methods discussed in this thesis work on a microscopic scale and can be categorized as Particle Methods. A particle method models fluids by simulate the movement and collision of every individual particle in a fluid.

Implementing these methods is easy because we only need to keep track the movement of every particle and change its velocity in case of collisions.

Moreover these methods are well suited for parallel computing because upon updating a particle we only need to know its velocity and the trajectory of its neighbours. It is however hard to analytically prove that these methods behave according to the Euler and Navier Stokes equations.

In the following subsections we discuss the Euler and Navier Stokes equa- tions. Thereafter we introduce the index notation that is used throughout this thesis. Both the Euler and the Navier Stokes equations are rewritten in this index notation. Lastly we introduce useful properties of weighted integrals of a spherical integrand.

This thesis’ main focus is chapter 2 where we discuss the Lattice Boltz- mann method. After deriving the Boltzmann equation we use its macro- scopic properties to derive the Euler equations. The Navier Stokes equa- tions are derived using perturbation theory. Lastly we discuss the numerical method itself and how to implement it along with various boundary condi- tions. The derivation of the Lattice Boltzmann method is mostly based on Erlend Magnus Viggen’s PHD thesis [10].

After discussing the Lattice Boltzmann Method we will introduce the 1

2 CHAPTER 1. INTRODUCTION AND PRELIMINARIES multiparticle collision (MPC) dynamics method in chapter 3. It is shown that the MPC method satisfies the conservation of mass, momentum and energy. Thereafter an outline is given on how to show that the MPC method

### 1.1 Fluid Dynamics

When modeling a fluid it is typically seen as a continuum where we assume that it can be seen as a continuous body and disregard the fact that it is actually made of individual atoms. By stating that it must conserve mass, momentum and energy we can derive the Euler equations. We describe a fluid by the macroscopic variables of fluid density ρ, fluid velocity u and the pressure p.

1.1.1 The Euler model

The Euler model was originally published by Leonhard Euler in 1757. It consists of a set of three equations which govern the conservation of mass, momentum and energy.

We express these equations using the material derivative, Dλ(x(t), t)

Dt = ∂λ(x(t), t)

∂t + (λ · ∇)λ, the equations are then written as [9],

Dρ

Dt = −ρ∇ · u, ρDu

Dt = −∇p + F , ρDe

Dt = −p∇ · u.

(1.1a) (1.1b) (1.1c) Here the external body force density is denoted as F and the internal energy is e. The first equation is known as the continuity equation. The second equation is the momentum equation and describes how a particle’s velocity is changed with respect to the change of pressure and the external body force. The Euler equations can be We will now see that the Euler equations are a simplified version of the Navier Stokes equations.

1.1.2 Navier Stokes model

The Navier Stokes model uses the same mass conservation equation, but as- sumes that the stress is influenced by a diffusing viscous term and a pressure term. The momentum equation is of the form,

ρDu

Dt = ∇ · σ + F ,

1.1. FLUID DYNAMICS 3 called the Cauchy momentum equation with Cauchy stress tensor σ. The Navier Stokes equations are

Dρ

Dt = −ρ∇ · u, ρDu

Dt = −∇p + ∇ · σ^{0}+ F ,
ρDe

Dt = −p∇ · u + σ_{αβ}∇ · u − ∇ · q.

(1.2a)
(1.2b)
(1.2c)
The Cauchy stress tensor is then defined as, σ = −pI + σ^{0} with σ^{0} being
the deviatoric stress tensor was found to equal,

σ^{0} = 2

3µ(∇ · u)I + µ(∇u + (∇u)^{T}).

The deviatoric stress tensor is responsible for the friction inside the fluid with the parameter µ being the dynamic shear viscosity. Another new quantity was introduced in the energy conservation equations, the heat flux q. Using Fourier’s heat conduction law the heat conduction can be expressed as,

q = −κ∇T

With T being the temperature. Using the ideal gas law we can relate the pressure in a fluid with that of the density, gas constant R and temperature,

p = ρRT

This identity will be used later on in 2.1.4 to derive the Navier Stokes formu- las from the Boltzmann equation. When both the deviatoric stress tensor and the heat flux are zero, then the Navier Stokes equations are reduced to the Euler equations.

1.1.3 Galilean Invariance

Galilean invariance states that laws of motions do not change when we change our inertial reference view. In this section we will proof that the Navier Stokes equations are Galilean invariant. The motivation behind this section is that when deriving numerical methods to solve the Navier Stokes equations we want those solutions to have the same properties as the ana- lytical solutions.

To prove that the Navier Stokes equations are Galilean invariant we consider a fixed coordinate system with time t and position x and a moving coordinate system with velocity V for which the following equations hold,

¯

x = x − V t, ¯t = t

4 CHAPTER 1. INTRODUCTION AND PRELIMINARIES Let ξ(x, t) be a solution to the Navier Stokes equations in the fixed co- ordinate system, then the solution in the moving system is derived as,

¯ξ(¯x, ¯t) = ξ(x, t) − V . We then derive the partial derivatives,

∇¯¯ξ = ∇ξ, ∂¯ξ

∂¯t = ∂ξ

∂t + V · ∇ξ.

From which we can find that the material derivative with respect to the fixed coordinate system is identical to the material derivative with respect to the moving coordinate system,

D¯ξ D¯t = ∂ ¯ξ

∂¯t + ¯ξ · ¯∇ ¯ξ = ∂ξ

∂t + V · ∇ξ + (ξ − V ) · ∇ξ = Dξ Dt

It follows that the velocity gradients and fluid acceleration are Galilean invariant, but the velocity and time are not invariant.

### 1.2 Index notation

Throughout the thesis we use index notations instead of the more commonly used vector notation. This ensures that equations involving higher order tensors are more understandable. Furthermore it is easier to implement computer code of equations written in index notations. We now illustrate the index notation by comparing Newton’s second law of motion written in vector notation and in index notation,

F =

Fx

Fy

F_{z}

= m

ax

ay

a_{z}

= ma The same equation can be written in index notation as,

F_{α}= ma_{α}

here α is a generic index which allows us to write the law of motion for all
three spatial directions, but with only one equation. We now show some
examples that illustrate the index notation. Using Einstein’s summation
convention^{1} the inner product of two vectors are written as,

a · b = axbx+ ayby+ azbz =X

α

aαbα

Divergence can then be written in index notation as,

∇ · ξ(x) = ∂ξα

∂xα

1Repeating an index twice in a single term implies summation over all possible values of that index.

1.2. INDEX NOTATION 5 Kronecker delta

The Kronecker delta symbol is a generalization of the identity and is written as,

δ_{αβ} =

(1 if α = β 0 if α 6= β

Note that we now have two indexes. The amount of indexes imply the tensor
order of an object. A third order tensor R could be written as R_{αβγ}. Using
the summation convention we can change the index of a quantity using the
Kronecker delta,

δ_{αβ}a_{β} = a_{α}.

Both the Euler and the Navier Stokes model can be rewritten into this index notation.

The Euler model

First we can rewrite the material derivative into index notation, Dλ

Dt = ∂

∂t+ uα

∂

∂x_{α}

λ

Thus the Euler equations in index notations are,

∂ρ

∂t + ∂ρuα

∂xα

= 0,

∂ρu_{α}

∂t +∂ρuαuβ

∂x_{β} = − ∂p

∂x_{α} + F_{α},

∂ρe

∂t +∂ρuαe

∂xα

= −p∂uα

∂xα

.

(1.3a) (1.3b) (1.3c)

Here we used the fact that with the continuity equation we have, ρDλ

Dt = ∂ρλ

∂t +∂ρuαλ

∂xα

.

The Navier Stokes model

For the Navier Stokes equations we write the deviatoric stress tensor and heat flux as,

σ_{αβ}^{0} = µ ∂u_{β}

∂xα

+ ∂u_{α}

∂xβ

− δ_{αβ}2
3

∂u_{γ}

∂xγ

, q_{α} = −κ∂T

∂xα

. (1.4)

6 CHAPTER 1. INTRODUCTION AND PRELIMINARIES Secondly the Navier Stokes equations written in index notation are,

∂ρ

∂t + ∂ρuα

∂x_{α} = 0,

∂ρuα

∂t +∂ρu_{α}u_{β}

∂x_{β} = − ∂p

∂x_{α} +∂σ_{αβ}^{0}

∂x_{β} + Fα,

∂ρe

∂t +∂ρuαe

∂x_{α} = −pδ_{αβ}∂u_{β}

∂x_{α} + σ_{αβ}^{0} ∂u_{β}

∂x_{α} − ∂qα

∂x_{α},

(1.5a) (1.5b) (1.5c) We can now see why the index notation is more useful to us than normal vector notation. Though the index notation is less compact than using del operators , it is more easy to see how the derivatives effect the equation.

The following two sections introduces some general properties of integrals which will be used throughout this thesis.

Weighted integral of a spherically symmetric integrand.

Let f (x) be a spherical symmetric function around x = 0. Then any integral
weighted with x_{α}x_{β} integrated over the entire volume of x can be written
as,

Z

x_{α}x_{β}f dx = 1
3δ_{αβ}

Z

|x|^{2}f dx.

When weighting over x_{α}x_{β}x_{γ}x_{δ} we get,
Z

xαx_{β}xγx_{δ}f dx = 1

15(δ_{αβ}δ_{γδ}+ δαγδ_{βδ}+ δ_{αδ}δ_{βγ})
Z

|x|^{4}f dx.

This can be generalized for any 2n repeated x_{i}x_{j}’s,
Z

x1x2...xn2−1x2n

| {z }

2n times

f dx = 1

(2n + 1)!!δ1,2,...2n−1,2n

Z

|x|^{2n}f dx. (1.6)

The (2n + 1)!! term is a double factorial which is defined as the product of all preceding odd integers up to 2n + 1. It can be written as,

(2n − 1)!! = Π^{n}_{i=1}(2i − 1) = (2n)!

2^{n}n!.

For example, 5!! is equal to 5 · 3 · 1 = 15. The delta term in (1.6) is the sum of all distinct permutations of Kronecker deltas. As an example, the delta term for n = 3 is,

δ_{αβγδ} = [δ_{αβ}δ_{γδ}δ_{ζ}+ δ_{αβ}δγδ_{δζ}+ δ_{αβ}δ_{γζ}δ_{δ}+ δαγδ_{βδ}δ_{ζ}+ δαγδ_{β}δ_{δζ}+ δαγδ_{βζ}δ_{δ}+
δαδδβγδζ+ δαδδβδγζ+ δαδδβζδγ+ δαδβγδδζ+ δαδβδδγζ+ δαδβζδγδ+
δ_{αζ}δ_{βγ}δ_{δ}+ δ_{αζ}δ_{βδ}δ_{γ}+ δ_{αζ}δ_{β}δ_{γδ}]

1.2. INDEX NOTATION 7 Gaussian integrals

In later chapters we need to determine the value of Gaussian integrals which are of the form,

Z ∞ 0

x^{2n}exp −ax^{2} dx
These values can easily^{2} be determined to equal,

Z ∞ 0

x^{2n}exp −ax^{2} dx =r π
a

(2n − 1)!!

a^{n}2^{n+1} . (1.7)

2This can be done by using for example parameter differentiation

8 CHAPTER 1. INTRODUCTION AND PRELIMINARIES

### Chapter 2

### The Lattice Boltzmann Method

Particle Methods such as the Lattice Boltzmann method work on a micro-
scopic scale by simulating the movement and collisions of the particles in a
fluid. This is done by creating a large set of particles and updating the posi-
tion of a particle according to its velocity. This means that for every particle
in the fluid we need to know 6 variables, the position x = (x, y, z) and its
velocity ξ = (ξ_{x}, ξ_{y}, ξ_{z}). We now face a problem. In one cubic centimeter of
water there are about 10^{23}atoms! Simulating all of these atoms would mean
that we need to store 6 ∗ 10^{23} variables in our computer’s memory and also
do computations on all of these particles. Even with the power of today’s
computers it will be impossible to do simulations with this many particles.

The Lattice Boltzmann method solves this problem by not working di- rectly on the particles of a fluid, but by using the statistical properties of the probability density function f (x, ξ, t). The distribution function can tell us the odds of finding a particle at a position x, with velocity ξ and at time t.

In the following section we will derive the Navier Stokes equations using this distribution function. We first explain some of the properties of the distribution function in section 2.1.1 after which we derive the boltzmann equation and its distribution function’s equilibrium distribution, called the Maxwell Boltzmann distribution. Thereafter we show that the boltzmann equation admits to all conservation equations on a microscopic scale. We then use Chapman Enskog expansion to show that these equations are also satisfied on a macroscopic scale.

Lastly we derive the Lattice Boltzmann method in 2.2 and show how to implement various Boundary Conditions in 2.2.3.

9

10 CHAPTER 2. THE LATTICE BOLTZMANN METHOD

### 2.1 Deriving the Navier Stokes equations

In this section we will derive both the Euler and the more accurate Navier Stokes equations using the distribution function f (x, ξ, t). In the first sub- section we show the macroscopic properties that must hold for a valid dis- tribution function f . We then derive the Boltzmann Equation and its equi- librium distribution which will later on be used to find the conservation equations. In doing so we assume that the fluid we want to model consists of identical particles meaning that every particle in the fluid has the same mass.

2.1.1 Moments

In the following sections we are interested in the macroscopic properties of the distribution function f . We can find these properties by looking at its moments. We define the moment of f as the integrand over the velocity space dξ weighted with a function of velocity. The mass density ρ(x, t) is defined as the zeroth moment of f ,

ρ(x, t) = Z

f (x, ξ, t)dξ. (2.1)

By using ξ as a weight we get ξf (x, ξ, t)dξ, which is the momentum density,

ρu(x, t) = Z

ξf (x, ξ, t)dξ. (2.2)

Using (2.1) and (2.2) we can find the fluid velocity u, which is the average velocity of the particles in the fluid.

Lastly we can find the energy density by taking |ξ|^{2} as a weight.

ρE(x, t) = 1 2

Z

|ξ|^{2}f (x, ξ, t)dξ. (2.3)
We will use these moments later on in section 2.1.3 to derive the conser-
vation laws. To find the momentum and energy conservation equations we
will need to expand ξ into a microscopic and a macroscopic term.

ξ = v + u (2.4)

Here v is called the peculiar velocity, which is the deviation of the particle velocity ξ from the fluid velocity u. We can think of v being the microscopic information of a particle and u as a macroscopic property of the fluid. It can easily be seen thatR vf dξ = 0.

Z

vf dξ = Z

ξf dξ − Z

uf dξ = ρu − ρu = 0, (2.5)

2.1. DERIVING THE NAVIER STOKES EQUATIONS 11 where we used (2.1) and the fact that u is not dependent on ξ.

It follows that we can describe (2.3) as the sum of a microscopic term and a macroscopic term,

ρE(x, t) = 1 2

Z

(|v|^{2}+ v · u + |u|^{2})f (x, ξ, t)dξ

= 1 2

Z

|v|^{2}f dξ +1
2

Z

|u|^{2}f dξ

= ρe(x, t) + 1
2|u|^{2}ρ.

(2.6)

Here e is the internal energy density of the fluid which satisfies the following moment equation,

ρe(x, t) = 1 2

Z

|v|^{2}f dξ. (2.7)

Pressure

The pressure at position x and time t can be derived by looking at the change of momentum of particles bouncing off a surface. By then dividing by an dt and the area of the surface we find the pressure. A more thorough explanation can be found at [10]. The resulting equation for the pressure can then be found to be,

p(x, t) = 1 3

Z

|v|^{2}f dξ = 2

3ρe (2.8)

From this we see that the pressure is proportional to the internal energy density. Meaning that if the internal energy or the density increases then the pressure will also increase.

2.1.2 The Boltzmann Equation and the Maxwell Boltzmann distribution

We now derive the Boltzmann equation and find its equilibrium distribution, called the Maxwell Boltzmann distribution. The Boltzmann equation will tell us how the distribution function evolves with time.

We have seen that the function is a function of x, ξ and t. We can therefore write its derivative as,

df dt = ∂f

∂t dt dt + ∂f

∂xα

dxα

dt + ∂f

∂ξα

dξα

dt .

Some of these terms can be replaced by using ^{dx}_{dt}^{α} = ξα and ^{dξ}_{dt}^{α} = Fα/ρ
where the last term came from Newton’s second law. We now rewrite the
equation into,

df dt = ∂f

∂t + ξ_{α} ∂f

∂xα

+F_{α}
ρ

∂f

∂ξα

. (2.9)

12 CHAPTER 2. THE LATTICE BOLTZMANN METHOD The differential df /dt of f can be seen as the rate of change of the distribution due to particle collisions. Aside the differential being 0 when there are no collisions, we can not say much about the effect of collisions without knowing more about the fluid’s properties. We do however know that collisions at a time t are dependent on the distribution of the position and velocity of each particle in the fluid. We therefore introduce the collision operator Ω(f ) which is dependent on the distribution function. Substituting the collision operator for the differential df /dt we find what is called the Boltzmann equation.

∂f

∂t + ξ_{α} ∂f

∂x_{α} +F_{α}
ρ

∂f

∂ξ_{α} = Ω(f ). (2.10)

Properties of the collision operator

The collision operator must satisfy the same conservation laws as the Boltz- mann distribution, namely mass, momentum and energy conservation. We can summarize these laws into the following equations,

Mass conservation:

Z

Ω(f )dξ = 0, (2.11a)

Momentum conservation:

Z

ξΩ(f )dξ = 0, (2.11b) Energy conservation:

Z

|ξ|^{2}Ω(f )dξ = 0, (2.11c)
The original collision operator introduced by Boltzmann was hard to
calculate on a computer and had its origin from an assumption he called

¨Stossazahlansatz¨, or ¨molecular chaos assumption¨. Later in 1954 P.L. Bhat- nagar, E.P. Gross and M. Krook introduced the BGK collision operator [2].

This operator has the advantage that it is much simpler to compute and im- plement in a numerical algorithm. The BGK operator uses the assumption that the distribution function tends to relax towards its equilibrium func- tion. In the next section we discuss the Maxwell Boltzmann distribution function. Thereafter we introduce the BGK operator.

The Maxwell Botzmann distribution

When the distribution function reaches an equilibrium, then the particles’

velocities ξ, are evenly distributed about the mean velocity u. Once equi- librium is reached we have f (ξ, x, t) = f (v, x, t), for all x and t. This means that the distribution function f (xi, x, t) in equilibrium form is only dependent on the magnitude of the peculiar velocity v. We can now derive the equilibrium distribution of the Boltzmann equation, called the Maxwell

2.1. DERIVING THE NAVIER STOKES EQUATIONS 13
Boltzmann distribution, f^{(eq)}(|v|) and because the equilibrium distribution
is independent of the on the direction we can write,

f^{(eq)}(|v|) = f_{x}^{(eq)}(vx)f_{y}^{(eq)}(vy)f_{z}^{(eq)}(vz)

This relationship is satisfied by Gaussian functions, thus we write,
f_{x}^{(eq)}(vx) = ae^{−bv}^{2}^{x}, f_{y}^{(eq)}(vy) = ae^{−bv}^{2}^{y}, f_{z}^{(eq)}(vz) = ae^{−bv}^{2}^{z}

f^{(eq)}(|v|) = ae^{−b|v}^{2}^{|}.

The constants a and b are independent of v. We derive these constants using
the density moment and energy moment of f^{(eq)} which is easy as we can use
equation 1.7 to find these integrals.

ρ = Z

f^{(eq)}(|v|)dξ = a
Z

e^{−b|v}^{2}^{|}dξ = a
Z

e^{−b(v}^{x}^{2}^{+v}^{y}^{2}^{+v}^{2}^{z}^{)}dξ

= a Z ∞

−∞

e^{−bv}^{2}^{x}dvx

Z ∞

−∞

e^{−bv}^{2}^{y}dvy

Z ∞

−∞

e^{−bv}^{2}^{z}dvz = a

π b

3/2

Thus we get,

a = ρ b π

3/2

To find the constant b we take the energy moment, ρe =

Z

|v|^{2}f^{(eq)}(|v|)dξ = ρ b
π

3/2Z

|v|^{2}e^{−b|v}^{2}^{|}dξ

We cannot expand this integral into three integrals like we did for the density
moment because there is a new |v|^{2} term within the integral. Instead we
use the spherical symmetry of the equilibrium distribution such that we can
write dξ = 4π|v|^{2}d|v|,

ρe = ρ 2

b π

3/2Z ∞ 0

|v|^{2}e^{−b|v}^{2}^{|}4π|v|^{2}d|v|

= ρ2π b π

3/2Z ∞ 0

|v|^{4}e^{−b|v}^{2}^{|}d|v|

This last integral can be found using Gaussian integrals and equals^{3}_{8}√

πb^{−5/2},
thus we get,

b = 3 4e = ρ

2p, a = ρ

3 4πe

3/2

, (2.12)

resulting in the Maxwell Boltzmann distribution,
f^{(eq)}(|v|) = ρ

3 4πe

3/2

exp

−3 4

|v|^{2}
e

. (2.13)

14 CHAPTER 2. THE LATTICE BOLTZMANN METHOD A general expression of the moments of the equilibrium function The integral in the previous section was found using Gaussian integrals. We now derive a general expressions of the moments of the equilibrium function,

Z

|v|^{2n}f^{(eq)}dξ =
Z

4π|v|^{2n+2}f^{(eq)}d|v|^{2}
Substituting (2.12) and using (1.7) gives,

4πρ b π

^{3}_{2} Z

|v|^{2n+2}e^{−b|v}^{2}^{|}d|v|^{2} = 4πρ b
π

^{3}_{2}r π
b

(2n + 1)!!

(2b)^{n}
1
2^{2}b
After canceling common terms we find,

Z

|v|^{2n}f^{(eq)}dξ = (2n + 1)!! ρ
(2b)^{n}

= (2n + 1)!! p^{n}
ρ^{n−1}

(2.14)

As a side note we can use (1.6) to find integrals of the form, Z

v1v2...v2n−1v2n

| {z }

2n times

f^{(eq)}dξ = δ1,2,...2n−1,2n

(2n + 1)!!

Z

4π|v|^{2n}f^{(eq)}d|v|^{2}

Substituting (2.14), Z

v_{1}v_{2}...v_{2n−1}v_{2n}

| {z }

2n times

f^{(eq)}dξ = ρ

(2b)^{n}δ1,2,...2n−1,2n

= p^{n}

ρ^{n−1}δ1,2,...2n−1,2n.

(2.15)

We will use these expressions to derive the Navier Stokes equations in sec- tion 2.1.4. One might have noticed that we did not determined the odd moments of the equilibrium function. This is because the resulting integrals are antisymmetric and integrate over the whole velocity space. Therefore those moments are all zero.

A simplified collision operator

The BGK uses the assumption that the distribution function relaxes towards the Maxwell Boltzmann distribution. Bhatnagar, Gross and Krook used this assumption to find a operator that models the relaxation directly instead of following the details of every collision like the original collision operator did. The BGK operator can be written as,

Ω(f ) = −1 τ

f − f^{(eq)}

. (2.16)

2.1. DERIVING THE NAVIER STOKES EQUATIONS 15 The relaxation time τ indicates how quickly the distribution relaxes towards its equilibrium. Because we are not dealing with any integrals in this opera- tor it is much simpler to compute. The conservation properties of the BGK operator are easily derived by using the density, momentum and energy mo- ments of both the distribution f and the Maxwell Boltzmann distribution.

We now have enough information to derive the conservation laws from the Boltzmann equation and its collision operator. In the next section we first derive the Euler equations. We will then see that only using traditional techniques is not enough to derive the Navier Stokes equations. In 2.1.4 we will use the Chapman Enskog expansion to derive the Navier Stokes equations.

2.1.3 Conservation laws

We have now discussed enough of the Boltzmann equation to derive the conservation equations with microscopic terms. In this section we use the moments defined in 2.1.1 and the conservation laws of the collision opera- tor (2.11) to find a set of equations representing the conservation of mass, momentum and energy for the Boltzmann distribution function.

The basic idea is that we take the zeroth, first and second moment of (2.10) in order to find the equations for mass, momentum and energy respectively.

Conservation of mass

Taking the first moment of (2.10) means that we integrate both sides over the velocity space.

Z ∂f

∂t + ξα

∂f

∂x_{α} +Fα

ρ

∂f

∂ξ_{α}

dξ =

Z

Ω(f )dξ.

The right hand side is zero due to the conservation of mass of the collision operator and we can move the derivatives of t and x out of the integral as they are not functions of ξ. Furthermore we can see that ξα ∂f

∂xα = ^{∂ξ}_{∂}^{α}^{f}

α such that we get,

∂

∂t Z

f dξ + ∂

∂x_{α}
Z

ξαf dξ + Fα

ρ

Z ∂f

∂ξ_{α}dξ = 0.

We have seen the values for the first and second term in 2.1.1 as 2.1 and 2.2. For the last term we must use Stokes’ multidimensional integration by parts to find its value,

Z ∂f

∂ξα

dξ = Z

∂ξ

f dξ = 0 (2.17)

Here we used the assumption that f becomes zero when the velocity ξ tends to infinity as it is physically impossible for a particle to have infinite velocity.

16 CHAPTER 2. THE LATTICE BOLTZMANN METHOD Inserting these moments into the previous equation results in the conti- nuity equation,

∂ρ

∂t +∂ρu_{α}

∂x_{α} = 0
Conservation of momentum

In a similar way we will take the first momentum over the Boltzmann equa- tion such that we will find the conservation of momentum equation,

Z ξα

∂f

∂t + ξ_{β} ∂f

∂x_{β} +Fβ

ρ

∂f

∂ξ_{β}

dξ =

Z

ξαΩ(f )dξ.

Using the same reasoning as before we get, Z

ξ_{α}∂f

∂tdξ + Z

ξ_{α}ξ_{β} ∂f

∂x_{β}dξ +
Z

ξ_{α}Fβ

ρ

∂f

∂ξ_{β}dξ =0.

∂

∂t Z

ξ_{α}f dξ + ∂

∂x_{β}
Z

ξ_{α}ξ_{β}f dξ + Fβ

ρ Z

ξ_{α} ∂f

∂ξ_{β}dξ =0.

(2.18)

The first term is known to be equal to ρu_{α}, but we do not yet know the second
and the third term. The third term can be found using multidimensional
integration by parts similar to findingR _{∂f}

∂ξαdξ, Z

ξ_{α} ∂f

∂ξβ

dξ =
Z ∂ξ_{α}

∂ξβ

f dξ = −ρδ_{αβ}.
Here we used that ^{∂ξ}_{∂ξ}^{α}

β = δαβ and (2.1). For the second term we use the
peculiar velocity and write ξ_{α}ξ_{β} = (u_{α}+ v_{α})(u_{β}+ v_{β}).

Z

ξαξβf dξ = Z

(uαuβ+ uαvβ+ uβvα+ vαvβ)f dξ

= u_{α}u_{β}
Z

f dξ + u_{α}
Z

v_{β}f dξ + u_{β}
Z

v_{α}f dξ +
Z

v_{α}v_{β}f dξ

= u_{α}u_{β}ρ +
Z

v_{α}v_{β}f dξ.

(2.19) If we now define

σ_{αβ} = −
Z

v_{α}v_{β}f dξ (2.20)

and insert the last equations into (2.18) then we get momentum conservation equation,

∂ρu_{α}

∂t +∂ρu_{α}u_{β}

∂xβ

= ∂σ_{αβ}

∂xβ

+ F_{β}. (2.21)

2.1. DERIVING THE NAVIER STOKES EQUATIONS 17 Here σαβ is the Cauchy stress tensor. Note that although the continuity equation is fully described in the macroscopic properties of f , the momentum equation is not. This is due to how we derived (2.20): it is defined by the distribution function and its peculiar velocity. In the next section we will use Chapman Enskog expansion to derive the Cauchy stress tensor in macroscopic terms of the distribution function.

Conservation of energy

For the last conservation equation we take the second moment over het Boltzmann equation.

Z

ξ_{α}ξ_{α} ∂f

∂t + ξ_{β} ∂f

∂x_{β} +Fβ

ρ

∂f

∂ξ_{β}

dξ =

Z

ξ_{α}ξ_{α}Ω(f )dξ.

∂

∂t Z

ξ_{α}ξ_{α}f dξ + ∂

∂x_{β}
Z

ξ_{α}ξ_{α}ξ_{β}f dξ + Fβ

ρ Z

ξ_{α}ξ_{α}∂f

∂ξ_{β}dξ =0, (2.22)
Similar to the conservation of momentum we need to find the result of
two momentum equations, R ξ_{α}ξα∂f

∂ξβdξ and R ξ_{α}ξαξ_{β}f dξ. The first term
can be found using partial integration,

Z

ξ_{α}ξ_{α} ∂f

∂ξ_{β}dξ = −

Z ∂ξαξα

∂ξ_{β} f dξ = −
Z

2ξ_{β}f dξ = −2ρu_{β}.
For the last term we used (2.2).

The second moment is harder to find. We first rewrite the equation in the same way that we did for momentum equation.

Z

ξ_{α}ξ_{α}ξ_{β}f dξ =
Z

(u_{α}u_{β}u_{β}+ u_{α}v_{β}v_{β}+ 2v_{α}v_{β}u_{β}+ v_{α}v_{β}v_{β})f dξ

= u_{α}|u|^{2}
Z

f dξ + u_{α}
Z

|v|^{2}f dξ + 2u_{β}
Z

v_{α}v_{β}f dξ +
Z

v_{α}|v|^{2}f dξ
Then we substitute the values of (2.1), (2.7) and (2.20) for the first three
integrals and define the heat flux q as,

qα= 1 2

Z

vα|v|^{2}f dξ (2.23)

To get, Z

ξ_{α}ξ_{α}ξ_{β}f dξ = u_{α}|u|^{2}ρ + u_{α}ρe + 2u_{β}σ_{αβ}+ 2q_{α}. (2.24)
Dividing (2.22) by 2 and substituting (2.6) for the first term and the
moments we just found for the other terms gives the energy conservation
equation,

∂ρe

∂t +∂ρu_{α}e

∂xα

+1 2

∂ρ|u|^{2}

∂t + ∂ρu_{α}|u|^{2}

∂xα

= ∂u_{β}σ_{αβ}

∂xβ

+ Fβuβ− ∂q_{α}

∂xβ

.

18 CHAPTER 2. THE LATTICE BOLTZMANN METHOD This equation can be further simplified if we subtract uα × (2.21). After using the product rule the right hand side will only contain the stress tensor and heat flux.

∂ρe

∂t +∂ρu_{α}e

∂x_{α} +1
2

∂ρ|u|^{2}

∂t +∂ρu_{α}|u|^{2}

∂x_{α}

−u_{α} ∂ρu_{α}

∂t +∂ρuαu_{β}

∂x_{β}

= σ_{αβ}∂u_{β}

∂x_{α}−∂q_{α}

∂x_{α}.
Using product rules and the continuity equation we can derive the following
equality,

1 2

∂ρ|u|^{2}

∂t +∂ρu_{α}|u|^{2}

∂xα

= uα

∂ρu_{α}

∂t +∂ρu_{α}u_{β}

∂xβ

,
ρu_{α} ∂u_{α}

∂t + u_{β}∂u_{β}

∂xα

+1

2|u|^{2} ∂ρ

∂t +∂ρu_{α}

∂xα

= ρu_{α} ∂u_{α}

∂t + u_{β}∂u_{α}

∂xβ

+ |u|^{2} ∂ρ

∂t +∂ρu_{β}

∂xβ

After switching indexes on the right hand side for the second and fourth term we find,

1

2|u|^{2} ∂ρ

∂t +∂ρuα

∂x_{α}

= 0

which is valid due to the continuity equation. This leaves us with the con- servation of energy equation.

∂ρe

∂t +∂ρu_{α}e

∂xα

= σαβ

∂u_{β}

∂xα

− ∂q_{α}

∂xα

(2.25)

The Euler model

The Euler equations can now easily be derived by assuming that the Boltz- mann distribution function is approximately equal to the Maxwell Boltz- mann distribution function,

f ≈ f^{(eq)}.

We can find the corresponding stress tensor σ and the heat flux q by deriving
the second and third moment of f^{(eq)}, which were determined in 2.1.2.

σ^{(eq)}_{αβ} = −
Z

vαvβf^{(eq)}dξ = −pδαβ

q^{(eq)}_{α} = 1
2

Z

vαv_{β}vγf^{(eq)}dξ = 0

Noting that if f ≈ f^{(eq)} then σ ≈ σ^{(eq)} and q ≈ q^{(eq)} then we can find
the Euler equations by substituting these into (2.21) and (2.25),

∂ρuα

∂t +∂ρu_{α}u_{β}

∂x_{β} = ∂p

∂x_{β} + Fα,

∂ρe

∂t +∂ρu_{α}e

∂xα

= −p∂u_{α}

∂xα

.

2.1. DERIVING THE NAVIER STOKES EQUATIONS 19 2.1.4 Deriving the Navier Stokes equations

In the previous section we derived the Euler and Navier Stokes equations where some terms were dependent on the microscopic properties of the fluid.

By assuming the distribution function f is approximated by its equilibrium
distribution f^{(eq)} we saw that we could write the Euler equations in terms
of only its macroscopic properties. The goal in this section is to derive the
Navier Stokes equations and also find the stress tensor σ_{αβ} and the heat
conduction q_{α} each only depending on the macroscopic properties of the
fluid. To achieve this goal we first have to introduce new mathematical
trick called the Chapman Enskog expansion. The first part of this section
talks about finding a dimensionless Boltzmann equation and how to apply
the Chapman Enskog expansion on this new equation. The second part
will use the newly derived equations from the Chapman Enskog expansion
to find more detailed description of the distribution function as well as the
stress tensor σαβ and the heat conduction qα.

In the Chapman Enskog expansion we expand the distribution function
f into a sum of perturbed functions f^{(n)} each corresponding to a specific
order of the Knudsen number. This Knudsen number Kn is a dimensionless
number and is defined as the ration between the mean free path of a particle
with respect to a characteristic length. Before we do this however we have to
find the dimensionless Boltzmann equation. For this a characteristic length
x0, characteristic velocity ξ0 and a characteristic time to = x0/ξ0 is used to
replace the dimensional variables into new dimensionless variables,

˜
x = _{x}^{x}

0

ξ =˜ _{ξ}^{ξ}

0

˜t = t_{x}^{ξ}^{0}

0

˜

τ = τ_{x}^{ξ}^{0}

mfp

F = F˜ _{ρξ}^{x}^{0}

0

f = f˜ ^{c}_{ρ}^{3}^{0}

Here we introduced the mean free path xmfp, mean free time tmfp of a par- ticle and also the speed of sound within the fluid. If we now replace the old dimensional variables with their respective dimensionless variables, we obtain the dimensionless Boltzmann equation,

x_{mfp}
x_{0} (∂ ˜f

∂˜t + ˜ξ_{α} ∂ ˜f

∂ ˜x_{α} + ˜F_{α}∂ ˜f

∂ ˜ξα

) = −1

˜

τ( ˜f − ˜f^{(eq)}). (2.26)
The dimensionless Boltzmann equation displays the previously mentioned
Knudsen number Kn defined as x_{mfp}/x0 or equivalently t_{mfp}/t0. When
we let the Knudsen number go to zero then the previous equation implies
that the distribution function approaches the equilibrium function. This
motivates us to explore what happens at different orders of the Knudsen
number. For this we expand the function f into a sum of perturbed functions
f^{(n)}(where f^{(0)}= f^{(eq)}),

f = f^{(eq)}+ f^{(1)}+ ^{2}f^{(2)}+ ...

= f^{(eq)}+ f^{(1)}+ O(^{2})

20 CHAPTER 2. THE LATTICE BOLTZMANN METHOD
The ^{n} terms describes the order on which the function f^{(n)} has an effect
on the distribution function. Mathematically this can be written as,

f^{(1)}/f^{(eq)}= O(Kn), f^{(n)}/f^{(eq)}= O(Kn^{n}).

It makes sense to assume that the different orders of the Knudsen number are all semi-independent, therefore we can derive a set of equations from the Boltzmann equation where each equation represents a different order of the Knudsen number. To find these equations we first write the Boltzmann equation in terms of the newly expanded distribution function,

∂

∂t+ ξα

∂

∂xα

+F_{α}
ρ

∂

∂ξα

(f^{(eq)}+ f^{(1)}+ O(^{2})) = 1

τ (f^{(1)}+ O(^{2}))
The extra ^{−1} comes from the fact that ˜ is of order . The set of equations
derived from this expansion will now be used to finally derive the Navier
Stokes equations from the Boltzmann equation. For this derivation it is
sufficient to look only at the equation resulting from the zeroth order of the
Knudsen number,

O(Kn) : ∂f^{(eq)}

∂t + ξα

∂f^{(eq)}

∂xα

+Fα

ρ

∂f^{(eq)}

∂ξα

= f^{(1)}

τ (2.27)

Furthermore by applying (2.20), (2.22) on the newly perturbated func- tion we find that we can expand both σαβ and qα,

σ = σ^{(0)}+ σ^{(1)}+ O(^{2})
q = q^{(0)}+ q^{(1)}+ O(^{2})
with,

σ_{αβ}^{(n)}= −
Z

vαv_{β}f^{(n)}dξ, q^{(n)}_{α} = 1
2

Z

vα|v|^{2}f^{(n)}dξ. (2.28)
The terms corresponding to the zeroth Knudsen order are the terms we
found in 2.1.3. In the next section we use (2.27) and (2.28) to find the stress
tensor and heat conduction with the first order Knudsen terms.

Analysis at the zeroth Knudsen number

One thing to note is that the density, momentum and energy moments for
f^{(n)} where n ≥ 1 are all zero. This follows because those moments of the
equilibrium function f^{(eq)} are identical to the same moments of f . These
moments will be used to derive an equation for f^{(1)} in terms of f^{(eq)}. It
will however not be enough to just take any of the mentioned moments over
(2.27) as then the right hand side will be zero and thus we will not be able

2.1. DERIVING THE NAVIER STOKES EQUATIONS 21
to find an equation for f^{(1)}. Instead we divide by f^{(eq)} and use a reverse
chain rule,

f^{(1)}

f^{(eq)} = − τ
f^{(eq)}

"

∂f^{(eq)}

∂t + ξα

∂f^{(eq)}

∂x_{α} +Fα

ρ

∂f^{(eq)}

∂ξ_{α}

#

= −τ

"

∂ ln f^{(eq)}

∂t + ξα

∂ ln f^{(eq)}

∂x_{α} +Fα

ρ

∂ ln f^{(eq)}

∂ξ_{α}

#

Once we have found the partial derivatives of ln f^{(eq)} with respect to time,
position and the velocity, we can find the higher order term f^{(1)} by,

f^{(1)}= −τ f^{(eq)}

"

∂ ln f^{(eq)}

∂t + ξα

∂ ln f^{(eq)}

∂x_{α} +Fα

ρ

∂ ln f^{(eq)}

∂ξ_{α}

#

The first order perturbation of σ_{αβ} and qα can then be found by applying
(2.28) to the newly found perturbed function f^{(1)}.

In the remaining part of this section we will derive the partial deriva-
tives of ln f^{(eq)}. Let us first recap on what the equilibrium function and its
logarithm looks like,

f^{(eq)}(|v|) = ρ

3 4πe

3/2

exp

−3 4

|v|^{2}
e

ln f^{(eq)}(|v|) = ln ρ +3
2ln 3

4π

− 3

2ln e − 3
4e|v|^{2}

Although we have so for written the equilibrium distribution only in
terms of the peculiar velocity, it is also dependent on the macroscopic vari-
ables ρ, u and e, all of which are dependent on both position and time. The
partial derivative of ln f^{(eq)} with respect to time and position can therefore
be written as,

∂ ln f^{(eq)}

∂t = ∂ ln f^{(eq)}

∂ρ

∂ρ

∂t +∂ ln f^{(eq)}

∂u_{β}

∂u_{β}

∂t +∂ ln f^{(eq)}

∂e

∂e

∂t, (2.29a)

∂ ln f^{(eq)}

∂xα

= ∂ ln f^{(eq)}

∂ρ

∂ρ

∂xα

+∂ ln f^{(eq)}

∂uβ

∂u_{β}

∂xα

+∂ ln f^{(eq)}

∂e

∂e

∂xα

, (2.29b)
Note that in these equations there is no partial derivative with respect to
the velocity as this term is not dependent on either time or space. Once
we have found the partial derivatives of ln f^{(eq)} with respect to the density,
mean velocity and energy, we can use the Euler conservation equations to
further simplify the expression for f^{(1)}.

The mentioned partial derivatives can easily be derived using a combi- nation of chain rules and product rules. For the derivative with respect to

22 CHAPTER 2. THE LATTICE BOLTZMANN METHOD

uβ we use the equivalence, ξβ = uβ+ vβ.

∂ ln f^{(eq)}

∂ρ = ∂

∂ρln ρ = 1 ρ,

∂ ln f^{(eq)}

∂u_{β} = −3
4e

∂

∂u_{β} (ξ_{α}ξ_{α}− 2ξ_{α}u_{α}+ u_{α}u_{α}) = − 3

4e(−2ξ_{β}+ 2u_{β}) = 3v_{β}
2e ,

∂ ln f^{(eq)}

∂e = ∂

∂e

−3

2ln e − 3
4e|v|^{2}

=

3

4e^{2}|v|^{2}− 3
2e

,

Plugging these into (2.30) results in,

∂ ln f^{(eq)}

∂t = 1 ρ

∂ρ

∂t −3v_{β}
2e

∂u_{β}

∂t +

3

4e^{2}|v|^{2}− 3
2e

∂e

∂t, (2.30a)

∂ ln f^{(eq)}

∂xα

= 1 ρ

∂ρ

∂xα

−3v_{β}
2e

∂u_{β}

∂xα

+

3

4e^{2}|v|^{2}− 3
2e

∂e

∂xα

, (2.30b)

Lastly the partial derivative of ln f^{(eq)}with respect to the velocity ξ_{β} is,

∂ ln f^{(eq)}

∂ξα

= − 3 4e

∂

∂ξα

(ξ_{α}ξ_{β}− 2ξ_{β}u_{β}+ u_{β}u_{β}) = −3

4e(2ξ_{α}− 2u_{α}) = −3v_{α}
2e
Inserting these derivatives into (2.29a) while writing u_{α} + v_{α} for ξ_{α} and
canceling common terms gives,

f^{(1)} = − τ f^{(eq)}[1
ρ

∂ρ

∂t + (uα+ vα) ∂ρ

∂x_{α}

+3v_{β}

2e

∂u_{β}

∂t + (uα+ vα)∂u_{β}

∂x_{α}

+

3

4e^{2}|v|^{2}− 3
2e

∂e

∂t+ (u_{α}+ v_{α}) ∂e

∂x_{α}

− 3

2ρeF_{α}v_{α}]

The equation can be further simplified by inserting the Euler equations (1.3)
f^{(1)} = − τ f^{(eq)}[1

ρ

−ρ∂u_{α}

∂xα

+ v_{α} ∂ρ

∂xα

+3v_{β}

2e

1

ρF_{β}−1
ρ

∂p

∂xβ

+ v_{α}∂u_{β}

∂xα

+

3

4e^{2}|v|^{2}− 3
2e

−p ρ

∂u_{α}

∂x_{α} + v_{α} ∂e

∂x_{α}

− 3

2ρeF_{α}v_{α}]

and by using (2.8) we can rewrite the derivative of p with respect to x_{β},

∂p

∂xβ

= 2 3

∂ρe

∂xβ

= 2 3

e ∂ρ

∂xβ

+ ρ ∂e

∂xβ

We now use this and cancel common terms to simplify the expression for
f^{(1)}. Furthermore every tensor term containing only α or β, but not both

2.1. DERIVING THE NAVIER STOKES EQUATIONS 23 can be written in terms of either α or β.

f^{(1)}= − τ f^{(eq)}[−∂uα

∂xα

+vα

ρ

∂ρ

∂xα

+ 3

2eρvβFβ−v_{β}
ρ

∂ρ

∂x_{β} −v_{β}
e

∂e

∂x_{β} + 3
2evαvβ

∂u_{β}

∂xα

+

3

4e^{2}|v|^{2}− 3
2e

−2e 3

∂uα

∂xα

+ vα

∂e

∂xα

− 3

2ρeFαvα]

= −τ f^{(eq)}[−∂uα

∂x_{α} −v_{β}
e

∂e

∂x_{β} + 3

2evαv_{β}∂u_{β}

∂x_{α} +

3

4e^{2}|v|^{2}− 3
2e

−2e 3

∂uα

∂x_{α} + vα

∂e

∂x_{α}

]

Lastly simplifying the expression and canceling common terms gives,

= −τ f^{(eq)} 3
2evαvβ

∂u_{β}

∂xα

+|v|^{2}
2e

3v_{α}
2e

∂e

∂xα

− ∂uα

∂xα

−5vα

2e

∂e

∂xα

Finally we rearrange terms such that on the left hand side we only have an odd number of occurrences of v, resulting in,

f^{(1)} = −τ f^{(eq)} 3|v|^{2}
4e^{2} − 5

2e

v_{α} ∂e

∂x_{α} + 3
2e

v_{α}v_{β}∂uβ

∂x_{α} −|v|^{2}
3

∂u_{α}

∂x_{α}

Using this perturbed distribution we can find the stress tensor and the heat conduction. The stress tensor can be found by finding the second moment of the perturbed function.

σ_{αβ}^{(1)}= −
Z

v_{α}v_{β}f^{(1)}dξ

= Z

v_{α}v_{β}τ f^{(eq)} 3|v|^{2}
4e^{2} − 5

2e

v_{γ} ∂e

∂x_{γ} + 3
2e

v_{γ}v_{δ}∂u_{δ}

∂x_{γ} −|v|^{2}
3

∂u_{γ}

∂x_{γ}

dξ

This can be split up into two integrals each being either an odd or an even moment,

σ_{αβ}^{(1)}
τ = ∂e

∂x_{γ}

Z 3|v|^{2}
4e^{2} − 5

2e

v_{γ}v_{α}v_{β}f^{(eq)}dξ + 3
2e

Z

v_{γ}v_{δ}∂u_{δ}

∂x_{γ} −|v|^{2}
3

∂uγ

∂x_{γ}

v_{α}v_{β}f^{(eq)}dξ

Because the equilibrium function f^{(eq)} is symmetric its odd moments are
zero. Therefore the first integral will vanish meaning that we will only have
to calculate the integral on the right hand side leaving us with,

σ^{(1)}_{αβ} = 3τ
2e

Z

vγv_{δ}∂uδ

∂x_{γ} − |v|^{2}
3

∂uγ

∂x_{γ}

vαv_{β}f^{(eq)}dξ

= 3τ 2e

∂u_{δ}

∂xγ

Z

vγvδvαvβf^{(eq)}dξ −∂u_{γ}

∂xγ

Z |v|^{2}

3 vαvβf^{(eq)}dξ