含转动非平衡效应Boltzmann模型方程统一算法与跨流域绕流问题模拟研究

李志辉，吴俊林，蒋新宇，张涵信

（1.中国空气动力研究与发展中心，绵阳 621000；2.国家计算流体力学实验室，北京 100191）

0 引 言

飞行器往返大气层先后经历连续流、近连续流、滑移过渡流以至高稀薄自由分子流［1］等，不同流域气体热力学性质及绕流状态互不相同，特别是位于连续流与自由分子流之间的过渡区流动，无论在实验技术还是数值方法、理论描述方面均是相当难于处理的流动。在稀薄过渡流区，由于连续介质假设失效，N－S方程已不能有效描述飞行器绕流流动现象；而基于微观分子运动与碰撞随机统计模拟的DSMC方法因受网格划分、时间步长限制，难以对低克努森Kn∞数高空近连续过渡流进行数值仿真［2－4］。如何准确模拟近连续滑移过渡流区高超声速气动力／热环境，一直是人们所关心的问题［4－5］，急需发展相关基础理论与数值计算方法。

气体分子运动论（气体动理学理论）的基本方程——Boltzmann（玻尔兹曼）方程［6］本身可描述各个流域的气体分子输运现象，通过求解Boltzmann方程可以研究各流域气体动力学问题。但该方程是一个高度非线性积分－微分方程，除Maxwell分布等少数几个解析解外，几乎不可能求出精确解。为避免求解复杂碰撞项所带来的困难，众多学者基于守恒不变量，由数学上较简单的统计和碰撞松弛模型代替Boltzmann方程碰撞项，提出许多类似Boltzmann方程各阶矩的气体运动论模型方程。最简单的BGK模型方程于1954年提出［7］，随后，人们发展了多种形式的Boltzmann模型方程，诸如 Holway的椭球统计模型［8］、Shakhov基于BGK方程的修正而得到的高阶推广模型方程［9］等。因此，求解从Boltzmann方程简化而得的气体运动论模型方程，是一个较为经济有效的解决稀薄气体流动问题的途径。基于文献［7－10］基础，我们［11－15］从研究描述各流域均适用的气体分子速度分布函数方程出发，吸收计算数学指数型积分求解原理，提出发展气体分子运动论离散速度坐标法，研制经改进的Gauss－Hermite积分方法等系列离散速度数值积分技术，消除原分布函数对速度空间的连续依赖性，将Boltzmann模型方程化为在各个离散速度坐标点处基于时间、位置空间具有非线性源项的双曲型守恒方程，利用二阶Runge－Kutta方法数值求解碰撞松弛源项方程，引入NND格式数值求解对流运动项。发展可用于气体分子速度分布函数有限差分求解的气体运动论边界条件数学模型并进行数值处理。由此建立求解一维、二维、三维Boltzmann模型方程的气体运动论统一算法（Gas－Kinetic Unified Algorithm，GKUA），并开展跨流域气体绕流问题模拟应用［16－18］。为了开展近空间飞行器跨流域非平衡绕流问题研究，近年来，基于对转动自由度松弛变化［19］、气体分子运动论 Rykov模型［20－21］研究，在气体分子速度分布函数演化更新求解中考虑转动自由度影响，把气体分子转动能作为分布函数的自变量，在统一算法理论计算框架下，提出考虑转动非平衡效应Boltzmann模型方程数值算法［22］。本文在此基础上，针对高、低不同马赫数（1.5≤Ms≤25）氮气一维激波结构和全局 Knudsen数（9×10－4≤Kn∞≤10）竖直平板、Ramp制动器跨流域绕流、尖双锥再入体绕流算例，进一步研究含转动非平衡效应Boltzmann模型方程在全飞行流域绕流问题中的模拟研究，计算剖析飞行器从高稀薄自由分子流到连续流跨越中间过渡带非平衡流动变化规律与复杂绕流现象。

1 控制方程及计算方法

