地外天體起飛羽流導流氣動力效應仿真

蘇楊，蔡國飆，*，舒燕，葉青，張明星，賀碧蛟

（1.北京航空航天大學 宇航學院，北京100083； 2.北京空間飛行器總體設計部，北京100094；3.北京航天長征飛行器研究所，北京100076）

探測器完成地外天體表面任務后，起飛器負責攜帶相關設備脫離地外天體。起飛器在起飛過程中，發(fā)動機噴出的羽流與起飛平臺表面作用后反向流動至起飛器表面，會對起飛器產生明顯的氣動力效應。這一羽流導流氣動力效應會對起飛過程中起飛器產生力矩作用，對其姿態(tài)控制和保持產生影響，可能導致起飛器無法正常工作或者工作質量下降。因此，對地外天體起飛過程中起飛器及起飛平臺可能受到的羽流導流氣動力效應進行研究十分必要[1-2]。

目前成功實現月面起飛的國家只有美國和俄羅斯，由于相關數據的保密性，僅有少數文獻對Apollo登月艙下降級的凹碗型導流裝置的導流效果進行研究[3]。國外對羽流導流效果研究方面多關注于大氣環(huán)境下火箭起飛過程，關注點主要集中于導流機構受到的熱效應及聲載荷影響等[4-8]。中國自2013年成功實現嫦娥三號探測器月面軟著陸及月面巡視以來，相繼開展了地外天體起飛羽流導流氣動力熱效應的實驗研究和仿真研究。賀衛(wèi)東等[9-12]在高超聲速低密度風洞中，利用氮氣為工質，設計了1 ∶10縮比模型的羽流導流實驗，并對平板、凹槽和錐面的羽流導流機構的導流效果進行了實驗和仿真驗證。張明星等[13]在真空羽流實驗系統中模擬了真空環(huán)境，針對起飛器與錐形導流機構在不同距離、不同偏轉角度工況下起飛器受到的羽流導流氣動力效應開展縮比實驗研究，并與仿真結果進行了對比。張萃等[14]對導流機構受到羽流沖擊后的熱載荷進行了分析。雖然研究人員對地外天體起飛羽流導流效應進行了相關的縮比實驗研究和仿真研究，但是對大尺寸起飛器所受到的羽流導流氣動力效應研究并不充分。此外，由于受到著陸過程振動或天體表面不平整影響，起飛器在起飛平臺表面可能存在初始偏角，或在起飛過程中需要進行姿態(tài)調整。因此，起飛器起飛過程中，隨著上升距離和偏轉角度的改變，需對引起的羽流導流氣動力效應變化規(guī)律進行深入研究。

本文利用計算流體力學/直接模擬蒙特卡羅（CFD/DSMC）耦合方法[15-17]，對圓錐導流結構作用下，地外天體起飛過程中起飛器受到的羽流導流氣動力效應進行了研究，分析了不同上升距離以及偏轉角度條件下，起飛器受到的力矩變化規(guī)律，并對起飛過程中可能出現的現象及原因進行了分析。

1 仿真模型

仿真模型如圖1所示，主體分為起飛器和起飛平臺。起飛器底部包含4個半球形機構、弧形底板和一個主發(fā)動機，主發(fā)動機軸心通過起飛器幾何中心且與起飛器軸心重合，如圖2所示。該發(fā)動機為雙組元發(fā)動機，推進劑為一甲基肼（MMH）/四氧化二氮（NTO）。起飛平臺包括一個圓錐導流結構和平面。

上升距離D定義為發(fā)動機出口距離起飛平面的高度，偏轉角度θ定義為以起飛器質心為中心，沿著Ypc負方向進行偏轉的角度。當發(fā)生偏轉時，半球形機構4靠近起飛平臺，半球形機構2遠離起飛平臺。如圖2所示。

