火星環繞器羽流效應仿真研究

陳 杰，賀碧蛟，蔡國飆

（北京航空航天大學宇航學院，北京 100191）

火星環繞器羽流效應仿真研究

陳 杰，賀碧蛟?，蔡國飆

（北京航空航天大學宇航學院，北京 100191）

火星探測器發動機羽流效應影響飛行任務的安全，采用差分求解NS方程與直接模擬蒙特卡洛耦合的方法對火星環繞器姿軌控發動機羽流力、熱以及污染效應進行了仿真研究。結果表明：發動機羽流對火星環繞器氣動力、熱及污染效應影響主要在發動機儲箱上，最大熱流值可以達到1789 W／m2，需要注意；發動機羽流對展開太陽能電池翼的氣動力和氣動力矩相對于發動機推力和推力力矩來說，小3～4個量級，影響較??；已分離著陸器在距發動機出口300 m以內時，羽流對已分離著陸器的影響較大，氣動熱流在100 W／m2以上；在100 m以內時，氣動熱流在1000 W／m2以上，需要對敏感設備采取一定的防護措施。不同熱適應系數氣動熱流仿真對比計算結果表明，隨著熱適應系數的增加，熱流密度值呈線性增加。數值模擬結果可以為火星環繞器工程設計提供參考。

火星探測；火星環繞器；羽流效應；DSMC

1 引言

2020年我國將進行首次火星探測，火星環繞器作為火星探測系統的一部分，將要完成地火轉移飛行和軌道修正、完成火星捕獲、進入火星停泊軌道、與巡視器可靠分離、攜帶載荷完成環火遙感等科學探測任務。由于發動機羽流會對火星環繞器及其上敏感設備、太陽能帆板等產生氣動力和力矩、氣動熱以及污染效應，對火星環繞器的正常工作造成影響。為了使火星環繞器能夠順利完成任務，需要對火星環繞器發動機羽流效應進行仿真研究。

直接模擬蒙特卡洛（DSMC）方法［1］是目前模擬發動機羽流效應較為成熟的方法?；贒SMC方法，北京航空航天大學的蔡國飆等人開發了PWS（Plume Work Station）軟件［2-3］，驗證算例表明軟件具有較高的置信度［4］。PWS軟件已成功應用在神舟飛船、嫦娥3號著陸器等航天器的羽流仿真研究中［5-6］。

本文采用差分求解NS（Navier-Stokes）方程與DSMC相耦合的方法［7］，使用PWS對不同工況下發動機羽流對環繞器表面及敏感元件的氣動力、熱和污染效應進行仿真分析，對比不同熱適應系數對環繞器氣動熱流仿真結果的影響。對于火星環繞器與火星著陸巡視器分離的過程，由于已釋放的著陸巡視器與環繞器的距離不斷改變，著陸巡視器與羽流作用是一個復雜的非穩態過程。本文采用將羽流流場與已釋放著陸巡視器解耦的辦法對羽流對已釋放著陸巡視器的影響進行估算。

2 仿真模型

PWS軟件采用二級構造法來構建復雜仿真模型［2］。首先在不影響羽流效應仿真分析結論的前提下，對實際模型進行適當簡化。然后對PWS軟件定義的幾種常用的平面（直線）和曲面（曲線）等基本元素，使用“與”和“或”運算最終組合成擁有復雜邊界的仿真模型。圖1為使用PWS軟件建立的火星環繞器羽流仿真簡化模型。

3 數值方法

考慮到發動機內流場以及羽流核心區流場具有較高的壓強和密度，使用DSMC方法求解，計算量過大。本文采用差分求解NS方程與DSMC相耦合的方法，對火星環繞器發動機羽流效應進行仿真研究。