基于转动松弛特性 Rykov模型［20－21］，采用转动惯量来描述气体分子的自旋运动，利用分子总角动量守恒作为一个新的碰撞不变量。基于气体分子速度分布函数f＝f（r，ξ，t，e），这里r、ξ分别是位置空间和速度空间的坐标，e为内能，在求解Boltzmann模型方程统一算法框架［11－18］下，基于权因子1和e对速度分布函数进行无穷积分，引入约化速度分布函数f0与f1，发展含转动非平衡效应的Boltzmann模型方程，其无量纲形式为：

宏观流动参数可以由f0和f1对速度空间积分求得：

粘性系数μ和碰撞次数Z可定义为：

这里，T∞和T＊均为有量纲的值，但B＝T∞／T＊是无量纲量［23］。对于氮气N2而言，T＊＝91.5K，1／δ＝1.55。

速度分布函数f0、f1是依赖于时间t、位置空间r和速度空间ξ的函数，需将速度空间离散降维，去掉分布函数对速度空间的连续依赖性。气体运动论离散速度坐标法［10－11，13－14，24］通过一套与积分规则相一致的离散速度坐标点分布，用具有N个元素的函数簇代替原分布函数对积分变量的连续依赖性，把模型方程转换为在各个离散速度坐标点处彼此独立的关于位置空间和时间的双曲型守恒方程，而宏观流动参数则可以通过相应于离散速度坐标法的积分求解规则动态确定。如经过速度空间离散化处理的速度分布函数方程在计算平面内矩阵形式为：

其中，U＝JF，E＝·U，G＝·U，H＝·U，

采用文献［11，13－14，24］所建立二阶有限差分格式，将时间分裂数值计算方法应用于方程（4）中位置空间的三个对流运动方程和 碰撞松弛源项方程数值求解，得到在每个离散速度坐标点处直接求解分子速度分布函数的气体运动论数值格式：

2 边界条件数值处理方法

对于流场边界条件数值处理，按照统一算法基于分子速度分布函数边界处理数学模型［14，17，24］，对于特征速度小于零的无穷远处来流条件，我们以无量纲的来流宏观参数n∞＝1、U∞、V∞、W∞和T∞＝1表征的平衡态分布函数确定：

而在流场出口边界上，认为分布函数趋于均匀变化而不存在变化梯度，即：

物面边界条件由气体分子与物体表面碰撞时，遵循无穿透通量守恒条件［13－15，17］和平动、转动能量平衡关系，推导得到：

其中，nw为无穿透通量守恒边界条件确定的物面反射气体分子数密度，

物体表面气体的平动温度和转动温度可表示为：

公式（10）、（11）中的适应系数αt、、αr和可以有两种赋值方法。第一种，对应于完全热适应边界条件：αt＝1，＝0，αr＝1，＝0。第二种，适应系数通过分析氮气中热传导实验数据得到［25］：αt＝0.8，＝0.1，αr＝0.5＝0.8。

3 数值计算结果与讨论

3.1 高、低不同马赫数定常正激波内流动

正激波相当于一个从均匀超声速上游到一个均匀亚声速下游流动过渡［26］。引入无量纲参数：

为检验考虑转动非平衡影响的Boltzmann模型方程统一算法对高、低不同马赫数激波内流动问题模拟能力，计算氮气中Ms＝1.53、3.2、10、25激波内流动。图1分别绘出马赫数Ms＝1.53、3.2、10、25激波结构无量纲密度ρ＊分布，图中实线“Cal.”表示统一算法结果，符号“○”为来自文献［27］实验数据，图1（a）、（b）中符号“▽”为来自文献［21］模拟值，图1（c）、（d）中符号“△”为来自文献［19］使用广义退化Boltzmann方程（GBE）计算结果。可看出，对于高、低不同马赫数1.53～25激波内流动密度分布统一算法结果均与实验数据、典型文献模拟值吻合很好。图1关于弱激波Ms＝1.53到强激波Ms＝25的计算，证实本文方法在计算分析双原子氮气高、低不同马赫数激波内流动问题方面的准确可靠性。图2给出了Ms＝7激波内流动平动温度和转动温度分布，其中符号“Ref.”为来自文献［21］模拟值，图2（b）标示出激波转动温度实验测试数据［28］。可看出，从整个激波结构的变化情况来看，不仅平动温度而且转动温度分布本文计算值均与考虑转动能自由度影响的Rykov模型结果相一致，且与转动温度分布实验测试数据相吻合，验证了统一算法程序求解考虑转动非平衡影响的一维流动问题可行性。

