梁 杰，阎 超，李志辉，李绪国

（1．北京航空航天大学 航空科学与工程学院，北京 100191；2．中国空气动力研究与发展中心，四川 绵阳 621000）

0 引 言

飞行器在飞行过程中因机动需要快速改变飞行姿态或轨道，通常采用反作用控制系统（RCS）来实现，分布在飞行器表面不同位置的RCS发动机提供的推力完成俯仰、偏航和滚转等姿态变化以及轨道的转移［1］。RCS发动机工作时产生的喷流气体与周围气流发生相互作用，对推力和飞行器的压心位置会产生附加影响，并改变飞行器表面的热环境。在连续流区域，经过几十年持续的风洞试验和数值模拟研究［2－4］，对不同类型的喷流干扰及流动机理有了非常细致的了解。但在稀薄大气环境下，尤其是飞行器处于稀薄过渡流域时，由于连续介质假设失效，飞行器表面出现了滑移效应，喷流与前方边界层（粘性层）相互作用引起的流动分离和流场结构均出现了不同程度的变化。无论是地面风洞试验所用的相似准则，还是数值求解连续介质假设下的NS方程，都已无法满足精确模拟稀薄过渡流区横向喷流干扰特征的要求。另外，在高空稀薄气流中喷流会形成羽流，它所涉及的问题还包括对周围环境的污染［5］，影响目标的识别。在文献［6］中，我们详细研究了三维平板模型在不同的喷流压力比环境下音速喷流与稀薄大气干扰的流动规律，计算得到的干扰流场结构、流动分离特征和分离长度与低密度风洞辉光放电和油流测量结果有较好的一致性，确认了建立的直接模拟Monte Carlo（DSMC）算法的可靠性和计算精度。由于旋成体外形的飞行姿态不同，可能存在攻角，也可能还存在侧滑角，飞行器外流场的复杂性导致旋成体上横向喷流干扰流场较平板的喷流干扰要复杂得多。更为重要的差别在绕旋成体流动的分离边界，很可能从上表面绕到下表面，产生所谓的包裹效应［4］。这些复杂的横向喷流流动现象难以进行风洞试验模拟，都需要依靠数值模拟开展细致的研究分析。本文在对平板模型喷流干扰研究的基础上，重点开展了细长钝双锥外形高超声速稀薄来流与超声速／高超声速横向喷流干扰效应的研究，计算分析了不同高度、不同速度、不同攻角、不同喷流推力下复杂流场结构和对气动力特性的影响规律。考察了RCS喷管出口参数不同（均匀／非均匀）对喷口附近分离涡和分离长度的影响。计算中采用了研究建立的两级直角网格及碰撞和取样网格自适应算法［7］，以及MPI并行环境下的基于静态随机负载平衡技术的DSMC并行算法。

1 DSMC数值方法

1．1 DSMC方法

DSMC方法［8］是用若干个模拟分子代替真实气体的分子，并在计算机上存贮模拟分子的位置坐标、速度分量以及内能，这些量随模拟分子的运动、与边界的碰撞以及分子之间的碰撞而改变，最后通过统计网格内模拟分子的运动状态实现对真实气体流动问题的模拟。该方法的关键是在时间步长Δt内将分子的运动与碰撞解耦。在Δt时间内让所有分子运动一定的距离并考虑在边界的反射，然后计算此Δt内具有代表性的分子间的碰撞。

1．2 计算网格及网格自适应方法

在DSMC方法中，流场中的网格是用来选取可能的碰撞分子对以及对宏观流动参数取样。流场采用均匀的直角坐标网格，追踪分子的效率非常高，计算区域内的模拟分子可以直接根据分子的位置坐标来确定分子所属的网格，而不必跟踪分子从一个网格运动到另一个网格［7］。其缺点是无法精确地描述物面边界。本文结合直角网格计算的高效率以及表面非结构网格对飞行器几何外形的精细描述，建立了混合网格结构的DSMC数值模拟方法。在描述物面几何形状的非结构网格建立以后，直接将其嵌入到直角网格的流场中，使DSMC计算对流场网格的依赖程度大大降低。同时通过判断分子运动轨迹方程和物面三角形面元上任一点的位置方程，唯一确定出分子与物面的碰撞点坐标，解决了这种混合网格流场分子运动与物面碰撞的难题。对分子在物体三角形面元上碰撞、反射前后的流场参数进行统计取样就可以获得飞行器的整体气动力特性以及表面力、热载荷分布。

