一种卫星双组元推力器羽流热影响分析方法

徐春生 李洋 肖应廷 曹鹏

（中国空间技术研究院，北京 100094）

1 引言

卫星推力器羽流会影响卫星的工作性能和寿命，在卫星布局阶段须考虑羽流对周围设备产生的影响。目前，分析卫星推力器羽流影响的主要方法包括理论计算和对试验数据的拟合。理论计算一般采用DSMC方法［1-2］，N-S方程和DSMC 方法相结合的方式［3］，或采用点源分析模型［4］。欧洲航天局（使用的Plume软件［5］）、北京航空航天大学［3］、上海交通大学［6］、中国空间技术研究院［1-2］等主要采用上述理论计算方法。理论计算方法在进行羽流分析时，对连续流使用N-S求解，对自由分子流采用二维或三维DSMC 方法进行分析。开展羽流影响分析的另一个途径是对试验数据的拟合。1983年，德国汉堡科技大学开展了MBB双组元10N 推力器的点火试验，综合考虑了推力器点火脉冲等因素，获得了羽流在各阶段的形状、热流量分布、干扰力和冲量分布等数据［7］。文献［8］中采用经过试验数据支持的分析方法，对国际海事卫星-2（Inmarsat-2）因10 N 推力器羽流而受到的冲击和热影响进行了分析。文献［9］在MBB 试验数据的基础上采用数据拟合的方法，得到了双组元10N 推力器的羽流动力学效应和热影响分析的公式。文献［10］和文献［11］在文献［9］的基础上，推导出了羽流对抛物面天线和太阳翼的动力学效应和热影响。

理论计算方法假设推力器连续稳态点火，不考虑点火脉冲的实际情况；而对试验数据拟合的分析方法，所得的结果与实际情况具有更好的一致性。另外，对试验数据进行拟合来计算羽流热影响的方法，计算量非常小，可以快速获得结果。不过，文献［10］的方法和Plume软件，均难以计算形状不规则或无法用解析公式表达的复杂曲面受到的羽流影响。例如：卫星上的通信天线，其反射面一般为赋形抛物面或赋形椭球面，与理论抛物面或理论椭球面相比存在几何上的差别。若用文献［10］的方法，即将赋形反射面简化为标准反射面，会存在结果上的偏差。本文以文献［9］中推导的公式为基础，给出一种分析星上设备受推力器羽流热影响的方法，直接从设备的三维模型提取出STL数据，不再受设备形状是否适合建模的束缚。相对于文献［10］的方法，本文提出的方法除对形状规则设备或能用解析公式表达设备的分析达到同样的效果外，还能分析文献［10］中方法不能分析的对象（无法用解析公式表达的设备）受到的羽流热影响，扩大了羽流热影响分析的应用范围。另外，本文所用方法只需要设备的三维模型，无需复杂的参数设定，计算过程较简单。

本文的方法适用于不同推力的双组元推力器。由于文中采用的羽流热流场分布模型来自对MBB双组元10N 推力器测试数据的拟合（即文献［9］所推导公式），因此只针对10N 双组元推力器的羽流热影响进行了仿真。

2 双组元推力器羽流冲击模型

对星上设备产生热影响的羽流分布主要与以下3个参数有关：从推力器出口到设备表面的距离矢量r，羽流场的分布角κ，以及羽流场与设备表面法线之间的夹角β，见图1［10］。

双组元推力器羽流热流场的空间分布可用函数H（·）表示，以10N 双组元推力器为例，其热流场模型为［9］

式中：r0为试验中所取距离，m；H（r0）为对应r0的热流量试验结果，W／cm2；l为曲线拟合时的参数，当κ＞40°时，取l＝12。

图1 双组元推力器羽流冲击模型Fig．1 Plume impingement model of bipropellant thruster

3 基于STL数据的羽流热影响分析

STL是美国3DSystem 公司1988年开发的数据交换标准［12］。符合该标准的数据文件，是一系列离散的三角形网格的集合。STL 数据文件记录了每个三角形网格的法向矢量和三角形的顶点位置，如图2所示。图2中，P1、P2、P3分别为STL 数据中第i个三角形网格的3个顶点（括号中为顶点坐标），n为三角形网格的法向矢量。

图2 STL数据三角形网格Fig．2 Triangle mesh of STL data

在进行羽流热影响分析时，首先，从设备的三维模型提取STL数据；然后，根据分析的需要，将STL数据文件中的有效三角形网格提取出来，剔除掉不需要分析的三角形网格；最后，以所提取的三角形网格为对象，建立其所受羽流热影响的分析模型。整个方法的流程见图3。

图3 基于STL数据的羽流热影响分析流程Fig．3 Flow chart of plume thermal effect analysis based on STL data

3．1 提取有效的STL数据