對于火星環繞器與火星著陸巡視器分離的過程，本文采用解耦的辦法對羽流對已釋放著陸巡視器的影響進行估算。即首先計算發動機羽流場，得到羽流場壓強、密度、溫度、速度分布，然后根據著陸巡視器與環繞器距離和與發動機推力夾角，確定火星探測器分離軌跡，并得到軌跡上羽流壓強、密度、溫度和速度參數，最后根據局地平衡假設及分子動力學理論，計算分離軌跡上的羽流質量流率、壓強和熱流密度，估算羽流效應。

1）差分求解NS方程

目前計算連續流域流動的最常用方法是數值求解NS方程的方法。矢量形式的NS方程如式（1）所示［8］，其中，Q為守恒變量，F、G為無粘通量，Fv、Gv為粘性通量，S為源項。

一般真實的流場都是以湍流狀態存在的，因此，控制方程中的粘性項和熱傳導項中的系數需由層流和湍流共同確定［9］，如式（2）、式（3）：

式中，μ、μl和μt分別為粘性系數、層流和湍流粘性系數，Pr、Prl和Prt分別為普朗特數、層流和湍流普朗特數。層流粘性系數一般隨溫度而變化，可由Sutherland公式較為精確地給出［9］。在求解雷諾平均控制方程時，為了使其封閉，必須引入計算湍流粘性系數的湍流模型［8］，計算中采用SST k-ω兩方程模型。

2）DSMC方法

由于DSMC計算域的流場密度和氣體溫度較低，計算時做以下假設［10］：（1）假設流場中分子間的碰撞為二體碰撞；（2）僅考慮分子的轉動內能；（3）氣體流動為定常流動；（4）假設分子為可變徑硬球模型（VHS）。

對于分子與壁面的作用模型，目前應用較為廣泛的是Maxwell反射模型。Maxwell反射模型采用熱適應系數α來表示分子在與壁面作用的過程中，鏡面反射與漫反射各占的比例［10］，計算公式如式（4）：

其中，qi和qr是來流和反射流的能流量，qw是反射溫度為壁溫時完全Maxwell漫反射下的能流量。在完全Maxwell漫反射時，α＝1，完全鏡面反射時，α＝0。熱適應系數表征了反射分子的溫度在多大程度上已經“適應”了壁面的溫度。

4 羽流效應計算

火星環繞器在不同階段將執行不同的飛行任務，需要不同數量、不同推力的發動機組合工作。對每個發動機工況，首先使用差分求解NS方程的方法得到發動機內流場及部分外流場，再使用得到的發動機出口截面參數作為發動機二維自由羽流場DSMC計算的入口條件，然后從得到的二維流場計算結果中截取流場參數作為三維羽流場入口條件，對三維羽流場進行仿真，得到不同工況下發動機羽流對火星環繞器及敏感元件的氣動力、熱及污染效應。采用不同熱適應系數對發動機儲箱的氣動熱流密度進行了對比計算。

1）入口條件

對火星環繞器進行全流場計算時，計算域較大，發動機推力較大以及發動機數量較多。在不影響計算結果的情況下，為減小計算量，將發動機噴管向外延伸一個母線傾角為15°的圓錐，以此圓錐面上的流場參數作為三維DSMC計算入口條件。

圖2為分別以發動機出口和發動機噴管向外延伸一個母線傾角為15°的圓錐為DSMC入口條件得到的羽流場壓力對比云圖。可以看出，兩種入口條件得到的羽流場基本吻合，只有在離發動機較遠處、壓力較小的流場區域有些許偏差，且偏差在工程允許范圍內。結果說明，將發動機噴管向外延伸一個母線傾角為15°的圓錐作為DSMC入口條件，對羽流流場計算結果影響不大，此做法較為合理。對120 N和25 N發動機羽流場分析可以得到類似的結論。

為了得到此圓錐面上的流場參數，首先差分求解NS方程，得到發動機內流場及部分外流場結果，然后以發動機出口截面參數作為入口條件用DSMC方法對二維發動機羽流場進行計算，最后在得到的二維羽流場中截取圓錐面母線上的流場參數，作為三維羽流場DSMC計算入口條件。

