神舟飞船平移发动机对太阳能帆板羽流效应数值模拟研究

蔡国飙 贺碧蛟

（北京航空航天大学宇航学院）

1 引言

飞船、空间站等大型航天器具有质量大、舱段多的特点，其上需要布置多块太阳能帆板才能满足航天器对电力的需求。由于帆板数量较多，对与部分发动机距离较近的帆板，发动机真空羽流会对其产生较大的气动力和力矩作用。为了对发动机真空羽流效应进行数值模拟，必须针对不同的流动状态建立相应的数学物理模型，并采用与之相应的数值模拟方法。目前能够比较成功地模拟真空羽流效应的研究方法是直接模拟蒙特卡罗方法[1]（DSMC）。

国外针对发动机羽流问题的研究开展得比较早，开发了一系列基于DSMC的模拟计算软件，其中比较典型的通用计算软件有：由美国约翰逊航天中心（Johnson Space Center）的 LeBeau[2-4]等人开发的DAC（DSMC Analysis Code）软件，由俄罗斯理论与应用机械研究所（the Institute of Theoretical and Applied Mechanics）的 Ivanov[5]等人开发的 SMILE（Statistical Modeling In Low-density Environment）软件，由美国Cornell大学的Dietrich和Boyd等人开发的MONACO 计算软件[6，7]等。

文章所介绍的基于 DSMC的 PWS（Plume WorkStation）羽流计算软件是依据模块化的原则，适当引入面向对象的编程思想，特别设计具有较强通用性的粒子边界处理模块，实现PWS软件在羽流模拟计算的通用化，并在国外“CUBRC”试验研究的基础上开展试验验证工作。使用PWS软件就神舟飞船120N平移发动机羽流对太阳能帆板气动力效应进行了数值模拟研究，并且对比了由不同帆板角度对作用在帆板中心点的力和力矩影响。

2 羽流数值模拟模型和方法

流场依据Knudsen（Kn）数划分为三大领域，即连续流区领域（Kn＜0.1）、过渡领域（0.1＜Kn＜10）和自由分子流领域（Kn＞10）[8]，Kn= λ/L，其中 λ 是平均分子自由程，L为流动特征长度。在Kn数较大，即气体变稀薄时，基于连续流介质的N-S方程失效。

真空羽流流场包括了全部流态，即连续介质流、过渡领域流和自由分子流等[1]。羽流核心区为连续流，外围为过渡领域流和自由分子流，羽流流动极为复杂，很难用一种方法来描述。对于连续流动，对应的数学方程是N-S方程；对于自由分子流，对应的数学方程是无碰撞项的Boltzmann方程；而对于过渡领域流的研究，目前尚无完善的理论，与之相对应的数学方程是完全的Boltzmann方程，它是一个非常复杂的非线性积分-微分方程，除了极其简单的情形，它的理论求解极为困难。

为了对羽流流场进行数值模拟，必须针对不同的流动状态建立相应的数学物理模型，并采用与之相应的数值模拟方法。目前能够比较成功地模拟自由分子流和过渡领域流的数值方法是DSMC。DSMC方法是20世纪60年代初由G.A.Bird提出并发展的直接模拟方法，在数学上已被证明与Boltzmann方程的相容。

G.A.Bird[1]提出的DSMC直接从物理实际出发，利用少量的模拟分子代替真实流场内数目众多的气体分子，用计算模拟代替物理过程，经统计平均获得宏观流动参数，达到求解稀薄气体问题的目的。典型的DSMC计算流程图如图1所示。

图1 DSMC计算流程图

由于羽流场DSMC计算域的密度值较低，分析时常作如下假设：（1）流场中分子的碰撞为二体碰撞；（2）由于气体温度不高，仅考虑分子的转动内能，忽略分子的振动内能和分子化学非平衡效应；（3）气体流动为定常流动；（4）将分子视为变径硬球分子（VHS）。

3 PWS软件及验证

DSMC使用大量的“模拟分子”以模拟真实的气体，为此需要较大的存储空间和较高计算速度，这极大地限制了DSMC的应用。DSMC的计算主程序主要有6个主模块：粒子模块（粒子参数存储）、网格模块、粒子进入模块、粒子运动模块、边界模块（主要指固体边界）和粒子碰撞模块。考虑到计算速度的需要，在6个主模块的基础上再增加一个并行模块。在羽流效应计算时，DSMC的各模块信息交换如图2所示。

