張 熇，吳成賡，劉立輝，孫 潔，賀碧蛟，蔡國飆

(1. 北京航空航天大學宇航學院，北京 102206；2. 北京空間飛行器總體設計部，北京 100094；3. 探月與航天工程中心，北京 100086)

0 引 言

月球作為距離地球最近的自然天體，因其獨特的空間位置以及豐富的科學信息，具有極高的探測價值，實現(xiàn)月球表面軟著陸是探月工程的重要目標之一。月球探測器月面著陸過程需要變推力發(fā)動機進行減速，發(fā)動機工作產生的高溫高壓燃氣進入真空環(huán)境迅速膨脹擴散形成羽流。著陸過程中探測器與月面距離不斷減小，發(fā)動機羽流與月面相互作用，不僅直接沖擊著陸緩沖機構與探測器底板，還會引發(fā)復雜的激波-邊界層干擾以及渦現(xiàn)象，惡化探測器下方熱環(huán)境，產生極大的羽流氣動熱效應，破壞探測器的正常結構，嚴重時甚至會導致著陸任務失敗。

為確保月球探測器月面安全穩(wěn)定著陸，必須精準評估發(fā)動機羽流的熱效應，分析不同狀態(tài)下羽流對不同部位的熱影響，以此為探測器著陸過程的關機策略制定提供支撐，為探測器關鍵部位的熱防護設計提供輸入。月球探測器發(fā)動機羽流場主要分為連續(xù)流區(qū)與稀薄流區(qū)兩部分，具體包括連續(xù)介質流、過渡領域流和自由分子流三種流態(tài)，涉及到連續(xù)流與稀薄流耦合，研究難度大。自20世紀50年代以來，計算流體力學的發(fā)展為連續(xù)流求解提供了較為完善的方法。在稀薄氣體動力學研究方面，國外基于直接模擬蒙特卡洛(DSMC)方法開發(fā)了一批通用的計算軟件，并針對探測器月面著陸過程的羽流效應分析進行了諸多應用，主要包括：美國約翰遜航天中心的DAC(DSMC Analysis Code)軟件、美國康奈爾大學的MONACO計算軟件、俄羅斯理論與應用機械研究所的SMILE(Statistical Modeling In Low-density Environment)軟件等，國內北京航空航天大學、上海交通大學、中國空氣動力研究與發(fā)展中心、中國航天空氣動力技術研究院等單位基于DSMC方法針對羽流的氣動作用特性也開展了一系列數(shù)值模擬研究。

月面著陸過程是一個復雜的高度非定常瞬態(tài)問題，需要量化評估探測器在距離月面不同高度下發(fā)動機羽流的熱效應，分析月面坡度的存在對羽流熱效應的附加影響。在探測器的觸月緩沖過程中，因為慣性發(fā)動機出口可能與月面逼近至更小距離而惡化底部熱環(huán)境，同時發(fā)動機若出現(xiàn)延時關機或拖尾效應，所產生的羽流熱影響也需要加以研判。因而有必要開展探測器高度、月面坡度和發(fā)動機狀態(tài)對羽流熱效應的影響研究。本文基于差分求解N-S方程與DSMC耦合的方法，數(shù)值模擬不同條件下的發(fā)動機羽流場，獲取探測器表面受羽流氣動作用后的熱流密度分布，研究探測器高度、月面坡度和發(fā)動機狀態(tài)對熱流密度分布的影響規(guī)律，最終評估月球探測器著陸過程的羽流熱效應。

1 數(shù)值模擬方法

流場的稀薄程度可由宏觀Knudsen () 數(shù)描述，=/，其中是平均分子自由程，為流場特征長度。當< 0.01時，流場可視為連續(xù)流，數(shù)越大則流場稀薄程度越高，基于連續(xù)介質假說的N-S方程將失效。月球探測器發(fā)動機羽流同時包含連續(xù)流與稀薄流的跨流域流動，且還需考慮羽流與月面的相互作用狀態(tài)，因此目前最為高效與實用的數(shù)值模擬方法是根據(jù)流場稀薄程度不同，對流場進行合理分區(qū)以及正確設置分界面，進而基于不同適配的數(shù)學物理模型分別求解。本文通過差分求解N-S方程計算探測器與月面特定相對高度、角度狀態(tài)下的發(fā)動機噴管內流場與近月面流場，在N-S方程求解結果的基礎上使用DSMC方法對外圍羽流場完成進一步計算。