圖3為從3000 N發動機出口開始的5 m×5 m的外流場壓力分布云圖，圖中紅線表示作為三維DSMC計算入口的圓錐面位置。圓錐面所在位置的粒子壓強和密度與發動機出口處相比，已大大降低，達到了減小計算量的目的。

2）邊界條件

計算域邊界：當粒子經計算域邊界運動出計算域時，按逃逸處理，即將粒子刪除。

模型壁面：當粒子撞擊到模型邊界時，將粒子運動按Maxwell反射處理。模型壁面的溫度設置為常溫300 K（模型壁面溫度的變化，對羽流仿真結果影響不大）［2］，發動機壁面溫度設置為1000 K，熱適應系數設定為1。在對比算例中，分別計算了熱適應系數為0、0.25、0.5、0.75和1的對比算例。

入口條件邊界：將入口條件圓錐面設置為吸收邊界條件，即在入口條件邊界內部不存在粒子，粒子運動穿過吸收邊界時即被刪除。

5 結果分析

5.1 差分求解NS方程仿真結果

首先對發動機內流場及部分外流場使用差分求解NS方程的方法進行仿真計算。

圖4為3000 N發動機內流場及部分外流場流場溫度分布云圖?？梢钥闯觯瑴囟仍谌紵覂茸罡撸瑖姽軘U張段為特形面，氣體充分擴張。高溫高速噴流噴入真空環境中，羽流邊界形成了拋物線形的膨脹波，噴管型面作用使噴流形成兩個連在一起的錐形激波。3000 N發動機流場壓強、密度、Ma數流場分布云圖也呈現出類似的規律。

5.2 羽流氣動力、熱以及污染效應的影響

在差分求解NS結果的基礎上，使用DSMC方法對發動機二維羽流場和三維羽流場進行仿真計算。對所有工況仿真結果分析總結，發現在3000 N+4×120 N+2×25 N發動機組合工作時，火星環繞器及其上設備所受到的羽流效應影響最大。

圖5為該工況下的三維流場切面溫度分布圖?？梢钥闯觯诎l動機儲箱上方流場溫度最高，這是由于3000 N與120 N發動機噴流相互阻滯的緣故，這與發動機儲箱上的熱流密度最大的仿真結果相一致。流場密度、壓強分布圖等也呈現出相似的規律。

圖6為該工況下環繞器底板及其上設備羽流影響分布云圖。可以看出，在該工況下，發動機羽流對發動機儲箱的羽流效應影響最大，儲箱上的最大熱流值可以達到1789 W／m2，CO2質量到達率可以達到9.12×10-5kg／m2／s。對其它結果云圖分析發現，羽流氣動力壓強、H2O質量到達率等分布呈現出類似的規律，各最大值均在發動機儲箱上；而對于離3000 N主發動機距離相對較遠的火星環繞器上的其它敏感元器件來說，其附近的羽流流場密度壓強、溫度、密度等變得相對較小，對敏感元器件的羽流氣動力、熱以及污染效應也因此變得較小。因此，在該火星環繞器設計方案下，主要考慮發動機羽流對發動機儲箱所造成的氣動力、熱及污染效應。

5.3 太陽能電池翼上羽流氣動力和力矩的影響

由于太陽能電池翼展開半徑較大，需要考慮各工況下展開太陽能電池翼上的羽流氣動力以及其對探測器中心點坐標（0，0，0.705）的力矩值。圖7為0°、45°和90°的太陽翼所受到的發動機羽流氣動力壓強分布云圖。

對不同角度的太陽能電池翼上的力和力矩的計算結果總結分析，得到不同角度（°）下的太陽能電池翼所受到的氣動力在坐標軸上的分力值如表1所示。

表1 不同角度太陽能電池翼上羽流氣動力Table 1 Aerodynamic force on solar wings with differrent angles

