设为首页 加入收藏 期刊导航 网站地图
  • 首页
  • 期刊
    • 数学与物理
    • 地球与环境
    • 信息通讯
    • 经济与管理
    • 生命科学
    • 工程技术
    • 医药卫生
    • 人文社科
    • 化学与材料
  • 会议
  • 合作
  • 新闻
  • 我们
  • 招聘
  • 千人智库
  • 我要投搞
  • 办刊

期刊菜单

  • ●领域
  • ●编委
  • ●投稿须知
  • ●最新文章
  • ●检索
  • ●投稿

文章导航

  • ●Abstract
  • ●Full-Text PDF
  • ●Full-Text HTML
  • ●Full-Text ePUB
  • ●Linked References
  • ●How to Cite this Article
Advances in Condensed Matter Physics 凝聚态物理学进展, 2013, 2, 88-96
http://dx.doi.org/10.12677/cmp.2013.24012 Published Online November 2013 (http://www.hanspub.org/journal/cmp.html)
Polar Coordinate Lattice Boltzmann Studying of
Shocking Process*
—Investigation of Non-Equilibrium Effects in Complex System
Chuandong Lin1, Aiguo Xu2#, Guangcai Zhang2, Yingjun Li1#
1State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining and Technology, Beijing
2National Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing
Email: #Xu_Aiguo@iapcm.ac.cn, #lyj@aphy.iphy.ac.cn
Received: Sep. 17th, 2013; revised: Oct. 8th, 2013; accepted: Oct. 17th, 2013
Copyright © 2013 Chuandong Lin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Abstract: The Polar Coordinate Lattice Boltzmann (PCLB) model is used to investigate the shocking processes. The
discrete velocity model used in this PCLB is the D2V97 presented by Watari and Tsutahara. The convective term is
treated with a modified Warming-Beam scheme. Two shocking processes are investigated. In one case, the evolution of
a shock wave travels inwards within an annular area is studied. In the other case, the evolution of a shock wave travels
inwards within a circular area is studied. The non-equilibrium effects around the shock fronts are investigated by ana-
lyzing the simulated results of velocity moments of local discrete distribution function. To further understand the inter-
face from a more fundamental level, we give the sketch of the actual distribution function in velocity space around the
shock front.
Keywords: Lattice Boltzmann Method; Shock Wave; Non-Equilibrium Effects
使用极坐标格子玻尔兹曼方法研究冲击过程*
—复杂系统中非平衡效应的探索
林传栋 1，许爱国 2#，张广财 2，李英骏 1#
1中国矿业大学(北京)深部岩土力学与地下工程国家重点实验室，北京
2北京应用物理与计算数学研究所计算物理重点实验室，北京
Email: #Xu_Aiguo@iapcm.ac.cn, #lyj@aphy.iphy.ac.cn
收稿日期：2013 年9月17 日；修回日期：2013年10月8日；录用日期：2013 年10 月17 日
摘 要：本文使用极坐标高速可压格子玻尔兹曼模型研究冲击过程。其中离散速度模型是由 Watari 和Tsutahara
提出的 D2V97，对流项采用修正的 Warming-Beam 格式处理。分别模拟研究了冲击波在环形区域和圆形区域中
向内传播过程中物理量的变化规律，并通过离散速度分布函数的速度矩分析了冲击波波阵面附近的非平衡效应。
为了从更基本层面理解界面行为，给出了波阵面处真实分布函数在速度空间的示意图。
关键词：格子玻尔兹曼方法；冲击波；非平衡效应
*资助信息：本研究得到了中国国家自然科学基金项目[批准号：11074300 和91130020]、中国国家重点基础研究发展计划[批准号：2013
CBA01504]的资助
#通讯作者。
Open Access
88
使用极坐标格子玻尔兹曼方法研究冲击过程
Open Access 89
1. 引言
最近几十年，格子玻尔兹曼(Lattice Boltzmann,
LB)方法作为一种联系宏观和微观行为的介观方法，
在众多复杂流体领域的模拟和应用方面取得了巨大
成功[1,2]。多数早期的LB 模型只能用于等温不可压缩
流动的模拟。而对于普遍存在的高速可压流体，比如
在航空航天(绕流湍流)、天体物理(射流)、冲击爆轰(瓦
斯爆炸)等领域，流体的压缩效应不可忽略。对于此类
问题的研究就需要借助于能够模拟高速可压流体的
LB 模型[3-20]。当前的高速可压LB 模型大都是基于笛
卡尔直角坐标系。实际上，基于极坐标、柱坐标或球
坐标的 LB模型将会更加方便处理一些物理系统，比
如惯性约束核聚变(Inertial Confinement Fusion, ICF)、
超新星爆炸、内爆外爆等。
Figure 1. Sketch of discrete velocity model
图1. 离散速度模型示意图



