Simulations and means of characterization of gravity-stabilized flow on a self-affine surface

characterization of gravity-stabilized flow on a self-affine surface

by

Beatrice Baldelli

Department of Physics

University of Oslo

Thesis submitted for the degree of Master of Science June 15, 2019

Self-affine surfaces are found in many physical systems, making it important to study flow over such surfaces.

This project develops a code to simulate, through the coupled Lattice Boltzmann Method, gravity-stabilized flow, with lighter fluid flowing in an system initially filled with heavier fluid, over a self-affine surface of Hurst exponent,H = 0.8.

Methods to qualitatively and quantitatively characterize the flow are introduced. To- gether with the expected stabilizing effect of gravity over long time scales, gravity is found to increase the initial flushing of the heavier fluid. Some explanations are given for this effect.

The rate of dissipation is presented as a means of characterizing the magnitude of bouyancy in the system and of comparing simulations to possible experiments.

3

1 Introduction 7

2 Theory 11

2.1 Fluid dynamics . . . 11

2.2 Lattice Boltzmann equation . . . 14

2.2.1 Kinetic theory . . . 14

2.2.2 BBGKY hierarchy . . . 16

2.2.3 Discretization in velocity space . . . 19

2.2.4 Discretization in time and space . . . 24

2.3 Lattice Boltzmann method . . . 27

2.3.1 D2Q9 . . . 27

2.3.2 Core time step algorithm . . . 28

2.3.3 Boundary conditions . . . 29

2.3.4 Advection-Diffusion LBM . . . 32

3 Method 35 3.1 Code outline . . . 35

3.2 Forcing . . . 36

3.3 Boundary conditions . . . 37

3.4 Self-affine surface . . . 38

3.5 Helmholtz minimum dissipation theorem . . . 39

3.6 Benchmark problem: Poiseuille flow . . . 42

4 Results and discussion 45 4.1 Helmholtz minimum dissipation theorem . . . 46

4.2 Scalar field . . . 47

4.2.1 Advection over diffusion . . . 48

4.3 Gravity . . . 50

4.4 Ideas for future development . . . 60

5 Conclusion 61

6 Bibliography 63

5

Introduction

Many naturally occurring surfaces on the planet, from landscapes such as mountain ranges and flatlands [1] to underground fractures [2, 3], are self-affine surfaces. It is not then surprising that much effort has been made to study the geometry of such sur- faces [4]. Many studies have also been conducted on the flow of fluids over a self-affine surface, or between two such surfaces in the case of fractures [5, 6, 7]. The study of the geometry of self-affine surfaces has led to the development of numerical methods to generate discrete approximations of said surfaces [8]. This is obviously a requirement for the study of flow over these natural surfaces via numerical simulations. Interestingly, these methods are also used in the entertainment industry to generate realistic-looking landscapes and planets.

The flow has also been studied in the presence of varying temperature [3] or concentra- tion of a solute [9, 10]. The two situations are, from a theoretical and computational point of view, similar. Temperature and solute concentration are represented by a scalar field which is advected by the fluid flow and follows a diffusion term proportional to the laplacian of the field itself. They both satisfy the advection-diffusion equation. In the studies that have been carried out [11, 9, 10], typically, heavier fluid is introduced through the inlet in a system filled with lighter fluid. Cold fluid, for example, is in- jected into a fracture whose surfaces are kept at a higher constant temperature [3], or the system is initially filled with fluid free of solute and a certain amount of solute is released at the inlet [10].

In this thesis we want to develop the tools to study the opposite process: lighter fluid flowing in a system initially filled with heavier fluid. The aim is to create a code which allows us to investigate the stabilizing effect of gravity, both for flow with varying tem- perature and flow in the presence of solute. The main difference between the two systems is given by the boundary conditions on the bottom surface. For the flow in the presence of solute, the surface is set to be impermeable to both fluid and solute, while, in the case of thermal flow, the surface is still impermeable to fluid, but the temperature is set to a constant value on it. This can be interpreted physically as the surface belonging to a body so large that it acts as a heat reservoir.

We developed the code for a 2Dsystem with a bottom self-affine solid surface (or more 7

precisely, profile) and a flat open top surface, representing, for example, contact with another fluid. For the flow of fluid with solute the system is initially filled with fluid of constant solute concentration, and solute-free fluid flows into the system for the duration of the simulation. This could be, for example, fresh water flushing out salt water. The flow goes from the inlet on the left (through which solute-free fluid enters the system) to the outlet on the right. The solute is heavier than the fluid, in order to achieve the light-heavy fluid process we are interested in. In the case of thermal flow, the system is initially filled with fluid at the same temperature as the one imposed on the bottom surface, and warmer fluid is injected at the inlet.

The self-affine profile is generated with a Fourier transform method [12]. The flow is sim- ulated via the coupled Lattice Boltzmann Method. This method solves two discretized Boltzmann equations which describe the evolution of two one-particle distribution func- tion. One of the Boltzmann equations approximates the Navier-Stokes equation and therefore simulates the fluid flow; the other approximates the advection-diffution equa- tion which controls the behaviour of the solute concentration. The two equations are coupled: the solute concentration is advected with the velocity of the fluid, the fluid velocity is in turn affected by gravity which is solute-concentration dependent. The method relies on the Boussinesq approximation, which holds in the regime of density variations small enough that their effect on inertia can be neglected, while still influenc- ing gravity. By this approximation the density entering the Navier-Stokes equation is that of the fluid in the absence of solute (or at a reference temperature) and the density variations given by the presence of solute (or changes in temperature) influence the flow only through a gravitational forcing. This idea of coupled Lattice Boltzmann equations was first introduced in the 1990s [13, 14] and has since been used in several articles in the field [9, 15].

We want to be able to extract information about the total amount of solute concentra- tion and the average temperature in the system, and to change the system parameters, such as the relative strength of gravity over fluid particles with different solute concen- tration (or temperature). In the simulation this is done by varying a parameter, and in physical terms it means changing the solute so that the difference between solute-free- fluid and solution densities vary, or changing the thermal expansion coefficient.