表2為工作的兩臺發動機A1和A2的各自推力在坐標軸上的分力值，以及兩臺發動機推力合力在坐標軸上的推力值。表3為不同角度下的太陽能電池翼所受到的羽流氣動力相對于參考點坐標的力矩值。表4為A1、A2發動機推力相對于參考點坐標的力矩值以及它們的和推力相對于參考坐標的力矩值。

表2 發動機推力Table 2 Thrust of motors

表3 太陽能電池翼上羽流氣動力相對于參考點的力矩Table 3 Torque of aerodynamic force on solar wings relative to the reference point

表4 發動機推力相對于參考點的力矩Table 4 Torque of motor thrust relative to the reference point

由表1和表2可以看出，發動機羽流對太陽能電池翼的羽流氣動力在0.001 N量級，而發動機的推力則在10 N量級，即與發動機推力相比，太陽能電池翼上的羽流氣動力影響可以忽略。

由表3和表4可以看出，發動機羽流對太陽能電池翼的羽流氣動力相對于參考點（0，0，0.705）的力矩在0.01 N·m量級，而發動機推力相對于參考點的力矩值在10 N·m量級，與發動機推力矩相比，羽流氣動力所造成的力矩影響也較小。

圖8為太陽能電池翼上的氣動力隨著角度變化的變化趨勢圖。圖9為太陽能電池翼上的氣動力相對于參考點（0，0，0.705）的氣動力矩隨著角度變化的變化趨勢圖。

由圖8可以看出，45°太陽能電池翼所受到的x軸向的氣動力最大，即45°太陽翼受到的繞火星環繞器中心軸的扭轉力最大。隨著太陽能電池翼偏轉角度的增加，太陽能電池翼上所受到的沿z軸的氣動力逐漸增大，即90°時最大，0°時最??；沿y軸的氣動力先減小后增大，在45°時最小，90°時最大。由圖9可以看出，沿x軸方向的羽流氣動力矩隨著太陽能電池翼轉角的增大而增大，這是因為隨著轉角的增大，太陽能電池翼受到的沿x軸方向的力越來越大，而力臂基本相同；沿y軸方向的羽流氣動力矩在各角度下都基本為0；沿z軸方向的力矩隨著太陽能電池翼轉角的增大先減小后增大，在轉角為45°時絕對值最大，方向與z軸正向相反；合力矩隨著太陽能電池翼轉角的增大而增大。

5.4 發動機羽流對已分離著陸器的影響

采用解耦的方法估算不同推力發動機羽流對距離發動機出口不同距離的已分離火星著陸器的影響。圖10為距發動機出口不同距離下，3000 N發動機羽流場氣動熱流分布對比圖?？梢钥闯?，在與發動機軸線夾角相同時，熱流密度隨著與發動機出口距離的增加而減小；而在距離發動機出口相同距離時，隨著與發動機軸線夾角度數的增加，熱流密度也呈下降趨勢。仿真結果符合羽流場結構特點，即隨著流場位置與發動機軸線夾角的不斷增大，流場逐漸由核心流域過渡到稀薄流域，密度、壓強、溫度等流場參數值不斷降低；而在與發動機軸線同一夾角下的位置上，噴流逐漸向真空膨脹，密度、壓強、溫度等流場參數值隨著噴流的膨脹而降低。對其它結果云圖分析發現，羽流場質量通量分布以及正應力分布也呈現出類似的規律。表5為3000 N發動機工作時，在已分離著陸器距發動機出口不同距離下的羽流對著陸器表面的最大質量通量、最大正應力和最大熱流密度值?？梢钥闯?，當已分離著陸器距發動機出口300 m時，質量通量達到10-5kg／m2量級，正應力達到0.01 Pa量級，熱流達到100 W／m2量級；當已分離著陸器距發動機出口100 m時，質量通量達到10-4kg／m2量級，正應力達到0.1 Pa量級，熱流達到1000 W／m2量級，即對于分離著陸器來說，發動機羽流對已分離著陸器的影響在300 m內應該考慮，在100 m內應該對相關敏感元器件采取一定的防護措施。