本文主要針對不同上升距離D和偏轉角度θ條件下起飛器受到的羽流導流氣動力效應進行數值模擬研究，具體工況位置如表1所示，其中D=200 mm，θ=3°和5°工況并未進行計算，因為D=200 mm 為起飛初始距離時，偏轉角度 θ不會過大。

2 仿真方法

CFD/DSMC耦合方法已經成熟地應用于羽流的研究中，在羽流的連續(xù)流區(qū)選用CFD方法進行仿真，在稀薄流區(qū)采用DSMC方法進行仿真研究。本文利用這一方法對地外天體起飛過程中起飛器受到的羽流導流氣動力效應進行仿真研究。

2.1 CFD方法

本次研究中選用計算流體學仿真軟件Fluent作為CFD方法的求解器，對連續(xù)流區(qū)進行仿真計算。Fluent采用基于密度的求解器，并選用SST k-ω模型計算湍流黏性系數，通量格式采用二階精度的Roe平均通量差分法（ROE-FDS），時間推進采用下上三角矩陣對稱Gauss-Seidel方法（LUSGS）。

連續(xù)流邊界條件如圖3所示，計算目的是為稀薄流區(qū)DSMC提供參數條件，考慮到羽流流場下游不會影響上游，因此進行連續(xù)流計算的模型中并未建立起飛器模型。在連續(xù)流計算過程中，所有偏轉角度θ=0°工況均采用二維軸對稱計算模型，所有偏轉角度θ不為0°的工況均采用三維計算模型。發(fā)動機入口設為壓強入口，總壓為0.8 MPa，總溫為3040 K。發(fā)動機壁面設置為無滑移絕熱壁面，導流機構壁面為恒溫壁面（300 K）。

為了同時保證計算精度和計算效率，對網格的質量進行了評價，圖4為3種不同網格結果對比，其中黑色對應200萬網格，紅色對應350萬網格，藍色對應500萬網格。從圖中可以看出，200萬網格計算結果與350萬網格的確存在一定差異，約為3.4%，而350萬網格與500萬網格計算結果無明顯差異，約為0.1%。結合計算能力，選取350萬網格即可得到滿意的網格質量和計算效率。

圖3 CFD計算邊界條件Fig.3 Computing boundary conditions of CFD

2.2 DSMC方法

求解DSMC選用的軟件為北京航空航天大學自主研發(fā)的PWS軟件，這一軟件也已經成功應用于羽流的仿真研究中，其精度已經經過多次驗證[18-19]。

DSMC計算所需入口邊界條件由CFD計算所得的結果選取，為了確保入口邊界有效，防止受到下游流場的干擾，入口邊界根據克努曾數Kn以及馬赫數Ma選取，保證Kn＜0.01，Ma＞1。DSMC入口邊界從CFD計算結果中選取三維網格坐標、密度、壓強、溫度、速度、馬赫數和組分摩爾分數等參數作為DSMC稀薄流計算邊界條件。

圖4 3種不同網格壓強結果對比Fig.4 Pressure comparison of three different grid results

圖5 DSMC計算邊界條件Fig.5 Computing boundary conditions of DSMC

稀薄流區(qū)仿真計算域如圖5所示，網格尺寸選擇為當地分子自由程的1/3，每個網格中粒子數不低于15個。計算過程中選取可變硬球模型（VHS）作為二體碰撞模型，并采用純擴散氣表面相互作用模型來達到足夠的表面粗糙度。所有壁面均采用恒溫壁面（300 K），壁面熱適應系數統一選取為1。

3 結果與分析

3.1 仿真校驗

