飛行器氣動(dòng)加熱燒蝕工程計(jì)算

張志豪，孫得川

(大連理工大學(xué)航空航天學(xué)院，工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024)

0 引言

飛行器在高速飛行時(shí)，會(huì)產(chǎn)生較大的氣動(dòng)熱。高溫氣體在飛行器飛行過(guò)程中持續(xù)向機(jī)體內(nèi)部傳熱，從而使機(jī)體尤其是頭錐、翼梢前緣等部位溫度升高，進(jìn)而影響結(jié)構(gòu)強(qiáng)度，造成設(shè)備失效等不良后果。對(duì)于氣動(dòng)加熱強(qiáng)烈的部位，常采用燒蝕材料進(jìn)行熱防護(hù)。熱防護(hù)材料通常由多層材料組成，其防、隔熱性能直接關(guān)系到能否完成既定任務(wù)，因此，在設(shè)計(jì)階段，預(yù)估氣動(dòng)熱環(huán)境和熱防護(hù)材料性能是必不可少的環(huán)節(jié)。

要實(shí)現(xiàn)對(duì)飛行器熱防護(hù)材料溫度場(chǎng)及燒蝕情況的快速預(yù)估，需要實(shí)現(xiàn)在一定效率和精度下的氣動(dòng)熱、燒蝕、多層材料溫度場(chǎng)的耦合計(jì)算。國(guó)內(nèi)已有研究人員對(duì)氣動(dòng)熱［1-4］、燒蝕機(jī)理［5-8］等方面展開了不同程度的研究，但仍然缺乏一套高效并具有一定精度的，能夠?qū)崿F(xiàn)氣動(dòng)熱環(huán)境下飛行器熱防護(hù)層性能快速預(yù)報(bào)及輔助設(shè)計(jì)的方法。

氣動(dòng)熱計(jì)算主要有兩類方法:一是直接對(duì)N-S方程進(jìn)行數(shù)值計(jì)算，二是采用工程算法。基于求解N-S 方程的氣動(dòng)熱算法雖然精度較高，但計(jì)算效率低，且對(duì)設(shè)計(jì)初期，氣動(dòng)外形還未完全確定的情況下，該方法并不適用。而工程算法具有較高的可信度和效率，只需給定特征外形和飛行參數(shù)即可預(yù)估氣動(dòng)熱。為此，本文將飛行彈道作為輸入?yún)?shù)，集成氣動(dòng)熱計(jì)算的工程算法(也保留通過(guò)計(jì)算流體力學(xué)或?qū)嶒?yàn)得到的熱流數(shù)據(jù)輸入接口)、材料燒蝕模型和傳熱計(jì)算，實(shí)現(xiàn)了具有較高效率和精度的氣動(dòng)熱環(huán)境下熱防護(hù)材料燒蝕傳熱的計(jì)算方法和軟件。可對(duì)輸入的不同飛行參數(shù)和多層材料組合進(jìn)行快速評(píng)估，達(dá)到輔助設(shè)計(jì)的目的。

1 計(jì)算方法

以飛行彈道為輸入?yún)?shù)的氣動(dòng)熱燒蝕計(jì)算主要包括3 部分:氣動(dòng)熱計(jì)算、材料燒蝕計(jì)算、多層材料瞬態(tài)溫度場(chǎng)計(jì)算。

1.1 氣動(dòng)熱工程計(jì)算

目前常用的高超聲速氣動(dòng)熱工程算法綜合考慮了外流場(chǎng)和實(shí)驗(yàn)數(shù)據(jù)，不僅計(jì)算效率高，而且對(duì)于飛行器的典型幾何特征位置(如鈍頭旋成體，翼梢前緣等部位)的駐點(diǎn)、球頭表面和錐體表面等具有較高的可信度。

駐點(diǎn)熱流的工程計(jì)算方法采用Fay-Riddell 公式［9］。該公式是依據(jù)相似性假定，對(duì)高溫氣體的邊界層方程進(jìn)行簡(jiǎn)化得到的，適用于來(lái)流總焓hs在1 549 ～24 158 kJ/kg 之間，壁溫Tw在300 ～3 000 K之間的計(jì)算。

