Transcript Gas-kinetic united algorithm for re

Sino-German Symposium on Modern Numerical Methods
for Compressible Fluid Flows and Related Problems
Gas-Kinetic Unified Algorithm for Reentry
Hypersonic Flows Covering Various Flow Regimes
Solving Boltzmann Model Equation in
Thermodynamic Nonequilibrium
Zhi-Hui Li
Hypervelocity Aerodynamics Institute , National Lab. for CFD
( China Aerodynamics Research & Development Center )
May 21--26, 2014
Beijing  China
Outline
1 Introduction
2 Unified Boltzmann Model Equation in Nonequilibrium
3 Development of Discrete Velocity Ordinate Method
4 Construct Gas-Kinetic Numerical Scheme for Solving
Velocity Distribution Function
5 Development of Discrete Velocity Quadrature
Methods for Macroscopic Flow Variables
6 Gas-Kinetic Boundary Conditions and Numerical
Procedures for the Velocity Distribution Function
7 Gas-Kinetic Parallel Algorithm for 3D Complex Flows
8 Numerical Study of Three-dimensional Complex Flows
Covering Various Flow Regimes
9 Concluding Remarks
1. INTRODUCTION
Kn   D
Continuum Slip flow
Rarefied transition
Free molecule
 To study the aerodynamics from various flow regimes,
the traditional way is to deal with different methods.
 Rarefied flow: DSMC
 Continuum flow: Euler, Navier-Stokes
 Two methods are totally different, and the computed
results are difficult to be linked up smoothly with
altitude.
 Problem: Engineering development of current or intending
spaceflight projects is closely concerned with complex
aerothermodynamics of hypersonic flows in the intermediate
range of Knudsen numbers, especially in the rarefied transition
and in the near-continuum flow regimes.
 Challenge: How to solve multi-scale non-equilibrium flows
over the whole flow regimes during spacecraft re-entry?
 Boltzmann Equation(BE): describe molecular transport
phenomena from full spectrum of flow regimes. The
difficulties encountered in solving BE are associated with
nonlinear multidimensional integral nature of collision term
 f
 

   4
' '
 Vx
 Vy
 Vz
 F          0 ( f f1  ff 1 ) V  V1   d  d V1
t
x
y
z
V
f
f
f
f
 From the kinetic theory of gases, numerous statistical or
relaxation kinetic model equations resembling to the original BE
have been put forward
 BGK, ES, Shakhov, polynomial, hierarchy kinetic, and McCormack
 Enlighten: Solve the nonlinear kinetic models simplified by BE
and probably finds a more economical and efficient approach to
deal with complex gas flows.
 BGK equation provides an effective and tractable ways:
 f
 V     BGK ( f  f M )
t
r
f
f M  n /  2  RT
3 2 exp[  c 2
( 2 RT ) ]
 Way 1: Lattice Boltzmann Method (LBM)
 Applying Lattice Gas Cellular Automata (LGCA) and DVM, LBM
methods have been developed for solving fluid dynamic problems,
such as continuum or near-continuum flow in low speed,
turbulence, microflow or porous medium models.
 Frisch, Pomean, Nie X, Doolen, Succi, Lee, Qian, Chen, Luo, Yong,
Guo, Yan, Zhong, Wang etc.
 Way 2: Gas-kinetic KFVS-, BGK-type schemes
 The Maxwellian distribution function is translated into
macroscopic flow variables in equilibrium, some gas-kinetic
methods are developed to solve inviscid gas dynamics.
 Beam, Pullin, Macrossan, Chen etc.: KFVS-type
 Applying the asymptotic expansion of velocity distribution
function to Maxwellian distribution on flux conservation at
cell interface, BGK-type schemes have been presented, such
as BGK-Euler,BGK-NS,BGK-Burnett,BGK-Super-Burnett.
 Prendergast, Kun Xu, Kim C, Tang, Li, Zhong etc.: BGK-type
 Way 3: Gas-kinetic numerical algorithm by directly solving