图1 高低不同马赫数下激波内流动密度分布计算比较Fig.1 Density distribution in shock wave of various Mach numbers

图2 Ms＝7激波结构内流动温度分布计算比较Fig.2 Temperature profiles of shock wave at Ms＝7

3.2 二维再入制动器全飞行流域绕流计算分析

为了检验与应用本文发展的统一算法计算再入不同Knudsen数二维高超声速绕流问题可靠性，拟定竖直平板绕流试验状态［29］：来流马赫数Ma∞＝8.3666，壁温与总温比Tw／T0＝0.6667，γ＝1.4。

表1绘出利用壁面完全热适应边界条件计算Kn∞＝0.1、0.01竖直平板绕流阻力与热传导系数与文献［29］数据比较，可看出，在高稀薄流Kn∞＝0.1与近连续流Kn∞＝0.01跨越不同流区上述竖直平板绕流物面力／热系数GKUA结果均与典型文献吻合很好，二者最大偏差不超过4.9%。

图3、图4分别绘出Kn∞＝0.01平板绕流流场平动温度Tt、转动温度Tr等值线分布与文献［29］数据比较，可看出两种结果变化趋势相当吻合。对此近连续流场，平板绕流前端的减速区内出现明晰的脱体激波，由于Knudsen数不是小量，流场稀薄效应明显存在，平板前端脱体激波较厚。图5描绘了该平板绕流Kn∞＝0.01、Ma∞＝8.3666流线分布，给出了气体绕流方向与流动变化特点，并与图5（b）显示的文献［29］中的流线图进行了对比。结果表明，GKUA流线分布（图5a）与文献［29］计算结果之间具有较好的一致性。对此近连续绕流Kn∞＝0.01状态，气体流经平板，在物体前部一定区域因跨越脱体激波而流向改变，绕过物体最高点，进入背风区，出现流动分离、回流涡区。从图3～图5统一算法计算得到的流场结构与文献［29］结果定性化比较看出，GKUA结果流场分辨率较好，加之表1物面力／热系数定量化比较，进一步证实本文考虑转动非平衡效应的Boltzmann模型方程统一算法计算再入不同流区高超声速流动问题准确可行性。

表1 跨流区竖直平板绕流阻力和热传导系数比较Table 1 Comparison of drag and heat transfer coefficients past erectplate covering various flow regimes

图3 平板绕流Kn∞＝0.01平动温度Tt等值线分布计算比较Fig.3 Translational temperature contours in erect plate flow with Kn∞＝0.01

图4 平板绕流Kn∞＝0.01转动温度Tr等值线分布计算比较Fig.4 Rotational temperature contours in erect plate flow with Kn∞＝0.01

图5 平板绕流Kn∞＝0.01、Ma∞＝8.3666流场流线分布Fig.5 Streamline structures past erect plate with Kn∞ ＝0.01

针对如图6所示的一个平板与25°倾角斜板构成二维Ramp再入制动器外形，计算Ma∞＝10的高马赫数绕流流场。设置：克努森数Kn∞＝6×10－4、来流气体温度T∞＝52K、密度ρ∞＝3.91×10－4kg／m3、速度V∞＝1477m／s、壁面温度Tw＝290K。

图6 Ramp制动器外形及尺寸（单位：cm）Fig.6 Ramp body geometries（unit：cm）