式中:qws為駐點(diǎn)熱流;Pr 為普朗特?cái)?shù);ρw為壁面密度;ρs為駐點(diǎn)密度;μs為駐點(diǎn)粘性系數(shù);μw為壁面粘性系數(shù);ue為邊界層速度;hD為空氣解離焓;hw為壁面焓;Le 為路易斯數(shù)。

球頭表面熱流密度按照Lees 簡(jiǎn)化公式［10］進(jìn)行歸一化計(jì)算。

式中:qwl為球頭表面熱流;θ 為球頭圓心角;γ∞為來(lái)流比熱比;Ma∞為來(lái)流馬赫數(shù);

錐體熱流密度qwr同樣通過(guò)歸一化計(jì)算，且僅與駐點(diǎn)熱流qws和來(lái)流參數(shù)有關(guān)。

式中:θc為錐角;為當(dāng)?shù)貜楏w半徑;Rn為球頭半徑;

在給定飛行彈道后，通過(guò)飛行馬赫數(shù)，攻角和大氣參數(shù)計(jì)算出來(lái)流參數(shù)，即可快速計(jì)算出特征部位的氣動(dòng)熱。

1.2 燒蝕計(jì)算

目前飛行器常用燒蝕防熱材料有硅基材料和碳基材料兩種。不同類型的燒蝕材料其燒蝕模型略有不同，但本質(zhì)仍然是依據(jù)質(zhì)量守恒和能量守恒。硅基材料的燒蝕模型如下。

本文算例采用的高硅氧/酚醛材料是典型的硅基材料。其燒蝕機(jī)理為:當(dāng)硅基材料開始燒蝕后，表面不斷向后退移，如圖1所示。原始材料在氣動(dòng)加熱下，溫度逐漸升高，當(dāng)溫度達(dá)到熱解溫度后，材料中的樹脂成分開始熱解、蒸發(fā)并進(jìn)入氣體邊界層;隨著樹脂材料的熱解和蒸發(fā)，剩余材料的溫度在熱流的作用下持續(xù)升高，到達(dá)某一溫度后開始碳化，形成碳化層。最靠近外部的壁面溫度最高，材料中的SiO2等成分處于液體態(tài)，并形成液態(tài)層。在發(fā)生燒蝕的情況下，高速氣流會(huì)對(duì)材料表面產(chǎn)生很高的剪切力，從而導(dǎo)致材料燒蝕表面被不斷吹除，使燒蝕材料表面不斷向內(nèi)退移。通過(guò)材料的熱容、熔化、蒸發(fā)、吸熱反應(yīng)、引射、熱阻塞效應(yīng)等吸熱或隔熱。

圖1 高硅氧/酚醛復(fù)合材料燒蝕模型［11］Fig.1 Ablation model of silica/phenolic composite［11］

硅基材料燒蝕機(jī)理主要由質(zhì)量損失機(jī)理和吸熱機(jī)理兩部分構(gòu)成。質(zhì)量損失機(jī)理是指:隨著溫度升高，材料中的樹脂成分首先發(fā)生熱解反應(yīng)，生成C和氣體，表面的C 與O2反應(yīng)燃燒，生成CO2等氣體帶走質(zhì)量;隨著溫度繼續(xù)升高，材料中的高硅氧纖維開始軟化、熔融，在材料表面形成一層薄SiO2液態(tài)層，一部分蒸發(fā)氣化被帶走，另一部分在氣動(dòng)剪切力的作用下被吹除。吸熱機(jī)理是指:材料熱容吸熱;樹脂成分熱解吸熱;高硅氧纖維熔化吸熱;SiO2液態(tài)層蒸發(fā)吸熱、熱解氣體和蒸發(fā)氣體引射引起的熱阻塞效應(yīng);C 燃燒放熱;材料表面熱輻射;材料被氣動(dòng)力吹除帶走熱量。