cos 21
sin 21
kir ki
ki ki
iN
iN






 





 




(2)
这里选取 24
i
N

，00


，11.5

，23.5


，
37.5


，412.5


。
本文使用极坐标LB 模型模拟研究在二维系统中
冲击波向内部传播的物理过程，并根据离散速度分布
函数的速度矩(包括矩和中心距)的数值结果重点分析
了冲击波波阵面附近系统的非平衡效应和表现，进而
恢复了真实分布函数在速度空间的部分主要特征。需
要说明的是，LB 模型秉承了玻尔兹曼方程可以描述
系统非平衡行为的功能和属性。因此，LB 模型可用
于模拟研究典型物理过程中的非平衡现象。本文所获
得的关于非平衡效应的结果对于流体动力学系统的
宏观建模(比如基于 Euler方程、Navior-Stokes 方程)
的改进具有重要的启发和指导作用。
局部平衡分布函数的表达式为：
24 2
2
2
23
4
11
22
8
12
26
24
eq ki
ki k
ki kiki kiki
ki kiki ki
u
uu u
fF EEE
E
uu uuu
u
E
EE
uuuu
E


  



 



 















+



(3)
式中




1
22 2222 2
123
4222322
43123 21
2222 2222
23 311123
00 1234
1
k kkkkkkk
kk kkk
kk kkkkk
F
BEBE B
EB E
FBFFFF
   
2

 
 

 

 
 

 







++
++ +
2. 模型简介
基于极坐标的 Bhatanger-Gross-Krook 近似下的有
限差分 LB 热学模型的演化方程为：
11
eq
ki kiki
kirkiki ki
ff f
f
f
trr



 
 

 

(1)
其中 024B

，1112B

，213B，32B

，
4
B16

。
其中 ki
f
与eq
ki
f
分别为局部离散分布函数和平衡分布
函数，

为决定恢复平衡态速度的松弛因子，


,r

和

,ki

kir


分别为物理空间和速度空间的坐标。
运用 Chapman-Enskog展开法，本模型在流体力
学条件下能够恢复到 Navier-Stokes 方程：

0
t



 
u (4)
2.1. 离散速度模型
本文采用由 Watari和Tsutahara[21-23]提出的离散
速度模型，该模型有 5组离散速度，速度大小为为
, 0,1,2,3,4
kk

，并且每组都个分速度，见图1。
i
N
 
 
T0
E
t



 



 

uIuu
uI uu
(5)
离散速度的数学表达式为：
使用极坐标格子玻尔兹曼方法研究冲击过程


22
2
11
2
22
10
2
EuE u
t
Eu
 
 
 

 
 





uu
uu uu
(6)
其中

为温度，

rr
uu


uee为速度，T为温度，

PT

为压强，


1ET



为单位质量的内能，
为绝热系数，

2




与 分别表示粘性系数和热传
导系数，


2P2



 。
2.2. 算法格式
对公式(1)中的时间项采用解析解，空间项采用修
正的 Warming-Beam格式处理，并引入符号

,, ,, ,,,
it
kiki irki iki

f
fffiriitfiriki


 ，可得 极
坐标中 LB 方程的演化形式[18]：
1
121
itt r r
ki 2
f
termterm term termterm


 (7)
其中
d1
dd
rkir ki
tt
CC
rr


d

，
dd
exp exp
tteqeq
ki kiki
tt
termf ff



 
 
 

 



,,1
1
,1 ,
if 0
if 0
r kiirkiirr
r
r kiirkiirr
Cf fC
term
Cf fC













2
2
2
2, ,1,2
2,1,,2
11 if
11 if
22
22
rr rrr
r
rr rrr
r kiir kiirkiir
rkiirkiir kiir
CCS SC
term
CCS SC
Sf ff
Sfff








 

 



 

