再入飛行器燒蝕熱防護一體化計算方法

周印佳，張志賢，付新衛，阿嶸

中國空間技術研究院 錢學森空間技術實驗室，北京 100094

燒蝕熱防護是再入飛行器最成功的防熱方法之一，但燒蝕熱防護系統響應的計算是也非常復雜的，涉及到平衡或者非平衡化學反應、燒蝕機理、燒蝕外形變化、燒蝕率和燒蝕量的預測以及結構熱響應等諸多因素[1-2]。國外使用最廣泛的燒蝕預測數值方法涉及到計算流體動力學(Computational Fluid Dynamics，CFD)的邊界層法(BLIMP程序)和材料熱響應的一維內部熱傳導、CMA程序(Aerotherm’s Charring Material Thermal Response and Ablation Program)。隨著計算機硬件的發展，近年來基于Navier-Stokes(N-S)方程方法耦合材料熱響應的預測逐漸增多。文獻[3]利用N-S方程求解軸對稱再入流場與材料熱響應的耦合，基于表面能量平衡方法計算出表面溫度和燒蝕率，流場計算和材料響應計算以準靜態或瞬態方式在每個時間步上進行耦合。計算中表面會發生后退和網格動態變化，并且流場網格在固定的彈道點上進行重畫以適應新的表面形狀。文獻[4-5]將軸對稱N-S程序與CMA程序耦合，采用松耦合的方法迭代得到CFD計算的熱流、CMA計算的表面溫度以及燒蝕率引起的表面溫度和流率變化低于1%。計算中允許燒蝕產物與非平衡氣體反應，但是表面不產生后退。文獻[6]采用類似的計算方法對軸對稱流場進行計算，并且計算中允許表面后退。文獻[7-8]將軸對稱熱化學非平衡流場計算和二維固體熱傳導計算耦合。流體和固體的耦合通過表面質量和能量平衡實現，并假設了一個包含表面化學非平衡的熱化學燒蝕模型。文獻[9]采用松耦合算法耦合了熱化學非平衡N-S程序和考慮表面燒蝕和形變的多維內部熱傳導。將N-S程序得到的氣動參數結果視為常數，材料響應程序按時間推進計算燒蝕、形變表面溫度和內部溫度。但是該方法中燒蝕過程的高度非線性導致熱傳導和燒蝕率計算的嚴重不穩定性，為了解決該問題，文獻[10]開發了一種迭代技術以消除前面所述的不穩定性，并顯著提高燒蝕過程的精度。

國外的氣動熱預測和材料響應計算通常采用在SACCARA[9]、GIANTS[11]、COYOTE II、FIAT等內部代碼的基礎上進行修改和應用的方式。而近年來得益于商用CFD和有限元軟件的快速發展，采用商業軟件耦合預測氣動熱和材料熱響應的研究也逐漸增多。文獻[12]利用Fluent與LS-DYNA商用軟件，以松耦合的方式求解了零攻角軸對稱飛行器的非平衡化學反應流動，并采用網格運動算法實現再入飛行器燒蝕導致的表面后退模擬。

以上數值方法通過求解N-S方程和高精度的CMA程序或有限元程序，提高了計算氣動參數和燒蝕響應的精度，但是該方法計算負擔較重，對計算能力要求較高，并且CMA高精度程序很復雜，尤其是考慮到程序輸入所需要的大量數據更是如此。另外，生成表面化學輸入和熱解氣體焓需要運行額外的平衡或者非平衡化學程序。化學程序本身需要大量的輸入，并需要材料組分以及熱解氣體組分的相關知識。為了完成這種分析，需要開發熱響應模型并與電弧風洞試驗數據相關聯，還需要一組可用的包含大量壁面條件的表面化學表，還要開發熱解氣體焓模型。如果這些模型都沒有，熱解氣體和熱防護材料的確切成分都未知，建立和運行像CMA這樣的高精度模型無論是從準確度還是從結果收斂性上都是有問題的。

