La ice Bol zmann lux sol e o simula ion o iscous hype sonic lows
V In e na ional Con e ence on Pa icle-based Me hods – Fundamen als and Applica ions
PARTICLES 2017
P. W igge s, M. Bischo , E. Oña e, D.R.J. Owen, & T. Zohdi (Eds)
LATTICE BOLTZMANN FLUX SOLVER FOR SIMULATION OF
HYPERSONIC FLOWS
Z. X. MENG¹, C. SHU², L. M. YANG3, W. H. ZHANG 4, F. HU5 AND S. Z. LI6
¹ College o Ae ospace Science and Enginee ing, Na ional Uni e si y o De ense Technology
Changsha, 410073, China
Co esponding au ho : [email protected]
² Depa men o Mechanical Enginee ing, Na ional Uni e si y o Singapo e
10 Ken Ridge C escen , Singapo e 119260, Singapo e
m[email p o ec ed]
3 Depa men o Mechanical Enginee ing, Na ional Uni e si y o Singapo e
10 Ken Ridge C escen , Singapo e 17576, Singapo e
yang[email p o ec ed]m
4 College o Ae ospace Science and Enginee ing, Na ional Uni e si y o De ense Technology
Changsha, 410073, China
zhang[email p o ec ed]
5 College o Ae ospace Science and Enginee ing, Na ional Uni e si y o De ense Technology
Changsha, 410073, China
hu [email protected]
6 College o Ae ospace Science and Enginee ing, Na ional Uni e si y o De ense Technology
Changsha, 410073, China
[email protected]
Key wo ds: Non- ee pa ame e D1Q4 model, la ice Bol zmann lux sol e , hype sonic lows,
ini e olume me hod
Abs ac : In his pape , a s able La ice Bol zmann Flux Sol e (LBFS) is p oposed o
simula ion o hype sonic lows. In LBFS, he ini e olume me hod is applied o sol e he
Na ie -S okes equa ions. One-dimensional La ice Bol zmann model is applied o econs uc
he in iscid lux ac oss he cell in e ace, while he iscous lux is sol ed by con en ional
smoo h app oxima ion unc ion. The p esen wo k ex ends he exis ing LBFS o calcula e
hype sonic low ield on he leewa d, which is ha d o ge con e gen esul s due o ex emely
low p essu e e ec s in his a ea. Simula ion o a biconics model is s udied. I is disco e ed ha
he ail a ea o double cone is ela ed o he maximum Mach numbe ha could be con e gen .
The la ge he diame e o ail a ea is, he smalle Mach numbe could be con e gen . Hence,
he low p essu e a ea behind double cone ail will ha e la ge e ec s du ing he LBFS simula ion
o hype sonic low. Two measu emen s a e applied in his pape o o e come he low p essu e
p oblem. The i s one is o apply a local block g id e inemen me hod based on he low
condi ions o imp o ing he s abili y. The second is o add a cons ain pa ame e o elimina e
nega i e alue and gi e ou a p ope one. Hence, LBFS is able o ge con e gen esul o he
hype sonic low ield on bo h windwa d and leewa d. Se e al nume ical examples a e es ed o
compa e he pe o mance o me hod p esen ed in his pape . Simula ion esul s show ha
me hod p esen in his pape is able o calcula e hype sonic low ield on he leewa d wi h bo h
ine accu a e and e icien .
203
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
2
1 INTRODUCTION
The hype sonic low ield is a esea ch ocus wi h i s special cha ac e is ics [1]: (1) s ong
shock e ec will cause ie ce comp ession a e he shock wa e, (2) iscosi y g ea ly changed
om wall su ace o a ield a ea, (3) o hype sonic ehicles lying in high al i ude,
con inuous medium assump ion is no longe app op ia e because o low densi y e ec . The
compu a ional luid dynamics (CFD) is becoming mo e and mo e popula in simula ing
hype sonic low ield. The ini e olume me hod (FVM) is widely used in exis ing nume ical
me hods [2]. This la gely a ibu es o i s spa ial disc e iza ion is ca ied ou di ec ly in he
physical space. Hence, s uc u al g id is sui able o disc e iza ion, which p o ides an e ec i e
way o sol ing complex geome y p oblems. In he p ocess o sol ing N-S equa ions by FVM,
lux sol e is he key o e alua e iscous and in icid luxes a cell in e ace. The iscous lux
is e alua ed by applying a smoo h unc ion app oxima ion and in iscid lux is calcula ed by
a ious upwind schemes such as Roe scheme [3], an Lee scheme [4] and AUSM (Ad ec ion
Ups eam Spli ing Me hod) scheme [5].
Bol zmann equa ion-based me hod is an al e na i e app oach o simula ing comp essible
low ield, including DVBE [6][8] (disc e e eloci y Bol zmann equa ion) and gas-kine ic
scheme [9]-[13]. Gas-kine ic scheme ha e be e e iciency han DVBE me hod. Howe e , bo h
o hem a e less e icien and mo e complica ed han con en ional Na ie - S ocks slo e s.
La ice Bol zmann Flux Sol e (LBFS) is a mo e e icien me hod which is p oposed by Ji
e al. [14] and o simula ing in iscid lows. LBFS is imp o ed by Shu and Yang e al. [15]-
[19] and has been p o ed ca ching s ong shock wa e and expansion wa e e y well. In LBFS,
Eule equa ions a e disc e ized by FVM and he in iscid lux a cell in e ace is econs uc ed
by local solu ion o 1-D comp essible La ice Bol zmann model. LBFS is mo e e icien and
easy o apply han DVBE and gas-kine ic me hod due o only 1-D La ice Bol zmann model is
applied and mac oscopic go e ning equa ions a e sol ed. Howe e , because o 1-D model is
applied along no mal di ec ion a cell in e ace, he angen ial e ec canno be p ope ly
conside ed. Viscous lux should be aken in o accoun o sol e iscous p oblems. F om
chapman-Enskog expansion analysis [19], he in iscid lux can be ully de e mined by
equilib ium non-equilib ium dis ibu ion unc ion a cell in e ace, and he non-equilib ium pa
can be ea ed as nume ical dissipa ion. Fo comp essible iscous lows, especially o
hype sonic lows, he nume ical dissipa ion should be con olled. A swi ch unc ion is p oposed
by Yang e al. [20] o weigh iscous nume ical dissipa ion and has been p o ed pe o med well
in iscous simula ion. Fo hype sonic low ield, physical alues such as p essu e in leewa d
side can be ex emely low. Hence i is e y easy o ge a nega i e alue in calcula ing p ocess,
which make i e y di icul o ge a con e gen esul .
In his wo k, a s able LBFS is p oposed o simula ion o 2-D comp essible hype sonic
iscous lows. The in iscid lux is calcula ed by LBFS and iscous lux is compu ed by smoo h
unc ion app oxima ion. Mul iblock g ids and local g id e inemen me hod a e applied o
imp o e calcula ion s abili y. Also a cons ain pa ame e is added o elimina e nega i e alue
and gi e ou a p ope one in calcula ing p ocess. Besides, he implici LU-SGS me hod is used o
speed up he con e gence a e. In he end, a biconics model in Ma=9.86 low ield is simula ed o
alida e he de eloped sol e . The esul s o sol e p esen ed in his wo k is compa ed wi h an
Lee and Roe scheme o e alua e he p ecision p oposed in his a icle.
204
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
3
2 METHODOLOGY
In LBFS, Na ie -S ockes equa ions a e sol ed in mac oscopic scale and local solu ion o
La ice Bol zmann equa ions a e applied in o cons uc in iscid lux sol e a he cell in e ace.
I has been p o ed ha pa icle po en ial ene gy is independen om la ice eloci y [15].
La ice eloci y would be decided by highe momen um conse a ion ela ion in his wo k
ins ead o by a i icial selec ion in adi ional me hod.
N-S equa ions o in e g al o m wi hou sou ce e m can be w i en as ollowing.
( )0
c
d ds

  