在三维软件中，将设备的三维模型转换为STL格式的数据文件（即图2所示三角形网格的集合）。在这个过程中，是将整个模型的外表面都网格化，而不能只将某些特定曲面进行网格化处理。推力器和设备的相对位置关系，决定了羽流只会对整个模型的某些特定曲面产生影响（不考虑热传导的影响）。因此，如何从整个STL数据文件中提取出受羽流影响曲面的三角形网格，成为羽流分析的前提。

在生成STL数据文件时，须选择基准坐标系，本文选择卫星的本体坐标系。推力器和设备的位置信息都在该坐标系下描述。根据这个特点，可以提取STL数据中每个三角形网格的法向矢量n，比较它与卫星本体坐标系某一轴的夹角，如图4所示。

图4 三角形网格法向和卫星本体坐标系坐标轴关系示意图Fig．4 Relation between normal vector of triangle mesh and coordinate system of satellite

从图4可以看出，若所要分析的曲面为设备的＋Z面，则在STL数据中提取出n（Z）＞0（即三角形网格法向和卫星本体坐标系Z轴夹角小于90°）的三角形网格，n（Z）为三角形网格的法向矢量在卫星本体坐标系Z轴上的分量。

3．2 羽流热影响分析建模

提取出所要分析表面的三角形网格后，开始根据羽流热流场分布函数进行羽流分析的建模。

设备表面到推力器出口的距离矢量

式中：P10N为所分析的推力器在卫星本体坐标系下的位置矢量；P为STL数据中每个三角形网格的几何中心在卫星本体坐标系下的位置矢量。

羽流场的分布角为

式中：E为推力器主轴在卫星本体坐标系下的方向矢量。

羽流场与设备表面法线之间的夹角为

由式（1）～（4），可以得到羽流热流场的空间分布如下。

按照式（5），求得每个三角形网格几何中心点处的羽流场热流密度，然后与三角形网格的面积相乘，累加后即得出所要分析的整个设备表面受到的热影响，见式（6）。

式中：N为所要分析设备表面的三角形网格总数；Pi为STL数据中第i个三角形网格的几何中心在卫星本体坐标系下的位置矢量；ri为STL 数据中第i个三角形网格几何中心到推力器出口的距离矢量；

4 应用

以3个应用实例对本文方法与文献［10］方法、Plume软件作一比较。

4．1 与文献［10］方法比较

图5是分别采用本文方法和文献［10］方法获得的卫星太阳翼由于推力器羽流场而受到热影响的分布图。因本例中分析的是受推力器羽流影响的太阳翼＋Z面，所以提取STL数据中n（Z）＝1的三角形网格作为有效网格。两种方法的输入参数一致，太阳翼网格化参数相当（文献［10］采用面积为238．3mm2的长方形网格，本文采用面积为237．7mm2的三角形网格）。从计算结果来看，热流密度分布结果一致性良好，文献［10］得到的最大热流密度为16．72 W／m2，整体所受热影响为199．5 W；本文方法得到最大热流密度为17．2W／m2，整体所受热影响为187．0 W。

图5 太阳翼上羽流热影响分布Fig．5 Plume thermal effect distribution on solar array

图6（a）是采用本文方法分析某赋形抛物面天线的反射面获得的羽流热影响分布图。文献［10］方法并不能分析这种无法用解析公式表达的不规则曲面，只能将其简化为标准抛物面进行分析，分析结果见图6（b）。本例中分析的是受推力器羽流影响的通信天线反射面的下表面，因此提取STL 数据中n（Z）＜0的三角形网格作为有效网格。从分析结果来看，本文方法得到的最大热流密度为1．242kW／m2，最大热流密度作用区域为（-18．5mm，1 470．5 mm，4 420．5mm）、（-12．7mm，1 473．8mm，4 423．0mm）、（-10．8mm，1 463．9mm，4 416．9mm）围成的面积为37．6mm2的三角形区域；整个反射面下表面受到羽流影响的热容量为174．8 W。文献［10］方法得到的最大热流密度为1．166kW／m2，整个反射面下表面受到羽流影响的热容量为106．1 W。

比较图6（a）和图6（b）可见，本文方法和文献［10］方法得到的热流密度分布相同，只是由于图6（a）所示模型中受到羽流热影响较大区域的形状不规则，比图6（b）所示的标准抛物面更靠近推力器，因此用前者得出的最大热流密度比后者的要大，更符合实际情况。

图6 通信天线反射面上羽流热影响分布Fig．6 Plume thermal effect distribution on communication antenna’s reflector

4．2 与Plume软件比较

图7是分别采用本文方法和Plume软件得到的某推力器羽流对天线展开臂转接头部位的热影响分布。该转接头结构较复杂，文献［10］方法无法建立羽流分析模型。由于Plume软件难以对复杂的模型建模，因此仅建立转接头下半部分（热流密度最大部位）的模型。图7（a）和图7（b）显示出热流密度分布相同，从热流密度最大值来看，本文方法所得结果（168kW／m2）和Plume软件得到的结果（172kW／m2）仅相差2．4%。这说明该天线展开臂转接头部位局部热流密度过大，须采取相应热控措施，或规划推力器工作时机来加以解决。

