李明鸿 王刚化 莫顺华
贵州贵飞飞机设计研究院有限公司，贵州省安顺市，561000
摘 要 本文应用BGK格式数值模拟高超音速下的钝体绕流问题，研究流场及壁面气动力/热的特征。
本文首先简单介绍BGK格式及其构造方法，然后将其应用到双椭球头高超音速绕流的数值模拟。计算结果与已有研究相吻合，验证了BGK格式的良好性能。
本文对高马赫数与高努森数条件下的流动进行了初步计算分析。当努森数较大时，边界无滑移条件已不适用，可以采用动理论壁面边界条件。
关键词 BGK格式，高超音速流动，钝体绕流
1 前言
近年来，高超音速已经越来越受到人们的关注，同时，高超音速飞行器也在以下方面得到应用：
1.对付时间敏感目标（导弹机动发射车）。即使是“飞毛腿”这种早期的战术弹道导弹，在国际上也能引发大量的政治问题。
2.可重复使用的极高速军用航天器。它能搭载宇航员进入太空，执行维护卫星这一类的航天任务。
3.无需完全进入太空的可重复使用的跨大气层飞行器，绕地球高速高空飞行，进行实时侦查，执行重要情报搜集任务。
由此可见，对于高超音速的研究是非常有必要的。
1.1高超音速流动的特点
高超音速流动一般是指自由流马赫数大于5的流动，具有重要的理论研究和工程应用价值。相比于传统飞行器，高超音速飞行器将面临如下气动关键技术问题[1,2]
1）高速飞行的气动加热导致飞行器表面温急剧升高，从而出现强烈的热力学、化学非平衡效应，给飞行器热防护造成很大的困难。
2）高马赫数与高飞行高度引起飞行器壁面边界层变厚，边界层内出现强烈的粘性相互作用（激波/边界层、激波/激波相互作用等），从而极大地影响壁面压力与热流分布；
3）对临近空间亚轨道飞行器及空天飞机等飞行高度较高的飞行器，低密度大气环境下的稀薄效应也是不容忽视的问题。事实上，对于更低一些的飞行高度，如70km以下，飞行器表面的某些特殊部位，如尖头附近，稀薄效应的影响仍然较大。因此，高超音速飞行器面临的是一个包含连续流、滑移流和部分过渡流的多尺度流动问题。
4）飞行器发射和返回过程中，甚至巡航中飞行器不同部位，气流速度跨越了从亚声速到高超音速这样极大的范围，如何统一求解也是一个亟待解决的问题。
5）此外，由于高超音速飞行器飞行过程中，特别是机动飞行时，大气环境和飞行姿态变化都很大，流场呈现出很强的非定常性，对CFD方法而言，意味着更高的计算精度要求和巨大的计算量。
总之，高超音速飞行器运行过程中将面临从低速到高超声速、从连续流到部分过渡流、高度非定常、复杂物理/化学过程等多种因素耦合在一起的复杂流动问题。准确预测它们对飞行器气动力/热特性的影响直接关系到飞行器的设计，同时也是CFD领域面临的严峻挑战。
2 气动BGK格式
2.1 气动BGK格式的基本原理
本文研究BGK格式在高超音速流动中的应用，这里将简单介绍该方法，细节参见文献[6]。从微观运动出发，利用统计力学的方法描写流体运动，建立关于微观粒子速度分布函数所满足的方程，即Boltzmann方程。对三维问题，在无外加作用力的情况下，方程具有如下形式：
.png)
于是对式（2-2）、（2-5）在相空间内对积分可得
.png)
这里n 代表不同的离散时间步，i, j,k分别代表离散网格单元在x,y,z方向上的编号。对有限体积法，求解点定义于单元的中心，和表示单元的边界。值得注意的是，相比于其他计算格式，如Roe 格式而言，这里的通量既包含了无粘通量，也包含了粘性通量，分别对应粒子的自由运动与相互碰撞。这种两部分通量的耦合使得多维格式的构造变得简单，同时也提高了格式对流场高阶物理量，如热流的模拟精度。
设在单元内为常数，则式（2-2）有形式上的通解
.png)
具体求解时，首先利用前一时间步已知的守恒量构造出单元内的初始分布函数，得到初始平衡态分布函数，利用通解求出单元边界上随时间变化的分布函数，再由该分布函数计算出边界上的守恒量通量，最后代入有限体积法格式得到下一时间步上的守恒量。为减少计算量，本文采用空间分
.png)
可见上述的Pr 数修正基本上不会增加计算量，它可以由前面计算出的守恒量和通量表示出来。
2.5边界条件
BGK 格式的边界条件提法很灵活，既可以同其他差分格式一样直接对宏观量进行边界处理，如无反射边界条件，固壁条件等[5]；也可通过靠近边界的网格点上的物理量构造出的分布函数来计算边界上的守恒量。后者对近连续流动中的固壁边界非常有用，它可以直接得到壁面速度滑移和温度跳跃条件。
2.6 定常流动中的加速收敛技术
对于定常流动计算，为加快计算的收敛，可以采用当地时间步长方法，也就是各网格点按照给定的CFL 数计算出局部时间步长，各单元流场计算按此局部时间步长向前推进。对于网格变形较大的计算，当地时间步长是一种十分有效的方法。另外，也可以采用隐式格式。本文中所涉及的均为定常流动，因此采用结合LU-SGS技术的隐式BGK格式。
3 高超音速双椭球体绕流场计算
3.1 研究意义
双椭球体在外形上与飞行器头部相类似，因此也更加能够通过计算此模型的高超音速下的流场从而解决更多的实际问题。同时，外形复杂度的增高，对网格质量的要求也会随之增高，对于网格的更高要求在这里有所介绍。
3.2 模型构造及网格划分
计算模型的几何参数都是依照[7]来给定。由于主要考虑的是双椭球体上表面的流场，由对称性，如下图1只计算半模：
.png)
在这里我们依然采用的是30×70×20网格（图2），即在X、Y、Z方向上的节点数为30、70、20，在Y方向上网格加密是为了在沿物面法向上贴近壁面处的网格足够密，从而使得结果更精确。
3.3 计算结果
下表1中，是来流的条件以及初始条件：
.png)
.png)
通过以上四幅图片，可以明显的看到在两个椭球的交界处，出现了二次激波。二次激波的出现，使得流场的性质有了新的变化，同时也要求在二次激波出现的地方需要进行网格的加密处理，但由于时间关系，该问题将有待于后续工作来解决。
有图7可以看出，热流密度与文献[6]吻合得不错。在中间部分最小值部分有些许误差主要原因是：由于二次激波处参数会发生巨大变化，所以需要更多的网格来描述激波附近的参数，这样才能精确地计算出结果，所以需要加密交界处上的网格。
.png)
4 结论
在这里通过研究更为复杂的双椭球体高超音速绕流问题，首先是验证了BGK格式对于复杂外形的物体绕流流场的计算的可靠性。其次是对于具有复杂外形的模型的网格划分进行了初步的讨论。同时对于高超音速下的气动力/热的研究也提供了一个高效的方法。
参考文献
[1] 瞿章华，刘伟，曾明，柳军，高超声速空气动力学，长沙：国防科技大学出版社，2001.
[2] K. Xu, A gas-kinetic BGK scheme ofr the Navier-Stokes equations, and its connection with artificial dissipation and Godunov method. J. Comput. Phys. 2001, 171, 289
[3] 李启兵，BGK格式的拓展及应用研究，[清华大学博士后出站报告]. 北京：清华大学，2004
[4] K. Xu, M. Mao, L. Tang, A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow. J. Comput. Phys. 2005, 203, 405
[5] Q.B. Li, S. Fu, K. Xu, Application of gas-kinetic scheme with kinetic boundary conditions in hypersonic flow, AIAA J. 2005, 43(10):2170-2176
[6] 贺立新，间断Galerkin有限元方法及其与有限体积混合计算方法研究，[博士学位论文]. 四川绵阳：中国空气动力研究与发展中心研究生部，2007
[7] 李志辉，从稀薄流到连续流的气体运动论统一数值算法研究，[博士学位论文]. 四川绵阳：中国空气动力研究与发展中心研究生部，2001
作者：李明鸿，1986年5月，男，贵州省都匀市，苗族，本科，工程师，飞机强度/载荷方向。