Deriving Navier Stokes with particle methods.

(1)

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.

(2)

ii

(3)

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

(6)

vi CONTENTS

(7)

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 equations. 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 macroscopic properties to derive the Euler equations. The Navier Stokes equations are derived using perturbation theory. Lastly we discuss the numerical method itself and how to implement it along with various boundary conditions. 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

(8)

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 ,

(9)

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 + ∇ · σ⁰+ F , ρDe

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

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

σ⁰ = 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

(10)

4 CHAPTER 1. INTRODUCTION AND PRELIMINARIES Let ξ(x, t) be a solution to the Navier Stokes equations in the fixed coordinate 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¹ 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.

(11)

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,

σ_αβ⁰ = µ ∂u_β

∂xα

+ ∂u_α

∂xβ

− δ_αβ2 3

∂u_γ

∂xγ

, q_α = −κ∂T

∂xα

. (1.4)

(12)

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_α +∂σ_αβ⁰

∂x_β + Fα,

∂ρe

∂t +∂ρuαe

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

∂x_α + σ_αβ⁰ ∂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|²f dx.

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

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

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

|x|⁴f dx.

This can be generalized for any 2n repeated x_ix_j’s, Z

x1x2...xn2−1x2n

| {z }

2n times

f dx = 1

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

Z

|x|²ⁿ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)!! = Πⁿ_i=1(2i − 1) = (2n)!

2ⁿ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,

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

(13)

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²ⁿexp −ax² dx These values can easily² be determined to equal,

Z ∞ 0

x²ⁿexp −ax² dx =r π a

(2n − 1)!!

aⁿ2ⁿ⁺¹ . (1.7)

2This can be done by using for example parameter differentiation

(14)

8 CHAPTER 1. INTRODUCTION AND PRELIMINARIES

(15)

Chapter 2 The Lattice Boltzmann Method

Particle Methods such as the Lattice Boltzmann method work on a microscopic 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 position 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²³atoms! Simulating all of these atoms would mean that we need to store 6 ∗ 10²³ 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 directly 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

(16)

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 distribution function f . We then derive the Boltzmann Equation and its equilibrium 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 |ξ|² as a weight.

ρE(x, t) = 1 2

Z

|ξ|²f (x, ξ, t)dξ. (2.3) We will use these moments later on in section 2.1.3 to derive the conservation 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)

(17)

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|²+ v · u + |u|²)f (x, ξ, t)dξ

= 1 2

Z

|v|²f dξ +1 2

Z

|u|²f dξ

= ρe(x, t) + 1 2|u|²ρ.

(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|²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|²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)

(18)

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