當高保真的熱燒蝕響應模型無法實現時，尤其是在需要進行大量快速計算的初樣設計階段，需要多次迭代以求得防熱層厚度，相對不那么復雜的計算方法顯然更實用。這些近似方法通常采用穩態燒蝕假設并采用燒蝕熱Q*來預測再入過程中的燒蝕后退量。這些方法常采用近似解析方程來確定材料內部的瞬態一維熱傳導。采用近似熱傳導方程求解在指定的熱防護系統結構界面處溫度條件下的所需絕熱層厚度。這些方法能夠給設計者提供大概的答案，但是其缺點在于這些方法的穩態燒蝕假設在飛行中通常是不成立的，一般會高估后退量。另外一個不足是采用的近似瞬態熱傳導方程一般只在半無限固體或平板假設成立時才有效。

本文所提出的計算方法旨在提供一種具有較高精度，并能夠快速應用的燒蝕熱防護一體化計算手段，適用于缺少高精度材料熱響應模型的初步設計階段。該方法利用定義好的飛行器幾何數據，采用工程算法估算駐點的對流熱流和輻射熱流，整合具有較高精度的燒蝕模型，通過Landau變換簡化燒蝕后退帶來的節點刪除過程并保證空間離散精度，最后求解瞬態有限差分熱傳導來獲得不同尺寸、形狀和材料的熱防護內部響應。

1 物理模型與計算方法

1.1 氣動熱環境

采用工程算法預測駐點的對流熱流和輻射熱流。Sutton-Graves冷壁對流熱流計算公式[13]為

(1)

式中：K為大氣化學的函數，針對地球大氣和火星大氣分別取為K=1.741 5×10-4和K=1.902 7×10-4;Rn為球頭半徑；ρ∞為自由來流密度；V為飛行速度。式(1)已經過眾多地面和飛行驗證，對于鈍頭體的熱流預測偏差通常在5%～10%。

通常氣體溫度在8 000 K以下時輻射熱流可以忽略不計，但是隨著再入速度的增大，輻射熱流所占比重逐漸增大，例如Apollo飛船再入輻射熱流可占總熱流載荷的40%～60%[14]，此時輻射熱流不可忽略。駐點的輻射熱流采用Tauber-Sutton公式[15]計算，該公式適用于地球或者火星再入。Tauber-Sutton熱流公式的表達式為

(2)

在彈道上的每個點，將對流熱流和輻射熱流相加得到總加熱熱流。

1.2 熱化學燒蝕計算模型

如圖1所示,燒蝕材料表面能量平衡方程可以表示為

(3)

(4)

圖1 燒蝕材料表面能量平衡及內剖面原理示意圖

qnet=qconv(1-hw/hr)φHALφblow+qrad-

(5)

式中：hr為恢復焓；ε為輻射系數；σ為斯忒藩-玻耳茲曼常數；Tref為參考溫度；φHAL為侵蝕系數，由Hove-Shih理論可知φHAL≥1[17]，對于清潔大氣飛行取φHAL=1[18]；φblow是由氣體進入邊界層導致的，根據邊界層的薄膜理論可以由下式計算[19-20]

(6)

式中：吹氣系數a′取值范圍為0.3～1.3[21]，對于碳-酚醛材料取a′=1.3，碳-碳材料取a′=1.5[18]。

正常質量流量的計算式為

(7)

未校正質量流量B′0的計算式為

(8)

未校正的熱傳遞系數gh0=(qconv/hr)φHAL，熱傳遞系數gh=gh0φblow。壁內焓hu的計算式為

(9)

式中：Cp為氣體比熱；C∞=2.3 J/(g·K)；D=800 K；Tref=300 K。

為了計算燒蝕速率，可以采用碳-空氣表面熱化學平衡計算的簡化擬合公式[22-23]:

(10)

(11)

在最低燒蝕溫度范圍內可以采用非平衡反應速率，碳質材料的速率控制氧化有多種表達形式：

(12)

(13)

式中：R為通用氣體常數。

將速率控制表達式和平衡表達式合并成總燒蝕率為

(14)