根據(jù)上述兩個(gè)機(jī)理得到燒蝕計(jì)算的能量守恒和質(zhì)量守恒方程

式中:ψ 為引射因子;δL為表面液態(tài)層厚度;qw為熱流;qrs為空氣離解熱;qe為輻射換熱;ρ 為材料密度;v-∞為燒蝕速率;Cp為材料熱容;Tp為樹脂碳化溫度;fp為樹脂質(zhì)量分?jǐn)?shù);fc為碳化分?jǐn)?shù);Δhp為樹脂蒸發(fā)熱;Δhc為碳化熱;Δhv為SiO2蒸發(fā)熱;τ 為剪切力;p 為壓力。

式中:λL為液態(tài)層導(dǎo)熱系數(shù);Mj為樹脂分子量;Mair為空氣分子量。

用迭代法對(duì)(7)式～(10)式進(jìn)行計(jì)算，得到燒蝕溫度和燒蝕速率，然后以Tw作為溫度邊界條件，v-∞作為邊界運(yùn)動(dòng)速度，根據(jù)材料內(nèi)部的傳熱計(jì)算，就可計(jì)算出瞬態(tài)溫度場(chǎng)。

1.3 多層材料溫度場(chǎng)計(jì)算

對(duì)于大面積防熱來(lái)說(shuō)，沿物面方向溫度梯度遠(yuǎn)小于厚度方向的溫度梯度，因此采用一維傳熱可獲得局部的傳熱情況。溫度場(chǎng)計(jì)算控制方程為一維非穩(wěn)態(tài)導(dǎo)熱方程:

其數(shù)值求解過(guò)程如下:首先是對(duì)求解區(qū)域進(jìn)行離散化，即對(duì)于m 層材料組成的組合壁面，將每層材料分別任意等分為N(i)，i=1，2，…，m 個(gè)節(jié)點(diǎn)，如圖2所示，然后對(duì)每個(gè)節(jié)點(diǎn)列出對(duì)應(yīng)的差分方程，最后聯(lián)立求解各個(gè)節(jié)點(diǎn)的方程組，就可以得到各點(diǎn)的瞬時(shí)溫度。

圖2 多層材料離散Fig.2 Difference of multi-layer heat transfer

方程求解采用FTCS 隱式差分格式，該格式構(gòu)造簡(jiǎn)單，計(jì)算效率高，且無(wú)條件穩(wěn)定。

對(duì)于任意一層內(nèi)部材料i，內(nèi)部節(jié)點(diǎn)j，差分方程寫為

外邊界條件包括熱流、輻射換熱和空氣電離產(chǎn)生的熱流。

需要注意的是，當(dāng)j 位于兩層材料的交界面時(shí)，控制體包含兩種材料，其差分方程為

2 計(jì)算驗(yàn)證

按照上述公式及數(shù)值方法編寫了相應(yīng)的耦合計(jì)算軟件，該軟件可以在輸入彈道(或熱流)后，計(jì)算飛行器特征部位任意指定位置的熱流、多層材料溫度場(chǎng)以及燒蝕率，并可實(shí)時(shí)修改防熱層的層數(shù)、厚度等相關(guān)參數(shù)，便于優(yōu)化設(shè)計(jì)。

首先將數(shù)值計(jì)算值與實(shí)驗(yàn)值進(jìn)行對(duì)比，驗(yàn)證計(jì)算的準(zhǔn)確性，再通過(guò)某導(dǎo)彈的假想飛行彈道，給出本文所述方法的具體應(yīng)用。

2.1 熱流計(jì)算

首先對(duì)飛行器3 個(gè)不同飛行高度H 和不同速度v 條件下的熱流進(jìn)行計(jì)算，并與參考實(shí)驗(yàn)數(shù)據(jù)［12］進(jìn)行對(duì)比。圖3給出了高度7.6 km、21.4 km 和36.6 km 時(shí)對(duì)應(yīng)不同速度的駐點(diǎn)熱流，可見本文的計(jì)算與實(shí)驗(yàn)數(shù)據(jù)符合較好。