After a working code was set up, we started looking at the qualitative characteristics of the flow. We first verified that Helmholtz theorem, stating that Stokes flow has the minimum rate of dissipation, held true for our simulations. The rate of dissipation is in fact a means of characterization of the flow in the presence of gravity. The rate of dissipation can be measured experimentally, and therefore, allows the comparison of simulations and experiments. It may be used to quantify buoyancy effects and to compare flows in different systems and over different surfaces.

We then started looking at the qualitative effects gravity has on the flow field, and the quantitative effects on the total solute concentration, or average temperature, as a func- tion of time. Gravity creates pools of high solute-concentration (or low-temperature) fluid in the valleys of the bottom surface, left behind after the solute-free (or high- temperature) fluid has flushed most of the solute (or cold fluid) out of the system. In the initial time steps of the simulation, however, gravity helps flush the solute (or the

cold fluid) at a higher speed than is found in its absence. Possible explanations are given for this process. The code could be used to conduct further studies on the system, such as exploring the parameter space, or changing the geometry of the bottom surface.

This thesis has the following outline. Chapter 2 covers the theoretical background;

specifically the macroscopic description of fluids with the Navier-Stokes equation (§

2.1), the mesoscopic description with the Boltzmann equation and its discretization (§

2.2) and the Lattice Boltzmann Method (§2.3). Chapter 3 goes over the methods used;

here is a description of the system studied in this project and the programs written for the simulations. In chapter 4 the results are presented and discussed; some ideas are given for the future progress of this project.

Theory

Fluids can be described both macroscopically as a continuous media and microscopically as collections of molecules or atoms. In the macroscopic description the characteristic length,l_M, is given by the distances over which macroscopic quantities, such as pressure or temperature, change. The characteristic length of the microscopic description is given by the molecular sizeld. Kinetic theory explores an intermediate, or mesoscopic, scale characterized by the mean free path l_{mf p}. Rather than considering the fluid as a continuous whole or accounting for the full newtonian motion of all particles, it describes the fluid through a distribution function f(x,p, t). Section 2.1 gives a brief account of two core equations of the macroscopic description of fluids: the continuity and Navier- Stokes equations. Section 2.2 delves into the mesoscopic description of fluids, deriving the Boltzmann equation, first in its continuous, then discrete form. Section 2.3 describes how the discrete Boltzmann equation has been adapted into a computational method to simulate fluid behaviour.

2.1 Fluid dynamics

The continuity and Navier-Stokes equation are derived in any fluid dynamics textbook;

this section follows the approach of Landau and Lifshitz [16].

To derive the continuity equation, which expresses the conservation of mass, we start by considering a volume of fluid,V₀. The mass contained in this volume is given by the integral of the mass density over it:

Z

V0

ρdV. (2.1)

The change in mass is expressed by the time derivative of said mass, _∂t^∂ R

V0ρdV. Since we take the volume V0 to be fixed the order of the derivative and the integral can be exchanged. Mass cannot be spontaneously created or destroyed, therefore the only change in mass in V₀ is given by mass flowing in or out through the surface bounding the volume, ∂V0. TakingdS to be an infinitesimal area on∂V0 andnto be the normal

11

to the surface pointing outwards, the mass flowing through ∂V0 is given by I

∂V0

ρu·ndS (2.2)

whereu is the fluid velocity. The divergence theorem can be used to express this as a volume integral, which can then be equated to the time derivative of mass:

Z

V0

∂ρ

∂tdV =− I

∂V0

ρu·ndS=− Z

V0

∇·(ρu)dV. (2.3)

The minus sign is there because fluid flowing out of the volume decreases the amount of mass in the volume and vice versa fluid flowing in through the volume boundary increases the amount of mass in the volume. Equation (2.3) must hold for any volume, which entails that the integrands are equal:

∂ρ

∂t =−∇·(ρu). (2.4)

This is the continuity equation.

The fluids considered in this thesis are incompressible, i.e. their density is constant in time, ^∂ρ_∂t = 0, and space, ∇ρ = 0. The right-hand side of the continuity equation can be written as −∇·(ρu) = −(ρ∇·u+u·∇ρ) so that, for incompressible fluids, the continuity equation takes the simple form

∇·u= 0. (2.5)

The Navier-Stokes equation describes the changes of the fluid momentum. It equates the time derivative of the fluid momentum to the forces acting on the fluid and the dissipation of momentum through friction. The total force with which the rest of the fluid acts on a certain fluid volume, V0, is given by the integral of pressure over the surface of the fluid volume

− I

∂V0

pndS=− Z

V0

∇pdV, (2.6) where the divergence theorem has been applied again to obtain a volume integral. The negative pressure gradient, -∇p, can be interpreted as a force per unit volume acting on the fluid. The motion of the fluid volume is then given by

ρdu

dt =−∇p. (2.7)

Here ^du_dt is the change in velocity of the volume as it moves through the fluid. The changeduof the volume velocity during a time dtis given by two terms: the change of the velocity at a fixed point in space in the timedt, ^∂u_∂tdt, and the difference in velocity between the point in space the volume starts in and the point it travels to during the timedt, (udt)·∇u. The derivative ^du_dt can therefore be written as

du dt = ∂u

∂t +u·∇u (2.8)

In addition to the force given by the rest of the fluid, an external potential may be present (e.g. gravity) which also contributes to the volume velocity change:

∂u

∂t +u·∇u=−1

ρ∇p+g, (2.9)

whereg is the external force per unit mass per unit volume. The last thing to consider is the dissipation of momentum caused by friction between neighbouring fluid elements.

This is taken into account by substituting the pressure with a more general stress tensor,

−pδ_ik → σ_ik=−pδ_ik+σ_ik⁰ , (2.10) where p is the pressure and σ⁰_ik is the viscous stress tensor. The form of σ_ik⁰ is found by making some physical considerations about the conditions it must satisfy and then finding the most general rank-2 tensor which also satisfies these conditions.

There is friction between neighbouring fluid elements when they have different veloc- ities, therefore σ_ik⁰ must depend on the velocity gradients and cannot have terms not containing them. Assuming the gradients to be small, we can take σ⁰_ik to be a linear function of the first derivatives only. The last observation is that there can be no vis- cous stress for a fluid in uniform rotation, so thatσ_ik⁰ must contain only the derivative combinations

∂ui

∂x_k + ∂uk

∂xi

, (2.11)

