高空飛行環境中液體運載火箭底部熱環境研究

周志壇，丁逸夫，樂貴高，梁曉揚，孫培杰，盛英華

(1.南京理工大學機械工程學院，南京 210094；2.上海宇航系統工程研究所，上海 201109)

0 引 言

近年來，隨著載人登月、空間站組建、深空探測等一系列重大航天活動相繼開展，我國加大了對重型運載火箭及其動力系統的研究工作，確定了液氧煤油發動機和液氫液氧發動機是運載火箭動力系統的最佳選擇[1-3]。液體運載火箭發動機在工作時，由于外界環境壓力過低，燃氣流進入外界環境后急劇膨脹，在火箭底部形成回流，并對箭體底部形成對流加熱效應，同時，高溫燃氣中噴射熾熱的CO2和H2O等混合氣流對火箭底部形成輻射加熱效應[4]。火箭底部相當于背風面，容易受發動機噴流回流造成的對流加熱與輻射加熱耦合作用，導致溫度迅速升高。對火箭底部熱環境估計過低會給火箭整體的安全性帶來極大威脅，甚至誘發爆炸等重大事故造成飛行失敗，而估計過高又會導致熱防護設備保守設計，增加發射成本[5-7]，因此，對液體運載火箭底部進行熱環境分析研究成為當前的重中之重。

國外已有大量學者對火箭底部熱環境開展研究工作。1970年，Kramer[8]對泰坦III底部熱環境進行試驗測試，結果顯示在剛起飛1 s內輻射熱流是總熱流的主要來源，隨著飛行高度上升，輻射熱流減少，對流熱流升高；1994年，Reijasse等[9]對阿里安5進行數值模擬和縮比試驗對比研究，得出影響對流熱流的主要因素是飛行高度和噴管數量以及噴管排列方式；2007年，Negishi等[10]對H-IIA底部熱環境進行數值模擬，結果顯示底部熱流峰值在30 km 高空附近。

近年來國內也有部分學者開展了噴流試驗以及數值仿真方面的工作。史亞男[11]采用有限體積法求解RANS方程，湍流模型選取RNGk-ε兩方程模型對火箭底部熱環境進行仿真；周一鵬等[12]利用反向蒙特卡羅法計算火箭發動機噴管內的輻射熱流；戴欣怡[13]采用標準k-ε兩方程模型分析了液體發動機噴流對底部熱環境的影響。

相比于國外對火箭底部熱環境研究程度，國內目前還處于早期階段，尤其是我國液體運載火箭飛行試驗數據極少，國內學者對火箭底部熱環境多為數值模擬方面的研究，缺乏對比分析導致仿真方法的有效性難以得到檢驗。并且由于計算資源有限，所以大多數學者對數值模型進行了不同程度簡化，如幾何模型只考慮四分之一、網格無邊界層、熱輻射計算中采用了精度較低的P-1模型等，這些都在很大程度上影響了計算精度。

隨著計算機運算速度和內存規模飛速發展，為開展含邊界層網格的結構化全箭三維流動數值模擬創造了有利條件。本文對液體運載火箭底部熱環境進行數值模擬研究，在與試驗測試完成對比驗證基礎上，通過大量研究分析不同飛行高度環境下，火箭底部熱流的變化規律。

1 飛行試驗

試驗測試對象為某三級液體運載火箭，火箭共分為三個子級，由于液氧煤油發動機密度比沖高，因此一子級推進劑選擇為液氧和煤油。在火箭底部放置兩個熱流傳感器，傳感器位置如圖1所示。

測試結果為該火箭按照實際軌道飛行過程中，箭體底部兩個熱流傳感器反饋的數據。

圖1 火箭底部熱流傳感器位置Fig.1 The heat flux sensor of rocket tail compartment

2 計算模型

2.1 計算網格

采用多重網格方法對三維模型進行網格劃分，該方法通過在不同疏密網格下求解原方程組，平順高、低頻誤差，加速收斂過程。

