徐春生 李洋 肖應(yīng)廷 曹鵬

（中國空間技術(shù)研究院，北京 100094）

1 引言

衛(wèi)星推力器羽流會影響衛(wèi)星的工作性能和壽命，在衛(wèi)星布局階段須考慮羽流對周圍設(shè)備產(chǎn)生的影響。目前，分析衛(wèi)星推力器羽流影響的主要方法包括理論計算和對試驗數(shù)據(jù)的擬合。理論計算一般采用DSMC方法［1-2］，N-S方程和DSMC 方法相結(jié)合的方式［3］，或采用點源分析模型［4］。歐洲航天局（使用的Plume軟件［5］）、北京航空航天大學［3］、上海交通大學［6］、中國空間技術(shù)研究院［1-2］等主要采用上述理論計算方法。理論計算方法在進行羽流分析時，對連續(xù)流使用N-S求解，對自由分子流采用二維或三維DSMC 方法進行分析。開展羽流影響分析的另一個途徑是對試驗數(shù)據(jù)的擬合。1983年，德國漢堡科技大學開展了MBB雙組元10N 推力器的點火試驗，綜合考慮了推力器點火脈沖等因素，獲得了羽流在各階段的形狀、熱流量分布、干擾力和沖量分布等數(shù)據(jù)［7］。文獻［8］中采用經(jīng)過試驗數(shù)據(jù)支持的分析方法，對國際海事衛(wèi)星-2（Inmarsat-2）因10 N 推力器羽流而受到的沖擊和熱影響進行了分析。文獻［9］在MBB 試驗數(shù)據(jù)的基礎(chǔ)上采用數(shù)據(jù)擬合的方法，得到了雙組元10N 推力器的羽流動力學效應(yīng)和熱影響分析的公式。文獻［10］和文獻［11］在文獻［9］的基礎(chǔ)上，推導(dǎo)出了羽流對拋物面天線和太陽翼的動力學效應(yīng)和熱影響。

理論計算方法假設(shè)推力器連續(xù)穩(wěn)態(tài)點火，不考慮點火脈沖的實際情況；而對試驗數(shù)據(jù)擬合的分析方法，所得的結(jié)果與實際情況具有更好的一致性。另外，對試驗數(shù)據(jù)進行擬合來計算羽流熱影響的方法，計算量非常小，可以快速獲得結(jié)果。不過，文獻［10］的方法和Plume軟件，均難以計算形狀不規(guī)則或無法用解析公式表達的復(fù)雜曲面受到的羽流影響。例如：衛(wèi)星上的通信天線，其反射面一般為賦形拋物面或賦形橢球面，與理論拋物面或理論橢球面相比存在幾何上的差別。若用文獻［10］的方法，即將賦形反射面簡化為標準反射面，會存在結(jié)果上的偏差。本文以文獻［9］中推導(dǎo)的公式為基礎(chǔ)，給出一種分析星上設(shè)備受推力器羽流熱影響的方法，直接從設(shè)備的三維模型提取出STL數(shù)據(jù)，不再受設(shè)備形狀是否適合建模的束縛。相對于文獻［10］的方法，本文提出的方法除對形狀規(guī)則設(shè)備或能用解析公式表達設(shè)備的分析達到同樣的效果外，還能分析文獻［10］中方法不能分析的對象（無法用解析公式表達的設(shè)備）受到的羽流熱影響，擴大了羽流熱影響分析的應(yīng)用范圍。另外，本文所用方法只需要設(shè)備的三維模型，無需復(fù)雜的參數(shù)設(shè)定，計算過程較簡單。

本文的方法適用于不同推力的雙組元推力器。由于文中采用的羽流熱流場分布模型來自對MBB雙組元10N 推力器測試數(shù)據(jù)的擬合（即文獻［9］所推導(dǎo)公式），因此只針對10N 雙組元推力器的羽流熱影響進行了仿真。

2 雙組元推力器羽流沖擊模型

對星上設(shè)備產(chǎn)生熱影響的羽流分布主要與以下3個參數(shù)有關(guān)：從推力器出口到設(shè)備表面的距離矢量r，羽流場的分布角κ，以及羽流場與設(shè)備表面法線之間的夾角β，見圖1［10］。

雙組元推力器羽流熱流場的空間分布可用函數(shù)H（·）表示，以10N 雙組元推力器為例，其熱流場模型為［9］

式中：r0為試驗中所取距離，m；H（r0）為對應(yīng)r0的熱流量試驗結(jié)果，W／cm2；l為曲線擬合時的參數(shù)，當κ＞40°時，取l＝12。

圖1 雙組元推力器羽流沖擊模型Fig．1 Plume impingement model of bipropellant thruster

