A Pract Practica ical l Intr Introduct oduction ion to the the Lattice Lattice Boltzmann Boltzmann Method Alexander J. Wagner Depar Dep artment tment of Physics North North Dakot Dakota a Sta State University
[email protected]
2
Contents 1 The Boltzmann equation
5
2 Derivation of hydrodynamics
2.1 2.2
2.3 2.4 2.5
7
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . The Boltzmann equation . . . . . . . . . . . . . . . . . . 2.2.1 The BGK approximation . . . . . . . . . . . . . 2.2. 2.2.22 Mo Mome men nts of the the equi equili libr briu ium m dist distri ribu buti tion on func functi tion on Mass conservation . . . . . . . . . . . . . . . . . . . . . Momentum conservation . . . . . . . . . . . . . . . . . . Energy conservation . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
3 Lattice Boltzmann
3.1 3.2 3.3 3.4 3.5 3.6
13
The lattice Boltzmann equation . Taylor expansion . . . . . . . . . One dimensional implementation Isothermal Lattice Boltzmann . . Non-ideal fluids . . . . . . . . . . Boundaries . . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
4 Applications of Lattice Boltzmann
4.1 4.2 4.3
5.3
Introduction . . . . . . . . No conserved parameter . 5.2.1 Magnetic Systems Conserved zeroth moment 5.3.1 Diffusion . . . . . 5.3.2 Electrostatics . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
13 14 14 15 17 18 21
Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Liquid-gas systams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Two-fluid systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Non-traditional Lattice Boltzmann methods
5.1 5.2
7 7 7 8 8 8 9
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
21 21 21 23
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
23 23 23 23 24 24
6 Possible Pro jects
25
A Einstein convention for vectors and tensors
29
B Functional Derivatives
31
C Lagrangian multipliers
33
D Evaluating Gaussian integrals
35
E Galilean transformations
37
F D1Q5 Lattice Boltzmann code
39
3
4
CONTENTS
Preface The lattice Boltzmann method is increasingly attracting researchers in many areas from turbulence to multi-phase multi-phase flow in porous media. Several Several textbooks have have been written written to address address the need of students to learn about this relatively new method. The aim of this introduction is to provide a succinct description of the field and to provide students with sample codes so that they can immediatly apply their knowledge to practical applicati plications ons.. Source Source codes, which which are provid provided ed under under the GNU copyri copyrigh ght, t, can be downl download oaded ed at http://www.physics.ndsu.nodak.edu/wagner/LB.html. This manuscript is still in a rough state and is continuously beeing improoved. Your comments and corrections are always welcome.
A.J. Wagner, Fargo, March 2005.
Chapter 1
The Boltzmann equation We have already seen 1 that the dynamics of the Boltzmann equation always mimimizes the H Functional given by H (t) =
dx dv f ( f (x, v, t)log(f )log(f (x, v, t)) .
(1.1)
So we can conclude that the equilibrium distribution function f 0 in a volume V for a given density n, mean momentum momentum nu and energy nǫ = 1/ 1 /2 nu2 + 3/2 nθ will minimize the H -functional. -functional. We can use Lagrangian Lagrangian multipliers multipliers to minimize minimize this functional. functional. With Lagrangian Lagrangian multiplie multipliers rs the H -functional -functional reads
H (t) =
dx dv f (x, v, t)log f ( f (x, v, t)
−
−λ1(nV dx dv f ( f (x, v, t)) −λ2α(nuαV − dx dv f ( f (x, v, t)vα ) −λ3(nǫV
−
dx dv f ( f (x, v, t)
v2
(1.2)
2
Now we calculate the functional derivative of the H -functional -functional with respect to the distribution function. function. Since the equilibrium equilibrium distribut distribution ion minimizes minimizes the H -functiona -functionall this derivativ derivativee has to 0 vanish when the distribution is the equilibrium distribution f = f . δH ( δH (t) 0= δf !
= 1 + log(f log(f 0 ) + λ1 + λ2α vα + λ3
f = f =f 0
v2
(1.3)
2
We can solve this for the equilibrium distribution: f 0 = exp
−1 −
λ2 λ1 + 2 2λ3
−
λ3 2
λ2 +v λ3
2
(1.4)
Introducin Introducingg a new set of Lagrangian Lagrangian multipliers multipliers a, bα and c we can also write this as 0
f = a exp
(b
− v) 2 c
(1.5)
Note that this expression does not depend on spatial variables. So we know that the solution will be homogeneous. We now find the Lagrange multipliers by invoking the conservation laws nV = 1
0
dx dv f
nuα V =
0
dx dv f vα
Kerson Kerson Huang, Huang, Statistical Statistical Mechanics, Mechanics, 2nd ed., chapter 4.
5
nu2 3 V + nθV = 2 2
dx dv f 0
v2
2
(1.6)
6
CHAPTER CHAPTER 1. THE BOLTZMANN BOLTZMANN EQUA EQUATION
Solving Solving these equations equations for a, bα and c we find that the equilibrium distribution is given by n f = exp (2πθ (2πθ))3/2 0
−
(v
− u)2 2θ
(1.7)
which is known as the Maxwell-Boltzmann distribution.
Problems 1.1: Which function φ(x) minimizes the functional Ψ[φ Ψ[ φ] =
b a
φ2 (x) + φ4 (x)dx? dx?
− − b a
φ2 (x) + φ4 (x)dx if the integral b a) is a consta constant nt?? Do you obtain obtain differen differentt kinds kinds of soluti solutions ons for certai certain n a φ(x)dx = c(b ranges of the constant c?
1.2: Which function φ(x) minimizes the functional Ψ[φ Ψ[ φ] =
−
1.3: Derive (1.7) [i.e. determine a, b, and c for equation (1.5)].
Chapter 2
Derivation of the hydrodynamic equations from the Boltzmann equation 2.1 2.1
Intr Introdu oduct ctio ion n
In this lecture we will examine the hydrodynamic limit of the Boltzmann equation and derive the transport equations for the macroscopic quantities from first principles. We will see that the macroscopic scopic equations equations of motion motion are simply the conservation conservation equations equations for continuous continuous fields. Because Because of the general concepts involved the transport equations we derive are not only applicable for dilute gases, which we require for the Boltzmann equation to apply, but also for much denser fluids. This is the reason that recently a numerical method called “lattice Boltzmann” has been develope veloped d for the simula simulatio tion n of fluids. fluids. We will will cove coverr the lattice lattice Boltzmann Boltzmann approac approach h in the next lecture.
2.2 2.2
The The Boltz Boltzma mann nn equa equatio tion n
The Boltzmann equation we derived in the last lecture is given by ∂ t f + v∂ x f + F ∂ v f =
dv1′ dv2′ dv2 (f 1′ f 2′
− f 1f 2)P 1212
12
→
′
′
(2.1)
Solving this equation analytically is very challenging and can only be done for special cases.
2.2.1 2.2.1
The The BGK BGK appro approxim ximati ation on
There Bhatnagar, Gross and Krook noticed, however, that the main effect of the collision term is to bring the velocity velocity distributio distribution n function function closer to the equilibrium equilibrium distribution. distribution. The equilibrium equilibrium distribution is given by (1.7) n f 0 (v) = e−(v −u) /2θ (2.2) (2πθ (2πθ))3/2 2
where n is the number density of particle, u is the mean velocity and θ = kT is the temperature. These macroscopic quantities are given by moments of the distribution function:
f = n
(2.3)
f vα
= nuα
(2.4)
f v2
= nu2 + 3nθ 3nθ
(2.5)
7
8
CHAPTER CHAPTER 2. DERIVA DERIVATION OF HYDRODYNAMICS HYDRODYNAMICS
The simplest way of approximating the the collision term is by using a single relaxation time approximation 1 0 dv1′ dv2′ dv2 (f 1′ f 2′ f 1 f 2 )P 12 (f f ) f ) (2.6) 12→1 2 τ One important feature of this approximation is that mass, momentum and energy are are still exactly conserved by the collision term. With this approximation the Boltzmann equation reads
−
′
∂ t f + v∂ x f + F ∂ v f =
′
≈
1 0 (f τ
−
− f ) f )
(2.7)
We can use this expression to approximate the probability density function by the equilibrium distribution and its derivatives: f = f 0
2.2.2
− τ ( τ (∂ t f + v∂ x f + F ∂ v f ) f )
(2.8)
Momen Moments ts of the equili equilibriu brium m distribu distribution tion functi function on
Because we can express the distribution function in terms of the equilibrium distribution function we will require require several velocity velocity mom moment entss of the equilibrium equilibrium distribution distribution function (1.7). They are
f 0
= n,
(2.9)
− uα) f 0 (vα − uα )(v )(vβ − uβ )
= 0,
(2.10)
= nθδαβ ,
(2.11)
− uα)(v )(vβ − uβ )(v )(vγ − uγ )
= 0,
(2.12)
= 5nθ2 δαβ .
(2.13)
f 0 (vα
f 0 (vα
f 0 (vα
− uα)(v )(vβ − uβ )(v − u)2
This is easily derived using Gaussian integrals (see appendix D).
2.3
Mass Mass conse conserv rvat atio ion n
Integrating over the Boltzmann equation we obtain for the mass conservation equation ∂ t
1 dvfvα + F dv∂ v f = τ ∂ t n + ∂ α (nuα ) = 0
dvf + ∂ α
⇔
dv( dv (f 0
− f ) (2.14)
which is the continuity equation.
2.4
Moment Momentum um conser conserv vation ation
Multiplying the Boltzmann equation with vα and integrating we obtain as a momentum conservation conservation equation ∂ t
dvfvα + ∂ β
dvfvα vβ + F β
∂ t (nuα ) + ∂ β
dv∂ vβ f vα
dvfvα vβ
− nF α
=
1 τ
dv( dv(f 0
− f ) f )v
= 0
We now need to use equation (2.8) to approximate
dvfvα vβ =
0
dvf vα vβ τ ( τ (∂ t
−
0
dvf vα vβ +∂ γ
)+ O(∂ 2 ) (2.15) dvf 0 vα vβ vγ +nF α uβ + nuα F β )+O
9
2.5. ENERGY CONSERV CONSERVATION
To first order in derivatives our conservation equation now reads ∂ t (nuα ) + ∂ β (nuα uβ ) =
−∂ α(nθ) nθ) + nF α .
(2.16)
This equation equation is known known as the Euler equation. equation. Using the continuit continuity y equation (2.14) (2.14) we can also write it as 1 ∂ t uα + uβ ∂ β uα = ∂ α (nθ) nθ) + F α . (2.17) n To calculate the equations of motion to second order in the derivatives we need to evaluate the higher order terms in eqn. (2.15).
−
∂ t (
f 0 vα vβ )
= ∂ t (nuα uβ + nθδαβ ) = ∂ t (nuα )uβ + nuα ∂ t uβ + ∂ t nθδαβ + n∂ t θδαβ = ∂ γ (nuα uγ )uβ ∂ α (nθ) nθ)uβ nF α uβ nuα uγ ∂ γ uβ uα ∂ β (nθ) nθ) nuα F β 2 ∂ γ (nuγ )θδ αβ nuγ ∂ γ (θδαβ ) ∂ γ uγ θ 3 = ∂ γ (nuα uβ uγ ) ∂ β (nθ) nθ)uα ∂ α (nθ) nθ)uβ n(F αuβ + uα F β ) 2 ∂ γ (nθuγ )δαβ nθ∂ γ uγ 3
− − − − −
− − − − −
− −
−
−
−
(2.18)
where we used equations equations (2.16),(2.17 (2.16),(2.17)) and (2.24). We need to com combine bine this with ∂ γ
f 0 vα vβ vγ
= ∂ β (nθuα ) + ∂ α (nθuβ ) + ∂ γ (nθuγ )δαβ + ∂ γ (nuα uβ uγ )
(2.19)
where we used equation equation (2.12). Combining Combining the terms of the two previous equations equations we have have ∂ t
0
dvf vα vβ + ∂ γ
= nθ( nθ(∂ α uβ + ∂ β uα )
dvf 0 vα vβ vγ
− 23 nθ∂ γ uγ
(2.20)
Now we can use this expression to obtain the first order momentum conservation from equation (2.15): 2 n∂ t uα + nuβ ∂ β uα = ∂ α (nθ) nθ) + nF α + ∂ β [η (∂ β uα + ∂ α uβ ∂ γ uγ δαβ )] (2.21) 3 where η = nθτ is the viscosity. This equation is known as the Navier-Stokes Equation.
−
2.5 2.5
−
Ener Energy gy cons conser erv vation ation
Multiplying the Boltzmann equation with (v ( v u)2 and integrating we obtain for the energy equation
−
dv∂ t f ( f (v
− u)
2
+
⇔ ⇔
− τ
∂ t (vα
dv∂ α f vα (v
∂ t +∂ α ∂ t
2
2
− u) + F α dv∂ v f ( f (v − u) dvf ( dvf (v − u)2 + dvf 2( dvf 2(vvα − uα )∂ t uα
dvfvα (v
dvf ( dvf (v
2
− u)
− u)2 + ∂ α
− uα)(v )(vβ − uβ )
+
dvfvα 2(v 2(vβ
dvfvα (v
∂ γ f 0 vγ (vα
α
− uβ )∂ αuβ
=
1 τ
= 0
− u)2 + 2nθ∂ 2nθ∂ α uα
− uα)(v )(vβ − uβ )
∂ α uβ
= 0
dv( dv(f 0
− f )( f )(vv − u)2
10
CHAPTER CHAPTER 2. DERIVA DERIVATION OF HYDRODYNAMICS HYDRODYNAMICS
3∂ t (nθ) nθ) + ∂ α
⇔
−τ ∂ αuβ
dvfvα (v
− u)2 + 2nθ∂ 2nθ∂ α uα
∂ α uβ + ∂ β uα
−
2 ∂ γ uγ δαβ 3
= 0
(2.22)
We now need to approximate the remaining integral using eqn. (2.8).
dvfvα (v
=
− u)2
0
2
dvf vα (v
+
− u) − τ [τ [ dv∂ v f vα (v − u)2 ]
= 3nθuα
− τ [τ [
0
dv∂ t f vα (v
dv∂ t f vα (v
− u)
2
+
−5nθF α] + O(∂ 2).
2
− u)
+
dv∂ β f vα vβ (v
dv∂ β f 0 vα vβ (v
− u)2
− u)2 (2.23)
So to zeroth order we can write the energy conservation equation as ∂ t θ + uα ∂ α θ =
− 23 ∂ αuαθ + O(∂ 2).
(2.24)
where we have again used the continuity equation (2.14). To obtain the first order equation we now need to evaluate the two integrals in equation (2.23).
dv∂ t f 0 vα (v
= ∂ t
0
f vα (v
− u)2 2
− u)
+
2(vγ f 0 vα 2(v
− uγ )∂ tuγ
= ∂ t (3nθu (3nθuα ) + 2nθ∂ 2 nθ∂ t uα
2 − ∂ α(nθ) nθ) + nF α ] + 3nu 3 nuα (−uβ ∂ β θ − ∂ β uβ θ) 3 1 +2nθ +2nθ[[−uβ ∂ β uα − ∂ α (nθ) nθ) + F α ] n ∂ β (−3θnuα uβ ) − 2nθ∂ β (uα uβ ) − 5θ∂ α (nθ) nθ) + 5nF 5 nF α
−
= 3θ[ ∂ β (nuα uβ )
=
(2.25)
For the second integral we get
= ∂ β = ∂ α
dv∂ β f 0 vα vβ (v
0
f vα vβ (v dvf 0 (vα
− u)2
− u)
2
+
f 0 vα vβ 2(v 2(vγ
− uγ )∂ β uγ
− uα)(v )(vβ − uβ )(v )(v − u)2 − ∂ β (3nθu (3nθuα uβ )
+2nθ +2nθ((uα ∂ β uβ + uβ ∂ β uα ) = ∂ α (5nθ (5nθ2 ) + ∂ β (3nθu (3nθuα uβ ) + 2nθ∂ 2 nθ∂ β (uα uβ )
(2.26)
Combining both integrals from equation (2.23) and the forcing term we obtain ∂ α (5nθ (5nθ2 )
− 5θ∂ α(nθ) nθ)
= 5nθ∂ α θ
(2.27)
So that we get the heat conduction equation ∂ t θ + uα ∂ α θ =
−
2 1 5nθ ∂ α uα θ + ∂ α ( ∂ α θ ) + τ ∂ α uβ ∂ α uβ + ∂ β uα 3 n 3
−
2 ∂ γ uγ δαβ 3
(2.28)
11
2.5. ENERGY CONSERV CONSERVATION
Problems 2.1: Derive (2.9–2.13). 2.2: Show that
f 0
f 0 vα
f 0 vα vβ
f 0 vα vβ vγ
= n, = nuα , = nuα uβ + nθδαβ , = nθ( nθ(uα δβγ + uβ δαγ + uγ δαβ ) + nuα uβ uγ
for the Maxwell-Boltzmann distribution of equation (1.7). 2.3: Derive the hydrodynamic equations for a BGK Boltzmann equation (2.7) with an equilibrium
distribution given by
f 0 (v) = n exp( π v2 )
−
(2.29)
where n = f dv. Hints: Consider Consider what quantities quantities (if any) are conserved conserved for this evolutio evolution n equation. What quantities are therefore appropriate fundamental variables? When deriving the hydrodynamic equation(s) for this model make sure that you replace all non-conserved quantities with appropriate approximations (in the same way that we used (2.8) to approximate f vvdv in (2.15).
12
CHAPTER CHAPTER 2. DERIVA DERIVATION OF HYDRODYNAMICS HYDRODYNAMICS
Chapter 3
Lattice Boltzmann In the last chapter we have seen that the Boltzmann equation describes a dynamics that includes Newtonian Newtonian hydrodynamic hydrodynamicss in the long waveleng wavelength th limit. This may at first sight seem surprising surprising since the Boltzmann Boltzmann equation equation was derived only in the limit for rare gases and fluids are dense. dense. And if we look closely we see that the transport coefficients, the viscosity ν , the heat conductivity κ and the speeds of sound are all closely related related in a way that is not true for all fluids. But the general general structure of the conservation laws is so general that they apply for most continuous media. This opens the way for an alternativ alternativee way to simulate simulate fluids. Instead Instead of trying trying to discretize discretize the continuity, Navier-Stokes and heat equations directly a simple discretization of the Boltzmann equation surfices. Let us briefly review what we required to derive the conservation equations. It is important to notice notice that we reduced all our calculations calculations to calculation calculationss over over the equilibrium equilibrium distribution. distribution. And we only required some basic moments of the equilibrium distribution which are given by eqns. (2.9– 2.13). 2.1 3). So we can conclu conclude de that any distribution function with these moments would lead to the same macroscopic equations.
3.1 3.1
The The latti lattice ce Bolt Boltzm zman ann n equa equati tion on
Let us now write down a simple discretization of the Boltzmann equation with BGK approximation(2.7) for the collision term f ( f (x + vi , vi , t + 1)
1 − f ( f (x, vi , t) + F ( F (vi ) = [f 0 (n,u,θ) n,u,θ) − f (x, vi , t)] τ
(3.1)
Here we have discretized velocity space to a finite number of velocity vectors vi , space to a lattice where we require that x + vi is again a lattice position and time only takes on integer values. Since the velocity vectors are fixed now we usually denote f ( f (x, vi , t) f i (x, t) and F ( F (vi ) F i . The force terms F i are defined as a generalization of the force of (2.7), i.e. F i F α ∂ vα f . f . In particu particular lar we demand that the moments match up:
≡
F i =
i
F i viα =
i
F i viα viβ =
i
F i viα viβ viγ =
i
dvF α ∂ vα f
↔
≡
= 0, 0,
(3.2)
dvF β ∂ vβ f vα
=
−nF α,
(3.3)
dvF γ ∂ vγ f vα vβ
=
−n(F αuβ + uαF β ),
(3.4)
dvF δ ∂ vδ f vα vβ vγ
=
−n[F α(θδβγ + uβ uγ ) + F β (θδαγ + uαuγ ) +F γ (θδαβ + uα uβ )]. )].
(3.5)
We now need to show that the hydrodynamic limit of this this discretized version of the Boltzmann are still the well known equations for fluid flow. 13
14
CHAPTER CHAPTER 3. LATTICE LATTICE BOLTZMANN BOLTZMANN
Problems 3.1.1: Show that the moments given in (3.2) to (3.5) are indeed what I claim they are.
3.2
Taylor ylor expa expans nsion ion
To determine what the macroscopic equations are that the lattice Boltzmann equation simulates we perform a Taylor expansion of equation (3.1). We obtain to second order ∂ t f i + viα ∂ α f i +
1 1 [∂ t (∂ t f i + viα ∂ α f i ) + ∂ β (∂ t f i viβ + viβ viα ∂ α f i )] + F i + O(∂ 3 ) = (f i0 2 τ
− f i) (3.6)
We notice that this is not quite the Boltzmann equation (2.8) that we set out to simulate because there are a large number of additional terms with the second derivative. These terms are discretization error because of the simple discretizati discretization on scheme we used. Howeve However, r, we will not be deterred deterred by that for the moment and we can still write f i = f i0
− τ ( τ (∂ t f + viα ∂ α f i + F i )
(3.7)
to express the f i in terms of the equilibrium distribution f i0 . Now expressing expressing everything everything (exccept (exccept the collision term) in terms of the equilibrium distribution we get ∂ t f i0
− τ ∂ t(∂ tf i + viα∂ αf i + F i) + viα∂ αf i0 + ∂ α(∂ tf i0viα + viαviβ ∂ β f i0 + F i) +
1 ∂ t (∂ t f i + viα ∂ α f i ) + ∂ β (∂ t f i0 viβ + viβ viα ∂ α f i0 ) + O(∂ 3 ) = 2
1 0 (f τ i
− f (3.8) i)
We note now that, through a lucky coincidence, the discretization errors are of exactly the same form as the higher order terms for the expression of the distribution function in terms of the equilibrium distribution function1 . We can now write ∂ t f i0
−
+ viα ∂ α f i + τ
1 2
∂ t (∂ t f i + viα ∂ α f i ) + ∂ β (∂ t f i0 viβ + viβ viα ∂ α f i0 ) + O(∂ 3 ) =
1 0 (f τ i
− f i)
(3.9) which is exactly the same equation we would have obtained for the Boltzmann equation, exccept that the relaxation time is renormalized to be τ 1/2. So if we choose an equilibrium distribution with the appropriate moments (2.9–2.13) we will automatically simulate the hydrodynamic equations to the same order that we derived the hydrodynamic limit.
−
3.3
One dimensi dimensional onal impleme implement ntatio ation n
So in order to implement this we only need to define an equilibrium distribution which fulfills the equivalent definition of (2.9–2.13 which are:
f i0
= n,
(3.10)
− uα)
= 0,
(3.11)
= nθδαβ ,
(3.12)
= 0,
(3.13)
= nθ2 δαβ
(3.14)
i
f i0 (viα
i
i
f i0 (viα
i
1
You also need to notice that
F i
− uα)(v )(viβ − uβ )
− uα)(v )(viβ − uβ )(v )(viγ − uγ )
f i0 (viα
i
f i0 (viα
− uα)(v )(viβ − uβ )(v − u)2 i
= O(∂ ), ), but that is clear from (3.6).
15
3.4. ISOTHERMAL ISOTHERMAL LATTICE LATTICE BOL BOLTZMANN
where the difference between (3.14) and (2.13) is due to the fact that we are considering a onedimensional model instead of a three dimensional one (see problem 3.1). Since these are 5 equations (in one dimension) we can expect that we will require at least, and probably exactly, a set of 5 velocities vi and corresponding equililibrium distribution terms f i0 . If we choose choose the symmet symmetric ric velocity set vi = 2c, c, 0, c, 2c (3.15)
{ } {− −
}
we obtain for the equilibrium equilibrium distributio distribution n f 00
=
f 10
=
0 f − 1
=
f 20
=
0 f − 2
=
n 4 c4 + 3 θ 2 + 6 θ u2 + u4 5 c2 θ + u2 4 c4 n 3 θ 2 + 4 c3 u 6 θ u2 u4 + 4 c2 θ + u2 c 3 θ u + u3 6 c4 n 3 θ 2 4 c3 u 6 θ u2 u4 + 4 c2 θ + u2 + c 3 θ u + u3 6 c4 n 3 θ2 2 c3 u + 6 θ u2 + u4 c2 θ + u2 + 2 c 3 θ u + u3 24 c4 2 3 2 4 n 3θ +2c u+6θu +u c2 θ + u2 2 c 3 θ u + u3 24 c4
− −
−
−
−
−
−
−
−
− −
− −
(3.16) (3.17) (3.18) (3.19) (3.20)
This just leaves us with the actual implementation of the lattice Boltzmann method defined in (3.1). One possible implementation, which also employs my GUI, is given in appendix F.
Problems 3.3.1: Show that for a one-dimensional system you obtain (3.14) instead of (2.13). 3.3.2: How would you write a lattice Boltzman code (say one dimensional for simplicity) for an
isothermal isothermal system? system? In such a system system the energy energy is not conserved conserved but instead we force the equilibrium equilibrium temperature temperature to be b e a constant constant θ0 . a) Which Which macroscopic macroscopic equations equations will you need to simulate simulate?? b) What is the form of the Navier-Stok Navier-Stokes es equation? equation? Hint: you will want to consider the derivation of (2.21). How does this derivation change for an isothermal system? ( difficult ) c) What are the required required moments moments of the equilibrium equilibrium distribution? distribution? d) Hence, Hence, what is the minimum number number of velocities velocities you are likely to require? require? e) Calculate Calculate the equilibrium distributio distribution n for this model. f ) Implement the the model in C. (See Appendix F for an example implementation of the D1Q5 model.) Please note that this is an enourmous simplification and the math required will be much simplified also.
3.4 3.4
Isoth Isother erma mall Latt Lattic ice e Boltz Boltzma mann nn
Most lattice Boltzmann simulations are used to only simulate the continuity and Navier-Stokes equations. equations. The temperature temperature is assumed to be constant constant and the equilibrium equilibrium distribution distribution will no longer longer conserve conserve energy; instead instead it serves serves as a thermostat. thermostat. This removes removes the requireme requirement nt for the equilibrium equations to fulfill equation (3.14). This moment was only needed to calculate the heat equations. equations. For simplicity simplicity let us now consider a one-dimensio one-dimensional nal LB model. For the full thermothermohydrodynamic model we needed 5 velocities. Now that we have dropped one constraint, you would
16
CHAPTER CHAPTER 3. LATTICE LATTICE BOLTZMANN BOLTZMANN
expect that we need a 4 velocity velocity model to fulfill the remaining remaining 4 constrain constraints. ts. But if there was a way way to reduce the number of required velocities further we could save some memory and computation time. If you now consider that you will not be interested in the absolute value of the temperature, you can use the determination of the temperature as an additional degree of freedom. You use the 0 0 4 equation to determine f − 1 , f 0 f 1 , and θ. This will certainly work but we need to remember that we also want θ to be a constant independent of n and u. We will now see that this nearly works. Using the D1Q3 velocity set vi = 1, 0, 1 it is easy to see that (3.10) to (3.12) require
{−
0 f − 1
f 00 f 10
}
1 n( u + θ + u2 ) 2 = n(1 θ u2 ) 1 = n(u + θ + u2 ) 2
− − −
=
(3.21) (3.22) (3.23) (3.24)
Using these solutions for the f i0 we can calculate θ from (3.13): θ=
1 3
2
− u3 .
(3.25)
We know that the velocity has to be much smaller than the lattice velocity c = 1 and θ is nearly constant constant.. Most standard standard lattice Boltzmann Boltzmann models use these smaller smaller velocity velocity sets. For models in an arbitrary number of dimensions this usually means that the third moment of the equilibrium distributio distribution n function function is modified to
f i0 (viα
i
− uα)(v )(viβ − uβ )(v )(viγ − uγ ) = nuα uβ uγ
(3.26)
and it it assumed that the term in u3 can be neglected. If it is not negligible this terms will lead to violations of Galilean invariance. Depending on the dimensionality of the space you want to simulate there are a variety of velocity sets and correspondin correspondingg equilibrium equilibrium distributions distributions that are frequently frequently used in the literature. literature. Equiblibrium distribution is given by: f i0
= nwi
3 9 1 + 2 u.vi + 4 (u.vi )2 c 2c
−
3 u.u . 2c2
(3.27)
The weights wi depend depend on the set of velocit velocites. es. The values values for com common monly ly used used models models are given given below. For D2Q9 we have 2
6
5
3
wi =
1 0
7
4
8
For D3Q15 the weights are:
4/9 i = 0 1/9 i = 1, 1 , 2, 3, 4 1/36 i = 5, 5 , 6, 7, 8
(3.28)
2/9 i = 0 1/9 i = 1 1/72 i = 7
−6 − 14
(3.29)
1/3 i = 0 1/18 i = 1 1/36 i = 7
−6 − 18
(3.30)
wi =
For D3Q19 the weights are: 2 6 3
1 5 4
wi =
17
3.5. NON-IDE NON-IDEAL AL FLUID FLUIDS S
For D3Q27 the weights are:
wi =
Note:
8/27 2/27 1/216 1/54
i =0 i =1 6 i = 7 14 i = 15 26
− − −
(3.31)
The moments moments of the equilibrium equilibrium distribution distribution function function (3.10) to (3.12) (3.12) and (3.26) (3.26) are not sufficient to determine the wi for large velocity velocity sets. You can see this easily if you consider consider that the weights weights to D3Q15 are a subset of D3Q27. D3Q27. You could simply set the weights weights such that you recover recover the D3Q15 model from D3Q27 and you would still have the same moments. These weights are determined by considerations that go beyond the Taylor expansion presented here.
Problems 3.4.1: Show that the constraint (3.26) is not Galilean invariant. 3.4.2: Calculate the weights wi in (3.27) for the D1Q3 model.
3.5 3.5
NonNon-id idea eall fluid fluidss
We discussed that the Boltzmann equation leads to the Navier-Stokes equation for an ideal gas. It would, of course, be much more useful if we could also simulate non-ideal systems. For such systems the Navier-Stokes equation is given by n∂ t uα + nuβ ∂ β uα =
2 −∂ β P αβ αβ + ∂ β [η (∂ β uα + ∂ α uβ − ∂ γ uγ δαβ )] 3
(3.32)
Comparing this to the ideal Navier-Stokes equation (2.21) we see that we can formally recover the non-ideal non-ideal behavior if we choose nF α = ∂ β (P αβ nθδαβ ) (3.33) αβ
−
−
for any given non-ideal pressure tensor P αβ αβ . These pressure tensors for an iso-thermal system can be derived from the Free energy of the system. You can implement a forcing term F by F i0
= nwi
3 9 F α viα + 4 (F α uβ + F β uα )viα viβ 2 c 2c
−
3 F α uα . c2
(3.34)
This description of the simulation of non-ideal fluids is quite superficial. In particular we have only considered the bulk terms of the pressure. There are additional terms which relate to the surface surface tension. Introducin Introducingg a forcing forcing of the form of (3.33) will lead to a surface surface tension, tension, but you can not derive what that surface tension is from the analysis presented here. The terms needed are formally of higher order (like 3 n). But since these gradiens are not small near the interface they are important for determining the shape of the interface and the surface tension. tension. Up to know know nobody seems to have have a way of consisten consistently tly deriving the hydrodynamic hydrodynamic terms in a way that includes all these terms and work in this area is ongoing. Notice:
∇
Problems 3.5.1: Show that (3.34) fulfills the equations (3.2)–(3.4). Hint:
Consider the moments of (3.27).
18
CHAPTER CHAPTER 3. LATTICE LATTICE BOLTZMANN BOLTZMANN
3.5.2: What are the forcing terms F i for a D1Q3 model. 3.5.3: Write a D1Q3 model for a Van der Waals gas and show that you can observe phase-
separation into a liquid and a gas phase. The Van der Waals equation is given by (V
−
a b) P + 2 V
= N kT
(3.35)
where a and b are parameters related to the attraction of the particles and the excluded volume volume respectively respectively.. We make the assumption assumption that we can define a local pressure that is related to the local density by n = N/V where N would be the number of particles at the lattice size and V the volume of the lattice size. a) Show Show that the Pressure is also given by nkT P = 1 nb/N
−
−
an2 N 2
(3.36)
b) Calculate Calculate the critical critical temperature temperature T c, the critical density nc and the critical pressure P c . 2 Hint: Consider that at the critical point ∂ V V P = 0 and ∂ v P = 0. c) Show that b/N = 1/(3n (3nc ) and a/N 2 = 9/ 9 /8 kT c /nc . Hence you can use nc and T c as the variable variabless in your program. (This makes makes it easier easier for you to decide decide whether you would expect phase-separ phase-separation ation in a given given system). Calculate Calculate (3.33) using the approximati approximation on ∂ x P = (P ( P [[x + 1] P [ P [x 1])/ 1])/2.
−
−
d) Calculate Calculate the numerical numerical phase diagram (i.e. nliquid (T/Tc) T/Tc) and ngas (T/Tc)) T/Tc)) using your program. What do you observe about the accuracy of the method? 3.5.4: You will have notice that the approach in the previous problem was not very accurate.
The reason for this is that we only used a first order accurate discretization for the pressure for our second order accurate accurate lattice Boltzmann Boltzmann method. What is a second second order accurate accurate definition of the first-order derivative of P? Use your new definition in your code and calculate the phase-diagram again.
3.6 3.6
Boun Bounda dari ries es
In many cases cases we will want want to implem implemen entt some some boundar boundaries ies in the fluid. It is usuall usually y assume assumed d that that there there is a non-sl non-slip ip boundar boundary y condit condition ion desired desired at those boundari boundaries. es. The simplest simplest way of implementing such a boundary is to draw the boundary and then mark all links that are cut by this boundary. boundary. Instead Instead of free streaming streaming on these these links the densities densities are “bounced-ba “bounced-back” ck” i.e. f i (x, t + 1) = f −i (x, t)
(3.37)
where the velocity index i is defined through v−i = vi . The effective boundary then lies halfway between the links. If the b oundary oundary is mov moving, ing, we clearly need to modify this b oundary oundary condition. condition. But what will it be modified modified to? The simple simplest st argume argument nt I could could come up with is the follow following ing:: we perform perform a Galilean transformation of the distribution into the rest frame of the boundary, then we perform the bounce-back operation, and then we transform the flow back into the original frame of reference 2 . Let us define a Galilean transform as
−
−
f ′ (vi ) = f i (vi ) + G(vi , U ) U )
(3.38)
We can then write the moving bounce back boundary condition as ′ f ( f (vi , t + 1) = f − U ) = f −i + G(vi , U ) U ) + G( vi , U ) U ) i + G( vi , U )
− −
2
− −
This derivation derivation follows roughly my argument from a paper on Lees-Edwa Lees-Edwards rds boundary conditions[2]. conditions[2].
(3.39)
19
3.6. BOUND BOUNDARIE ARIES S
For the moments of the Galilean transformation we obtain
(G(vi , U ) U ) + G( vi , U ) U )) = 0
− −
i
− −
(G(vi , U ) U ) + G( vi , U )) U ))vviα
i
− −
(G(vi , U ) U ) + G( vi , U )) U ))vviα viβ
i
− −
(G(vi , U ) U ) + G( vi , U )) U ))vviα viβ viγ
i
(3.40)
= 2nU α
(3.41)
= 0(?)
(3.42)
= ?
(3.43)
The usual choice is
6 wi nU α viα c2 So for each density crossing the moving boundary the streaming step is replaced with
− −
G(vi , U ) U ) + G( vi , U ) U ) =
f i (x, t + 1) = f −i (x, t) +
6 wi nU α viα . c2
(3.44)
(3.45)
As a practical matter it may often be more efficient to perform a streaming step for all densities and then and additional additional step that corrects corrects for the solid boundary boundary conditions. conditions. Assume Assume that a link from point (x, (x, y ) to (x (x + vix , y + viy ) is a boundary link associated with a boundary moving with velocity U. If f ˆi are the densities after the streaming step you will want to perform a swap of the densities and add the appropriate velocity terms: f i (x + vi , t + 1) = f −i (x, t + 1) = Need to double check this formula for signs.
6 f ˆ−i (x, t) + 2 wi nU α viα , c 6 f ˆi (x + vi , t) wi nU α viα . c2
−
(3.46) (3.47)
20
CHAPTER CHAPTER 3. LATTICE LATTICE BOLTZMANN BOLTZMANN
Chapter 4
Applications of Lattice Boltzmann 4.1 4.1
Turbu urbulen lence ce
4.2 4.2
Liqu Liquid id-g -gas as syst systam amss
4.3 4.3
TwoTwo-flu fluid id syst system emss
21
22
CHAPTER CHAPTER 4. APPLICATIONS APPLICATIONS OF LATTICE LATTICE BOLTZMANN BOLTZMANN
Chapter 5
Non-traditional Lattice Boltzmann methods 5.1 5.1
Intr Introdu oduct ctio ion n
While lattice Boltzmann methods were derived in the context of Fluid Mechanics the general algorith gorithm m can be used used to simula simulate te a variet ariety y of other other physi physical cal phenomen phenomena. a. We will focus here on methods that do not conserve momentum and therefore lead to a different set of equations.
5.2 5.2
No cons conser erv ved para parame meter ter
Let us first consider what we obtain if we use a lattice Boltzmann methods that has no conservation laws. A physical realization of such a system may be given by a magnetic system. We will identify the magnetization as
f i (x, t) = m(x, t)
(5.1)
i
If we now choose an equilibrium distribution with
f i0 = m +
i
M,
f i0 vi = j,
i
f i0 vi vi = ψ1,
(5.2)
i
we obtain for the first moment
∇− − ∇
∂ t m + j
5.2.1 5.2.1
τ
1 2
2
ψ=
M + O(∂ 3).
(5.3)
Magne Magnetic tic Syste Systems ms
The typical evolution equation for a magnetic system (also known as modelA [ ?]) is given by .∂ t m =
−µ
(5.4)
A simple form for µ for a ferromagnetic system is given by µ(m) = Am + Bm 3
− κ∇2m
(5.5)
This suggests a choice of j = 0, = (Am+ Am + Bm 3 ) and ψ = km/( km/(τ 0.5). This actually works! This should now be compare compare to Ising model simulation simulations. s. This also provides provides a great great opportunity opportunity to introduce the noise terms so that we can compare the results to a well studied system.
M −
−
23
24
5.3
CHAPTER CHAPTER 5. NON-TRADITIONA NON-TRADITIONAL L LATTICE LATTICE BOLTZMANN BOLTZMANN METHODS METHODS
Conse Conserv rved ed zero zeroth th momen momentt
= 0. The resulting equation of motion is
If we want to conserve the first moment we need to set
∇
∂ t m + j =
5.3. 5.3.1 1
Diffu Diffusi sion on
5.3.2 5.3.2
Electr Electrost ostati atics cs
− ∇ τ
1 2
2
ψ + O(∂ 3 ).
(5.6)
Chapter 6
Possible Projects Think of any problem you would would like to solve. solve. We can then decide whether whether and how you can use what you have learned to simulate the system you are interested in. There is a very large number of problems that can be addressed! Past Past projects include:
Jay Ihry Mixing of Milk in coffee.
25
26
CHAPTER CHAPTER 6. POSSIBLE PROJECTS PROJECTS
Bibliography [1] “Statistic “Statistical al Mechanics”, Mechanics”, 2nd edition, edition, K. Huang, Huang, Wiley. Wiley. [2] “Lees-Edw “Lees-Edwards ards boundary conditions conditions for lattice lattice Boltzmann Boltzmann”, ”, A.J. Wagner and I. Pagonaba Pagonabarrraga, J. Stat. Phys. 107, 521 (2002), [also cond-mat/0103218].
27
28
BIBLIOGRAPHY
Appendix A
Einstein convention for vectors and tensors We can write a vector in terms of its components
x1 x2 x3
x=
(A.1)
The scalar product of two vectors can then be written as c = x.y =
x1 x2 x3
y1 y2 y3
.
3
= x1 y1 + x2 y2 + x3 y3 =
xi yi
(A.2)
i=1
The product of a two dimensional tensor (also known as a matrix) A and a vector v is defined as w = A.v =
a11 a21 a31
a12 a22 a32
a13 a23 a33
v1 v2 v3
.
=
a11 x1 + a12 x2 + a13 x3 a21 x1 + a22 x2 + a23 x3 a31 x1 + a32 x2 + a33 x3
The product of two two-dimensional tensors A and B is given by C = A.B =
a11 a21 a31
a12 a22 a32
a13 a23 a33
.
b11 b21 b31
b12 b22 b32
b13 b23 b33
=
a1i bi1 i a2i bi1 i a3i bi1 i
a1i bi2 i a2i bi2 i a3i bi2 i
.
(A.3)
i a1i bi3 i a2i bi3 i a3i bi3
(A.4) As we progress to higher dimensional tensors it becomes more and more cumbersom to write down down all these componen components. ts. There is an easier way of writing this in terms of components. components. For this we use Greek Greek letters letters to enumera enumerate te the spatial spatial dimensions dimensions.. For (A.1) we write xα
(A.5)
where the free index α indicates that this is a vector expression. For (A.2) we get 3
c=
xα yα = xα yα
(A.6)
α=1
where we used the Einstein convention which tells us that a sum over repeated indices is implied. Note that this expression is a scalear expression, i.e. there are no free indices. indices. For the product of a matrix and a vector (A.3) we again obtain a vector equation 3
wα =
aαβ vβ = aαβ vβ
β =0
29
(A.7)
30
APPENDIX APPENDIX A. EINSTEIN EINSTEIN CONVENTION CONVENTION FOR VECTORS VECTORS AND TENSORS
The product of two matrices is a matrix equation ( i.e. it has two free indices) and (A.4) can be written written as cαβ = aαγ bγβ = aαγ bγβ (A.8)
γ
One important important matrix is the identity identity matrix. matrix. This is represeted represeted by the Kroneck Kronecker er delta in tensor tensor notation: 1 if α = β δαβ = (A.9) 0 othe otherw rwise ise
This often is appears in tensor equations. Consider for instance aαγ bǫβ δνγ δǫν = aαν bνβ .
(A.10)
You may want to point out that the vector notation is at least as compact as the tensor notation. That is correct, but it becomes difficult to manipulate expressions like tr ( u + ( u)T +
∇
∇
∇.u1).σ
.
(A.11)
(∂ α uβ + ∂ β uα + ∂ γ uγ δαβ )σαβ
(A.12)
It is much easier to manipulate the equivalent
if one is familiar with the tensor notation. notation. And for tensors of more than two dimensions dimensions vector vector notation is no longer used!
Problems A.1: Express v.A.w, where v and w are vectors and A is a two dimensional tensor, in tensor
notation. A.2: Are the following a scalar, vector, two-dimensional tensor or higher dimensional tensor ex-
pressions? a) aαβ bβγ cγδ b) aαβ bα cγ c) aαβγ bαδ cγδ d) aα bγδǫ cν eδ A.3: Simplify ∂ γ uα δβγ + ∂ γ uβ δαγ .
Appendix B
Functional Derivatives Recall that the derivative of a function f ( f (x) is defined as df (x) f ( f (x + dx) dx) = lim dx→0 dx dx
− f ( f (x)
(B.1)
For a functional, i.e. a function of a function we can define a functional derivative in a similar way: δ Ψ[f Ψ[f ((x)] Ψ[f Ψ[f ((x) + δf ( δf (y )] = l im δf ( δf (y) δf ( δf (y ) δf ( δf (y )→0
− Ψ[f Ψ[f ((x)]
(B.2)
where the variation δf ( δf (y) is only different from zero at x = y . These functionals are usually integrals and in the simplest case where Ψ is only a functional of f and not its derivatives we have b
Ψ[f Ψ[f ]] =
ψ(f (x))dx ))dx
(B.3)
a
The functional derivative is then simply given by b
δ Ψ[f Ψ[f ((x)] = δf ( δf (y)
a
δf ( δf (x) ψ (f ( f (x)) dx = δf ( δf (y) ′
b
ψ′ (f ( f (x))δ ))δ(x
a
′
− y)dx = ψ (y)
(B.4)
where ψ′ is the ordinary derivative of ψ. The functional derivative of functionals containing spatial derivativ derivatives es is a little more complicated, complicated, and we need to perform a partial partial integration. integration. For Ψ[ f ] f ] = b f (x). f ( f (x) we obtain a
∇
∇
δ Ψ[f Ψ[f ((x)] = δf ( δf (y)
b
a
δ f ( f (x) 2 f ( f (x) dx = δf ( δf (y )
∇
∇
b
− ∇ a
2
f ( f (x)δ(x
− y)dx = −∇2f ( f (y )
(B.5)
where we have made use of the fact that the boundary term vanishes because of the δ-function. Other derivative terms have to be treated similarly.
Problems B.1: What is the functional derivative of
a) Ψ[f Ψ[f ]] =
b 2 (f (x) a
− f 5(x) + f 2(x)∇2 f (x))dx ))dx ?
31
32
APPENDIX APPENDIX B. FUNCTIONAL FUNCTIONAL DERIVA DERIVATIVES
Appendix C
Lagrangian multipliers If you want to minimize a functional while you have other constraints you can often use Lagrangian multipliers to remove the contraints. b Assume Assume you want want to minimi minimize ze the functi functiona onall Ψ[φ Ψ[φ(x)] = a ψ (φ(x))dx ))dx with the constrain constraintt b ˆ as φ(x) = c. We then write a new functional Ψ a
ˆ = Ψ + λ( Ψ
b
φ(x)dx
a
− c)
(C.1)
ˆ if the constraint where λ is known known as a Lag Lagran range ge multipli multiplier. er. Note Note that that Ψ = Ψ constraint is fulfilled. Now, ˆ however, we relax the constraint and minimize Ψ by calculating 0=
δΨ(φ Ψ(φ) = ψ ′ (φ(y )) + λ δφ( δφ(y )
(C.2)
which we can invert to find φ(y, λ). These These soluti solutions ons will, in genera general, l, not fulfill fulfill our constra constraint int.. However, we can introduce a “cost” of not obeying the constraint by setting different values for λ. In particular we can determine the value of λ where this “cost” is just right so that φ(x, λ) does fulfill the constraing. We do this by inserting φ(x, λ) into the constraint b
φ(y, λ) = c.
(C.3)
a
We can solve this to determine the Lagrangian multiplier λ for which our solution obeys the constraint.
Problems C.1: Minimize
a
− a
−
a
1 2 1 1 f (x) + f 4 (x) + 2 4 2
with the constraint −a f (x) = 0. tanh(x/ 2). Solution: f ( f (x) = tanh(x/
√
33
df ( f (x) dx
2
(C.4)
34
APPENDIX APPENDIX C. LAGRANGIAN LAGRANGIAN MULTIPLIERS MULTIPLIERS
Appendix D
Evaluating Gaussian integrals First let us review some standard one dimensional Gaussian integrals that we will be using:
∞
π a 1 π a−3/2 2 3 π a−5/2 4
2
e−ax dx =
−∞
√
∞
2
x2 e−ax dx =
−∞
√
∞
2
x4 e−ax dx =
−∞
Now let us look at some of the more relevant three dimensional integrals that we can derive using the above equations. Integrating the Gaussian we get ∞
∞
∞
−∞
−∞
∞
dx dy dz e
a(x2 +y2 +z 2 )
−
=
−∞
∞
ax2
−
e
dx
−∞
π a
= To calculate the temperature we need ∞
∞
∞
e
−∞
3/2
−∞
.
(D.1)
2
∞
ax2
−
x e
e
dy
−∞
∞
e
ax
−
y e
2
−∞
−∞
∞
2
e−ay dy
−∞
1 π 2a a
3/2
2
e−az dz
dy
∞
e−ax dx
+
ay
−
2
−∞
∞
2
e−az dz
∞
2
dx
−∞
∞
2
−∞
= 3
∞
ay 2
−
dx
−∞
+
+y 2 +z 2 )
−∞
∞
=
2
z 2 e−az dz
−∞
.
(D.2)
To calculate the heat equation we need to evaluate terms of the form ∞
∞
∞
2
dx dy dz x2 (x2 + y2 + z 2 )e−a(x
−∞
−∞
−∞
∞
=
∞
4
x e
ax2
−
dx
−∞
+
∞
e
ax2
−
dx
−∞
−∞
∞
∞
2
x2 e−ax dx
−∞
dy
2
e−az dz
−∞
∞
x e
+
ay2
−
−∞
∞
2
2
e−az dz
−∞
dx dy dz (x2 + y 2 + z 2 )e−a(x
2
dy
∞
−∞
ay 2
−
2
y e
ay2
−
∞
2
e−az dz
dy
−∞
∞
2
e−ay dy
−∞
−∞
35
2
z 2 e−az dz
+y2 +z 2 )
36
APPENDIX APPENDIX D. EVALUA EVALUATING TING GAUSSIAN INTEGRALS
1 (2a (2a)2 1 = 5 (2a (2a)2 = 3
π a π a
3/2
3/2
+
1 (2a (2a)2
π a
3/2
+
1 (2a (2a)2
π a
3/2
(D.3)
Note that in the Gaussian integrals we perform a = 1/(2θ (2θ).
Problems D.1: Show what the value of the integrals (D.1), (D.2), and (D.3) are in a D dimensional space.
Appendix E
Galilean transformations Since the velocity set is fixed we clearly can not perform a true Galilean transformation which would simply require f ′ (vi ) = f ( f (vi + U ). U ). What we can require, however, is that the transformed moments agree agree with the Galilean Galilean transform transformation, ation, i.e. f i′ viN = f i (vi + U ) U )N for an appropriate number of powers N . N . This gives us
f i′
=
f i′ viα
=
f i′ viα viβ
=
f i′ viα viβ viγ
=
f i′ viα viβ viγ
f i
(E.1)
f i viα + nU α
(E.2)
f i viα viβ + n(uα U β + U α uβ + U α U β )
(E.3)
f i viα viβ viγ + n[(u [(uα uβ + θδαβ )U γ + (u (uα uγ + θδαγ )U β + (u (uβ uγ + θδ βγ )U α
+U α U β U γ ]
=
(E.4)
f i viα viβ viγ + n[(u [(uα uβ + θδαβ )U γ + (u (uα uγ + θδαγ )U β + (u (uβ uγ + θδ βγ )U α
+U α U β U γ ]
(E.5)
37
38
APPENDIX APPENDIX E. GALILEAN TRANSFORMA TRANSFORMATIONS TIONS
Appendix F
D1Q5 Lattice Boltzmann code 1
/* An implem implement entati ation on of the D1Q5 D1Q5 model model for non-id non-ideal eal system systems s */
2 3 4 5 6
#include
#include #include "mygraph.h" #define #define xdim 100
7 8 9 10 11 12
double f0[xdim],f1[xdim],f2[xdi f0[xdim],f1[xdim],f2[xdim],f3[xdim],f m],f3[xdim],f4[xdim],n[xdi 4[xdim],n[xdim]; m]; double double n0=1, n0=1, T0=0.333 T0=0.3333333 33333, 3, Amp=0.01 Amp=0.01, , omega=1; omega=1; double ug[xdim],pg[xdim],tg[xdi ug[xdim],pg[xdim],tg[xdim]; m]; int ugreq=0,pgreq=0,tgreq=0; ugreq=0,pgreq=0,tgreq=0; int next=0,pause=1,done=0,Rep next=0,pause=1,done=0,Repeat=1,iterati eat=1,iterations; ons;
13 14 15 16 17 18 19 20 21 22 23 24 25
void init(){ init(){ int int i; iterations=0; for (i=0;i
26 27 28 29
void iteration(){ double u,T,tmp,tmp2; u,T,tmp,tmp2; int int i;
30 31 32 33 34 35 36 37 38 39 40
iterations++; for (i=0;i
39
40
} tmp=f1[0]; memmove(&f1[0],&f1[1],(xdi memmove(&f1[0],&f1[1],(xdim-1)*sizeof(d m-1)*sizeof(double)); ouble)); f1[xdim-1]=tmp; tmp=f2[xdim-1]; memmove(&f2[1],&f2[0],(xdi memmove(&f2[1],&f2[0],(xdim-1)*sizeof(d m-1)*sizeof(double)); ouble)); f2[0]=tmp; tmp=f3[0]; tmp2=f3[1]; memmove(&f3[0],&f3[2],(xdi memmove(&f3[0],&f3[2],(xdim-2)*sizeof(d m-2)*sizeof(double)); ouble)); f3[xdim-1]=tmp2; f3[xdim-2]=tmp; tmp=f4[xdim-1]; tmp2=f4[xdim-2]; memmove(&f4[2],&f4[0],(xdi memmove(&f4[2],&f4[0],(xdim-2)*sizeof(d m-2)*sizeof(double)); ouble)); f4[0]=tmp2; f4[1]=tmp;
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
APPENDIX APPENDIX F. D1Q5 LATTICE LATTICE BOLTZMANN BOLTZMANN CODE
}
19 20 21
void GUI(){ GUI(){ static static int xdimi=xd xdimi=xdim; im;
22
DefineGraphN_R("n",&n[0],& DefineGraphN_R("n",&n[0],&xdimi,NULL); xdimi,NULL); DefineGraphN_R("u",&ug[0], DefineGraphN_R("u",&ug[0],&xdimi,&ugreq &xdimi,&ugreq); ); DefineGraphN_R("p",&pg[0], DefineGraphN_R("p",&pg[0],&xdimi,&pgreq &xdimi,&pgreq); ); DefineGraphN_R("T",&tg[0], DefineGraphN_R("T",&tg[0],&xdimi,&tgreq &xdimi,&tgreq); ); StartMenu("D1Q3",1); DefineInt("iterations",&it DefineInt("iterations",&iterations); erations); DefineDouble("omega",&omega); DefineDouble("Amp",&Amp); DefineDouble("n0",&n0); DefineDouble("T0",&T0); DefineFunction("init",&init); DefineGraph(curve2d_,"Graphs"); DefineInt("Repeat",&Repeat); DefineBool("next",&next); DefineBool("pause",&pause); DefineBool("done",&done); EndMenu();
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
}
41 42 43
void GetData( GetData(){ ){ int int i;
44 45 46 47 48 49 50 51 52 53 54 55
if (ugreq|| (ugreq||tgre tgreq) q) { for (i=0;i
41 }
1 2
}
3 4 5 6
int main(int main(int argc, char *argv[]){ *argv[]){ int newdata= newdata=1; 1; int int i;
7
init(); GUI();
8 9 10
while (done==0){ Events(newdata); GetData(); DrawGraphs(); if (next|| (next|| !pause){ !pause){ newdata=1; next=0; for (i=0;i
11 12 13 14 15 16 17 18 19 20 21 22 23 24
return return 0;
25 26
}