表5 3000 N發動機各計算結果的最大值Table 5 Maximum value of plume effects of 3000 N motor

5.5 熱適應系數對羽流氣動熱流仿真結果的影響

圖11為不同熱適應系數下發動機儲箱上最大熱流密度值的變化趨勢圖?？梢钥闯觯S著熱適應系數的增加，發動機儲箱上最大熱流密度值基本呈線性增加。

當熱適應系數為0時，發動機儲箱上的熱流密度為0，即此時粒子與發動機儲箱表面沒有發生能量交換，仿真結果符合理論實際。當熱適應系數由0.25增大到1時，發動機儲箱上熱流密度分布基本相同，只是熱流密度值隨著熱適應系數的增大而增大。發動機儲箱上的熱流密度最大值由熱適應系數α＝0.25時的464 W／m2，增大到熱適應系數α＝1時的1789 W／m2。

6 結論

對火星環繞器在入軌及器器分離過程中的羽流效應進行了仿真分析；對同一發動機工況下不同角度太陽能電池翼羽流氣動熱流進行了對比計算；對熱適應系數對羽流氣動熱流仿真結果的影響進行了研究，結果表明：

1）受羽流氣動力、熱及污染影響較大的為發動機儲箱，氣動熱流值可以達到1789 W／m2；羽流對環繞器表面及其他敏感元件的氣動力、熱及污染影響相對較小，氣動熱流值基本都在100 W／m2量級和10 W／m2量級。

2）發動機羽流對太陽能電池翼的氣動力以及氣動力相對于參考點的力矩數值較小，氣動力在0.001 N量級；相對于參考點的氣動力矩值在0.01 N·m量級，比發動機推力和發動機推力相對于參考點的力矩小3～4個量級；太陽能電池翼上羽流氣動力和力矩相對于發動機合力和合力矩可以忽略。

3）對于分離著陸器來說，發動機羽流對距發動機出口100 m的已分離著陸器的質量通量達到10-4kg／m2量級，正應力達到0.1 Pa量級，熱通量達到1000 W／m2量級；對距發動機出口300 m的已分離著陸器質量通量達到10-5kg／m2量級，正應力達到0.01 Pa量級，熱通量達到100 W／m2量級，即對于分離著陸器來說，發動機羽流對已分離著陸器的影響在300 m內應該考慮，在100 m內應該對相關敏感元器件采取一定的防護措施。

4）當太陽能電池翼處于不同角度時，45°太陽能電池翼所受到的x軸向的氣動力最大。隨著太陽能電池翼偏轉角度的增加，太陽能電池翼上所受到的沿z軸的氣動力逐漸增大；沿y軸的氣動力先減小后增大，在45°時最小。對于相對于參考點的氣動力矩來說，沿x軸方向的羽流氣動力矩隨著太陽能電池翼轉角的增大而增大；沿y軸方向的羽流氣動力矩在各角度下都基本為0；沿z軸方向的力矩隨著太陽能電池翼轉角的增大先減小后增大，在轉角為45°時絕對值最大，方向與z軸正向相反；合力矩隨著太陽能電池翼轉角的增大而增大。

5）熱適應系數對羽流氣動熱流值的仿真結果的影響是線性的，即隨著熱適應系數的增大，羽流氣動熱流值仿真結果呈線性增加。

（References）

［1］ Bird G.Molecular Gas Dynamics［M］.Oxford：Clarendon Press，1976.

［2］ 蔡國飆，賀碧蛟.PWS軟件應用于探月著陸器羽流效應數值模擬研究［J］.航天器環境工程，2010，27（01）：18-23.Cai G B，He B J.PWS software applied to numerical simulation on plume effect of lunar exploration lander［J］.Spacecraft Environment Engineering，2010，27（01）：18-23.（in Chinese）

