Lattice Gas Cellular Automata And Lattice Boltzmann Models Chapter7

6. Appendix

D.A. Wolf-Gladrow: LNM 1725, pp. 247–270, 2000. c Springer-Verlag Berlin Heidelberg 2000


248

6. Appendix

6.1 Boolean algebra George Boole (1815 - 1864), English mathematician. Boolean algebra is a mathematical structure B = (B, ∪, ∩,¯) consisting of a set B and two binary operations called union (∪) and intersection (∩) and one unitary operation called complementation (¯). The system which originally inspired this collection of laws is the algebra of sets with the familiar operations of set union and intersection. The following relations1 hold true: 1. Closure: (i) The union of two elements in B yields a unique element in B. a, b ∈ B → (a ∪ b) ∈ B (ii) The intersection of two elements in B yields a unique element in B. a, b ∈ B → (a ∩ b) ∈ B 2. Commutativity: ∀ a, b ∈ B: (i) a ∪ b = b ∪ a (ii) a ∩ b = b ∩ a 3. Associativity: ∀ a, b, c ∈ B: (i) (a ∪ b) ∪ c = a ∪ (b ∪ c) (ii) (a ∩ b) ∩ c = a ∩ (b ∩ c) 4. Distributivity: ∀ a, b, c ∈ B: (i) a ∩ (b ∪ c) = (a ∩ b) ∪ (a ∩ c) (ii) a ∪ (b ∩ c) = (a ∪ b) ∩ (a ∪ c) 5. The idempotent laws: ∀a ∈ B (i) a ∪ a = a (ii) a ∩ a = a 6. Identity elements 0 and 1: (i) In B there is the unique element 0 with the following properties: a ∪ 0 = a and a ∩ 0 = 0 ∀a ∈ B. 0 is the identity element with respect to union. (ii) In B there is the unique element 1 with the following properties: a ∪ 1 = 1 and a ∩ 1 = a ∀a ∈ B. 1 is the identity element with respect to intersection. Remark: If you consider the ordinary union and intersection of sets the 0-element corresponds to the empty set and the 1-element corresponds to the whole set. 1

We do not speak of axioms because some of relations are redundant.

6.1 Boolean algebra

249

7. Complementation: ∀ a ∈ B there exists a unique element a such that (i) a ∪ a = 1 (ii) a ∩ a = 0. a is called the complement of a. The elements a, b ∈ B obey the DeMorgan laws: (iii) a ∪ b = a ∩ b (iv) a ∩ b = a ∪ b Law of involution: ∀a ∈ B (v) (a) = a Principle of duality: The substitution (∪, ∩, 0, 1)

→

(∩, ∪, 1, 0)

transforms true expressions of the Boolean algebra into true expressions. Proof: All laws listed above obey the principle of duality. Law 7(v) is selfdual. Exercise 6.1.1. (**) Show that the structure Btf = (B = {true, f alse}, AN D, (inclusive)OR, N OT ) defines a Boolean algebra. Show that XOR can be expressed in terms of AN D, OR and N OT . Why is N = (B = {true, f alse}, AN D, XOR, N OT )(XOR = exclusive or) not a Boolean algebra? Exercise 6.1.2. (*) Consider the set B = {0, 1 ∈ IN} and the operations addition modulo 2 (+) and multiplication (·). Show: this structure obeys the laws of Boolean algebra if complementation is appropriately defined (how?). Exercise 6.1.3. (**) Prove that for all elements a and b in the set B of a Boolean algebra (B, ∪, ∩,¯): (a ∩ b) ∪ (a ∩ b) = a Exercise 6.1.4. (*) Consider the set of integers and the operations addition (+) and multiplication (·). Show that this structure is not a Boolean algebra.

250

6. Appendix

6.2 FHP: After some algebra one finds ... In this appendix the Lagrange multipliers of the equilibrium distributions for the FHP-I and the HPP model will be calculated (compare Subsection 3.2.5). For vanishing flow velocity the occupation numbers equal each other u=0

=⇒

Ni =

ρ =d b

(6.2.1)

(b number of cells per node, b = 6 for FHP-I, b = 4 for HPP; d density per cell) and therefore Ni (ρ, 0) =

1 =d 1 + exp[h(ρ, 0)]

(6.2.2)

and q(ρ, 0) = 0.

(6.2.3)

Because of invariance of the occupation numbers under the parity transform u

→

−u,

ci

→

−ci

(6.2.4)

it follows that h(ρ, −u) = h(ρ, u),

(6.2.5)

q(ρ, −u) = −q(ρ, u).

(6.2.6)

and The Lagrange multipliers h und q will be expanded up to second order in u: h(ρ, u)

