分形理論在碳化材料三維燒蝕熱防護(hù)計(jì)算中的應(yīng)用

楊偉斌， 朱慶勇

(中山大學(xué)航空航天學(xué)院,廣東深圳 518107)

引 言

在航空航天領(lǐng)域,燒蝕是高超聲速飛行器最重要的熱防護(hù)形式.當(dāng)燒蝕材料受熱時(shí),材料發(fā)生熱解,釋放出熱解氣體.燒蝕過(guò)程中,熱解氣體進(jìn)入邊界層,降低碳化層熱流,起到了熱阻塞的作用;同時(shí)熱解氣體流經(jīng)碳化層引起溫度的降低,起到冷卻劑的作用.在這類熱防護(hù)系統(tǒng)設(shè)計(jì)過(guò)程中,溫度場(chǎng)的預(yù)測(cè)格外重要.因溫度場(chǎng)預(yù)測(cè)不準(zhǔn)而不得不采用過(guò)大的防熱冗余,將嚴(yán)重影響飛行器性能;而不足的防熱設(shè)計(jì)將有可能導(dǎo)致飛行器損毀.

燒蝕材料通常用作熱防護(hù)材料,可以保護(hù)航天器免受大氣摩擦產(chǎn)生的大量熱量[1].針對(duì)材料的燒蝕機(jī)理以及高超聲速飛行器的燒蝕流場(chǎng)特性,國(guó)內(nèi)外展開了廣泛的研究[2-3].Gokcen等用一維模型計(jì)算了碳化物在燒蝕過(guò)程中的內(nèi)部熱響應(yīng)[4].李鴻權(quán)等對(duì)熱解氣在碳化層內(nèi)流動(dòng)進(jìn)行了研究，發(fā)現(xiàn)熱解氣對(duì)熱防護(hù)結(jié)構(gòu)的溫度分布有很大的影響[5].劉驍?shù)妊芯苛颂蓟療g材料的三維溫度場(chǎng)[6].多孔材料碳化過(guò)程中熱解氣體產(chǎn)生的機(jī)理已被許多研究者廣泛研究[7-8].用質(zhì)量守恒、動(dòng)量守恒和能量守恒方程描述物質(zhì)守恒,可以分析熱解氣體在模型中的壓力分布和溫度分布,氣體在多孔介質(zhì)中的流動(dòng)遵循Darcy定律,朱慶勇等通過(guò)從Darcy定律導(dǎo)出的數(shù)學(xué)模型來(lái)探討流體在多孔介質(zhì)中的擴(kuò)散,建立了毛細(xì)管內(nèi)流體以質(zhì)量和能量平衡形式擴(kuò)散的控制方程[9-10].

考慮到熱解氣體的運(yùn)動(dòng),氣體在碳化物中的滲透率對(duì)燒蝕過(guò)程和熱防護(hù)有著重要的影響,對(duì)于具有自相似性的多孔介質(zhì),可以用分形維數(shù)來(lái)評(píng)價(jià)其滲透率[11-14].分形滲透率模型與3個(gè)參數(shù)密切相關(guān),即多孔介質(zhì)的迂曲度分形維數(shù)、孔徑分形維數(shù)和最大孔徑.Yu等建立了滲透率模型,指出分形分析對(duì)評(píng)價(jià)多孔介質(zhì)的輸運(yùn)性質(zhì)具有重要意義,通過(guò)對(duì)輸運(yùn)性質(zhì)的研究,得到了雙重多孔介質(zhì)中氣體輸運(yùn)模型[13].Li等建立了考慮氣體壓縮性和滑移效應(yīng)的分形致密多孔介質(zhì)中的氣體輸運(yùn)模型[15].王登科等針對(duì)氣體滲流中存在的Klinkenberg效應(yīng)研究了煤體瓦斯中的滲流規(guī)律[16].Cao等采用二元混合格子Boltzmann模型(lattice Boltmann method,LBM)模擬了分?jǐn)?shù)Brown運(yùn)動(dòng)模型產(chǎn)生的隨機(jī)分形多孔介質(zhì)中的擴(kuò)散過(guò)程[17].