為了驗證第2節(jié)所述計算方法的精度，對文獻[13]在真空羽流實驗系統中進行的120 N雙組元發(fā)動機縮比起飛器羽流效應實驗進行了仿真計算，120 N發(fā)動機燃燒室參數總壓為0.8 MPa，總溫為3 040 K。計算過程中，連續(xù)流區(qū)和稀薄流區(qū)的參數設置與本次計算參數一致，即連續(xù)流區(qū)發(fā)動機壁面設置為無滑移絕熱壁面，導流機構壁面為恒溫壁面（300 K）；稀薄流區(qū)計算過程中選取VHS作為二體碰撞模型，并采用純擴散氣表面相互作用模型來達到足夠的表面粗糙度。所有壁面均采用恒溫壁面（300 K），壁面熱適應系數統一選取為1。

利用上述方法計算得到的120 N發(fā)動機羽流效應結果與實驗結果進行了對比，圖6為實驗中起飛器底部壓強測點位置以及4條仿真對比曲線位置，圖7為實驗結果與仿真結果對比，其中，s為圖6中各位置與中心點的平面距離，p為壓強。

圖7中，Line 1和Line 3曲線位于4個半球形機構位置，壓強較高，仿真結果與實驗結果偏差在10%左右，Line 2和Line 4曲線位于弧形底板位置，此處壓強較低，仿真結果與實驗結果偏差在20%左右，這一偏差主要是由實驗安裝誤差以及發(fā)動機工作狀態(tài)與理想狀態(tài)的偏差導致的。

圖6 4條仿真曲線與實驗壓強測點位置Fig.6 Four simulation curves and experimental pressure measuring point position

圖7 實驗與仿真結果對比Fig.7 Comparison of experimental and simulation results

3.2 連續(xù)流場仿真結果

本文利用Fluent軟件對羽流連續(xù)流區(qū)進行計算，針對18個計算工況，選取部分具有代表性的算例進行流場分析。圖8為Case 1（200 mm，0°）和Case 7（400 mm，0°）工況下羽流連續(xù)流區(qū)壓強云圖。

從對內流場的計算結果中發(fā)現，當上升距離D為200 mm時，由于導流裝置與噴管出口較近，導流裝置對噴管內流動產生較強的影響，發(fā)動機燃氣并未完全擴張至發(fā)動機出口，而是在發(fā)動機擴張段中間位置形成正激波，隨著上升距離D的增加，到達400 mm距離附近，正激波位置逐漸擴張至發(fā)動機出口。定義β為激波至發(fā)動機喉部距離與發(fā)動機擴張段長度的比值；φ為計算得到的發(fā)動機推力與額定推力的比值。β和φ與上升距離D的關系如圖9所示，可以看出，隨著上升距離D的增加，激波位置逐漸由發(fā)動機內部擴張至發(fā)動機外部；當激波位于發(fā)動機內時，發(fā)動機產生的推力較小，隨著激波位置外移，推力逐漸增大，并在激波到達發(fā)動機出口后趨于穩(wěn)定。

圖8 連續(xù)流場壓強云圖Fig.8 Pressure contour of continuous flow field

圖9 β和φ隨上升距離的變化Fig.9 Variation ofβandφwith rising distance

3.3 稀薄流場仿真結果

羽流稀薄流場采用北京航空航天大學自主研發(fā)的PWS軟件進行計算。針對18個計算工況，選取部分具有代表性的算例進行流場分析。圖10為Case 5（300 mm，3°）和Case 13（500 mm，3°）流場云圖。對DSMC計算所得的稀薄流場云圖進行分析可以看出，圓錐導流結構的羽流導流效果較好，但仍有少部分羽流氣體返流至起飛器底面，尤其在半球形機構表面形成局部高壓區(qū)域。

圖11為Case 5（300mm，3°）和Case 13（500mm，3°）起飛器底面云圖。從2組云圖中可看出，起飛器半球形機構和弧形底板局部區(qū)域明顯受到羽流返流影響。