0
0


,,1
1
,1 ,
if 0
if 0
ki ikii
ki iki i
Cf fC
term
Cf fC
 




















2
,,1,2
,1 ,,2
11 if
11 if
22
22
i
i
i kii kiikii
ikiikii kii
CCS SC
term
CCS SC
Sf ff
Sf ff
 

 















 

0
0


,1,2
,,1
,2 ,1
,1,
if 0
if 0
ki irki ir
r
ki irkiir
r
ki irkiir
r
ki irki ir
ff C
ff
ff C
ff













,1 ,2
,,1
,2 ,1
,1 ,
if 0
if 0
ki iki i
ki ikii
ki iki i
ki ikii
ff C
ff
ff C
ff




















需要说明的有四点：第一，it 是演化时刻，dt 是
时间步长， ir 与i

是极坐标系中的网格坐标， 与dr
d

是空间步长；第二，当r

与


的分母趋向于零时，
其r

与


的绝对值趋向于无穷大；第三，公式(7)可
以看做是一阶精度的算子分裂格式，并且没有考虑数
值计算的演化顺序；第四，

r
S

与


S


是开关函
数，如果令




r

1

SS


，那么 ，
这样公式(7)就变成了一阶迎风格式；第五，需要强调
的是开关函数要满足
2
r
rmte 2
0te rm





112S


，符合此条
件的函数有很多，比如：

12
12
1
,1,1
1
SS
SS
S


 3
S

 
 (8)
或




0
0
0
Exp1 1
Exp 1
S
SS
S





， (9)
不同的开关函数对计算精度有影响，选择什么样
开关函数最接近精确解，这项工作正在探索中。本文
的模拟过程采用的是一阶迎风格式。
2.3. 边界条件
本文模拟研究的物理系统所在的空间是扇形区
域，角度范围为 02
s
i
NN


 ，其中 1, 2,,
s
i
NN。
当
s
i
NN

时，扇形就变成了圆。本文选取 1
s
N

。
2.3.1. 径向边界
径向网格节点编号


1, 2,,1,
r
irN Nr
，在物
理区域的内侧和外侧分别增加两层虚拟网格


2
r
1, 0,1,
r
irN N

，这样网格节点便拓展为


1,0, 1,,,1,
rr
NN 2
r
N

。
环形区域：内、外半径分别为 与，径向步
长
1
R2
R




21 1
r
drR RN

，节点到原点的距离为


11rR irdr

。假定边界虚拟点上的物理系统处
于热力学平衡状态，那么分布函数就可以通过平衡分
布函数公式计算得到： eq
ki ki
f
f；对于非平衡微观边
界条件，偏离热力学平衡的部分 可以通过

ki 

eq
ki
ff
Open Access
90
使用极坐标格子玻尔兹曼方法研究冲击过程
插值格式计算得出。虚拟网格点上的宏观物理量可以
对边界内的物理场使用插值格式处理。
圆形区域：半径为 ，径向步长Rr
drR N，节
点到原点的距离为 12rirdr

。其径向外边界的处
理方式与环形区域的处理方法一样，而内边界虚拟网
格点 上的分布函数为[18]

1, 0ir 
 
 
,,,1,,, 12if1224
,,,1,,,12 if1224
firikifirikii
firikifirikii


 


 


2.3.2. 切向边界
径向网格节点编号


1, 2,,1,iNN





i

1,0, 1,,
，切向
角度步长dθ = 2π/(24Nθ)，每个节点所在的角度为θ =
iθ × dθ。在切向两侧分别增加两层虚拟网点
，得到新的网格编号

1, 0,
,N


1, 2NN



1, 2NN

。在虚拟点的分布函数计算方法如下：



,,,1if1
,,, ,,,23if
firNi kii
firi ki
f
ir iNk iiN





 

 



其中

1if 23
11 if24
ii
ii





，

24 if1
23 1if1
i
iii





3. 非平衡效应的宏观表现
在公式(3)中，局部平衡分布函数满足以下七个矩
关系：

0
ki
ki
f


 (10)

0
ki ki
ki
f


vu
(11)


2
02
ki ki
ki
f
E


vu (12)


0
kiki ki
ki
fE


