On the Lattice Boltzmann Method: Implementation and Applications
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classified information.
Kang Jin
Certificate of Approval:
Paul G. Schmidt
Professor
Mathematics and Statistics
Amnon J. Meir, Chair
Professor
Mathematics and Statistics
Wenxian Shen
Professor
Mathematics and Statistics
Jay Khodadadi
Professor
Mechanical Engineering
George T. Flowers
Interim Dean
Graduate School
On the Lattice Boltzmann Method: Implementation and Applications
Kang Jin
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
December 19, 2008
On the Lattice Boltzmann Method: Implementation and Applications
Kang Jin
Permission is granted to Auburn University to make copies of this dissertation
at its discretion, upon the request of individuals or institutions and
at their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Kang Jin was born in Shanghai, China in 1979. He did all his undergraduate study
in Shanghai. He graduated from East China Normal University in 2001 with a Bachelors
Degree in Mathematics. He then went to the United States in 2002. He accepted an
o?er from Auburn University, and began his graduate study as well as being a Teaching
Assistant. He got the Master of Science in Mathematics in 2005, and Then continue his
study as a Ph.D student in Auburn University.
iv
Dissertation Abstract
On the Lattice Boltzmann Method: Implementation and Applications
Kang Jin
Doctor of Philosophy, December 19, 2008
(B.S., East China Normal University, 2001)
125 Typed Pages
Directed by Amnon J. Meir
We studied the development and di?erent types of the Lattice Boltzmann method.
We gave several implementations. Then we presented two moving boundary treatments
for the Lattice Boltzmann method, the second one is new. We also gave an incompresi
bility enhancement for the Lattice Boltzmann method in order to better simulate some
problems using the moving boundary. Finally we gave a MHD solution using the Lattice
Boltzmann method.
v
Acknowledgments
I thank Auburn University Mathematics and Statistics department for o?ering me
the Graduate Teaching Assistant position. This is very important to me. Without this
o?er I could not be where I am today. I also thank the faculty in Math department
who give me lots of help. I thank my committee members Dr. Wenxian Shen, Dr. Paul
Schimdt and Dr. Jay Khodadadi. Special thank to my advisor, Dr. Amnon J. Meir. His
support and advice are the key to my work.
vi
Style manual or journal used Journal of Approximation Theory (together with the
style known as ?aums?). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package T
E
X (specifically
L
A
T
E
X) together with the departmental stylefile aums.sty.
vii
Table of Contents
List of Figures x
1Introduction 1
1.1 Background................................... 1
1.2 TheFHPmethod................................ 3
1.2.1 Basicmodel............................... 3
1.2.2 Macroscopicquantities......................... 7
1.3 TheLatticeBGKmethod........................... 7
1.3.1 TheD2Q7method........................... 7
1.3.2 Other frequently used Lattice Boltzmann Methods . . . . . . . . . 22
1.4 Applications................................... 25
2 Moving Boundary Problems 26
2.1 Background................................... 26
2.2 Moving boundary for LBM . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1 Moving boundary without mass conservation . . . . . . . . . . . . 30
2.2.2 Moving boundary with mass conservation . . . . . . . . . . . . . . 34
2.3 Sharp boundary definitiondetails....................... 39
2.3.1 Line................................... 40
2.3.2 Arc ................................... 41
2.4 Implementation................................. 41
2.5 ConclusionandFutureWork ......................... 42
3 Incompressibility 45
3.1 Background................................... 45
3.2 AnincompressibilityenhancedschemeforLBGK.............. 47
3.2.1 Thescheme............................... 47
3.2.2 Anexample............................... 48
3.3 Incompressibility with moving boundary . . . . . . . . . . . . . . . . . . . 48
3.4 Conclusion ................................... 54
4MHDWithConstantB 55
4.1 Backgound ................................... 55
4.2 Implementations ................................ 56
4.2.1 ExampleI................................ 56
4.2.2 ExampleII ............................... 59
4.3 Conclusion ................................... 61
viii
5Conclusion 62
Bibliography 64
Appendices 66
A FHP collision lookup table 67
B Partial Matlab code I: d2q9_rotating_polygon.m 69
C Partial matlab code II: d2q9_rotating_polygon_collision.m 82
ix
List of Figures
1.1 The hexagonal grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 FHPcollisionrule ............................... 6
1.3 Thed2q7_lattice. ............................... 8
1.4 Left boundary of a hexagonal grid. . . . . . . . . . . . . . . . . . . . . . . 15
1.5 DrivencavityatRe =10. ........................... 17
1.6 DrivencavityatRe = 100. .......................... 18
1.7 DrivencavityatRe = 200. .......................... 18
1.8 DrivencavityatRe = 400. .......................... 19
1.9 GhiaGhiaShin?sresult............................. 19
1.10 Driven cavity at Re = 800. .......................... 20
1.11 Dimensionless xvelocity profileatthegeometrycenter. .......... 20
1.12 Dimensionless yvelocity profileatthegeometrycenter........... 21
1.13 A uniform flowpastacylinder ........................ 23
1.14TheD2Q9lattice. ............................... 24
2.1 Noslip Boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2 Slip Boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.3 Above: boundary on the halfway; Below: boundary on the node . . . . . 29
2.4 Cases when q>1/2 or q<1/2. . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 Thechangeofstateillustration......................... 33
2.6 A halfway bounceback boundary. . . . . . . . . . . . . . . . . . . . . . . 35
x
2.7 Bounceback boundary when q<1/2. .................... 37
2.8 Bounceback boundary when q>1/2. .................... 38
2.9 An example of a boundary node. . . . . . . . . . . . . . . . . . . . . . . . 39
2.10 The change of state in this bounceback based scheme. . . . . . . . . . . . 40
2.11 The definition of a straight line boundary. . . . . . . . . . . . . . . . . . . 41
2.12 The definition of an arc boundary. . . . . . . . . . . . . . . . . . . . . . . 42
2.13Therotatingtriangleexample. ........................ 43
3.1 AtestusingtheincompressibleLatticeBGKmethod............ 49
3.2 The velocity contour of the damper with the standard Lattice BGK method. 50
3.3 The velocity contour of the damper with the incompressible Lattice BGK
method...................................... 51
3.4 The density of the fluid in the damper with the standard Lattice BGK
method...................................... 52
3.5 The density of the fluid in the damper with the incompressible Lattice
BGKmethod................................... 53
4.1 The boundary conditions of ?. ........................ 57
4.2 TheresultofthecubicboxMHDtest..................... 58
4.3 Theexampleofamagneticpump. ...................... 59
4.4 Thevelocitycontourofthemagneticpumpexample............. 60
xi
Chapter 1
Introduction
1.1 Background
The Lattice Gas model and Lattice Boltzmann method are methods for simulating
fluid flows. The flow of incompressible fliuds can be described by the NavierStokes
equation
?nullu
?t
+(nullu??)nullu = ??P + ??
2
nullu (1.1)
and the continuity equation
??nullu =0 (1.2)
where nullu is the velocity, P = p/?
0
the kinematic pressure, p the pressure, ?
0
the constant
mass density and ? = ?/?
0
the kinematic shear viscosity, and ? the dynamic shear
viscosity.
The Lattice Boltzmann method is derived from the Lattice Gas method. These
two methods are di?erent from other methods such as finite di?erence, finite volume,
and finite elements which are based on the discretization of partial di?erential equations
(topdown models [1]). These two method are based on a discrete microscopic model
which conserves desired quantities (such as mass and momentum), the partial di?erential
equations can then be derived by multiscale analysis (bottomup models).
First introduced in 1973, by Hardy, de Pazzis and Pomeau (HPP) [2], the HPP
method is a Lattice Gas method. It simulates the microscopic behavior of the fluid
by utilizing a square grid. The basic idea was to create a simple Cellular Automaton
1
obeying nothing but conservation laws at a microscopic level that was able to reproduce
the complexity of real fluid flows. Fluid particles of identical mass are only allowed to
travel on the lattice at unit speed. All lattice sites, which are the intersections of the
lattice, are exclusive. This means that only one particle is allowed to travel at each of
the four directions of one site. This gives a maximum of 4 particles at each site. Each
site can only take a finite number of states, actually 2
4
=16states. At each time step,
a collision occurs at each site, according to a collision rule which conserves the density
and the momentum. Then particles travel in straight lines (free streaming) until they
meet some other particle or the boundary.
This method is computer friendly, since only a 4bit variable is needed, and only
logical operations are repuired. Also, only the information from the four neighbours
are needed at each collision and streaming, this method is easy to parallelize. Simple
calculations as the HPP required, however, it leads to a macroscopical anisotropical
NavierStokes equation. This defect prevents the HPP from being widely used to model
fluid problems. In 1986 Frisch, Hasslacher and Pomeau (FHP) [3] introduced a lattice
gas method based on a hexagonal grid. This grid change made the FHP method exhibit
isotropy. Details of the FHP are discussed in the next section.
Similar to the HPP, logical operations made the FHP method easy to implement on
computers. Now the biggest problem of the cellular automata is noise, since it is based
on a FermiDirac distribution of the equilibrium population because of the exclusion
principle. The FermiDirac distruibution is a distribution that applies to particles called
fermions. Fermions have halfintegral values of the quantum mechanical property called
spin and are ?antisocial? in the sense that two fermions cannot exist in the same state.
Protons, neutrons, electrons, and many other elementary particles are fermions. The
2
results of the FHP is very noisy. Ensemble average and space average should both be
used. This may result in a grid size thousand times larger than the orighnal problem.
For example, if the final solution on a 100 ? 100 grid is needed, averaging on 10 ? 10
cells and ensemble averaging on 10 experiments, then the size of the calculation is 10
times the size that is on a 1000 ? 1000 grid! The lack of Galilean invariance is another
bigproblemoftheFHP.Thecollisionrulescanbewritteninalookuptable.Forthe
FHP method, this table should have a size of 2
7
?7. For multidimensional simulations,
huge lookup table associated with the collision rule made this almost impossible.
The Lattice Boltzmann method overcomes these defects very well. Instead of using
boolean variables at each site, the Lattice Boltzmann method uses real numbers. The
first method, proposed by McNamara and Zanetti [4], replaced the boolean variables
with their ensemble average. The statistics noise is greatly reduced. After that the
linear collision operator [5] came into being and then the enhanced collision rule [6].
The breakthrough was the single relaxation time approximation, known as the Lattice
BGK method, named after Bhatnagar, Gross, and Krook [7]. This method dramatically
reduced computation and gives pretty good results in various applications. In this article,
we will discuss in detail the Lattice BGK method. Simulations up to Reynolds number
1000 are presented.
1.2 The FHP method
1.2.1 Basic model
The Lattice Gas Cellular Automata simulates molecular collision in a discretized
fashion. Consider a hexagonal grid shown in Figure 1.1. Each site is surrounded by 6
3
Figure 1.1: The hexagonal grid. One site can contain a maximum of 7 particles.
neighbours, connected by unit vectors
nulle
i
=(sin(
?
2
(i?1)),cos(
?
2
(i?1))),i=1,...,6. (1.3)
Exclusion principles allow a maximum of 7 particles at one site, one moving particle
in each of the 6 directions, together with a rest particle in the middle. Here we use
an occupation boolean variable n
i
(nullx,t),i=0,1,...,6 (0 stands for the rest particle) to
indicate particle presence or nonpresence at the i
th
direction of site nullx at time t.Thusa
7bit variable is enough to carry the information at one site. All particles have the same
mass and the same speed.
One time step consists of a collision and a streaming. Collision only occurs at
the sites, while streaming takes place on the connection between each two sites. the
collision rules should conserve mass and momentum. Figure 1.2 shows the basic set of
collision rules [8]. The left column are the instates and the right are the outstates. In
state means that particles are moving towards the center of the site. After the collision
4
follows the outstate, particles then move away from the center of the site and begin
streaming. So a full time step is:
instate =? collision =? outstate =? streaming =? instate.
By rotating, flipping, and combining these 7 rules, one can get a full set of 128 collision
rules. Notice that some instates will lead to two equivalent outstates. It is not necessary
to pick an outstate randomly at every site. Notice that picking a random number is
very time consuming. Instead, one can pick a single random boolean variable for all the
sites at one time step.
The 6direction discretization makes this method lack degrees of freedom on speed
directions, yet it can display all the complexities of fliud phenomema [9]. This is the
simplest isotropic model. By limiting the types of collision, the FHP can be divided into
three types. The FHPI only allows collisions of type (a) and (b). No rest particle is
present. Types of collision other than (a) and (b) are replaced by simple streamig as if
the particles don?t collide at all. FHPII adds the rest particle, together with collision
type (c), (d), and (e). FHPIII includes all types of collisions.
A noslip boundary condition can be imposed by a bounce back scheme, which
means particles that hit the boundary at any angle should move back (bounce back) in
the opposite direction. Reflection will lead to a slip boundary. The Dirichelet boundary
condition can be imposed as a random variable on the boundary with a probability
distuibution indicating the value at the boundary, then applying the collision followed
by the bounce back scheme.
5
Figure 1.2: FHP collision rule
6
1.2.2 Macroscopic quantities
Noise is the biggest problem of the FHP method. Hence both space average and
ensemble average should be used. The space average is achieved by averaging on small
cells, for example, 16 ? 16 cells. The density is given by ? =
P
i
n
i
.And?
0
is the mean
density, which is the average on the whole grid.
Here are some methoddependent quantities derived by Frisch et al. [3] for the three
types of FHP methods.
FHPI FHPII FHPIII
d ?
0
/6 ?
0
/7 ?
0
/7
c
s
1
?
2
q
3
7
q
3
7
?
1
12
1
d(1?d)
3
?
1
8
1
28
1
d(1?d)
3
1
1?4d/7
?
1
8
1
28
1
d(1?d)
1
1?8d(1?d)/7
?
1
8
where d is the mean density per link, c
s
is the speed of sound, ? is the kinematic viscosity.
1.3 The Lattice BGK method
1.3.1 The D2Q7 method
The scheme
Consider again the hexagonal lattice showed in figure 1.3. This time, the occupation
number is replaced by its ensemble average value, or, the particle distribution function
f
i
(nullx,t). The meaning of this function is the probability of finding a particle moving in
the i
th
direction of the site nullx at time t. The collision rules in FHP are replaced with
a collision operator ?
i
, and the particle distribution function should satisfy the Lattice
7
Figure 1.3: The d2q7_lattice.
Boltzmann equation
f
i
(nullx + e
i
,t+1)?f
i
(nullx,t)=?
i
. (1.4)
This collision operator has lots of forms. Here we talk about the simplest one, the
BGK single relaxation time model. Introduce the single relaxation parameter ?,andthe
equilibrium distribution function f
eq
i
(nullx,t). By assuming that the particle distribution
function approaches the equilibrium state at a constant rate, we should get
?
i
= ?
1
?
(f
i
?f
eq
i
). (1.5)
This gives us the equation
f
i
(nullx + e
i
,t+1)=(1?w)f
i
(nullx,t)+wf
eq
i
(nullx,t) (1.6)
where the weight w =
1
?
. The equilibrium distribution function has the form
f
eq
i
(nullx,t)=?(nullx,t)
?
1?z
6
+
D
6c
2
(nulle
i
?nullu)+
D(D +2)
12c
4
(nulle
i
?nullu)
2
?
Dnullu
2
12c
2
?
,i=1,...,6 (1.7)
f
eq
0
(nullx,t)=?(nullx,t)(z?
u
2
c
2
) (1.8)
8
where ?(nullx,t)=
P
i
f
i
is the density. Here z is a parameter, we choose z =
1
2
.AlsoD is
the dimension, here D =2. c = e
i
,herec =1. And the speed of sound c
s
is controlled
by the parameter z by
c
s
=
r
1?z
2
. (1.9)
The kinematic viscosity can be adjusted by choosing a proper relaxation parameter ?,
and the relation is
? =
c
2
D +2
?
? ?
1
2
?
. (1.10)
Recovering the NavierStokes Equations
The conservation laws From the definition of the unit vectors e
i
,onecangetthe
following equations [11]
X
i
e
i?
=0,(1.)
X
i
e
i?
e
i?
=
c
2
b
D
?
??
2
X
i
e
i?
e
i?
e
i?
=0 3
X
i
e
i?
e
i?
e
i?
e
i?
=
c
4
b
D(D +2)
(?
??
?
??
+ ?
??
?
??
+ ?
??
?
??
),(1.4)
and
X
i
e
i?
e
i?
e
i?
e
i?
e
inull
=0,
where e
i?
stands for the ? direction component (one of the i,j directions on the 2
dimentional plane) of the unit vector nulle
i
.Usingthefirst two, one can obtain the moments
9
of the equilibrium distribution function. First sum the equilibrium distribution function
and get conservation of mass and momentum
X
i
f
eq
i
= ?,(1.5)
and
X
i
f
eq
i
e
i?
= ?u
?
. (1.16)
Also, from the rest of the equations, one gets
X
i
f
eq
i
e
i?
e
i?
=
?(1?z)c
2
D
?
??
+ ?u
?
u
?
,(1.7)
and
X
i
f
eq
i
e
i?
e
i?
e
i?
=
?c
2
D +2
(u
?
?
??
+ u
?
?
??
+ u
?
?
??
). (1.18)
The ChapmanEnskog expansion The distribution funcition can be expanded as
follows
f
i
= f
(0)
i
+ nullf
(1)
i
+ null
2
f
(2)
i
+ ... (1.19)
where null?1. Here we can use the Knudsen number Kn as null. The Knudsen number is
defined as
Kn =
?
L
where ? isthemeanfreepath,andL is the characteristic length. One can think this
expension of the distuibution function f as an equilibrium distribution function f
(0)
together with some pertubations f
(1)
,f
(2)
,..., of higher order in null.Wealsoexpandnullx and
10
t as
nullx =
nullx
1
null
+ ..., t =
t
1
null
+
t
2
null
2
+ ... (1.20)
where nullx
1
= o(null),t
1
= o(null),t
2
= o(null
2
). In this case, one can get
?
?x
?
= null
?
?x
1?
+ ..., (1.21)
and
?
?t
= null
?
?t
1
+ null
2
?
?t
2
+ ... .
Now we perform a Taylor expension on the Lattice Boltzmann equation (1.4) in both
space and time, we obtain
"
?
?
?t
+ e
i?
?
?x
?
?
+
1
2
?
?
?t
+ e
i?
?
?x
?
?
2
#
f
i
(nullx,t)=?
i
. (1.22)
Notice that Einstein summation is used. So e
i?
?
?x
?
=
P
?=1,2
e
i?
?
?x
?
. Using the expensions
of f,
?
?x
?
,
?
?t
, together with equation (1.5), we get
"
?
null
?
?t
1
+ null
2
?
?t
2
+ nulle
i?
?
?x
1?
?
+
1
2
?
null
?
?t
1
+ null
2
?
?t
2
+ nulle
i?
?
?x
1?
?
2
#
(1.23)
?
?
f
(0)
i
+ nullf
(1)
i
+ null
2
f
(2)
i
?
= ?
1
?
(f
(0)
i
+ nullf
(1)
i
+ null
2
f
(2)
i
?f
eq
i
).
Set the 0
th
order approximation f
(0)
i
to be the equilibrium distribution function f
eq
i
.
The conservation of mass and momentum require that
P
i
f
(k)
i
=0and
P
i
f
(k)
i
e
i?
=0,
11
for k =1,2. From these equations to firstorder in null we get
?
?t
1
f
(0)
i
+ e
i?
?
?x
1?
f
(0)
i
= ?
f
(1)
i
?
. (1.24)
Summing over i and from equation (1.15) and (1.16) we get
??
?t
1
+
?
?x
1?
?u
?
=0. (1.25)
Now multiply equation (1.24) by the unit vectors e
i?
and again sum over i, using equation
(??)gives
?
?t
1
?u
?
+
?
?x
1?
?u
?
u
?
?
?
?x
1?
?(1?z)c
2
D
?
??
=0. (1.26)
From equation (1.23) to secondorder in null and by equation (1.24) we get
?
?
?t
2
+
1
2
?
?t
1
?
?
?t
1
+ e
i?
?
?x
1?
?
+
1
2
e
i?
?
?x
1?
?
?
?t
1
+ e
i?
?
?x
1?
??
f
(0)
i
+
?
?
?t
1
+ e
i?
?
?x
1?
?
f
(1)
i
= ?
1
?
f
(2)
i
. (1.27)
Summing over i and using equation (1.25), all
?
?t
1
+ e
i?
?
?x
1?
vanished, and one gets
?
?t
2
? =0.
Again multiplying the equation by a unit vector e
i?
gives the following
?
?
?t
2
e
i?
+
1
2
?
?t
1
?
?
?t
1
e
i?
+ e
i?
e
i?
?
?x
1?
?
+
1
2
e
i?
e
i?
?
?x
1?
?
?
?t
1
+ e
i?
?
?x
1?
??
f
(0)
i
+
?
?
?t
1
e
i?
+ e
i?
e
i?
?
?x
1?
?
f
(1)
i
= ?
1
?
f
(2)
i
. (1.28)
12
By multiplying equation (1.24) by e
i?
e
i?
?
?x
1?
, one can rewrite the second term of f
(1)
i
as
?
?x
1?
e
i?
e
i?
f
(1)
i
= ??
?
?
?t
1
?
?x
1?
e
i?
e
i?
f
(0)
i
+
?
?x
1?
?
?x
1?
e
i?
e
i?
e
i?
f
(0)
i
?
. (1.29)
Combining this term with the third term of f
(0)
i
, one gets
?
?
?t
2
e
i?
+
1
2
?
?t
1
?
?t
1
e
i?
+ e
i?
e
i?
?
?x
1?
?
?
? ?
1
2
??
?
?t
1
?
?x
1?
e
i?
e
i?
+
?
?x
1?
?
?x
1?
e
i?
e
i?
e
i?
??
f
(0)
i
+
?
?t
1
e
i?
f
(1)
i
= ?
1
?
f
(2)
i
. (1.30)
Summing over i, the righthand side is 0. The second term of f
(0)
i
is 0 by equation (1.26).
The term of f
(1)
i
is 0 by the conservation of momentum requirement. The third term
of f
(0)
i
can be obtained from equation (??) to the order O(u) andthenconvertingtime
derivatives into spatial derivatives using equation (1.25), we get
?
?t
2
?u
?
=
?
? ?
1
2
??
?
?t
1
?
?x
1?
?(1?z)c
2
D
?
??
+
?
?x
1?
?
?x
1?
?c
2
D +2
(u
?
?
??
+ u
?
?
??
+ u
?
?
??
)
?
=
?
? ?
1
2
??
?
?x
1?
?
?x
1?
?c
2
D +2
u
?
+
?
?x
1?
??
2c
2
D +2
?
(1?z)c
2
D
?
?
?x
1?
?u
?
??
. (1.31)
13
By setting the kinematic shear viscosity ? =
?
? ?
1
2
?
c
2
D +2
and the kinematic bulk
viscosity ? =
?
? ?
1
2
?
?
2c
2
D +2
?
(1?z)c
2
D
?
, the above equation becomes
?
?t
2
?u
?
= ?
?
2
?x
2
1?
?u
?
+
?
?x
1?
?
?
?
?x
1?
?u
?
?
. (1.32)
Using all these equations (provided above), one can show that the NavierStokes equa
tion, the momentum equation
??
?t
u
?
+
?
?x
?
?u
?
u
?
= ?
?
?x
?
?
?(1?z)c
2
D
?
??
?
+ ?
?
2
?x
2
?
?u
?
+
?
?x
?
?
?
?
?x
?
?u
?
?
(1.33)
and the continuity equation
??
?t
+
?
?x
?
?u
?
=0 (1.34)
are satisfied. For an incompressible flow, these two equations reduce to equation (1.1)
and (1.2).
Boundary and initial conditions
The bounce back scheme is still good for the noslip boundary condition. Bounce
back boundary conditions only give firstorder accuracy. A Dirichlet boundary condition
can be achieved by solving the system of equations at the boundary sites
f
1
+
1
2
f
2
?
1
2
f
3
?f
4
?
1
2
f
5
+
1
2
f
6
= ?u
x
, (1.35)
?
3
2
(f
2
+ f
3
?f
5
?f
6
)=?u
y
, (1.36)
14
Figure 1.4: Left boundary of a hexagonal grid.
wherenullu =(u
x
,u
y
) is the velocity vector. For example, in the simulation of a driven cavity
flow, assume the top boundary has the speed nullu =(u
0
,0).Thef
2
,f
3
(up directions) are
from the inside sites, and you can keep f
1
,f
4
(horizontal directions) as constants. Only
f
5
,f
6
(down directions) are unknowns. This is a linear system of equations involving
two unknowns. Notice that when nullu =(0,0), it?s a bounce back scheme.
The left boundary shown in Figure 1.4 is not smooth in a microscopic view because
of the hexagonal structure of the grid. But it is smooth enough in a macroscopic view.
One can take the macroscopic boundary as the average of this boundary.
For initial conditions, one may use the equilibrium distribution from the given val
ues of ? and nullu. Bad initial distribution, for example, far away from the equilibrium
distribution, will make the method unstable, and eventually lead to blow up.
15
Implementation
Consider a driven cavity flow again. This time, an array of 7 floatingpoint variable
is needed for the information at one site. So we create a T ? 7 matrix M,whereT is
the total number of sites. We number the sites the same way as in FHP. The program
structure is also pretty much the same as in FHP, except that in the LBGK method, the
collision and the streaming are combined together by the Lattice Boltzmann equation.
The equilibrium distribution function is a mapping from the given ? and nullu to the matrix
M.Sofirst we can use this to initialize M. In the collision, we calculate the ? and nullu by
?(nullx,t)=
X
i
M(n,i), (1.37)
u
x
(nullx,t)=
null
M(n)?nulle
ix
, (1.38)
and
u
y
(nullx,t)=
null
M(n)?nulle
iy
, (1.39)
where n is the number of the site corresponding to nullx,andnulle
i
=(sin(
?
2
(i?1)),cos(
?
2
(i?
1))), which is the link to the 6 neighbour sites. Then apply the weighted equilibrium
distribution function (equation 1.6) with a proper ?.
Results and data analysis
Driven Cavity Here we present a driven cavity example again. In figures 1.5, 1.6,
1.7, 1.8, and 1.10, we give the velocity vectors (left) and velocity contour (right) at the
steady state for Reynolds number 10, 100, 200, 400 and 800. We give the result of Ghia,
Ghia, and Shin [12] for Reynolds number at 400 for comparison (The computations were
16
Figure 1.5: Driven cavity at Re =10.
performed using the timemarching capabilities of WIND to approach the steadystate
flow starting from the freestream conditions). We also give the velocity profiles for u
and v through the geometirc center of the cavity. For comparison, refer to Shuling Hou
and Qisu Zou et al [13].
Flow past a cylinder
This is an example of a uniform flow past a cylinder. This example is done on a
360 ? 1000 grid. The speed of the uniform flow coming from the left is 0.5. And the
flow is free to flow out at the right end. The cylinder was placed at the center of the left
inlet with a diameter 120. The Reynolds number is 400. The top and bottom are noslip
boundaries. It is well know that at a Reynolds number greater than 100, the flow past a
cylinder will give a Von Karman vortex street. Here we give the figures of both velocity
17
Figure 1.6: Driven cavity at Re = 100.
Figure 1.7: Driven cavity at Re = 200.
18
Figure 1.8: Driven cavity at Re = 400.
Figure 1.9: Ghia Ghia Shin?s result. The plot of the velocity contour with a Reynolds
number of 400
19
Figure 1.10: Driven cavity at Re = 800.
Figure 1.11: Dimensionless xvelocity profile at the geometry center of the cavity for
Reynolds number 10, 100, 200, 400, 800. Dashed line is the result from Ghia Ghia Shin
at Reynolds number 400.
20
Figure 1.12: Dimensionless yvelocity profile at the geometry center of the cavity for
Reynolds number 10, 100, 200, 400, 800. Dashed line is the result from Ghia Ghia Shin
at Reynolds number 400.
21
contour and vorticity contour at time step 5000. Both figures show the back half of the
cylinder.
1.3.2 Other frequently used Lattice Boltzmann Methods
D2Q9 method
D2Q9methodisthemostcommonlyused2dimensionalLBM.Ithas3discretegrid
speed: 0, c and
?
2c. It has 9 directions: 8 moving directions, plus the rest particle. This
can be discribed by the link vectors that point to the neighbors, which read
e
i
=(0,0) i =0, (1.40)
e
i
=(?c,0),(0,?c),i=1,2,3,4, (1.41)
and
e
i
=(?c,?c),i=5,6,7,8. (1.42)
Figure 1.14 shows the D2Q9 Lattice. The equilibrium dsitribution is
f
eq
i
(x,t)=W
i
?(x,t)
?
1+3
e
i
?
??
u
c
2
+
9
2
(e
i
?
??
u )
2
c
4
?
3
2
??
u
2
c
2
?
,i=0,...,8, (1.43)
where
W
i
=
4
9
,i=0, (1.44)
W
i
=
1
9
,i=1,2,3,4, (1.45)
W
i
=
1
36
,i=5,6,7,8, (1.46)
22
Figure 1.13: An example of a uniform flow past a cylinder for Re = 400.Above:velocity
contour. Below: vorticity contour.
23
Figure 1.14: The D2Q9 lattice.
and the kinematic viscosity is given by
? =
c
2
3
?
? ?
1
2
?
. (1.47)
D3Q19 method
The D3Q19 is a 19direction 3dimensional LBM. The link vectors are
e
i
=(0,0,0) i =0, (1.48)
e
i
=(?c,0,0),(0,?c,0),(0,0,?c),i=1,...,6, (1.49)
and
e
i
=(?c,?c,0),(?c,0,?c),(0,?c,?c) i =7,...,18. (1.50)
Note that just like the D2Q9 LBM, D3Q19 also features 3 discrete speed: 0, c,and
?
2c,
so it shares the same form of equilibrium distribution with D2Q9, with a di?erent set of
24
W
i
s, which read
W
i
=
1
3
,i=0, (1.51)
W
i
=
1
18
,i=1,...,6, (1.52)
W
i
=
1
36
,i=7,...,18. (1.53)
1.4 Applications
In the next 3 Chapters, we will give several applications of the Lattice Boltzmann
method. In Chapter 2, we will talk about the moving boundary problems using the
Lattice Boltzmann method. We present two di?erent moving boundary treatments. In
Chapter 3, we give an incompressibility enhancement of the Lattice Boltzmann method.
In Chapter 4, we give an application that couples the Lattice Boltzmann method with a
Poisson?s equation to solve a MHD flow problem.
25
Chapter 2
Moving Boundary Problems
2.1 Background
Basic types of problems
The moving boundary problems are a type of CFD problems which have time de
pendent solid moving boundaries. Numerous physical phenomena involve solidfluid
interaction. This is in contrast to a static boundary where the boundary nodes are fixed
on a fixed mesh. These types of problems include deformation of structures, moving solid
objects, and boundaries that evolve in time. The main methods of solving these moving
boundary problems can be brieflyclassified into two major categories: the Lagrangian
methods, and the Eulerian methods.
On one hand, the Lagrangian methods maintain a mesh that follows the movement
of the solid. The solidfluid interface is explicitly and accurately tracked. And boundary
conditions can be applied at the exact location of the interface at each time step. How
ever, the regeneration of the mesh is time consuming. Also, in some cases, the resulting
grid may be unevenly distributed. This degrades the accuracy on the boundary.
The Eulerian methods, on the other hand, use a mesh that remains fixed. These
method can roughly be categorized into di?use type and sharp type. For the di?use
type, the exact location of the solidfluid interface is unknown. The methods accuracy
is to the order of the grid scale. The boundary conditions are applied in the sense of the
governing equations of the fluid and the solid. The obvious advantage is that there is
no need to regenerate the mesh at every time step. However, due to the unclarity of the
26
exact boundary location, when the interface is arbitrarily shaped, improved resolution
is di?cult to obtain.
Perhaps the more commonly used method is the immersed moving boundary method
(Peskin 1977). A fixedCartesiangridisused,andanexplicitsharpinterfaceisdefined.
This method is a mixture of Eulerian and Lagrangian framework. Clearly, this method
shares the virtures of both the Eulerian and the Lagrangian method: accurate boundary
condition can be applied without regenerating the mesh. At the same time, however,
several problems arise. One issue that needs to be dealt with is keeping track of those
nodes that go from the solid state to fluid state, or fluid state to solid state. This is
an issue that is not encountered in a pure Lagrangian or pure Eulerian method. A
steep change will result in a huge discontinuity on the boundary. One method is to
use a fractionalstep scheme (H. S. Udaykumar), which is based on a onedimensional
interpolation.
Basic boundary types for LBM
For Lattice Boltzmann method, methods that are close to the immersed moving
boundary can be used. This is due to the nature of the lattice boltzmanm method,
that is, it uses a fixed Cartisian grid, and it is easy to manage local behaviors. In the
next section, we will present two methods for treating moving boundaries in the LBM.
The first one is a widely used method which based on a quadratic interpolation. This
method is second order accurate, but at the cost of violating mass conservation. The
second method we will present is a newly developed method, which is closer in spirit
to the fractionalstep scheme. It conserves mass, and the change of state is gradual,
without a steep discontinuity.
27
Figure 2.1: Noslip Boundary
Before continuing, we describe several basic boundary types for LBM. The two most
common types of boundary conditions are: Dirichlet boundary condition and Newmann
boundary condition. In LBM, the simple realizations are of noslip boundary and slip
boundary (Figure 2.1, 2.2). The noslip boundary uses a bounceback scheme, without
considering the angle of incidence of the particle. This makes the noslip boundary very
easy to implement, since one doesn?t need to know the exact shape of the boundary.
The slip boundary is of reflection type. The angle of incidence should equal to the angle
of reflection. This makes the slip boundary much more di?cult and complicated to
implement.
The moving particles in LBM travell to the neighbor node in one time step. In order
to make sure that all particles are on the nodes, one could either set the boundary at
the nodes, or in the middle of two neighbor nodes. Figure 2.3 shows the two normally
used boundary positions. One of the di?culties of the moving boundary is to find a way
so that the boundary can be set in any intermediate position.
28
Figure 2.2: Slip Boundary
Figure 2.3: Above: boundary on the halfway; Below: boundary on the node
29
2.2 Moving boundary for LBM
2.2.1 Moving boundary without mass conservation
The scheme
A moving boundary treatment for LBM was proposed by P. Lallemand and L.S. Luo
[14]. This method is based on the simple bounceback boundary scheme and quadratic
interpolations, yielding a noslip boundary condition. As stated in the previous section,
the di?cultyliesindefining the boundary at a position other than the halfway between
nodes and on the nodes. For simplicity, we just study the case of one particle direction.
Figure 2.4 shows two cases when the boundary is not at a standard position. Here q
is the distance between the boundary and the closest node, assuming that the distance
between any two neighboring nodes is 1.
For the case q<1/2, the problem is how to define the left direction particle on node
r
j
, since it should came from the right direction particle traveling between r
j
and r
0
j
at
the previous time step. The idea is to implement a quadratic interpolation at position
x using the information on r
j
, r
0
j
,andr
00
j
. For the case q>1/2, position x is outside
of r
j
, r
0
j
,andr
00
j
. To avoid using extrapolation, one can do the streaming (propagation)
first, and then implement a quadratic interpolation on node r
j
using the information on
x, r
0
j
,andr
00
j
.
Also to be considered is the extra term F
i
introduced by the fluidsolid interac
tion. By considering the mass conservation
P
F
i
=0and the momentum conservation
P
??
e
i
F
i
= ?
??
u
w
, one can get F
i
= w
i
(
??
e
i
?
??
u
w
),where
??
u
w
isthespeedofthemoving
boundary, and w
i
is the weight of the mass on i direction.
30
Figure 2.4: Cases when q>1/2 or q<1/2.
31
The actual formulas are: For q<1/2
f
i
(r
j
,t)=q(1 + 2q)
b
f
i
(r
j
,t)+(1?4q
2
)
b
f
i
(r
0
j
,t)?q(1?2q)
b
f
i
(r
00
j
,t)+w
i
(
??
e
i
?
??
u
w
), (2.1)
And for q>1/2
f
i
(r
j
,t)=
1
q(2q +1)
b
f
i
(r
j
,t)+
2q?1
q
f
i
(r
0
j
,t)?
2q?1
2q +1
f
i
(r
00
j
,t)+
w
i
q(2q +1)
(
??
e
i
?
??
u
w
), (2.2)
where
b
f
i
is the distribution from the previous time step, and for D2Q9 LBM,
w
i
=
?
?
?
?
?
2/3,i=1,2,3,4
1/6,i=5,6,7,8
. (2.3)
Note in the case q>1/2, there is a correction of 1/[q(2q +1)][14]. This is obtained
by considering the analytic solution for the Couette flow. These two formulas coincide
when q =1/2, and reduce to a simple bounceback scheme.
Change of state treatment
For immersedmovingboundary problems with aCartisiangridand asharpinterface,
the change of state (a node goes either from fluid state to solid state, or from solid state
to fluid state) should be carefully considered. Figure 2.5 shows nodes that change state
as a sharp interface moves. The nodes that changed their state are marked with triangles.
The two di?erent states are marked with squares and circles.
32
Figure 2.5: The change of state illustration.
On the one hand, the change of a node from the fluid state to solid state is simplely
considered a nonissue, since the fluid motion is not calculated in solid. And the treat
ment is to set the distribution on those nodes to zero. On the other hand, the change of
a node from the solid state to fluid state is not trivial. Some method has to be used to
"create" the particle distribution of the LBM at those nodes.
Several methods could be used. A commonly used method is to use the information
from neighbors and extrapolate the distributions. One should use the direction which is
closest to
??
n, the outnormal vector of the moving interface. That is, to pick the direction
??
e
i
which maximizes the quantity
??
n ?
??
e
i
. Then, use a linear or quadratic extrapolation, or
some other method, to get the distributions on the newly changed state nodes. Another
method is to use the velocity of the moving interface and the average density of the
system, or a local average density, to get the distributions.
33
An example
2.2.2 Moving boundary with mass conservation
An obvious compromise of the quadratic interpolation moving boundary is that
the mass is not conserved. Perhaps the total mass conservation in the whole system
is not violated, but locally the problem could be serious in some cases. In Figure 2.5
consider the moving boundary as a slim bar. What the quadratic interpolation does
is, it continuously eleminate particles that are in front of the bar along the moving
direction, and creates particles behind it. This results in an e?ect that fluid particles
are continuously infiltrating through the bar, this reduces the presure di?erence between
the two sides of the bar. If this issue is critical to the whole system, then some other
method should be used. This is the motivation for us to develop a new method for LBM
that conserves mass.
The scheme
Let?s consider again a halfway bounceback boundary shown in Figure 2.6. And
first let?s only consider this in one direction. Then at the next time step, due to the
streaming, the following changes will be performed.
f
1
= f
0
1
,
f
0
2
= f
2
,
and
f
2
= f
1
.
34
Figure 2.6: A halfway bounceback boundary.
Inthecasewhenq<1/2, the boundary nodes no longer occupy a whole unit space,
but just a part of it. The first graph in Figure 2.7 shows this case. In order to explain
things, the second graph in Figure 2.7 reflects the left direction distribution to the right
side. This makes all the particles to be right direction particles. After the streaming
process, all particles move to the right by one unit distance, as indicated in the third
graph of Figure 2.7. Then the space occupied by f
1
, f
2
,andf
0
2
(second graph of figure
2.7) has been replaced by f
0
1
, f
1
,andf
2
(third graph of figure 2.7). By comparing the
space being occupied, one gets the transition formulas
f
1
=(1/2+q)f
0
1
, (2.4)
f
0
2
= f
2
+
(1/2?q)
(1/2+q)
f
1
, (2.5)
and
f
2
=
2q
(1/2+q)
f
1
+(1/2?q)f
0
1
. (2.6)
35
Similarly, as shown in Figure 2.8, in the case when q>1/2, the transition formulas
read
f
1
=
(q?1/2)
(q +1/2)
f
1
+ f
0
1
, (2.7)
f
0
2
=
f
2
q +1/2
, (2.8)
and
f
2
=
q?1/2
q +1/2
f
2
+
1
q +1/2
f
1
. (2.9)
In either of these two cases, the occupied space of the boundary nodes does not
equal to a unit space, thus the f
i
s can not represent the true distribution of a LBM
on that node. One can use the normalized distribution
e
f
i
= f
i
/(1/2+q) to get macro
quantities like density and velocity on the boundary nodes. For a D2Q9 method, one
needs to find the q for all directions. Figure 2.9 shows a typical example of a boundary
node. The occupied space on each direction is di?erent. But in most of the cases, for
opposite directions, they are the same. The collision process includes normalizing the
distribution, collision, then unnormalizing. In this process, the mass is not perfectly
conserved. However, if the speed of the moving boundary is much slower than the grid
speed, this imperfection can be neglect.
The change of state
In this scheme, the change of state transitionissmoothandrequiresnoextrawork.
As shown in Figure 2.10, the upper graph shows a situation when the boundary is very
close to a node, which is the case q<1/2. Although the occupied space is small on the
boundary node, the normalized distribution should not be too di?erent from the neighbor
36
Figure 2.7: Bounceback boundary when q<1/2.
37
Figure 2.8: Bounceback boundary when q>1/2.
38
Figure 2.9: An example of a boundary node.
node. Consider the case that after a time step, this node is coverd by the boundary, as
shown in the lower graph of Figure 2.10. Now the node is gone, but most of the volume
is still there. And the neighbor node becomes the boundary node with q>1/2.Sovery
naturally, it will take over this part that the previous boundary node left. Given that
the boundary moves slowly enough, the sum of the occupied space of the previous two
node will not di?er too much from the occupied space of the later boundary node. Thus
a smooth transition is achieved.
2.3 Sharp boundary definition details
The implementation of the moving boundary requires the knowledge of the exact
location of the sharp boundary. We demostrate the definition of two simple and widely
used boundaries: line and arc. Before locating the exact boundary, one needs to find the
39
Figure 2.10: The change of state in this bounceback based scheme.
boundary layer first. This is a layer that contains nodes that are closest to the boundary
in all directions.
2.3.1 Line
The definition of a line can be done by using a normal vector of this line. The
advantage of this is that, by taking the projection along the normal vector, it is very
easy to get the perpendicular distance to the line. Figure 2.11 shows a straight line
boundary, which lies between v and v0. n is the normal vector to the line, and l is the
shortest distance between the line and the center point c.Wewishtofind q, the distance
between v and the line along v0?v. Assume that v?v0 =1, q is given by
q =
v?n?l
(v?v
0
)?n
.
40
Figure 2.11: The definition of a straight line boundary.
2.3.2 Arc
The arc is of course defined by a partial circle. Figure 2.12 shows an arc boundary,
where r is the radius of the arc. Let e
i
= v
0
?v, the distance of v to the arc along e
i
is
given by
q = v?
e
i
e
i