= h0 (ρ) + h2 (ρ)u2 + O(u4 )

(6.2.7)

q(ρ, u)

= q1 (ρ)u + O(u ).

(6.2.8)

3

All other low order terms vanish because of the parity constraints. It is remarkable that h2 and q1 are scalars instead of tensors of rank 2. This fact is a consequence of the isotropy of lattice tensors (compare Section 3.3) of rank 2. The expansions (6.2.7) and (6.2.8) will be inserted into the Fermi-Dirac distribution (3.2.28). Then the distribution is expanded in a Taylor series with respect to u at u = 0 up to second order.

Ni (u)

∂Ni ∂Ni ·u+ ·v ∂u ∂v 2 2 1 ∂ Ni 2 ∂ Ni 1 ∂ 2 N1 2 + · u + · v + O(u3 ), · u · v + 2 ∂u2 ∂u∂v 2 ∂v 2

= Ni (u = 0) +

6.2 FHP: After some algebra one finds ...

Ni (u) =

1 , 1 + exp[x(u)]

x(u) = h0 + h2 u2 + q1 uci , x(0) = h0 , 1 = d, Ni (h0 ) = (1 + exp[h0 ])

→

1−d , d 1−d ln d ∂Ni ∂x · ∂x ∂uα

exp[h0 ]

=

h0

=

∂Ni ∂uα

=

∂Ni ∂x

=

−

=

2h2 uα + q1 ciα → q1 ciα , % & ∂Ni ∂x ∂ · ∂uα ∂x ∂uα ∂x ∂Ni ∂ 2 x ∂ 2 Ni · + · ∂x∂uα ∂uα ∂x ∂u2α !2 ∂ 2 Ni ∂x ∂Ni ∂ 2 x · + ∂x2 ∂uα ∂x ∂u2α

∂x ∂uα ∂ 2 Ni ∂u2α

= = =

∂ 2 Ni ∂x2

∂2x ∂u2α

→  u=0

exp[x] (1 + exp[x])2

d(d − 1)q1 ciα , →  u=0

d−1 2 d = d(d − 1) d

→ d(d − 1)(2d − 1)q12 c2iα + d(d − 1)2h2 , exp[x](1 + exp[x])2 − exp[x] · 2(1 + exp[x]) exp[x] = − (1 + exp[x])4 exp[x](exp[x] − 1) = (1 + exp[x])3 1−d 1−d ( − 1)d3 = d(d − 1)(2d − 1), → d d =

2h2 .

For α = β: ∂ 2 Ni ∂uα ∂uβ

%

∂Ni ∂x · ∂x ∂uα

&

=

∂ ∂uβ

=

∂ 2 Ni ∂x ∂Ni ∂2x · + · ∂x∂uβ ∂uα ∂x ∂uα ∂uβ

251

252

6. Appendix

∂2x ∂ 2 Ni ∂x ∂x ∂Ni · · · + 2 ∂x ∂uα ∂uβ ∂x ∂uα ∂uβ

=

→ d(d − 1)(2d − 1)q12 ciα ciβ , ∂2x ∂uα ∂uβ

Ni (u)

=

0,

= d + d(d − 1)q1 ci · u 1 + d(d − 1)(2d − 1)q12 c2iα u2α + d(d − 1)h2 u2 . 2

(6.2.9)

At this point the coefficients h2 and q1 are not known yet. The expression for Ni , however, looks much simpler now: a polynomial instead of a rational function with an exponential function in the denominator. Next insert the Ni according to eq. (6.2.9) into the definitions of mass and momentum density and use the moments relations (3.2.3 - 3.2.5):  Ni ρ = i

=



d+

i



d(d − 1)q1 ci · u

i

      =ρ =0 1 d(d − 1)(2d − 1)q12 c2iα u2α + 2 i    1 = 3d(d − 1)(2d − 1)q12 u2 2  + d(d − 1)h2 u2 . i

   = 6d(d − 1)h2 u2 From this a relation between h2 and q1 follows: h2 ρu

1 (1 − 2d)q12 , 4  = Ni ci

=

i

=

 i

dci +

   =0

 

i

d(d − 1)q1 (ci · u)ci  = 3d(d − 1)q1 u



6.2 FHP: After some algebra one finds ...

+

1 

+

 

q1

=

h2

=

Ni (u)

i

2

d(d − 1)(2d − 1)q12 c2iα u2α ci  =0

d(d − 1)h2 u2 ci ,

i

 =0

with G(ρ) =



2 , d−1 1 − 2d , (d − 1)2

= d + 2dci · u + 2d =