［3］ 侯鳳龍，蔡國飆.應用于DSMC方法的直角網格技術［J］.北京航空航天大學學報，2010，36（06）：694-699.Hou F L，Cai G B.Cartesian grid techmology applied in DSMC method［J］.Journal of Beijing University of Aeronautics and Astronautics，2010，36（06）：694-699.（in Chinese）

［4］ He B，He X，Zhang M，et al.Plume aerodynamic effects of cushion engine in lunar landing［J］.Chinese Journal of Aeronautics，2013，26（2）：269-278.

［5］ 蔡國飆，賀碧蛟.神舟飛船平移發動機對太陽能帆板羽流效應數值模擬研究［J］.載人航天，2011，17（04）：48-53.Cai G B，He B J.DSMC simulation for thee plume effect of the shenzhou spacecraft attitude control thrusters on the solar board［J］.Manned Spaceflight，2011，17（04）：48-53.（in Chinese）

［6］ 張熇，蔡國飆，許映喬，等.嫦娥三號著陸器軟著陸過程中羽流仿真分析及試驗研究［J］.中國科學：技術科學，2014，44（04）：344-352.Zhang H，Cai G B，Xu Y Q，et al.Simulation and experimental study of the plume during the Chang’E-3 lunar landing［J］.Scientia Sinica Technologica，2014，44（04）：344-352.（in Chinese）

［7］ 唐振宇，賀碧蛟，蔡國飆.解耦N-S／DSMC方法計算推力器真空羽流的邊界條件研究［J］.推進技術，2014，35（07）：897-904.Tang Z Y，He B J，Cai G B.Investigation on boundary conditions in decoupled N-S／DSMC method for vacuum plume simulation of thrusters［J］.Journal of Propulsion Technology，2014，35（07）：897-904.（in Chinese）

［8］ 約翰D.安德森.計算流體力學基礎及應用［M］.吳頌平，劉趙淼，譯.北京：機械工業出版社，2015.JohnD.Anderson.Foundation and Application of Computational Fluid Dynamics［M］.Wu S P，Liu Z M，translate.Beijing：China Machine Press，2015.（in Chinese）

［9］ 李斌.超燃沖壓發動機進氣道的數值模擬研究［D］.湘潭：湘潭大學，2013.Li B.Numerical Simulation of Scramjet Inlet／Isolator［D］.Xiangtan：Xiangtan University，2013.（in Chinese）

［10］ 沈青.稀薄氣體動力學［M］.北京：國防工業出版社，2003.Shen Q.Rarefied Gas Dynamics［M］.Beijing：National Defense Industry Press，2003.（in Chinese）

Simulation Study of Plume Effect on Mars Probe

CHEN Jie，HE Bijiao?，CAI Guobiao
（School of Astronautics，Beihang University，Beijing 100191，China）

To mitigate the threat of the motor plume on the flight mission of the Mars probe，the hybrid NS equations／DSMC method was used to simulate the plume effect on the Mars probe.The simulation results showed that the plume of the motor mainly affected the storage box and the maximum heat flow reached 1789 W／m2.The effects of the aerodynamic force and the torque relative to the reference point on the expended solar wing were 3～4 orders’smaller than that of the motor.The plume effect on the separated Mars lander within 300 m was large，and the value of the aerodynamic heat flux was more than 100 W／m2.Within 100 m，the value of the aerodynamic heat flux was more than 1000 W／m2.So some protective measures for the sensitive equipment was necessary.The simulation with different thermal accommodation coefficients showed that with the increase of the thermal accommodation coefficient，the heat flow value increased linearly.The numerical simulation results can provide references for the engineering design of the Mars probs.

Mars exploration；Mars orbiters；plume effect；DSMC

V476.4

A

1674-5825（2017）06-0743-08

2017-04-17；

2017-10-02

陳杰，男，碩士研究生，研究方向為姿軌控發動機羽流效應仿真。E-mail：fyaycj＠163.com

?通訊作者：賀碧蛟，男，博士，講師，研究方向為航天器真空羽流效應。E-mail：hbj＠buaa.edu.cn

（責任編輯：龐迎春）