vvIuu



(13)


0
kiki ki ki
ki
f
E
 




vvv
ue ee ueeeuuuu (14)


022
kikiki ki
ki
fE


vvvu uu

2 (15)


0222
kikiki ki ki
ki
fE


vvvvuuI uu

E (16)
在公式(10)~(12)中，eq
ki
f
可以用 ki
f
替换，而在公
式(13)~(16)中如果用 ki
f
替换 eq
ki
f
，那么这四个等式将
不再成立，此时这四个公式两侧的差异就是系统偏离
热力学平衡造成的。这里给出两种物理量：矩 与
中心距 ，
m
M
*
m
M





2
3
3,1
4,2
2
2
kiki ki ki
ki
kiki ki ki ki
ki
kiki kiki ki
ki
kiki ki kiki ki
ki
ff
ff
ff
ff











Mvv
M vvv
Mvvv
Mvvvv
(17)





***
2
* ***
3
***
3,1
***
4,2
2
2
kiki ki ki
ki
kiki ki ki ki
ki
kiki kiki ki
ki
kiki ki kiki ki
ki
ff
ff
ff
ff













Mvv
M vvv
Mvvv
Mvv
*
**
vv
(18)
其中 *
ki ki

vvu，下标“3,1”表示三阶张量缩并为一
阶张量，“4,2”表示四阶张量缩并为二阶张量。矩


1 3,1,
M
3,


M
1,
e
是矢量，包含两个组分 和
3,1,r
M
3,
M

；矩


22,
M


eM
2,rr 2,r
M
e是二阶张量，包含四个组
分 、M

、2, r
M

、2,
M


，其中 2,
M2,rr
M



，
也是就是说 有三个独立变量、
2
M2,rr
M2,r
M

、2,
M


；
同理，


4,2,
M
4,2


M
2,rr 4, 2 ,r
M
ee也只有三个独立变量
、
4,
M

、4,2,
M


；

33
M,

 
eeeM是三
阶张量，含有八个组分，其中只有四个独立变量
、
3,
Mrrr 3,rr
M

、3,r
M


M、3,


。对于中心距 在
数学上与矩 有类似的分析结果。
*
m
M
m
M
在概率论里，对于一维分布函数

f
v，中心距
 
3
*
3d
M
fv vuv
描述的是分布函数


f
v的“ 偏
斜度”，中心距
 
4d
*
4
M
fv v uv
描述的是分布
函数


f
v

的“扁平度”。对于高斯分布函数

2
12exp 2fvv u




，可以算出 *
4
M3

。
如果一个分布函数对应的且 ，那么这个
分布函数在中心位置将比高斯分布函数尖锐；如果一
个分布函数对应的
*
4
M
*
43M
3*
21M

且 ，那么这个分布函
数在中心位置将比高斯分布函数平缓。
*
2
M1
在物理意义上， 、 描述的是粒子涨落不同
自由度的2次和3次关联； 是粒子涨落在不同自由
度3次关联的一次缩并； 是粒子涨落在不同自由
度4次关联的一次缩并。矩 和中心距 的迹都是
守恒量，其迹的数值大小与温度密切相关；矩 和
2
M3
M
3,1
M
4,2
M
2
M*
2
M
2
M
Open Access 91
使用极坐标格子玻尔兹曼方法研究冲击过程
中心距 的斜对角上的组分与流体的剪切作用有
关。矩 和描述了两部分物理量的叠加：系统
宏观流动引起的能量通量和系统非平衡过程的“微观
涨落能流”1。中心距 和描述了系统非平衡过
程(比如热扩散、粘性等过程)的“微观涨落能流”。
处于非平衡态的系统存在“微观涨落能流”，矩 、
和中心距 、都不为0。
*
2
M
3
M

fv

fv
u
3,1
M
*
3
M

v

v
mm

**
mm


*
3
M
*
3,1
M

ki
M

ki
M
*
3,1
M

m
M

*
m
M
M



3
M
*
m
M
3,1
M

f

f
*
m

另外，平衡态系统的粒子 速度 分布函 数满足
，而非平衡态系统的粒子速度分布函数
。为了描述非平衡态系统粒子速度分布
函数偏离对称性的程度，可以通过矩和中心距进一步
定义关系式：

