航空动力学报
JOURNAL OF AEROSPACE POWER
1999年　第14卷　第2期　Vol.14　No.2　1999



真空羽流场的DSMC并行数值模拟*
蔡国飙**　刘世俭　王慧玉　庄逢甘
　　【摘要】　对DSMC方法进行小喷管真空羽流数值模拟的计算软件进行了并行化改造。计算表明在多处理器（节点）上实现并行计算能够大幅度地提高运算速度，缩短运行时间，增加计算规模，计算结果与串行计算一致。并对目前串行计算不能完成的真空羽流计算的算例成功地进行了DSMC并行计算，计算结果与理论分析一致。
　　主题词：　喷管气流　蒙特卡罗法　数值模拟　计算　　　　　　　自由词：　并行计算
　　分类号：　V439
1　引言
　　真空羽流场包括羽流核心区的连续流，外围区的自由分子流和过渡区的过渡流。其流动极为复杂，很难用一种方法来描述，特别是对于过渡区的流动，目前尚无完善的理论，与之相对应的数学方程是完全的Boltzman方程，除特殊情况外求解极为困难。比较成功地能够模拟真空羽流场的数值技术是直接模拟的蒙特卡洛方法（DSMC），这种方法由Bird G.A.等人在60年代提出，几十年来此方法被证明是成功的，并得到深入发展。但由于其计算量过大，计算时间过长，对内存的要求很高，而使得其应用受到限制［1］。
　　计算机技术的迅速发展，特别是网络技术在计算机领域的广泛应用，使得并行计算机在科学计算中的优势越来越明显。从理论上讲，并行计算规模可无限增大，并行计算速度可无限提高。真空羽流场DSMC实现并行计算，将能够满足其计算规模与计算速度的要求。
　　80年代末，Furlani T.R.和Lordi J.A.等人［2］，采用三种并行方法实现过渡流场DSMC并行计算，在32节点上计算速度最快可提高17倍。90年代初，Wilmmth R.G.和Carlsm A.B.等［3］在异构并行环境实现DSMC的并行计算。他们利用Sun工作站、Gray Ymp巨型机和Intel Ipsc/860系列机所构成的较为复杂的异构并行环境进行DSMC并行计算，得出与串行计算较为一致的结果，计算速度有更为明显的提高，计算规模能够基本上满足实际需要。
2　并行计算和PVM系统
　　并行计算机已有很长的发展历程，其种类较多，可分为单指令流多数据流(SIMD)和多指令流多数据流(MIMD)两大类。SIMD并行机按结果和运行机理可分为阵列机、流水线机和向量机。MIMD并行机按存储器的组成也可分为存储器分布式MIMD和存储器共享式MIMD并行机。
　　并行计算机的种类不同，相应的并行算法也截然不同，SIMD并行机需要对计算问题进行微观并行处理，对并行算法的依赖性强，串行占优的计算问题其并行效率很难提高。存储器分布式MIMD并行机可由识别串行算法的内在并行性加以改造，而使整个计算实现并行，对并行算法的依赖性较弱，通过适当选取并行对象能不同程度地包容串行计算过程。
　　本文应用的并行计算机是基于信息传递机制的松散耦合的大规模并行计算机系统(MPP)，属存储器分布式MIMD并行机。
　　存储器分布式MIMD并行机，能够同时对多条指令及其各自相关的数据进行处理。但其对系统软件要求比较高。在算法运行过程中，需要协调进程之间的同步问题以满足它们的数据往来或控制依赖要求。80年代末，美国橡树岭国家实验室（ORNL）推出的PVM系统能够较好地完成这些任务。PVM系统意即并行虚拟机（Parallel Virtual Machine），是一个并行系统软件包，它能够将异构并行网络和其它系列机用作同一个并行计算环境来操作，完成大规模的计算任务。
　　并行计算的目的是要克服运行速度和内存空间的不足，使计算顺利完成。因此实现并行化计算的原则应为：(1)尽可能提高计算速度；(2)能够处理更大规模的计算任务；(3)并行计算过程不影响计算结果的准确性。
　　并行算法与并行程序设计的优劣将直接导致并行计算效率有很大的差异，并行效果可以由并行加速比与可扩展性来衡量。并行加速比是表示采用多个处理器并行计算时计算速度所能得到的加速倍数。设tseq表示用串行计算机求解某个计算问题所需的时间，tp是用P个处理器求解问题的时间。在理论情况下，程序的每一部分若能完全并行，P个处理器的加速比应该等于P，但在一般情况下，这是不可能的。定义加速比Sp和并行效率Ep为［4］：
Sp＝tseq／tp　　　　　　　　Ep＝Sp／P
　　加速比Sp一般比P小，Ep也小于1，且随着处理器个数的增多Ep将减小。
3　DSMC并行计算及真空羽流场的DSMC并行计算过程
3．1　DSMC方法简介
　　DSMC方法是根据稀薄气体流动的数值模拟研究的需要发展起来的一种新的计算气体流动方法，以跟踪分子运动为依据，运用统计力学的方法来获得气体流动参数［1］。在流场中总有一些分子发生碰撞，而另一些分子自由运动，DSMC方法采取在一定时间间隔内将分子自由运动和分子碰撞解耦。模拟每个网格中的分子间的碰撞时，忽略分子的具体位置及其位移，而模拟分子自由运动和与壁面碰撞时忽略分子间的相互碰撞。利用少量的模拟分子代替真实流场内数目众多的气体分子，用计算模拟实际物理过程，即对模拟分子跟踪，记录它们的位置和速度变化，计算分子之间、分子和物面间的碰撞，再经统计平均获得所需的各种流动参数。
3．2　DSMC并行处理
　　由前所述，DSMC法进行真空羽流流场数值模拟可分为四部分：(1)分子的自由运动及与壁面碰撞；(2)分子由来流进入区域；(3)分子间的碰撞；(4)统计处理流场参数。主要计算量在前三部分，因此对第四部分计算不考虑并行处理。
　　分子自由运动和与壁面碰撞是各自独立无关的，因此对此过程并行处理可以分子为对象，最佳并行为各节点处理等量的分子数。
　　分子由来流进入区域与边界来流状态有关，边界上不同区域来流一般不同。对此过程并行处理应在考虑来流状态计算量的基础上，按空间区域并行，最佳并行为各节点处理等计算量的区域。
　　分子碰撞过程是DSMC法计算量最大的部分。由于位置不同其分子碰撞概率也不同，等网格数区域内的分子碰撞过程的计算量一般有较大的差异，因此并行处理时并行对象为区域，而各节点处理网格数的多少则应考虑碰撞概率的不同进行分配。
　　根据并行化原则的要求，各节点计算时在每个计算过程应保持同步以减少等待。同时由于信息通信量的大小严重影响运行速度，因此一个计算问题的整体并行需兼顾计算量的均匀分布和降低通讯量两个方面。DSMC真空羽流场计算并行化处理以计算量最大的分子碰撞过程为主要因素。根据碰撞概率的不同划分网格数，概率大的区域网格数少，概率小的区域网格数多，在此基础上按区域分子数众寡调整区域大小使并行计算接近同步。
3．3　二次过滤法
　　内存分布式MIMD并行计算机系统在进行并行计算时一般存在信息传递问题，在并行化处理时应尽可能减少信息传递量，恰当选取并行对象，充分利用各节点独立计算和信息处理过程中串行计算自依赖性特点，降低节点间信息依赖度。用DSMC法直接模拟流场计算时，由于分子运动是随机的，无论节点计算区域怎样划分，任何两个节点均存在信息传递问题。在节点较多时，信息传递方式则会严重影响并行加速比。