本文计算得到该Ramp制动器绕流流场马赫数等值线与来自文献［30］的Boltzmann－ES椭球统计模型结果比较，见图7所示，可看出二者吻合较好，流场变化规律较为一致，计算再现了上表面水平板前端尖前缘产生斜激波与斜板上部压缩斜激波相互干扰绕流现象，得到了较高分辨率的流场计算结果，验证了本文方法程序计算的正确性。

该状态下Ramp外形绕流上部流场x＝0.1m和x＝0.2m截线处气流密度、温度与马赫数随z轴横向分布本文计算与文献［30］模拟结果比较情况，见图8、图9所示。可看出在流场不同位置本文计算值与文献［30］椭球统计模型结果变化趋势吻合较好，揭示了Ramp凹面体外形前、后部流场出现不同强度压缩斜激波流动变化规律，验证了含转动非平衡效应Boltzmann模型方程数值算法对二维高超声速连续流区Kn∞＝6×10－4、Ma∞＝10绕流问题计算适应性。

图8 Ramp上部流场x＝0.1m截线处密度、温度与马赫数横向分布Fig.8 Density（ρ），temperature（T）and Mach（M）distribution in the upside field of Ramp at x＝0.1m

图9 Ramp上部流场x＝0.2m截线处密度、温度与马赫数横向分布Fig.9 Density（ρ），temperature（T）and Mach（M）distribution in the upside field of Ramp at x＝0.2m

图10绘出该二维Ramp制动器上表面压力系数Cp和摩阻系数Cf沿x轴线分布本文计算值（Cal.）与文献［30］中Boltzmann－ES模型结果（Ref.）比较情况，可看出两者相当吻合，证实所发展含转动非平衡效应Boltzmann模型方程统一算法对飞行器物面气动特性计算的可靠性与精度。

图10 Ramp上表面压力、摩阻系数分布计算比较Fig.10 Distribution of pressure and frictional drag coefficients on the upper surface of Ramp

为了检验求解转动非平衡效应Boltzmann模型方程统一算法对二维Ramp制动器外形再入各流域绕流问题计算能力，拟定高稀薄自由分子流到连续介质绕流状态Kn∞＝10～0.0009、Ma∞＝1.8进行计算分析。图11给出了跨流区Kn∞＝10、0.8、0.01、0.0009四个不同克努森数下计算得到Ramp制动器再入绕流流场尾部矢量流线与流动局部结构。可看出，对于Kn∞＝10的高稀薄自由分子流，气流完全附着物面流动；而当流动进入稀薄过渡流，如Kn∞＝0.8，在飞行器尾部上端尖角后部绕流区开始出现较小的回流涡结构；对于Kn∞＝0.01的近连续流区，飞行器尾部绕流出现双涡结构，Kn∞＝0.0009的连续流区，尾部流场又褪变为存在一个巨大的回流涡流动结构，该图直观展示了二维制动器再入过程跨越飞行各流域所出现的不同绕流变化规律。

图11 Ramp制动器再入过程Kn∞＝10、0.8、0.01、0.0009绕流结构（Ma＝1.8，CFL＝2）Fig.11 Flow structures of Ramp re－entering for all flow regimes with Kn∞ ＝10，0.8，0.01and 0.0009（Ma＝1.8，CFL＝2）

3.3 双维外形体再入稀薄流区绕流计算分析

为了检验考虑转动非平衡影响的Boltzmann模型方程统一算法对三维飞行器再入绕流问题的模拟能力，拟定双锥再入体外形及尺寸如图12所示，为了与适于高稀薄流模拟的GBE解算器［31］计算比较，拟定来流马赫数Ma∞＝3、Kn∞＝1尖双锥绕流状态，图13绘出统一算法计算该状态驻点线密度分布与文献［31］结果（符号“Δ”）对比情况，可看出两种方法计算结果很好一致，证实考虑转动非平衡效应的统一算法用于计算三维再入飞行器绕流问题可行性。