W FF
(1)
In which he conse a i e a iables
W
, in iscid lux
c
F
and iscous lux
F
a e gi en by Eq. (2)
in 2-D low ield.
u
E










W
,
()
n
nx
c
ny
n
U
uU n p
U n p
E pU














F
,
0
x xx y xy
x yx y yy
xx yy
nn
nn
nn









 


F
(2)
The exp essions o no mal eloci y
n
U
and o al ene gy
E
a e shown as Eq. (3). Whe e
/ [( 1) ]ep

 
is he po en ial ene gy o mean low.
22
1()
2
nx y
U nu n
Ee u
 

(3)
The in eg al o luxes in Eq. (1) can be disc e ized by FVM and he app oxima ed summa ion
o m can be w i en as
1
1()
N
I
ci i i
i
I
dF FS
d 



W
(4)
In which I ep esen s he con ol olume index,
I

is he olume and
N
is numbe o aces o
con ol olume I. In his pape , he in iscid lux
c
F
in Eq. (4) will be sol ed by LBFS wi h non-
ee pa ame e D1Q4 model [18] and iscous lux
F
will be sol ed by cen al di e ence
me hod.
2.1 Non- ee pa ame e D1Q4 la ice Bol zmann model
The dis ibu ion o disc e e la ice eloci ies o D1Q4 model is shown in Fig. 1. This model
con ains 4 equilib ium dis ibu ion unc ions
1
eq
,
2
eq
,
3
eq
,
4
eq
and 2 la ice eloci ies
1
d
,
2
d
(shown as Eq. (5)). The de i a ion p ocess de ails a e shown in e e ences [18] .
2
43
-d
1
-d
2
d
1
d
2
1
205
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
4
Figu e 1: Dis ibu ion o disc e e la ice eloci ies o D1Q4 model
2 2 2 23 2
12 2 1 1
122
11 2
2 2 2 23 2
12 2 1 1
222
11 2
2 2 2 23 2
12 1 2 1
322
21 2
2 2 2 23 2
12 1 2 2
422
21 2
2
1
( 3)
2( )
( 3)
2( )
( 3)
2( )
( 3)
2( )
eq
eq
eq
eq
d d d u d u d c u uc
dd d
d d d u d u d c u uc
dd d
d d d u d u d c u uc
dd d
d d d u d u d c u uc
dd d
du




    

    

   

   

  2 22 4