图1　二次过滤法
　　在存储器分布式MIMD并行机上，用PVM系统进行节点间通信是通过信息“打包”发送和接受“解包”实现的。针对DSMC并行计算的特点，本文提出了一种较为实用的通信方法—“二次过滤法”。如图1，每个计算节点处理的分子信息通过第一次过滤将不属于自己的分子信息过滤出来“打包”广播发送给其它节点，其它节点均接收到信息包后通过第二次过滤将属于本节点处理的信息留下。其特点是：每个节点只需要发送一个包就可将信息发送出去，同时每个节点可只考虑自己处理区域的边界情况而不必保存更多的信息，通过减少打包次数、降低信息通讯时间，当节点增多时可避免因频繁发送信息包所可能出现的梗塞现象，各节点不保存其它区域的边界信息，能够节约更多的内存。
3．4　真空羽流场的DSMC并行计算过程
　　DSMC法进行真空小喷管羽流场并行计算采用主－从模式完成，程序设计如图2所示。主控机处理主控程序，各节点完成从程序计算。
　　主控程序的作用在于启动从程序，将从程序加载到各节点上去，并随时进行监控管理，发出相应的指令，调整计算范围。对循环计算作出判断，接收各节点计算结果，并最终完成流场的参数统计计算。
　　从程序的作用是完成直接模拟过程，依据并行化原则将流场划分为多个计算区域，每个节点处理其中的一个区域。各节点的计算内容是相似的，均包括分子自由运动，分子逃逸注销，来流分子加入，分子碰撞等过程。在分子碰撞计算之前，节点间进行信息传递，将不属于本计算区域的分子信息送出去，接收来自其它区域的分子信息，使模拟区域始终保持完备统一。
4　算例计算
4．1　基本假设
　　卫星和空间站的调姿发动机的喷管一般为轴对称拉伐尔喷管，喷管扩张段为锥形，由于喷管外的环境为真空环境，因而喷管羽流为自由膨胀流动，没有激波的出现。由于喷管内外流动均为较低密度流动，分析时作如下基本假设：(1)流场中分子的碰撞均为二体碰撞；(2)不考虑分子内自由度激发和化学非平衡效应；(3)假设喷管壁面为恒温壁面；(4)气体流动为定常流动；(6)将分子视为硬球分子，即分子间作用力忽略不计。