1 − 2d 2 2 1 − 2d 2 c u −d u 1 − d iα α 1−d

ρ ρ + ci · u + ρG(ρ)Qiαβ uα uβ , 6 3

1 6 − 2ρ 3 6−ρ

1 and Qiαβ = ciα ciβ − δαβ . 2

253

254

6. Appendix

6.3 Coding of the collision operator of FHP-II and FHP-III in C

/* collision with 1 rest particle (FHP-II) 17.5.89 dwg */ for(ix=0; ix < IXM; ix++) { for(iy=0; iy < IYM; iy++) { /* i’s -> register a,...,f a b c d e f r

= = = = = = =

*/

i1[ix][iy]; i2[ix][iy]; i3[ix][iy]; i4[ix][iy]; i5[ix][iy]; i6[ix][iy]; rest[ix][iy];

/* three body collision <-> 0,1 (bits) alternating <-> triple = 1 */ triple = (aˆb)&(bˆc)&(cˆd)&(dˆe)&(eˆf); /* two-body collision <-> particles in cells a (b,c) and d (e,f) no particles in other cells <-> db1 (db2,db3) = 1 */ db1 = (a&d&˜(b|c|e|f)); db2 = (b&e&˜(a|c|d|f)); db3 = (c&f&˜(a|b|d|e)); /* rest particle and 1 particle ra rb rc rd re rf

= = = = = =

*/

(r&a&˜(b|c|d|e|f)); (r&b&˜(a|c|d|e|f)); (r&c&˜(a|b|d|e|f)); (r&d&˜(a|b|c|e|f)); (r&e&˜(a|b|c|d|f)); (r&f&˜(a|b|c|d|e));

/* no rest particle and 2 particles (i,i+2) */

6.3 Coding of the collision operator of FHP-II and FHP-III in C

ra2 rb2 rc2 rd2 re2 rf2

= = = = = =

255

(f&b&˜(r|a|c|d|e)); (a&c&˜(r|b|d|e|f)); (b&d&˜(r|a|c|e|f)); (c&e&˜(r|a|b|d|f)); (d&f&˜(r|a|b|c|e)); (e&a&˜(r|b|c|d|f));

/* change a and d <-> three-body collision triple=1 or two-body collision db1=1 or two-body collision db2=1 and eps=1 (- rotation) or two-body collision db3=1 and noeps=1 (+ rotation) <-> chad=1 */ eps = irn[ix][iy]; noeps = ˜eps;

/* random bits */

cha=(triple|db1|(eps&db2)|(noeps&db3)|ra|rb|rf|ra2|rb2|rf2); chd=(triple|db1|(eps&db2)|(noeps&db3)|rd|rc|re|rd2|rc2|re2); chb=(triple|db2|(eps&db3)|(noeps&db1)|rb|ra|rc|rb2|ra2|rc2); che=(triple|db2|(eps&db3)|(noeps&db1)|re|rd|rf|re2|rd2|rf2); chc=(triple|db3|(eps&db1)|(noeps&db2)|rc|rb|rd|rc2|rb2|rd2); chf=(triple|db3|(eps&db1)|(noeps&db2)|rf|ra|re|rf2|ra2|re2); chr=(ra|rb|rc|rd|re|rf|ra2|rb2|rc2|rd2|re2|rf2); /* change: a = a ˆ chad k1[ix][iy] = k2[ix][iy] = k3[ix][iy] = k4[ix][iy] = k5[ix][iy] = k6[ix][iy] = rest[ix][iy]

*/

i1[ix][iy]ˆcha; i2[ix][iy]ˆchb; i3[ix][iy]ˆchc; i4[ix][iy]ˆchd; i5[ix][iy]ˆche; i6[ix][iy]ˆchf; ˆ= chr;

/* collision finished (except at the boundaries) */ }}

256

6. Appendix

/* i======================================================i i i i c o l l i s i o n FHP-III i i ------------------i i i i ( two, three and four body collisions) i i 162 bit-operations i i======================================================i

collision with 1 rest particle (FHP-III) 24.6.91 Armin Vogeler

*/

