Transcript Document

Sino-German Symposium on Advanced Numerical Methods for Compressible Fluid Mechanics
and Related Problems, May 21-27, 2014, Beijing, China
Discrete unified gas-kinetic
scheme for compressible flows
Zhaoli Guo
(Huazhong University of Science and Technology, Wuhan, China)
Joint work with Kun Xu and Ruijie Wang (Hong Kong University of Science and Technology)
Outline
• Motivation
• Formulation and properties
• Numerical results
• Summary
Motivation
Non-equilibrium flows covering different flow regimes
Re-Entry Vehicle
10-3
Inhalable particles
Chips
10-2
10-1
100
10
Challenges in numerical simulations
Modern CFD:
• Based on Navier-Stokes equations
• Efficient for continuum flows
• does not work for other regimes
Particle Methods: (MD, DSMC… )
• Noise
• Small time and cell size
• Difficult for continuum flows / low-speed non-equilibrium flows
Method based on extended hydrodynamic models :
• Theoretical foundations
• Numerical difficulties (Stability, boundary conditions, ……)
• Limited to weak-nonequilibrium flows
Lockerby’s test (2005, Phys. Fluid)
 = const
the most common high-order continuum equation sets (Grad’s 13
moment, Burnett, and super-Burnett equations ) cannot capture
the Knudsen Layer, Variants of these equation families have,
however, been proposed and some of them can qualitatively
describe the Knudsen layer structure … the
agreement
slight
quantitative
with kinetic theory and DSMC data is only
A popular technique: hybrid method
Limitations
Numerical rather than physical
Artifacts
MD
NS
Time coupling
Dynamic scale changes
Hadjiconstantinou
Int J Multiscale Comput Eng 3 189-202, 2004
Hybrid method is inappropriate for problems
with dynamic scale changes
Efforts based on kinetic description of flows
# Discrete Ordinate Method (DOM) [1,2]:
•
•
•
•
Time-splitting scheme for kinetic equations (similar with DSMC)
dt (time step) <  (collision time)
dx (cell size) <  (mean-free-path)
numerical dissipation  dt
Works well for highly non-equilibrium flows, but encounters difficult for
continuum flows
# Asymptotic preserving (AP) scheme [3,4]:
• Consistent with the Chapman-Enskog representation in the continuum limit
(Kn  0)
• dt (time step) is not restricted by  (collision time)
• at least 2nd-order accuracy to reduce numerical dissipation [5]
Aims to solve continuum flows, but may encounter difficulties for free molecular flows
[1] J. Y. Yang and J. C. Huang, J. Comput. Phys. 120, 323 (1995)
[2] A. N. Kudryavtsev and A. A. Shershnev, J. Sci. Comput. 57, 42 (2013).
[3] S. Pieraccini and G. Puppo, J. Sci. Comput. 32, 1 (2007).
[4] M. Bennoune, M. Lemo, and L. Mieussens, J. Comput. Phys. 227, 3781 (2008).
[5] K. Xu and J.-C. Huang, J. Comput. Phys. 229, 7747 (2010)
Efforts based on kinetic description of flows
# Unified Gas-Kinetic Scheme (UGKS) [1]:
• Coupling of collision and transport in the evolution
• Dynamicly changes from collision-less to continuum according to the local
flow
• The nice AP property
A dynamic multi-scale scheme, efficient for multi-regime flows
In this report, we will present an alternative
kinetic scheme (Discrete Unified Gas-Kinetic
Scheme), sharing many advantages of the UGKS
method, but having some special features .
[1] K. Xu and J.-C. Huang, J. Comput. Phys. 229, 7747 (2010)
Outline
• Motivation
• Formulation and properties
• Numerical results
• Summary
# Kinetic model (BGK-type)
¶ t f + x ×Ñ f = Wº -
Distribution function
1é
f - f eq ù
ë
û
t
f = f (x , x, t )
Particel velocity
Equilibrium:
f eq = f eq [x,W (x , t ), J (x , t ),...]
Flux
Conserved variables
Example:
f
f
eq
eq
= f
= f
M
S
é c2
r
=
exp ê(3 + K )/ 2
êë 2R T
(2p R T )
= f
f eq = f ES =
M
é
c ×q
ê1 + (1 - P r)
êë
5pR T
ù
ú
ú
û
æ c2
öù
÷ú
çç
÷ú
çè R T - 5 ÷
ø
û
é 1
ù
exp ê- c ×L ×c ú
êë 2
ú
û
det (2p L )
r
æ
1 ö
÷s
L ij = RT dij + ççç1 ÷
ij
è
ø
Pr ÷
Maxwell (standard BGK)
Shakhov model
ES model
ér ù
ê ú
W = êêr u ú
ú=
êr E ú
êë ú
û
Conserved variables
òj
(x)fd X
é1 ù
ê
ú
ú
j ( x) = ê
x
ê
ú
ê1 x 2 ú
ê2
ú
ë
û
Conservation of the collision operator
ò j Wfd X =
1
j [f - f eq ]fd X = 0
ò
t
A property: for any linear combination of f and f eq , i.e.,
f ¢ = (1 - b ) f + b f eq
The conservation variables can be calculated by
W =
ò j (x)f ¢d X
# Formulation: A finite-volume scheme
¶ t f + x ×Ñ f = Wº -
1é
f - f eq ù
ë
û
t
j+1/2
Point 1: Updating rule for cell-center distribution function
1. integrating in cell j:
j
f jn + 1 - f jn +
j+1
D t é n + 1/ 2
Dt é n + 1
x ëêf j + 1/ 2 - f jn- +1/1/22 ù
=
Wj + Wnj ù
ú
ë
û
û
Dx
2
Mid-point
Trapezoidal
2. transformation:
Dt
f%= f W
2
Dt
2t - D t %
2D t
f%+ = f +
W=
f +
f eq
2
2t + D t
2t + D t
3. update rule:
D t é n + 1/ 2
+ ,n
n+1
f%
= f%
+
x êëf j - 1/ 2 - f jn++1/1/22 ù
j
j
ú
û
Dx
Key: distribution function at cell interface
12
Point 2: Evolution of the cell-interface distribution function
f jn++1/1/22
How to determine
Again
j+1/2
1
¶ t f + x ×Ñ f = Wº - éëf - f eq ù
û
t
s j+
1. integrating along the characteristic line
f jn++1/1/22 - f jn =
h é n + 1/ 2
n
ù
êëWj + 1/ 2 + W (x j + 1/ 2 - xh) ú
û
2
h =
Dt
2
j
j+1
explicit
Implicit
2. transformation:
h
W
2
f = f -
So
f j n+ +1/1/2 2 = f
+ ,n
(x j + 1/ 2
f
+
= f +
ìï f
ïï
- xh ) = ïí
ïï
ïïî f
+ ,n
+ ,n
h
2t - h
2h
W=
f +
f eq
2
2t + h
2t + h
(x j ) + 21 (x j + 1/ 2 - x j - xD t )s j+ ,
x³ 0
(x j + 1 ) + 21 (x j + 1/ 2 - x j + 1 - xD t )s j++ 1,
x< 0
3. original:
f jn++1/1/22 =
2t
h
f j n+ +1/1/2 2 +
f eq (W jn++1/1/22 )
2t + h
2t + h
W jn++1/1/22 =
ò y (x)f
n + 1/ 2
j + 1/ 2
(x)d x
s j+ = Ñ f j+
Slope
# Boundary condition
Bounce-back
W x ·u
f (x w , xi , t + h ) = f (x w , - xi , t + h ) + 2r w i i w ,
wi R T
xi ×n > 0
é
2
r w = êê1 RT
êë
- 1
ù
å W i xi ×u w úú
xi ×n > 0
ú
û
n
é
ù
ê
ú
ê å wi f ( xi ) + 2 å wi f ( xi ) ú.
xi ×n < 0
êë xi ×n = 0
ú
û
Diffuse Scatting
n
f (x w , xi , t + h) = f eq (xi ; r w , u w ),
å
(xi ×n ) f eq ( xi , r w , u w ) =
xi ×n > 0
å
xi ×n < 0
rw
xi ×n > 0
(xi ×n ) f ( x w , xi , t + h)
# Properties of DUGKS
1. Multi-dimensional
• It is not easy to device a wave-based multi-dimensional scheme based on
hydrodynamic equations
• In the DUGKS, the particle is tracked instead of wave in a natural way (followed by its
trajectory)
2. Asymptotic Preserving (AP)
(a) time step (t) is not limited by the particle collision time ():
(b) in the continuum limit (t >> ):
f jn++1/1/22 =
2t
h
f j n+ +1/1/2 2 +
f eq (W jn++1/1/22 )
2t + h
2t + h
in the free-molecule limit: (t << ):
f jeq+ 1/ 2 - t Df jeq+ 1/ 2 + h¶ t f jeq+ 1/ 2
Chapman-Ensokg expansion
f jn++1/1/22 = f n (x j + 1/ 2 - xh )
(c) second-order in time; space accuracy can be ensured by choosing linear or
high-order reconstruction methods
# Comparison with UGKS
Unified GKS (Xu & Huang, JCP 2010)
Starting Point:
¶ t f + x ×Ñ f = Wº ¶ tW + Ñ ×F = 0
1é
f - f eq ù
ë
û
t
F =
ò xy f d x
Macroscopic flux
Updating rule:
f jn + 1 - f jn +
1
Dx
W jn + 1 - W jn +
ò
tn + 1
tn
Dt é
x éëf j + 1/ 2 (t ) - f j - 1/ 2 (t ) ù
dt
=
Wj ( f jn + 1,W
û
ë
2
n+1
) + Wj ( f jn ,W n ) ù
û
j+1/2
tn + 1
1
F j + 1/ 2 (t )dt = 0
| V j | òtn
j
j+1
If the cell-interface distribution f(t) is known, the update both f and W can
be accomplished
Unified GKS (cont’d)
j+1/2
1
¶ t f + x ×Ñ f = Wº - éëf - f eq ù
û
t
Key Point:
j
Integral solution:
j+1
1 t eq
f j + 1/ 2 (t ) = ò f (x ¢, t ¢)e - (t - t ¢)/ t dt ¢+ e - (t - tn )/ t f (x j + 1/ 2 - x(t - t n ),t n )
t tn
Free transport
Equilibrium
After some algebraic, the above solution can be approximated as
f j + 1/ 2 (t ) = (1 - e- t / t
)
fˆjeq+ 1/ 2 (tn ) + e- t / t fˆj + 1/2 (tn )
Dt ? t
fˆjeq+ 1/ 2 (t n )
Chapman-Enskog expansion
Dt = t
fˆj + 1/ 2 (t n )
Free-transport
DUGKS vs UGKS
(a)Common:
• Finite-volume formulation;
• AP property;
• collision-transport coupling
(b) Differences: in DUGKS
• W are slave variables and are not required to
update simultaneously with f
• Using a discrete (characteristic) solution instead
of integral solution in the construction of cellinterface distribution function
# Comparison with Finite-Volume LBM
ci
Lattice Boltzmann method (LBM)
¶ t fi + c i ×Ñ fi = -
fi
eq
1
[fi - fieq ]
t
é c i ×u (c i ×u )2 u ×u
ù
= wi r ê1 +
+
+
L
ú
RT
2(RT )2 2RT
ú
ëê
û
Standard LBM: time-splitting scheme
Collision
¶ t fi = -
1
[fi - fieq ]
t
fi ¢(x , t ) = fi (x , t ) -
Free transport ¶ t fi + c i ×Ñ fi = 0
Evolution equation:
Viscosity:
fi (x + c i D t , t + D t ) = fi ¢(x , t )
fi (x + c i D t , t + D t ) - fi (x , t ) = -
(
rn = t -
Dt
p
2
)
Dt
[fi (x , t ) - fieq (x , t ) ]
t
Dt
[fi (x ,t ) - fieq (x ,t ) ]
t
Numerical viscosity is absorbed into the physical one
Limitations:
1. Regular lattice
2. Low Mach incompressible flows
# Comparison with Finite-Volume LBM
Finite-volume LBM (Peng et al, PRE 1999; Succi et al, PCFD 2005; )
1
[fi - fieq ]
t
¶ t fi + c i ×Ñ fi = -
f j ,i (t n + 1 ) - f j ,i (t n ) +
j+1/2
D t éˆ
ˆj - 1/ 2,i (t n ) ù= - 1 [fi (t n ) - fieq (t n ) ]
f
(
t
)
f
n
j
+
1
/
2
,
i
ú
û
D x ëê
t
j
Micro-flux is reconstructed without considering collision effects
Viscosity:
rn = t p
Numerical dissipation cannot be absorbed
Limitations (Succi, PCFD, 2005):
1. time step is limited by collision time
D t < 2t
Difference between DUGKS and FV-LBM:
2. Large numerical dissipation
DUGKS is AP, but FV-LBM not
j+1
Outline
• Motivation
• Formulation and properties
• Numerical results
• Summary
Test cases
• 1D shock wave structure
• 1D shock tube
• 2D cavity flow
Collision model:
f
eq
= f
S
= f
M
é
c ×q
ê1 + (1 - P r)
êë
5pR T
æ c2
öù
÷ú
çç
÷ú
çè R T - 5 ÷
ø
û
Shakhov model
1D shock wave structure
Parameters: Pr=2/3,  = 5/3, Tw
Left: Density and velocity profiles; Right: heat flux and stress (Ma=1.2)
DUGKS agree with UGKS excellently
Again, DUGKS agree with UGKS
excellently
DUGKS as a shock capturing scheme
Density (Left) and Temperature
(Right) profile with different grid
resolutions (Ma=1.2, CFL=0.95)
1D shock tube problem
Parameters: Pr=0.72,  = 1.4, T0.5
(r , u , p )L = (1.0, 0.0, 1.0)
(r , u , p )R = (0.125, 0.0, 0.1)
Domain: 0  x  1;
Mesh: 100 cell, uniform
Discrete velocity : 200 uniform gird in [-10 10]
Reference mean free path
By changing the reference viscosity at left boundary, the flow can
changes from continuum to free-molecular flows
m0 = 10, 1, 0.1, 10- 3 , 10- 5
l 0 = 12.77, 1.277, 0.1277, 1.277 ´ 10- 3 , 1.277 ´ 10- 5
1
UGKS
present
2.2
0.8
2