图2 DSMC计算模块信息交换关系图

根据模块化和面向对象的软件设计方法，PWS软件对DSMC计算程序进行如下改进：

（1）使用松散的软件架构保证各个计算模块的相对独立性，即改进其中一个模块时，基本上不需要改动其他模块；

（2）对各个模块的功能实现采用通用性设计，即在使用单/多组分分子、不同分子模型（碰撞模型和内能模型）、不同固体边界条件时不需要改动源代码。

如图2所示，基于DSMC的PWS羽流计算通用软件主要包括网格模块、粒子参数模块、粒子进入模块、粒子运动及边界模块、粒子碰撞模块、并行模块模块6大模块。在这6大基本模块的基础上，针对各种实际问题（单组分、多组分、多种碰撞模型、多种边界条件），PWS羽流计算软件可以方便地搭建数值模拟应用平台。

软件计算结果的验证是软件编制不可缺少的一环。本文使用了美国Holden[9]等人做的“CUBRC”第7次试车数据作为软件计算验证的对比数据。“CUBRC”实验是一个超音速稀薄气流对双圆锥体的气动力和气动加热效应的测量试验，试验用的双圆锥体构型如图3所示。图4为利用PWS软件对双圆锥体表面压强分布的计算结果与试验测量数据对比图，图5为利用PWS软件对双圆锥体表面热流密度分布的计算结果与试验测量数据对比图，轴向距离X指从双圆锥体的顶点为原点，计算点在双圆锥体的中心线上投影的坐标值（m）。

图3 “CUBRC”试验双圆锥体结构图（单位：mm）

图4 双圆锥体表面压强分布计算值与试验值对比

图5 双圆锥体热流密度分布计算值与试验值对比

从图4、图5可以看出双圆锥体表面压强和热流密度计算值与试验值相符，最大值和最小值的位置一致；除了2个点之外，其他各点的压强计算值相对试验值的偏差都在±20%以内；热流密度计算值相对试验值的偏差都在-30%～+10%范围以内，特别是在双圆锥连接处附近符合得相当好。经过与实验值比对，可以认为PWS计算软件具有较高的置信度。

4 太阳能帆板羽流效应计算

神舟飞船平移发动机推力相对较小，故可直接将发动机出口截面作为DSMC计算的模拟分子入口条件，使用PWS软件直接对发动机三维羽流场进行数值模拟，得到发动机单机/双机工作时羽流对太阳能帆板在0°、45°和90°三种状态下产生的气动力效应值，具体计算算例如表1所示。发动机与太阳能帆板的相对位置如图6所示。

表1 对比算例

图6 发动机与太阳能帆板相对位置

4.1 计算条件

（1）粒子入口条件

粒子入口条件为：以发动机出口截面作为PWS软件计算外喷流和后流区的粒子入口截面条件；“模拟粒子”选取燃烧产物主要组分 N2、H2O、CO、CO2和H2等5种气体分子，质量分数如表2所示。

表2 发动机燃烧产物主要组分质量分数

粒子入口截面物理相对值分布如图7所示，图中各参考量的的定义为：参考密度ρ0=0.005kg/m3，参考压强P0=2000Pa参考温度T0=3000K，参考轴向速度u0=3500m/s，参考径向速度v0=500m/s，发动机出口半径R0=0.039m。

图7 粒子入口物理参数分布图

（2）边界条件

真空边界：粒子经过真空边界后逃逸（即删除）；

固体边界：粒子撞击到固体边界按照Maxwell反射[1]处理，Maxwell反射需确定固体边界壁面温度和热适应系数，发动机壁面温度设定为800 K，飞船和太阳能帆板表面温度设定为300 K，热适应系数统一设定为1。

4.2 计算结果与分析

图8至图 13为算例1至算例6太阳能帆板所受的发动机羽流气动压强分布图，图中P指的是太阳能帆板表面正压强（Pa）。

图8 算例1太阳能帆板受到的羽流气动压强分布图