圖3 駐點(diǎn)熱流密度Fig.3 Heat flux of stagnation point

圖4給出了彈頭球頭半徑Rn=6.6 mm 的球頭熱流密度分布，在相同狀態(tài)下與文獻(xiàn)［13］的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，可見當(dāng)θ ＜65°符合較好，而在θ ＞70°后，計(jì)算值誤差較大。

圖4 球頭熱流密度Fig.4 Heat flux of bulb

2.2 燒蝕模型驗(yàn)證

燒蝕模型驗(yàn)證以實(shí)驗(yàn)數(shù)據(jù)作為對(duì)比，計(jì)算中，熱流、恢復(fù)焓、參考?jí)毫Φ扰c文獻(xiàn)［14］的實(shí)驗(yàn)一致，高硅氧/酚醛部分物性參數(shù)取值于文獻(xiàn)［15-18］. 表1對(duì)比了計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)，其中材料燒蝕率計(jì)算值與實(shí)驗(yàn)值的平均誤差為3.6%，說(shuō)明計(jì)算和實(shí)驗(yàn)符合較好。

表1 燒蝕數(shù)值計(jì)算與電弧加熱實(shí)驗(yàn)對(duì)比Tab.1 Comparison of calculated results and arc heater experimental results of the ablation properties

3 應(yīng)用算例

以某導(dǎo)彈的假想飛行彈道計(jì)算氣動(dòng)熱和熱防護(hù)材料在熱環(huán)境下的燒蝕情況和溫度場(chǎng)分布。導(dǎo)彈的初速v0=2 km/s，初始飛行高度為H =40 km，終止飛行高度H=0 km，飛行時(shí)長(zhǎng)37 s. 在飛行過(guò)程中，俯仰角和攻角變化不大，俯仰角θ =30°，攻角α =0°. 圖5所示為飛行速度和高度，由于飛行速度高，飛行高度基本呈線性變化。

圖5 飛行速度和高度Fig.5 Flight velocity and altitude

導(dǎo)彈為鈍頭體，球頭半徑Rn=6.6 mm，半錐角θc=20°. 彈頭表面熱防護(hù)材料為4 mm 高硅氧/酚醛復(fù)合材料，彈體材料為鉻鋼，材料物性參數(shù)見表2，初溫25 ℃. 計(jì)算中，大氣參數(shù)根據(jù)飛行高度由標(biāo)準(zhǔn)大氣模型確定。

表2 材料物性參數(shù)Tab.2 Physical parameters and thickness of thermal protection layers

通過(guò)本文集成的計(jì)算方法和程序，計(jì)算了彈頭指定位置的熱流、熱防護(hù)材料燒蝕厚度和溫度場(chǎng)分布隨時(shí)間的變化。以下給出駐點(diǎn)、球頭θ =30°和錐段(當(dāng)?shù)貜楏w半徑)R1=10 mm 處的計(jì)算結(jié)果。

圖6為以上3 個(gè)部位沿飛行彈道的熱流密度，可見隨著飛行速度的加快和高度下降(大氣壓力及密度升高)，熱流密度也隨之逐漸升高。飛行速度在飛行28 s 左右達(dá)最大值，隨后由于阻力升高，速度迅速下降，熱流密度也隨之快速下降。

圖6 計(jì)算部位熱流密度Fig.6 Heat fluxes at calculated locations

圖7為對(duì)應(yīng)位置的燒蝕情況。駐點(diǎn)的燒蝕厚度δ 約為1. 54 mm，球頭θ = 30°處燒蝕厚度為0. 65 mm，錐段燒蝕厚度約為0.2 mm. 由圖7可知，在飛行的初始階段，由于飛行高度較高，空氣稀薄，熱流強(qiáng)度不足以使表面材料升高到燒蝕溫度(約1 400 ℃);從20 s 左右開始，熱防護(hù)材料發(fā)生燒蝕，隨著時(shí)間推移，燒蝕厚度不斷增加，在27 s 左右燒蝕率最大;隨著飛行速度下降，熱流密度減小，燒蝕率又逐漸降低，在35 s 左右，錐段燒蝕基本停止，隨后駐點(diǎn)和球頭的燒蝕也逐漸停止。對(duì)應(yīng)熱流密度可知，燒蝕發(fā)生在熱流密度約4 500 kW/m2以上的熱環(huán)境下。