1.1 差分求解N-S方程

采用有限差分、有限體積方法對N-S方程進行求解，控制方程主要包括：

連續(xù)方程：

(1)

動量方程：

(2)

能量方程：

(3)

式中：為密度，為時間，為速度矢量，為壓強，為黏性應力張量，為總能量，為熱傳導系數(shù)。

湍流模型選用適應性較好的二方程模型中的SST-湍流模型，以有限體積方法離散微分方程組。時間推進采用LU-SGS隱式方法，每一單元采用局部時間步長?？臻g推進采用AUSMP迎風格式。

1.2 DSMC方法

DSMC方法直接從物理實際出發(fā)，利用少量的模擬分子代替真實流場內數(shù)目眾多的氣體分子模擬計算物理過程，經(jīng)統(tǒng)計平均獲得宏觀流動參數(shù)，達到求解稀薄氣體問題的目的。

基于DSMC方法針對稀薄流區(qū)進行求解時，因為流場密度、溫度相對較低的特點，可作如下假設：流場中分子的碰撞均為二體碰撞；僅考慮分子轉動內能，忽略分子振動內能和分子化學非平衡效應；氣體流動為定常流動；分子統(tǒng)一視為變徑硬球分子(VHS)。DSMC方法以跟蹤統(tǒng)計粒子、實現(xiàn)粒子相對時間的隨機演化為特征，具有繼承性強、可靠性高、通用性好等諸多優(yōu)點，是解決稀薄氣體流動問題的有效方法，目前已經(jīng)廣泛應用于工程實踐。

北京航空航天大學基于DSMC方法自主研發(fā)了羽流效應計算軟件PWS，軟件設計過程中依據(jù)模塊化原則引入面向對象的編程思想，創(chuàng)新性地提出了流場與邊界解耦的高精度網(wǎng)格處理方法、區(qū)分粒子碰撞和運動的可變時間步長方法、針對多組分氣體的可變粒子權重方法，有效解決了常規(guī)DSMC方法在復雜邊界描述困難、計算量大等方面的一系列難題，成功實現(xiàn)了PWS軟件在羽流效應數(shù)值模擬上的通用化應用。

1.3 數(shù)值模擬校驗

基于美國卡爾斯本大學巴法羅研究中心(CUBRC)在高能激波風洞中試車的實驗數(shù)據(jù),校驗PWS軟件的數(shù)值模擬精度。CUBRC實驗測量了稀薄超聲速來流對雙圓錐體的氣動力熱數(shù)值，其中來流工質為氮氣，來流壓強為2.23 Pa。

圖1為PWS數(shù)值模擬得到的雙圓錐體流場馬赫數(shù)和流線分布圖。結果顯示，超聲速來流在雙圓錐體前形成了多道斜激波，激波的交匯處(即雙圓錐面的交接處附近)流動較為復雜且形成了明顯渦流，氣動熱影響較為嚴重。

圖1 數(shù)值模擬校驗馬赫數(shù)及流線分布圖

PWS數(shù)值模擬得到的雙圓錐體表面熱流密度數(shù)值結果與CUBRC實驗結果如圖2所示。對比可得，數(shù)值模擬與實驗結果變化趨勢一致、大小基本相符，偏差在-26%～+10%范圍內，在流動復雜的渦流區(qū)域符合較好。因此，本文認為PWS軟件的羽流熱效應數(shù)值模擬結果具有較高的精度。

圖2 雙圓錐體表面熱流密度數(shù)值模擬結果與實驗數(shù)據(jù)

2 探測器羽流熱效應數(shù)值模擬

2.1 連續(xù)流區(qū)數(shù)值模擬

