高超声速混合流动的粒子模拟耦合算法*

姜婷婷，钱 坤，陈伟芳

(浙江大学 航空航天学院，杭州 310027)

高超声速混合流动的粒子模拟耦合算法*

姜婷婷，钱 坤，陈伟芳

(浙江大学 航空航天学院，杭州 310027)

通过采用基于三维非结构网格的子网格技术，对适用于高超声速混合流动的自适应时间步长粒子模拟耦合算法(Improved Hybrid Particle Simulation Method, IHPSM)进行了改进，在保证算法计算效率的同时降低数值误差。通过对三维双曲钝锥外形的数值仿真及与DSMC(Direct Simulation Monte Carlo)计算结果的对比分析，证明改进的IHPSM算法具有较高的计算精度，且能较大幅度提高计算效率。基于改进的IHPSM算法，文中针对双曲钝锥外形进行了稀薄气体效应和飞行马赫数对高超声速流动影响规律的研究。结果表明，气体稀薄程度的增加会减缓双曲钝锥前端流场宏观物理量的变化梯度，使流场激波结构变厚；来流马赫数的增大会使激波明显增强，但对激波厚度与结构的影响较小。

子网格技术；高超声速混合流动；粒子模拟耦合算法；IHPSM

0 引言

钱学森[1]根据Kn数的数值大小，将流动按照气体的稀薄程度划分为连续流区、滑移流区、过渡流区与自由分子流区。在连续流区，通常通过数值求解N-S(Navier-Stokes)方程研究飞行器的气动特性；在气体稀薄程度很大的区域，N-S方程失效，通常通过数值求解简化的Boltzmann方程或者DSMC(Direct Simulation Monte Carlo)方法来研究飞行器的气动特性。然而，高超声速飞行器在近连续流区以及过渡流区飞行时，高速自由来流与飞行器控制面或翼面上的激波、边界层相互作用，会出现连续流和稀薄流共存的现象[2-4]，即混合流动现象。在混合流动中，局部稀薄流的存在使得基于连续介质模型假设的经典连续性方程在理论上失效，而局部高密度区的存在，使得DSMC方法因其巨大的计算资源耗费量亦失去应用价值。近年来，N-S/DSMC耦合算法在求解这类流动问题中受到学术界广泛关注[5-6]。但由于N-S方程与DSMC方法模拟本质的巨大差异，N-S/DSMC耦合算法在分区并行计算、界面信息传递、复杂外形适应性等方面仍面临巨大挑战。

不同于N-S/DSMC耦合算法的思想，Macrossan[7]提出了一种粒子模拟耦合算法(Hybrid Particle Simulation Method, HPSM)。该方法采用DSMC方法计算流场的局部稀薄流区，采用基于DSMC方法发展而来的平衡粒子模拟方法(Equilibrium Particle Simulation Method, EPSM)[8]计算流场的局部连续流区。与DSMC方法类似，EPSM对仿真分子进行数值模拟。但不同的是，EPSM省去了复杂的碰撞过程，而是在遵守能量守恒和动量守恒的基础上，通过从Maxwell平衡分布中随机抽样来得到仿真分子碰撞之后的状态，因此提高了计算效率。基于Macrossan的思想，陈伟芳[9-10]、朱荣丽[11]等研究者进一步在计算模型、计算效率等方面发展了粒子模拟耦合算法，并通过相关高超声速典型算例验证了算法的有效性。由于在整个计算区域内均采用粒子模拟方法进行数值求解，HPSM具有物理模型清晰、区域边界条件自然、界面信息传递简单等优点。然而，尽管如此，目前国内外的研究者针对HPSM所做的研究仍不充分，若要实现其广泛的应用，仍需解决诸多难题。首先，传统的连续介质假设失效判据在HPSM中并不适用，因此亟需建立一套合理的划分DSMC计算区域和EPSM计算区域的流场区域划分判据。其次，虽然HPSM已经在一定程度上较DSMC方法提高了计算效率，但距离工程应用的需求仍有较大差距。此外，粒子模拟方法(无论是DSMC方法还是EPSM)对计算网格的尺寸较为敏感[12-13]，当计算网格尺寸远大于当地平均分子自由程时(这种情况在高密度的局部连续流区不可避免)，会导致该方法产生较大的数值误差。在前期研究工作[14]中，本文作者针对前两个问题，构建了适用于HPSM算法的平衡态失效判据(Equilibrium Breakdown Parameter)和自适应时间步长方法，从而建立了改进的HPSM算法(Improved Hybrid Particle Simulation Method, IHPSM)。验证算例表明，IHPSM在保证计算精度的同时较大程度的提高了计算效率。