(15)

1.3 Landau變換

對于常見的鈍頭熱防護罩外形，一維模型就能夠滿足要求。軸對稱幾何的一維熱傳導方程為

(16)

式中：λ為固體熱導率；Cp,s為固體比熱;r為徑向。

Landau變換的最初發展是用于解決一維平板的相變問題，這里采用Landau坐標系[24]實現燒蝕的動網格。該坐標系重新定義了空間坐標，從而使在Landau空間內剩余區域的無量綱厚度始終為1，如圖2所示并且在求解過程中某一點的Landau坐標始終為常值。圖2中z坐標系的原點固定在初始燒蝕界面處，而x坐標系與瞬時燒蝕界面相關聯。任意一點的Landau坐標[25]為

圖2 一維平板模型坐標系示意圖

(17)

式中：η為Landau坐標；L(t)為隨時間變化的防熱層厚度。

對于軸對稱幾何外形，z坐標系和r坐標系的關系可以表示為

r=Ro-z,z=Ro-r

(18)

而η坐標系和r坐標系的關系可以表示為

(19)

利用式(19)和微積分中的鏈式法則，熱傳導方程可以轉換成以下形式:

(20)

(21)

式中：j為節點號；n為時間節點；Δt為時間步長。

因此，任意節點的移動速度可以表征為與表面后退率相關的公式

(22)

燒蝕問題的網格方案通常有平移網格和收縮網格2種方案，平移網格方案的節點位移關系相對來說要更簡單，但是會在背壁面位置的網格間產生很大的不連續，如圖3所示，可能會導致在空間離散精度上產生負數值結果。而本文采用的基于Landau變換的收縮網格方案避免了復雜的節點刪除過程，并且在整個求解過程中可以保證空間離散精度，因為該方案的節點數和無量綱單元厚度Δzj(t)/L(t)是常數[26]。

圖3 均布收縮網格方案對比

1.4 有限差分計算

變換后的熱傳導方程采用有限差分法進行隱式離散化，對時間的偏導數采用后向差分格式：

(23)

對空間坐標的偏導數全部采用中心差分格式：

(24)

(25)

(26)

對于平板模型有式(20)中的m=0，將式(23)～式(26)代入到熱傳導方程中，得到

(27)

外表面(η=1)為熱流邊界條件qnet(nΔt)，并假設內表面(η=0)給定絕熱邊界條件，最終平板形式熱防護的Landau坐標系下的熱傳導方程可化為如式(28)的矩陣形式，同理亦可得到圓柱和球頭形式熱防護在Landau坐標系下的熱傳導方程。公式推導過程比較冗長，此處不再詳細給出。

(28)

1.5 計算流程

計算的有限差分網格采用均勻網格，初始給定一個防護層厚度以及所需的最大節點間距Δx，并計算出能夠保持節點間距略低于或等于所需最大間距的節點數量，然后將節點編號存儲在一個數組內。

另外，該方法也可以利用割線法進行迭代，用來預估防熱層厚度。只要給定飛行軌跡和背壁溫度約束，就可以根據防熱層厚度的初始預估值按照圖4所示流程進行熱響應計算；然后根據割線法的需要，防熱層厚度比初始預估值增加0.1 cm再進行一輪計算，如果背壁溫度與之前相比沒有顯著差異，則將初始預估值減小一半，并重復計算過程。一旦這些初始計算完成，割線法就可以按部就班地校正防熱層厚度，直到計算所得的背壁溫度與要求的背壁溫度約束相差小于0.1 ℃。該過程通常經過10～15次迭代即可收斂。

圖4 計算流程圖

2 結果與討論

為了驗證本文提出的計算方法，以典型鈍頭體的地球再入和PICA電弧風洞模擬算例分別針對相對高密度和中低密度的材料進行了計算和對比驗證。

2.1 碳-碳材料鈍頭體地球大氣再入