which vanish when the velocity can be expressed as the tensor product of angular velocity and radius,u=Ω×r.

The most general rank-2 tensor which satisfies all of these conditions is σ_ik⁰ =µ

∂ui

∂x_k +∂uk

∂xi

−2 3δik

∂ul

∂x_l

+ζδik

∂ul

∂x_l, (2.12)

where η and ζ are the dynamic viscosity and second viscosity coefficients respectively, and are both positive and independent of velocity. The viscosity coefficients are often assumed to be constant, so that inserting (2.12) in (2.9) gives

∂u

∂t +u·∇u=−1

ρ∇p+g+µ ρ∆u+

ζ+1

3µ

∇(∇·u). (2.13) This is the Navier-Stokes equation. For incompressible fluids (∇·u= 0) the equation simplifies to

∂u

∂t +u·∇u=−1

ρ∇p+g+ν∆u, (2.14)

whereν = ^µ_ρ is the kinematic viscosity.

2.2 Lattice Boltzmann equation

A system of N particles can be described by the probability distribution function f(x1, ...,xN;p₁, ...,p_N;t) expressing the probability that at time t the particles are at positionsx₁, ...,x_N with momentum p₁, ...,p_N. Since f(x₁, ...,x_N;p₁, ...,p_N;t) is a probability density, its evolution follows a continuity equation (§ 2.2.1). This equation is exact, but quite unmanageable, becausef(x1, ...,xN;p₁, ...,p_N;t) depends on 6N+ 1 variables, whereN is of order O(10²³). The strategy is to transform the one equation into a system of N coupled equations: the BBGKY hierarchy (§ 2.2.2). The first of these equations controls the one-particle distribution functionf1(x1;p₁;t), which gives the probability of finding at time t a particle at position x₁ with momentum p₁. The second equation controls the two-particle distribution function, f₂(x₁,x₂;p₁, ...,p₂;t), and so on up to the N-particle distribution function, fN(x1, ...,xN;p₁, ...,p_N;t). This does not seem to be an improvement, but the advantage is that, under certain assump- tions, it is possible to close first of these equations, which controls the evolution of the one-particle distribution function, thus obtaining the Boltzmann equation. In order to render it suitable for a computational method, the Boltzmann equation is discretized over velocity space (§ 2.2.3) and over space and time (§2.2.4).

The derivation of the continuous Boltzmann equation follows the lectures by D.Tong [17] and the approach of G.M.Kremer [18].

2.2.1 Kinetic theory

The state at time t occupied by a system of N molecules can be specified in phase space, by 3N position coordinates, x_i, and 3N momentum coordinates, p_i. These satisfy Hamilton’s equations:

∂p_i

∂t =−∂H

∂xi

, ∂x_i

∂t = ∂H

∂p_i (2.15)

whereH is the Hamiltonian which we can write in the form H =

N

X

i=1

p²_i 2m +

N

X

i=1

V(xi) +X

i<j

U(xi−xj). (2.16) The first term is the kinetic energy andV is an external potential that is felt equally by all particles. The third term accounts for the interaction between particles with U(xi−xj) being the potential energy between particlesiand j. We are interested not in the evolution of the coordinates of a specific instance of a system, but in the evolution of the distribution functionf(x_i,p_i, t). The distribution function gives the probability that at time t the system is in a state specified by (x_i,p_i). The distribution function being a probability implies that it must satisfy a continuity equation since probability is locally conserved, and that it is normalized to

Z ^N Y

i=1

d³pid³xi

!

f(xi,p_,t) = 1. (2.17)

Probability is locally conserved, which means that if we take a small region of phase space with volumeV, the change of the probability in the volume will be balanced by the flow of probability through the surface S of the region. This is expressed by

Z

V

∂f

∂tdV + I

S

uf ·ndS= 0 (2.18)

where n is the normal to the surface pointing outwards and u is a 6N dimensional vector representing the flow in phase space. We can rewrite the time derivatives as

uα= ∂x^a_i

∂t , i= 1, ..., N; a= 1,2,3;

α=ia= 1, ...,3N (2.19)

in the 3N position directions and u_α= ∂p^a_i

∂t , i= 1, ..., N; a= 1,2,3;

α= 2N+ia= 3N+ 1, ...6N (2.20) in the 3N momentum directions.

Using the divergence theorem, we can transform the surface integral into a volume integral and get

0 = Z

V

∂f

∂tdV + Z

V

∇ ·(uf)dV

= Z

V

∂f

∂tdV +

N

X

i=1

Z

V

∂

∂x_i · ∂xi

∂t f

+ ∂

∂p_i · ∂p_i

∂t f

dV. (2.21)

This equation has to hold for any volume V, which implies that it holds for the inte- grands. Using Hamilton’s equations (2.15), we can write:

0 = ∂f

∂t +

N

X

i=1

∂

∂x_i · ∂H

∂p_if

− ∂

∂p_i · ∂H

∂x_if

= ∂f

∂t +

N

X

i=1

∂f

∂xi

·∂H

∂p_i − ∂f

∂p_i ·∂H

∂xi

(2.22) where in the second step we have used the fact that the Hamiltonian (2.16) doesn’t contain terms that mix momentum and position variables. Equation (2.22) can be written as the total derivative of the distribution function

df dt = ∂f

∂t +

N

X

i=1

∂f

∂xi

·∂x_i

∂t + ∂f

∂p_i ·∂p_i

∂t

= 0 (2.23)

expressing the fact that probability doesn’t change along a trajectory in phase space.

Equation (2.23) is a closed equation for the evolution of the distribution function de- scribing the system. The problem with this is that in a typical system we want to study N is very large, leading to a function of O(10²³) variables.

The path to solve this problem is to go from a function of 6N variables to a system ofN coupled equations forN functions, and try to close the first of these equations through some sensible approximation.

2.2.2 BBGKY hierarchy

Fromf(xi,ξ_i, t) we can derive over distribution functions by integrating other a subset of the particles. For example:

fn(x1, ...,xn;p₁, ...,p_n;t) =

Z ^N Y

i=n+1

dxidp_i

!

f(x1, ...,xN;p₁, ...,p_N;t) (2.24) is the probability fornparticles to have at timetcoordinates in an infinitesimal region around (xi,ξ_i). To obtain an equation for the evolution of fn we integrate equation (2.23) over the otherN −nparticles:

