黄 飞，陈 智，程晓丽，沈 清

（中国航天空气动力技术研究院，北京 100074）

0 引 言

DSMC方法自出现以来已成功应用于航天、微机电、材料加工等领域［1－5］。在DSMC模拟方法中，网格的主要作用有两个方面，一是为了采样统计出宏观流场变量，另一方面主要是为了在模拟过程中选取碰撞对，即两个相互碰撞的分子必须在同一网格中。研究表明，DSMC模拟中必须保证网格尺寸在当地分子自由程的三分之一量级才能有效保证模拟结果的正确性，即在DSMC模拟过程中碰撞对的选取在空间上必须满足一定的限制，这一要求给DSMC模拟相对较高的密度区域（自由程相对较小）带来了很大难度，并且在全流场模拟过程中很难保证全流场网格尺寸满足这一要求。

为了降低模拟过程中对网格尺寸的限制，同时在不损失精度的前提下提高计算效率，Bird［6］引入了子网格的概念，该方法将流场中的计算网格在不同方向分为若干子网格，流场中的碰撞对在子网格中选取，而背景网格用来进行流场的采样统计，这一方法的成功应用被许多DSMC研究者所采纳。然而，这一静态子网格的方法很难保证子网格中的粒子至少在两个以上，尤其是在对高密度区域进行较多子网格划分时，为了解决大量静态子网格中可能少于两个模拟粒子的问题（尤其在三维情况下），Bird［7］引入了动态子网格方法，该方法根据背景网格中的平均分子数动态分配子网格数，这种方法大大改善了静态子网格的缺陷，但仍然难以消除此类问题出现的可能性；另一方面，无论是静态子网格还是动态子网格都一定程度上增加了碰撞对选取时此类特殊情况下的逻辑判断。鉴于此，LeBeau［8］等人引入了虚拟子网格方法，该方法从网格中首先选取碰撞对中的第一个分子，然后遍历其它所有分子，选择离第一个分子最近的分子作为碰撞对中的另一个分子，此种选取方式有效避免了Bird子网格方法中的缺陷，然而该方法需要对网格中的分子数N进行有效控制，不能太大，否则其在网格中进行的遍历搜索，引入了N2次方量级的搜索次数及距离计算，大大增加了计算量。国内关于DSMC模拟技术方面的工作主要集中于IP方法［9］、并行算法［10］、位置元方法［11］、网格处理［12］及 LD－DSMC［13］等方面，尚未发现关于权因子方面的研究报道。

综上分析可以发现，无论是子网格还是虚拟子网格，其本质都是尽可能的降低选取碰撞对时两个分子之间的空间距离，尽可能保证相对较近的分子得到有效碰撞。基于此，本文在虚拟子网格的方法上引入自适应的碰撞距离思想，将相对较粗的网格中的碰撞对限制在合理的距离范围内以满足近距离的碰撞要求，与LeBeau虚拟子网格方法不同的是，本文在选取碰撞对时不需要遍历网格中所有的分子，因而也无需对网格中的模拟分子数进行多余的控制，从而能够在满足计算精度的要求下放宽对网格的限制。

1 计算方法

1.1 DSMC 方法介绍［6］

DSMC方法的主要思想是采用有限个模拟分子代替真实分子，在计算机中进行分子与分子之间、分子与物面之间的相互碰撞模拟，整个模拟中分子之间及其与物面之间不断交换动能与内能，待模拟足够的时间步以后，采样统计给出每个网格中的宏观流场结果。计算中分子之间采用VHS碰撞模型，能量交换采用L－B模型［14］，按照Bird的能量按自由度分配原则采用取舍法进行抽样分配物面边界条件采用完全漫反射条件。

1.2 碰撞对选取的主要关键技术点

碰撞距离的求解方法：

由于本文计算代码为非结构网格，故而主要以非结构四面体网格为例进行介绍。从前面分析可以发现，本文的方法关键在于自适应碰撞距离d的确定，d的合理选取直接关系到模拟结果的准确性。计算求解中d由下式确定：

其中，L＝Lmax／（／n），λi为网格中的局部自由程，由流场参数确定，Lmax为四面体中的最大边长，n为每个等效虚拟子网格中的最小分子数，为可调参数，为采样一定步数后网格中平均分子数，即为：

从以上可以发现，碰撞距离d由两个参数L、λi确定，λi／3的限制保证了碰撞距离限制在三分之一局部自由程范围之内，并且λi随流场宏观参数的变化而改变，本文在DSMC采样前流场每隔500步更新一次，同时对相应的自由程λi、碰撞距离d都进行更新，从而能够有效做到碰撞距离随流场的变化进行自适应更新；L主要是为了解决在模拟过程中局部高密度网格中模拟分子数相对较少时碰撞对难以选取的问题。

DSMC模拟中碰撞对选取的主要步骤：