图14绘出Ma∞＝3绕流Kn∞＝1和Kn∞＝0.1尖双锥流场密度等值线分布，图中看出，这两种状态克努森数较高，流动均呈现较强的稀薄流特征，不存在明显的激波，流动参数分布没有强间断现象，而且图14（a）对应Kn∞＝1的高稀薄绕流物体前方扰动区域较图14（b）对应更低克努森数Kn∞＝0.1的稀薄过渡区绕流要大得多；Kn∞＝0.1对应的气体流动趋近物体前端出现了激波强扰动现象。

图12 双锥再入体外形Fig.12 Bicone geometries

图13 驻点线密度分布Fig.13 Density distribution along the stagnation line for Ma∞＝3，Kn∞＝1

图14 不同克努森数尖双锥Ma＝3绕流场密度分布Fig.14 Density distribution around bicone with different Kn

图15直观显示Kn∞＝1和Kn∞＝0.1两种状态下尖双锥绕流流线，反映出两个状态的绕流面貌呈现明显不同的流动结构，Kn∞＝1对应的高稀薄近自由分子流，气流完全附着物面流动过锥体后，平滑过渡到背风下游流动状态，呈现明显的高稀薄流特征；而Kn∞＝0.1对应的更低克努森数稀薄过渡区绕流，则在锥体后端面背风区开始出现流动分离，形成上下两个回流涡，鞍结点拓扑结构，开始呈现出近连续流绕流面貌。

图15 不同克努森数尖双锥绕流（Ma∞＝3）流线结构Fig.15 Streamline distribution around bicone at different Kn with Ma∞＝3

为了剖析跨流区流动的非平衡效应，图16、图17分别绘出了Kn∞＝1、0.1两种状态对应的绕流流场温度分布。在Kn∞＝1对应的高稀薄流区，平动温度与转动温度分布彼此差别很大，计算得到的物体前部绕流区平动温度最大值约为3.2，而转动温度最大为2.4左右；对应Kn∞＝0.1的稀薄过渡区绕流，平动温度与转动温度之间的差异明显减小，物体前部高温区平动温度最大值下降为2.8，转动温度升高约为2.65；而从图3、图4所绘竖直平板绕流看出，对于Kn∞＝0.01的近连续绕流，流场平动温度与转动温度分布趋于一致，流场等值线结构差别不大，流动属于近平衡态流动。这体现了高Kn∞数稀薄过渡流区存在严重的流动非平衡，且克努森数越大非平衡效应越明显，说明了在稀薄流区考虑转动非平衡效应的必要性。

图16 Ma∞＝3，Kn∞＝1尖双锥绕流温度分布Fig.16 Temperature contours around bicone at Ma∞＝3，Kn∞＝1

图17 Ma∞＝3，Kn∞＝0.1尖双锥绕流温度分布Fig.17 Temperature contours around bicone at Ma∞＝3，Kn∞＝0.1

4 结 论

基于对转动自由度松弛变化、气体分子运动论Rykov模型研究，应用稀薄流到连续流的气体运动论统一算法数值求解考虑转动非平衡影响的Boltzmann模型方程。针对高、低不同马赫数定常正激波结构、从自由分子流到连续流二维竖直平板、Ramp制动器及三维尖双锥外形体再入稀薄流区绕流问题计算分析，得到如下结论：

（1）通过对弱激波到强激波1.53≤Ms≤25考虑转动非平衡影响的氮气激波结构内流动计算结果与实验数据、典型文献模拟值比较分析，验证了含转动非平衡效应Boltzmann模型方程统一算法对求解强、弱激波一维流动问题的准确可靠性。

（2）通过对二维竖直平板、Ramp再入制动器从自由分子流到连续流跨流域绕流流场转动非平衡效应模拟研究，验证了统一算法对全飞行流域模拟计算的一致适用性。

（3）通过对三维尖双锥外形再入稀薄气体绕流状态数值计算，对比分析绕流流场平动温度与转动温度分布，揭示了来流Kn∞数越高，过渡区非平衡效应越严重，飞行器绕流流动非平衡影响越明显。证实考虑转动非平衡影响的Boltzmann模型方程统一算法对三维复杂飞行器再入绕流问题模拟能力。