for(ix=0; ix < IXM; ix++) for(iy=0; iy < IYM; iy++) {

a = b = c = d = e = f = r = s = eps

ns h1 h2 h3 h4 h5 h6 /* h0 /* z1

i1[ix][iy]; i2[ix][iy]; i3[ix][iy]; i4[ix][iy]; i5[ix][iy]; i6[ix][iy]; rest[ix][iy]; sb[ix][iy]; = irn[ix][iy];

= = = = = = = 0 = 2 =

˜s; /* no solid bit */ a&c&e; /* 3 particles a,c,e */ b&d&f; /* 3 particles b,d,f */ aˆcˆe; /* 1 or 3 particles in a,c,e */ bˆdˆf; /* 1 or 3 particles in b,d,f */ a|c|e; /* at least 1 particle in a,c,e */ b|d|f; /* at least 1 particle in b,d,f */ particles in a,c,e or!! in b,d,f and no solid */ ns&(h5ˆh6); particles in a,c,e or b,d,f and no solid */ ns&((˜h3&h5)ˆ(˜h4&h6));

/* three-body collisions */

6.3 Coding of the collision operator of FHP-II and FHP-III in C

257

c3 = (h1ˆh2)&h0; /* head-on collisions with spectator */ c2s = z1&((a&d)ˆ(b&e)ˆ(c&f)); /* two- and four-body collisions */ c24 = ns&((fˆa)|(aˆb))&(˜((fˆc)|(aˆd)|(bˆe))); /* rest particle and 1 particle collisions */ r1 = r&((h3&˜h1)ˆ(h4&˜h2))&h0; /* no rest particle and 2 particles collisions */ r2 = z1&h0&˜r; /* no no le ri

s,c24, r1, r2 and c2s collision */ = ˜(c24|r1|r2|c2s)&ns; = c24&eps; /* c24 collision and left rotation */ = c24&˜eps; /* c24 collision and right rotation */

/* change bitfield: */ /*|---| |---------| |---------| |------| |--------| |-----|*/ /*| | | | | |----|| | | ||---| | | |*/ k1[ix][iy] =(s&d)|(ri&bˆle&f)|(no&(aˆc3))|(b&f&r2)|((bˆf)&r1)|(˜d&c2s); k2[ix][iy] =(s&e)|(ri&cˆle&a)|(no&(bˆc3))|(a&c&r2)|((aˆc)&r1)|(˜e&c2s); k3[ix][iy] =(s&f)|(ri&dˆle&b)|(no&(cˆc3))|(b&d&r2)|((bˆd)&r1)|(˜f&c2s); k4[ix][iy] =(s&a)|(ri&eˆle&c)|(no&(dˆc3))|(c&e&r2)|((cˆe)&r1)|(˜a&c2s); k5[ix][iy] =(s&b)|(ri&fˆle&d)|(no&(eˆc3))|(d&f&r2)|((dˆf)&r1)|(˜b&c2s); k6[ix][iy] =(s&c)|(ri&aˆle&e)|(no&(fˆc3))|(e&a&r2)|((eˆa)&r1)|(˜c&c2s); krest[ix][iy] = rˆ(r1|r2); } /*-----------

end of collision

---------- */

258

6. Appendix

6.4 Thermal LBM: derivation of the coefficients Constraints from the definitions of mass and momentum. Mass: ρ=



Fieq = A0 + 6 (A1 + A2 ) + (3C1 + 12C2 + D0 + 6D1 + 6D2 ) u2

i

→ constraints 1 and 2 3C1 + 12C2 + D0 + 6D1 + 6D2 = 0

(6.4.1)

A0 + 6 (A1 + A2 ) = ρ

(6.4.2)

Momentum: %  eq j= ci Fi = u 3B1 + 12B2 + i

! & 9 E1 + 36E2 + 3G1 + 12G2 u2 4

→ constraints 3 and 4 3B1 + 12B2 = ρ 9 E1 + 36E2 + 3G1 + 12G2 = 0 4

(6.4.3) (6.4.4)

Conservation of mass, momentum, and energy. The expansions (5.2.25) will be substituted into the conservation equations for mass, momentum and energy   1   0 =  ci  [Fi (x + ci , t + 1) − Fi (x, t)] i c2i /2 which lead to 0

=  (5.2.25)

=



 1  1 ' ciα   Fi (x, t) + ∂t Fi + ciβ ∂xβ Fi + ∂t ∂t Fi 2 i ciα ciα /2 ( 1 + ∂xβ ∂xγ ciβ ciγ Fi + ciβ ∂t ∂xβ Fi + O(∂ 3 Fi ) − Fi (x, t) 2   1

  (1) (0) 2 (1) (1) 2 (2) (0) ciα   ∂t Fi +  ∂t Fi +  ∂t Fi i ciα ciα /2 1 (0) (1) (1) (1) (0) +ciβ ∂x(1) Fi + ciβ 2 ∂x(1) Fi + 2 ∂t ∂t Fi β β 2 & 1 (0) (1) 2 (1) (1) (0) 3 + 2 ∂x(1) ∂ c c F + c  ∂ ∂ F + O( ) iβ iγ iβ t xγ xβ i i β 2

6.4 Thermal LBM: derivation of the coefficients