?
s
r
2
?
?
?
?
?
v?(v?
e
i
e
i

)
e
i
e
i

?
?
?
?
2
.
2.4 Implementation
Here we present an example of this moving boundary scheme. This example is a
rotating triangle in a square domain full of fluid. All boundaries are noslip boundaries.
This example is done in a 300300 grid using D2Q9 Lattice Boltzmann method. Figure
41
Figure 2.12: The definition of an arc boundary.
2.13 shows the velocity contour with streamlines. This picture is taken after 2 full cycles
(720 degree).
2.5 Conclusion and Future Work
In this chapter we presented two moving boundary treatment for the Lattice Boltz
mann method. Both define sharp boundaries, and both require that the speed of the
moving boundary is much slower than the grid speed c.Thefirst uses quadratic interpo
lation. It is of second order accurate, but at the cost of voilating the mass conservation.
The second is of first order accurate, and conserves mass.
We would like to point out a future work of the second massconserved moving
boundary treatment. One draw back of this scheme is that it lacks the shear stress. The
reason is that it only consider the normal direction movement of the boundary, but not
42
Figure 2.13: The rotating triangle example.
43
the tangential direction movement. The first treatment is also the same, but since it
count the momentum transfer, the e?ectisnotsevere.Thismayexplainwhythepart
of the streamlines that connect to the boundary in figure 2.13 is almost perpendicular
to the boundary. So a future work of this moving boundary treatment is to consider the
shear stress that is given by the tangential direction movement of the boundary.
44
Chapter 3
Incompressibility
3.1 Background
In the previous chapters, we used Lattice Boltzmann method on incompressible
flows. However, we need to point out that an incompressible flow is an ideal flow.
It only exists in theorey. In practice, any flow is compressible to some extent. The
Lattice Boltzmann method, which is based on the Boltzmann equation for gas, simulates
compressible fluids with some finite speed of sound c
s
. When the fluid speed is su?ciently
small compared with c
s
, we should get a solution that converges to the incompressible
limit.
In recent years, several improvement were made to the Lattice Boltzmann method
in order to better approximate the incompressibility. This can be roughly categorized
into two groups: the improved singlerelaxation time model (Lattice BGK method), and
the multirelaxation time model [15]. The unique nature of the multirelaxation time
model makes it a better method for simulating incompressible flows. However, due to
the simplicity of implementation and being 30% more e?cient than the multirelaxation
time model, the Lattice BGK method is the favorite method for many. Unfortunately,
theoretically it is impossible to maintain a constant density in the Lattice BGK method.
In this chapter, we focus on the Lattice BGK method.
In the real world, for an incompressible fluid, the density ? = ?
0
+ ?? ,where?
0
is a constant and ?? is the density fluctuation, which should be of the order O(M
2
)
(M is the Mach number, and M ?? 0).X.HeandLSLuo[16]improvedtheLattice
45
BGK method by substituting ? = ?
0
+?? into the equilibrium distribution function f
eq
i
,
neglecting the terms of the order O(M
2
) and higher. For a D2Q9 Lattice BGK method,
the result is a new equilibrium distribution function, which reads
f
eq
i
(x,t)=w
i
?
? + ?
0
?
3
(
??
e
i
?
??
u )
c
2
+
9
2
(
??
e
i
?
??
u )
2
c
4
?
3
2
??
u
2
c
2
??
. (3.1)
This improved method can simulate both steady and unsteady flow problems.
Another improved Lattice BGK method is provided by Guo, Shi and Wang [17].
This scheme redefined the distribution when the fluidisatrest(
??
u =0). The new
equilibrium distribution function (again, for D2Q9) is
g
eq
i
= v
i
p
c
2
+ w
i
?
3
(
??
e
i
?
??
u )
c
2
+
9
2
(
??
e
i
?
??
u )
2
c
4
?
3
2
??
u
2
c
2
?
, (3.2)
and the macroscope velocity and pressure are given by
u =
8
X
i=1
c
??
e
i
g
i
(3.3)
and
p =
c
2
?v
i
"
8
X
i=1
g
i
?
3
2
??
u
2
c
2
#
. (3.4)
This scheme is an artificially incompressible Lattice BGK method due to the "neg
ative rest particle distribution". It also can simulate both steady and unsteady flow
problems.
46
3.2 An incompressibility enhanced scheme for LBGK
3.2.1 The scheme
In a traditional Lattice BGK method, pressure and density are proprotional to each
other. Theoretically, in an incompressible flow, pressure disturbances travel at infinite
speed. But in Lattice BGK method, pressure disturbances travel at the grid speed c.To
improve the incompressibility, it is reasonable to increase the speed at which pressure
disturbances travel. Practically, an incompressible flow is a flow in which the density
fluctuation ?? will dissipate in a very short time, or more precisely, the system will reach
a constant density in a very short time.
To achive a faster "dissipation", we implement an extra "propagation" after every
propagation step. Each node should "push out" or "pull in" the amount that equals to
the density fluctuation. The new distribution function then reads
b
f
i
(nullx + e
i
,t)=f
i
(nullx + e
i
,t)+
???
0
?
f
i
(nullx,t), (3.5)
where ? =
P
f
i
,andis?
0
the "supposed" constant density. In practice, it can be the
constant initial density. This is a redistribution of mass that is due to the density
fluctuation of the Lattice BGK method. It would be natural to redistribute the extra
mass to the neighbors. Since the redistributed amount is proportional to f
i
,thevelocity
is kept the same.
47
3.2.2 An example
We give an example using the new incompressible Lattice BGK method. Consider a
closed channel with a density fluctuation that is initially set as 1.05 at the left side, and
gradually decreases to 0.95 at the right side. We compare the results with a standard
Lattice BGK method. We were not supprised to see that our new proposed incom
pressible Lattice BGK method reached (almost) constant density much faster than the
standard Lattice BGK method due to the extra propagation.
3.3 Incompressibility with moving boundary
An interestingexample is obtained if we combine the mass conserving moving bound
ary with the incompressible Lattice BGK method. As a matter of fact, it makes more
sense now, since the continuity equation
??
?t
+??(?
??
u )=0
is satisfied. The incompressibility
D?
Dt
=
??
?t
+???
??
u =0
leads to
??
??
u =0.
We consider an example of a hydraulic rotation damper. A damper is a device that
restrains or depresses motion. It transfers the energy of the motion to kinetic energy
(heat) due to friction so that a sudden strong motion that might be dangerous to a
48
Figure 3.1: A test using the incompressible Lattice BGK method. From top to bottom:
1. velocity contour with normal Lattice BGK; 2. velocity contour with incompressible
Lattice BGK; 3. density plot with normal Lattice BGK; 4. density plot with incom
pressible Lattice BGK.
49
Figure 3.2: The velocity contour of the damper with the standard Lattice BGK method.
system is prohibitted. In this example, we consider a half cylinder system which holds
ahydraulicfluid. A devider is attached to the center of the half cylinder and allowed
to rotate. Between the divider and the cylinder wall is a gap that allows the fluid to go
from one side of the divider into the other. Because of the incompressibility of the fluid,
the rotation of the divider will force the fluid to go throught the small gap, which create
a counter force, thus damps the rotation.
50
Figure 3.3: The velocity contour of the damper with the incompressible Lattice BGK
method.
51
Figure 3.4: The density of the fluid in the damper with the standard Lattice BGK
method.
52
Figure 3.5: The density of the fluid in the damper with the incompressible Lattice BGK
method.
53
3.4 Conclusion
In this chapter we present an incompressibility enhancement for the Lattice Boltz
mann method. Since the pressure travel at an infinite speed, we perform an extra
propagation of the density disturbances. The result is that, it gives about the same ve
locity solution as the normal Lattice Boltzmann method. But in some extream situations
where the density disturbances are too strong, the incompressibility enhancement will
keep the particle distribution in a "safe" mode so that the Lattice Boltzmann method
can continue without encounter any problem.
54
Chapter 4
MHD With Constant B
4.1 Backgound
We present an application of the Lattice Boltzmann method to Magnetohydrody
namics (MHD). MHD is the theory of macroscopic interaction of electrically conducting
fluids and electromagnetic fields. Examples of such fluids include plasmas, liquid metals,
and salt water. The ideal MHD is that, in a moving conducting fluid, the magnetic field
will induce current, the interaction of the current and magnetic field creates a force in
this field and alters the velocity of the fluid, and also changes the magnetic field itself.
Assuming that the fluid is an incompressible viscous fluid, the governing equation of a
MHD system are the NavierStokes equation
nullu
t
+(nullu??)nullu = ??P + ??
2
nullu +
??
J ?
??
B +
??
F, (4.1)
and the Ohm?s law
??
J = ?(??? +
??
E +
??
u ?
??
B), (4.2)
together with the continuity equations
??
??
u =0 and ??
??
J =0, (4.3)
55
where nullu is the velocity,
??
J is the current,
??
B is the magnetic field,
??
E is the electric field,
P is the pressure, ? is the electric potential, ? is the kinematic viscosity, and ? is the
fluid conductivity.
We implement a LBM for the fluid. To solve Ohm?s law (4.2), we take the divergence
of both sides. Since the divergence of
??
J and
??
E are equal to zero (assuming there is no
net electric charge in the field), we get a poisson?s equation
4? = ?? (
??
u ?
??
B).
Inonetimestep,wesolvethepoisson?sequationfor? with a boundary condition that
yields a wellposed problem, then use ? to update
??
J ,thenupdate
??
u using LBM. For
simplicity, we neglect the induced magnetic field. So in this case,
??
B is a constant external
magnetic field.
4.2 Implementations
4.2.1 Example I
In this example, we consider a cubic box which contains electrically conducting
incompressible fluid. The boundary conditions of ? are shown in Figure 4.1. The little
arrows indicate dirichlet boundary, and neumann boundary elsewhere. The front side of
the dirichlet boundary is set to 1, and the back side is set to 0. The magnetic field
??
B is
set to 1 all over the cubic box, the direction is shown in the figure. The external electric
field
??
E is set to 0. We used D3Q19 for the fluid.
Figure 4.2 shows the result.
56
Figure 4.1: The boundary conditions of ?.
57
Figure 4.2: The result of the cubic box MHD test.
58
Figure 4.3: The example of a magnetic pump.
4.2.2 Example II
This is a magnetic pump example. In this example, we consider a loop which
contains electrically conducting incompressible fluid. The boundary conditions of the
electric potential ? are shown in Figure 4.3. The boundary of ? at the front shaded area
is set to 0, and at the back shaded area is set to 1. Other boundaries are all Neumann
boundaries. The magnetic field
??
B is set to 1 on z direction and 0 on x and y directions.
The external electric field
??
E is set to 0. The electic conductivity ? of the fluid is set to
1. We used D3Q19 for the fluid. Figure 4.4 gives the result at t = 500.
59
Figure 4.4: The velocity contour of the magnetic pump example.
60
4.3 Conclusion
In this chapter we coupled the Lattice Boltzmann method with the Ohm?s law. We
presented two examples. They showed that the 3dimension D3Q19 is a very fast solver
for the NavierStokes equation.
61
Chapter 5
Conclusion
In the last chapter, we give an overall conclusion of the Lattice Boltzmann method
and the applications.
The Lattice Boltzmann is a fast solver of the NavierStokes equation. Several dif
ferent discretizations are available. Take the 2dimension methods as an example: the
D2Q7 is the fastest 2dimension method due to the least number of discrete speed, while
the D2Q9 is a little slower but gives a little more accurate result. So one can pick the
right method to keep the best balance between performance and accuracy.
The biggest advantage of the Lattice Boltzmann method over other methods is that
it is very easy to handle complicated noslip boundaries. This is due to the fact that to
implement the bounce back boundary condition, one doesn?t need to know the angle of
incidence. So on programming using the Lattice Boltzmann method, one just need to
specify the fluid nodes and the solid nodes.
By using the moving boundary treatments we presented in Chapter 2, the Lattice
Boltzmann method can be used on even more complicated situations without compro
mising the performance. This is due to the fact that Lattice Boltzmann method uses
a fixed grid and an immersed moving boundary. We presented two moving boundary
treatments in Chapter 2: the second order accurate moving boundary, and the first order
massconserved moving boundary. As we discussed in the conclusion section of Chapter
2, one future work of the mass conserved moving boundary is to consider the shear stress.
62
The Lattice Boltzmann method is also easy to couple with other method to simulate
problems that are related to fluid. We gave examples of MHD flows in Chapter 4.
Overall, the Lattice Boltzmann method can be widely used to di?erent types of
problems, and it can do very well.
63
Bibliography
[1] Dieter A. WolfGladrow, LatticeGas Cellular Automata and Lattice Boltzmann
Models, An Introduction, SpringerVerlag Berlin Heidelberg 2000. ISSN 00758434
ISBN 3540669736
[2] J. Hardy, Y. Pomeau & O. de Pazzis, Time evolution of twodimensional model
system. I. Invariant states and time correlation functions, J. Math. Phys. 14 (1973),
pp. 17461759.
[3] U. Frisch, B. Hasslacher, and Y. Pomeau. Latticegas automata for the NavierStokes
equation. Physical Review Letters, 56:15051508, 1986.
[4] G.R.McNamaraandG.Zanetti.Use of the Boltzmann equation to simulate lattice
gas automata. Physical Review Letters, 61:23322335, 1988.
[5] F. J. Higuera and J. Lim?nez. Boltzmann approach to lattice gas simulations.Euro
physics Letters, 9 (7):663668, 1989.
[6] F. J. Higuera, S. Succi, and R. Benzi. Lattice gas dynamics with enhanced collisions.
Europhysics Letters, 9 (4):345349, 1989.
[7] P. L. Bhatnagar, E. P. Gross, and M Krook. A model for collision processes in
gases. I. Small amplitude processes in charged and neutral onecomponent systems.
Physical Review, 94 (3):511525, 1954.
[8]J.Buick,W.Easson,andC.Greated.Simulation of wave motion using a lattice gas
model. International Journal for Numberical Methods in Fluids, 22:313321, 1996.
[9] Sauro Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,
Clarendon Press, Oxford, 2001. ISBN 0198503989
[10] D. R. Noble, S. Chen, J. G. Georgiadis, and R. O. Buckius. A consistent hydro
dynamic boundary condition for the lattice Boltzmann method. Physics of Fluids, 7
(1):203209, 1995.
[11] S. Wolfram. Cellular automaton fluids 1: Basic theory. Journal of Statistical Physics,
45(3/4): p471529, 1986.
[12] Ghia, Ghia, and Shin, High Resolutions for incompressible flow using the navier
stokes equations and a multigrid method, Journal of Computational Physics, Vol.
48, p387411, 1982.
64
[13] Shuling Hou and Qisu Zou, Shiyi Chen, Gary Doolen, Allen C. Cogley. Simulation
of Cavity Flow by the Lattice Boltzmann Method. Journal of Computational Physics
118, p329347 1995.
[14] Pierre Lallemand, LiShi Luo. Lattice Boltzmann method for moving boundaries.
Journal of Computational Physics 184 (2003) 406421
[15] Paul J. Dellar, Incompressible limits of lattice Boltzmann equations using multiple
relaxation times, Journal of Computational Physics 190 (2003) 351370
[16] Xiao He and LiShi Luo, Lattice Boltzmann Model for the Incompressible Navier
Stokes Equation, Journal of Statistical Physics, Vol.88, Nos. 3/4, 1997
[17] Zhaoli Guo, Baochang Shi, and Nengchao Wang, Lattice BGK Model for Incom
pressible NavierStokes Equation, Journal of Computational Physics 165, 288306
(2000)
65
Appendices
66
Appendix A
FHP collision lookup table
This is the collision lookup table. It gives the two outstates corresponding to each
of the 128 instates (right in bracket).
0 [0 0]
1 [1 1]
2 [2 2]
3 [3 3]
4 [4 4]
5 [66 66]
6 [6 6]
7 [7 7]
8 [8 8]
9 [36 18]
10 [68 68]
11 [38 69]
12 [12 12]
13 [74 22]
14 [14 14]
15 [15 15]
16 [16 16]
17 [96 96]
18 [9 36]
19 [98 37]
20 [72 72]
21 [42 42]
22 [13 74]
23 [102 75]
24 [24 24]
25 [52 104]
26 [84 44]
27 [45 54]
28 [28 28]
29 [90 108]
30 [30 30]
31 [110 110]
32 [32 32]
33 [33 33]
34 [65 65]
35 [35 35]
36 [18 9]
37 [19 98]
38 [69 11]
39 [39 39]
40 [80 80]
41 [81 50]
42 [21 21]
43 [83 101]
44 [26 84]
45 [54 27]
46 [77 86]
47 [87 87]
48 [48 48]
49 [49 49]
50 [41 81]
51 [51 51]
52 [104 25]
53 [105 114]
54 [27 45]
55 [107 107]
56 [56 56]
57 [57 57]
58 [116 89]
59 [117 117]
60 [60 60]
61 [122 122]
62 [93 93]
63 [63 63]
67
64 [64 64]
65 [34 34]
66 [5 5]
67 [67 67]
68 [10 10]
69 [11 38]
70 [70 70]
71 [71 71]
72 [20 20]
73 [100 82]
74 [22 13]
75 [23 102]
76 [76 76]
77 [86 46]
78 [78 78]
79 [79 79]
80 [40 40]
81 [50 41]
82 [73 100]
83 [101 43]
84 [44 26]
85 [106 106]
86 [46 77]
87 [47 47]
88 [88 88]
89 [58 116]
90 [108 29]
91 [109 118]
92 [92 92]
93 [62 62]
94 [94 94]
95 [95 95]
96 [17 17]
97 [97 97]
98 [37 19]
99 [99 99]
100 [82 73]
101 [43 83]
102 [75 23]
103 [103 103]
104 [25 52]
105 [114 53]
106 [85 85]
107 [55 55]
108 [29 90]
109 [118 91]
110 [31 31]
111 [111 111]
112 [112 112]
113 [113 113]
114 [53 105]
115 [115 115]
116 [89 58]
117 [59 59]
118 [91 109]
119 [119 119]
120 [120 120]
121 [121 121]
122 [61 61]
123 [123 123]
124 [124 124]
125 [125 125]
126 [126 126]
127 [127 127]
68
Appendix B
Partial Matlab code I: d2q9_rotating_polygon.m
%
% Test of a D2Q9 Lattice BGK model on a square area w/ a rotating polygon.
% An implementation of the moving boundary with mass conservation.
% An implementation of the incompressible enhancement.
%
%% Define the general global variables
global cont_mode;
cont_mode = 0
global dimension; % The square domain
dimension = [1 1];
global n_poly; % Number of the sides of the polygon
n_poly = 3;
global center;
center = [.5 .5]; % Center of the rotating polygon
global N; % Number of fineness
N=100;
global diameter; % The diameter of the rotating object
diameter = 1/N + 0.9;
global D; % The hydraulic diameter
if n_poly == 3
H=N;
69
S = N*sqrt(3)/2*diameter;
D = 2*(H^2sqrt(3)/4*S^2)/(4*H+3*S)
elseif n_poly == 4
H=N;
L = H/sqrt(2);
D = 2*(H^2L^2)/(4*H+4*L)
end
global x y; % Define the domain length, width (including the ghost layer)
x = dimension(1)*N + 3;
y = dimension(2)*N + 3;
global speed; % Rotation angle speed
T=3000;
speed = pi/T;
global r; % Density
r=1;
global Re; % Reynolds number
Re=31;
global omega;
viscosity = r*D^2*speed/Re
omega = 1/(viscosity*3+1/2)
global F; % Distribution function matrix
F = zeros(x*y, 9);
global U; % Velocity field
U = zeros(x*y, 2);
% global SequenceU;
% SequenceU = zeros(x*y,2,500);
global Rho; % Density
70
Rho = zeros(x*y, 1);
global Fc; % Distribution function after the collision
Fc=F;
global c; % c is the unit speed
c=1;
global e;
e=[00]; %i=1
e = [e; 1 0; 0 1; 1 0; 0 1]; % i = 2,3,4,5
e = [e; 1 1; 1 1; 1 1; 1 1]; % i = 6,7,8,9
% 736
%\/
% \/
% 4  1 2
% /\
%/\
% 859
global w;
w = [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36];
global dx;
dx = 1; % the increment on positive x direction
global dy;
dy = x; % the increment on positive y direction
global opposite;
opposite = [1 4 5 2 3 8 9 6 7];
global q2 q3 q4 q5 q6 q7 q8 q9;
71
q2=[];
q3=[];
q4=[];
q5=[];
q6=[];
q7=[];
q8=[];
q9=[];
global i_polygon_update2 i_polygon_update3 i_polygon_update4 i_polygon_update5
i_polygon_update6 i_polygon_update7 i_polygon_update8 i_polygon_update9;
i_polygon_update2 = [];
i_polygon_update3 = [];
i_polygon_update4 = [];
i_polygon_update5 = [];
i_polygon_update6 = [];
i_polygon_update7 = [];
i_polygon_update8 = [];
i_polygon_update9 = [];
global A B C;
A=0;B=0;C=0;
%% Numbering order: x > y
% An example of this code with polygon = 3
% Starting from (1, 1) which is the lower left corner
%
%  y 
72
%  /\
%   / \  < rotating triangle
%  x  /____\  (counterclockwise rotation)
% (x+1) 
% (1)(x)
%
% _______________________________________________________________________
%  The index should start with i_ 
%  Those without the i_ are binary index that are for logical operation 
% ______________________________________________________________________
%
% actual boundary
%
% boundary ghost
% node node
%
%x x x x x

