摘要：研究了二维非结构网格高超声速稀薄流DSMC方法中气体与壁面相互作用的数值模型。发展了一种基于辐射平衡的壁面温度边界条件，与恒温壁面边界条件相比，该温度边界条件可以克服恒温壁面边界条件的自身缺陷并且能够在流场模拟中适时给出更贴近真实情况的壁面温度；采用了由完全漫反射和纯镜面反射模型组合而成的Maxwell壁面反射模型，该壁面反射模型能够更好地描述气体模拟分子在壁面散射的真实情况。采用本文方法对钝头体外形进行了数值模拟，结果表明本文采用的气体壁面相互作用模型能够提高钝头体背风面后部流场对温度变化的敏感度，并且随着飞行高度提升，敏感度有减弱的趋势。
关键词：辐射平衡壁面边界；Maxwell壁面模型；DSMC
中图分类号：O356 文献标志码：A 文章编号：1005-2615（2015）03-0348-07
Gas-Surface Interaction Models for DSMC Method
Qu Cheng，Wang Jiangfeng
Abstract：Gas-surface interaction models in direct simulation Monte Carlo （DSMC） method for hypersonic rarefied flow is studied based on unstructured grids. A wall boundary condition based on radiative equilibrium is developed， which can overcome the defects of isothermal boundary condition and provide more precise wall temperature distribution during numercial simulations. The Maxwell reflection boundary model， coupled with diffuse reflection and specular reflection， is applied. Numerical results of a blunt body configuration arc obtained. Test cases show that the developed models can improve the sensitivity of the temperature variation in afterbody field， and the sensitivity decreases with the simulated flying altitude increases.
Key words： radiative equilibrium boundary； Maxwell boundary model； direct simulation Monte Carlo
到目前为止，人类在过渡领域的飞行和风洞试验充满危险并且代价高昂，数值模拟虽然不能完全替代物理试验，却能够极大降低此类试验的数量，DSMC（Direct simulatlon Monte Carlo）方法是Bird直接从流动的物理模拟出发发展的一种直接模拟蒙特卡洛方法[1]。在模拟过渡领域三维真实气体流动问题方面成功的案例很多，例如美国国家航天局（NASA）利用DSMC方法成功模拟了AFE飞船、欧洲Hermes航天飞机DSMC模拟结果与飞行测量数据的出色相符、中国长征-4B火箭末子级推进剂在轨排放扰动力矩的成功预测等[2-4]。近年来，国内关于DSMC的研究工作主要集中于并行计算技术、自适应时间步长[5]以及羽流应用[6-7]等方面。在气体与壁面相互作用模型的研究方面，贺群武[8]等人在逆温度抽样算法的基础上发展了一种DSMC方法中壁面处给定热流边界条件的实施方法，黄飞[9]等人采用一种统计方法实现了DSMC方法在热力学非平衡碰撞时绝热壁条件的抽样。
基金项目：“十二五”国防基础科研资助项目。
收稿日期：2013-11-04；修订日期：2014-12-09
DSMC方法中通常采用的是恒温壁面边界，气体模拟分子与壁面作用后散射进入流场的能量信息由初始设定的壁面温度给出[5]，但是，这种壁面温度条件会使流场中初始与壁面作用散射出的模拟分子能量过高[10]，与实际情况并不相符。气体模拟分子与壁面相互作用最常采用的是纯镜面反射模型和完全漫反射模型，但是试验研究证明，无论是纯镜面反射模型还是完全漫反射模型都不能很好地描述分子在物体表面散射的真实情况[3]。针对上述问题，本文主要研究了稀薄流问题中的气体模拟分子与壁而相互作用的数值模型。首先，发展了一种基于辐射平衡的壁面温度边界条件，在数值模拟中能够适时给出更贴近真实情况的壁面温度。随后，采用由完全漫反射和纯镜面反射模型组合而成的壁面反射模型，控制完全漫反射和镜面反射分子的比例，从而能更好地描述模拟分子在壁面散射的真实情况。最后，编写了非结构网格二维DSMC程序，对过渡流域钝头体外形再入绕流流场进行了数值模拟与分析。
1 计算方法
1.1 DSMC方法介绍
DSMC法使用大量的模拟分子模拟真实的气体，用一个模拟分子代表一定数量的真实气体分子，该方法的本质是在时间步长△t内，对分子的运动和分子间的碰撞进行解耦。整个模拟中分子之间以及分子和物面之间不断进行能量交换，在充足的模拟时间后，通过采样统计得到每个网格的宏观流场结果。DSMC方法中，计算网格起两种作用[11]。第一种是对流场宏观流动参数进行空间离敞，第二种是促进碰撞分子碰撞对的选择，使其满足基本的几何接近。本文采用非结构贴体网格，计算中采用可变硬球（VHS）分子模型，利用Bird非时间计数器（NTC）法进行碰撞对的选取，使用Larsen-Borgnakke唯象论模型处理模拟分子间的能量交换[5]。 1.2 壁面边界条件
1.2.1 恒温壁面边界条件
在DSMC方法中，根据流场初始给定的来流密度和速度，恒温壁面边界中的初始壁面温度可以通过求解斯特藩定律（Stefan-Boltzmann law）[12]得出
式中：ε为表面发射系数，σSB为斯特藩-玻耳兹曼常数。在实际问题中，一般假设壁面为黑体表面，令ε=1。
1.2.2 辐射平衡壁面边界条件
恒温壁面边界条件假定物面温度初始给定且保持不变，然而，DSMC方法模拟中时间参数与真实流动中的物理时间等同，所有计算都是非定常的[1]，实际壁面温度会随流场进程发生改变，因此，恒温壁面边界并不能很好地符合真实情况。文献[12]在处理连续流高超声速非定常流动问题时采用了基于Stefan-Boltzmann law定律及热传导的Fourier定律的辐射壁温度边界条件，通过Newton迭代适时给出流场中的壁面温度.取得了良好的计算结果。
本文对文献[12]中连续流辐射壁温度边界条件进行了推广，发展了稀薄流DSMC方法基于辐射平衡的壁面边界条件，该边界条件能够克服恒温壁面边界条件的自身缺陷并且在流场模拟中适时给出更贴近真实情况的壁面温度，主要思想如下：
（1）根据给定的流场计算条件.将初始壁面温度设为利用斯特藩定律得到的温度。
（2）假设壁面上编号为i的单元表面的热流密度是根据斯特藩定律完全由该单元表面辐射的，这样，通过第n时间步的热流密度就可以得到第n+1时间步的壁面温度，计算公式如下
式中热流密度为单位时间内单位面积上入射分子与反射分子的总能量之差。即
式中：下标in，out分别表示入射分子和反射分子；Q表示每个分子携带的能量，包含平动能和转动能；△tn表示循环到第n 步的作用时间；Si表爪编号为i的单元表面积。
通过以上处理方法，在程序中不断更新飞行器的壁面温度，在实际处理中，认为飞行器的舱内温度为300 K，当通过式（2）计算得到的温度低于300 K时，将壁面温度设为舱内温度。
1.3 壁面反射模型
纯镜面反射模型实现相对容易，但是不能形成壁面附近实际存在的附面层流动，完全漫反射模型则相对复杂，也能够模拟附面层流动，但是该模型依然不能很好地满足真实情况，为了较好地描述模拟分子在壁面散射的真实情况，本文采用由纯镜面反射模型和完全漫反射模型组合而成的Maxwell壁面反射模型，假设模拟分子在壁面发生完全漫反射和纯镜面反射的比例为α：（1―α），实际问题中通常控制α=0. 8～0.9[13]，在本文中，取α=0.85。
本文对纯镜面反射模型和完全漫反射模型的处理方法如下：
（1）镜面反射模型是一种弹性碰撞模型，假定来流分子在物体表面的反射与光滑弹性球在光滑的完全弹性表面上的反射相同，即模拟分子与壁面碰撞后，其速度平行于壁面的切向分量保持不变，而法向速度分量大小不变，方向与碰撞前相反。
（2）完全漫反射模型是一种非弹性碰撞模型[3]，假定离开物面的分子以平衡的Maxwell分布散射，平衡条件是表面温度，Maxwell分布的温度与来流的静温相同。漫反射分子的分布函数满足
1.4 远场边界的处理
高超声速稀薄流数值求解问题中遇到的流动一般都不是封闭的，在远场边界处会不断有新的分子进入流场[14]。如果来流速度C0与远场边界网格单元的外法矢成一定的夹角θ，那么，每个时间步长从该单元进入流场的模拟分子数可以表示为
式中：n∞为自由来流的分子数密度；为来流分子的最概然速度；Si为远场边界网格i的来流面面积；△t为时间步长；f为一个模拟分子代表的真实分子数；s为分子速度比，是来流宏观速度c0与之比，即；erf为误差函数。crf（α）定义为
在数值模拟时，计算得到的Ni是实数，实际每个时间步长进入的分子数是整数，本文在程序中将截断部分的分子个数保存起来，加到下一时间步长进入边界网格的分子数中。
2 数值算例及结果分析
2.1 计算条件
本文用非结构网格DSMC方法对过渡流域钝头体外形高超声速绕流流场进行了数值模拟。在每一高度下来流气体的组分仅考虑了N2和O2，组分摩尔百分比近似取为N2：O2=0.763：0.237，不考虑来流气体组分、摩尔比值等随高度的变化。来流速度取为V∞=7.5 km/s，钝头体几何及部分非结构网格如图1所示，表1给出了本文数值模拟的计算状态及相应的气体壁面相互作用模型。
2.2 结果比较与分析
图2给出了计算状态1和计算状态2的流场总温分布等值线及云图，观察图2（a）和图2（b）可以看出非常明显的弓形激波，激波层后流场温度升高明显，在激波过渡区内出现了流场最高温度，钝头体迎风面前部流场温度等值线及云图分布基本一致。对比图2（a）和图2（b）可以发现，图2（b）钝头体背区域温度明显低于图2（a）钝头体背风区域温度，图2（b）采用的气体壁面相互作用模型对温度的变化更为敏感，这可能是因为从壁面散射出去的分子能量信息是由壁面温度给出，在此计算高度下，采用恒温壁面条件会使得初始与钝头体壁面发生碰撞散射进入流场的模拟分子带有过高能量，而
计算状态2采用的壁面条件降低了初始设定的壁面温度，并且假设辐射能从流场中直接耗散而不是被散射模拟分子吸收。
图3给出了计算状态3和计算状态4的流场总温分布等值线及云图，观察图3（a）和图3（b）同样可以发现非常明显的弓形激波，钝头体迎风面前部流场温度等值线及云图分布规律与图2（a）和图2（b）相符，对比图3（a）和图3（b）可以发现图3（b）钝头体背风面温度明显低于图3（a）钝头体背风面温度，同样与对比图2（a）和图2（b）观察到的现象相符。 图4给出了计算状态9和计算状态10的流场总温分布等值线及云图，观察图4（a）和图4（b）可以看出非常明显的弓形激波，钝头体迎风面前部流场温度等值线及云图分布基本一致。对比图4（a）和图4（b）可以发现，图4（b）钝头体背风区域温度稍微高于图4（a）钝头体背风区域温度，这可能是因为在此计算高度下，虽然计算状态10采用的壁面条件降低了初始设定的壁面温度，但是随着流场的发展，钝头体背风面的热流密度进一步增加，采用本文发展的壁面条件能够依据热流密度适时更新壁面温度，使背风面温度高于初始设定的恒定温度，进而导致流场背风区域温度高于计算状态9。
图5给出了计算状态1～2和计算状态5～10的后体驻点线总温分布。观察图5（a）可以发现相比于计算状态1，计算状态2采用的壁面条件使流场背风区总温降低明显，观察图s（b）可以发现相比于计算状态5，计算状态6采用的壁面条件使流场背风区总温降低，但是降低幅度明显低于图5（a），观察图5（c）可以发现相比于计算状态7，计算状态8采用的壁面条件反而使流场背风区总温有所提升，观察图5（d）可以发现相比于计算状态9，计算状态10采用的壁面条件使流场背风区总温有所提升，与图5（c）结果类似。
表2给出了不同计算状态后体驻点线在X=1.5 m和X=3 m处的总温以及两坐标处的总温之差。观察表2可以发现，计算状态1和计算状态2的△T相差最大，表明在此计算高度下本文采用的气体壁面相互作用模型对温度变化的敏感度最高，进一步观察计算状态5～10可以发现，两种模型△T有接近的趋势，这表明随着计算高度上升本文采用的气体壁面相互作用模型对温度变化的敏感度有减弱的趋势。
图6给出了计算状态1和计算状态2的驻点线气动参数分布，观察图6（a～d）可以发现两种计算状态下的驻点线密度、平动温度、转动温度、速度分布基本保持一致，但计算状态2采用的气体壁面相互作用模型使激波离驻点距离略微缩短。
3 结束语
本文开展了基于二维非结构网格的高超声速稀薄流DSMC气体与壁面相互作用数值模型的研究，采用了基于辐射平衡的壁面温度条件，该边界条件能够克服恒温壁面条件的自身缺陷，在DSMC模拟过程中适时给出更贴近真实情况的壁面温度。并且，采用由完全漫反射和纯镜面反射模型组合而成的壁面反射模型，能够更好地描述模拟分子在壁面散射的真实情况。数值算例表明本文采用的气体壁面相互作用模型能够提高钝头体背风面后部流场对温度变化的敏感度，并且随着计算高度的提升，钝头体背风面后部流场对温度变化的敏感度有减弱的趋势。而在钝头体迎风面，驻点线气动参数与采用恒温壁面条件和完全漫反射模型状态的结果基本保持一致。
参考文献：
[1] 沈青.稀薄气体动力学[M].北京：国防工业出版社，2003.
Shen Qing. Rarefied gas dynamics[M]. Beijing： National Defence Industry Press， 2003.
[2] 樊菁.稀薄气体动力学：进展与应用[J].力学进展，2013，43（2）：185-201.
Fan Jing. Rarefied gas dynamics： Advances and applications[J]. Advances in Mechanics， 2013，43 （2）：185-201.
[3] 叶品.高超声速稀薄流中化学反应的DSMC数值模拟研究[D].西安：西北工业大学，2007.
Ye Pin. DSMC method for the simulation of chemical reaction of hypersonic rarefied gas flow[D]. Xi'an： Northwestern Polytechnical University， 2007.
[4] Celenligil M C， Moss J N. Application of the DSMC method to hypersonic flow about a delta wing[R]. AIAA-95-2053，1995.
[5] 王学德.高超声速稀薄气流非结构网格DSMC及并行算法研究[D].南京：南京航空航天大学，2006.
Wang Xuede. DSMC method on unstructured grids for hypersonic rarefied gas flow and its parallelization [D]. Nanjing： Nanjing University of Aeronautics and Astronautics，2006.
[6] 蔡国飙，郭景晶，庄逢甘，等，考虑化学反应的真空羽流的数值模拟[J].宇航学报，2001，22（4）：13-19.
Cai Guobiao， Guo Jingjing， Zhuang Fenggan， et al. Numerical simulation ofvacuum plume with chemical reaction[J]. Journal of Astronautics， 2001， 22（4）：13-19.
[7] 焦子龙，庞贺伟，海洋卫星水色仪羽流污染DSMC模拟研究[J].宇航学报，2007，28（6）：1529-1532.
Jiao Zilong， Pang Hewei. DSMC， simulation on plume contamination to ocean colorscanner of HY-1 satellite[J]. Journal of Astronautics， 2007，28（6）：1529-1532. [8] 贺群武，王秋旺，索晓娜.DSMC方法下给定热流边界条件的实施[J].工程热物理学报，2006， 27（4）：684-688.
He Qunwu， Wang Qiuwang， Sou Xiaona. Implementation of heat flux specified boundary conditions in DSMC， method[J]. Journal of Engineering Thermophysics， 2006，27（4）：684-688.
[9] 黄飞，赵波，程晓丽，等，热力学非平衡状态下DSMC绝热壁边界的抽样方法[J].宇航学报，2013，34 （11）：145l-1455.
Huang Fei， Zhao Bo， Cheng Xiaoli， et al. A new sampling method of Adiabatic boundary condition in DSMC under thermodynamic nin-equilibrium [J]. Journal of Astronautics， 2013，34 （11）： 1451-1455.
[10] Padilla J F. Assessment of gas-surface interaction models for computation of rarefied hypersonic flows [D]. Michigan： The University of Michigan，2008.
[11] Gallis M A， Torczynski J R， Rader D J， et al. Convergence behavior of a new DSMC algorithm [J]. Journal of Computational Physics， 2009（228）： 4532-4548.
[12] 周欣欣，吴颂平.高超声速非定常流动的数值模拟与气动热计算[J].计算力学学报，2008，25（S）：23-28.
Zhou Xinxin， Wu Songping. Numerical simulation of unsteady hypersonic flow and acrodynamic heating calculation [J]. Chinese Journal of Computational Mechanics， 2008，25（S）：23-28.
[13] Boyd I D. Computation of hypersonic flows using the direct simulation monte carlo method [C]//21st AIAA Computational Fluid Dynamics Conference. San Diego， CA：[s.n.]，2013.
[14] Bird G A. Molecular gas dynamics and the direct simulation of gas flow [M]. Oxford： Clarendon Press， 1994.