259

and finally sorted according to orders in    1  /

  (1) (0) (1) (0) 0 = ciα    ∂t Fi + ciβ ∂xβ Fi i ciα ciα /2 % 1 (1) (1) (0) (1) (1) (2) (0) (1) Fi + ∂t ∂t Fi (6.4.5) +2 ∂t Fi + ∂t Fi + ciβ ∂x(1) β 2 &  1 (1) (1) (0) (1) (1) (0) 3 + O( ) + ∂xβ ∂xγ ciβ ciγ Fi + ciβ ∂xβ ∂t Fi 2

Terms of first order in : mass. To first order in  eq. (6.4.5) yields:   (1) (0) (0) ∂t Fi + ciα ∂x(1) Fi 0 = α i

or

(1)

∂t ρ + ∂x(1) j =0 α α

(continuity equation),

(6.4.6)

→ no further constraints from mass conservation. Terms of first order in : momentum.  

(1) (0) (0) ciα ∂t Fi + ∂x(1) c c F 0 = iα iβ i β 0

=

i (1) ∂t (ρuα )

whereby (0)

Pαβ :=

(0)

+ ∂x(1) Pαβ . β



(6.4.7)

(0)

ciα ciβ Fi

i

is the momentum flux tensor with components ! 3 (0) Pxx = 3A1 + 12A2 + C1 + 12C2 (3u2 + v 2 ) + (3D1 + 12D2 )u2 4 ! 9 C1 + 36C2 + 3D1 + 12D2 u2 = 3A1 + 12A2 + 4 ! 3 + C1 + 12C2 + 3D1 + 12D2 v 2 4

260 (0) Pxy (0) Pyy

6. Appendix

! 3 (0) C1 + 24C2 uv = Pyx (6.4.8) 2 !
 3 C1 + 12C2 u2 + 3v 2 + (3D1 + 12D2 ) u2 = 3A1 + 12A2 + 4 ! 3 C1 + 12C2 + 3D1 + 12D2 u2 = 3A1 + 12A2 + 4 ! 9 C1 + 36C2 + 3D1 + 12D2 v 2 + 4 =

The momentum flux tensor should yield # (0) Pαβ

= ρuα uβ + p δαβ = ρ

u2 uv

uv v2

$ + p δαβ .

(6.4.9)

Comparison of (6.4.8) and (6.4.9) leads to: 3A1 + 12A2 9 C1 + 36C2 + 3D1 + 12D2 4 3 C1 + 12C2 + 3D1 + 12D2 4 3 C1 + 24C2 2

= p = ρ =

0

= ρ.

This results in the three independent constraints 5 to 7:

and

3A1 + 12A2 = p 3 C1 + 24C2 = ρ 2

(6.4.10)

6D1 + 24D2 = −ρ

(6.4.12)

(6.4.11)

The first order terms of the moment equation lead to the Euler equation (1)

(ρuβ ) ρ∂t uα − uα ∂x(1) = β  (6.4.6) = (1)

∂t (ρuα )

(1)

(1)

ρ∂t uα + uα ∂t ρ

= =

(0)

Pαβ −∂x(1) β

= −∂x(1) (ρuα uβ + pδαβ ) β  (6.4.9) =

−ρuβ ∂x(1) uα − uα ∂x(1) (ρuβ ) − ∂x(1) p α β β

(6.4.13)

6.4 Thermal LBM: derivation of the coefficients

and therefore

(1)

uα − ∂x(1) p ρ∂t uα = −ρuβ ∂x(1) α β

261

(6.4.14)

Terms of first order in : energy.  1  (1) (0) (0) ∂t ciα ciα Fi + ∂x(1) 0= c c c F β iα iα iβ i 2 i 



 1 (1)  1 (1) (0) (0) (1)  1 ∂t ρuα uα + p  ciα ciα Fi = ∂t Pαα = ∂t    2 2 2   = i ρε I = ρεK where p is to be identified with the internal energy, i.e. p = ρεI , and 12 ρuα uα is the kinetic energy. 1 (1)  (0) ∂ ciα ciα ciβ Fi 2 xβ i 



  ! !  3   B1 + 24B2 uβ + 9 E1 + 72E2 + 3 G1 + 24G2 uα uα uβ  = ∂x(1)  β  8 2  2       =: f1 (ρ, εI ) =: f2 (ρ, εI ) = ∂x(1) [f1 (ρ, εI )uβ + f2 (ρ, εI )uα uα uβ ] β

(6.4.15)

and thus (1)

∂t (ρεK + ρεI ) = −∂x(1) [f1 (ρ, ε)uβ + f2 (ρ, ε)uα uα uβ ] β (1)