本文针对第3个问题，在IHPSM算法的基础上引入三维非结构子网格技术，以减小算法对网格尺寸的依赖性。同时，本文通过对超声速双曲钝锥外形的数值计算来验证引入子网格技术的IHPSM的可行性。此外，本文采用IHPSM算法开展了稀薄气体效应和飞行速度对高超声速双曲钝锥外形流场结构的影响规律研究。

1 自适应时间步长IHPSM算法

IHPSM[14]通过平衡态失效判据划分计算区域，然后在满足连续介质模型假设且达到局部平衡态分布的流动区域使用EPSM进行计算，其他区域使用DSMC方法进行计算[15]，从而实现对整个流场的数值求解。

1.1 平衡态失效判据

在N-S/DSMC耦合算法中，通常采用连续介质模型假设失效判据来划分N-S方程与DSMC计算区域。而EPSM适用于气体状态达到近平衡态分布时的情况，它相当于通过数值模拟仿真分子的运动与碰撞过程来求解Euler方程。因此，EPSM的适用条件更为严苛。在IHPSM中，首先通过连续介质假设失效判据来划分局部连续流区和局部稀薄流区，然后在局部连续流区通过“平衡态失效参数[14]”判断网格内分子状态是否处于局部平衡态分布，即是否满足EPSM的使用条件。

在此，采用的连续介质假设失效判据为当地克努森数[2]，其定义为

(1)

式中λ为平均分子自由程；ρ为气体密度；T为平动温度；u为速度；a为声速。

当网格内KnGLL,max>0.05时，为稀薄流动区域，采用DSMC方法来模拟该区域的碰撞过程；当KnGLL,max<0.05时，为连续流动区域，此时需要继续判断该网格是否达到平衡分布状态。基于Bird[16-17]的研究，选取网格中的分子碰撞数作为判断网格内的流动是否达到平衡态的准则。分子碰撞数f的定义：

(2)

式中 Δt为计算时间步长；γ为分子碰撞频率。

对于VHS模型来说，分子碰撞频率γ的表达式：

式中dref为气体分子在参考温度为0 ℃时的参考直径；n为气体分子数密度；kB为Boltzmann常数；ms为气体分子质量；Tref为参考温度；ω为粘性指数。

在一个时间步长内，若网格中的平均分子碰撞数f≥10，则网格内的分子运动状态将达到平衡态，即可采用EPSM进行计算。

1.2 自适应时间步长方法

为了提高计算效率，在IHPSM中，每个网格的时间步长不再保持统一，而是在计算过程中不断调整，以适应流场物理量的变化。

在自适应时间步长方法[14]中，定义3种与时间步长相关的变量：基时间步长Δt、时间区Ti(Ti≥1，且为整数，i为网格单元编号)以及网格时间步长Δti。基于DSMC方法对时间步长的要求，选取初始基时间步长[11]Δt：

(3)

式中i为网格编号；Δxi、Δyi、Δzi为网格i在3个方向的参考长度；u∞、v∞、w∞为来流气体运动速度；γ∞为来流气体碰撞频率。

Δt将作为整个流场推进过程中的最大时间步长，且不随着推进过程发生改变。在高超声速流场中，大多数网格的分子碰撞频率随着时间推进过程而增加。因此，这个基时间步长对于大多数网格来讲不会出现因时间步长过小而降低计算效率的情况。

确定了基时间步长以后，网格内当地时间步长随计算过程的自动调整步骤可描述如下：

(1)将初始时间区设为Ti=1；

(2)每隔一个基时间步长Δt，计算网格 的分子碰撞数fi，同时将时间区更新为Ti=int(fi)+1；

(3)令

(4)

1.3 子网格计算技术

1.4 IHPSM算法的主要计算步骤

IHPSM算法的主要计算步骤描述如下：

(1)按照式(3)求出初始基时间步长Δt，并赋予各网格初始时间区Ti=1，则各网格的初始时间步长均为Δt。