Z ^N Y

i=n+1

dxidp_i





∂f

∂t +

N

X

j=1

∂f

∂xj

·∂x_j

∂t + ∂f

∂p_j ·∂p_j

∂t



= 0. (2.25) In the time-derivative term we can simply switch the integral and the derivative and obtain

∂

∂t Z ^N

Y

i=n+1

dx_idp_if = ∂f_n

∂t . (2.26)

To deal with the terms in the sum we can divide them between the ones involving particles we are integrating over and the other particles. Let us first look at the terms involving position coordinates:

Z ^N Y

i=n+1

dxidp_i





n

X

j=1

∂f

∂xj

·∂xj

∂t +

N

X

j=n+1

∂f

∂xj

·∂xj

∂t





=

n

X

j=1

∂x_j

∂t · ∂

∂xj

Z ^N Y

i=n+1

dx_idp_if

+

N

X

j=n+1

Z dp_j

N

Y

i=n+1 i6=j

dx_idp_i I

S

∂x_j

∂t f·ndS

=

n

X

j=1

∂x_j

∂t ·∂f_n

∂x_j (2.27) where in the first sum we have followed the same procedure as for the time derivative term and in the second sum we have used the divergence theorem. This second term goes to zero: if our system is in a container, the surface S corresponds to the walls of the container and thereforef is zero on the surfaceS; if there is no container and the system is free to explore the whole space, the surfaceS is at infinity andf must go to zero at infinity in order for it to be normalized to one. We look now at the momentum derivative terms:

Z ^N Y

i=n+1

dxidp_i





n

X

j=1

∂f

∂p_j ·∂p_j

∂t +

N

X

j=n+1

∂f

∂p_j ·∂p_j

∂t



 (2.28)

The terms in the second sum can be transformed through the divergence theorem and go to zero by the same arguments discussed above. To discuss the terms in the first sum we insert the Hamiltonian (2.16) in the first of Hamilton’s equations (2.15):

Z ^N Y

i=n+1

dxidp_i









n

X

j=1

∂f

∂p_j ·









∂V(xj)

∂xj

+

N

X

k=1k6=j

∂U(xj−xk)

∂xj

















=

n

X

j=1









∂V(x_j)

∂xj

+

n

X

k=1k6=j

∂U(x_j−x_k)

∂xj









· ∂

∂p_j Z ^N

Y

i=n+1

dx_idp_if

+

n

X

j=1 N

X

k=n+1

Z

dx_kdp_k









∂U(xj−xk)

∂x_j · ∂

∂p_j Z ^N

Y

i=n+1 i6=k

dxidp_if









=

n

X

j=1

∂V(x_j)

∂xj

+

n

X

k=1k6=j

∂U(x_j−x_k)

∂xj









·∂f_n

∂p_j

+ (N −n)

n

X

j=1

Z

dx_kdp_k

∂U(xj−xk)

∂x_j ·∂fn+1

∂p_j

. (2.29) We can now insert all of these results back into (2.25) and find an equation for the evolution of the n-particles distribution function

∂f_n

∂t +

n

X

j=1









∂xj

∂t ·∂f_n

∂x_j −









∂V(xj)

∂x_j +

n

X

k=1k6=j

∂U(xj−x_k)

∂x_j









·∂f_n

∂p_j









= (N−n)

n

X

j=1

Z

dx_kdp_k

∂U(x_j−x_k)

∂xj

·∂f_n+1

∂p_j

(2.30) We now have N equations for the N distribution functions, but the variables we are interested in, to determine the macroscopic behavior of the system can be derived from the one-particle distribution functionf1(x,p, t), as shown below.

For convenience we now switch to a description in terms of velocity coordinatesξ, rather than momentum ones. We therefore take the one-particle distribution function to be f1(x,ξ, t). We also take it to be normalized to the particle mass m, rather than unity.

Z

dxdξf1(x,ξ, t) =m (2.31)

Variables like the average mass density, ρ(x, t), velocity, u(x, t) and kinetic energy density per unit mass,(x, t) are given by the first velocity moments of the one-particle distribution function:

ρ(x, t) = Z

dξf1(x,ξ, t) (2.32) ρu(x, t) =

Z

dξξf₁(x,ξ, t) (2.33) ρ(x, t) =

Z

dξ|ξ−u|²

2 f₁(x,ξ, t) (2.34) The challenge is now to obtain a closed equation for the evolution of the one-particle distribution function: the Boltzmann equation. The equation we have now depends on the two-particle distribution function:

∂f₁

∂t +∂x

∂t ·∂f₁

∂x −∂V(x)

∂x ·∂f₁

∂ξ

| {z }

streaming

= (N−1) Z

dx⁰dξ⁰∂U(x⁰−x)

∂x⁰ ·∂f₂

∂ξ⁰

| {z }

collision

(2.35)

Through some assumptions the collision integral can be simplified to a form that only depends on the one-particle distribution function [17, 18]. Even though it closes the equation for the evolution off₁, it is a rather complicated collision operator to work with.

It has been proposed that other collision operators, J(f), can be used in place of the original Boltzmann collision operator as long as they possess the following properties, which are true for the Boltzmann collision operator [19]:

1. J(f₁) preservers the five collision invariants, Ψ_k(ξ):

Z

dξdxJ(f)Ψk(ξ), k= 1, ...,5 (2.36) where

Ψ_k=







1, k= 1 ξ, k= 2,3,4 ξ², k= 5

(2.37)

2. J(f₁) leads to an evolution off₁ such that it satisfies the H-theorem.

One such collision operator was proposed by Bhatnagar, Gross and Krook [20], hence its name of BGK collision operator.

J_BGK(f₁) =−1

τ(f₁−f₁⁰) (2.38)

whereτ is the relaxation time andf₀ is the local equilibrium distribution.

This is the operator we will use in our simulations. It assumes that we are close enough to the local equilibrium,f0, that it is possible to linearize the collision operator around

f0. We now have BGK approximation to the Boltzmann equation which we now need to discretize in all its variables.

∂f

∂t +ξ·∂f

∂x+F ·∂f

∂ξ =−1

τ f−f⁰