連續(xù)流區(qū)包括發(fā)動機噴管內流場及近月面流場，數(shù)值模擬結果主要用于截取粒子入口截面，以作為稀薄流區(qū)的粒子布入條件。本文以某型月球探測器為參考，設定探測器在著陸過程中主要采用變推力發(fā)動機的2500 N工作模式進行減速，發(fā)動機混合比為1.572，燃燒室總壓為0.267 MPa。通過熱力計算確定總溫為2836.5 K，獲得燃燒產物的組分質量分數(shù)列于表1。

表1 發(fā)動機2500 N工作模式主要燃燒產物質量分數(shù)

綜合以上條件作為連續(xù)流區(qū)數(shù)值模擬的壓強入口，發(fā)動機內壁面屬性設置為無滑移絕熱壁面，其余邊界設置為壓強大小為0的壓強出口，通過差分求解N-S方程得到包括發(fā)動機噴管內流場與噴流近月面流場的連續(xù)流區(qū)數(shù)值模擬結果，其中探測器相對月面無偏角工況采用二維軸對稱計算模型，有偏角工況采用三維計算模型，圖3為典型高度連續(xù)流區(qū)的求解結果示例。

圖3 典型高度連續(xù)流區(qū)求解結果馬赫數(shù)分布圖

在連續(xù)流區(qū)設置分界面截取入口條件時，為確保流場計算結果有效，同時避免受到下游流場干擾，分界面所在位置流場需滿足<0.01以及>1，并且盡量避開發(fā)動機噴管下方核心區(qū)的復雜波系。

對于發(fā)動機噴管出口與月面間距離較小工況，為避免月面反射羽流對入口條件產生影響，分界面選取為噴管出口向軸線方向延伸至月面所得的一個母線傾角為30°的圓臺側面(圖3中右側的深色虛線)。對于發(fā)動機噴管出口與月面間距離較大工況，出口位置基本不受月面反射羽流的影響，分界面可直接選取為噴管出口。

2.2 稀薄流區(qū)數(shù)值模擬

稀薄流區(qū)數(shù)值模擬采用PWS軟件三維計算完成。DSMC計算域以粒子入口截面為起始面，沿粒子主流擴散方向進行擴展設置。經(jīng)進一步添加探測器幾何模型和流場空間體網(wǎng)格，即可完成發(fā)動機羽流場的求解，最終獲得羽流與月面、探測器相互作用后探測器表面的熱流密度分布。

探測器幾何模型包括探測器主體、著陸緩沖機構與發(fā)動機，在仿真中作為幾何包絡邊界模擬探測器的真實表面。幾何模型表面主體采用結構化的直角面網(wǎng)格，物面邊界連接區(qū)域采用非結構化的三角面網(wǎng)格。面網(wǎng)格主要參與同模擬羽流分子的相互作用計算過程，并記錄存儲產生的熱流密度，以此分析著陸過程中探測器重點位置受羽流氣動熱效應后的具體影響。

流場空間計算體網(wǎng)格配置以涵蓋探測器主要關注部位以及控制計算量為原則，體網(wǎng)格類型選擇正交化直角網(wǎng)格，結構劃分以粒子分布密度為標準，在數(shù)密度大的區(qū)域進行合理加密。對于無偏角的對稱工況，計算域可結合探測器主體的對稱特性，單獨劃分探測器的1/4區(qū)域進行數(shù)值模擬。圖4為網(wǎng)格劃分及入口截面的示意圖，該工況下發(fā)動機噴管出口與月面間距離較小，在連續(xù)流區(qū)設置分界面截取得到的粒子入口條件即為圖中右側的深色表面。進行1/4對稱設置后的計算域大小為(1.53 m×2.47 m×2.47 m)，流場體網(wǎng)格尺寸與附近探測器面網(wǎng)格尺寸基本保持一致，在DSMC入口截面和月面附近進行了適當加密，最小網(wǎng)格尺度為5 mm。

圖4 數(shù)值模擬網(wǎng)格劃分及入口條件示意圖

2.3 典型數(shù)值模擬結果