f

f
eq
ki
f

(19)
eq
ki
f (20)
通过伽利略不变性可知，矩和 同时包含宏
观流动 和微观热力学运动的信息，而中心距和
则单纯描述了微观分子的热力学涨落过程。
m m

4. 数值模拟与研究
4.1. 雨贡纽关系
当冲击波稳定传播时波前和波后的物理量满足
雨贡纽(Hugoniot)关系：


00
00


0 0
000
11 2
H
Du
uuu
P


 

H
H




HH
D
P


Du
PP
EE


(21)
其中下标 0和H分别表示冲击波波前和波后的物理
量，D为冲击波的传播速度。如果波前和波后的状态
方程满足理想气体状态方程，

1

EP ，那 么
公式(21)中的能量方程可写成


0
P
0
00
11
2H
PP

1
1
H
P



 








 (22)
4.2. 数值模拟
4.2.1. 环形区域
在内外半径分别为 ，的环形
物理区域内，冲击波由外向内传播，波速D = 2，初
始时刻的物理场为：
1
R200022002R






, ,,1,1,0,1
,,,1.50000,0.66667, 0, 2.33333
ri
ro
uu P
uu P









这里下标 i表示 ，该部分为冲击
波波前区域，下标 表示 ，该部分
为波后区域 2，波阵面前后的物理量满足雨贡纽关系
(21)。计算网格为
2000 2001.6r
o2001.6 r
2000 3
r
NN

2002

。
图2为冲击波在环形区域向内传播过程中物理量
沿径向的一维分布图，模拟得到的波后物理量的大小
为




,,,1.50031,1.55593, 2.33437,0.66686
r
TPu

，
即相对于初始时刻的波后物理量的数值分别增加了


0000
00 00 00 00
0.2 ,0.2 ,0.4 ,0.3。考虑到环形区域的几何
效应，冲击波波后区域的物质随着时间演化将不断聚
集压缩，数值模拟得到的波后参数应该比初始条件高
一些。
4.2.2. 圆形区域
在半径为 1R

的圆形区域，冲击波以波速 2D

向内传播，初始时刻的物理场为：






, ,,1,1,0,1
,,,1.50000,0.66667, 0, 2.33333
ri
ro
uu P
uu P









这里下标 i表示区域 00r.8

，下标 表示区域
o
0.8 1r

。计算网格为 。图 3为冲
击波在圆形区域传播过程物理量沿径向分布的演化
1000 3
r
NN


Figure 2. Physical quantities (density, temperature, pressure, ve-
locity) versus radius in the progress of a shock travelling inwards
within an annular area
图2. 冲击波在环形区域向内传播时物理量沿径向的分布图
2初始条件给出的是强间断的物理场，而实际的冲击波波阵面是一
个光滑过渡层，因此在数值模拟的初始阶段物理量分布不光滑。在
数值模拟一段时间后，冲击波将稳定传播，物理量也会光滑分布。
然后，我们可以截取冲击波稳定传播时的物理场作为初始条件再次
进行数值模拟。
1目前所使用的物理学概念大都是描述平衡态系统的物理参数，而
这里用于描述非平衡现象的高阶矩在流体力学中没有对应概念。本
文使用的新概念与流体语言中的传统概念在微观机理上可能不同。
Open Access
92
使用极坐标格子玻尔兹曼方法研究冲击过程
Open Access 93
(a) (b)
(c) (d)
Figure 3. Physical quantities versus radius in the evolution of a shock travelling within a circular area: (a) Density, (b) Temperature,
(c) Velocity, (d) Pressure
图3. 冲击波在圆形区域内传播时的物理量沿径向分布的演化图：(a) 密度，(b) 温度，(c) 速度，(d) 压强
图，选择的时刻分别是t = 0.000，0.200，0.345，0.500。
在初始阶段，冲击波向内传播，波后区域的物质
具有方向向内的宏观速度，由于惯性作用波后不断有
物质涌入，该区域的物质不断被挤压，密度、温度和
压强不断上升；同时，被挤压的物质将相互排斥，在
斥力的作用下波后区域的这部分物质的速度将不断
下降。需要说明的是，在冲击波向内过程，波阵面附
近的温度不断上升，冲击波的波速不断加快，见图
3(c)；
当t = 0.345 时，冲击波到达原点，此时的密度、
温度、压强的大小在径向单调分布，极点附近的数值
最大，向外依次递减，见图3(a)、(b)和(d)；
冲击波向外传播时，波后区域的物质具有向外的
宏观速度，在原点附近的速度最小，在波阵面处的速
度最大。波前区域物质的宏观速度指向内，并且该部
分物质的速度不断减小，而密度、温度、压强不断上
升。
4.3. 非平衡效应的分析
4.3.1. 非平衡效应的宏观表现
例如图 2，在冲击波稳定传播过程，波前的物理
系统处于平衡态，波后的物理系统也将很快达到新的
平衡态，而波阵面附近的非平衡效应则表现的比较明
显。图 4给出了对应于图 2所示冲击波波振面附近的
中心距 在径向的分布图。这里只给
出了各个中心距的独立变量。从图 4可以看出以下几
点：
*** *
23314
,, ,MMMM
，，
2
第一，图 4(b)、(d)、(f)、(g)中显示
的数值都比较小，数量级是10−3，这说明我们所研究
的非平衡系统偏离平衡态的程度并不是很大；
*** *
233142
,, ,
，，
ΔΔΔ Δ
第二，图 4(c)中的 与图(d)中的 的图像一样，
图(e)中的 与图(f)的图像一样。这是因为在理
论上
*
3
M
*
31
，
Δ
*
3
Δ
*
31
M，