高彈道系數(95 760 N/m2)的鈍頭再入飛行器再入地球大氣，高度在30 s內從92 km降低到海平面。再入初始速度為6 000 m/s，駐點半徑Rn=0.05 m。熱防護材料為碳-碳材料，密度為1 900 kg/m3，材料的熱物性參數如表1所示。來流的密度和速度與時間的關系為

表1 碳-碳材料物性

(29)

利用以上再入條件將本文方法與熱平衡積分(Heat Balance Integral，HBI)算法的計算結果[18]進行了對比分析，圖5和圖6分別為壁面溫度和燒蝕后退量的結果對比。2種方法計算的峰值溫度基本一致，預測的燒蝕后退量相差約0.2 mm。可見本文算法與經典的熱平衡積分算法計算結果吻合較好，能夠滿足工程需要。

圖5 表面溫度計算結果對比

圖6 表面燒蝕后退量計算結果對比

2.2 PICA材料電弧風洞燒蝕模擬

為了測試星塵號(Stardust)飛船的熱防護材料PICA的燒蝕性能，NASA開展了一系列電弧風洞試驗，同時開展了一系列相應的仿真分析[28]，對實際燒蝕表面溫度與試驗中熱量計測量結果的差異進行對比分析。試驗和仿真針對4種不同工況條件，試驗中氣體為干燥空氣加7%質量的氬氣，氬氣的濃度是試驗條件與實際飛行條件最大的不同之處。通過測量銅制熱量計的傳熱速率得出氣體最小焓值為30 MJ/kg，而通過分析激波層輻射的頻譜得出焓值為40.6±1.4 MJ/kg。

在所有工況計算中，與文獻[28]保持一致，PICA材料表面溫度均設為3 000 K，PICA材料參數為：密度ρs=240 kg/m3[28]，熱擴散系數為 7×10-7m2/s[29]。表2給出了利用本文方法計算的燒蝕速率和燒蝕后退率與文獻計算結果的對比。

表2 計算結果對比

在熱流和壓力相對較低的工況1和工況2條件下，本文方法預測結果與文獻結果都吻合較好，燒蝕后退率和燒蝕速率的偏差分別為18%、13.8% 和9.1%、8%。而在壓力更高(60.8 kPa)的工況3和工況4條件下，后退量預測出現了較大差異。本文方法預測的結果明顯低于文獻計算結果，燒蝕后退率和燒蝕速率的偏差分別達到了25%、25.5%和47.7%、47.4%。研究表明[30]，當前采用的燒蝕模型對于非碳化燒蝕材料具有較高精度，而對于碳化燒蝕材料則具有一定局限性，需要進一步完善。PICA材料以碳化燒蝕為主，在高壓力高熱流條件下，會發生大量的碳化和熱解，燒蝕模型需要更精細地考慮材料碳化、熱解能量吸收、熱解氣體通過固體時的對流以及表面化學作用的影響。另外，PICA這樣的中低密度燒蝕材料的燒蝕性能對壓力高度敏感[31]。因此，針對此類具有高壓力高熱流的中低密度燒蝕材料的燒蝕計算在實際應用之前需要經過大量的驗證。

3 結 論

1)本文建立的再入飛行器燒蝕熱防護系統燒蝕與瞬態溫度耦合響應一體化預測方法，能夠有效預測平板、圓柱和球體熱防護材料包括氣動熱、燒蝕后退、瞬態溫度響應在內的動態響應變化。

2)通過仿真對比和已有研究成果表明，對于具有較高密度的材料其燒蝕后退較小，本文方法預測的燒蝕及熱響應比較準確。而對于以碳化燒蝕為主的PICA材料，當前采用的燒蝕模型具有局限性，需要進一步完善。隨著壓力和熱流的增大，低密度材料(例如燒蝕性能對壓力高度敏感的PICA材料)的預測偏差逐漸增大。

3)在使用未經驗證的材料響應模型時有必要保持謹慎。在碳-碳材料算例中，本文方法與文獻數據吻合較好。對于其他具有不同熱防護系統材料的飛行器還需要進一步驗證，尤其是采用低密度材料熱防護系統的飛行器需要進一步研究。