根據(jù)上述數(shù)值模擬結果獲取著陸緩沖機構與探測器底板的熱流密度數(shù)值，進而分析探測器著陸過程羽流氣動熱效應的影響。

圖5為發(fā)動機出口中心距離月面0.434 m時，在發(fā)動機2500 N工作模式羽流影響下的探測器表面熱流密度分布圖。結果表明，著陸緩沖機構上受到的羽流氣動熱影響顯著嚴重于探測器底板，且不同部位熱流密度的分布特征區(qū)別明顯。在此條件下，著陸緩沖機構上從支柱到足墊的熱流密度總體呈逐漸減小的變化趨勢，探測器底板上從中心到邊緣的熱流密度總體呈逐漸增大的變化趨勢。

圖5 探測器表面氣動熱流密度分布

3 結果與討論

3.1 著陸高度影響分析

為評估著陸過程中,不同時刻探測器受發(fā)動機羽流熱效應的影響程度，設置8個工況數(shù)值模擬并對比著陸緩沖機構與探測器底板的熱流密度數(shù)值。本節(jié)8個工況發(fā)動機出口中心與月面間距離分別為3.08 m,2.68 m,2.28 m,1.88 m,1.48 m,1.08 m,0.757 m與0.434 m；發(fā)動機工作模式為2500 N；月面無坡度；月面溫度為300 K。

著陸緩沖機構主要分為支柱與足墊兩部分，如圖6所示。在支柱上重點關注關鍵點P,P,P,P,P，在足墊上設置關鍵點P,P,P。

圖6 著陸緩沖機構關鍵點分布示意圖

圖7為發(fā)動機以2500 N模式在不同高度工作時，著陸緩沖機構上各關鍵點的熱流密度隨發(fā)動機出口與月面間距離的變化關系。當發(fā)動機出口中心距離月面在2.28 m以上時，著陸緩沖機構上的熱流密度均在10 kW/m以下。當距離減小到1.48 m時，著陸緩沖機構上的熱流密度開始突破100 kW/m。

圖7 著陸緩沖機構關鍵點氣動熱流密度隨距離變化關系

在探測器下降過程中,隨著距離減小，支柱表面(P～P)的熱流密度近似呈指數(shù)增加，當發(fā)動機出口中心距離月面0.434 m時，支柱上熱流密度達到最大值，約為429 kW/m。在探測器下降過程中隨著距離減小至1.48 m，足墊表面(P～P)的熱流密度增加規(guī)律與支柱相似，但隨著距離進一步減小，呈現(xiàn)出減小的趨勢。當發(fā)動機出口中心距離月面在0.757 m～1.48 m時，足墊表面各點熱流密度相繼達到最大水平，其中最大值約為165 kW/m。

在探測器底板上設置關鍵點P,P,P,P,P,P,P，分布圖如圖8所示。

圖8 探測器底板關鍵點分布示意圖

圖9為發(fā)動機以2500 N模式在不同高度工作時，探測器底板上各關鍵點的熱流密度隨發(fā)動機出口與月面間距離的變化關系。當發(fā)動機出口中心距離月面在2.28 m以上時，探測器底板上的熱流密度均在1 kW/m以下。當距離減小到1.48 m時，底板上的熱流密度開始突破10 kW/m。

圖9 探測器底板關鍵點氣動熱流密度隨距離變化關系

在探測器下降過程中,隨著距離減小，探測器底板上各關鍵點的熱流密度變化趨勢基本一致，總體上近似呈指數(shù)增加，當發(fā)動機出口中心距離月面0.434 m時，熱流密度達到最大值，約為47 kW/m。在相同高度下，探測器底板邊緣位置關鍵點的熱流密度相比靠近中央位置的關鍵點更大。針對所取關鍵點按P～P,P～P,P～P三組分別進行對比，結果顯示越接近探測器底板邊緣的關鍵點熱流密度更高。

為探究隨著陸高度減小,探測器不同位置受發(fā)動機羽流熱效應影響程度的變化規(guī)律，需結合流場圖進一步分析，圖10為探測器在三個典型高度下的羽流場壓強分布圖。