滲透率隨燒蝕過(guò)程的變化而變化,為了準(zhǔn)確評(píng)估燒蝕過(guò)程,文章結(jié)合全局坐標(biāo)和局部坐標(biāo)的關(guān)系,推導(dǎo)出了材料的滲透率張量,并將新模型引入燒蝕模型中,得到了模型溫度隨燒蝕過(guò)程的變化,并將溫度分布與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了比較,驗(yàn)證了方法的有效性.模擬計(jì)算了不同孔徑分形維數(shù)和不同迂曲度分形維數(shù)材料三維球頭結(jié)構(gòu)的燒蝕過(guò)程.

1 燒蝕模型

1.1 物理模型

燒蝕結(jié)構(gòu)受到外部熱流的影響,發(fā)生化學(xué)熱解反應(yīng),并同時(shí)受擴(kuò)散因素的影響.當(dāng)溫度較低時(shí),燒蝕速率由化學(xué)動(dòng)力學(xué)(Arrhenius方程)和擴(kuò)散速度共同決定.根據(jù)壁面組元質(zhì)量守恒和相容條件,表面燒蝕率可計(jì)算為[18]

(1)

熱解氣包含O2,O,N2,N,CO,CN,C2N,C1,C2,C3,C4,C5等13種組分.各組分之間反應(yīng)如下[19-20]

(2)

(3)

其中,Mi是組分i的分子質(zhì)量,Mav是氣體平均分子質(zhì)量.燒蝕率可由式 (1)～ (3)計(jì)算得出.

1.2 滲透率模型

圖1為掃描電子顯微鏡下的纖維材料表面形貌.圖像顯示碳化纖維材料的橫截面,經(jīng)過(guò)圖像處理運(yùn)用盒子計(jì)數(shù)法可以得到如圖2所示的分形維數(shù)、最大孔徑和孔隙率等材料屬性,從而得出熱解氣體在碳化層的滲透率[21].

圖1 燒蝕材料的掃描電鏡圖像Fig.1 SEM image of ablative material

在計(jì)算多孔介質(zhì)微孔中的氣體擴(kuò)散時(shí),須考慮氣體的滑移效應(yīng).對(duì)于氣體流動(dòng),可以分為黏滯流、分子流和黏滯-分子流[22].哪種流動(dòng)起主要作用,主要由Knudsen數(shù)決定.

其中,λ為微孔尺寸,l為氣體平均分子自由程.平均分子自由程可以表示為

其中,μ為動(dòng)力黏性系數(shù),P為壓強(qiáng),M為分子質(zhì)量,T為材料溫度.熱解氣的平均分子自由程可以計(jì)算為1.643×10-8m.燒蝕材料的孔徑大約為1 μm,進(jìn)而計(jì)算得出Knudsen數(shù)為0.033,在計(jì)算中要考慮氣體的滑移效應(yīng).

朱慶勇等綜合考慮模型傾斜角度、材料分形特征、氣體滑移邊界條件等,推導(dǎo)出氣體在多孔介質(zhì)中的滲透率張量形式[22]為

其中

在實(shí)際工程當(dāng)中,熱流往往沿著一個(gè)特定的方向,因此可僅考慮氣體在一個(gè)方向上的擴(kuò)散,方程表示為

其中，DpT和DtL分別表示材料模截面上的孔徑分形維數(shù)和長(zhǎng)度方向的迂曲度分形維數(shù).

(a)Tortuosity fractal dimensions

1.3 燒蝕模型[18]

當(dāng)發(fā)生熱解之前,熱傳遞方式主要是材料的熱傳導(dǎo).溫度達(dá)到一定程度,材料熱解為熱解氣和碳化結(jié)構(gòu).熱解氣經(jīng)由碳化層的孔隙擴(kuò)散到邊界層.

通過(guò)上述分析,文章做出以下簡(jiǎn)化和假設(shè)[2]:

(1)在燒蝕過(guò)程中,燒蝕材料由外到內(nèi)可分為碳化層-熱解面-原始材料3層結(jié)構(gòu);

(2)熱解反應(yīng)僅發(fā)生在熱解面上;

(3)熱解氣體與碳化層之間換熱充分,二者無(wú)溫差;

材料的熱解過(guò)程主要發(fā)生在一定的溫度范圍內(nèi),一般在熱解面附近.首先,材料受到外部熱流的影響逐漸被加熱,其中表面溫度最先升高,同時(shí)由于材料自身熱傳導(dǎo)的影響,原始材料層溫度也會(huì)升高.當(dāng)溫度達(dá)到熱解溫度后,材料發(fā)生熱解吸熱,形成熱解氣體和多孔碳化層.同時(shí)由于航天器使用環(huán)境使得邊界層熱流密度一般很大,熱解層很薄,厚度可以忽略不計(jì),近似熱解面.熱解氣通過(guò)多孔碳化層注入熱防護(hù)機(jī)構(gòu)邊界層帶走部分熱量.