本文工作仅是含转动非平衡效应Boltzmann模型方程统一算法在跨流域绕流问题方面的初步应用与计算检验，将其应用于跨流区复杂高超声速气动热力学问题以及考虑转动能与振动能、混合气体非平衡效应适应性研究，均有待进一步研究。

致谢：本文部分三维绕流计算在国家超级计算天津中心进行，特此感谢。

［1］TSIEN H S.Superaerodynamics，mechanics of rarefied gases［J］.JournalofAeronauticsScience，1946，13（12）：653－664.

［2］KOPPENWALLNER G，LEGGE H.Drag of bodies in rarefied hypersonic flow［R］.AIAA Paper 85－0998.Progress in Astronautics and Aeronautics：Thermophysical Aspects of Reentry Flows.New York，1985，103：44－59.

［3］LI Z H，WU Z Y.DSMC simulation on rarefied aerodynamics of Apollo－CM［J］.ACTAAerodynamicaSinica，1996，14（2）：230－233.（in Chinese）李志辉，吴振宇.阿波罗指令舱稀薄气体动力学特征的蒙特卡罗数值模拟［J］.空气动力学学报，1996，14（2）：230－233.

［4］YOON S，GNOFFO P A，WHITE J A，et al.Computational challenges in hypersonic flow simulations［R］.AIAA 2007－4265.

［5］ZHUANG F G，CUI E J，ZHANG H X.Some development of future spacecrafts and aerodynamics tasks［C］.Proc.of First Aerodynamics and Aerothermodynamics，2006，8：1－12.（in Chinese）庄逢甘，崔尔杰，张涵信.未来空间飞行器的某些发展和空气动力学的任务［C］.中国第一届近代空气动力学与气动热力学会议论文集，2006，8：1－12.

［6］CHAPMANN S，COWLING T G.The mathematical theory of non－uniform hases［M］.3rded.Cambridge，1990.

［7］BHATNAGAR P L，GROSS E P，KROOK M.A model collision processes in gases［J］.Phys.Rev.，1954，94：511－525.

［8］HOLWAY Jr L H.New statistical models for kinetic theory，methods of construction［J］.Phys.Fluids，1966，9：1658－1673.

［9］SHAKHOV E M.Kinetic model equations and numerical results［M］.In：Oguchi H Ed.Proceedings of 14thInternational Symposium on Rarefied Gas Dynamics，Tokyo：University of Tokyo Press，1984，pp137－148.

［10］YANG J Y，HUANG J C.Rarefied flow computations using nonlinear model Boltzmann equations［J］.J.ofComput.Phys.，1995，120：323－339.

［11］LI Z H，ZHANG H X.Pilot study on gas－kinetic unified numerical algorithm for flows from rarefied transition to continuum［J］.ACTAAerodynamicaSinica，2000，18（3）：251－259.（in Chinese）李志辉，张涵信.稀薄流到连续流的气体运动论统一数值算法初步研究［J］.空气动力学学报，2000，18（3）：251－259.

［12］LI Z H，ZHANG H X.Study on gas kinetic unified algorithm for flows from rarefied to continuum［J］.ACTAAerodynamica Sinica，2003，21（3）：255－266.（in Chinese）李志辉，张涵信.稀薄流到连续流的气体运动论统一算法研究［J］.空气动力学学报，2003，21（3）：255－266.

［13］LI Z H，ZHANG H X.Study on gas kinetic algorithm for flows from rarefied transition to continuum using boltzmann model equation［J］.ActaMechanicaSinica，2002，34（2）：145－155.（in Chinese）李志辉，张涵信.稀薄流到连续流的气体运动论模型方程算法研究［J］.力学学报，2002，34（2）：145－155.

［14］LI Z H，ZHANG H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］.J.ofComput.Phys.，2004，193：708－738.

［15］LI Z H.Gas－kinetic unified algorithm for re－entering complex problems covering various flow regimes by solving Boltzmann model equation［M］.Advances in Spacecraft Technologies.In－Tech Publisher，2011：273－332.