(2.39) where F = −∇V /m is the acceleration given by the external potential and we have omitted the subscript since it is clear we are talking about the one-particle distribu- tion function. When working with numerical methods it is convenient to work with dimensionless equations. It is also a good approach because the law of similarity states that the physics of fluid systems is characterized by the value of certain dimensionless numbers (i.e. Reynolds number, P´eclet number...). If two systems share the same di- mensionless numbers, the physics of one can be recovered from the other through a simple scaling.

We can choose reference quantities for three independent variables; the others will be determined by the relations between the variables. We can choose for example the ref- erence mass, m₀, the reference time, t₀ and the reference velocity, v₀. The reference length will then be l₀ =v₀t₀. We can now write our variables as function of the refer- ence quantities: ξ =v0ξ^∗, t=t0t^∗, ...; where the starred quantities are dimensionless.

The local equilibrium distribution, which is given by the Maxwell distribution, can be written in terms of dimensionless quantities as

f⁰ = ρ^∗ 2πθ^∗exp

−v^∗2 2θ^∗

, (2.40)

whereθ^∗ =θ/θ0, with the characteristic temperature given by θ0=m0v₀²/k_B. 2.2.3 Discretization in velocity space

In order to get the LBE, we need to discretize the Boltzmann equation (2.39) in velocity space, space and time. For the velocity space discretization we will follow the Hermite series expansion as in the Shanet al. [21]. The expansion of the distribution function f(x,ξ, t) can be carried out in general in D-dimensional space, but since we are simu- lating a 2-dimensional fluid system, we will limit the derivation to the case D=2.

The distribution function is expanded in terms of the Hermite polynomials in velocity space ξ:

f(x,ξ, t) =ω(ξ)

∞

X

n=0

1

n!a⁽ⁿ⁾(x, t) H⁽ⁿ⁾(ξ) (2.41) where thenth-order Hermite polynomialH⁽ⁿ⁾ is annth rank symmetric tensor defined as:

H⁽ⁿ⁾(ξ) = (−1)ⁿ

ω(ξ) ∇ⁿ_ξω(ξ) (2.42)

and ω(ξ) is the weight function:

ω(ξ) = 1 2πe^−ξ

2

2 . (2.43)

The tensors H⁽ⁿ⁾ are a generalization of the usual Hermite polynomials to higher dimensions. The first three are:

H⁽⁰⁾(ξ) = 1, (2.44a)

H_i⁽¹⁾(ξ) =ξi (2.44b)

H_ij⁽²⁾(ξ) =ξiξj−δij. (2.44c) Finding an expression forH⁽ⁿ⁺¹⁾ in terms ofH⁽ⁿ⁾ and H⁽ⁿ⁻¹⁾ is straightforward:

H_ii⁽ⁿ⁺¹⁾_1...in = (−1)ⁿ⁺¹

ω(ξ) ∂i1...∂in∂iω(ξ) = (−1)ⁿ⁺¹

ω(ξ) (−1)∂_i1...∂in(ξiω(ξ))

= (−1)ⁿ

ω(ξ) ∂i1...∂in(δiinω(ξ) +ξi∂inω(ξ))

= (−1)ⁿ

ω(ξ) ξ_i∂_i1...∂_inω(ξ) +

n

X

k=1

δ_iik∂_i1...∂_ik−1∂_ik+1...∂_inω(ξ)

!

. (2.45) We can now recognizeH⁽ⁿ⁾ andH⁽ⁿ⁻¹⁾.

H_ii⁽ⁿ⁺¹⁾_1...in =ξiH_i⁽ⁿ⁾_1...in−

n

X

k=1

δii_kH_ii⁽ⁿ⁻¹⁾_{1...ik−1ik+1...in} (2.46) The Hermite polynomials are orthonormal with respect to the inner product defined in terms of the weighing function as hf, gi = R

ωf gdξ. This means that they satisfy the following relation:

Z

ω(ξ)H_i^(m)_1...imH_j⁽ⁿ⁾_1...jn =δ_mnδⁿ_i

1,...,in;P(j1,...,jn) (2.47) whereδ_iⁿ

1,...,in;P(j1,...,jn) is unity if and only if i1, ..., in is a permutation of j1, ..., jn. The coefficients of expansion a⁽ⁿ⁾(x, t) are found through integration on the whole velocity space of the productfH⁽n):

a⁽ⁿ⁾(x, t) = Z

f(x,ξ, t)H⁽ⁿ⁾(ξ). (2.48) Observe that they are tensors of rank n, just like the Hermite polynomials. Looking at the expression for the Hermite polynomials (2.42), it is clear that the expansion coefficients, a⁽ⁿ⁾, are linear combinations of the velocity moments of the distribution function,f. Substituting (2.44) in (2.48) the first few coefficients are found to be:

a⁽⁰⁾ = Z

f dξ (2.49)

a⁽¹⁾_i = Z

f ξ_idξ (2.50)

a⁽²⁾_ij = Z

f(ξiξj−δij)dξ (2.51)

We can approximate the distribution function toNth order by truncating the expansion:

f(x,ξ, t)≈f^N(x,ξ, t) =ω(ξ)

N

X

n=0

1

n!a⁽ⁿ⁾(x, t) H⁽ⁿ⁾(ξ). (2.52) Since we are interested in studying the macroscopic behaviour of the fluid system, we do not need the full detail of the distribution function, but rather we are interested in its first few velocity moments. These give the thermodynamic variables: mass density, ρ(x, t), fluid velocity,u(x, t) and pressure,p(x, t).

ρ(x, t) = Z

f dξ (2.53)

ρu(x, t) = Z

fξdξ (2.54)

Pij = Z

f vivjdv (2.55)

wherev =ξ−u(x, t) is the departure from the average velocity.

The full second moment P_ij is also interesting in that it gives us the stress tensor through:

σ_ij =P_ij −pδ_ij. (2.56)

It is easy to see that we express these macroscopic variables in terms of the expansion coefficients:

ρ=a⁽⁰⁾ (2.57)

ρui=a⁽¹⁾_i (2.58)

P_ij =a⁽²⁾_ij −ρ(u_iu_j−δ_ij) (2.59) In order to find a discretized expression for the expansion coefficients, we truncatef to order N in the formula (2.48) and express the integrand as

