Transcript Artificial Compressibility Method - Beijing Computational Science

Artificial
Compressibility
Method and
知彼知己者
百戰不殆
Lattice
Boltzmann
Method
If you know your
enemies and know
yourself,
you can win a hundred battles without a single loss
Similarities and Differences
Taku Ohwada
(大和田 拓）
Department of Aeronautics & Astronautics, Kyoto University
（京都大学大学院工学研究科航空宇宙工学専攻）
Collaborators : Prof. Pietro Asinari, Mr. Daisuke Yabusaki
May 4, 2011,
Spring School on the lattice Boltzmann Method
Beijing Computational Science Research Center
May 2-6
1
0. What is a good numerical method ?
• Performance
• Cost (CPU)
• Education (Human CPU)
2
Many paths to the summit
1. INTRODUCTION
4
Kinetic methods for fluid-dynamic equations
Gas Kinetic Scheme
Lattice Boltzmann Method
Why Kinetic ?
Boltzmann Eq.  Euler, Navier-Stokes
The path is INDIRECT !
5
An indirect method is not
always your best choice.
6
Extraction of essence
Subtraction rather than addition
7
Why Kinetic ?
Case of compressible flows
Kinetic gadget
fM
(i  ui ) 2
fM 
exp(
)
3/ 2
(2RT )
2RT

f M yields the flux for the Euler equations.
u


 1 




2
u  p

       f M d
 [  (e  u 2 / 2)  p ]u 
 2 / 2 




The discontinuities at cell-interfaces
produce numerical dissipation, which
suppresses spurious oscillations around
shock waves……….. Shock Capturing !
Riemann Problem
Riemann Problem
Riemann Problem
Kinetic Flux Splitting
f
f

t
x
Characteristics :
x   t  Const.
f
0
 1 


   2  f d
 / 2 



Prerequisite of Gas Kinetic Scheme is
Taylor expansion !!!!
Experiment using undergraduate students.
11
Gallery of Undergraduate Students ‘ Works
Euler
Blasius flow
Navier-Stokes:
Any asymptotic
method is NOT
employed.
Pressure distribution along the wall
Present
13
Case of incompressible flows
Gas kinetic Scheme
Lax-Wendroff
M 0
No kinetic ingredient !!!!
14
Case of incompressible flows
LBM
f (t , xi ) (  1,2,...,Q)
f (t  1, xi）＝ f (t , xi  V )

  g (h(t , xi  V ))  f (t , xi  V )

h  u, v, P
f (t , xi )  u(t , xi ), v(t , xi ), P(t , xi )
15
f (t  1, xi）＝ f (t , xi  V )

  g (h(t , xi  V ))  f (t , xi  V )

 1
f (t  1, xi）＝ g (h(t , xi  V ))
LBM -> Lattice Kinetic Scheme (LKS)
No kinetic ingredient !!
LKS ->
Lattice Scheme (LK)
-> ACM
16
Incompressible case:
More difficult than compressible case !!!
ui
 0,
xi
ui
ui P
 2u i
uj


t
x j xi
x 2j
u j ui
2P

2
xi
xi x j
17
1
A
LBM !
Poisson Free…
• Poisson Free !!!
• 2nd order accurate BB
• Small time step t ~ (x) 2
• Parallel computation !!!
Bouzidi, Firdaouss, Lallemand (2001)
Ginzburg, d’Humières (2003)
18
LBM solves INSE via
Artificial Compressibility Equations
Prof. Asinari ’s morning lecture
Chapman-Enskog expansion
Hilbert expansion (diffusive scale)
P ui
3(x)

0
t
xi
ui
ui P
 2ui
uj

 2
t
x j xi
x j
2
LBM
19
Artificial Compressibility Method (ACM)
P ui
k

0
t
xi
(Chorin,1967)
(Témam, 1969)
ui
ui
 2 ui
P
uj


t
x j
xi
x 2j
usually
k  0.1~10
LBM
k  3(x)
2
20
Considering the fact that the lattice Boltzmann method
starts with the kinetic theory and has been derived to
conserve high-order isotropy, the lattice Boltzmann
method should be more accurate than the artificial
compressibility method in capturing pressure waves.
He, Doolen, Clark (JCP2002)
ACM:
LBM:
Macroscopic
(356 papers)
Kinetic
(4053 papers)
21
Devil’s Project LBM-ACM
LBM
Chapman-Enskog
Expansion
ACM
• Lattice Structure
• Collocated Grid
LBM
ACM
Finite Difference
(Finite Volume)
Kinetic
22
2. Numerical Computation of ACM
Cartesian Grid
D2Q9
Time Step
Finite Difference
Space : Central Difference
Time : Semi-Implicit
23
ui
t
P
t
0
x
24
P
1  u v 


(  x)
0
2 
t
   x x 
u
u
u P
 2u  2u
u
v

 ( 2 
)
2
t
x
y x
x
y
v
v
v P
 2v  2v
u
v

 ( 2 
)
2
t
x
y y
x
y
h (i ) ( j )
h(( i  1) , j )  h(( i  1) , j )
( x , y ) ~  x h(i , j ) 
x
2
h (i ) ( j )
h(i , ( j  1) )  2h(i , j )  h(i , ( j  1) )
( x , y ) ~  yy h(i , j ) 
2
y
2
25
Basic form of ACM
P
1  u v 


(  x)

0
t
  2  x x 
u
u
u P
1  2u  2u
u
v


(

)
t
x
y x
Re x 2 y 2
v
v
v P
1  2v  2v
u
v


(

)
t
x
y y
Re x 2 y 2
uijn  u(nt, i , j )
uijn 1  uijn
t
n 1
n
vij  vij
t
1
n
u  u  v  u  P 
( xx   yy )uij
Re
n
ij
n
x ij
n
ij
n
y ij
n
x ij
1
u  v  v  v  P 
( xx   yy )vijn
Re
n
ij
n
x ij
Pijn1  Pijn
t

n
ij
n
y ij
n
y ij
2
 u
 v
1

n 1
x ij
n 1
y ij
 0
26
4th order accurate div u for compact stencil
 xu ~  xu 
 yv ~  yv 
2
6
2
6
 xu   y v  O( )
 yv ~  yv 
D(u, v)   xu 
2
6
 yyyv  O( 4 )
 xxxu   xxyv  O( 2 )
2
 xu ~  xu 
 xxxu  O( 4 )
2
6
2
6
 xx y v  O( 4 )
 yy xu  O( 4 )
 xx y v   y v 
2
6
 yy xu   xu   y v  O( 4 )
Test Problems
• Generalized Taylor-Green Vortex (2D, 3D)
• Circular Couette Flow (2D)
• Flow Past a Cylinder in a Channel (2D)
• Lid-driven Cavity Flow (2D, 3D)
• Flow Past a Sphere in Uniform Flow (3D)
28
2D Generalized Taylor-Green
Periodic
Boundary
29
Time history of L1 error
  1 Re  0.001
30
Convergence Rate (t=100)
31
3D Generalized Taylor-Green (3D-GTG)
32
Circular Couette Flow
R
5R
33
Error of Velocity Uθ (ACM)
34
Error of Velocity Uθ (MRT)
35
Error of Pressure P (ACM)
36
Error of Pressure P (MRT)
37
Comparison
Uθ
P
ACM
MRT
ACM
MRT
38
Convergence rate
39
The Flow Past a Cylinder
M. Schäfer, S. Turek, (1996)
Non-Slip Boundary
Poiseuille Flow
2.1D
D
2D
0th-order
extrapolation
Non-Slip Boundary
40
t=100
|div u| (Re=100)
41
|div u| (Re=100)
t=100
42
Re=100 (unsteady)
43
* M. Schäfer, S. Turek, (1996)
LBM: Mussa, Asinari, Luo,
JCP 228 (2009)
44
Adaptive Mesh Refinement
(Re=100)
D: the diameter of the cylinder
Poiseuille Flow
D
0th-order
extrapolation
Simple Interpolation
45
t=100
Velocity u
46
t=100
Velocity v
47
t=100
Pressure
48
2D Lid-driven Cavity Flow
Top boundary
49
2D Lid-driven Cavity Flow
CPU Time
(129×129, 100000step)
(intel Corei7, openMP)
ACM
MRT
55.01 sec 405.51 sec
×7.37
356 papers
4053 papers
50
3D Lid-driven Cavity Flow
Grid : 100³
Lo D. C. , Murugesan K. , Young D. L. , Int. J. Numer. Meth. Fluids (2005)
51
Flow past a sphere
Re=300
52
Conclusion
Comparisons of LBM and ACM
• Taylor-Green Vortex Flow
• Circular Couette Flow
• Flow past a Cylinder
• Lid-driven Cavity Flow
LBM – ACM is NOT decisively positive !!!
Performance of ACM in 3D problems
• Lid-driven Cavity Flow
• Flow past a Cylinder
ACM is capable of practical 3D simulation
53
54
Intermission
Re=2000
55
Re=500
56
3. Theory of ACM
57
Convergence of the solution of ACE
P ui
k

0
t
xi
ui
ui 1 uk
 2ui
P
uj

ui 
 2
t
x j 2 xk
xi
x j
to the solution of INSE in the limit of k  0 .
Témam (1969)
How fast ?
A mathematical analysis shows
(Moise and Ziane, Journal of Dynamics and Differential Equations, Vol. 13, No. 2, 2001)
| uik  uiNS |H1  Ck 1/ 2
Numerical observations show that the error is O(k) !!
59
Diffusive mode of ACE solution
P ui
k

0
t
xi
0  k  1
ui
ui P
 2ui
uj

 2
t
x j xi
x j
Assumption
ui ui
ui ~
~
~ O(1)
x j
t
P P
~
~ O(1)
x j t
60
ui  u  ku  k u  
0
i
1
i
2
2
i
P  P  kP  k P  
0
1
2
ui0
 0,
xi
0
i
2
0
0
i
2
j
ui1
Oseen type
P 0

,
xi
t
u
u
P
1 u
0 u
uj
uj


0
t
x j
x j xi
x
1
i
0
i
INSE
O(1)
u
u
P
0 u
uj


0
t
x j xi
x
0
i
2
1
i
1
2 1
i
2
j
O(k )
The intrinsic error of ACM O(k )
k ~ (x) ~ t
2
Acoustic mode of ACE solution
P ui
k

0
t
xi
k
ui
ui P
 2ui
uj

 2
t
x j xi
x j
P ui

0
t
xi
ui P

0
t xi
2
2P
2  P
c
t 2
x 2j
c
1
k
Characteristic time ~ Characteristic length / Speed of wave
 1/ c

k
62
Initial data for ACE = Initial data for INSE
(ui (t  0), P(t  0))  (u (t  0), P (t  0))
0
i
wi  ui  ui0 , q  P  P0
0
wi (t  0)  0, q(t  0)  0
wi
P 0 q
k(

)
0
t
t
xi
0
2
wi

w

u

w

wi

q
0
i
i
i
uj
 wj
 wj


t
x j
x j
x j xi
x 2j
P 0
t
acts as an initial impact !!!
q
P 0

t
t
wi
0
t
Magnitude of acoustic mode
q
P 0

~ O(1)
t
t
t  k
q
~q


wi
P 0 q
k(

)
0
t
t
xi
q ~ O( k )
wi ~ O(k )
64
q
ui (0,.)  ui0 (0,.)
ui (0,.)  ui0 (0,.)  ku i1 (0,.)
P (0,.)  P 0 (0,.)
P (0,.)  P 0 (0,.)  kP 1 (0,.)
q
4. Improvement of ACM
66
4.1 Suppression of acoustic mode
q wi
k

0
t
xi
wi q
 2wi

 2
t xi
x j
  t / k Q  q / k Wi  wi / k
 2Q  2Q
 3Q
 2  k 2
2

x j
x j 
1. Bulk viscosity
q wi
k

0
t
xi
wi q
 2wi
 w j

 2  
t xi
x j
xi x j
 2 Q  2Q
 3Q
 2  k (   ) 2
2

x j
x j 
67
2. Dashpot
wi
q
k(
  q) 
0
t
xi
wi q
 2wi

 2
t xi
x j
 2Q
 2Q
Q
 3Q
 (1    k ) 2   k
 k 2
2

x j

x j 
Dashpot
68
1. Bulk viscosity
k
P ui

0
t
xi
ui
ui
 2 ui
P
uj


t
x j xi
x 2j
k
P ui

0
t
xi
ui
ui
 2 ui
P
 ui
uj




t
x j xi
x 2j
xi xi
2. Dashpot
k
P ui

0
t
xi
ui
ui
 2 ui
P
uj


t
x j xi
x 2j
k(
u
P
  P)  i  0
t
xi
ui
ui
 2 ui
P
uj


t
x j xi
x 2j
69
Time history of L1 error
  1 Re  0.001
 1
70
4.2 Suppression of Checkerboard Instability
1
Acoustic mode killer
2
Checkerboard killer
3
stability
71
v
Staggered grid
P
u
u
v
u, v, P
u, v, P
u, v, P
u, v, P
Collocated grid
72
2D-LID Re=5000
P-contour
73
P-contour
ACM (s=0,2)
(256×256)
C.-H. Bruneau*
(2048×2048)
*C.-H. Bruneau & M. Saad, Comput. & Fluids 35, 326-348 (2006)
74
Stream line
Pressure
75
Stream line
Pressure
76
4.3 Enhancer of stability
1
Acoustic mode killer
2
Checkerboard killer
3
stability
77
Linear Stability Analysis
Pressure mode
Velocity modes
S=0
S=2
Checkerboard
killer
~
k
The normalized wave number
78
S=0
S=2
79
Stability in Lid-Driven
ACM
ACM
MRT
LBGK
Grid
(s=0)
(s=2)
32×32
×
○
×
×
64×64
×
○
×
×
96×96
○
○
○
×
128×128
○
○
○
×
256×256
○
○
○
○
(χ=0.0378) (χ=0.0002)
80
Richardson Extrapolation in the Mach Numbers
Diffusive mode
Linear Combination
81
High order accurate ACM
Spatial discretization
4th order accuracy
(5point central finite difference)
Time-marching
2nd order accuracy
(Semi-implicit Runge-Kutta)
82
4th order accurate div u for compact stencil
 xu ~  xu 
 yv ~  yv 
2
6
2
6
 xu   y v  O( )
 xu ~  xu 
D(u, v)   xu 
2
6
 yyyv  O( 4 )
 xxxu   xxyv  O( 2 )
2
 yv ~  yv 
 xxxu  O( 4 )

2
2
6
6
 xx y v  O( 4 )
 yy xu  O( 4 )
 xx y v   y v 
2
6
 yy xu   xu   y v  O( 4 )
ui0
 0,
xi
nd order
4th
2
order accuracy
accuracy
u
u
P
0 u
uj


0
t
x j xi
x
0
i
0
i
2
0
0
i
2
j
O(1)
O( )
ui1
P 0
 
xi
t
2
0
1
2 1
1
ui1

u

u

ui

P
1
0
i
i
uj
uj


0
2
t
x j
x j xi
x j
ui2
P1
 
xi
t
u
u
P
2 u
0 u
1 u
uj
uj


 u j
t
x j
x j xi
x
x j
2
i
0
i
O( )
4
2
i
2
2 1
i
2
j
1
i
3D Generalized Taylor-Green (3D-GTG)
85
Curved Solid Boundary
LBM
e.g. Interpolation Bounce-Back
Bouzidi, Firdaouss, Lallemand (2001)
Ginzburg, d’Humières (2003)
ACM:
Macroscopic data
Interpolation or extrapolation
86
y
N
Fluid
Quadratic
extrapolation
W
B
P
W
E
B
Solid
Body
S
x
P
E
x
at B
87
y
N
W
B
P
Fluid
Quadratic
extrapolation
E
W B
Solid
Body
P
E
x
at B
S
x
at P
88
Why not ACM?
Since the ACM does not employ any kinetic theory gadget,
it is much easier than the LBM. Up to now, any decisive
inferiority of ACM to LBM has not been found. Conversely,
superiority of ACM over LBM has been found in some
fundamental test problems.
Therefore, it is highly recommended to master ACM
before learning LBM.
89
Richardson extrapolation = 累遍増約術
Its usefulness for practical computations can hardly be
overestimated.
Birkoff & Rota, Ordinary differential equations (1978).
Lewis Fry Richardson, ``Approximate arithmetical solution by finite
differences of physical problems including differential equations, with
an application to the stresses in a masonry dam ‘’, Phil. Trans. Royal
Soc. London, Series A 210: 307–357 (1910).
谢谢您为您的关注
91
Finite Volume ACM
93
94
95
96
Excitation of acoustic mode
P ui

0
t
xi
ui  u  ku  k u  
k
P  P  kP  k P  
ui
ui P
 2ui
uj

 2
t
x j xi
x j
0
i
0
1
i
2
1
2
0
The initial data for (ui
The initial data for
2
i
2
, P0 ) must be compatible with INSE.
(uin , Pn ) (n  1,2,3,...)
ui1
P 0

,
xi
t
0
1
2 1
1
ui1

u

u

ui

P
1
0
i
i
uj
uj


0
2
t
x j
x j xi
x j
troublesome