，，，

北京航天長征飛行器研究所 高超聲速飛行器防隔熱中心，北京 100076

高超聲速飛行器在臨近空間長時間、高機動飛行時，高溫真實氣體效應、化學非平衡效應、激波/邊界層干擾，以及轉(zhuǎn)捩等特征明顯，給氣動熱環(huán)境預示帶來一系列復雜問題,特別是更難以準確預測邊界層轉(zhuǎn)捩。地面風洞試驗目前僅能做到小于Ma12的模擬，因此，直接測量高超聲速飛行器表面熱流密度這一關(guān)鍵參數(shù)，對摸清真實飛行條件下高超聲速飛行器熱環(huán)境和飛行器熱防護設(shè)計都有著極為重要的意義。國外一些研究者在該方面開展了許多研究[1-5]。典型的是德國SHEFEX-Ⅱ測得Ma9.3飛行熱流數(shù)據(jù)，受傳感器使用溫度不能超過700℃的限制，僅測得飛行器尾部熱流。國內(nèi)在地面風洞試驗技術(shù)和熱流測量技術(shù)方面開展了許多研究[6-9]，獲得了Ma10以下大量風洞熱流測量數(shù)據(jù)。樹脂基材料由于密度小、導熱系數(shù)低，被廣泛應用于高超聲速飛行器熱防護結(jié)構(gòu)。在高溫氣動加熱條件下，通過材料的燒蝕/熱解/碳化及對邊界層的“熱阻塞”作用，實現(xiàn)對高超聲速飛行器熱防護目的。國內(nèi)外學者針對樹脂基材料、蜂窩材料溫度場計算方法進行了大量的理論和試驗研究[10-23]，但未研究樹脂基復合材料導熱微分方程中加入導熱系數(shù)隨溫度變化的燒蝕溫度場計算方法，也未定量給出碳化層厚度及計算結(jié)果與實測結(jié)果的比較。

本文采用在飛行器表面開孔安裝長時耐高溫熱流傳感器直接測量熱流密度方法，國內(nèi)首次獲得Ma12以上飛行器表面臨近空間熱流密度真實變化歷程和邊界層轉(zhuǎn)捩特征，驗證了熱環(huán)境預示方法的正確性。在考慮熱解吸熱的樹脂基常規(guī)固體導熱微分方程中加入材料導熱系數(shù)隨溫度變化項后的導熱微分方程，計算了飛行器熱結(jié)構(gòu)分層溫度和燒蝕碳化層厚度，并與實測飛行器熱結(jié)構(gòu)燒蝕碳化量和分層溫度進行了比較，驗證了修改計算方法的合理性和正確性。

1 錐面熱流密度

再入飛行器以高超聲速在大氣層飛行時，受氣體粘性作用影響飛行器表面相對速度為零，在緊貼近飛行器表面薄層里，速度由零迅速增至來流速度，此薄層即為邊界層。在邊界層內(nèi)氣體動能部分轉(zhuǎn)化為熱能，存在很大的溫度梯度，由此溫度梯度形成了對高超聲速飛行器表面的氣動加熱。采用參考焓方法計算有攻角再入飛行器錐面熱流密度。在轉(zhuǎn)捩區(qū)內(nèi)熱流密度用層流值和湍流值加權(quán)平均方法。

軸向X位置的錐面熱流密度qw：

(1)

式中：qx為不考慮焓值影響的錐面熱流，按照層流、湍流、轉(zhuǎn)捩3種狀態(tài)分別計算；hre為恢復焓；hw為壁焓。層流狀態(tài)下：

(2)

(3)

式(2)中帶“*”參考值由當?shù)匚锩鎵毫蛥⒖检蚀_定，h*參考焓和hre恢復焓分別為：

h*=0.19hre+0.23he+0.58hw

(4)

(5)

湍流狀態(tài)下：

當Re<107，

(6)

當Re≥107，

(7)

h*=0.22hre+0.28he+0.50hw

(8)

轉(zhuǎn)捩狀態(tài)下：

qxtr=(1-Γ)qxL+ΓqxT

(9)

式中：qxtr為轉(zhuǎn)捩熱流密度；qxL為層流熱流密度；qxT為湍流熱流密度；Γ為加權(quán)函數(shù)。

2 熱結(jié)構(gòu)燒蝕碳化熱解模型

2.1 燒蝕表面控制方程

高超聲速飛行器熱防護通常采用燒蝕熱防護方式，其燒蝕防熱機理是在氣動熱環(huán)境作用下，依靠材料自身物理、化學變化吸收氣動加熱，同時燒蝕產(chǎn)物蒸發(fā)氣體和樹脂熱解氣體向邊界層引射，起到對邊界層高溫氣體的阻塞，從而減少了部分氣動加熱。常用熱結(jié)構(gòu)防熱材料為硅基類復合材料。該類材料特點是在一定高溫條件下飛行器表面形成一層熔融的液態(tài)層，其燒蝕表面控制方程包括：