1）给定碰撞对中第二个分子选取的循环次数限制数m及第一个分子重新选取的循环次数限制数p。

2）首先在网格中随机选取碰撞对中的第一个分子A。

3）在网格中随机选取第二个分子B，计算两个分子之间的距离di。

今天，我们主要来了解一下中国古代最具代表性的十大乐器，分别是 ：琵琶、二胡、编钟、箫、笛、瑟、琴、埙（）、笙、鼓。

4）判断0＜di≤d是否满足，满足则停止搜索进行分子之间的碰撞；不满足接着判断循环次mi数是否满足mi＜m，满足返回至3），不满足则继续判断pi是否满足pi＜p，满足则返回至2），不满足则停止搜索，进行分子之间的碰撞。

1.3 气动力参数定义

2 算例验证一：平板外形

2.1 计算模型及计算条件

本文采用不同网格尺寸下的平板外形对计算方法进行了系统的计算分析。其中计算模型及尺寸如图1所示，板长L＝1.2043m，图中边界面A′ABB′、C′CD′D及 A′C′FE为对称边界，ABDC及 ACC′A′为远场来流条件，EFD′B′为壁面边界，BDD′B′为真空出口边界；来流气体为氮气，速度1412.5m／s，温度300K，密度4.65×10－6，壁温500K。

图1 计算模型示意图Fig.1 Model in simulation

2.2 结果分析

为了对本文所发展的虚拟子网格方法进行验证计算，本文在三维计算域中共选取了三种不同尺寸的网格进行了计算分析，其三组网格数分别为10850、74897及334275个单元，计算中的全场模拟分子数约分别为25万、120万及350万。其中第一组网格为最粗的网格，对其分别采用虚拟子网格方法与不采用该方法两种工况进行了计算分析。

图2为计算给出的流场压力等值线云图，从中可以看出超声速来流在平板前形成了斜激波，激波边界层的相互干扰使得波后压力升高，激波后在靠近前缘点附近压力升至最高，之后沿壁面逐渐降低，在出口处受真空边界的影响，压力减小变快。图3为不同网格下的物面压力系数分布结果。从中可以看出，物面压力系数的整体分布规律与前面所述规律类似，在前缘附近先升至最高，之后沿板面逐渐降低，在板尾部降低最快。从中还可发现，相对较密的后两组网格所得结果与Bird结果吻合较好，而10850个计算单元下的粗网格在不采用虚拟子网格方法时的计算结果明显偏离Bird所得计算结果。当加入子网格方法后此粗网格亦能得到与Bird一致的计算结果。此外，在板前缘压力系数结果吻合较差，这主要是由于在板前缘流场压力梯度相对较大，相对较粗的网格统计结果难以描述这种梯度变化。

图2 流场压力等值线云图（334275cell）Fig.2 Pressure contours

图3 不同网格下的物面压力系数分布Fig.3 Pressure coefficient distribution along the wall

图4、图5分别为不同网格下的物面摩擦系数及物面热流系数的分布结果。从中可以看出，在前缘处摩阻系数、热流系数达最大值，随后沿壁面降低，在尾缘处由于气体真空膨胀加速，使摩擦系数、热流系数在尾缘处有所增加。从中亦可发现，相对较密的后两组74897、334275个单元的网格计算结果都与Bird计算结果吻合很好，与压力系数的吻合规律类似，10850个计算单元下的粗网格的计算结果在采用本文所发展的虚拟子网格方法之后亦能得到满意的结果。

图4 不同网格下的物面摩擦系数分布Fig.4 Friction coefficient on the surface

图5 不同网格下物面热流系数分布Fig.5 Heat flux coefficient on the surface

3 算例验证二：钝锥外形

3.1 计算模型及计算条件

本文算例二主要针对钝锥外形进行了计算分析，其中计算网格共采用了三种不同的尺寸。计算模型尺寸如图6所示，钝锥头部半径为R＝35mm，半锥角5°，长1m。来流气体为空气，来流速度为5000m／s，飞行高度90km，壁温300K。

图6 计算模型示意图Fig.6 Model in simulation

3.2 结果分析

计算分析中采用的三种不同尺寸的网格数分别为654114、128454及34173个单元，其中最后一组最粗的网格采用本文所发展的虚拟子网格与不采用该方法两种工况进行了计算分析。图7、图8分别给出了流场压力云图及不同网格尺寸下的表面压力分布结果，从中可以看出在钝锥头部形成了明显的激波结构，在驻点区域压力较高，随后沿锥面降低，在驻点下游附近由于膨胀作用，压力有所下降，而在锥尾部区域由于出口边界为真空边界条件，流动在此区域加速膨胀使得压力在尾部有所下翘；从中还可发现，在粗网格下采用本文所发展的虚拟子网格方法能够得到与密网格一致的结果，而不采用该方法时，较粗网格下的结果明显偏离其它密网格的结果。