通過對羽流流場和起飛器底部云圖進行分析，可以發(fā)現當起飛器上升距離D較低，偏轉角度θ較小時，靠近起飛平臺的半球形機構4和弧形底板受到羽流作用明顯（Case 5半球形機構4），其表面壓強高于遠離起飛平臺一側（Case 5半球形機構2）；當上升距離D達到一定值，并偏轉較大角度時，遠離起飛平臺一側的半球形機構2和弧形底板受到羽流作用明顯（Case 13半球形機構2），其表面壓強高于靠近起飛平臺一側（Case 13半球形機構4）。

圖10 稀薄流場壓強云圖Fig.10 Pressure contour of rarefied flow field

圖11 起飛器底面壓強云圖Fig.11 Pressure contour of bottom of ascender

3.4 起飛器所受力矩

本次仿真起飛器受到的沿Ypc方向的力矩結果如圖12所示。由于偏轉方向為Ypc負向，因此在上升距離D較低，偏轉角度較小時，力矩Typc為正值。從圖12（a）中可以看出，在偏轉角度 θ=0°時，起飛器受到的Ypc方向力矩基本為0 N·m，這主要是由于起飛器屬于對稱結構；在偏轉角度θ=1°時，起飛器受到的Ypc力矩隨著上升距離D的增加首先為正向力矩，且逐漸減小，但是當上升距離到達500 mm附近，力矩從正向轉向負向力矩，并逐漸負向增加；在偏轉角度 θ為3°和5°時，力矩規(guī)律與偏轉角度為1°時相似，在上升距離D為450 mm左右，力矩逐漸從正向力矩轉變?yōu)樨撓蛄?，然后隨著上升距離D的增加先負向增加，后負向減小。從圖12（b）中可以看出，上升距離D=200～400 mm時，Ypc方向的力矩隨著起飛器偏轉角度θ的增加而正向增大；但是當上升距離D＞500 mm時，Ypc方向的力矩隨著起飛器偏轉角度θ的增加而負向增加。

圖12 Y pc方向力矩變化趨勢Fig.12 Torque variation trend of Y pc direction

由此可見，當上升距離D和偏轉角度θ逐漸增加時，起飛器受到的力矩出現反向增加現象，由糾正偏轉力矩逐漸變?yōu)榧觿∑D的力矩。

3.5 仿真結果分析

圖13 Case 1工況下羽流密度場膨脹波和壓縮波示意圖Fig.13 Schematic diagram of expansion and compression waves of plume density field of Case 1

計算得到的壓強極值位置反向以及力矩反向現象可以通過流場流動特點進行解釋。圖13為Case 1（200 mm，0°）工況下連續(xù)流區(qū)和稀薄流區(qū)的羽流密度場云圖（ρ為羽流密度），從圖中可以看出，由于起飛器與起飛平臺結構復雜，在羽流從噴管流出直至作用于起飛器的過程中，存在較為復雜的波系。在這一復雜的波系中，促使羽流作用于起飛器表面的主要為如下2點：①由于噴管內壓強較高，從而在噴管出口附近形成的膨脹波；②羽流作用于圓錐導流結構后產生的壓縮波。當羽流經過上述膨脹波和壓縮波后，速度大小和方向均會發(fā)生改變，如圖14所示。

圖14中，v為速度，l為羽流與起飛平臺作用位置距離起飛平臺軸線的距離，L為圓錐導流結構底部半徑，α為羽流與起飛平臺作用后壓縮波與導流機構表面形成的偏轉角度。圖15為Case 15（700 mm，0°）工況下馬赫云圖，從圖中可以看出羽流與著陸器表面作用位置即膨脹波與起飛平面作用位置，因此l選取方式較為明確，而壓縮波為曲面，α選取為羽流作用點與壓縮波曲面切線角度。

圖14 經過壓縮波和膨脹波后的羽流速度變化Fig.14 Variation of plume velocity after expansion and compression waves

圖15 Case 15工況下l與α選取示意圖Fig.15 Schematic diagram of l andαselection under working condition of Case 15