為確保激波能夠完整描述，對噴管及噴流區域網格進行加密處理，整個計算域均采用結構網格以保證網格的正交性及光順性。在數值計算同時，為防止邊界處的溫度、壓力等物理量變化劇烈，本文對邊界層計算域的網格進一步加密，確保壁面附近網格Y+值小于3，具體如圖2所示。

圖2 計算網格Fig.2 Computational mesh

2.2 數值方法

2.2.1控制方程

基于液體運載火箭多組分燃氣滿足連續介質假設和理想氣體狀態方程的特點，建立三維多組分可壓縮Navier-Stokes方程，包括質量、動量與能量輸運的控制方程，形式如下。

多組分輸運方程為：

(1)

式中：Yl為組分l的質量分數，Rl為組分l在化學反應后的凈生成率，Sl為自定義源項的離散相引起的生成率。

組分擴散通量Jl可表示為：

(2)

式中：Dl，m為組份l的質量耗散系數，DT，l為組份l的熱擴散系數。

在直角坐標系中，單一組分l的可壓縮N-S方程的守恒形式為：

(3)

(4)

(5)

(6)

式(3)～(6)中:U為流動變量；F,G,H為氣流通量矢量;Fv,Gv,Hv為黏性通量矢量;K為熱傳導系數；T為環境溫度；p,ρ,e,τ,μ分別為壓力、密度、比動能、應力、黏性系數;u,v,w分別為速度在x,y,z方向上的分量。

2.2.2空間離散格式

本運載火箭采用了大推力發動機設計方案，燃燒室總壓、總溫高，噴流速度為典型超聲速，高空階段箭體以超聲速飛行，來流與噴流強烈干擾而形成復雜的發動機噴流效應，二階迎風TVD格式收斂速度快，數值穩定性高，由于設計了魯棒性強的通量限制器，不僅能求得高分辨率的復雜波系結構，而且抑制強間斷解的非物理震蕩，擬選用二階迎風TVD格式計算全箭這類大型噴流模型。

將含箭體的整個計算域劃分為六面體的結構化網格，采用有限體積法對每個體網格單元進行數值積分，流動變量U在體單元中心的體積平均為：

(7)

式中:Vi,j,k為計算體網格單元，采用二階迎風TVD格式對對流通量進行離散。

(8)

其中，數值流通量分別為：

(9)

式中:λk，μk，νk為線性化替換矩陣的特征值;αk,βk,γk為線性化替換矩陣的展開項系數;ek(A),ek(B),ek(C)為線性化替換矩陣的特征向量，黏性通量采用二階中心差分離散。

2.2.3湍流模型

湍流模型采用Realizablek-ε兩方程模型，與標準k-ε兩方程模型相比較，湍動耗散率ε方程中的產生項不再包含有湍動能k方程中的產生項Gk，并且湍動黏度μt中的系數Cμ不是常數，而是與應變速率相關，這樣的形式能更好地表示能量轉換。

湍動能k方程為：

(10)

式中:μ為混合物黏性，Gk為平均速度梯度引起的湍動能k的產生項，常數系數σk=1.0。

湍動耗散率ε方程為：

(11)

式中:常數系數σε=1.2；C1=1.44，C2=1.9，E為平均特征應變率。

2.2.4輻射模型

對輻射熱流采用離散坐標法(DOM)進行求解，該方法具有易于處理散射問題，易于和流動方程聯立求解以及計算精度較高的特點。該方法最早由Chandrasekhar[14]提出，由輻射傳遞方程沿著s方向進行離散得到：

在r位置沿著s方向的輻射傳遞方程為：

(12)

由式(12)可得光譜強度Iλ(r,s)的輻射傳遞方程為：

(13)

式中:r為位置矢量，s為方向矢量，s′為散射方向矢量，δ為氣體層厚度，a為吸收系數，n為折射率，σ為黑體輻射常數，σs為散射系數，I為光譜輻射強度，T為黑體的熱力學溫度，Φ為散射相函數，Ω′為立體角，λ為波長，Ibλ為普朗克方程給定的黑體強度。

氣體的折射率n會受環境溫度Tδ以及波長λ影響，具體如下：

(14)