∂t (ρεK )

= = =  (6.4.13),(6.4.14)

(6.4.16)

(1) 1 ∂t ( ρuα uα ) 2 1 1 (1) (1) uα ∂t (ρuα ) + ρuα ∂t uα 2 2 1 − uα [ρuβ ∂x(1) uα + uα ∂x(1) (ρuβ ) + ∂x(1) p α β β 2

+ρuβ ∂x(1) uα + ∂x(1) p] α β =

−ρuα uβ ∂x(1) uα − uα ∂x(1) p α β 1 − uα uα ∂x(1) (ρuβ ) β 2

(6.4.17)

262

6. Appendix

(1)

ρ∂t εI

(1)

(1)

∂t (ρεI ) − εI ∂t ρ

=

(1)

[f1 uα + f2 uα uβ uβ ] + εI ∂x(1) (ρuα ) = −∂t (ρεK ) − ∂x(1) α α  (6.4.6), (6.4.16) 1 uα + uα ∂x(1) p + uα uα ∂x(1) (ρuβ ) = ρuα uβ ∂x(1) α β β  2 (6.4.17) −∂x(1) [f1 uα + f2 uα uβ uβ ] + εI ∂x(1) (ρuα ). α α

(6.4.18)

Substitution of p = ρεI and expansion of all terms leads to (1)

ρ∂t εI

= ρuα uβ ∂x(1) uα + ρuα ∂x(1) ε + uα εI ∂x(1) ρ α I α β          (2)

(1)

(3)

1 1 + ρuβ uβ ∂x(1) uα + uα uβ uβ ∂x(1) ρ α α 2 2       (4)

− f1 ∂x(1) uα α

  

(5)

− uα ∂x(1) f α 1

  

(6)

(7)

− uα uβ uβ ∂x(1) f α 2 





− f2 uβ uβ ∂x(1) uα − 2f2 uα uβ ∂x(1) uβ   α    α 

(8)

(9)

+ ρεI ∂x(1) uα α 



(11)



(10)

+ uα εI ∂x(1) ρ α 



(6.4.19)



(12)

which should yield (1)

ρ∂t εI = − ρuα ∂x(1) εI − ρεI ∂x(1) uα   α   α  (13)

(6.4.20)

(14)

Terms (1) and (10), (4) and (9), and (5) and (8) cancel each other when f2 = ρ/2. The sum (2) + (3) + (7) + (12) gives (13) when f1 = 2ρεI . Finally, (6) + (11) gives (14). Thus we obtain the constraints 8 and 9: 3 B1 + 24B2 = 2ρεI (= f1 ) 2

(6.4.21)

6.4 Thermal LBM: derivation of the coefficients

9 3 ρ E1 + 72E2 + G1 + 24G2 = (= f2 ) 8 2 2

263

(6.4.22)

Calculation of the coefficients. B1 and B2 are constrained by (6.4.3) and (6.4.21) 3B1 + 12B2 = ρ 3 B1 + 24B2 = 2ρεI 2 which lead to the unique solution 4 ρ(1 − εI ) 9 ρ (4εI − 1). B2 = 36 B1 =

(6.4.23) (6.4.24)

A1 and A2 are constrained only by (6.4.2) and (6.4.10) and A0 only by (6.4.2): A0 + 6A1 + 6A2 = ρ 3A1 + 12A2 = p (= ρεI ). In order to obtain a unique solution one may require (compare a similar constraint in Section 5.4, Eq. 5.4.4) that A1 /B1 = A2 /B2 . This leads to 5 A0 = ρ(1 − εI + 2ε2I ), 2

4 A1 = ρ (εI − ε2I ), 9

A2 = ρ

1 (−εI + 4ε2I ), 36

e.g. the expressions given by Alexander et al. (1993). The Cν and Dν are constrained by (6.4.1), (6.4.11), and (6.4.12): 3C1 + 12C2 + D0 + 6D1 + 6D2 = 0 3 C1 + 24C2 = ρ 2 and

6D1 + 24D2 = −ρ.

The choice of Alexander et al. is consistent with these constraints. Eν and Gν are constrained by (6.4.4) and (6.4.22): 9 E1 + 36E2 + 3G1 + 12G2 = 0 4 3 ρ 9 E1 + 72E2 + G1 + 24G2 = . 8 2 2 Alexander et al. (1993) have chosen G1 = 0 = G2 . Then E1 and E2 are uniquely given by 4 (6.4.25) E1 = − ρ 27 ρ (6.4.26) E2 = 108

264

6. Appendix