在地外天體起飛的過程中，隨著起飛器上升距離D和偏轉角度θ的改變，l與α均隨之改變。偏轉過程中，靠近起飛器一側羽流與起飛平面作用位置始終在圓錐導流結構上，導流效果較好。力矩反向等現象主要受到遠離起飛器一側羽流流動變化影響。圖16給出了遠離起飛器一側羽流與起飛平臺作用點相對位置變化規(guī)律，其中l(wèi)/L為作用點位置與圓錐半徑的比值，當l/L=1時，代表圓錐半徑位置。從圖16中可以看出，隨著上升距離D和偏轉角度θ的增加，遠離起飛器一側的羽流作用點逐漸脫離錐面進入平面內。

圖17給出了遠離起飛器一側壓縮波角度α變化規(guī)律，可以看出：D≤300 mm時，α隨著偏轉角度θ增加而減?。划擠 =400～500 mm時，α隨著偏轉角度θ增加先減小后急劇增加，最后趨于平緩，這主要是由于隨著θ的增加，羽流作用點位置逐漸由圓錐導流結構移動至平面附近；當D=700 mm附近，α隨著偏轉角度 θ增加逐漸增大，當θ達到3°后逐漸降低，這主要是由于700 mm時，羽流作用點在1°時已經位于圓錐邊緣位置附近，并且隨著θ增加逐漸遠離圓錐。

圖16 遠離起飛器一側羽流作用點相對位置變化規(guī)律Fig.16 Variation of relative position of plume impact point on the side far away from ascender

圖17 遠離起飛器一側壓縮波角度α變化規(guī)律Fig.17 Variation of compression wave angleα on the side far away from ascender

通過上述分析可以發(fā)現，上升距離D的增加和偏轉角度θ的增加都將導致發(fā)動機羽流與導流機構的作用點發(fā)生改變，尤其在距離起飛器較遠距離一側，作用點在較大距離和角度下可能由圓錐表面移動至平面，導致壓縮波角度α出現較大幅度增加，根據圖14所示羽流速度變化與角度α的關系，可以看出α較大時，羽流更容易返流至起飛器底面，從而形成壓強極值位置反向以及力矩反向等現象。

4 結 論

本文利用CFD/DSMC耦合的方法，對圓錐導流結構作用下，地外天體起飛過程中起飛器受到的羽流導流氣動力效應進行了研究。分析了不同的上升距離和偏轉角度情況下，起飛器可能受到的羽流導流氣動力影響，得到如下結論：

1）當上升距離較低時，發(fā)動機燃氣會在噴管擴張段中部形成正激波，無法完全擴散至噴管出口，發(fā)動機推力較小，隨著上升距離的增加，激波位置逐漸移至發(fā)動機噴管出口外，推力增加并趨于穩(wěn)定。

2）起飛器底部受到羽流返流影響較嚴重區(qū)域隨著上升距離和偏轉角度的增加逐漸從靠近起飛平臺一側變?yōu)檫h離起飛平臺一側。

3）隨著上升距離和偏轉角度的增加，起飛器受到的羽流導流氣動力矩發(fā)生方向的轉變，如在負向偏轉角度下，由最初的正向力矩逐漸轉變?yōu)樨撓蛄亍?/p>

4）通過對羽流流場中壓縮波和膨脹波的位置變化及其對羽流流動影響進行分析，發(fā)現力矩反向等現象產生的原因為發(fā)動機羽流與起飛平臺的作用點從圓錐導流結構表面移動至平面，從而使起飛平臺表面壓縮波角度增加，羽流作用于起飛平臺后的流動方向由貼近起飛平臺向側面流動急劇轉變?yōu)榉磸椫疗痫w器底面方向流動。

研究結果指出了地外天體起飛過程中圓錐導流結構可能引起的羽流效應影響，為后續(xù)的研究及設計工作提供有效的參考。根據仿真結果，建議在起飛過程中，利用額外姿控發(fā)動將起飛器姿態(tài)控制在偏角較小范圍內。