圖7 計(jì)算部位燒蝕厚度Fig.7 Ablation thicknesses at calculated locations

圖8為導(dǎo)彈飛行37 s 后，彈頭的溫度分布。其中，駐點(diǎn)溫度最高，約1 700 ℃，球頭30°處約1 600 ℃，錐段則為1 500 ℃左右。因?yàn)閺楏w材料為鉻鋼，其導(dǎo)熱系數(shù)較大，在厚度方向上基本為等溫分布，當(dāng)熱防護(hù)層的初始厚度為4 mm 時(shí)，彈體內(nèi)壁溫度升高到約300 ℃，滿足不了熱防護(hù)的要求。因此，將熱防護(hù)材料厚度增加到6 mm，重新進(jìn)行計(jì)算，得到同一時(shí)刻的彈體內(nèi)壁溫度約為170 ℃，可滿足設(shè)計(jì)要求。

圖8 計(jì)算部位溫度場(chǎng)分布Fig.8 Temperature fileds at calculated locations

4 結(jié)論

本文集成了氣動(dòng)熱工程算法、材料燒蝕計(jì)算和傳熱計(jì)算，形成了一種飛行器氣動(dòng)熱及多層防熱材料燒蝕溫度場(chǎng)預(yù)報(bào)的方法。通過(guò)算例驗(yàn)證了該方法具有較高的效率和可信度。采用該方法可以快速進(jìn)行給定彈道條件下的氣動(dòng)熱、熱防護(hù)材料燒蝕性能和彈體溫度場(chǎng)計(jì)算，為飛行器熱防護(hù)層設(shè)計(jì)提供依據(jù)。

References)

［1］呂麗麗.高超聲速氣動(dòng)熱工程算法研究［D］. 西安:西北工業(yè)大學(xué)，2005.LYU Li-li. Research for aerodynamic heat methods of hypersonic［D］.Xi’an:Northwestern Polytechnical University，2005.(in Chinese)

［2］楊凱，高效偉. 高超聲速氣動(dòng)熱環(huán)境工程算法［J］. 導(dǎo)彈與載人航天運(yùn)載技術(shù)，2010(4):19 -23.YANG Kai，GAO Xiao-wei. Engineering algorithm for aeroheating environment of hypersonic aircrafts［J］. Missiles and Space Vehicles，2010(4):19 -23.(in Chinese)

［3］薛鵬飛，龔春林，谷良賢.基于軸對(duì)稱比擬的氣動(dòng)熱計(jì)算方法研究［J］. 計(jì)算機(jī)仿真，2012，29(6):97 -100.XUE Peng-fei，GONG Chun-lin，GU Liang-xian. Study of method for aerothermal heating calculating based on axisymmetric analogue［J］.Computer Simulation，2012，29(6):97 -100.(in Chinese)

［4］李建林，唐乾剛，霍林，等.復(fù)雜外形高超聲速飛行器氣動(dòng)熱快速工程估算［J］. 國(guó)防科技大學(xué)學(xué)報(bào)，2012，34(6):89 -93.LI Jian-lin，TANG Qian-gang，HUO Lin，et al. The rapid engineering aero-heating calculation method for complex shaped hypersonic vehicles［J］. Journal of National University of Defense Technology，2012，34(6):89 -93.(in Chinese)

［5］歐東斌，陳連忠，陳海群，等. 硅基復(fù)合材料燒蝕特性試驗(yàn)研究［J］.宇航材料工藝，2009(1):85 -90.OU Dong-bing，CHEN Lian-zhong，CHEN Hai-qun，et al. Ablation properties of silicon based coposites［J］. Aerospace Materials＆ Technology，2009(1):85 -90.(in Chinese)