6.5 Schl¨ afli symbols Regular polytopes can be characterized by Schl¨ afli2 symbols instead of listing, for example, the coordinates of the whole set of vertices. Coxeter (1963, p.126/7) defines a polytope “as a finite convex3 region of n-dimensional space enclosed by a finite number of hyperplanes”. A polytope is characterized by its ensemble of vertices. Two-dimensional polytopes are called polygons. Threedimensional polytopes are called polyhedra. The part of the polytope that lies in one of the hyperplanes is called a cell (each cell is a (n − 1)-dimensional polytope; example: consider the cube where the cells are squares). The cells of polyhedra are called faces; they are polygons bounded by edges or sides. Edges join nearest-neighbor vertices. Thus a four-dimensional polytope Π4 has solid cells Π3 , plane faces Π2 (separating two cells), edges Π1 , and vertices Π0 . A polygon with p vertices is said to be regular if it is both equilateral (all sides are equal) and equiangular (all angles between nearest neighbor vertices are equal). If p > 3 a polygons can be equilateral without being equiangular (a rhomb, for example), or vice versa (a rectangle). Regular polygons are denoted by {p} (the Schl¨ afli symbol = number of vertices put in cranked brackets); thus {3} is an equilateral triangle, {4} is a square, {5} is a regular pentagon, and so on. A polyhedron is said to be regular if its faces are regular and equal, while its vertices are all surrounded alike. If its faces are {p}’s (i.e. regular polygons), q surrounding each vertex, the polyhedron is denoted by the Schl¨ afli symbol {p, q}. In three dimensions there exist only five regular polyhedra, namely the Platonic solids (compare Section 3.4). Consider, for instance, the cube. The faces are squares (4 edges) and each vertex is surrounded by 3 faces. Accordingly the cube is denoted by the Schl¨ afli symbol {4, 3}. The Schl¨ afli symbols for the other Platonic solids read: tetrahedron {3, 3}, octahedron {3, 4}, dodecahedron {5, 3}, icosahedron {3, 5}. Please note that the dual polytope to {p, q} has the Schl¨ afli symbol {q, p}. A polytope Πn (n > 2) is said to be regular if its cells are regular and there is a regular vertex figure4 at every vertex (‘regular surrounded’). It can be shown that as a consequence of this definition all cells are equal ({p, q} for n = 4) and the vertex figures are all equal ({q, r} for n = 4). A regular polytope Π4 is denoted by the Schl¨ afli symbol {p, q, r} where the cells are {p, q} and r is the number of cells that surround an edge. The three regular polytopes in four dimensions 2 3 4

After the Swiss mathematician Ludwig Schl¨ afli (1814-95). “A region is said to be convex if it contains the whole of the segment joining every pair of its points.” (Coxeter, 1963, p.126) “If the mid-points of all the edges that emanate from a given vertex O of Πn lie in one hyperplane ..., then these mid-points are the vertices of an (n − 1)dimensional polytope called the vertex figure of Πn at O.” (Coxeter, 1963, p. 128)

6.5 Schl¨ afli symbols

{3, 3, 3},

{3, 3, 4},

{4, 3, 3},

265

(6.5.1)

are bounded by tetrahedrons ({3, 3}) or cubes ({4, 3}). Π4 = {4, 3, 3} is the hypercube (not FCHC!). Similarly, a regular polytope Π5 whose cells are {p, q, r} must have vertex figures {q, r, s}, and thus will be denoted by Π5 = {p, q, r, s}.

(6.5.2)

It can be shown (Coxeter, 1963) that the parameters of the Schl¨ afli for regular polyhedra are constrained by “Schl¨ afli’s criterion” which reads 1 1 1 + > p q 2 and sin

for {p, q}

π π π sin > sin p r q

(6.5.3)

for {p, q, r}.

(6.5.4)

This and p, q, r ≥ 3 leads to {3, 3},

{3, 4},

{4, 3},

{3, 5},

{5, 3}

(6.5.5)

and {3, 3, 3},

{3, 3, 4},

{4, 3, 3},

{3, 4, 3},

{3, 3, 5},

{5, 3, 3}.

(6.5.6)

{3, 4, 3} is the face-centered hypercube (FCHC). Since Schl¨ afli’s criterion is merely a necessary condition, it remains to be proved that the corresponding polytopes actually exist (this can be very laboriously!; result: all above mentioned polytopes exist). Further reading: Rothman and Zaleski (1997, p. 265-270).

266

6. Appendix