0.6
0.4
1.8
0.2
1.6
0
0
0.5

=10: Free-molecular flow
1
1.4
0
1
0.8
0.5
1
0.5
1

UGKS
present
0.6


UGKS
present
0.4
0.2
0
0
UGKS
present
1
UGKS
present
2.2
2

0.6
0.4
1.8
0.2
1.6
0
0
0.5

1
1.4
0
=1: transition flow
0.5
1
0.5
1

1
0.8
UGKS
present
0.6


0.8
0.4
0.2
0
0
UGKS
present
1
UGKS
present
2.2
2

0.6
0.4
1.8
0.2
1.6
0
0
0.5

1
1.4
0
=0.1: low transition flow
0.5

1
1
0.8
UGKS
present
0.6


0.8
0.4
0.2
0
0
0.5
1
1
2.2
UGKS
present
0.8
2

0.6
1.8
0.4
1.6
0.2
0
0
0.5

1
1.4
0
0.8
1
0.5
1

1
=0.001: slip flow
0.5
UGKS
present
0.6


UGKS
present
0.4
0.2
0
0

1
2.2
0.8
UGKS
present

2
0.6
1.8
0.4
1.6
0.2
0
0
0.5

1
1.4
0
0.5
1
0.5
1

1
=1.0e-5: continuum flow
0.8
UGKS
present
0.6


UGKS
present
0.4
0.2
0
0
2D Cavity Flow
Parameters: Pr=2/3,  = 5/3, T0.81
Domain: 0  x, y  1;
Mesh: 60x60 cell, uniform
Discrete velocity : 28x28 Gauss-Hermite
Kn=0.075
Temperature.
White and background: DSMC
Black Dashed: DUGKS
Kn=0.075
Heat Flux
Kn=0.075
Velocity
UGKS: Huang, Xu, and Yu, CiCP 12 (2012)
Present DUGKS
Temperature and Heat Flux
Kn=1.44e-3; Re=100
Comparison with LBM
Stability:
Re=1000
LBM becomes unstable on 64 x 64 uniform mesh
UGKS is still stable on 20 x 20 uniform mesh
80 x 80 uniform mesh
LBM becomes unstable as Re=1195
UGKS is still stable as Re=4000 (CFL=0.95)
Velocity
DUGKS
LBM
DUGKS
LBM
Pressure fields
Summary
• The DUGKS method has the nice AP property
• The DUGKS provides a potential tool for compressible flows in
different regimes
Thank you for your attention!