图2　真空羽流场并行计算程序流程图
　　由于该流动问题包含连续介质流、过渡流和分子自由流三种流动，因此需用分子气体动力学的理论来分析。为了数值模拟上的方便，在上述基本假设下，我们统一用DSMC方法来数值模拟真空下低密度羽流场。喷管喉部流动参数由给定条件通过差分求解N-S方程数值计算喷管流场给出［5］。
4．2　计算结果及分析
　　用DSMC方法并行模拟真空下发动机喷管的羽流场，模拟时采用R-S-F方法［6］。并行采用宏观无重叠区域分解法，程序流程如图2。为了验证并行计算结果的正确性，首先计算一个与串行计算相一致的算例，考虑单计算机CPU和内存的限制，算例1的喷管尺寸和原始数据按表1所示。网格划分喷管内为轴向×径向＝300×200，喷管外为150×400。初始时每个网格内放置4个分子。当模拟区域分子总数变化很小时，认为流动达到稳定状态，此后取时间平均使统计涨落减小到要求的范围内即可。
　　并行计算与串行计算结果对比如图3，由此可看出并行计算与串行计算结果是一致的。
表1　算例1喷管几何尺寸及流动参数

喉部直径D面积比Ae／At (出口／喉部) 扩张角θ气体总温T *总压P*喷管壁温Tw
1.0×10-3 m100 20°N21000 K9200 Pa1000 K

　
图3　串、并行计算轴向压强曲线对比　　　图4　并行加速比和并行效率与节点数的关系
　　算例1的计算验证了并行计算结果的正确性，在模拟计算过程中，同时得出并行程序单位计算量的计算时间、并行加速比及并行效率等与计算节点数的关系，如图4所示。由图4可以看出，通过对DSMC法的并行处理，计算速度明显加快，计算时间减少，同时还可以看出计算节点增多并行加速比始终在增大，呈现出较好的可扩展性。


图5　算例2等压强线、等温度线、等马赫数线及速度矢量图
　　在算例1验证的基础上，不改变喷管的几何尺寸和形状，根据现有能够使用的并行机条件，对气体较高总压、较大区域的真空羽流场进行并行数值模拟，所给算例2气体总压为92000 Pa，是算例1的10倍，轴向计算范围由算例1的0.025 m增加到0.073 m，计算网格数增加到260000个，计算结果分别如图5所示。从图5可以看出计算结果与真空羽流规律是一致的［7］。
5　结束语
　　本文实现了真空小喷管羽流场的并行数值模拟。在分布式并行机上进行并行计算，通过兼顾计算规模和计算速度的协调提高，能够使DSMC法进行真空羽流场的数值模拟具有更强的适应性，为真空羽流场的进一步研究打下良好的基础。本文提出的并行化原则与二次过滤法是对并行工作的经验总结，可为大型计算问题在分布式并行机上实施并行化处理提供参考。并行计算能够通过牺牲计算速度来换取大的计算规模。随着更为先进的并行计算机的推广使用，羽流场的并行数值模拟将具有更强的适应性。
*本文系国家863高技术资助项目
作者简介：**男　31岁　博士　副教授　北京航空航天大学宇航学院　100083 
作者单位：北京航空航天大学
参考文献
　　1　Bird G A.Molecular Gas Dynamics.Clarendon Press,Oxford,1976
　　2　Furlani T R,Lordi J A.A Comparison of Parallel Algorithms for the Direct Simulation Monte Carlo MethodⅡ：Application to Exhaust Plume Flowfields.AIAA-89-1667
　　3　Wilmoth R G,Carlson A B.DSMC Analysis in a Heterogeneous Parallel Computing Environment.AIAA 92-2861
　　4　孙家昶，张林波，迟学斌，汪道柳．网络并行计算与分布式编程环境．北京：科学出版社，1996
　　5　蔡国飙．真空羽流问题的数值模拟研究：［博士论文］．北京航空航天大学，1995
　　6　樊菁，沈青．直接模拟蒙特卡罗方法中的一个新算法—随机取样频率法．第六届全国计算流体力学会议论文集，北京：科学出版社，1992
　　7　蔡国飙，王慧玉，祖国君，张振鹏．真空小喷管羽流场的Monte Carlo直接模拟．推进技术，1997，8(2)：44
（责任编辑　杨再荣）
1998年6月收稿；1998年10月收到修改稿。