6.6 Notation, symbols and abbreviations General remarks: Latin indices refer to the lattice vectors and run from 0 or 1 to l where l is the number of non-vanishing lattice velocities. Greek indices assign the cartesian components of vectors and therefore run from 1 to D where D is the dimension. If not otherwise stated (Einstein’s) summation convention is used, i.e. sum l mation is performed over repeated indices (ni ci = i=1 ni ci ). No summation l will be done over primed indices (nj
vj
= j
=1 nj
vj
). Table 6.6.1. Notation (miscellaneous symbols) Symbol

Meaning

∇ ∂ & | ∧ ∼ ∪ ∩ ¯ ◦

nabla operator partial derivative AND (Boolean operator) OR (inclusive or; Boolean operator) XOR (exclusive or; Boolean operator) NOT (Boolean operator) union (Boolean algebra) intersection (Boolean algebra) complementation (Boolean algebra) composition (of two elements; group theory)

6.6 Notation, symbols and abbreviations Table 6.6.2. Notation (Latin letters) Symbol

Meaning

A A(s → s ) Ass
Aν , Bν ... (t) ai B B b C C ci ciα cs cv D DkQb d E E EA F F F Fm f f0 G Gα1 α2 ...αn G(ρ) g(ρ) g H Hi h I I(P ) J(f ) J(a, b) j Kn k kB ks

Lagrange multiplier transition probability transition matrix free parameters of equilibrium distributions state of cell i at time t Boolean algebra Lagrange multiplier number of lattice velocities (’bits’) collision operator number of corners lattice vectors, lattice velocities cartesian component of the ith lattice velocity speed of sound heat capacity at constant volume (spatial) dimension lattice notation (k=dimension, b=number of lattice velocities) density per cell evolution operator number of edges Ekman number operator that interchanges particles and holes body force number of faces distribution functions Coriolis parameter Coriolis parameter at ϕ0 isometric group generalized lattice tensor of rank n g-factor (breaking Galilean invariance) = G(ρ)/2; g-factor (breaking Galilean invariance) group element Boltzmann’s H (= - entropy) Zanetti invarinats Lagrange multiplier identity operator measure of information collision operator (BGK) Jacobi operator momentum density Knudsen number number of states (CA) Boltzmann constant friction coefficient

267

268

6. Appendix

Table 6.6.3. Notation (Latin letters; continued) Symbol

Meaning

L L Lα1 α2 ...αn Ma IN IN0 Ni ni O(2 ) Oαβ P P Pαβ p pi Q Q(f, f ) Qiαβ q IR Re Re,g Ro r rj S S T (MA) Tαβγδ Tx,y t U u = (u, v, w) WM Wi wi x = (x, y, z) ZZk Z

a lattice characteristic length scale lattice tensor of rank n Mach number the set of natural numbers (integers) the set of non-negative integers mean occupation number (real variable) occupation number (Boolean variable) on the order of 2 orthogonal transformation matrix discrete probability distribution kinematic pressure momentum flux tensor pressure probabilities set of possible automata states collision integral Q-tensors (FHP) (vectorial) Lagrange multiplier the set of real numbers Reynolds number grid Reynolds number Rossby number range (CA) cartesian coordinates of nodes streaming (propagation) operator entropy temperature momentum advection tensor (MAT) components of the wind stress time characteristic speed velocity Munk scale global equilibrium distributions weights (generalized lattice tensors) cartesian coordinates residue class (integers modulo k) discrete set

6.6 Notation, symbols and abbreviations

Table 6.6.4. Notation (Greek letters) Symbol

Meaning

β Γ ∆1 ∆t ∆x δαβ δαβγδ  αβγ εI εK θ κ κ λ ν ξ ξ ρ σ τ ψ ψn Ω Ω Ωi ω

gradient of the Coriolis parameter phase space collision function time step spatial step size Kronecker symbol generalized Kronecker symbol expansion parameter Levi-Civita symbol internal energy density kinetic energy density = kB T temperature in energy units diffusion coefficient magnitude of wave number mean free path shear viscosity bulk viscosity random Boolean variable mass density collision cross section collision time stream function collision invariants set of events angular velocity of the Earth collision operator SOR or viscosity parameter

269

270

6. Appendix

Table 6.6.5. Abbreviations Acronym

Meaning

BC BGK CA EFD FCHC FHP HPP LBM LBGK LGCA MD ODE PCLBM PDE PI MSC NSE q.e.d. SOR

boundary conditions Bhatnagar, Gross, Krook cellular automata explicit finite difference face-centered hypercube Frisch, Hasslacher, Pomeau Hardy, Pazzis, Pomeau lattice Boltzmann model lattice BGK models lattice-gas cellular automata molecular dynamics ordinary differential equation pressure corrected lattice Boltzmann model partial differential equation pair interaction multi-spin coding Navier-Stokes equation quot erat demonstrandum successive over-relaxation