图7 流场压力云图（654114cell）Fig.7 Pressure contours

图8 不同网格下的物面压力分布Fig.8 Pressure on the surface at different grid sizes

图9、图10分别给出了不同网格尺寸下表面热流及摩阻分布，从中可以看出，相对较密的前两组网格654114、128454个单元的网格结果一致，而采用虚拟子网格方法34173个单元的粗网格亦能得到与密网格一致的热流、摩阻分布结果，同压力类似，不采用虚拟子网格方法的粗网格单元很难与密网格结果吻合。通过钝锥外形的算例分析进一步证明了本文所发展的方法能够在相对较粗的网格下取得相对较为满意的结果。

图9 不同网格下的物面热流分布Fig.9 Heat flux on the surface at different grid sizes

图10 不同网格下的摩阻分布Fig.10 Friction on the surface at different grid sizes

4 结 论

本文针对DSMC模拟过程中的碰撞特征，经过计算分析给出了一种可进一步提高DSMC方法模拟效率的自适应虚拟子网格方法。通过在粗网格下的计算分析可以看出，本文所发展的方法能够在相对较大的粗网格下得到与密网格一致的结果。本文的方法一定程度上放宽了DSMC模拟中的网格尺度，有效降低了全流场的模拟分子数，从而节省了内存，提高了计算效率。

［1］SHEN Q.DSMC method and the calculation of rarefied gas flow［J］.AdvanceinMechanics，1996，26（1）：1－13.（in Chinese）沈青.DSMC方法与稀薄气体计算的发展［J］.力学进展，1996，26（1）：1－13.

［2］WU Q F，CHEN W F.DSMC method of nonequilibrium thermodynamics and chemical reaction in the high temperature rarefied gas flow［M］.Changsha：National University of Defense Technology Press，1999.（in Chinese）吴其芬，陈伟芳.高温稀薄气体热化学非平衡流动的DSMC方法［M］.长沙：国防科技大学出版社，1999.

［3］RAULT D F G.Study of shuttle body flap effectiveness at high altitudes using DSMC［R］.AIAA－94－2021，1994

［4］HUANG F，ZHAO B，CHENG X L，et al.Real gas effects of reentry vehicle aerodynamics under hypersonic［J］.Journalof Astronautics，2014，35（3）：283－290.（in Chinese）黄飞，赵波，程晓丽，等.返回器高空稀薄气动特性的真实气体效应研究［J］.宇航学报，2014，35（3）：283－290.

［5］LVANOV M，GIMELSHEIN S.Influence of real gas effects on control surfaces efficiency at high flight altitudes［R］.AIAA 93－5116.

［6］BIRD G A.Molecular gas dynamics and direct simulation of gas flow［M］.London：OxfordUniversity Press，1994，438－451.

［7］BIRD G A.Forty years of DSMC，and Now？［A］.22ndInternational Symposium of Rarefied Gas Dynamics［C］.Sydney Australia，2000.

［8］LEBEAU G J，BOYLES K A，LUMPKIN F E.Virtual subcells for the direct simulation Monte Carlo method［R］.AIAA 2003－1031.

［9］FAN J，SHEN C.Micro－Scale gas flows［J］.AdvanceinMechanics，2002，32（3）：321－336.（in Chinese）樊菁，沈青.微尺度气体流动［J］.力学进展，2002，32（3）：321－336.

［10］WANG X D，WU Y Z，XIA J，et al.A paprallel algorithm of 3Dunstructured DSMC method and its application［J］.Journal ofAstronautics，2007，28（6）：1500－1505.（in Chinese）王学德，伍贻兆，夏健，等.三维非结构网格DSMC并行算法及应用研究［J］.宇航学报，2007，28（6）：1500－1505.

［11］FAN J，LIU H L，SHEN C，et al.A molecular reflection deterministic criterion used in the position element algorithm of direct statistical simulation［J］.ACTAAerodynamicaSinica，2000，18（2）：180－187.（in Chinese）樊菁，刘宏立，沈青等.直接统计模拟位置元方法中的分子表面反射确定论判据［J］.空气动力学学报，2000，18（2）：180－187.

［12］LIANG J，YAN C，DU B Q.An algorithm study of three－dimensional DSMC simulation based on two－level cartesian coordinates grid structure［J］.ACTAAerodynamicaSinica，2010，28（4）：466－471.（in Chinese）梁杰，阎超，杜波强.基于两级直角网格结构的三维DSMC算法研究［J］.空气动力学学报，2010，28（4）：466－471.

［13］SU W，HE X Y，CAI G B.Extension of the low diffusion particle method for near－continuum two－phase flow simulations［J］.ChineseJournalofAeronautics，2013，26（1）：37－46.

［14］BORGANOFF C，LARSEN P S.Statistical collision model for Monte Carlo simulation of polyatomic gas mixture［J］.Journal ofComputationalPhysics，1975，18：405－420.