圖10 不同著陸高度下探測器羽流場壓強分布圖

不同著陸高度下的流場結果顯示，在探測器下降過程中，隨著發(fā)動機出口與月面間距離不斷減小，羽流與月面作用反射產生的返流可逐漸擴散至著陸緩沖機構的支柱上，且作用強度不斷增強，使得支柱上各關鍵點熱流密度總體呈單調上升趨勢。隨著探測器高度不斷減小，羽流相對足墊表面的入射角度和作用強度不斷變化，當發(fā)動機出口距離月面在1.08 m左右時，羽流在足墊上的滯止作用最劇烈，而隨著陸高度進一步減小，足墊周圍羽流與月面更為接近，受到月面的冷卻作用更強，羽流的能量特性進一步降低。綜合以上原因，隨著陸高度減小，足墊上關鍵點的熱流密度變化呈現(xiàn)先增大后減小的趨勢。

不同著陸高度下的流場結果表明，在探測器下降過程中，隨著發(fā)動機出口與月面間距離不斷減小，探測器底板下方羽流場壓強逐漸增大，羽流作用強度不斷增強，因此底板上各關鍵點的熱流密度隨距離減小總體上也呈單調增大趨勢。由于探測器下方的中央位置靠近發(fā)動機為高壓區(qū)，邊緣位置因外部真空環(huán)境影響為低壓區(qū)，結合流場圖可以判斷，羽流與月面作用后的向上反射過程具有強烈的向外擴散特征，底板邊緣下方相對中央位置流場壓強更大，因此受到的氣動熱效應影響也更為嚴重，印證了高度相同時底板上靠外側關鍵點熱流密度相比內側更大的規(guī)律。

3.2 月面坡度影響分析

實際著陸過程中,落月位置可能存在坡度，為評估月面坡度對探測器著陸過程發(fā)動機羽流的熱效應影響，在發(fā)動機出口軸線與月面有無偏角的兩種狀態(tài)下共設置8個工況，分析著陸緩沖機構與探測器底板上羽流氣動熱流密度的最大值大小。偏角大小由著陸過程中可能出現(xiàn)的極限角度確定為12°，考慮到發(fā)動機出口距離月面2.28 m以上時探測器受羽流熱效應影響較小，因此只設置距離不大于2.28 m的工況，圖11為有無偏角狀態(tài)下著陸緩沖機構與探測器底板熱流密度最大值隨距離的變化關系。

圖11 探測器相對月面不同姿態(tài)著陸緩沖機構與底板熱流密度最大值隨距離變化關系

數(shù)值模擬結果顯示，當發(fā)動機出口中心與月面間距離相同時，著陸緩沖機構與探測器底板上的熱流密度最大值在有偏角狀態(tài)下均小于無偏角狀態(tài)。針對羽流的流動狀態(tài)進行分析，當落月位置存在坡度時，探測器底板與月面不平行，探測器下方一側區(qū)域空間更大，具有更強的導流作用，有利于羽流的擴散膨脹。存在坡度時，發(fā)動機羽流與月面的作用過程將由近似垂直正入射變?yōu)樾比肷?，有助于減弱羽流的壓縮減速與滯止作用，減小羽流導致的氣動熱效應，最終改善探測器周圍的流場熱環(huán)境。

本文研究結果表明，月球探測器著陸過程中月面坡度的存在對羽流氣動熱效應無加劇作用，落月位置無坡度狀態(tài)相對有坡度狀態(tài)，發(fā)動機羽流對探測器的氣動熱更惡劣，因此可采用無坡度狀態(tài)的氣動熱計算結果指導探測器的熱防護設計。

3.3 緩沖過程慣性及延時拖尾效應影響分析