［6］孫冰，林小樹，劉小勇，等. 硅基材料燒蝕模型研究［J］. 宇航學(xué)報(bào)，2003，24(5):282 -286.SUN Bing，LIN Xiao-shu，LIU Xiao-yong，et al. Study on ablation model of silica-reinforced composites［J］. Journal of Astronautics，2003，24(5):282 -286.(in Chinese)

［7］張俊華，李錦文，魏化震，等. 高性能酚醛樹脂基燒蝕復(fù)合材料的研究［J］. 纖維復(fù)合材料，2009(1):15 -18.ZHANG Jun-hua，LI Jin-wen，WEI Hua-zhen，et al. Research on high performance phenolic resin ablative composites［J］. Fiber Composites，2009(1):15 -18.(in Chinese)

［8］易法軍，孟松鶴，梁軍，等. 防熱復(fù)合材料高溫體積燒蝕模型［J］. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2001.33(6):725 -728.YI Fa-jun，MENG Song-he，LIANG Jun，et al. A high-temperature volumetric ablation model for thermal protective composites［J］.Journal of Harbin Institute of Technology，2001，33 (6):725 -728.(in Chinese)

［9］Fay J A，Riddell F R. Theory of stagnation point heat transfer in dissociated air［J］. Journal of Aeronautic Science，1958，25(2):73 -85.

［10］Lees L. Laminar heat transfer over blunt nosed bodies at hypersonic flight speed［J］. Jet Propulsion，1956，26(4):259 -269.

［11］時(shí)圣波. 高硅氧/酚醛復(fù)合材料的燒蝕性能與熱力學(xué)行為研究［D］.哈爾濱:哈爾濱工業(yè)大學(xué)，2013.SHI Sheng-bo. Ablation properties and thermomechanical behavior of silica/phenolic composites［D］.Harbin:Harbin Institute of Technology，2013.(in Chinese)

［12］Rose P H，Stark W I. Stagnation point heat transfer measurements in dissociated air［J］. Journal of Aerospace Science，1958，25(2):86 -97.

［13］Kemp N H，Rose R H，Detra R W. Laminar heat transfer around blunt bodies in dissociated air［J］. Journal of the Aerospace Science，1959，26(7):421 -430.

［14］姜貴慶，李仲平，俞繼軍. 硅基復(fù)合材料粘性系數(shù)的參數(shù)辨識(shí)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2008，26(4):425 -455.JIANG Gui-qing，LI Zhong-ping，YU Ji-jun. Prarmeter identification of silicon based composites＇viscous coefficient［J］. Acta Aerodynamica Sinica，2008，26(4):425 -455.(in Chinese)

［15］姜貴慶，劉連元. 高速氣流傳熱與燒蝕熱防護(hù)［M］. 北京:國(guó)防工業(yè)出版社，2003:52 -66.JIANG Gui-qing，LIU Lian-yuan.Heat transfer of hypersonic gas and ablation thermal protection［M］. Beijing:National Defense Industry Press，2003:52 -66.(in Chinese)

［16］張志成，潘海林，劉初平，等.高超聲速氣動(dòng)熱和熱防護(hù)［M］.北京:國(guó)防工業(yè)出版社，2003:222 -248.ZHANG Zhi-cheng，PAN Hai-lin，LIU Chu-ping，et al. Hypersonic aerodynamic heat and thermal protection［M］. Beijing:National Defense Industry Press，2003:222 -248.(in Chinese)

［17］Henderson J B，Verma Y P，Tant M R，et al. Measurement of the thermal conductivity of polymer composites to high temperatures using the line source technique［J］. Polymer Composites，1983，4(4):219 -224.

［18］Henderson J B，Tant M R. A study of the kinetics of high-temperature carbon-silica reactions in an ablative polymer composite［J］. Polymer Composites，1983，4(4):233 -237.

［19］Bruneton E，Tallaron C，Gras-Naulin N，et al. Evolution of the structure and mechanical behaviour of a carbon foam at very high temperatures［J］. Carbon，2002，40(11):1919 -1927.