f^N(x,ξ, t) H⁽ⁿ⁾(ξ) =ω(ξ)q(x,ξ, t), (2.60) whereq is a polynomial inξ. Let us look at the degree ofq, by inserting the truncated expansion forf in the above expression.

ω(ξ)

N

X

m=0

1

m!a^(m)(x, t) H^(m)(ξ) H⁽ⁿ⁾(ξ) =ω(ξ)q(x,ξ, t) (2.61) Clearly the greatest degreeq can have is 2N, since the largest value that bothmandn is given by the order at which the expansion is truncated,N. We can now use Gaussian quadrature to evaluate thenth coefficient, by choosing the set of discrete values for the

velocityξ_a, a= 1, .., dsuch that:

a⁽ⁿ⁾= Z

ω(ξ)q(x,ξ, t) =

d

X

a=1

w_aq(x,ξ_a, t)

=

d

X

a=1

wa

ω(ξ_a)f^N(x,ξ_a, t) H⁽ⁿ⁾(ξ_a) (2.62) wherewa are the weights of the quadrature.

The fundamental theorem of Gaussian quadrature states that the optimal discrete ve- locity values are the roots of thenth corresponding orthogonal polynomial [22], which in our case is thenth Hermite polynomial, and that the weights are given by

wa= hH⁽ⁿ⁻¹⁾, H⁽ⁿ⁻¹⁾i

H⁽ⁿ⁻¹⁾(ξ) H⁰⁽ⁿ⁾(ξ), (2.63)

where

H⁰⁽ⁿ⁾= dH⁽ⁿ⁾

dξ =ξH⁽ⁿ⁾−H⁽ⁿ⁺¹⁾. (2.64)

Using the recurrence relation for Hermite polynomial (2.46) and their orthonormality (2.47), we find the weights to be

w_a= n!

nH⁽ⁿ⁻¹⁾(ξ_a)2. (2.65)

This applies to integrals in one dimension, while there is not a general theory for higher dimensions. However some velocity sets can be found through ’production’ formulas.

Without going into the details of their derivation, the higher dimension velocity values, ξ_a, can be found as combinations of the values in one dimension, and their weights are products of the one-dimension weights. Here we are interested in the D2Q9 velocity set, which will be used in the numerical simulations. We are using the convention DdQq, where d is the dimension and q is the number of discrete velocity values. The D1Q3 set contains the values ξ1 = 0, ξ2 = √

3, ξ3 = −√

3 with the corresponding weights w₁ = 2/3, w₂ = w₃ = 1/6, which can be derived from formula (2.65) setting n = 3.

From this we find the D2Q9 set:

ξ₀ = (0,0) w₀ = ⁴₉

ξ_1,3 = (±√

3,0), ξ_2,4 = (0,±√

3) w1,2,3,4 = ¹₉ ξ_5,6,7,8= (±√

3,±√

3) w5,6,7,8 = ₃₆¹

.

We would like to have velocities that lead us from one site to another on a square lattice of unitary length. For this reason, rather than using these values,ξ_a, we will change to:

ca=csξ_a, (2.66)

wherecs= 1/√ 3.

We are now able to obtain the full distribution function up to Nth order from the set ofdvalues f^N(x,ca, t). This in turn means that knowing the value off^N for a discrete set of velocity values allows us to study the system with the same precision as the thermodynamic variable given by the firstN velocity moments would. Defining

fa(x, t) = wa

ω(c_a)f(x,ca, t), (2.67) we now have d functions of space and time, for which we need to find equations to describe their evolution. To do so is easy by calculating the BGK Boltzman equation at the values ca and multiplying it bywa/ω(ca).

∂fa(x, t)

∂t +c_a·∂fa(x, t)

∂x + wa

ω(c_a)g·

∂f(x,c, t)

∂c

c=ca

| {z }

f orceterm

=−1 τ











fa(x, t)− wa

ω(c_a)f⁰(x,ca, t)

| {z }

equilibrium term











. (2.68)

The expansion for the velocity-derivative, ∂_cf, in the force term is easily derived:

∂f(x,c, t)

∂c =

∞

X

n=0

1

n!a⁽ⁿ⁾(x, t)(−1)ⁿ∂ⁿ⁺¹ω(c)

∂cⁿ⁺¹

=−ω(c)

∞

X

n=0

1

n!a⁽ⁿ⁾(x, t)H⁽ⁿ⁺¹⁾(c). (2.69) Let us call the force term in (2.68),F_a. Inserting (2.69) in the expression for the force term we find its expansion to be

F_a≡ w_a ω(c_a)g·

∂f(x,c, t)

∂c

c=ca

=w_a

∞

X

n=0

1

n!a⁽ⁿ⁾H⁽ⁿ⁺¹⁾(c_a)

≈waρ[g·ca+ (g·ca)(u·ca)−g·u], (2.70) where we have used (2.49) and kept only the first few terms of the expansion.

Let us now look at the equilibrium term. We have to expand the equilibrium distribution using the same technique we used for f. First we rewrite the equilibrium distribution (2.40) in terms of the weighing function (2.43):

f⁰= ρ θω

ξ−u

√ θ

(2.71)

remembering thatv=ξ−uis the departure from the average velocity. The expansion coefficients for the equilibrium distribution are then given by an expression similar to (2.48):

a⁽ⁿ⁾₀ (x, t) = Z

f⁰(x,ξ, t)H⁽ⁿ⁾(ξ)dξ =ρ Z

ω(η)H⁽ⁿ⁾(

√

θη+u)dη (2.72) where, in the last step, we inserted (2.71) and changed integration variable toη= (ξ− u)/√

θ. Inserting the Hermite polynomials from (2.44), we get the first few equilibrium expansion coefficients:

a⁽⁰⁾₀ =ρ (2.73a)

a⁽¹⁾_0;i =ρu_i (2.73b)

a⁽²⁾_0;ij =ρu_iu_j, (2.73c)

where we have adopted the isothermal assumption and takenθ= 1.

We can now rewrite the equilibrium term like we did for the force term:

f_a⁰≡ wa

ω(ca)f⁰(x,ca, t) =wa

∞

X

n=0

1

n!a⁽ⁿ⁾H⁽ⁿ⁾ (2.74)

≈waρ

1 +ca·u+1 2

(ca·u)²−u²