eq
f
*
30

M，


0
eq
f
*
3,1
M，故


**
33
f

MΔ，


f
**
31

，，
Δ
31
M；
第三，图 4(b)中*
2,rr

与*
2,


的大小相等，正负号
相反，即 
**
2, 2
rr ,0


。这是因为


*
2,
eq
rr
f
M与


*
2,
eq
f

M的物理量纲分别对应 与r

自由度上的内
能，而 与
*
2,
rr
*
2,


分别表示非平衡系统在与r

自由
度上偏离能量均分状态的程度3。
第四，在图 4(d)、(f)、(h )中 的峰
值最大，而
***
3,3,1,4,2,
,,
rrr rrr

***
3,3,1,4,2,
,,

 
的数值接近于零。这是
因为它们的物理量纲都与“能量涨落能流”有关。由
3能量均分定理描述的是平衡态系统的物理规律，处于非平衡态的
物理系统在各个自由度的能量不一定相等。
使用极坐标格子玻尔兹曼方法研究冲击过程
Figure 4. Central moments and their corresponding non-equilibrium manifestations in the progress of a shock wave travelling inwards
within an annular area. Figures (a)-(h) are for , respectively
MMM M
** ** ****
22 3331314242
ΔΔ ΔΔ,
，， ，，
，， ，，， ，
**
22
，，
图4. 冲击波在环形区域向内传播时 M*与在径向的数值分布。图(a)~(h)依次对应
*
ΔMMM M
***** *
33 31314242
ΔΔ ΔΔ,
，， ，，
，，， ，
于冲击波的冲击作用，波阵面附近的物质粒子在速度
空间分布的对称性将得到破坏，系统在 自由度上将
存在较大的“能量涨落能流”，而在
r

自由度上基本
上没有“能量涨落能流”。
4.3.2. 非平衡效应的微观机制
在冲击波波阵面附近，系统非平衡效应的演化过
程可以解释如下：冲击波是一个能量传播的过程，在
冲击波附近会产生较大的宏观物理量的梯度。在冲击
波作用下，物质粒子在 自由度上的能量(包括动能、
r
内能)首先得到增加，然后这部分变化的内能将通过分
子间的碰撞传递到

自由度。同时， r自由度的内能
在冲击波作用下继续升高，如此往复，直到宏观物理
量的梯度完全消失，系统在各个自由度的能量均分，
系统将达到新的热力学平衡。
根据 *
m

的模拟数值，我们可以分析得 知真实分
布函数在速度空间的一些主要特点，分析步骤如下：
第一步：图4(b)显示，在波阵面附*
2, 0
rr

，
*
2, 0


。这说明


r
f
v比麦克斯韦分布函数“扁平”，


f
v

比麦克斯韦分布函数“瘦高”，并且


r
f
v

的
峰值低于麦克斯韦分布函数的峰值，而
f
v