the velocity distribution function
 Applying discrete ordinate technique and reduced
velocity distribution functions, finite difference explicit
and implicit methods, and discrete-velocity models of
conservation and dissipation of entropy for solving
one- and two-dimensional BGK-Boltzmann model
equations have been set forth for high Mach flows.
 Chu, Shakhov, Morinishi, Chung, Yang, Aoki, Tcheremissine,
Mieussens, Kolobov, Aristov , Titarev etc.
 Specially, Gas-Kinetic Unified Algorithm (GKUA) since 1999:
 Boltzmann model equation can be described by modifying the
BGK model from two sides on molecular collision relaxing
parameter and local equilibrium distribution function
 Discrete velocity ordinate method (DVOM) and numerical
integration techniques are developed to dynamically evaluate
macroscopic flow variables
 Gas-kinetic numerical schemes are constructed to directly capture
the evolution and update of velocity distribution function
 Unified gas-surface interaction model is token by the velocity
distribution function
 The GKUA has presented and applied from low speed to
extremely hypersonic flows for MEMS and spacecraft reentry,see
IJNMF2003; JCP2004; JCP2009; CMA2011; Adv. Space. Tech. 2011.
 Recently, a unified gas-kinetic scheme is developed from the
combination of BGK and DVOM by Xu and Huang, Chen etc.
 How to solve BE models involving thermodynamics
non-equilibrium effect ?
 Relaxation kinetic models involving rotational degrees of
freedom
Morse, Rykov etc.
 Inelastic relaxing phenomena model is presented with
experience treatment and the expansion of ChapmanEnskog with small disturbance
 C.S.Wang-Chang, G.E.Uhlenbeck etc.
 Rykov model is applied to simulate hypersonic flow around plate
Titarev etc.
 This study is aimed at extending the GKUA to solve BE
models involving thermodynamics non-equilibrium effect
for possible engineering applications to spacecraft re-entry.
2. Unified Boltzmann Model Equation in
Nonequilibrium Effect
 Based on Rykov model, relaxation effect of rotational degrees
of freedom is considered into the evolution and update of VDF
 The spin movement of diatomic molecule is described by
moment of inertia, and the conservation of total angle
momentum is taken as a new Boltzmann collision invariant



 
 
K  m r  V  M  m r  V  J




 Internal energy is equally distributed in each degree of
freedom by introducing energy model partition function.


1  mc 2

kT  kT t     
 e  fded V
2
2
n  2

5
3
 Internal energy is taken as the independent variables of VDF
2
f
(0)
(t , r ,  , M )  n (
m
2  kT
3
)
2
 m (  U ) 2 