稀薄过渡流区的喷流干扰流场结构非常复杂，流场中的气体密度、压力有几个量级的变化。喷流核心区和激波干扰区的密度非常高，其它大部分区域的密度相对较低。为了解决流场中密度急剧变化的流动特征，采用了网格自适应的方法。即在背景网格的基础上，根据流场中密度梯度的变化分别对碰撞网格和取样网格进行细化。由于流场的梯度沿各个方向的变化是不同的，梯度变化大的方向网格细分的更密一些，因此沿三个坐标方向是各自独立地进行自适应，碰撞分子则是在自适应后最小的亚网格内选取，保证了计算的空间精度。

碰撞网格自适应的目的是使亚网格的尺度与当地平均自由程一致，取样网格自适应是为了捕捉当地的流动梯度，正常情况下取样网格的数量远远小于碰撞网格的数量。为了快速定位自适应后分子所属的网格，采用如下的搜索方法：

对于流场中均匀的矩形网格，模拟分子在背景网格中的位置可以用下面的编号表示：

式中ij是在（x，y，z）三个坐标方向的序号，rj是分子坐标，Δj＝Lj／Nj是网格间距，Lj是沿相应坐标方向的流场大小，Nj是沿相应坐标方向划分的背景网格数，j0是背景网格的起始坐标。则模拟分子在背景网格中的序号为：

每个背景网格根据当地的流动梯度沿三个坐标方向分别划分成lx，ly，lz个二级碰撞或取样网格，二级网格在一级（背景）网格中的起始编号为：

因为二级网格也是均匀的矩形网格，因此二级网格的序号也可以按照公式（2）同样的方式获得：

1．3 随机负载平衡的DSMC并行算法

采用DSMC方法模拟三维RCS喷流与稀薄大气的干扰流场，由于喷流核心区的气流密度相当高，导致模拟分子数和网格划分引起的计算量仍然非常大，需要开展并行计算。DSMC并行算法采用区域分解，根据计算的处理器数（或CPU核心数）的多少将计算区域划分为等量的子区域，每个处理器在其分配的子区域内部独立地计算模拟分子的碰撞和迁移，离开子区域的模拟分子把携带的信息传递给对应子区域的处理器。并行计算总的时间包括每个处理器计算碰撞、迁移的时间、处理器之间的通讯时间以及各处理器之间为同步而等待的时间。提高并行计算的效率主要是通过减少通讯和同步等待的时间来实现。

如果模拟分子在一个时间步长之后，离开原来所在的网格迁移进入相邻的网格，并且改变后的网格隶属于其它的处理器，这样的分子及其携带的所有信息就被记录下来并放入缓冲数组，当所有的模拟分子都被确定以后，放入缓冲数组的分子就在相应的处理器之间进行信息交换。

DSMC方法的一个显著特征是模拟分子在流场中的分布是不断变化的，计算时间依赖于子区域内的分子个数、分子的碰撞次数以及分子与物面的碰撞次数。因为模拟前不可能预测到不同子区域内的分子个数，因此在区域分解时将各子区域的计算负载划分得比较均衡就特别困难。如果简单地将流场均匀划分成若干同等大小的子区域，就有可能造成处理器之间的负载分布不均匀，负载最少的处理器就要化费较长的时间等待负载最大的处理器，造成同步等待的时间显著增加。这种方式并不适合DSMC并行模拟。

为了减少不同处理器计算时间的差别，通常采用负载平衡技术来实现［9－10］，本文根据稀薄大气条件下飞行器横向喷流干扰特点，发展了一套静态随机负载平衡方法用于解决不同处理器之间的计算时间同步问题，该方法基于概率近似原理，将计算区域的全部网格平均分配给指定的所有处理器。由于采用相同的概率随机选取流场网格，按照均分后的数量分配给每个处理器，因此当计算区域的网格数量较多时，每个处理器包含近似相等的物体边界、高密度流动区域和稀薄气体区域的网格数，这样每个处理器的计算负载也非常接近。当用于并行计算的CPU数量不多时（Np＜32），这种方式得到的并行算法加速比和并行效率都非常高［6］。

2 计算结果分析