的峰值
Open Access
94
使用极坐标格子玻尔兹曼方法研究冲击过程
高于麦克斯韦分布函数的峰值；根据图4(d)中*
3

和(f)
中的分布，我们可以得知
*
3,1


f
v

是对称分布的，
而


r
f
v是非对称分布的，并且

r

f
v在的部分
比在 的部分“胖”。我们称之为“负偏斜”。
图5给出了在冲击波波振面处分布函数在一维速度空
间
0
r
v
0
r
v
r


与


示意图，其中长划线代表

r
f
v，短划 线代
表

f
v

，实线线代表麦克斯韦分布函数 eq
f
。需要
指出的是，麦克斯韦分布函数描述的是平衡态系统的
粒子速度分布律，此时系统的粒子速度分布具有各向
同性，也就是系统满足能量均分定理，所以图5中的
麦克斯韦分布函数

eq
r
f
v和


eq
f
v

是重合的，这里
统一用 eq
f
表示。
第二步：由图4(b)中 可知冲击波波振面附
近的分布函数在速度空间

*
2, 0
r


,
r




的轮廓投影应该关
于r

或


成轴对称分布。根据图5给出的分析结果可
知分布函数

f
v

在速度空间


是对称分布的，因此
分布函数在速度空间的轮廓投影的对称轴就是坐标
轴r

。图 4(h)可以给出一致的结论。图6给出了冲击波
波振面处的真实分布函数在二维速度空间

,
r



的
轮廓投影示意图。
第三步：联立图5和图6的分析结果就可以得到真
实速度分布函数在速度空间


,
r



的三维示意图，见
图7。
5. 总结
本文使用极坐标高速可压格子玻尔兹曼模型研
究了內爆冲击过程。所用离散速度模型是由Watari和
Tsutahara提出的D2V97，空间项使用修正的Warming-
Figure 5. Sketch of the Maxwellian and actual distribution func-
tions versus velocity r

and


, respectively. The long dashed
line is for distribution function


r
f

, the shot dashed one is for
distribution function

f



, and the solid line is for e
q
f.
图5. 真实分布函数和麦克斯韦分布函数在速度空间 r

与


的一
维示意图。

r
f

、

f


与e
q
f分别用长划线、
短划线和实线表示。
v
θ
v
r
Figure 6. The sketch of contours of the actual distribution func-
tions in velocity space

,
r




,

图6. 真实分布函数在速度空间

r


的投影轮廓示意图
Figure 7. The 3-dimensional sketch of the actual distribution func-
tions in velocity space

,
r





图7. 真实分布函数在速度空间

,
r


