稀薄气体中高超速圆柱绕流仿真研究

2013-10-21黄海龙

哈尔滨商业大学学报（自然科学版） 2013年4期

黄海龙

(1.哈尔滨工业大学土木学院，哈尔滨 150090;2.北京宇航系统工程研究所，北京 100076)

随着飞行高度逐渐增加，当飞行器高度达到或超过70 km 时，气体的稀薄效应开始明显.该区域内，Kn较大，气体处于非连续区，此时，流体连续介质假设失效，传统的连续流体运动方程(NS 方程)和相应的CFD 方法已经失效.针对这种稀薄气体效应，当采用适当的研究方法，例如采用分子动力学模型和相应的数值方法，其中直接模拟蒙特卡洛方法(DSMC)等方法就是该领域内典型的研究方法.

直接模拟Monte carlo(DSMC)方法是依赖物理的概率模拟方法，在求解过渡领域流动的众多解析、数值和模拟方法中，目前只有DSMC 方法可以模拟三维复杂真实稀薄气体流动［1－6］.直接模拟蒙特卡洛方法可以模拟飞行器穿越高度为70～128 km 范围内的过渡区以及飞行器在70 km 以上高度巡航时的流场稀薄效应，其所涉及的可研究领域包括高超声速巡航飞行器、空间站、卫星、返回舱等多种飞行器.

本文将用连续流体运动方程(NS 方程)和DSMC 方法对四种典型Kn下的圆柱绕流进行计算，分析比较不同Kn时，流场中的间断结构变化情况以及圆柱体后方的分离区变化情况.其中，Kn=0.001 时，流动环境处于连续介质区，可采用连续流体运动方程(NS 方程)求解，其余三种绕流流动采用DSMC 方法研究.

1 DSMC 方法概述

Bird［2］提出的DSMC(Direct Simulation Monte－Carlo)方法是依赖物理的概率模拟方法，它来源于分子动力学方法，采用几率论方法判断分子间是否发生碰撞.DSMC 方法从微观角度出发，利用较少量的模拟分子，代表真实流体的大量分子，用计算机模拟由于气体分子运动、碰撞而引起动量和能量的输运、交换，产生气动力和气动热这一宏观物理过程.在求解过渡领域流动的众多解析、数值和模拟方法中，DSMC 方法是惟一取得巨大成功的方法，从宏观参量到细观速度分布函数的水平上，该方法均能得到实验的支持［3－7］.

DSMC 的模拟过程主要包括两个部分:分子的微观运动和分子碰撞.一般，在实际情况中，分子的运动和碰撞是耦合在一起的，而数值模拟方法中，很难将二者解耦，这给计算过程造成了一定困难.因此，在使用DSMC 方法对稀薄气体流动进行模拟计算时，必须对分子运动和分子碰撞加以特殊对待.在微观领域下，考虑限制DSMC 时间步长小于平均碰撞时间，这时，在平均意义上，单个分子在一个时间步长内最多只有一次碰撞的机会;同时，再要求网格尺寸小于分子平均自由程，则碰撞只会在同一网格中的分子间发生.这样就可在数值计算中，将分子的运动和碰撞实现解耦.在模拟初始，根据物理条件，在计算域内布置模拟分子并给定其位置和速度，在每一个时间步内，分别计算运动和碰撞.处理运动时还包括分子进入或逸出计算域(对应于进出口边界条件)、与壁面的相互作用;碰撞的处理主要是碰撞对的选择、碰撞截面积的计算.流动的宏观量则由网格内模拟分子的运动参数统计平均得到.

2 计算条件和软件设置

2.1 计算条件

本文计算了来流马赫数为5，以圆柱半径为特征长度的Kn分别为0.001、0.01、0.1、1 四种条件下的圆柱绕流.

图1 流动分区示意图

图1 是流动分区示意图.随着Kn的增加，气体流动可分为连续流、滑移流、过渡流和自由分子流四个流区，许多情况下，将滑移流区和过渡流区统称为过渡区，除了连续流区外的其他区域统称为稀薄流区.从图1 可以看出，Kn=0.001 处于连续流区和滑移流区的边界，这时气体比较稠密，连续介质假设仍然成立，因此，此区域内的流体流动常用NS 方程进行计算.而如果用DSMC 计算的话，计算量太大，对CPU和内存的需求都很大，计算时间也很长，计算结果也与求解NS 方程的结果一致，因此，本文采用连续流方程(NS 方程)计算Kn=0.001 时的圆柱绕流流场，可利用商软Fastran 进行仿真.

对于Kn为0.01，0.1、1 的三种情况，本文将采用基于DSMC 方法的程序进行计算.下面将以Kn=0.01为例进行详细介绍.