2．1 平板模型喷流干扰计算与试验对比

为了验证本文算法的精确度，计算模拟了三维平板模型横向喷流干扰的分离距离并与低密度风洞试验结果进行了对比。风洞试验气体为氮气，前室总温593K、总压400kPa，流场马赫数12．99。音速喷流的气体为氦气，喷嘴前室总温为300K，模拟的三个试验状态的总压变化见表1。在DSMC数值模拟中，采用VHS分子模型，L－B内能交换模型，壁面采用完全漫反射模型（下同）。平板模型长120mm，宽68mm，高15mm。音速喷嘴出口直径2mm，在平板中心距前端70mm的位置。

喷流干扰的分离长度是通过测量喷管出口中心到主分离线的径向距离得到的。用DSMC方法计算的三个试验状态的分离长度，试验测量的分离长度都列在了表1中。计算与试验总体上吻合较好，只是计算得到的喷流干扰流场分离长度比试验测量的分离长度偏小些，最大偏差3mm（Case3），最大相对偏差10%（Case1）。这是由于喷嘴尺寸较小，为简化起见，计算中音速喷嘴出口参数设置为均匀流，即出口只有轴向速度，径向速度为0。这样会减弱喷流沿径向的膨胀速度，可能会造成由于逆向流动减弱而引起计算的分离距离偏小，喷管出口参数不同对分离长度的影响在下节的算例中还会进一步分析说明。

表1 DSMC计算模拟的试验状态及分离距离Table 1 Test conditions and separation distance

模拟的激波／边界层干扰、流动分离特征以及表面极限流线与低密度风洞辉光放电和表面油流图谱的对比分析详见文献［6］。这些结果验证了本文算法的可靠性以及对稀薄过渡流区复杂横向喷流干扰流场的模拟能力。

2．2 钝双锥横向喷流与稀薄大气干扰特性数值模拟

下面选取总长度2．1m的细长钝双锥外形作为研究对象，喷流中心距头部1m，方向与Y轴正向相同。这样就可以考察头部斜激波与喷流弓形主激波干扰引起的复杂流动结构。RCS发动机分别采用推力为150N和10N的型面喷管，150N喷管的出口马赫数近似为4．4，燃气组分为 H2、N2、NH3和 H2O，10N喷管的出口马赫数6．3，燃气组分同上。模拟的大气参数组分为N2和O2。

2．2．1 飞行高度对喷流干扰流场的影响

图1分别给出了喷流推力150N时85km、90km和100km高度、飞行速度6000m／s的干扰流场密度分布等值线云图。飞行高度较低时，飞行器头部产生较强的斜激波，头激波与喷流弓形主激波相交产生激波／激波空间干扰，喷流与前方边界层（粘性层）内气体相互作用形成三维分离流动区。喷流弓形激波后为高压力和高密度的强剪切层区，由于喷流阻碍了来流，在喷流的下游存在较宽的低密度和低压区。

随着飞行高度的增加，来流气体的稀薄度增大，头部斜激波变弱，喷流与稀薄来流的干扰也越来越弱，喷流自身的影响区域却越来越大，在100km以上高空，喷流的影响已经覆盖整个飞行器。如果喷流的燃气组分中含有较强污染效应的成分（如NH3和H2O），会粘附在飞行器表面，从而影响传感器的正常工作。另外，随着飞行高度的增加，喷流与粘性层内气体的干扰减弱，引起的分离也逐渐减小，100km的计算状态已经没有分离的形成。

图1 不同高度150N推力横喷干扰流场密度分布Fig．1 Density distribution of lateral jet interaction flowfield of 150Nthrust in different altitudes

2．2．2 喷管出口参数分布对喷流干扰流场影响

为方便起见，前面的计算中喷管出口参数采用的是均匀分布，沿径向的气体速度为零，但实际上燃气在喷管内流动到达出口时分布并不是均匀的。图2、图3分别给出了采用真实喷管非均匀出口参数（求解NS方程获得）和均匀分布计算的90km高度喷流干扰分离区附近的流线局部放大图，不论是主分离涡还是马蹄涡均有明显差别，采用真实喷管型面出口参数计算的分离区明显大于采用均匀出口参数计算的分离区。

