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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

在進(jìn)行羽流熱影響分析時(shí)，首先，從設(shè)備的三維模型提取STL數(shù)據(jù)；然后，根據(jù)分析的需要，將STL數(shù)據(jù)文件中的有效三角形網(wǎng)格提取出來，剔除掉不需要分析的三角形網(wǎng)格；最后，以所提取的三角形網(wǎng)格為對(duì)象，建立其所受羽流熱影響的分析模型。整個(gè)方法的流程見圖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)格的集合）。在這個(gè)過程中，是將整個(gè)模型的外表面都網(wǎng)格化，而不能只將某些特定曲面進(jìn)行網(wǎng)格化處理。推力器和設(shè)備的相對(duì)位置關(guān)系，決定了羽流只會(huì)對(duì)整個(gè)模型的某些特定曲面產(chǎn)生影響（不考慮熱傳導(dǎo)的影響）。因此，如何從整個(gè)STL數(shù)據(jù)文件中提取出受羽流影響曲面的三角形網(wǎng)格，成為羽流分析的前提。

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

圖4 三角形網(wǎng)格法向和衛(wèi)星本體坐標(biāo)系坐標(biāo)軸關(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)星本體坐標(biāo)系Z軸夾角小于90°）的三角形網(wǎng)格，n（Z）為三角形網(wǎng)格的法向矢量在衛(wèi)星本體坐標(biāo)系Z軸上的分量。

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

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

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

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

羽流場(chǎng)的分布角為

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

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

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

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

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

4 應(yīng)用

以3個(gè)應(yīng)用實(shí)例對(duì)本文方法與文獻(xiàn)［10］方法、Plume軟件作一比較。

4．1 與文獻(xiàn)［10］方法比較

圖5是分別采用本文方法和文獻(xiàn)［10］方法獲得的衛(wèi)星太陽(yáng)翼由于推力器羽流場(chǎng)而受到熱影響的分布圖。因本例中分析的是受推力器羽流影響的太陽(yáng)翼＋Z面，所以提取STL數(shù)據(jù)中n（Z）＝1的三角形網(wǎng)格作為有效網(wǎng)格。兩種方法的輸入?yún)?shù)一致，太陽(yáng)翼網(wǎng)格化參數(shù)相當(dāng)（文獻(xiàn)［10］采用面積為238．3mm2的長(zhǎng)方形網(wǎng)格，本文采用面積為237．7mm2的三角形網(wǎng)格）。從計(jì)算結(jié)果來看，熱流密度分布結(jié)果一致性良好，文獻(xiàn)［10］得到的最大熱流密度為16．72 W／m2，整體所受熱影響為199．5 W；本文方法得到最大熱流密度為17．2W／m2，整體所受熱影響為187．0 W。

圖5 太陽(yáng)翼上羽流熱影響分布Fig．5 Plume thermal effect distribution on solar array

圖6（a）是采用本文方法分析某賦形拋物面天線的反射面獲得的羽流熱影響分布圖。文獻(xiàn)［10］方法并不能分析這種無法用解析公式表達(dá)的不規(guī)則曲面，只能將其簡(jiǎn)化為標(biāo)準(zhǔn)拋物面進(jìn)行分析，分析結(jié)果見圖6（b）。本例中分析的是受推力器羽流影響的通信天線反射面的下表面，因此提取STL 數(shù)據(jù)中n（Z）＜0的三角形網(wǎng)格作為有效網(wǎng)格。從分析結(jié)果來看，本文方法得到的最大熱流密度為1．242kW／m2，最大熱流密度作用區(qū)域?yàn)椋?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ū)域；整個(gè)反射面下表面受到羽流影響的熱容量為174．8 W。文獻(xiàn)［10］方法得到的最大熱流密度為1．166kW／m2，整個(gè)反射面下表面受到羽流影響的熱容量為106．1 W。

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

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

4．2 與Plume軟件比較

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

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

5 結(jié)論

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

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

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

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

（References）

［1］孫威，閆少光，田東波．衛(wèi)星發(fā)動(dòng)機(jī)羽流污染模型的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ā)動(dòng)機(jī)羽流對(duì)敏感器熱效應(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］蔡國(guó)飆，賀碧蛟．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］程曉麗，王強(qiáng)，閻喜勤．衛(wèi)星姿控發(fā)動(dòng)機(jī)高空羽流場(chǎng)工程分析［J］．空間科學(xué)學(xué)報(bào)，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］恒潤(rùn)科技公司．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］王平陽(yáng)，鹿洪芳，程惠爾．雙組元姿控發(fā)動(dòng)機(jī)反流特性的DSMC模擬［J］．工程熱物理學(xué)報(bào)，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 推力器羽流試驗(yàn)數(shù)據(jù)建立羽流場(chǎng)數(shù)學(xué)模型［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