(1)表面能量平衡方程

fcΔhc+αSiO2fSiO2Δhv]

(10)

式中：Ψ為引射因子；ε為材料輻射系數(shù)；σ為斯蒂芬-波爾茲曼常數(shù)；Tw為表面溫度；T0為初始溫度；Cp為材料比熱；mt為質(zhì)量燒蝕速度；fp為樹脂含量；Δhp為樹脂分解熱；fc為碳化分數(shù)；Δhc為碳燃燒熱；αSiO2為二氧化硅蒸發(fā)分數(shù)；fSiO2為二氧化硅質(zhì)量分數(shù)；Δhv為二氧化硅蒸發(fā)熱。左邊第一項為表面氣動加熱，包括質(zhì)量引射對氣動加熱的引射因子Ψ，第二項為表面向環(huán)境輻射熱；右邊第一項為材料初始時刻升至燒蝕溫度時本身吸收熱量，第二項為樹脂熱解吸熱，第三項為碳燃燒熱，第四項為SiO2蒸發(fā)吸熱。

(2)質(zhì)量守恒方程

(11)

式中：B為積分因子；ρ為材料密度；δL為液態(tài)層厚度；P為邊界層外緣壓力；μw為液態(tài)層粘度；τ為剪切力。

(3)引射因子

(12)

式中：層流下β=0.62，湍流下β=0.2；層流下η=0.26，湍流下η=0.33；Ma為空氣分子量；Mj為引射氣體分子量。Ψ反映了質(zhì)量引射時氣動熱流密度與無質(zhì)量引射時氣動熱流密度之比,質(zhì)量燒蝕速度mt可通過迭代方法求得。

2.2 熱結(jié)構(gòu)溫度響應/熱解碳化多元耦合控制方程

高超聲速飛行器熱結(jié)構(gòu)溫度場通常采用一維固體內(nèi)部導熱微分方程計算。求解區(qū)域為除燒蝕層的熱結(jié)構(gòu)內(nèi)部溫度和承力結(jié)構(gòu)內(nèi)壁溫度。常規(guī)固體導熱微分方程求解樹脂基熱結(jié)構(gòu)內(nèi)部溫度響應時不考慮材料物性參數(shù)隨溫度變化和樹脂熱解碳化過程的吸熱，計算的承力結(jié)構(gòu)內(nèi)壁溫度往往偏高。偏差主要原因是氣動熱通過表面?zhèn)飨騼?nèi)部時，防熱材料中的樹脂在高溫作用下會發(fā)生熱解碳化，一般300℃時開始分解，800℃左右分解碳化完畢，樹脂基材料密度由100%降至60%左右。這期間材料熱物理性能發(fā)生了明顯變化，熱物性參數(shù)已不是常溫狀態(tài)下的熱物性參數(shù)，同時材料熱解碳化過程吸收一定熱量，特別是長時間飛行條件下，該作用更為明顯。為合理計算高溫樹脂基熱結(jié)構(gòu)內(nèi)部溫度響應，應在考慮了材料熱解吸熱的常規(guī)導熱微分方程中加入材料熱物性參數(shù)隨溫度變化項。

考慮材料動態(tài)熱解吸熱、材料熱物性參數(shù)隨溫度變化條件下，飛行器熱結(jié)構(gòu)內(nèi)部溫度響應的控制方程為：

(13)

根據(jù)修正的阿侖尼烏斯(S.Arrhenius)熱動力學方程，熱解氣體質(zhì)量生成率為：

(14)

式中：ρ0為原始材料密度；ρc為材料完全碳化的極限密度；K0為頻率因子；E為活化能；R為氣體常數(shù)；n為反應級數(shù)。

材料導熱系數(shù)隨溫度變化關(guān)系為：

λ=a+bT+cT2

(15)

式中：a，b，c為通過試驗回歸計算所得系數(shù)。對式(13)所示的控制方程采用時間和空間有限差分離散，以式(10)燒蝕移動方程求得的Tw為邊界條件，通過迭代方法即可求得考慮物性參數(shù)變化和熱解吸熱的飛行器熱結(jié)構(gòu)內(nèi)部溫度響應和碳化層厚度。

3 計算結(jié)果與討論

3.1 實測表面熱流與理論計算比較

圖1為國內(nèi)首次實測Ma12以上飛行器0°子午線(迎風母線)，X/L=0.53處(X為測點距飛行器原點距離，L為飛行器總長)熱流密度真實變化歷程和理論預示結(jié)果的對比。由圖1可知，實測飛行器表面熱流密度變化趨勢與理論預示熱流密度變化趨勢一致，實測熱流密度略高于理論計算值，結(jié)合飛行器熱流傳感器地面試驗結(jié)果分析認為屬于熱流傳感器系統(tǒng)測量誤差?？偟膩砜蠢碚擃A測熱流密度結(jié)果與實測值偏差較小，兩者偏差小于20%。熱環(huán)境工程設(shè)計允許的誤差在30%。