图9 算例2太阳能帆板受到的羽流气动压强分布图

图10 算例3太阳能帆板受到的羽流气动压强分布图

图11 算例4太阳能帆板受到的羽流气动压强分布图

图12 算例5太阳能帆板受到的羽流气动压强分布图

从图8至图 13可以看出T1和T2平移发动机羽流对太阳能帆板的气动压强分布：T2发动机单机工作时羽流对0°状态太阳能帆板气动压强最大值达到12Pa，且压强值在1Pa以上的作用面积较小；T2发动机单机工作时羽流对45°状态太阳能帆板气动压强最大值达到10Pa，且压强值在1Pa以上的作用面积较大；T2发动机单机工作时羽流对90°状态太阳能帆板气动压强最大值仅达到2.4Pa，但压强值在1Pa以上的作用面积较大；T1和T2发动机双机工作时羽流对0°状态太阳能帆板气动压强最大值依然为12Pa，且压强值在1Pa以上的作用面积与单机工作相比有所增加，但还是较小；T1和T2发动机双机工作时羽流对45°状态太阳能帆板气动压强最大值达到10Pa，且压强值在1Pa以上的作用面积很大；T1和T2发动机双机工作时羽流对90°状态太阳能帆板气动压强最大值达到7.5Pa，且压强值在1Pa以上的作用面积很大。这些说明，最大气动力效应工况应该出现在双机工作的45°或90°位置。表 3为六种工况下作用点为太阳能帆板的中心点的力和力矩的值。

图13 算例6太阳能帆板受到的羽流气动压强分布图

从表3可以看出，发动机羽流对太阳能帆板影响最大的工况出现在算例5，即双机工作时帆板处于45°角位置时，其作用在帆板中心点上的合力为约为24N，合力矩约为7.7Nm。

5 结束语

PWS软件是一款基于DSMC的羽流计算通用软件，并且具备了可扩展特性，松散的架构为软件的进一步完善提供了方便的接口，软件的计算值与国外“CUBRC”试验测量值符合得很好，认为PWS软件具有较高的置信度。

使用PWS软件对神舟飞船平移发动机羽流效应进行计算表明：发动机羽流对着太阳能帆板羽流气动压强最大值出现在的算例1和算例4，即单机和双机工作时帆板处于0°角位置时，但作用力不大；发动机羽流对太阳能帆板影响最大的工况出现在算例5，即双机工作时帆板处于45°角位置时，其作用在帆板中心点上的合力约为24N，合力矩约为7.7Nm。

平移发动机对处于羽流直接撞击区的太阳能帆板有较大的力效应作用，需要结合飞船质心进行进一步分析，数值模拟结果可以为神舟飞船平移发动机使用模式提供一定参考作用。 ◇

［1］Bird G A.Molecular gas dynamics［M］.Oxford:Clarendon Press，1976

［2］Wilmoth R G，Carlson A B，LeBeau G J.DSMC grid methodologies for computing low-density.hypersonic flows about reusable launch vehicles，AIAA 1996-1812［R］

［3］Boyles K A，LeBeau G J，Lumpkin III F E.The use of virtual subcells in DSMC analysis of orbiter aerodynamics at high altitudes upon reentry.AIAA 2003-1030［R］

［4］LeBeau G J，Boyles K A，Lumpkin III F E.Virtual sub-cells for the direct simulation Monte Carlo method.AIAA 2003-1031［R］

［5］Ivanov M S，Markelovf G N，Gimelshein S F.Statistical simulation of reactive rarefied flows:numerical approach and applications.AIAA 1998-2669［R］

［6］Dietrich S，Boyd I D.A scalar optimized parallel implementation of the DSMC method.AIAA 1994-355［R］

［7］Dietrich S，Boyd I D.Parallel implement at ion on the IBM SP-2 of the Direct Simulation Monte Carlo method.AIAA 1995-2029[R]

［8］沈青.稀薄气体动力学［M］.北京：国防工业出版社，2003

［9］Holden M，Wadhams T.Code validation study of laminar shock/boundary layer and shock/shock interactions in hypersonic flow.Part A:experimental measurements，AIAA 2001-1031［R］