M
1
exp  
)
 (4  JkT ) exp( 
2
kT
2
JkT


 Boltzmann Model Eq. in non-equilibrium effect is presented in
the framework of GKUA covering various flow regimes.
 f
 V    v r ( f r  f )  vt ( f t  f )
t
r
f
vr 
Pt
t

1
vt 
Z
f  N ( e )(
m
t
f
r
 n(
2 kTt
m
2 kT
)
3 2
)
Pt
t
3 2
1
 (1 
Pt
)
t
Z

4 ( 5  2  ( 7  2  )

5 (  1)(   2 )
R
2
 1 2

T
n

1


n
T  1
mc i 5 mc
q mc i (  e )
exp( 
)[1 

( 
)  (1   ) i

]
2 kTt
15 nkT t kTt 2 2 kTt
kTt

mc
mc
2
2
2
t
2
qi
t
qi
2
r
mc i 5 mc
q mc i
e
exp( 
)( kT ) exp( 
)[ 1  

( 
) i
 (1 
)]
2 kT
kT
P kT 2 2 kT
PkT
kT
1
e
r
 All macroscopic flow variables are determined by moments of
the distribution function over the velocity space.

   
n  r , t           f  dedV x dV y dV z

   
nU i ( r , t )          V i f  dedV x dV y dV z

2

m
c
   
nkT t ( r , t )         
f  dedV x dV y dV z
2
2

1    
kT r ( r , t )          ef  dedV x dV y dV z
n

   
 ij ( r , t )          mc i c j f  dedV x dV y dV z  P  ij
3

m    
q it ( r , t ) 
        c 2 c i f  dedV x dV y dV z
2

   
q ir ( r , t )          c i efdedV x dV y dV z
q i  q it  q ir
 Focus:
f
5
2
T 
3
2
Tt  T r
How to solve? Seven independent variables

V

、r e、t
 The energy VDF is integrated by the weight factors of 1 and e
on the internal energy. Two energy-level reduced distribution
functions are introduced to remove the continuous dependence
of Boltzmann Model Eq. on internal energy:
 
f 0 ( t , r , V )   fde
 
f1 ( t , r , V )   efde
 The unified and reduced VDF equations in non-equilibrium
effect is obtained for various flow regimes.
f 0
 f 0
r
t

V


v
(
f

f
)

v
(
f
 f0 )
i
r
0
0
t
0
  t
xi

 f1
 f1
r
t

 Vi 
 v r ( f1  f1 )  v t ( f1  f1 )
  t
xi
 All flow variables are evaluated and updated by the reduced
non-equilibrium VDFs of f 0 、f1 over the velocity space.
3. Development of Discrete Velocity Ordinate Method
 The VDFs f 0 f1 remain with probability density distribution on the
principle of probability statistics, VDF possesses the property of
exponential function exp(-c2), not Maxwellian distribution.
 VDF is confined to the finite region
0 .8
-1 0
0 .6
0 .7 5
-5
0 .5 5
0
0 .5
10
0 .4 5
0 .4
0 .3 5
0 .3
0 .2 5
0 .2
0 .1 5
0 .1
0 .0 5
M s  3 .8
-2 0
0 .7
-5
0 .6
0
0 .5 5
10
0 .5
Ms 8
-3 0
-1 0
0 .6 5
 
0 .5 5
R e duce d distribution: g(x,V x )
M s  1 .4
-2 0
0 .6 5
0 .6
 
0 .8
R e duce d distribution: g(x,V x )
0 .7
R e duce d distribution: g(x,V x )
0 .8 5
 
0 .7 5

V
0 .4 5
0 .4
0 .3 5
0 .3
0 .2 5
0 .2
0 .1 5
0 .1
0 .5
-5
0 .4 5
0
0 .4
5
10
0 .3 5
0 .3
0 .2 5
0 .2
0 .1 5
0 .1
0 .0 5
0 .0 5
0
0
0
-0 .0 5
-3
-2
-1
0
1
V x /C m re f
2
3
4
-0 .0 5
-5
-4
-3
-2
-1
0
1
V x /C m re f
2
3
4
5
6
7
-0 .0 5
-1 0
-5
0
5
10
15
V x /C m re f
Bimodal
distribution
 Discrete Velocity Ordinate Method (DVOM) can be developed

to discretize the finite velocity region removed from U and to
replace continuous dependency of VDF on velocity space.
 The selection of DVO points is optimized and corresponded
with the evaluation points and weight coefficients of the
integration rule in a way that the approximation is exact.

b
aW
N
(V ) f (V ) dV   W f (V )
 1
 Boltzmann’s H-theorem and conservation condition are
guaranteed at each of DVO points with self-adaption.



b
N
(m )


V
2

1 
V
f

f

d
V

0
a


(1 )
(2)
(3)
2
 The VDF Eq. is transformed into hyperbolic conservation laws
with nonlinear source terms at each of (V  , V  , V  )
x
Q
t

F
x
x

 f 0  , , ( t , x , y , z ) 
Q  J 

 f1 , , ( t , x , y , z ) 
F
y
 V y Q
F
y

y
y
F
z
z
z
F
x
 V x Q
F
z
 V z Q
vr ( f r
 f 0  , , )  v t (
0  , ,
S  J 
r
 v r ( f1 , ,  f1 , , )  v t (
 S
f 0t , ,  f 0 , , ) 

t
f1 , ,  f1 , , ) 

4. Construct Gas-Kinetic Numerical Scheme for
Solving Velocity Distribution Function
 The finite-difference method from CFD can be extended to
directly solve the discrete VDFs.
U
t

F

 f 0  , , ( t , x , y , z ) 
U  J 

f
(
t
,
x
,
y
,
z
)
 1 , ,


G

F UU

H

G  VU
 S
H WU
t
v ( f r


f
)

v
(
f
r
0  ,  ,
0  ,  ,
t
0  ,  ,  f 0  ,  , )
S  J 

r
t
v
(
f

f
)

v
(
f

f
)
 r 1 , ,
1 ,  ,
t
1 ,  ,
1 ,  , 

U  V x  x  V y  y  V z  z
V  V x  x  V y  y  V z  z
W  V x  x  V y  y  V z  z
J   ( x , y , z )  (  , ,  )
 The finite difference second-order scheme for solving the
discrete velocity distribution functions are constructed as
U
n 1
 LS (
t
2
) L (
t
2
) L (
t
2
) L (  t ) L (
t
2
) L (
t
2
) LS (
t
2
 This can be split into four steps.
U
*
U
**
U
***
 L s (  t )U
n
U
n
  t (1    t / 2 ) S
 L (  t )U
 L (  t )U
***
2
2

b t
 1  b  t   
2
2


U

**

 ***
a t
 1  a  t   
  2 U
2


2
U
2
*
**
2
U
t

 *
c t
 L (  t )U  1  c  t   
  2 U
2


2
n 1
n
U
t
U
t
U
t
 S



H

G

F

0
0
0
)U
n
 Second-order Runge-Kutta method solves non-linear colliding
relaxation source term:



 tU
*
  1   t   S (U )
2


n
U U   t   t U
*
 tU
U
n
*



*
  1   t   S (U )
2


**
n 1
U 
n
t
2
( t U
*
  tU
**
)
 NND scheme with primitive variables solves convective term:
U
U
*
U
n 1

n
1

[U
t

n
n
n
(Q i  1 / 2  Q i 1 / 2 )
U
*
2

t
*
*
( Q i  1 / 2  Q i  1 / 2 )]



Q i  1 / 2  H i  1 / 2 (U L )  H i  1 / 2 (U R )
U L ,i 1 / 2  U
U R ,i 1 / 2  U
p ,i

p ,i 1
1
min mod(  U
2

1
2
p , i 1 / 2 ,  U p , i  1 / 2 )
min mod(  U
p ,i 1 / 2 ,  U p ,i  3 / 2 )
5．Development of Discrete Velocity Quadrature
Methods for Macroscopic Flow Variables
 Once discrete VDFs are solved, macroscopic flow moments in
the physical space are updated by the appropriate discrete
velocity quadrature method.
 The new Gauss-type integration methods and self-adaptive
technique are presented to simulate hypersonic flows.
 Bell-type Gauss quadrature formula:

0
2

1/ 2
N
2
exp( V ) f (V ) dV   W  f (V  )
 1
W (V ) 
2

1/ 2
2
exp( V )
p  (V )  (V  b ) p  1 (V )  g  p   2 (V )
 The macroscopic flow variables can be evaluated by the
integrating summation with the weight function
6. Gas-Kinetic Boundary Conditions and Numerical
Procedures for the Velocity Distribution Function
 The gas-kinetic algorithm is based on time evolution of the
VDF, the interaction of gas/surface and aerothermodynamics
are expressed by directly acting on the VDF.
 Escape the statistical fluctuation of DSMC.
 Obviate the difficulties in expressing rarefied effect by
macroscopic N-S solvers.
 Molecules hitting surface must be reflected back to the gas,
the reflected VDFs are
f 0w , ,
 (V x   U w ) 2  (V y   V w ) 2  (V z   W w ) 2 

exp  

( T at ) 3 / 2
T at


nw
f 1w  f 0 w  T a
r
t
T 
t
a
Ei
Ni
t
 ( t   )(
r
t
Ei
Ni
 Tw ) 
1
2
r
 (2
t
r
Ei
Ni
 Tw )
 
c n  0
r
T 2
r
a
Ei
Ni
r
 ( r   )( 2
t
r
Ei
Ni
t
 T w )  2 (
r
t
Ei
Ni
 Tw )
 The number density n w of molecules diffusing from the surface
are determined from mass balance condition on the surface.

w
c n  0 c n f M dV
nw
( T w ) 3 2
nw  2
 
c n  0

 c  0  c n f 0 dV
n
c n  0 c n  e


V x2  V y2  V z2
Tw
dV x dV y dV z  c  0  c n  f 0 dV x dV y dV z
n
  c n  0 c n  f 0 dV x dV y dV z
t
Ta

cn  (cn  cn ) / 2
 If
, molecules don’t strike on the surface, the discrete
VDFs at wall cells are solved by using second-order upwinddifference approximations from adjacent grids in flow field.
7． Gas-Kinetic Parallel Algorithm for
Three-Dimensional Complex Flows
 For three-dimensional flow, GKUA computing space relates to
discrete velocity, physical and energy-level multi-dimensional
space, and the GKUA requires six-dimensional arrays to access
discrete VDFs at all discrete ordinate points.
① Parallel domain decomposition of discrete velocity space
 Computing space decompose as physical space
velocity space  ( ,  ,  )
 r (i , j , k )
and
V
   r  V
Np
i  
   i
i 1
j
 
(i  j )
 Data from sub-space  map to corresponding processors
decompose  as N subspace  in block-layout manner
i
V
p
 i   r   vi
Vi
Np
 V    Vi
i 1
PE i
 Variable dependency relations of GKUA: For  domain
decomposition, complete parallelization without data
communication arise from velocity space during solving the
discrete VDFs. Data communication arise in the sum-reduction
computation  V during evaluating macroscopic flow variables.
② Data communication analysis
V
 Evaluate data traffic by two fork tree parallel reduction:
C V  14  N
p
N
p
N
p
NiN jNk
 domain decomposition is carried by three-dimensional way:
V
N p N p N p
 V      V , , ,
 1  1  1
  , ,   V , , ,   r
 Variables   , , map to processors PE
accordance with processor arrays
i, j ,k
,distribute N
p
N
p
N
p
in
C V  14  V 
N
p

N
N
p

N
N
p
V  N i N j N k N N N
N
 To get smaller C V ,take N  N , N  N , N  N
 For spacecraft reentry with high Mach numbers,
decomposition strategy  V is suitable for many DVO points
p
③

p

p

parallel scalability
The computation load produced from parallel reduction
sum of discrete velocity numerical integration only
account for 1 5 of total workload of the algorithm.
The number of parallel processors can reach up to the
maximum N  N  N so that the parallel scalability can
be effectively enlarged.
p
p
p
 For spacecraft simulation with Mach 25, the DVO points is up to
N   Nalgorithm
50  4solve
 10
, the
the BE model can realize high  N   100  80  to
performance
computation with thousands and tens of
thousands, even more massive parallelism.
2
S pe e d u p
S pe e d up R a tio
5
Id e a l S p e e d u p
1 .9
R ea l S pe e dup
1 .8
16
1
Id e a l S p e e d u p
15
R e al S pe e dup
14
13
0 .9
12
1 .7
P a ra lle l e fficie n cy
11
1 .6
10
9
1 .5
8
1 .4
7
6
1 .3
5
256~512CPU
1 .2
64~1024CPU
4
0 .8
0 .7
0 .6
3
1 .1
2
1
256
320
384
448
512
1
128
256
N u m be r o f P ro ce sso rs
Speed-up ratio
384
512
640
768
896
1024
N u m b e r of P ro ce ssors
0 .5
0
8
16
24
P ro ce sso rs
Parallel efficiency
 Parallel speed-up goes up as near-linearity with high
parallel efficiency.
32
1024~8000CPU
4950~20625CPU
 The gas-kinetic parallel algorithm possesses quite high parallel
efficiency and expansibility with good load balance, which
makes it possible to solve three-dimensional complex hypersonic
problems covering various flow regimes.
8．Numerical Study of Three-dimensional Complex
Flows Covering Various Flow Regimes
① Inner flows of shock wave structures in nitrogen
M
M s  1 . 53
M s  3 .2
M s  10
s
 1 . 53 ~ 25
M s  25
 Computed profiles agree with experimental data and solutions of GBE well
M s  10
② Sphere flows in perfect gas and nonequilibrium effect
Kn   0 . 0126
Ma   2 Re   236 . 6
T  273 . 15 K
Tw T  1 . 8 Tw T0  1
 The computation spends 47 seconds/step for perfect gas, however 78s/step for
non-equilibrium model, the computational load and memory increases two times
 The front shock for thermodynamic non-equilibrium effect is closer 14% to the
stagnation surface than that of perfect gas
③ Bicone reentry flows from high rarefied to continuum
flow regimes
Kn  1
Ma  3
Kn  1
Ma  3
76  42  61
Kn  0 . 1
Kn  0 . 0001
Stagnation line density profiles
 Computed results match solutions of GBE and theoretic analysis well
H  95 km
 Density and pressure contours around tine bicone with different Kn=1, 0.1
H  83 km
H  95 km
 Mach number contours and streamlines around tine bicone with Kn=1, 0.1
H  83 km
Ma  3
Kn  1
 Effect of thermodynamic non-equilibrium on translational and rotational
temperature contours around tine bicone with Kn=1, 0.1
Ma  3
Kn  0 . 1
Ma  3
Kn  1
 Overall temperature and heat flux contours around tine bicone with Kn=1, 0.1
Ma  3
Kn  0 . 1
H  32 . 6 km
Continuum flow around tine bicone with Kn=0.0001
Kn   0 . 00065
H  45 km
76  42  61
Re   26179
Ma   11 . 39
Continuum flow around tine bicone with Kn=0.0001
H  44 km
Ma   9 . 59
Kn   0 . 000564
Re   25684
9. Concluding Remarks
 The GKUA for re-entry hypersonic flows of spacecraft
has been developed for the whole range of flow regimes
by solving Boltzmann model equation involving nonequilibrium effect.
 The computations of hypersonic flows and aerodynamic
phenomena around sphere, double-cone and spacecraft
covering various flow regimes have indicated both high
resolution of the flow fields and good agreement with the
relevant theoretical, DSMC and experimental results.
 The GKUA provides a feasible way to simulate spacecraft
re-entry hypersonic aerothermodynamics with highperformance massively parallel computation.
 Further investigation on cargo re-entry capsule and large
spacecraft involving real-gas effect with internal energy.
Acknowledgements
 This work are supported by NSFC under Grants Nos.
91130018, 11325212, and National Key Basic Research
Program (2014CB744100)
 Massively parallel computations are run by National
Supercomputer Center in Tianjin and Jinan, and
National Parallel Computing Center in Beijing
 Joint work with my postgraduates Aoping Peng, Junlin
Wu, Xinyu Jiang, Ming Fang and Qiang Ma.
• Thanks for organizing committee and Sino-German
Center, specially to Prof. Jiequan Li and Song Jiang etc.
Thank you for your kind attention!
Wish you have a good time in Beijing！