［16］LI Z H，ZHANG H X，FU S，et al.A gas kinetic algorithm for flows in microchannel［J］.InternationalJournalofNonlinear SciencesandNumericalSimulation，2005，6（3）：261－270.

［17］LI Z H，ZHANG H X.Gas－kinetic numerical studies of threedimensional complex flows on spacecraft re－entry［J］.J.Comput.Phys.，2009，228：1116－1138.

［18］LI Z H，ZHANG H X.Study on the unified algorithm for threedimensional complex problems covering various flow regimes using Boltzmann model equation［J］.Sci.China－Phys.Mech.＆Astron.，2009，52（1）：124－138.（in Chinese）李志辉，张涵信.基于Boltzmann模型方程各流域三维复杂绕流问题统一算法研究［J］.中国科学，2009，39（3）：414－427.

［19］AGRWAL R K，CHEN R，CHEREMISIN F G.Computation of hypersonic shock wave flows of a diatomic using generalized Boltzmann equation［R］.AIAA 2007－4541.39th AIAA Thermophysics Conference.Miami，2007.

［20］RYKOV V A.Model kinetic equation of a gas with rotational degrees of freedom［J］.FluidDynamics，1975，10：959－966.

［21］RYKOV V A，TITAREV V A，SHAKHOV E M.Shock wave structure in a diatomic gas based on a kinetic model［J］.FluidDynamics，2008，43（2）：316－326.

［22］LI Z H，WU J L，PENG A P，ZHANG H X.Study of the unified algorithm solving Boltzmann model equation with the effect of rotational non－equilibrium［C］.The 8th Cross－straits Conference on Aeronautics and Astronautics.Beijing，2012.（in Chinese）李志辉，吴俊林，彭傲平，等.考虑转动非平衡效应Boltzmann模型方程统一算法应用研究［C］.第八届海峡两岸航空航天学术研讨会论文集，北京，2012.

［23］LARINA I N，RYKOV V A.Diatomic gas flow past a sphere based on kinetic equations［J］.Dokl.Akad.NaukSSSR，1976，227：60－62.

［24］LI Z H.Study on the unified algorithm from rarefied flow to continuum［D］.［PhD thesis］.Mianyang：China Aerodynamics Research and Development Center，2001.（in Chinese）李志辉.从稀薄流到连续流的气体运动论统一数值算法研究［D］.［博士学位论文］.中国空气动力研究与发展中心，2001.

［25］LARINA I N，RYKOV V A.The boundary condition on the body surface for a diatomic gas［C］.Proceedings of 15thInternational Symposium on Rarefied Gas Dynamics，Stuttgart，1986，1：635－643.

［26］LI Z H.Applications of gas kinetic unified algorithm for flows from rarefied transition to continuum［D］.［Post－doctor disser－tation］.Beijing：Tsinghua University，2003.（in Chinese）李志辉.稀薄流到连续流气体流动问题统一算法应用研究［D］.［博士后研究报告］.清华大学，2003，9.

［27］ALSMEYER H.Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam［J］.J.Fluid.Mech，1976，74：497－513.

［28］ROBBEN F，TALBOT L.Experimental study of the rotational distribution function of nitrogen in a shock wave［J］.Physicsof Fluids，1966，9：653－662.

［29］RYKOV V A，TITAREV V A，SHAKHOV E M.Numerical study of the transverse supersonic flow of a diatomic rarefied gas past a plate［J］.ComputationalMathematicsandMathematical Physics，2007，47（1）：136－150.

［30］ANDRIES P，BOURGAT J F.Numerical comparison between the Boltzmann and ES－BGK models for rarefied gases［J］.ComputerMethodsinAppliedMechanicsandEngineering，2002，191（31）：3369－3390.

［31］WILSON C D，AGARWAL R K，TCHEREMISSINE F G.Computation of hyersonic flow of a diatomic gas in rotationa nonequilibrium past 3Dblunt bodies using the generalized Boltzmann equation［R］.AIAA 2009－3836.