式中:T0為標準溫度，n0為氣體在標準溫度下的折射率，Tδ表示環境溫度，η為線脹系數，可由文獻[15]查得。

吸收系數a及散射系數σδ可表示為：

(15)

式中:kλ為光譜減弱系數，βλ為散射衰減系數，ρ為氣體密度。

對散射角為θ的散射相函數定義為Φ(θ)，離散坐標法求解散射光強正是利用勒讓德多項式展開，對Φ(θ)做近似處理，展開如下：

(16)

式中:Φg是第g階勒讓德函數，cg是第g階展開項的系數，可表示如下：

(17)

由式(14)可知，要獲得精確的散射相函數需要無窮多階勒讓德多項式展開，在實際計算時，由于只在前向周圍很窄的區域內存在顯著的前向峰，因此，將散射相函數近似表示為前向部分和截斷相函數之和，即：

(18)

式中:f為截斷因子，代表前向散射所占比例，M為近似的階數，cg*為第g階近似展開項系數。

3 數值結果

液體運載火箭從地面發射后，高度不斷提升，隨海拔變化的主要參數包括環境壓強、環境溫度以及飛行馬赫數。

為研究液體運載火箭底部熱環境變化趨勢，本文選取8個典型高度的工況進行數值模擬，表1為這8個工況的自由流條件參數，噴管喉部總溫為3800 K，總壓為18.66 MPa，氣體組分包括H2O,CO2,CO,H2以及空氣等。

3.1 不同高度下的馬赫數場分布

由圖3可知，在5 km～25 km高度，噴管尾焰中馬赫波系結構較為明顯，特別是在5 km和10 km 時，燃氣流場中可見多個噴流波節。由于發動機總壓遠大于外部環境壓強，因此燃氣在噴管出口附近形成桶形壓縮波，在馬赫數云圖中可觀察到清晰的尾焰噴流核心區、發展段以及射流邊界。

表1 自由流條件參數Table 1 Parametric free stream conditions

圖3 不同高度下的馬赫數云圖Fig.3 Contour of Mach number at different altitudes

隨著高度逐漸爬升，大氣環境壓強降低，噴流對外界環境的擾動區域變大，噴流馬赫數峰值增大，噴流波節數減少，桶形激波與噴管軸線的交點坐標向下游移動。在35 km～40 km，噴流膨脹半徑及膨脹距離劇增，桶形激波與噴管軸線不相交。

燃氣在拉瓦爾噴管中處于超臨界工作狀態，在擴張段形成超聲速膨脹波并進一步加速，該膨脹波在噴管尾部由于受到大氣壓縮作用呈桶形，在噴流激波和噴流邊界下游與噴管軸線相交。隨著高度上升，環境壓強和溫度呈下降趨勢，外部伴隨氣流對噴流的壓縮作用逐漸減弱，噴流的膨脹程度繼續加強，這與圖3情況一致。

3.2 不同高度下的底部溫度場分布

選取4個典型工況的底部溫度場進行對比展示，圖4為10 km,20 km,30 km和40 km這4個工況的箭體底部及噴管附近溫度場等值線分布圖。

圖4 不同高度下的溫度等值線圖Fig.4 Contour of temperature at different altitudes

從圖4中火箭底部溫度等值線的顏色及疏密程度可以看出，隨著飛行高度提升，火箭底部溫度整體為逐漸上升的趨勢。隨著海拔增加，噴管尾部噴流膨脹加強，膨脹后的噴流邊界與超聲速來流相互作用形成剪切層，噴管噴出的高溫燃氣在剪切層流動作用下回流至火箭底部，形成對流熱流，自由來流、發動機噴流激波與箭體底部及噴管壁面構成一個封閉區域，對流熱流在封閉的剪切層區域內形成再循環系統，對箭體底部持續進行對流傳熱，造成底部對流熱流值在前期隨著海拔上升逐漸升高。但隨著海拔繼續上升，尾部噴流激波膨脹劇烈，致使所構成的封閉區域減小，回流至箭體底部的高溫燃氣量減小，這也是后期箭體底部對流熱流下降的原因。

3.3 不同高度下的底部對流熱流分布