另外，从图4、图5分别绘出的采用喷管出口非均匀与均匀参数计算得到的表面极限流线分布，也可以看出两者的不同。在图4中数字1代表主分离线，2代表再附线，3是马蹄涡引起的二次分离线，与图2中的空间流线相对应。图5中仅可看到主分离线和再附线，同时主分离线和再附线之间的分离区范围也减小了。因此，喷管出口参数分布对分离区的大小和形状有明显的影响。

2．2．3 来流速度及喷管推力对干扰流场的影响

为了比较飞行器来流速度不同对飞行器横向喷流干扰流场的影响，图6给出了来流速度4000m／s、高度85km、喷流推力150N的干扰流场密度分布，与图1（a）相比，来流速度的降低使头部斜激波减弱，喷流主激波向上游扩展增强，喷流产生的逆向流动引起的分离距离和空间分离涡更大一些。图7、图8绘出了喷流推力为10N时来流速度分别为4000m／s和6000m／s的喷流干扰流场密度分布，飞行高度为85km，喷流马赫数6．3。从图中可以看出，尽管喷流马赫数较高，但由于推力小，喷流产生的逆流与粘性层内气体的干扰较弱，并没有引起明显的流动分离。而且，对于小推力喷流来说，来流速度的增加，使喷流主激波向下游倾斜得更大，对喷口下游区域的再压缩增强，使下游飞行器表面压力增大，对力矩产生的附加影响有所加大。计算表明，对于推力比较小的喷流，在来流参数不变的条件下，来流对喷流的压制作用要明显高于大推力喷流。

2．2．4 攻角变化对喷流干扰流场的影响

旋成体细长双锥外形的飞行姿态不同，如正攻角和负攻角所对应的干扰流场特性是不同的。10°攻角时（高度90km、来流速度6000m／s、喷流推力150N），由于头激波远离喷口，对喷口影响较小，与喷流弓形激波的作用位置上移，粘性层内气体的流动减弱，喷流的逆向流动引起的空间分离涡相比0°攻角时有所增强（图9）。同时下表面的气流向上方卷起，使分离产生的包裹效应减弱，表面分离形成的高压区减小。

图6 150N推力横喷干扰流场密度分布（V＝4000m／s）Fig．6 Density distribution of lateral jet interaction flowfield of 150Nthrust（V＝4000m／s）

图7 10N推力横喷干扰流场密度分布（V＝4000m／s）Fig．7 Density distribution of lateral jet interaction flowfield of 10Nthrust（V＝4000m／s）

图8 10N推力横喷干扰流场密度分布（V＝6000m／s）Fig．8 Density distribution of lateral jet interaction flowfield of 10Nthrust（V＝6000m／s）

图9 10°攻角时喷口上游流线及分离涡Fig．9 Stream line and seperation vortex of 10°AOA

图10 －10°攻角时喷口上游流线分布Fig．10 Stream line and seperation vortex of－10°AOA

图10 是－10°攻角的情况。这时的头激波靠近喷口，并产生较强的压制作用。喷口附近的当地压力比减小，喷流干扰引起的喷口上游空间的分离基本消失，但由于气流是向下表面方向流动，空间流动引起的包裹效应在表面引起的高压力区域有所增大。

2．2．5 喷流干扰对气动力特性的影响

对于旋成体上的侧向喷流干扰引起的流动分离，很可能从上表面绕到下表面跨过φ＝90°的分界线（定义φ＝0°为上表面迎风母线），假如喷流干扰的压力正增量区（在主分离线之后的相当大区域）包裹在旋成体上，那么在φ＜90°范围内的正增量投影将放大喷管推力，对应φ＞90°的干扰区，虽然表面压力仍具有正增量，但其在推力方向的投影为负增量，因此产生所谓的包裹效应。通常用放大系数来定义干扰引起的力特性变化量［4］：

推力放大系数 KF＝（T＋ΔT）／T

力矩放大系数 KM＝（M0＋ΔM）／M0

ΔT与ΔM表示由于喷流干扰引起的推力与力矩增量，本文中力矩参考点为头部顶点。