的三维分布示意图
Beam 格式处理，时间项使用解析算法。模拟研究了
冲击波在环形区域和圆形区域中向内传播时宏观物
理量(密度、温度、压强和速度)在径向的分布情况。
通过分析离散速度分布函数的高阶矩的空间分布特
征，研究了冲击波波阵面附近系统的非平衡效应，并
给出了真实分布函数在速度空间的主要特征。这些研
究结果有助于从更基本的层面上理解界面附近的非
平衡行为，也为传统数值模拟中物理建模的修正提供
了重要参考。
参考文献 (References)
[1] Succi, S. (2001) The Lattice Boltzmann Equation for Fluid Dy-
namics and Beyond. Oxford University Press, New York.
[2] 许爱国, 张广财, 李华, 朱建士等. (2011) 材料动力学的介观
模拟(北京应用物理与计算数学研究所讲义)[Z]. 北京.
[3] Xu, A., Zhang, G., Gan, Y., Chen, F. and Yu, X. (2012) Lattice
Boltzmann modeling and simulation of compressible flows.
Frontiers of Physics, 7, 582-600.
[4] Pan, X., Xu, A., Zhang, G. and Jiang, S. (2007) Lattice Boltz-
mann approach to high-speed compressible flows. International
Journal of Modern Physics C, 18, 1747-1764.
[5] Gan, Y., Xu, A., Zhang, G. and Li, Y. (2008) Finite-difference
Lattice Boltzmann scheme for high-speed compressible flow:
Two-dimensional case. Communications in Theoretical Physics,
50, 201-210.
[6] Gan, Y., Xu, A., Zhang, G., Yu, X. and Li, Y. (2008) Two-dimen-
sional Lattice Boltzmann model for compressible flows with
Open Access 95
使用极坐标格子玻尔兹曼方法研究冲击过程
Open Access
96
high mach number. Physica A, 387, 1721-1732.
[7] Gan, Y., Xu, A., Zhang, G. and Li, Y. (2011) Flux limiter Lattice
Boltzmann scheme approach to compressible flows with. Com-
munications in Theoretical Physics, 56, 490-498.
[8] Gan, Y., Xu, A., Zhang, G. and Li, Y. (2011) Lattice Boltzmann
study on Kelvin-Helmholtz instability: Roles of velocity and
density gradients. Physical Review E, 83, Article ID: 056704.
[9] Chen, F., Xu, A., Zhang, G., Gan, Y., Tao, C. and Li, Y. (2009)
Lattice Boltzmann model for compressible fluids: Two-dimen-
sional case. Communications in Theoretical Physics, 52, 681-
694.
[10] Chen, F., Xu, A., Zhang, G., Li, Y. and Succi, S. (2010) Multi-
ple-relaxation-time Lattice Boltzmann approach to compressible
flows with flexible specific-heat ratio and prandtl number. Eu-
rophysics Letters, 90, Article ID: 54003.
[11] Chen, F., Xu, A., Zhang, G. and Li, Y. (2010) Three-dimensional
Lattice Boltzmann model for high-speed compressible flows.
Communications in Theoretical Physics, 54, 1121-1128.
[12] Chen, F., Xu, A., Zhang, G. and Li, Y. (2011) Multiple-relaxa-
tion-time lattice Boltzmann model for compressible fluids. Phys-
ics Letters A, 375, 2129-2139.
[13] Chen, F., Xu, A., Zhang, G. and Li, Y. (2011) Flux limiter Lattice
Boltzmann for compressible flows. Communications in Theo-
retical Physics, 56, 333-338.
[14] Chen, F., Xu, A., Zhang, G. and Li, Y. (2011) Prandtl number
effects in MRT lattice Boltzmann models for shocked and un-
shocked compressible fluids. Theoretical and Applied Mechan-
ics Letters, 1, Article ID: 052004.
[15] Chen, F., Xu, A., Zhang, G. and Li, Y. (2011) Multiple-relaxa-
tion-time Lattice Boltzmann approach to Richtmyer-Meshkov.
Communications in Theoretical Physics, 55, 325-334.
[16] Gan, Y., Xu, A., Zhang, G. and Yang, Y. (2013) Lattice BGK
kinetic model for high speed compressible flows: Hydrodynamic
and nonequilibrium behaviors. EPL, 103, Article ID: 24003.
[17] Chen, F., Xu, A., Zhang, G. and Wang, Y. (2013) Two dimen-
sional MRT LB model for compressible and incompressible
flows. Frontiers of Physics, in press. arXiv: 1305.4736.
[18] Lin, C., Xu, A., Zhang, G., Li, Y. and Succi, S. (2013) Polar
coordinate lattice Boltzmann modeling of compressible flows.
e-print arXiv: 1302.7104.
[19] Yan, B., Xu, A., Zhang, G., Ying, Y. and Li, H. (2013) Lattice
Boltzmann model for combustion and detonation. Frontiers of
Physics, 8, 94-110.
[20] Lin, C., Xu, A., Zhang, G. and Li, Y. (2013) Polar coordinate
lattice Boltzmann modeling of detonation phenomena. e-print
arXiv: 1308.0653.
[21] Watari, M. (2011) Rotational slip flow in coaxial cylinders by
the finite-difference lattice Boltzmann methods. Communica-
tions in Computational Physics, 9, 1293-1314.
[22] Watari, M. (2010) Relationship between accuracy and number of
velocity particles of the finite-difference lattice Boltzmann
method in velocity slip simulations. Journal of Fluids Engineer-
ing, 132, Article ID: 101401.
[23] Watari, M. and Tsutahara, M. (2003) Two-dimensional thermal
model of the finite-difference lattice Boltzmann method with
high spatial isotropy. Physical Review E, 67, Article ID: 036306.

版权所有：汉斯出版社 (Hans Publishers) Copyright © 2012 Hans Publishers Inc. All rights reserved.