图7 天线展开臂转接头上羽流热影响分布Fig．7 Plume thermal effect distribution on antenna’s deployed arm connector

5 结论

本文以文献［9］推导出的羽流热流场分布函数为基础，提出一种基于STL数据的羽流热影响分析方法，主要利用对推力器试验数据的拟合实现，计算量较小。通过与文献［10］方法以及Plume软件的对比，可以得出如下结论。

（1）本文方法与文献［10］方法和Plume软件相比，在分析同样设备受到的羽流热影响时，具有良好的一致性。

（2）本文方法除了能分析形状规则设备或能用解析公式表达设备受到的羽流热影响，还能分析形状不规则设备或不能用解析公式表达设备受到的羽流热影响。这比文献［10］和Plume软件具有一定的优势。

本文所介绍的方法，包括但不限于双组元10N推力器的羽流热影响分析。只要提供类似文中式（1）的羽流热流场分布函数，即可分析其他类型推力器的羽流对星上设备的热影响。

（References）

［1］孙威，闫少光，田东波．卫星发动机羽流污染模型的DSMC方法研究［J］．航天器环境工程，2007，24（2）：95-98

Sun Wei，Yan Shaoguang，Tian Dongbo．The DSMC method for plume contamination modeling of satellite［J］．Spacecraft Environment Engineering，2007，24（2）：95-98（in Chinese）

［2］张健，尚志．航天器发动机羽流对敏感器热效应仿真研究［J］．航天器工程，2009，18（3）：59-65

Zhang Jian，Shang Zhi．Simulation of plume thermal effect on spacecraft sensitive devices［J］．Spacecraft Engineering，2009，18（3）：59-65（in Chinese）

［3］蔡国飙，贺碧蛟．PWS软件应用于探月着陆器羽流效应数值模拟研究［J］．航天器环境工程，2010，27（1）：18-23

Cai Guobiao，He Bijiao．PWS software applied to numerical simulation on plume effect of lunar exploration lander［J］．Spacecraft Environment Engineering，2010，27（1）：18-23（in Chinese）

［4］程晓丽，王强，阎喜勤．卫星姿控发动机高空羽流场工程分析［J］．空间科学学报，2006，26（1）：59-63

Cheng Xiaoli，Wang Qiang，Yan Xiqin．Engineering analysis on a satellite attitude control thruster plume at high altitudes［J］．Chinese Journal of Space Science，2006，26（1）：59-63（in Chinese）

［5］恒润科技公司．SYSTEMA——航天器空间任务及在轨分析工具［EB／OL］．［2012-05-12］．http：／／www．hirain．com／user-community／download／technical-information／systema

Hirain Technologies．SYSTEMA-analysis tool for spacecraft space operation on-orbit［EB／OL］．［2012-05-12］．http：／／www．hirain．com／user-community／down-load／technical-information／systema（in Chinese）

［6］王平阳，鹿洪芳，程惠尔．双组元姿控发动机反流特性的DSMC模拟［J］．工程热物理学报，2004，25（4）：640-642

Wang Pingyang，Lu Hongfang，Cheng Huier．DSMC simulation of plume backflow from bipropellant attitude control engine ［J］．Journal of Engineering Thermophysics，2004，25（4）：640-642（in Chinese）

［7］Trinks H．Experimental exhaust plume analysis with MBB 10N thruster［C］／／Proceedings of AIAA／SAE／ASME 19thJoint Propulsion Conference．Washington：AIAA，1983

［8］Trinks H，Fox S．Evaluation of the exhaust plume of the bipropellant MBB 10Nthruster for the Inmarsat 2 satellite ［C］／／Proceedings of AIAA／ASME／SAE／ASEE 25thJoint Propulsion Conference．Washington：AIAA，1989

［9］孙宝祥．利用10N 推力器羽流试验数据建立羽流场数学模型［J］．航天控制，2005，23（3）：26-29

Sun Baoxiang．Modeling plume field using plume test data for 10Nthruster［J］．Aerospace Control，2005，23（3）：26-29（in Chinese）

［10］Li Yang，Li Xuelin．Analyses and simulations of exhaust plume effects on satellite appendages［C］／／Proceedings of IEEE ICCSE．New York：IEEE，2010：184-188

［11］Xia Yongquan，Yin Pu．Simulation of thruster’s plume effects on satellite［C］／／Proceedings of IEEE IACSIT ICCEA．New York：IEEE，2011：604-608

［12］Kai C C，Jacob G G，Mei T．Interface between CAD and rapid prototyping systems，part 1：a study of existing interfaces［J］．International Journal of Advanced Manufacturing Technology，1997，13（8）：566-570