表2给出了喷管推力为150N时喷流干扰引起气动力系数的增量以及计算的推力和力矩放大系数，来流速度为6000m／s。表中数据表明，在相同的飞行状态下，随着飞行高度的增加，喷流干扰引起的气动力系数的增量在数值上是不断增大的，但推力和力矩放大系数却是减小的趋势，主要是由于来流密度降低使动压减小，造成推力和力矩的增量减小。另外，正攻角引起推力和力矩放大系数的减小，负攻角引起推力放大系数的增大，力矩放大系数的减小。喷管推力10N时喷流干扰对气动力特性的影响，与150N推力喷流干扰的变化规律类似，只是当推力比较小时，推力和力矩放大系数有所增大，主要是处于分母的推力和力矩比较小造成的，与气动中心的Φ1m高超声速风洞的侧向喷流干扰试验有相同的结论。

表2 喷流干扰对气动力特性的影响（150N喷管）Table 2 The effects of jet interaction on aerodynamic properties（150Nthrust）

3 结 论

综合分析上述数值模拟结果，可以得出以下几点初步研究结论：

（1）在喷流参数保持不变的情况下，随着来流稀薄度的增加，因喷流干扰引起的空间分离涡不断减弱直至消失，稀薄来流对喷流的影响越来越弱，而喷流自身的影响区域却越来越大。

（2）采用真实RCS喷管出口参数计算的空间分离涡和表面极限流线与采用均匀出口参数有明显的不同，喷管出口参数分布对分离区的大小、形状以及表面压力分布都有一定的影响。

（3）高超声速来流对小推力的RCS喷流有较强的压制作用，喷口附近的干扰流场结构及分离特征与当地流动条件密切相关。

（4）来流速度的增加，使喷流主激波更向下游倾斜，并对喷口下游的流动再压缩增强。

（5）飞行姿态不同，如正攻角和负攻角所对应的干扰流场特性明显不同，正攻角引起推力和力矩放大系数的减小，负攻角引起推力放大系数的增大，力矩放大系数的减小。

（6）在相同飞行状态下，随着飞行高度的增加，喷流干扰引起的气动力系数的增量在数值上是不断增大的，但推力和力矩放大系数却是减小的趋势。

（7）在稀薄流域，RCS喷流干扰引起的推力和力矩放大系数总体上是比较小的。

［1］ GIMELSHEIN S F，ALEXEENKO A A，LEVIN D A．Modeling of chemically reacting flows from a side jet at high altitudes［R］．AIAA 2002－0212，2002．

［2］ BRANDEIS J，GILL J．Experimental investigation of side jet steering for missiles at supersonic and hypersonic speeds［R］．AIAA 95－0316，1995．

［3］ SRIVASTAVA B．CFD analysis and validation of lateral jet control of a missile［R］．AIAA 96－0288，1996．

［4］ 李素循．激波与边界层主导的复杂流动［M］．北京：科学出版社，2006：162－190．（LI suxun．Shock wave and boundary layer dominant complicated flow［M］．Beijing：Science Press，2006：162－190（in Chinese））

［5］ GIMELSHEIN S F，ALEXEENKO A A，LEVIN D A．Modeling of the interaction of a side jet with a rarefied atmosphere［R］．AIAA 2001－0503，2001．

［6］ 梁杰，阎超，杨彦广，等．过渡区侧向喷流干扰的并行DSMC数值模拟研究［J］．宇航学报，2011，32（5）：1012－1018．（LIANG Jie，YAN Chao，YANG Yanguang，et al．Parallel DSMC simulation of lateral jet interaction in rarefied transitional region［J］．Journal of Astronautics，2011，32（5）：1012－1018（in Chinese））

［7］ 梁杰，阎超，杜波强．基于两级直角网格结构的三维DSMC算法研究［J］．空气动力学学报，2010，28（4）：466－471：466－471．（LIANG Jie，YAN Chao，DU Boqiang．An algorithm study of three－dimensional DSMC simulation based on two－level Cartesian coordinates grid structure［J］．Acta Aerodynamica Sinica，2010，28（4）：466－471（in Chinese））

［8］ BIRD G A．Molecular gas dynamics and the direct simulation of gas flows［M］．Oxford：Clarendon Press，1994．

［9］ IVANOV M，MARKELOV G，TAYLOR S，et al．Parallel DSMC strategies for 3Dcomputation［A］．Proc．Par－allel CFD＇96［C］，1996：485－492．

［10］LEBEAU G J．A parallel implementation of the direct simulation Monte Carlo method［J］．Comput．Methods Appl．Mech．Engrg．，1999，174：319－337．