2 2 22 4
2
34 6
34 6
c uc c
d u c uc c

  
(5)
In which
/c Dp


ep esen s pa icula eloci y o pa icles and
D
is space dimension. I has
been p o ed by Yang [18] ha physical conse a ion laws (Eq. (6)) can e i i y he N-S
equa ions by applying ela ions in Eq. (5).
i

means pa icle eloci y in i-di ec ion, o example
1 12 13 2
,,d dd
 
 
and
42
d

 
. Fo highe dimensional p oblems, he D1Q4 model should
be applied along he no mal di ec ion o cell in e ace [18] as is shown in Fig. 2 o he 2D case.
The no mal eloci y
n
U
(Eq. (3)) and angen ial eloci y
(,)
xy n
U uu U
 
  un
in Fig. 2 will
eplace
u
in Eq. (6). Finally, we will ge
nx x
u Un u

 
.
4
1
4
1
4
22
1
4
32
1
3
eq
i
i
eq
ii
i
eq
i ii
i
eq
i iii
i
u
uc
u uc


  
  












(6)
in e ace
U
τ
Un
u
Figu e 2: Applica ion o 1D model o 2D case
By applying Eq. (6) o Eq. (2), we will ge he conse a ion a iable
W
and in iscid lux
c
F
206
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
5
o 2-D low ield as ollowing.
2
2
()
()
( /2 ) /2
nx x
ny y
n
Un u
Un u
U eU

















W
,
2
2
2
2
()
()
( ( /2 ) ) /2
n
n x nx
c
n y ny
nn
U
U pn Uu
U pn Uu
U e pU U















 

F
(7)
2.2 In iscid lux model
Suppose ha cell in e ace is loca ed a
, 1/2 0
cj
x
, hen he in iscid lux a in e ace
c
F
is
decided by no mal eloci y and can be w i en as Eq. (8). In which he momen s
2
1, , / 2 T
a ii p
e


 

φ
and
(0, )
i
is he dis ibu ion unc ion a cell in e ace. Gene ally speaking.
(0, )
i
is summa ion o equilib ium pa
(0, )
eq
i
and non equilib ium pa
(0, )
neq
i
.
4
2
1
2
(0, )
( ( /2 ) )
n
c n i ai
i
nn
U
Up
U e pU






 





Fφ
(8)
To eco e N-S equa ions by Bol zmann equa ion om Chapman-Enskog analysis
[19][21][22][23][24], he non-equilib ium pa
(0, )
neq
i
can be w i en as ollowing.
(0, )
(0, )
neq ii
ii
x







(9)
Then applying Taylo se ies expansion in ome and physics space, Eq. (9) can be simpli ied as
ollowing.
 
(0, ) (0, ) ( , ) ( )
neq
i i ii
   

   
(10)
A he cell in e ace he equillib ium dis ibu ion unc ion is
(0, ) (0, )
eq
ii

.
(,)
ii
 

is
equilib ium dis ibu ion unc ion a su ouding poin o he cell in e ace. Subs i u ing Eq. (9)
in o Eq. (8) as ollowing.
 