下面的算例中，取圆柱取直径D=100 mm，来流温度，壁面温度.

根据Kn的定义有:

其中:λ为分子运动的平均自由程;L为流动特征长度，本文中取为圆柱半径D/2.

由可以求得:

DSMC 算法要求，划分网格时，每个网格尺度Δx、Δy、Δz均要小于当地的分子平均自由程λ 的1/3，实际计算中可以放宽到小于λ.

分子运动的平均自由程还有如下关系:

其中:Mref为分子参考质量，对空气而言取Mref=48.1×10－27kg;Dref为分子参考直径，流动介质为空气时，取Dref=4.19×10－10m;ρ为空气密度.

再结合气体分子数密度n 与气体密度ρ 的关系式:

可以得到对应Kn=0.01 下的气体分子数密度n:

分子热运动的平均速度c为:

其中，k=1.3806×10－23J/K为波尔兹曼常数;T为当地温度.代入来流温度T∞，并结合分子平均碰撞时间的定义，可以得到:

DSMC 算法要求，时间步长要小于分子平均碰撞时间Dt 的1/3.

来流温度选为T∞=300K，来流马赫数Ma∞=5，可以得到来流速度为:

2.2 网格划分

本文使用CGNS 格式的结构网格.采用单区或者多区的CGNS 格式的结构网格，同时实现分区并行计算.在计算前，由于流场中分子平均自由程λ的分布未知，因此，需要对流场进行预估，比如在来流部分满足小于来流分子平均自由程λ 的1/3，在激波后，网格要适度加密，在尾流区网格可以适度放宽.

下面给出Kn=0.01 的算例所用的网格在圆柱壁面附近的网格大约尺寸见表1.

表1 网格尺寸(Kn=0.01)

该网格的计算域体积Vcell约为5.5×10－5m3，总网格数Ncell为540 000，因此，每个模拟分子代表的真实分子数Nn为:

在计算过程中，可以根据具体条件对Nn进行适当的调整，比如，当计算条件不足时，需要增大Nn，流场中的模拟分子数就会减少，但此时的计算精度会下降，因此，对Nn的调整需谨慎.

该网格分了四个Block，可以四个线程并行计算.利用MPI 进行分区并行计算，为了提高计算效率，可以对计算域进行合理分区，使每个区的计算量大致相同，这样才不会因为某个区内的计算量过大，导致整体计算速度降低，见图2.

图2 网格

3 结果分析

下面给出Kn分别为0.001、0.01、0.1、1 四种条件下的圆柱绕流的计算结果.其中连续介质区(Kn=0.001)的工况由Fastran 完成，非连续介质区(Kn=0.01、Kn=0.1、Kn=1)的三个工况由DSMC方法计算得到.

3.1 Kn=0.001

图3 给出了Kn=0.001 工况下的马赫数云图和等值线.从图3 可以明显的看到，圆柱前方有清晰的脱体弓形激波，圆柱后缘遇到内折再次产生激波，圆柱正后方有一条狭长的条状低速区，低速区两侧有两条滑移线.图4 给出了Kn=0.001 工况下的压力云图和等值线，弓形激波后有明显的高压区.图5 是圆柱附近的流线图，可以看出，圆柱后方出现了分离，有很大的漩涡区.

图3 马赫数云图(Kn=0.001)

图4 压力云图(Kn=0.001)

图5 圆柱附近流线(Kn=0.001)

3.2 Kn=0.01

图6 给出了Kn=0.01 工况下的马赫数云图和等值线.从图6 中还可以明显的看到圆柱前方有清晰的脱体弓形激波，圆柱后缘遇到内折再次产生激波，圆柱正后方有一条狭长的条状低速区，低速区两侧有两条滑移线.但相较Kn=0.001 的结果，激波和滑移线都明显变粗，弓形激波的位置前移.图7 给出了Kn=0.01 工况下的压力云图和等值线，弓形激波后的高压区较Kn=0.001 的结果减弱.图8 是圆柱附近的流线图，圆柱后方仍存在分离区，但分离区范围减小.

图6 马赫数云图(Kn=0.01)

图7 压力云图(Kn=0.01)

图8 圆柱附近流线(Kn=0.01)

3.3 Kn=0.1

图9 给出了Kn=0.1 工况下的马赫数云图和等值线.从图9 中可以看出，弓形的脱体激波变得更厚，圆柱后缘的激波已经不明显，低速区两侧的两条滑移线已经消失.图10 给出了Kn=0.1 工况下的压力云图和等值线，弓形激波后的高压区进一步减弱.图11 是圆柱附近的流线图，圆柱后已没有分离区，分离区消失.