%  inside >
%
%  domain >
%
%
% Since we are using onthe node bounceback boundary, the inside nodes
% number will be (N  1)
% center of each space.
%
73
% ....
% ghost 0 1 ghost
%
% nodes showed 1234 N1NN+1
%
% < domain >
% < inside >
%
%
%% Define the domain
global domain;
domain = zeros(x*y,1);
% Each node will be assign a number:
% 0 for ghost_boundary
% 1 for boundary
% 2 for inside
%
% Domain = boundary + inside
%
for j = 2:y1
74
for i = 2:x1
domain((j1)*x+i) = 1;
end
end
for j = 3:y2
for i = 3:x2
domain((j1)*x+i) = 2;
end
end
%% Define the rotating polygon
global radius; % The radius of the inscribed circle
radius = diameter/2*sin( (pi  2*pi/n_poly)/2 );
global theta; % Initial polygon position
theta = pi*3/2;
% Define the sides of the polygon by the normal vector
global v_poly;
v_poly = zeros(n_poly, 2);
for i = 1:n_poly
v_poly(i,:) = [cos(theta+2*pi/n_poly*(i1)) sin(theta+2*pi/n_poly*(i1))];
end
global polygon;
polygon = zeros(x*y, 1);
for i = 2:x1
75
for j = 2:y1
v = [((i2)/N  center(1)) ((j2)/N  center(2))];
if v_poly*v? 0);
global i_boundary;
i_boundary = find( ( domainpolygon.*2 ) == 1 );
global v_poly_update;
v_poly_update = zeros(n_poly, 2);
for i = 1:n_poly
v_poly_update(i,:) = [cos(theta+2*pi/n_poly*(i1)) sin(theta+2*pi/n_poly
*(i1))];
end
polygon_update = zeros(x*y,1);
for i = 2:x1
for j = 2:y1
v = [((i2)/N  center(1)) ((j2)/N  center(2))];
76
if v_poly_update*v? 0);
i_non_domain = find( ( domainpolygon_update.*2 ) <= 0);
polygon_update2 = polygon_update;
polygon_update3 = polygon_update;
polygon_update4 = polygon_update;
polygon_update5 = polygon_update;
polygon_update6 = polygon_update;
polygon_update7 = polygon_update;
polygon_update8 = polygon_update;
polygon_update9 = polygon_update;
polygon_update2(1+dx:end) = polygon_update(1:enddx);
polygon_update3(1+dy:end) = polygon_update(1:enddy);
polygon_update4(1:enddx) = polygon_update(1+dx:end);
polygon_update5(1:enddy) = polygon_update(1+dy:end);
polygon_update6(1+dx+dy:end) = polygon_update(1:enddxdy);
polygon_update7(1+dy:enddx) = polygon_update(1+dx:enddy);
84
polygon_update8(1:enddxdy) = polygon_update(1+dx+dy:end);
polygon_update9(1+dx:enddy) = polygon_update(1+dy:enddx);
polygon_update2 = xor(polygon_update, polygon_update2)
& ( domainpolygon_update.*2 );
polygon_update3 = xor(polygon_update, polygon_update3)
& ( domainpolygon_update.*2 );
polygon_update4 = xor(polygon_update, polygon_update4)
& ( domainpolygon_update.*2 );
polygon_update5 = xor(polygon_update, polygon_update5)
& ( domainpolygon_update.*2 );
polygon_update6 = xor(polygon_update, polygon_update6)
& ( domainpolygon_update.*2 );
polygon_update7 = xor(polygon_update, polygon_update7)
& ( domainpolygon_update.*2 );
polygon_update8 = xor(polygon_update, polygon_update8)
& ( domainpolygon_update.*2 );
polygon_update9 = xor(polygon_update, polygon_update9)
& ( domainpolygon_update.*2 );
i_polygon_update2_backup = i_polygon_update2;
i_polygon_update3_backup = i_polygon_update3;
i_polygon_update4_backup = i_polygon_update4;
i_polygon_update5_backup = i_polygon_update5;
i_polygon_update6_backup = i_polygon_update6;
i_polygon_update7_backup = i_polygon_update7;
i_polygon_update8_backup = i_polygon_update8;
85
i_polygon_update9_backup = i_polygon_update9;
i_polygon_update2 = find(polygon_update2==1);
i_polygon_update3 = find(polygon_update3==1);
i_polygon_update4 = find(polygon_update4==1);
i_polygon_update5 = find(polygon_update5==1);
i_polygon_update6 = find(polygon_update6==1);
i_polygon_update7 = find(polygon_update7==1);
i_polygon_update8 = find(polygon_update8==1);
i_polygon_update9 = find(polygon_update9==1);
q2_backup = q2;
q3_backup = q3;
q4_backup = q4;
q5_backup = q5;
q6_backup = q6;
q7_backup = q7;
q8_backup = q8;
q9_backup = q9;
q2 = []; % q2 is the boundary length from the polygon to the nearest node
on direction 2
q3 = [];
q4 = [];
q5 = [];
q6 = [];
86
q7 = [];
q8 = [];
q9 = [];
q2 = q_rotating_polygon(i_polygon_update2, 2);
q3 = q_rotating_polygon(i_polygon_update3, 3);
q4 = q_rotating_polygon(i_polygon_update4, 4);
q5 = q_rotating_polygon(i_polygon_update5, 5);
q6 = q_rotating_polygon(i_polygon_update6, 6);
q7 = q_rotating_polygon(i_polygon_update7, 7);
q8 = q_rotating_polygon(i_polygon_update8, 8);
q9 = q_rotating_polygon(i_polygon_update9, 9);
%% Find and fill in newly released zone
Fc = F;
released = xor((polygon  polygon_update), polygon_update);
i_released = find(released==1);
count = numel(i_released);
for k = 1:count
coordi = [mod(i_released(k),x) floor(i_released(k)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
87
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
temp = [0 0 0 0 0 0 0 0 0];
for i = 2:9
temp1 = find(i_domain == i_released(k)+e(i,:)*[dx;dy]);
% Check if the previous node is in domain
temp2 = find(i_released == i_released(k)+e(i,:)*[dx;dy] );
% Check if the previous node is also a released node
temp3 = find(i_domain == i_released(k)e(i,:)*[dx;dy] );
% Check if the next node is in domain
temp4 = eval( [?find(i_polygon_update? num2str(i) ? == i_released(k))?] );
% Check if q value is available
if temp1 & isempty(temp2) & isempty(temp3) & temp4
temp(1) = temp(1) + 1;
temp(i) = 1;
% For F
F(i_released(k), i) = Fc(i_released(k)+e(i,:)*[dx;dy], i);
F(i_released(k), opposite(i)) = Fc(i_released(k)+e(i,:)*[dx;dy],
opposite(i));
end
88
if temp1 & temp3
% For F
F(i_released(k), i)=(Fc(i_released(k)+e(i,:)*[dx;dy], i)
+Fc(i_released(k)e(i,:)*[dx;dy], i) )/2;
F(i_released(k), opposite(i)) = ( Fc(i_released(k)+e(i,:)*[dx;dy],
opposite(i))
+Fc(i_released(k)e(i,:)*[dx;dy], opposite(i)) )/2;
end
end
% Fill in the rest particle
for i = 2:9
if temp(i) == 1
% For F
F(i_released(k), 1) = F(i_released(k), 1) + Fc(i_released(k)+e(i,:)
*[dx;dy], 1);
end
end
% For F
if temp(1)~=0
F(i_released(k), 1) = F(i_released(k), 1)/temp(1);
89
elseif temp(1)==0
fori=2:9
F(i_released(k), 1) = Fc(i_released(k), 1) + Fc(i_released(k)+e(i,:)
*[dx;dy], 1);
end
F(i_released(k), 1) = F(i_released(k), 1)/8;
end
end
%% Find and clear newly covered zone
covered = xor((polygon  polygon_update), polygon) & domain;
i_covered = find(covered==1);
count = numel(i_covered);
for k = 1:count
coordi = [mod(i_covered(k),x) floor(i_covered(k)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
temp = [0 0 0 0 0 0 0 0 0];
for i = 2:9
temp1 = find(i_domain == i_covered(k)+e(i,:)*[dx;dy] );
% Check if the next node is in domain
90
temp2 = find(i_covered == i_covered(k)+e(i,:)*[dx;dy] );
% Check if the next node is also a covered node
temp4 = eval( [?find(i_polygon_update? num2str(i) ?_backup
== i_covered(k))?] );
% Check if q value is available
if temp1 & isempty(temp2) & temp4
temp(i) = 1;
end
end
for i = 2:9
if temp(i) == 1
F(i_covered(k)+e(i,:)*[dx;dy], i) = Fc(i_covered(k)+e(i,:)*[dx;dy], i)
+ Fc(i_covered(k), i) * ( eval([?q? num2str(i) ?_backup( find(i_polygon_update?
num2str(i) ?_backup == i_covered(k) ) )? ]) + 1/2);
F(i_covered(k)+e(i,:)*[dx;dy], opposite(i)) = Fc(i_covered(k)
+e(i,:)*[dx;dy], opposite(i)) + Fc(i_covered(k), opposite(i))
* ( eval([?q? num2str(i) ?_backup( find(i_polygon_update? num2str(i) ?
_backup == i_covered(k) ) )? ]) + 1/2);
% boundary uncorrection
F(i_covered(k)+e(i,:)*[dx;dy], i) = F(i_covered(k)+e(i,:)*[dx;dy], i)
/ (3/2 + eval([?q? num2str(i) ?_backup( find(i_polygon_update? num2str(i) ?
_backup == i_covered(k) ) )? ]) );
F(i_covered(k)+e(i,:)*[dx;dy], opposite(i)) = F(i_covered(k)
91
+e(i,:)*[dx;dy], opposite(i)) / (3/2 + eval([?q? num2str(i) ?_backup( find
(i_polygon_update? num2str(i) ?_backup == i_covered(k) ) )? ]) );
end
end
end
B = B + sum(sum(F))/numel(i_domain);
%% Calculate Rho and U
Rho(:,:) = 0;
for i = 1:9
Rho(i_domain,:) = Rho(i_domain,:) + F(i_domain,i);
end
U(:,:) = 0;
for i = 1:9
U(i_domain,:) = U(i_domain,:) + F(i_domain,i)*e(i,:);
end
U(i_domain,1) = U(i_domain,1)./Rho(i_domain);
U(i_domain,2) = U(i_domain,2)./Rho(i_domain);
%% Calculate equilibrium distribution
92
for i = 1:9
Fc(i_domain,i) = Rho(i_domain,:)*w(i).*( 1 + 3/c^2*U(i_domain,:)*e(i,:)?
+ 9/(2*c^4)*(U(i_domain,:)*e(i,:)?).^2  3/(2*c^2)*(U(i_domain,1).^2
+ U(i_domain,2).^2) );
end
%% Collision
Fc(i_domain,:) = (1omega)*F(i_domain,:) + omega*Fc(i_domain,:);
%% Streaming (propagation)
%ForF
for i = 2:9
F(i_domain + e(i,:)*[dx; dy], i) = Fc(i_domain, i);
end
F(i_non_domain,:) = 0;
%% Boundary correction on F
93
%ForF
F(i_polygon_update2_backup, 2) = Fc(i_polygon_update2_backup, 2)
.*(1/2+q2_backup);
F(i_polygon_update2_backup, 4) = Fc(i_polygon_update2_backup, 4)
.*(1/2+q2_backup);
F(i_polygon_update4_backup, 4) = Fc(i_polygon_update4_backup, 4)
.*(1/2+q4_backup);
F(i_polygon_update4_backup, 2) = Fc(i_polygon_update4_backup, 2)
.*(1/2+q4_backup);
F(i_polygon_update3_backup, 3) = Fc(i_polygon_update3_backup, 3)
.*(1/2+q3_backup);
F(i_polygon_update3_backup, 5) = Fc(i_polygon_update3_backup, 5)
.*(1/2+q3_backup);
F(i_polygon_update5_backup, 5) = Fc(i_polygon_update5_backup, 5)
.*(1/2+q5_backup);
F(i_polygon_update5_backup, 3) = Fc(i_polygon_update5_backup, 3)
.*(1/2+q5_backup);
F(i_polygon_update6_backup, 6) = Fc(i_polygon_update6_backup, 6)
.*(1/2+q6_backup);
F(i_polygon_update6_backup, 8) = Fc(i_polygon_update6_backup, 8)
.*(1/2+q6_backup);
F(i_polygon_update8_backup, 8) = Fc(i_polygon_update8_backup, 8)
.*(1/2+q8_backup);
94
F(i_polygon_update8_backup, 6) = Fc(i_polygon_update8_backup, 6)
.*(1/2+q8_backup);
F(i_polygon_update7_backup, 7) = Fc(i_polygon_update7_backup, 7)
.*(1/2+q7_backup);
F(i_polygon_update7_backup, 9) = Fc(i_polygon_update7_backup, 9)
.*(1/2+q7_backup);
F(i_polygon_update9_backup, 9) = Fc(i_polygon_update9_backup, 9)
.*(1/2+q9_backup);
F(i_polygon_update9_backup, 7) = Fc(i_polygon_update9_backup, 7)
.*(1/2+q9_backup);
% Correction on newly released zone
count = numel(i_released);
for k = 1:count
coordi = [mod(i_released(k),x) floor(i_released(k)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
for i = 2:9
temp1 = find(i_domain == i_released(k)+e(i,:)*[dx;dy]);
% Check if the previous node is in domain
temp2 = find(i_released == i_released(k)+e(i,:)*[dx;dy] );
% Check if the previous node is also a released node
95
temp3 = find(i_released == i_released(k)e(i,:)*[dx;dy] );
% Check if the next node is a released node
temp4 = eval( [?find(i_polygon_update? num2str(i) ?
== i_released(k))?] );
% Check if q value is available
if temp1 & isempty(temp2) & isempty(temp3) & temp4
% For F
F(i_released(k), i) = F(i_released(k), i) * ( eval([?q? num2str(i)
?_backup( find(i_polygon_update? num2str(i) ?_backup == i_released(k)
+e(i,:)*[dx;dy]) )? ])  1/2);
F(i_released(k), opposite(i)) = F(i_released(k), opposite(i))
* ( eval([?q? num2str(i) ?_backup(find(i_polygon_update? num2str(i) ?
_backup == i_released(k)+e(i,:)*[dx;dy]) )? ])  1/2);
F(i_released(k)+e(i,:)*[dx;dy], i) = F(i_released(k)+e(i,:)
*[dx;dy], i) / (1/2 + eval([?q? num2str(i) ?_backup( find(i_polygon_update
? num2str(i) ?_backup == i_released(k)+e(i,:)*[dx;dy] ) )? ]) );
F(i_released(k)+e(i,:)*[dx;dy], opposite(i)) = F(i_released(k)
+e(i,:)*[dx;dy], opposite(i)) / (1/2 + eval([?q? num2str(i) ?_backup(
find(i_polygon_update? num2str(i) ?_backup == i_released(k)+e(i,:)
*[dx;dy] ) )? ]) );
end
end
96
end
% Correction on newly covered zone
count = numel(i_covered);
for k = 1:count
coordi = [mod(i_covered(k),x) floor(i_covered(k)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
temp = [0 0 0 0 0 0 0 0 0];
for i = 2:9
temp1 = find(i_domain == i_covered(k)+e(i,:)*[dx;dy] );
% Check if the next node is in domain
temp2 = find(i_covered == i_covered(k)+e(i,:)*[dx;dy] );
% Check if the next node is also a covered node
temp4 = eval( [?find(i_polygon_update? num2str(i) ?_backup
== i_covered(k))?] );
% Check if q value is available
if temp1 & isempty(temp2) & temp4
temp(i) = 1;
end
end
97
for i = 2:9
if temp(i) == 1
% For F
F(i_covered(k)+e(i,:)*[dx;dy], i) = F(i_covered(k)+e(i,:)*[dx;dy], i)
* (3/2 + eval([?q? num2str(i) ?_backup( find(i_polygon_update? num2str(i) ?
_backup == i_covered(k) ) )? ]) );
F(i_covered(k)+e(i,:)*[dx;dy], opposite(i)) = F(i_covered(k)+e(i,:)
*[dx;dy], opposite(i)) * (3/2 + eval([?q? num2str(i) ?_backup( find
(i_polygon_update? num2str(i) ?_backup == i_covered(k) ) )? ]) );
end
end
end
F_backup = F;
%% Boundary correction on the polygon for q<1/2
U_speed = U;
U_speed(:,:) = 0;
98
count = numel(q2);
for i = 1:count
if q2(i)<1/2
coordi = [mod(i_polygon_update2(i),x) floor(i_polygon_update2(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update2(i),:) = u?;
F(i_polygon_update2(i), opposite(2)) = (1/2+q2(i))*F_backup
(i_polygon_update2(i)+dx, opposite(2));
F(i_polygon_update2(i), 2) = (2*q2(i)/(1/2+q2(i)))*F_backup
(i_polygon_update2(i), opposite(2)) + (1/2q2(i))*F_backup
(i_polygon_update2(i)+dx, opposite(2));
F(i_polygon_update2(i)+dx, 2) = F_backup(i_polygon_update2(i), 2)
+ (1/2q2(i))/(1/2+q2(i))*F_backup(i_polygon_update2(i), opposite(2));
end
end
% Direction 3
count = numel(q3);
for i = 1:count
if q3(i)<1/2
coordi = [mod(i_polygon_update3(i),x) floor(i_polygon_update3(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
99
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update3(i),:) = u?;
F(i_polygon_update3(i), opposite(3)) = (1/2+q3(i))*F_backup
(i_polygon_update3(i)+dy, opposite(3));
F(i_polygon_update3(i), 3) = (2*q3(i)/(1/2+q3(i)))*F_backup
(i_polygon_update3(i), opposite(3)) + (1/2q3(i))*F_backup
(i_polygon_update3(i)+dy, opposite(3));
F(i_polygon_update3(i)+dy, 3) = F_backup(i_polygon_update3(i), 3)
+ (1/2q3(i))/(1/2+q3(i))*F_backup(i_polygon_update3(i), opposite(3));
end
end
% Direction 4
count = numel(q4);
for i = 1:count
if q4(i)<1/2
coordi = [mod(i_polygon_update4(i),x) floor(i_polygon_update4(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update4(i),:) = u?;
100
F(i_polygon_update4(i), opposite(4)) = (1/2+q4(i))*F_backup
(i_polygon_update4(i)dx, opposite(4));
F(i_polygon_update4(i), 4) = (2*q4(i)/(1/2+q4(i)))*F_backup
(i_polygon_update4(i), opposite(4)) + (1/2q4(i))*F_backup
(i_polygon_update4(i)dx, opposite(4));
F(i_polygon_update4(i)dx, 4) = F_backup(i_polygon_update4(i), 4)
+ (1/2q4(i))/(1/2+q4(i))*F_backup(i_polygon_update4(i), opposite(4));
end
end
% Direction 5
count = numel(q5);
for i = 1:count
if q5(i)<1/2
coordi = [mod(i_polygon_update5(i),x) floor(i_polygon_update5(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update5(i),:) = u?;
F(i_polygon_update5(i), opposite(5)) = (1/2+q5(i))*F_backup
(i_polygon_update5(i)dy, opposite(5));
F(i_polygon_update5(i), 5) = (2*q5(i)/(1/2+q5(i)))*F_backup
(i_polygon_update5(i), opposite(5)) + (1/2q5(i))*F_backup
(i_polygon_update5(i)dy, opposite(5));
101
F(i_polygon_update5(i)dy, 5) = F_backup(i_polygon_update5(i), 5)
+ (1/2q5(i))/(1/2+q5(i))*F_backup(i_polygon_update5(i), opposite(5));
end
end
% Direction 6
count = numel(q6);
for i = 1:count
if q6(i)<1/2
coordi = [mod(i_polygon_update6(i),x) floor(i_polygon_update6(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update6(i),:) = u?;
F(i_polygon_update6(i), opposite(6)) = (1/2+q6(i))*F_backup
(i_polygon_update6(i)+dx+dy, opposite(6));
F(i_polygon_update6(i), 6) = (2*q6(i)/(1/2+q6(i)))*F_backup
(i_polygon_update6(i), opposite(6)) + (1/2q6(i))*F_backup
(i_polygon_update6(i)+dx+dy, opposite(6));
F(i_polygon_update6(i)+dx+dy, 6) = F_backup(i_polygon_update6(i), 6)
+ (1/2q6(i))/(1/2+q6(i))*F_backup(i_polygon_update6(i), opposite(6));
end
end
102
% Direction 7
count = numel(q7);
for i = 1:count
if q7(i)<1/2
coordi = [mod(i_polygon_update7(i),x) floor(i_polygon_update7(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update7(i),:) = u?;
F(i_polygon_update7(i), opposite(7)) = (1/2+q7(i))*F_backup
(i_polygon_update7(i)dx+dy, opposite(7));
F(i_polygon_update7(i), 7) = (2*q7(i)/(1/2+q7(i)))*F_backup
(i_polygon_update7(i), opposite(7)) + (1/2q7(i))*F_backup
(i_polygon_update7(i)dx+dy, opposite(7));
F(i_polygon_update7(i)dx+dy, 7) = F_backup(i_polygon_update7(i), 7)
+ (1/2q7(i))/(1/2+q7(i))*F_backup(i_polygon_update7(i), opposite(7));
end
end
% Direction 8
count = numel(q8);
for i = 1:count
if q8(i)<1/2
103
coordi = [mod(i_polygon_update8(i),x) floor(i_polygon_update8(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update8(i),:) = u?;
F(i_polygon_update8(i), opposite(8)) = (1/2+q8(i))*F_backup
(i_polygon_update8(i)dxdy, opposite(8));
F(i_polygon_update8(i), 8) = (2*q8(i)/(1/2+q8(i)))*F_backup
(i_polygon_update8(i), opposite(8)) + (1/2q8(i))*F_backup
(i_polygon_update8(i)dxdy, opposite(8));
F(i_polygon_update8(i)dxdy, 8) = F_backup(i_polygon_update8(i), 8)
+ (1/2q8(i))/(1/2+q8(i))*F_backup(i_polygon_update8(i), opposite(8));
end
end
% Direction 9
count = numel(q9);
for i = 1:count
if q9(i)<1/2
coordi = [mod(i_polygon_update9(i),x) floor(i_polygon_update9(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
104
% For F
U_speed(i_polygon_update9(i),:) = u?;
F(i_polygon_update9(i), opposite(9)) = (1/2+q9(i))*F_backup
(i_polygon_update9(i)+dxdy, opposite(9));
F(i_polygon_update9(i), 9) = (2*q9(i)/(1/2+q9(i)))*F_backup
(i_polygon_update9(i), opposite(9)) + (1/2q9(i))*F_backup
(i_polygon_update9(i)+dxdy, opposite(9));
F(i_polygon_update9(i)+dxdy, 9) = F_backup(i_polygon_update9(i), 9)
+ (1/2q9(i))/(1/2+q9(i))*F_backup(i_polygon_update9(i), opposite(9));
end
end
%% Boundary correction on the polygon for q>=1/2
% F_backup2 = F;
% Direction 2
count = numel(q2);
for i = 1:count
if q2(i)>=1/2
coordi = [mod(i_polygon_update2(i),x) floor(i_polygon_update2(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
105
% For F
U_speed(i_polygon_update2(i),:) = u?;
F(i_polygon_update2(i), opposite(2)) = (q2(i)1/2)/(q2(i)+1/2)
*F_backup(i_polygon_update2(i), opposite(2)) + F_backup(i_polygon_update2(i)
+dx, opposite(2));
F(i_polygon_update2(i), 2) = (q2(i)1/2)/(q2(i)+1/2)*F_backup
(i_polygon_update2(i), 2) + 1/(q2(i)+1/2)*F_backup(i_polygon_update2(i),
opposite(2));
F(i_polygon_update2(i)+dx, 2) = 1/(q2(i)+1/2)*F_backup
(i_polygon_update2(i), 2);
end
end
% Direction 3
count = numel(q3);
for i = 1:count
if q3(i)>=1/2
coordi = [mod(i_polygon_update3(i),x) floor(i_polygon_update3(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update3(i),:) = u?;
F(i_polygon_update3(i), opposite(3)) = (q3(i)1/2)/(q3(i)+1/2)
*F_backup(i_polygon_update3(i), opposite(3)) + F_backup(i_polygon_update3
106
(i)+dy, opposite(3));
F(i_polygon_update3(i), 3) = (q3(i)1/2)/(q3(i)+1/2)*F_backup
(i_polygon_update3(i), 3) + 1/(q3(i)+1/2)*F_backup(i_polygon_update3(i),
opposite(3));
F(i_polygon_update3(i)+dy, 3) = 1/(q3(i)+1/2)*F_backup
(i_polygon_update3(i), 3);
end
end
% Direction 4
count = numel(q4);
for i = 1:count
if q4(i)>=1/2
coordi = [mod(i_polygon_update4(i),x) floor(i_polygon_update4(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update4(i),:) = u?;
F(i_polygon_update4(i), opposite(4)) = (q4(i)1/2)/(q4(i)+1/2)
*F_backup(i_polygon_update4(i), opposite(4)) + F_backup(i_polygon_update4(i)
dx, opposite(4));
F(i_polygon_update4(i), 4) = (q4(i)1/2)/(q4(i)+1/2)*F_backup
(i_polygon_update4(i), 4) + 1/(q4(i)+1/2)*F_backup(i_polygon_update4(i),
opposite(4));
107
F(i_polygon_update4(i)dx, 4) = 1/(q4(i)+1/2)*F_backup
(i_polygon_update4(i), 4);
end
end
% Direction 5
count = numel(q5);
for i = 1:count
if q5(i)>=1/2
coordi = [mod(i_polygon_update5(i),x) floor(i_polygon_update5(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update5(i),:) = u?;
F(i_polygon_update5(i), opposite(5)) = (q5(i)1/2)/(q5(i)+1/2)
*F_backup(i_polygon_update5(i), opposite(5)) + F_backup(i_polygon_update5(i)
dy, opposite(5));
F(i_polygon_update5(i), 5) = (q5(i)1/2)/(q5(i)+1/2)*F_backup
(i_polygon_update5(i), 5) + 1/(q5(i)+1/2)*F_backup(i_polygon_update5(i),
opposite(5));
F(i_polygon_update5(i)dy, 5) = 1/(q5(i)+1/2)*F_backup
(i_polygon_update5(i), 5);
end
108
end
% Direction 6
count = numel(q6);
for i = 1:count
if q6(i)>=1/2
coordi = [mod(i_polygon_update6(i),x) floor(i_polygon_update6(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update6(i),:) = u?;
F(i_polygon_update6(i), opposite(6)) = (q6(i)1/2)/(q6(i)+1/2)
*F_backup(i_polygon_update6(i), opposite(6)) + F_backup(i_polygon_update6(i)
+dx+dy, opposite(6));
F(i_polygon_update6(i), 6) = (q6(i)1/2)/(q6(i)+1/2)*F_backup
(i_polygon_update6(i), 6) + 1/(q6(i)+1/2)*F_backup(i_polygon_update6(i),
opposite(6));
F(i_polygon_update6(i)+dx+dy, 6) = 1/(q6(i)+1/2)*F_backup
(i_polygon_update6(i), 6);
end
end
% Direction 7
count = numel(q7);
109
for i = 1:count
if q7(i)>=1/2
coordi = [mod(i_polygon_update7(i),x) floor(i_polygon_update7(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update7(i),:) = u?;
F(i_polygon_update7(i), opposite(7)) = (q7(i)1/2)/(q7(i)+1/2)
*F_backup(i_polygon_update7(i), opposite(7)) + F_backup(i_polygon_update7(i)
dx+dy, opposite(7));
F(i_polygon_update7(i), 7) = (q7(i)1/2)/(q7(i)+1/2)*F_backup
(i_polygon_update7(i), 7) + 1/(q7(i)+1/2)*F_backup(i_polygon_update7(i),
opposite(7));
F(i_polygon_update7(i)dx+dy, 7) = 1/(q7(i)+1/2)*F_backup
(i_polygon_update7(i), 7);
end
end
% Direction 8
count = numel(q8);
for i = 1:count
if q8(i)>=1/2
coordi = [mod(i_polygon_update8(i),x) floor(i_polygon_update8(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
110
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
U_speed(i_polygon_update8(i),:) = u?;
F(i_polygon_update8(i), opposite(8)) = (q8(i)1/2)/(q8(i)+1/2)
*F_backup(i_polygon_update8(i), opposite(8)) + F_backup(i_polygon_update8(i)
dxdy, opposite(8));
F(i_polygon_update8(i), 8) = (q8(i)1/2)/(q8(i)+1/2)*F_backup
(i_polygon_update8(i), 8) + 1/(q8(i)+1/2)*F_backup(i_polygon_update8(i),
opposite(8));
F(i_polygon_update8(i)dxdy, 8) = 1/(q8(i)+1/2)*F_backup
(i_polygon_update8(i), 8);
end
end
% Direction 9
count = numel(q9);
for i = 1:count
if q9(i)>=1/2
coordi = [mod(i_polygon_update9(i),x) floor(i_polygon_update9(i)/x)+1];
v = [((coordi(1)2)/N  center(1)) ((coordi(2)2)/N  center(2))];
u = [v(2) v(1)];
u = u*speed*norm(v)/norm(u);
% For F
111
U_speed(i_polygon_update9(i),:) = u?;
F(i_polygon_update9(i), opposite(9)) = (q9(i)1/2)/(q9(i)+1/2)
*F_backup(i_polygon_update9(i), opposite(9)) + F_backup(i_polygon_update9(i)
+dxdy, opposite(9));
F(i_polygon_update9(i), 9) = (q9(i)1/2)/(q9(i)+1/2)*F_backup
(i_polygon_update9(i), 9) + 1/(q9(i)+1/2)*F_backup(i_polygon_update9(i),
opposite(9));
F(i_polygon_update9(i)+dxdy, 9) = 1/(q9(i)+1/2)*F_backup
(i_polygon_update9(i), 9);
end
end
polygon = polygon_update;
%% Bounce back boundary condition
% This is mainly for the outer circular cylinder wall
%ForF
Fc(i_boundary,:) = F(i_boundary,:);
F(i_boundary,2) = Fc(i_boundary,4);
F(i_boundary,4) = Fc(i_boundary,2);
F(i_boundary,3) = Fc(i_boundary,5);
112
F(i_boundary,5) = Fc(i_boundary,3);
F(i_boundary,6) = Fc(i_boundary,8);
F(i_boundary,8) = Fc(i_boundary,6);
F(i_boundary,7) = Fc(i_boundary,9);
F(i_boundary,9) = Fc(i_boundary,7);
%% Boundary uncorrection on F
F(i_polygon_update2, 2) = F(i_polygon_update2, 2)./(1/2+q2);
F(i_polygon_update2, 4) = F(i_polygon_update2, 4)./(1/2+q2);
F(i_polygon_update4, 4) = F(i_polygon_update4, 4)./(1/2+q4);
F(i_polygon_update4, 2) = F(i_polygon_update4, 2)./(1/2+q4);
F(i_polygon_update3, 3) = F(i_polygon_update3, 3)./(1/2+q3);
F(i_polygon_update3, 5) = F(i_polygon_update3, 5)./(1/2+q3);
F(i_polygon_update5, 5) = F(i_polygon_update5, 5)./(1/2+q5);
F(i_polygon_update5, 3) = F(i_polygon_update5, 3)./(1/2+q5);
F(i_polygon_update6, 6) = F(i_polygon_update6, 6)./(1/2+q6);
F(i_polygon_update6, 8) = F(i_polygon_update6, 8)./(1/2+q6);
F(i_polygon_update8, 8) = F(i_polygon_update8, 8)./(1/2+q8);
F(i_polygon_update8, 6) = F(i_polygon_update8, 6)./(1/2+q8);
F(i_polygon_update7, 7) = F(i_polygon_update7, 7)./(1/2+q7);
F(i_polygon_update7, 9) = F(i_polygon_update7, 9)./(1/2+q7);
F(i_polygon_update9, 9) = F(i_polygon_update9, 9)./(1/2+q9);
113
F(i_polygon_update9, 7) = F(i_polygon_update9, 7)./(1/2+q9);
114