圖1 熱流密度實測值和計算值比較Fig.1 Comparison between measured and calculated values of the heat flux density

3.2 理論計算轉(zhuǎn)捩熱流與實測比較

飛行器在臨近空間飛行時，隨著高度的降低，Re逐漸增大，當飛行器降低到一定高度，Re達到某個值時，飛行器表面的邊界層出現(xiàn)轉(zhuǎn)捩，流態(tài)由層流逐漸變?yōu)橥牧?。一直以來，邊界層轉(zhuǎn)捩問題在飛行器熱防護設(shè)計中備受關(guān)注，因為邊界層一旦發(fā)生轉(zhuǎn)捩，將對氣動力、氣動熱產(chǎn)生重要影響。圖2給出了實測飛行器0°子午線，X/L=0.53處邊界層轉(zhuǎn)捩真實歷程和理論計算結(jié)果的比較。由圖2可知熱流密度在某一時刻出現(xiàn)了突增，表明在該時刻邊界層出現(xiàn)轉(zhuǎn)捩。理論計算和實測熱流密度的轉(zhuǎn)捩起始時間點基本一致，曲線變化規(guī)律一致且實測值略高于理論值，主要原因是轉(zhuǎn)捩是非常復雜的流態(tài)變化過程，準確計算轉(zhuǎn)捩熱流值目前還很困難，計算方法有待于完善。

圖2 轉(zhuǎn)捩熱流實測值和計算值比較Fig.2 Comparison between the measured and calculated values of the transition heat flux density

3.3 熱結(jié)構(gòu)內(nèi)部分層溫度實測結(jié)果與理論計算比較

圖3、圖4是飛行器熱結(jié)構(gòu)X/L=0.53，0°子午線處距離防熱層表面4 mm、8 mm處理論計算值與實測分層溫度的比較。由圖3、圖4可知，理論計算防熱結(jié)構(gòu)分層溫度響應規(guī)律與遙測結(jié)果一致，距離防熱結(jié)構(gòu)表面越近則溫度變化速率越大，主要是防熱結(jié)構(gòu)表面熱傳遞是非穩(wěn)態(tài)過程，距表面越近則溫度梯度越大，沿厚度方向溫度梯度逐漸減小，表現(xiàn)溫度變化率趨于減緩。從計算值和實測值對比還可知，不考慮樹脂熱解吸熱和材料物性參數(shù)變化，理論計算值明顯高于實測值，圖3最大偏差275℃，圖4是320℃，而考慮熱解吸熱和物性參數(shù)隨溫度的變化，則理論計算值與實測值最大偏差小于70℃，可見在樹脂基常規(guī)導熱微分方程中考慮樹脂熱解吸熱和導熱系數(shù)隨溫度變化因素更符合高溫氣動加熱條件下樹脂基防熱材料物理化學變化過程，計算結(jié)果更為合理。

3.4 熱結(jié)構(gòu)燒蝕碳化實測結(jié)果與理論計算比較

圖5為飛行器0°和45°子午線，X/L=0.53兩個測點處防熱層碳化層厚度實測值和理論計算結(jié)果比較?？芍碚撚嬎愫蛯崪y結(jié)果變化趨勢一致，理論值略高于實測值，0°迎風子午線防熱層燒蝕碳化厚度均高于45°子午線，是由于0°迎風子午線的熱環(huán)境高于45°子午線，熱環(huán)境高的防熱層表面溫度也高，表面溫度越高，燒蝕碳化越嚴重。

圖4 距離表面8 mm層間溫度理論值與實測值比較Fig.4 Comparison between measured and calculated values the stratified temperature for 8 mm from the surface

圖5 熱結(jié)構(gòu)碳化層理論計算和實測結(jié)果比較Fig.5 Measured and calculated values comparison between the thickness of the carbonized layer

4 結(jié)束語

國內(nèi)首次測得了Ma12以上高超聲速飛行器表面熱流密度真實變化歷程和邊界層轉(zhuǎn)捩特征，本文通過對比分析再入飛行器熱流密度/燒蝕碳化/分層溫度實測與理論結(jié)果，得到以下結(jié)論：

1)熱流密度理論計算與實測結(jié)果對比表明，理論計算與實測結(jié)果較為吻合，兩者偏差小于20%。證明目前采用的高超聲速飛行器氣動熱理論預示方法合理可信；

2)研究表明，若不考慮樹脂的動態(tài)熱解特性和導熱系數(shù)隨溫度變化情況，計算的防熱結(jié)構(gòu)內(nèi)部溫度響應遠高于實測值；采用考慮樹脂動態(tài)熱解特性和導熱系數(shù)隨溫度變化的固體導熱微分方程方法計算的防熱結(jié)構(gòu)內(nèi)部溫度響應更為合理，計算方法具有一定的通用性。