圖3 燒蝕結(jié)構(gòu)示意圖Fig.3 Schematic diagram of ablation structure

原始層的傳熱方程

碳化層的傳熱控制方程

燒蝕邊界邊界條件

熱解面邊界條件

其中，qN為熱流密度;ρ表示密度,cp表示比熱,κ表示導(dǎo)熱系數(shù),h表示汽化潛熱.下標(biāo)1,2分別表示原始層和碳化層

熱解面上的熱解氣體質(zhì)量流量由熱解溫度決定,為

其中，下標(biāo)Ep表示熱解面所在單元.

碳化層為多孔介質(zhì),當(dāng)熱解氣通過(guò)碳化層流動(dòng)時(shí)服從氣體連續(xù)方程和Darcy定律[6]

其中，ε是孔隙率,[K]是滲透率張量,μ是熱解氣動(dòng)力黏性系數(shù).碳化層外邊界條件為

P=Pe

其中，Pe為外壁面壓強(qiáng).熱解面邊界條件為

2 計(jì)算結(jié)果和分析

根據(jù)前面闡述的方法進(jìn)行系列數(shù)值計(jì)算,將實(shí)驗(yàn)結(jié)果與參考文獻(xiàn)[24]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比.計(jì)算條件為:熱流密度為279.7 kW/m2,材料厚度為30 mm,燒蝕時(shí)間為600 s.材料物理參數(shù)和初始條件如表1所示，材料迂曲度分形維數(shù)為1.1，孔徑分形維數(shù)為1.8.

表1 材料的物理參數(shù)Table 1 Detailed parameters of calculation model and test material

圖4給出不同位置的溫度變化情況.對(duì)比可以看出,數(shù)值計(jì)算的溫度變化與實(shí)驗(yàn)結(jié)果基本一致.實(shí)驗(yàn)結(jié)束時(shí),近壁點(diǎn)溫度為1 206 K，與實(shí)驗(yàn)溫度 1 236 K 擬合良好.從中可以看出本文模型能夠很好地計(jì)算地?zé)岱雷o(hù)機(jī)構(gòu)的防護(hù)效果.

圖4 燒蝕結(jié)構(gòu)內(nèi)部溫度計(jì)算結(jié)果與實(shí)驗(yàn)對(duì)比Fig.4 Comparison of temperature increase at different thickness locations with test results

建立三維球頭模型計(jì)算材料孔徑分形維數(shù)和迂曲度分形維數(shù)不同對(duì)燒蝕的影響.球頭模型的半徑為1.9 cm.具體模型參見參考文獻(xiàn)[23].

2.1 不同孔徑分形維數(shù)材料的燒蝕

圖5給出了不同孔徑分形維數(shù)材料制成的熱防護(hù)結(jié)構(gòu)燒蝕70 s時(shí)的材料密度分布情況.圖6顯示的是不同孔徑分形維數(shù)的熱防護(hù)機(jī)構(gòu)燒蝕70 s時(shí)的熱解氣密度分布情況.結(jié)合溫度分布可以看出,在球頭位置,溫度最高,相應(yīng)地在密度上反映的是材料密度減少,形成熱解氣,造成球頭位置熱解氣密度較高.70 s 時(shí),球頭位置密度分別是1.719 g/cm3(Dp=1.7)，1.717 g/cm3(Dp=1.8)和1.714 g/cm3(Dp=1.9).可以看出隨著孔徑分形維數(shù)的增加,固體密度變小.相應(yīng)地,孔徑分形維數(shù)為1.7,1.8和1.9的球頭結(jié)構(gòu)在 70 s 時(shí)的氣體密度分別是0.187，0.189，0.192 g/cm3.

(a)Dp=1.7

(a)Dp=1.7

分別計(jì)算3種孔徑分形維數(shù)材料熱防護(hù)層的防護(hù)效果.具體溫度分布如圖7所示:計(jì)算70 s時(shí),孔徑分形維數(shù)1.7,1.8和1.9的材料組成的球頭前點(diǎn)溫度分別是1 901,1 909和1 919 K,而相應(yīng)地?zé)岱雷o(hù)層內(nèi)層的溫度分別是388,393和398 K.