3 基于STL數(shù)據(jù)的羽流熱影響分析

STL是美國3DSystem 公司1988年開發(fā)的數(shù)據(jù)交換標準［12］。符合該標準的數(shù)據(jù)文件，是一系列離散的三角形網(wǎng)格的集合。STL 數(shù)據(jù)文件記錄了每個三角形網(wǎng)格的法向矢量和三角形的頂點位置，如圖2所示。圖2中，P1、P2、P3分別為STL 數(shù)據(jù)中第i個三角形網(wǎng)格的3個頂點（括號中為頂點坐標），n為三角形網(wǎng)格的法向矢量。

圖2 STL數(shù)據(jù)三角形網(wǎng)格Fig．2 Triangle mesh of STL data

在進行羽流熱影響分析時，首先，從設(shè)備的三維模型提取STL數(shù)據(jù)；然后，根據(jù)分析的需要，將STL數(shù)據(jù)文件中的有效三角形網(wǎng)格提取出來，剔除掉不需要分析的三角形網(wǎng)格；最后，以所提取的三角形網(wǎng)格為對象，建立其所受羽流熱影響的分析模型。整個方法的流程見圖3。

圖3 基于STL數(shù)據(jù)的羽流熱影響分析流程Fig．3 Flow chart of plume thermal effect analysis based on STL data

3．1 提取有效的STL數(shù)據(jù)

在三維軟件中，將設(shè)備的三維模型轉(zhuǎn)換為STL格式的數(shù)據(jù)文件（即圖2所示三角形網(wǎng)格的集合）。在這個過程中，是將整個模型的外表面都網(wǎng)格化，而不能只將某些特定曲面進行網(wǎng)格化處理。推力器和設(shè)備的相對位置關(guān)系，決定了羽流只會對整個模型的某些特定曲面產(chǎn)生影響（不考慮熱傳導(dǎo)的影響）。因此，如何從整個STL數(shù)據(jù)文件中提取出受羽流影響曲面的三角形網(wǎng)格，成為羽流分析的前提。

在生成STL數(shù)據(jù)文件時，須選擇基準坐標系，本文選擇衛(wèi)星的本體坐標系。推力器和設(shè)備的位置信息都在該坐標系下描述。根據(jù)這個特點，可以提取STL數(shù)據(jù)中每個三角形網(wǎng)格的法向矢量n，比較它與衛(wèi)星本體坐標系某一軸的夾角，如圖4所示。

圖4 三角形網(wǎng)格法向和衛(wèi)星本體坐標系坐標軸關(guān)系示意圖Fig．4 Relation between normal vector of triangle mesh and coordinate system of satellite

從圖4可以看出，若所要分析的曲面為設(shè)備的＋Z面，則在STL數(shù)據(jù)中提取出n（Z）＞0（即三角形網(wǎng)格法向和衛(wèi)星本體坐標系Z軸夾角小于90°）的三角形網(wǎng)格，n（Z）為三角形網(wǎng)格的法向矢量在衛(wèi)星本體坐標系Z軸上的分量。

3．2 羽流熱影響分析建模

提取出所要分析表面的三角形網(wǎng)格后，開始根據(jù)羽流熱流場分布函數(shù)進行羽流分析的建模。

設(shè)備表面到推力器出口的距離矢量

式中：P10N為所分析的推力器在衛(wèi)星本體坐標系下的位置矢量；P為STL數(shù)據(jù)中每個三角形網(wǎng)格的幾何中心在衛(wèi)星本體坐標系下的位置矢量。

羽流場的分布角為

式中：E為推力器主軸在衛(wèi)星本體坐標系下的方向矢量。

羽流場與設(shè)備表面法線之間的夾角為

由式（1）～（4），可以得到羽流熱流場的空間分布如下。

按照式（5），求得每個三角形網(wǎng)格幾何中心點處的羽流場熱流密度，然后與三角形網(wǎng)格的面積相乘，累加后即得出所要分析的整個設(shè)備表面受到的熱影響，見式（6）。

式中：N為所要分析設(shè)備表面的三角形網(wǎng)格總數(shù)；Pi為STL數(shù)據(jù)中第i個三角形網(wǎng)格的幾何中心在衛(wèi)星本體坐標系下的位置矢量；ri為STL 數(shù)據(jù)中第i個三角形網(wǎng)格幾何中心到推力器出口的距離矢量；

4 應(yīng)用

以3個應(yīng)用實例對本文方法與文獻［10］方法、Plume軟件作一比較。

4．1 與文獻［10］方法比較