月球探測器著陸過程因落月慣性、關機信號延時、發(fā)動機拖尾等因素影響可能存在特殊工況，在羽流熱效應數(shù)值模擬研究中，需要考慮的主要包括緩沖過程發(fā)動機出口與月面達到極限距離以及發(fā)動機關機的延時拖尾效應兩類情況。為探究上述因素的影響程度，共設置12個工況。其中，選取發(fā)動機2500 N模式距離月面0.234 m條件下的工作狀態(tài)研究落月慣性的影響，選取發(fā)動機1000 N, 500 N兩種模式分別距離月面1.08 m, 0.757 m, 0.434 m與0.234 m條件下的工作狀態(tài)研究發(fā)動機延時拖尾效應的影響。

發(fā)動機在2500 N, 1000 N與500 N工作模式下著陸緩沖機構與探測器底板熱流密度的最大值隨距離的變化曲線如圖12所示。

圖12 發(fā)動機不同工作模式探測器關鍵部位氣動熱流密度最大值隨距離變化關系

數(shù)值模擬結果表明，隨著發(fā)動機出口距離減小，著陸緩沖機構與探測器底板上熱流密度最大值增加。當距離由0.434 m減小到0.234 m時，著陸緩沖機構上熱流密度最大值急劇升高至888 kW/m，約為0.434 m狀態(tài)下的兩倍。探測器底板上熱流密度最大值增幅較小，由0.434 m狀態(tài)下的78 kW/m略微增長至81 kW/m。由此可見，探測器在觸月緩沖過程中，若因慣性導致發(fā)動機出口與月面間距離減小至極限狀態(tài)(0.234 m)后，著陸緩沖機構附近熱環(huán)境將進一步惡化，而探測器底板受到的氣動熱效應變化較小。

數(shù)值模擬結果顯示，當發(fā)動機處于1000 N工作狀態(tài)時，于0.434 m高度在探測器上產生的最大熱流密度為154 kW/m，于0.234 m高度產生的最大值為399 kW/m。相同高度下，發(fā)動機以1000 N工作模式在著陸緩沖機構與探測器底板上產生的熱流密度最大值基本處于2500 N模式下的40～60%范圍內。當發(fā)動機處于500 N工作狀態(tài)時，于0.434 m高度在探測器上產生的最大熱流密度為71 kW/m，于0.234 m高度產生的最大值為167 kW/m。相同高度下，發(fā)動機以500 N工作模式在著陸緩沖機構與探測器底板上產生的熱流密度最大值約為2500 N模式下的20%。

綜上所述，落月慣性導致探測器發(fā)動機出口與月面間距離逼近至0.234 m后，羽流與月面相互作用將進一步加劇探測器受到的氣動熱效應；發(fā)動機關機過程的延時及拖尾效應依然對探測器具有一定程度的熱影響，但產生的熱流密度最大值明顯小于以2500 N模式工作時的對應數(shù)值，且大致比例與推力大小近似相關。

4 結 論

本文針對月球探測器著陸過程羽流熱效應進行了數(shù)值模擬研究，通過差分求解N-S方程與DSMC耦合的方法完成了發(fā)動機噴管內流場、近月面流場與外圍羽流場的數(shù)值模擬，獲得了探測器表面不同條件下受羽流氣動作用后的熱流密度分布，研究結果表明：

1) 探測器受到的羽流氣動熱效應總體隨著陸高度的減小急劇增強，在著陸高度為3.08 m至0.434 m范圍內，著陸緩沖機構和探測器底板關鍵點的熱流密度最大值分別為429 kW/m和47 kW/m；當著陸高度大于2.28 m時，著陸緩沖機構和探測器底板受羽流氣動熱效應的影響較小，熱流密度均小于10 kW/m。

2) 當月面存在坡度時，因導流作用及羽流入射角度改變，探測器受羽流氣動熱效應影響相對同高度無坡度狀態(tài)更小。

3) 因落月慣性導致探測器發(fā)動機出口與月面間距離逼近至0.234 m后，著陸緩沖機構與探測器底板熱流密度最大值分別升至888 kW/m與81 kW/m；當發(fā)動機關機過程出現(xiàn)延時及拖尾效應時，探測器將受到一定程度熱影響，但熱流密度最大值明顯小于同條件2500 N模式工作狀態(tài)。