|ξ|²Ω(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 implement in a numerical algorithm. The BGK operator uses the assumption that the distribution function tends to relax towards its equilibrium function. 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 equilibrium 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

(19)

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²^x, f_y^(eq)(vy) = ae^−bv²^y, f_z^(eq)(vz) = ae^−bv²^z

f^(eq)(|v|) = ae^−b|v²^|.

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²^|dξ = a Z

e^−b(v^x²^+v^y²^+v²^z⁾dξ

= a Z ∞

−∞

e^−bv²^xdvx

Z ∞

−∞

e^−bv²^ydvy

Z ∞

−∞

e^−bv²^zdvz = a

π b

3/2

Thus we get,

a = ρ b π

3/2

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

Z

|v|²f^(eq)(|v|)dξ = ρ b π

3/2Z

|v|²e^−b|v²^|dξ

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

ρe = ρ 2

b π

3/2Z ∞ 0

|v|²e^−b|v²^|4π|v|²d|v|

= ρ2π b π

3/2Z ∞ 0

|v|⁴e^−b|v²^|d|v|

This last integral can be found using Gaussian integrals and equals³₈√

π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|² e

. (2.13)

(20)

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|²ⁿf^(eq)dξ = Z

4π|v|²ⁿ⁺²f^(eq)d|v|² Substituting (2.12) and using (1.7) gives,

4πρ b π

³₂ Z

|v|²ⁿ⁺²e^−b|v²^|d|v|² = 4πρ b π

³₂r π b

(2n + 1)!!

(2b)ⁿ 1 2²b After canceling common terms we find,

Z

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

= (2n + 1)!! pⁿ ρⁿ⁻¹

(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|²ⁿf^(eq)d|v|²

Substituting (2.14), Z

v₁v₂...v_2n−1v_2n

| {z }

2n times

f^(eq)dξ = ρ

(2b)ⁿδ1,2,...2n−1,2n

= pⁿ

ρⁿ⁻¹δ1,2,...2n−1,2n.

(2.15)

We will use these expressions to derive the Navier Stokes equations in section 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)

(21)

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 operator it is much simpler to compute. The conservation properties of the BGK operator are easily derived by using the density, momentum and energy moments 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 operator (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.

(22)

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

∂ρ

∂t +∂ρu_α

∂x_α = 0 Conservation of momentum

In a similar way we will take the first momentum over the Boltzmann equation 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)

(23)

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|² Z

f dξ + u_α Z

|v|²f dξ + 2u_β Z

v_αv_βf dξ + Z

v_α|v|²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|²f dξ (2.23)

To get, Z

ξ_αξ_αξ_βf dξ = u_α|u|²ρ + 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|²

∂t + ∂ρu_α|u|²

∂xα

= ∂u_βσ_αβ

∂xβ

+ Fβuβ− ∂q_α

∂xβ

.

(24)

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|²

∂t +∂ρu_α|u|²

∂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|²

∂t +∂ρu_α|u|²

∂xα

= uα

∂ρu_α

∂t +∂ρu_αu_β

∂xβ

, ρu_α ∂u_α

∂t + u_β∂u_β

∂xα

+1

2|u|² ∂ρ

∂t +∂ρu_α

∂xα

= ρu_α ∂u_α

∂t + u_β∂u_α

∂xβ

+ |u|² ∂ρ

∂t +∂ρu_β

∂xβ

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

1

2|u|² ∂ρ

∂t +∂ρuα

∂x_α

= 0

which is valid due to the continuity equation. This leaves us with the conservation 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α

.

(25)

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⁽ⁿ⁾ 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

˜

τ = τ_x^ξ⁰

mfp

F = F˜ _ρξ^x⁰

0

f = f˜ ^c_ρ³⁰

Here we introduced the mean free path xmfp, mean free time tmfp of a particle 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₀ (∂ ˜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⁽ⁿ⁾(where f⁽⁰⁾= f^(eq)),

f = f^(eq)+ f⁽¹⁾+ ²f⁽²⁾+ ...

= f^(eq)+ f⁽¹⁾+ O(²)

(26)

20 CHAPTER 2. THE LATTICE BOLTZMANN METHOD The ⁿ terms describes the order on which the function f⁽ⁿ⁾ has an effect on the distribution function. Mathematically this can be written as,

f⁽¹⁾/f^(eq)= O(Kn), f⁽ⁿ⁾/f^(eq)= O(Knⁿ).

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⁽¹⁾+ O(²)) = 1

τ (f⁽¹⁾+ O(²)) The extra ⁻¹ 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⁽¹⁾

τ (2.27)

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

σ = σ⁽⁰⁾+ σ⁽¹⁾+ O(²) q = q⁽⁰⁾+ q⁽¹⁾+ O(²) with,

σ_αβ⁽ⁿ⁾= − Z

vαv_βf⁽ⁿ⁾dξ, q⁽ⁿ⁾_α = 1 2

Z

vα|v|²f⁽ⁿ⁾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⁽ⁿ⁾ 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⁽¹⁾ 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

(27)

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

f⁽¹⁾

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⁽¹⁾ by,

f⁽¹⁾= −τ 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⁽¹⁾.

In the remaining part of this section we will derive the partial derivatives 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|² e

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

4π

− 3

2ln e − 3 4e|v|²

Although we have so for written the equilibrium distribution only in terms of the peculiar velocity, it is also dependent on the macroscopic variables ρ, 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_β

∂t +∂ ln f^(eq)

∂e

∂t, (2.29a)

∂ ln f^(eq)

∂xα

= ∂ ln f^(eq)

∂ρ

∂xα

+∂ ln f^(eq)

∂uβ

∂u_β

∂xα

+∂ ln f^(eq)

∂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⁽¹⁾.

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

(28)

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|²

=

3

4e²|v|²− 3 2e

,

Plugging these into (2.30) results in,

∂ ln f^(eq)

∂t = 1 ρ

∂ρ

∂t −3v_β 2e

∂u_β

∂t +

3

4e²|v|²− 3 2e

∂e

∂t, (2.30a)

∂ ln f^(eq)

∂xα

= 1 ρ

∂ρ

∂xα

−3v_β 2e

∂u_β

∂xα

+

3

4e²|v|²− 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⁽¹⁾ = − τ f^(eq)[1 ρ

∂ρ

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

∂x_α

+3v_β

2e

∂u_β

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

∂x_α

+

3

4e²|v|²− 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⁽¹⁾ = − τ f^(eq)[1

ρ

−ρ∂u_α

∂xα

+ v_α ∂ρ

∂xα

+3v_β

2e

1

ρF_β−1 ρ

∂p

∂xβ

+ v_α∂u_β

∂xα

+

3

4e²|v|²− 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⁽¹⁾. Furthermore every tensor term containing only α or β, but not both

(29)

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

f⁽¹⁾= − τ f^(eq)[−∂uα

∂xα

+vα

ρ

∂ρ

∂xα

+ 3

2eρvβFβ−v_β ρ

∂ρ

∂x_β −v_β e

∂e

∂x_β + 3 2evαvβ

∂u_β

∂xα

+

3

4e²|v|²− 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²|v|²− 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|² 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⁽¹⁾ = −τ f^(eq) 3|v|² 4e² − 5

2e

v_α ∂e

∂x_α + 3 2e

v_αv_β∂uβ

∂x_α −|v|² 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.

σ_αβ⁽¹⁾= − Z

v_αv_βf⁽¹⁾dξ

= Z

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

2e

v_γ ∂e

∂x_γ + 3 2e

v_γv_δ∂u_δ

∂x_γ −|v|² 3

∂u_γ

∂x_γ

dξ

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

σ_αβ⁽¹⁾ τ = ∂e

∂x_γ

Z 3|v|² 4e² − 5

2e

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

Z

v_γv_δ∂u_δ

∂x_γ −|v|² 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,

σ⁽¹⁾_αβ = 3τ 2e

Z

vγv_δ∂uδ

∂x_γ − |v|² 3

∂uγ

∂x_γ

vαv_βf^(eq)dξ

= 3τ 2e

∂u_δ

∂xγ

Z

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

∂xγ

Z |v|²

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

Deriving Navier Stokes with particle methods.