圖5是分別采用本文方法和文獻［10］方法獲得的衛(wèi)星太陽翼由于推力器羽流場而受到熱影響的分布圖。因本例中分析的是受推力器羽流影響的太陽翼＋Z面，所以提取STL數(shù)據(jù)中n（Z）＝1的三角形網(wǎng)格作為有效網(wǎng)格。兩種方法的輸入?yún)?shù)一致，太陽翼網(wǎng)格化參數(shù)相當（文獻［10］采用面積為238．3mm2的長方形網(wǎng)格，本文采用面積為237．7mm2的三角形網(wǎng)格）。從計算結(jié)果來看，熱流密度分布結(jié)果一致性良好，文獻［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］方法并不能分析這種無法用解析公式表達的不規(guī)則曲面，只能將其簡化為標準拋物面進行分析，分析結(jié)果見圖6（b）。本例中分析的是受推力器羽流影響的通信天線反射面的下表面，因此提取STL 數(shù)據(jù)中n（Z）＜0的三角形網(wǎng)格作為有效網(wǎng)格。從分析結(jié)果來看，本文方法得到的最大熱流密度為1．242kW／m2，最大熱流密度作用區(qū)域為（-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的三角形區(qū)域；整個反射面下表面受到羽流影響的熱容量為174．8 W。文獻［10］方法得到的最大熱流密度為1．166kW／m2，整個反射面下表面受到羽流影響的熱容量為106．1 W。

比較圖6（a）和圖6（b）可見，本文方法和文獻［10］方法得到的熱流密度分布相同，只是由于圖6（a）所示模型中受到羽流熱影響較大區(qū)域的形狀不規(guī)則，比圖6（b）所示的標準拋物面更靠近推力器，因此用前者得出的最大熱流密度比后者的要大，更符合實際情況。

圖6 通信天線反射面上羽流熱影響分布Fig．6 Plume thermal effect distribution on communication antenna’s reflector

4．2 與Plume軟件比較

圖7是分別采用本文方法和Plume軟件得到的某推力器羽流對天線展開臂轉(zhuǎn)接頭部位的熱影響分布。該轉(zhuǎn)接頭結(jié)構(gòu)較復(fù)雜，文獻［10］方法無法建立羽流分析模型。由于Plume軟件難以對復(fù)雜的模型建模，因此僅建立轉(zhuǎn)接頭下半部分（熱流密度最大部位）的模型。圖7（a）和圖7（b）顯示出熱流密度分布相同，從熱流密度最大值來看，本文方法所得結(jié)果（168kW／m2）和Plume軟件得到的結(jié)果（172kW／m2）僅相差2．4%。這說明該天線展開臂轉(zhuǎn)接頭部位局部熱流密度過大，須采取相應(yīng)熱控措施，或規(guī)劃推力器工作時機來加以解決。

圖7 天線展開臂轉(zhuǎn)接頭上羽流熱影響分布Fig．7 Plume thermal effect distribution on antenna’s deployed arm connector

5 結(jié)論

本文以文獻［9］推導(dǎo)出的羽流熱流場分布函數(shù)為基礎(chǔ)，提出一種基于STL數(shù)據(jù)的羽流熱影響分析方法，主要利用對推力器試驗數(shù)據(jù)的擬合實現(xiàn)，計算量較小。通過與文獻［10］方法以及Plume軟件的對比，可以得出如下結(jié)論。

（1）本文方法與文獻［10］方法和Plume軟件相比，在分析同樣設(shè)備受到的羽流熱影響時，具有良好的一致性。

（2）本文方法除了能分析形狀規(guī)則設(shè)備或能用解析公式表達設(shè)備受到的羽流熱影響，還能分析形狀不規(guī)則設(shè)備或不能用解析公式表達設(shè)備受到的羽流熱影響。這比文獻［10］和Plume軟件具有一定的優(yōu)勢。

本文所介紹的方法，包括但不限于雙組元10N推力器的羽流熱影響分析。只要提供類似文中式（1）的羽流熱流場分布函數(shù)，即可分析其他類型推力器的羽流對星上設(shè)備的熱影響。

（References）

［1］孫威，閆少光，田東波．衛(wèi)星發(fā)動機羽流污染模型的DSMC方法研究［J］．航天器環(huán)境工程，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］張健，尚志．航天器發(fā)動機羽流對敏感器熱效應(yīng)仿真研究［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軟件應(yīng)用于探月著陸器羽流效應(yīng)數(shù)值模擬研究［J］．航天器環(huán)境工程，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］程曉麗，王強，閻喜勤．衛(wèi)星姿控發(fā)動機高空羽流場工程分析［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——航天器空間任務(wù)及在軌分析工具［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］王平陽，鹿洪芳，程惠爾．雙組元姿控發(fā)動機反流特性的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，F(xiàn)ox 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 推力器羽流試驗數(shù)據(jù)建立羽流場數(shù)學模型［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