選取4個典型工況的底部熱流云圖進行對比展示，從圖5可以看出，20 km處對流熱流大于10 km處及30 km處，而40 km處對流熱流值最小。

圖5 不同高度下的對流熱流云圖Fig.5 Contour of convection heatflux at different altitudes

隨著海拔增加，環境壓力下降，噴流膨脹加強，引起更多回流，火箭底部的對流換熱增加，另一方面，海拔升高時，前方來流馬赫數(飛行速度)增加，對箭體底部回流的阻滯效果加強，因而減弱了火箭底部的對流加熱。在低海拔時，環境壓強對對流熱流的影響大于來流速度，所以隨著海拔提升對流熱流表現為增加的趨勢；在高海拔地區，環境壓強對對流熱流的影響小于來流速度，所以隨著海拔提升對流熱流表現為降低的趨勢。

4 數值模擬與飛行試驗對比

在對某三級液體運載火箭數值模擬過程中，選擇了與熱流傳感器相同的位置布置兩個監測點，通過數值計算得到了監測點的熱流參數，圖6顯示了兩個監測點的數值計算總熱流值與飛行實際測量總熱流值變化曲線，可以看出兩者趨勢以及數值大小吻合度都較高?？偀崃鞒霈F了兩次峰值，分別是72 s 和99 s，對應的正是10 km附近和20 km附近。

圖6中顯示的計算總熱流值可分解為圖7中的對流熱流和輻射熱流。圖7表明,72 s附近的總熱流峰值主要來源于輻射加熱，99 s附近的總熱流峰值主要來源于對流加熱。對流熱流在50～99 s逐漸增加，并在99 s達到最大值，然后在99～133 s逐漸減小，輻射熱流在剛開始達到最大值，然后開始降低。

圖6 總熱流對比Fig.6 Comparison of total heating rate

圖7 數值計算的對流熱流以及輻射熱流Fig.7 Computed convection and radiation heating rate

從圖7還能看出，在較低飛行高度時，輻射熱流占總熱流的比重較大，而在更高的海拔時，輻射熱流逐漸減小，對流熱流在火箭底部熱環境中扮演更加重要的角色。因為液體運載火箭輻射傳熱介質主要為H2O和CO2，隨著飛行高度上升，環境含氣量降低，只能由燃氣排放的H2O和CO2作為介質進行輻射傳熱，所以輻射熱流逐漸減小，同時，隨著海拔升高，火箭燃氣噴流轉變為高度膨脹，來自噴管附近的高溫氣體回流至火箭底部，對底部產生對流加熱，因而對流熱流增加，當海拔繼續爬升時，外界來流馬赫數增加，來流對回流的阻滯效果加強，削弱了火箭底部的對流換熱，致使對流熱流下降。

5 結 論

本文對液體運載火箭的高空底部熱流問題進行了數值模擬，并與實際飛行數據進行了對比分析，主要得出以下結論：

1)基于合理假設，建立含燃氣-空氣多組分Navier-Stokes方程、Realizablek-ε兩方程模型和熱輻射模型的火箭燃氣噴流模型，采用高效的二階TVD格式離散對流通量，離散坐標法求解熱輻射方程，中心差分求解黏性引起的耗散項。對液體運載火箭高空底部熱流進行數值研究，仿真結果與測試數據吻合良好，校驗了本文所采用的數值方法的精度和有效性。

2)本文的三級液體運載火箭底部總熱流出現了兩次峰值，第一次峰值在10 km附近，主要受到輻射熱流影響所致；第二次峰值在20 km附近，主要由對流熱流所致。

3)分別對8個典型飛行高度的工況的熱流變化規律進行研究，火箭底部的輻射熱流在剛起飛階段達到最大值，然后隨著海拔上升逐漸降低，對流熱流表現為先上升后下降的趨勢，并在20 km附近達到最大值。

4)在分析液體運載火箭底部的熱環境時，需要考慮到對流/輻射耦合換熱的影響，本文所形成的火箭高空噴流模型和熱流數值模擬方法對液體運載火箭底部熱防護安全性設計具有一定的參考價值。