(2.75) where we kept only the first few terms of the expansion.

We can now write the velocity-discretized Boltzmann equation:

∂f_a(x, t)

∂t +c_a·∂f_a(x, t)

∂x +F_a(x, t) =−1

τ f_a(x, t)−f_a⁰(x, t)

(2.76) wherea= 1, ..., d.

Let us notice that we can apply the Gaussian quadrature to the velocity moments as well and express them as finite sums:

ρ=

d

X

a=1

f_a (2.77a)

ρu=

d

X

a=1

faca (2.77b)

P +ρuu=

d

X

a=1

f_ac_ac_a (2.77c)

2.2.4 Discretization in time and space

We have discretized the Boltzmann equation in velocity space and we now have d equations (2.76) for functions of xand twhere the velocitiesca are constants.

The next step is to discretize in space and time. This can be done with different techniques, but we are going to follow the method of characteristics. We transform the partial differential equations (2.76) in an ordinary differential equation by parametrizing the variablesxand t:

x=x(ζ), t=t(ζ). (2.78) By doing this we can express (2.76) in terms of the total derivative of f_a(ζ):

df_a

dζ = ∂f_a

∂t dt dζ +∂f_a

∂x ·dx

dζ =−f_a−f_a⁰

τ −F_a. (2.79)

Comparing (2.79) to (2.76 we get dt

dζ = 1 and dx

dζ =ca (2.80)

from which we see that the trajectory we are solving along isx=c_at+x₀, wherex₀ is an arbitrary constant. The starting point is set to (x(0), t(0)) ≡(x0, t0) and we want to integrate over a time step ∆t to the point (x(∆t), t(∆t) = (x0 +ca∆t, t0 + ∆t).

Since the starting point (x₀, t₀) is arbitrary, solving (2.79) along this trajectory gives an expression forfa(x+ca∆t, t+ ∆t) as a function of fa(x, t).

Equation (2.79) can be solved by variation of parameters. The solution to the corre- sponding homogeneous equation is:

f_a^C(ζ) =fa(ζ0)e⁻

ζ−ζ0

τ . (2.81)

The particular solution to the original equation (2.79) is then set to:

f_a^P(ζ) =C(ζ)e^{−ζ−ζτ0}. (2.82) Inserting this is (2.79) we find:

C(ζ) = Z ζ

ζ0

f_a⁰(ζ⁰)

τ −F_a(ζ⁰)

e^{ζ0−ζτ 0}dζ⁰. (2.83) Summing the general homogeneous solution and the particular inhomogeneous solution we get the general solution for (2.79):

f_a(ζ) =f_a^C(ζ) +f_a^P(ζ)

=e⁻

ζ−ζ0 τ

f_a(ζ₀) + Z _ζ

ζ0

f_a⁰(ζ⁰)

τ −F_a(ζ⁰)

e^{ζ0−ζτ 0}dζ⁰

. (2.84)

As discussed previously, we want to integrate over a time step. We therefore set ζ = ζ₀+ ∆tand transform back toxand t:

fa(x+ca∆t, t+ ∆t)

=e^−∆tτ [ fa(x, t) + Z t+∆t

t

f_a⁰(x+ca(t−t⁰), t⁰) τ

−F_a(x+c_a(t−t⁰), t⁰) e^t0−tτ

dt⁰ (2.85)

where we have substituted (x0, t0) with (x, t) since the starting point is arbitrary.

The integral in (2.85) can be discretized in many ways. The simplest is a rectangular discretization, where the integral is approximated by only one point.

f_a(x+c_a∆t, t+ ∆t) =

1−∆t τ

f_a(x, t) +∆t

τ f_a⁰(x, t)−∆tF_a(x, t) (2.86) where the exponential has been expanded and we have kept terms up to first order in

∆t. Using the trapezoidal approximation for the integral, the equation is fa(x+ca∆t, t+ ∆t)−fa(x, t) =

−∆t 2τ

fa(x+ca∆t, t+ ∆t)−f_a⁰(x+ca∆t, t+ ∆t)−τ Fa(x+ca∆t, t+ ∆t) +fa(x, t)−f_a⁰(x, t)−τ Fa(x, t) +O(∆t³)

. (2.87)

This does not look promising;f_a⁰(x+c_a∆t, t+ ∆t) depends onf_a(x+c_a∆t, t+ ∆t) and therefore (2.87) is time-implicit. Fortunately this problem can be solved by introducing a new set of populations

f¯a=fa+∆t 2τ

fa−f_a⁰−τ Fa

(2.88)

and inserting them into (2.87), which leads to f¯_a(x+c_a∆t, t+ ∆t)−f¯_a(x, t) =− ∆t

τ + ∆t/2

f¯_a(x, t)−f_a⁰(x, t)−τ F_a(x, t

. (2.89) The new populations have moments closely related to the original ones:

X

a

f¯_a=X

a

f_a+∆t 2

X

a

f_a−f_a⁰

τ −∆t

2 X

a

F_a=ρ (2.90)

X

a

f¯_ac_a=X

a

f_ac_a+∆t 2

X

a

f_a−f_a⁰

τ c_a−∆t 2

X

a

F_ac_a=ρu− ∆t

2 F, (2.91) whereP

aFa= 0 and P

aFaca=F, both of which can be seen to be true from (2.70), were used. The sums involving the collision operator, (f_a−f_a⁰)/τ, are zero since by construction the operator preserves collision invariants (see 2.36). Velocity and density can now be found as

ρ=X

a

f¯_a (2.92)

u=X

a

f¯_ac_a

ρ +∆tF

2 (2.93)

where the new populations satisfy f¯_a(x+c_a∆t, t+∆t)−f¯_a(x, t) =−∆t

¯

τ [ ¯f_a(x, t)−f_a⁰(x, t)]+

1−∆t

2¯τ

F_a(x, t)∆t. (2.94)

A new relaxation time ¯τ =τ+ ∆t/2 has been introduced. It can be shown through the Chapman-Enskog analysis [23], that the viscosity can be expressed asν =c²_s(¯τ−∆t/2).

For convenience,τ and fa are written in place of ¯τ and ¯fa henceforth.

The Chapman-Enskog analysis is an expansion of the one-particle distribution function around the local equilibrium. The parameter used for the expansion is the Knudsen number

Kn= l_{mf p} lM

, (2.95)

where l_{mf p} is the molecular mean free path and l_M is the characteristic length for the gradients of macroscopic quantities. Through the Chapman-Enskog analysis, the macroscopic equations of fluid dynamics can be recovered from the Boltzmann equation [24].

2.3 Lattice Boltzmann method

The Lattice Boltzmann method solves the discretized Boltzmann equation by evolving thedone-particle distribution functions,fa, for thedvalues of velocity. The discretized one-particle distribution functions,f_a, are often called populations in the Lattice Boltz- mann Method, and that is the term we will use from now on. The populations are stored in an array of rank 3 and size nx, ny, d, wherenx×ny is the size of the lattice anddis the number of values in the velocity set. This is for a system in two dimensions, and therefore a two-dimensional lattice. In general the populations array will have rank D+ 1, where D is the dimension of the system. The most common discrete velocity set used for 2D simulations is the D2Q9, presented in section 2.3.1. The equilibrium populations are held in a similarly shaped array, while density and the components of velocity are stored in a array each of rank 2 and size of the lattice, (nx, ny). At each time step the values of the populations are evolved according to the equations (2.86). This is the core time-step algorithm, detailed in section 2.3.2. Boundary conditions are a very important part of LBM simulations, and many methods to implement them, exist [23].

Here two variations of the bounce-back method are described (§ 2.3.3). The discussion up to this point is about a single LBM sets of populations to solve the Navier-Stokes equation. Section 2.3.4 shows how a second set of populations can be introduced to simulate the advection and diffusion of the solute (or temperature).

The outline of the Lattice Boltzmann Method in this section follows the description by T. Kr¨ugeret al. [23].

2.3.1 D2Q9

As was mentioned in section 2.2.3, the D2Q9 velocity set is derived from the D1Q3 set through a production formula. It is preferable to have velocities that, starting from one lattice site, after one time step end at another lattice site, rather than somewhere in the middle. For this reason, the rescaling 2.66 is adopted. The velocities in the D2Q9

are

c0= (0,0) w0 = ⁴₉

c_1,3= (±1,0), c_2,4= (0,±1) w_1,2,3,4= ¹₉ c_5,6,7,8 = (±1,±1) w_5,6,7,8= ₃₆¹, where thewa are the weights. Figure 2.3.1 shows the D2Q9 set velocities.

c₁ c₃

c2

c4

c5

c6

c7 c8

c₀

Figure 2.1: The D2Q9 velocity set. The velocity values are scaled so that in a time step the populations are carried from one lattice node to another.

2.3.2 Core time step algorithm

The algorithm for one time step can be divided into substeps. These need to be in a specific order, but one can start from any of the substeps.

In the moment update the velocity moments of the distribution function, usually just density and velocity, are calculated through:

ρ=

d

X

a=1

f_a, ρu_i =

d

X

a=1

c_a,if_a+∆t

2 F_i (2.96)

where i = x, y. In the case of the D2Q9 velocity set, d= 9 and the formulas can be written explicitly.

ρ=f₀+f₁+f₂+f₃+f₄+f₅+f₆+f₇+f₈ (2.97) ux= [(f1+f5+f8)−(f3+f6+f7) + ∆tFx/2]/ρ (2.98) uy = [(f2+f5+f6)−(f4+f7+f8) + ∆tFy/2]/ρ (2.99) From these we can find the new equilibrium values for the distribution function corre- sponding to the new macroscopic state. This is done through the expression we found for the equilibrium term, f_a⁰, when expanding in terms of Hermite polynomials:

f_a⁰ =w_aρ

1 +c_a·u+ 1 2

(c_a·u)²−u²

, (2.100)

where we have kept only terms up to second order in the discrete velocity. In particular for D2Q9 we find the following formulas for the nine equilibrium populations:

f₀⁰ = ^2ρ₉ 2−3u² ,

f₁⁰ = ₁₈^ρ 2 + 6ux+ 9u²_x−3u²

, f₅⁰= ₃₆^ρ

1 + 3(ux+uy) + 9uxuy + 3u² f₂⁰ = ₁₈^ρ 2 + 6u_y+ 9u²_y−3u²

, f₆⁰= ₃₆^ρ

1−3(u_x−u_y)−9u_xu_y + 3u² f₃⁰ = ₁₈^ρ 2−6ux+ 9u²_x−3u²

, f₇⁰= ₃₆^ρ

1−3(ux+uy) + 9uxuy + 3u² f₄⁰ = ₁₈^ρ 2−6u_y+ 9u²_y−3u²

, f₈⁰= ₃₆^ρ

1 + 3(u_x−u_y)−9u_xu_y + 3u² . (2.101) We can now calculate the new populations through the discretized Lattice Boltzmann equation:

f_a(x+c_a∆t, t+∆t) =

1− ∆t τ

f_a(x, t)+∆t

τ f_a⁰(x, t)+

1−∆t

2¯τ

F_a(x, t)∆t. (2.102) This equation is solved in two steps: first the right hand side (the collision term) is determined, then the population values are moved to new lattice sites following their velocity,ca. This last step, called streaming, can be carried on in several ways, but the simplest is probably to copy the post-collision populations on a temporary array and then copy them back onto the population array in their new position. To help visualize the streaming process, let us look at an example. Let us assume that the value, f₆^∗, of the post-collision population in the direction of c6 at timet at the lattice site (i, j), whereiindicates the column and j the row, isf:

f₆^∗(i, j, t) =f. (2.103)

At time t+δtthe value of the populationf6 at the lattice site (i−1, j+ 1) will be set tof:

f6(i−1, j+ 1, t+ ∆t) =f₆^∗(i, j, t) =f. (2.104) The use of a temporary array is just to prevent writing over values that have not yet been copied.

2.3.3 Boundary conditions

The treatment of boundary conditions is a very important aspect of Lattice Boltzmann simulations. Many different methods have been explored. Among the simplest ones is the bounce-back scheme. This method sets the boundary to be approximately half a lattice constant outside of the last fluid lattice site. Let us consider the boundary condition to be a wall, on which we have a no-slip condition. The momentum of a fluid particle hitting the wall is reversed: ρu → −ρu. In the full-way bounce-back (FBB)