• No results found

Deriving Navier Stokes with particle methods.

N/A
N/A
Protected

Academic year: 2021

Share "Deriving Navier Stokes with particle methods."

Copied!
53
0
0

Bezig met laden.... (Bekijk nu de volledige tekst)

Hele tekst

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

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

(4)

iv CONTENTS

(5)

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

(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],

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

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

(10)

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

Fz

= m

 ax

ay

az

= 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 convention1 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,

σαβ0 = µ ∂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α +∂σαβ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|2f dx.

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

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

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

|x|4f dx.

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

x1x2...xn2−1x2n

| {z }

2n times

f dx = 1

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

Z

|x|2nf 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)!! = Πni=1(2i − 1) = (2n)!

2nn!.

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

x2nexp −ax2 dx These values can easily2 be determined to equal,

Z 0

x2nexp −ax2 dx =r π a

(2n − 1)!!

an2n+1 . (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 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 1023atoms! Simulating all of these atoms would mean that we need to store 6 ∗ 1023 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

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

|ξ|2f (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)

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

= 1 2

Z

|v|2f dξ +1 2

Z

|u|2f 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|2f 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|2f 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

∂ξα

α

dt .

Some of these terms can be replaced by using dxdtα = ξα and 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

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

(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|) = fx(eq)(vx)fy(eq)(vy)fz(eq)(vz)

This relationship is satisfied by Gaussian functions, thus we write, fx(eq)(vx) = ae−bv2x, fy(eq)(vy) = ae−bv2y, fz(eq)(vz) = ae−bv2z

f(eq)(|v|) = ae−b|v2|.

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

e−b(vx2+vy2+v2z)

= a Z

−∞

e−bv2xdvx

Z

−∞

e−bv2ydvy

Z

−∞

e−bv2zdvz = a

π b

3/2

Thus we get,

a = ρ b π

3/2

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

Z

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

3/2Z

|v|2e−b|v2|

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|2d|v|,

ρe = ρ 2

 b π

3/2Z 0

|v|2e−b|v2|4π|v|2d|v|

= ρ2π b π

3/2Z 0

|v|4e−b|v2|d|v|

This last integral can be found using Gaussian integrals and equals38

π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)

(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|2nf(eq)dξ = Z

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

4πρ b π

32 Z

|v|2n+2e−b|v2|d|v|2 = 4πρ b π

32r π b

(2n + 1)!!

(2b)n 1 22b After canceling common terms we find,

Z

|v|2nf(eq)dξ = (2n + 1)!! ρ (2b)n

= (2n + 1)!! pn ρ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|2nf(eq)d|v|2

Substituting (2.14), Z

v1v2...v2n−1v2n

| {z }

2n times

f(eq)dξ = ρ

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

= pn

ρ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)

(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 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.

(22)

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)

(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

β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|2f dξ + 2uβ Z

vαvβf dξ + Z

vα|v|2f 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|2f 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β

.

(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|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α

.

(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(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 = x00 is used to replace the dimensional variables into new dimensionless variables,

˜ x = xx

0

ξ =˜ ξξ

0

˜t = txξ0

0

˜

τ = τxξ0

mfp

F = F˜ ρξx0

0

f = f˜ cρ30

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,

xmfp x0 (∂ ˜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 xmfp/x0 or equivalently tmfp/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)+ 2f(2)+ ...

= f(eq)+ f(1)+ O(2)

(26)

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(Knn).

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|2f(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

(27)

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



− 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

(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|2



=

 3

4e2|v|2− 3 2e

 ,

Plugging these into (2.30) results in,

∂ ln f(eq)

∂t = 1 ρ

∂ρ

∂t −3vβ 2e

∂uβ

∂t +

 3

4e2|v|2− 3 2e

 ∂e

∂t, (2.30a)

∂ ln f(eq)

∂xα

= 1 ρ

∂ρ

∂xα

−3vβ 2e

∂uβ

∂xα

+

 3

4e2|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

4e2|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

4e2|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

(29)

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

4e2|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

4e2|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 4e2 − 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)

= Z

vαvβτ f(eq) 3|v|2 4e2 − 5

2e

 vγ ∂e

∂xγ + 3 2e



vγvδ∂uδ

∂xγ −|v|2 3

∂uγ

∂xγ



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

σαβ(1) τ = ∂e

∂xγ

Z  3|v|2 4e2 − 5

2e



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

Z 

vγvδ∂uδ

∂xγ −|v|2 3

∂uγ

∂xγ



vαvβf(eq)

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)

= 3τ 2e

 ∂uδ

∂xγ

Z

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

∂xγ

Z |v|2

3 vαvβf(eq)



Referenties

GERELATEERDE DOCUMENTEN

Magnetic specific heat of the nearly-one-dimensional system tetramethyl ammonium nickel trichloride (TMNC).. Citation for published

These sign types were then compared with counterparts in six potential lexifier sign languages, American Sign Language (ASL), British Sign Language (BSL), Irish Sign Language (ISL),

The most significant differences in TB-associated obstructive pulmonary disease (TOPD) were more evidence of static gas (air) trapping on lung function testing and quantitative

Where traditionally a common scenario would be that the residential electricity system would simply consume energy from the grid and then be billed for this energy,

The use of a lower solids loading of 13% solids after water-impregnated steam pretreatment resulted in an ethanol concentration of 39 g/L, which was very close

transforr.tations for non-active plans have been considered or.. it is not worthwile to transform these plans any further. Ue will discuss now the behaviour of

Paulette Talliard NDLTD / ETD/ Workshop 16/ 17 September 2011 Stellenbosch University, Stellenbosch, South Africa.. Populating an

Resultantly, researchers and readers will only know whether or not the event had a significant impact over an entire event window that could potentially be as long as months,