0
(0, ) (0, ) (0, ) ( , ) ( )
eq
i i i ii
    
    
(11)
In which
0/
 

is he dimensionless collision ime and

d is he s eaming ime s ep and
ep esen s he physical iscous o N-S equa ions. The con ibu ion o non-equlib ium pa is
always ea ed as nume ical dissipa ion in LBFS. The e o e,
0

can be ega ded as he weigh
o nume ical dissipa ion. We ha e
 
0max ,
LR
 

in his wo k. In which
   
1, 1,
max , max
L R
LR
jj
jN jN
  
 
 
and
, anh
LR
jLR
pp
Cpp







.
L
N
and
R
N
a e he numbe o con ol
olume on le and igh side o cell in e ace. The de ails can be ound in [25].
Subs i u ing Eq. (11) in o Eq. (8) we will ge in iscid lux
c
F
which is decided by
angen ial eloci y as ollowing.
207

Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
6
44 4
0
11 1
0
(0, ) ( , ) (0, )
)
ci ai i ai i i ai
ii i
c cc
     

  


   


  
 
Fφ φ φ
F (F F
Ⅰ Ⅱ Ⅰ
(12)
To al in isicid lux a cell in e ace conside ing angen ial eloci y con ibu ion is shown
as ollowing.
0
)
cc c c


  F F (F F
Ⅰ Ⅱ Ⅰ
(13)
c

F
Ⅰ
ep esen s con ibu ion o equilib ium dis ibu ion unc ion
(0, )
eq
i
a cell in e ace and
c
F
Ⅱ
means equilib ium dis ibu ion unc ion
(,)
ii
 

a su ouding poin o he cell in e ace.
g
4L
in e ace
R
L
g
4R
g
2Rg1L
g
3L
in e ace
R
L
g
2L
g
1L
g
3L
g2Rg1Rg3R
g4R
s eaming
Figu e 3: S eaming p ocess based on D1Q4 model a cell in e ace
Supposing ha a local Riemann p oblem is o med a cell in e ace, which is shown in Fig.
3. Hence, equilib ium dis ibu ion unc ion
(,)
ii
 

a su ouding poin o cell in e ace
can be decided by he loca ion o
i


, which is shown as Eq. (14).
0
(,) 0
L
ii
ii R
ii

  


 



(14)
Mo e speci ically, Eq. (14) can be w i en as ollowing o D1Q4 model.
1, 3
(,) 2,4
L
i
ii R
i
i
i
 


 



(15)
By applying he ela ionship o
(0, ) (0, ) (0, )
eq neq
ii i
 
, he conse a ion lux which is
decided by angen ial eloci y can be w i en as ollowing.
44
11
2
(0, )+ (0, )
( /2 )
Meq neq
n ai ai
ii
n
U
Ue



 


 





Wφ φ
(16)
208
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
7
The subsc ip M means pa ame e on cell in e ace. The non-equilib ium pa has no
con ibu ion o calcula e conse a ion a iables acco ding o he compa ibili y condi ion [19],
which is shown as ollowing.
 
44
0
11
(0, ) (0, ) ( , ) 0
neq
ai a i i i
ii
  
 
   

φ φ
(17)
Subs i u ing Eq.( 17) in o Eq. (16) we will ge
44
11
(0, ) ( , )
M
ai ai i
ii
 
 
  

Wφ φ
(18)
By applying Eq. (15) in o Eq. (18), he conse a ion lux which is decided by angen ial eloci y
M
W
can be w i en as
1,3 2,4
MLR
ai ai
ii
 
 

Wφ φ
(19)
The e o e, densi y, no mal eloci y and p essu e a cell in e ace a e ob ained om Eq. (19).
Then he angen ial eloci y a cell in e ace can be calcula ed by he ollowing equa ion.
1,3 2,4
M LL RR
ii
ii
U U U
 

 
 

()
(20)
In which
M
U

,
L
U

and
R
U

a e angen ial eloci y a cell in e ace, on he le and igh side o
cell in e ace espec i ely. Based on Eq. (19) and Eq. (20), all a iables a cell in e ace such
as
M

,
M
p
M
U

and
M
n
U
can be ob ained. Subs i u ing hese a iables o Eq. (5) and we can ge
he equlib ium dis ibu ion unc ion
(0, )
i
a cell in e ace.
The
c

F
Ⅰ
in Eq. (13) including wo pa s, which a e
c

F
Ⅰ
decided by no mal eloci y and he pa
decided by angen ial eloci y. Then
c

F
Ⅰ
can be w i en as ollowing by applying a iables a
cell in e ace.
2
2
2
2
()
()
( ( /2 ) ) /2
n
n x nx
c
n y ny
n nn
U
U pn Uu
U pn Uu
U e pU U U















 

F
Ⅰ
(21)
Also
c
F
Ⅱ
is consis s o
c
F
Ⅱ
decided by no mal eloci y and he pa decided by angen ial
eloci y. Subs i u ing equilib ium dis ibu ion unc ion
(,)
ii
 

a su ouding poin o cell
in e ace o Eq. (8) we will ge
c
F
Ⅱ
as ollowing.
1,3 2,4
LR
c ii ii
ii


 
 

F
Ⅱ
 
(22)
Then he angen ial eloci y a cell in e ace can be calcula ed by he ollowing equa ion.
209
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
8
1,3 2,4
2 22
1,3 2,4
M LL RR
n ii ii
ii
M LL RR
n ii ii
ii
UU U U
UU U U
 
 
 
 
 
 
 
 


()
()
(23)
Now
c
FⅡ
can be exp essed as ollowing. In which
(,)
i ii
 

.
4
1
44
11
44
11
44
2
11
11
()
22
ii
i
M
iii x ii x
ii
c
M
iii y ii y
ii
M
i ii p i ii
ii
n U
n U
e U




 
 
  

 
 
 



















F
Ⅱ
(24)
By subs i u ing Eq. (21) and Eq. (24) o Eq. (13) we ge in iscid lux
c
F
on cell in e ace.
2.3 Viscous lux model
In his wo k, cen al di e ence scheme is applied o sol e he isous lux
F
. S ess enso

and lux ec o

in Eq. (2) a e gi en by
 
 
 
2
=+ 3
T
kT



   


   
u u uI
u
(25)
In which

is dynamic iscosi y and decided by Su he land’s law and u bulence model.
k
is
he mal conduc i i y and
I
is uni ma ix.
Fo physical a iables on cell in e ace, such as
andu, , w k
, can be calcula ed by a i hme ic
mean alue on le and igh side o cell in e ace, shown as ollowing.
1()
2
M LR
 
 
(26)
Whe e

ep esen s any physical a iables. The de i a i es in Eq. (25) a e calcula ed by ini e
di e ence scheme and de ils a e shwon in e e ence [26].
2.4 Me hod o imp o e nume ical s abili y
Mul iblock g ids a e applied o e ine local g id quali y. G ids local e inemen me hod is
used nea wall su ace and block in e ace. Hence, scheme s abli y has been imp o ed by
educing spa ial s ep, especially in he leewa d side o low ield. Besides, a cons ain
judgemen is added in ou code o a oid nega i e alue in calcula ing p ocess. In his wo k,
implici LU-SGS me hod [27] is applied o speed up he con e gence a e and keep calcula ion
obus .
210
Z. X. Meng, C. Shu, L. M. Yang, W. H. Zhang and S. Z. Li
9
3 NUMERICAL SIMULATION
Biconics model is used o e i y he p esen sol e . The shape o biconics model is shown
as Fig. 4. The head cu a u e adius is 3.84mm, semi-cone angle is 12.84 deg ee in on and 7
deg ee in back, while he leng h o on cone is 69.55mm and 122.24mm o ally. The ee
s eam has a p essu e o
59.92PaP


， empe a u e o
48.88KT


and Mach numbe
9.86Ma


.
Fig. 5 shows mach numbe dis ibu ion ob ained by p esen sol e in his a icle. I can be
seen clea ly ha LBFS can cap u e s ong shock wa es and low inside he bounda y laye
exac ly. The maximum Mach numbe is 16.93Ma.
Fig. 6 shows he p essu e con ou s bibonics model wi hou ail low ield, and Fig. 7 is he
one has ail ield. Fig. 8 and Fig. 9 show empe a u e dis ibu ions. One block g ids a e applied
o low ield wi hou ail simula ion. I is easie o calcula ion con e gen wi hou low p essu e
e ec in back ail low.
Figu e 4: Measu emen s o biconics model
Figu e 5: Mach numbe con ou s ob ained by LBFS
a Ma=9.86
Figu e 6: P essu e con ou s ob ained by LBFS
wi hou back ail low ield
Figu e 7: P essu e con ou s ob ained by LBFS wi h
back ail low ield
211