(2)在一个时间步长内，求出每个网格中的当地克努森数KnGLL，max和碰撞数fi，同时定义碰撞数累计值CNi=CNi+fi。

(3)当KnGLL，max≥0.05时：若fi≤1，则令网格属性IDi=1；若fi>1，则令IDi=-1。

当KnGLL，max<0.05时：若fi≤1，则令IDi=1；若1

为了更加清晰地描述该步骤，图1给出了较详细的计算流程。

(4)每个基时间步长之后，对仿真分子进行重新排序，同时统计仿真分子的物理信息，以便得到宏观物理信息。

(5)重复步骤(2)～(4)，直到计算结束。

2 数值模拟结果及分析

2.1 高马赫数双曲钝锥外形验证算例

为了验证自适应时间步长IHPSM计算程序，对高超声速过渡流条件下的三维双曲钝锥缩比模型进行数值计算，并与DSMC计算结果进行对比。其中，DSMC计算结果由课题组内部已经经过验证的DSMC计算软件[18]计算得到。三维双曲钝锥外形双曲面半角为42.5°，头部半径为0.013 62 m。来流介质为空气，计算条件：

图3为IHPSM方法计算得到的流场马赫数云图、压力云图和温度云图。从流场马赫数云图可看出，气流经过激波后开始减速，从超声速流动减为亚声速流动。流场压力云图表明，激波后压力增大，在驻点处压力达到最大值。激波后的温度变化则是先增大而后在驻点处又稍微降低，这主要是由等温壁条件所致。

图4给出了IHPSM方法、DSMC方法驻点线温度、密度的分布曲线，通过对比可看出，IHPSM方法获得了与DSMC计算方法基本一致的计算结果，表明IHPSM方法可较精确地描述高超声速混合流动的流场结构。

IHPSM方法和DSMC方法的测试计算机配置均为2 GHz的AMD处理器，内存24 G。2种方法均采用单核进行计算。流场达到稳定状态时，DSMC方法所耗费的计算时间约为720 h，而本文改进的IHPSM方法花费的计算时间仅为310 h。因此，在相同的仿真分子运动时间内，IHPSM方法在取得较高预测精度的前提下，显著地提高了计算效率。

2.2 双曲钝锥外形流场结构的变化规律研究

本节采用IHPSM算法开展气体稀薄效应与飞行速度对三维双曲钝锥外形流场结构的影响规律的研究。

首先，在来流马赫数为20、物面温度为500 K的条件下，取飞行高度分别为75、80、85、90 km的计算状态进行数值仿真，并对流场结果进行对比分析。不同飞行高度下计算得到的驻点线密度和温度分布曲线如图5所示。

从图5中各物理量变化趋势可看出，随着飞行高度的增大，气体越来越稀薄，钝锥前端流场中密度、温度的变化梯度逐渐减小，流场激波结构越来越厚。这是因为，随着来流气体稀薄程度的增大，分子的碰撞频率减小，粘性作用也因此而减小，从而导致流场物理量变化率减缓，宏观表现便是激波的厚度增加。

此外，在飞行高度80 km、物面温度为500 K的条件下，通过改变来流速度来研究来流马赫数变化对流动的影响。来流马赫数分别取6、10、15、20。不同马赫数状态下计算得到的驻点线密度和温度分布曲线如图6所示。从图6中可看出，随着来流马赫数增大，流体的动能在头部驻点区转化为气体内能，激波明显增强，导致激波后驻点区密度越来越大，而温度急剧升高。但从密度与温度分布来看，来流马赫数对激波厚度的影响不大。

3 结论

(1)在IHPSM算法中加入三维非结构子网格技术，降低了算法由于对计算网格的依赖性而产生的数值误差。

(2)通过对三维双曲钝锥外形的数值仿真以及与DSMC计算结果的比对分析，表明改进的IHPSM方法能同时满足高超声速混合流动问题对计算精度和计算效率的需求。

(3)针对双曲钝锥外形,进行了稀薄气体效应和来流马赫数对高超声速流动的影响规律研究。结果表明,随着来流气体稀薄程度的增加，钝锥前端流场宏观物理量的变化梯度会减缓，流场激波结构越来越厚；随着来流速度的增大，激波会明显增强，但激波的厚度与结构的变化不大。

[1] Tsien H S.Superaerodynamics,mechanics of rarefied gases[J].Journal of the Aeronautical Sciences,1946,13(12):653-664.