(a)Dp=1.7

從圖8可以明顯看出,模擬計(jì)算的前50 s,材料背面的溫度基本沒有變化,相當(dāng)于初始溫度.50 s 時(shí),3種不同孔徑分形維數(shù)材料熱防護(hù)結(jié)構(gòu)的背面溫度均為300.14 K,50 s僅上升0.14 K.50 s后球頭背面溫度快速上升,到70 s時(shí),3種材料(孔徑分形維數(shù)分別為1.7,1.8和1.9)的溫度已到達(dá)388,393和398 K,分別上升了88,93和98 K.可以明顯看出孔徑分形維數(shù)越小,熱防護(hù)效果越好.

圖8 燒蝕70 s時(shí)不同孔徑分形維數(shù)材料燒蝕結(jié)構(gòu)背面的溫度變化Fig.8 Temperature change in 70 s of ablative material with different pore area fractal dimensions

2.2 不同迂曲度分形維數(shù)材料的燒蝕

圖9給出了不同迂曲度分形維數(shù)的材料制成的球頭燒蝕70 s時(shí)的材料密度分布情況,圖10顯示的是70 s時(shí)的熱解氣體密度分布情況.在球頭位置,材料受到熱流的直接影響,溫度最高,相應(yīng)的材料熱解程度最高,形成更多的熱解氣,造成球頭位置熱解氣密度較高.

(a)Dt=1.05

(a)Dt=1.05

70 s時(shí),3種不同材料球頭位置密度分別是1.707 g/cm3(Dp=1.05),1.717 g/cm3(Dp=1.1)和1.724 g/cm3(Dp=1.15).可以看出隨著迂曲度分形維數(shù)的增加,固體密度變大.相應(yīng)地,3種材料的球頭結(jié)構(gòu)在70 s時(shí)的氣體密度分別是0.200,0.189,0.181 g/cm3.

圖11，12分別為不同迂曲度材料燒蝕結(jié)構(gòu)溫度場(chǎng)分布和背面的溫度變化.從圖11，12可以看出:球頭溫度分布受到燒蝕材料孔隙迂曲度分形維數(shù)的影響.當(dāng)燒蝕材料的迂曲度分形維數(shù)分別是1.05,1.1和1.15時(shí),計(jì)算70 s時(shí)球頭結(jié)構(gòu)的前端溫度分別是1 952,1 909 和1 878 K.相應(yīng)地,球頭背面溫度分別是418,393和378 K.

(a)Dt=1.05

圖12 不同迂曲度分形維數(shù)材料燒蝕結(jié)構(gòu)背面的溫度變化Fig.12 Temperature change in 70 s with different tortuosity fractal dimensions

幾種不同迂曲度的材料在燒蝕前期,熱防護(hù)效果差別不大.在計(jì)算50 s時(shí),球頭背面溫度均為300.1 K,沒有差別.隨著燒蝕的進(jìn)行,差距逐漸出現(xiàn),計(jì)算70 s結(jié)束時(shí),三者(迂曲度分形維數(shù)分別是1.05,1.1和1.15的球頭背面)的溫度分別是418,393和378 K,溫差達(dá)到40 K,出現(xiàn)明顯的差別.

3 結(jié)論

本文在三維球頭燒蝕模型中引入多孔介質(zhì)的分形滲透率張量模型,通過(guò)與實(shí)驗(yàn)對(duì)比驗(yàn)證了模型的可靠性.并分別研究了孔徑分形維數(shù)和迂曲度分形維數(shù)對(duì)熱解氣在三維球頭燒蝕過(guò)程中擴(kuò)散的影響,得出如下結(jié)論:(1)孔徑分形維數(shù)越小的材料熱防護(hù)效果越好.燒蝕70 s時(shí),分形維數(shù)分別為 1.7,1.8和1.9材料制成球頭的內(nèi)部溫度分別為388，393 和398 K.(2)迂曲度分形維數(shù)越大的材料熱防護(hù)效果越好.燒蝕70 s時(shí),分形維數(shù)分別為 1.05,1.1和1.15材料制成球頭的內(nèi)部溫度分別為418,393和378.通過(guò)研究不同燒蝕材料對(duì)結(jié)構(gòu)熱防護(hù)效果的影響,有助于更好地研究發(fā)展航空航天器的熱防護(hù)材料和熱防護(hù)結(jié)構(gòu).