周印佳，張志賢，付新衛(wèi)，阿嶸

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

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

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

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

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

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

1 物理模型與計算方法

1.1 氣動熱環(huán)境

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

(1)

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

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

(2)

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

1.2 熱化學燒蝕計算模型

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

(3)

(4)

圖1 燒蝕材料表面能量平衡及內(nèi)剖面原理示意圖

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

(5)

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

(6)

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

正常質(zhì)量流量的計算式為

(7)

未校正質(zhì)量流量B′0的計算式為

(8)

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

(9)

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

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

(10)

(11)

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

(12)

(13)

式中：R為通用氣體常數(shù)。

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

(14)

(15)

1.3 Landau變換

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

(16)

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

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

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

(17)

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

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

r=Ro-z,z=Ro-r

(18)

而η坐標系和r坐標系的關(guān)系可以表示為

(19)

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

(20)

(21)

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

因此，任意節(jié)點的移動速度可以表征為與表面后退率相關(guān)的公式

(22)

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

圖3 均布收縮網(wǎng)格方案對比

1.4 有限差分計算

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

(23)

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

(24)

(25)

(26)

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

(27)

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

(28)

1.5 計算流程

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

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

圖4 計算流程圖

2 結(jié)果與討論

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

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

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

表1 碳-碳材料物性

(29)

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

圖5 表面溫度計算結(jié)果對比

圖6 表面燒蝕后退量計算結(jié)果對比

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

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

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

表2 計算結(jié)果對比

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

3 結(jié) 論

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

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

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