[2] Boyd I D,Chen G,Candler G V.Predicting failure of the continuum fluid equations in transitional hypersonic flows[J].Physics of Fluids,1995,7(1):210-219.

[3] Gnoffo P A.CFD validation studies for hypersonic flow prediction [R].AIAA 2001-1025.

[4] Hash D B,Hassan H A.Assessment of schemes for coupling monte carlo and navier-stokes solution methods[J].Journal of thermophysics and Heat Transfer,1996,10(2):242-249.

[5] Schwartzentruber T E,Boyd I D.A modular particle-continuum numerical method for hypersonic nonequilibrium gas flows[J].Journal of Computational Physics,2007,225:1159-1174.

[6] Lian Y Y,Wu J S,Cheng G,et al.Development of a parallel hybrid method for the DSMC and NS solver[R].AIAA 2005-0435.

[7] Macrossan M N.A particle-only hybrid method for near-continuum flows[C]//American Institute of Physics,2001,388-395.

[8] Pullin D I.Direct simulation methods for compressible ideal gas flow[J].Journal of Computational Physics,1980,34:231-244.

[9] 陈伟芳,吴明巧,任兵.DSMC/EPSM混合算法研究[J].计算力学学报,2003,20(3):274-278.

[10] 吴明巧,陈伟芳,任兵.超声速平头圆柱绕流DSMC/EPSM混合算法研究[J].空气动力学学报,2002,20(2):206- 210.

[11] 朱荣丽,曹义华.高超声速飞行器DSMC/EPSM自适应当地时间步长混合算法 [J].空气动力学学报,2008,26(3):394-399.

[12] Macrossan M N.The equilibrium flux method for the calculation of flows with non-equilibrium chemical reactions[J].Journal of Computational Physics,1989,80(1):204-231.

[13] Breuer K S,Piekos E S,Gonzales D A.DSMC simulations of continuum flows[R].AIAA 1995-2088.

[14] Jiang T T,Xia C C,Chen W F.An improved hybrid particle scheme for hypersonic rarefied-continuum flow[J].Vacuum,2016,124:76-84.

[15] Laux M.Local time stepping with automatic adaptation for the DSMC method[R].AIAA 1998-2670.

[16] Bird G A.Near continuum impact of an under-expanded jet[C]//Proc.AIAA Computational Dynamics Conference,1973:103.

[17] Bird G A.Molecular gas dynamics and the direct simulation of gas flow[M].Clarendon:Oxford,1994.

[18] 吴雄,陈伟芳,石于中.过渡区球头气动力特性DSMC仿真研究[J].航天返回与遥感,2002,23(3):1-5.

(编辑：吕耀辉)

Study on hybrid particle simulation method for hypersonic mixed flows

JIANG Ting-ting,QIAN Kun,CHEN Wei-fang

(School of Aeronautics and Astronautics,Zhejiang University,Hangzhou 310027,China)

A hybrid all-particle scheme(Improved Hybrid Particle Simulation Method,IHPSM),which is intended for hypersonic mixed rarefied-continuum flows,is extended for improved efficiency and wider applicability by employing the three-dimensional unstructured sub-cell technique.To assess the validity of the improved hybrid particle simulation scheme,a typical test case for three dimensional hyperbolic blunt cone is simulated.The results of IHPSM method for flow field are in good agreement with the standard DSMC results.Notably,the hybrid procedures allow a large reduction in overall simulation expense relative to DSMC.In addition,the relationship between the hypersonic flow field and the rarefied gas or Mach number is studied by simulating a series of test cases of the hyperbolic blunt cone.As the inflow gas becomes rarer,the gradient of the macroscopic properties in front of the hyperbolic blunt cone becomes smaller,and the shock wave is thicker.The increase of Mach number increases the values of the macroscopic properties in the shock wave,but has little influence on the shock wave structure.

sub-cell technique;hypersonic mixed flows;hybrid all-particle scheme;IHPSM

2016-03-09；

2016-03-21。

国家重点基础研究发展计划(2014CB340201)；国家自然科学基金(51575487)。

姜婷婷(1989—)，女，博士生，流体力学专业。E-mail：jiangyan8140@163.com

V411

A

1006-2793(2017)03-0283-06

10.7673/j.issn.1006-2793.2017.03.003