石衛(wèi)波，黨雷寧，羅 躍，孫海浩，黃 潔

（中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所，綿陽 621000）

引 言

近地小行星撞擊地球是人類生存和發(fā)展的潛在威脅之一。近代以來著名的近地天體撞擊事件有1908年通古斯大爆炸[1]、2003年車里雅賓斯克流星事件等[2]。學術(shù)界已經(jīng)基本認同：6 500萬年前1顆直徑10 km的小行星撞擊現(xiàn)墨西哥尤卡坦半島位置，導致包括恐龍在內(nèi)的全球物種大滅絕，史稱“K-T事件”[3]。

小行星以極高速（12～70 km/s）進入地球大氣層，產(chǎn)生巨大的氣動加熱，引起材料大規(guī)模熔融、汽化和燒蝕，導致嚴重的質(zhì)量損失、外形變化甚至多次解體，產(chǎn)生的沖擊波和熱輻射向地面?zhèn)鞑?，破壞地面設(shè)施和生物[4-5]。小行星極高速進入的燒蝕問題，是小行星進入地球大氣的科學問題之一，是研究進入軌跡、光輻射、解體乃至空中爆炸現(xiàn)象的重要基礎(chǔ)[4]。

小行星材料大致分為3類[6]：鐵隕石、石鐵隕石、石隕石。目前，學術(shù)界主要通過小行星進入事件觀測、地面試驗和理論分析3種手段研究小行星的燒蝕問題。在小行星進入事件觀測中，可通過觀測到的軌跡分析燒蝕系數(shù)[7]。然而由于鐵質(zhì)小行星占比少[8]或缺乏觀測數(shù)據(jù)，鐵質(zhì)小行星的燒蝕系數(shù)數(shù)據(jù)較少。在地面試驗方面，美國國家航空航天局（National Aeronautics and Space Administration，NASA）艾姆斯研究中心（Ames Research Center）2018年在60 MW電弧風洞中開展了鐵質(zhì)小行星材料進入地球大氣層的燒蝕試驗[9]，試驗?zāi)M熱流32.9 MW/m2、滯止壓力126.5 kPa，相當于30 m直徑的小行星在65 km海拔高度以20 km/s的速度飛行，試件材料為IAB-MG型Campo Del Cielo鐵隕石。中國空氣動力研究與發(fā)展中心超高速研究所（簡稱CARDC超高速研究所）2021年在20 MW電弧加熱器上開展了球錐外形鐵隕石模型的燒蝕試驗[10]。理論分析方面的工作大多以石質(zhì)小行星為研究對象，針對鐵質(zhì)小行星燒蝕的工作較少。Girin[11]基于熔融層的失穩(wěn)提出熔融噴濺的燒蝕模型，計算了6 mm直徑鐵質(zhì)微流星的燒蝕，但缺乏與觀測或地面試驗的對比。?apek等[12]提出另一種熔融噴濺的燒蝕理論，?apek等[13]對該理論進行了改進。該理論認為小行星表面熔融層在溫度達到熔點后立即脫落，分散成液滴在大氣中繼續(xù)減速。該理論對質(zhì)量0.006～1.1 g微流星燒蝕的預測結(jié)果與光曲線觀測結(jié)果符合較好。

從研究現(xiàn)狀看，國內(nèi)外對鐵質(zhì)小行星燒蝕問題研究較少，且主要針對微流星，不僅在行星防御范疇之外，建立的燒蝕理論也缺乏地面試驗詳細測量支撐。本文基于CARDC超高速研究所電弧加熱器鐵隕石燒蝕試驗中觀察到的現(xiàn)象，建立了鐵質(zhì)小行星材料燒蝕的熔融層剪切燒蝕模型，并采用具有移動邊界的氣動熱、燒蝕與內(nèi)部熱傳導耦合求解技術(shù)，對燒蝕試驗進行了分析。

1 鐵質(zhì)小行星材料燒蝕模型與計算方法

CARDC超高速研究所于2021年在20 MW電弧加熱器上，針對如圖1所示的球錐外形，開展了3個狀態(tài)（見表1）的燒蝕試驗[10]。電弧加熱器產(chǎn)生的高溫氣體經(jīng)噴管加速形成高焓高速氣流，對試驗?zāi)Ｐ颓岸嗣鏇_擊加熱，通過調(diào)整加熱器狀態(tài)參數(shù)模擬小行星極高速進入大氣層時的駐點熱流、壓力、焓值，研究材料的燒蝕特性，試驗和測試方法示意如圖2所示[14]。支座安裝在試驗艙中的快速送進系統(tǒng)上，前端安裝試驗?zāi)Ｐ?。試驗過程中的參數(shù)測量：模型背面設(shè)置熱電偶測量背溫；高清攝像機記錄流場結(jié)構(gòu)和局部材料瞬時變化過程；比色高溫計測量模型駐點附近溫度；紅外攝像儀獲取模型表面溫度分布；光柵光譜儀測量流場中的發(fā)射光譜；雙目視覺系統(tǒng)在線實時測量模型縱截面輪廓變化。試驗后的參數(shù)測量主要包括燒蝕質(zhì)量、線燒蝕率。試驗視頻截圖（氣流自左向右）如圖3所示，可看到熔融層受剪切向后流動并覆蓋試驗?zāi)Ｐ腿砻妗?/p>

表1 燒蝕試驗?zāi)M狀態(tài)Table 1 Simulation state of ablation test

圖1 鐵隕石燒蝕試驗?zāi)Ｐ虵ig.1 Iron meteorite ablation test model

圖2 試驗和測試方法示意圖Fig.2 Diagram of experimental set-up

圖3 鐵石燒蝕過程的視頻截圖Fig.3 Video screenshot of iron stone ablation process

上述試驗現(xiàn)象得到的啟發(fā)：把熔融、蒸發(fā)和剪切作為鐵質(zhì)小行星材料燒蝕的主要因素，熔融是材料在氣動加熱下由固態(tài)到液態(tài)的相變，蒸發(fā)是由液態(tài)到氣態(tài)的相變，剪切是高速氣流對熔融層的機械剝蝕作用。

為評估熔融、蒸發(fā)和剪切在鐵質(zhì)小行星材料燒蝕中的作用，本文分別建立熔融燒蝕模型和熔融剪切燒蝕模型，模擬電弧射流中的球錐外形鐵隕石燒蝕情況。熔融燒蝕模型假設(shè)表面達到熔點后熔融層立即去除。熔融剪切燒蝕模型則考慮了表面熔融層的蒸發(fā)以及在氣流剪切下的流失。

燒蝕計算需要和氣動熱、材料內(nèi)部熱傳導耦合起來，是一個非定常問題（圖4）。氣動熱是熱傳導計算的邊界條件，內(nèi)部熱傳導結(jié)果是燒蝕計算的輸入條件，燒蝕引起的外形變化對氣動熱產(chǎn)生影響，燒蝕引起的物面后退影響熱傳導物面邊界條件。其中內(nèi)部熱響應(yīng)和材料燒蝕計算采用緊耦合，其它模塊采用松耦合。

圖4 鐵隕石模型氣動熱/燒蝕耦合計算流程圖Fig.4 Flowchart of aerodynamic heat/ablation coupling calculation of iron meteorite model

1.1 氣動熱計算方法

電弧加熱器試驗條件下模型氣動熱計算根據(jù)模型表面不同區(qū)域和邊界層不同流動狀態(tài)選取，駐點熱流采用Fay-Riddell公式[15]；層流熱流采用Zoby等[16]給出的參考焓方法計算簡捷，具有較高的精度；湍流熱流采用布拉休斯平板表面摩阻關(guān)系式、Eckert的參考焓壓縮性修正以及Colburn的雷諾比擬關(guān)系相結(jié)合的方法[17]。

1.2 熔融燒蝕模型與計算方法

對于極高速進入地球大氣層的鐵質(zhì)小行星，氣動加熱作用時間短，表面熱流只影響局部厚度的溫度響應(yīng)，尚未深入到物體內(nèi)部，因此，可將鐵質(zhì)小行星的瞬態(tài)傳熱看作半無限大物體的瞬態(tài)傳熱問題。

一維燒蝕過程的導熱微分方程為[18]

其中：T為r處的溫升；a為介質(zhì)的熱擴散系數(shù)，a=λ/(ρ·c)；t為導熱時間。

初始條件為

邊界條件為

其中：S(t)為移動界面；L為熔解熱（J·kg–1）；傳入熱量qw等于金屬熔化吸收和向小行星內(nèi)部導熱。燒蝕時表面在移動，采用相對動坐標 ξ，ξ=r-S(t)，導熱方程經(jīng)坐標變換為

其中：T是 ξ的函數(shù)，偏微分方程轉(zhuǎn)化為常微分方程。由于假定在小的時間步長內(nèi)小行星表面熱流密度為常數(shù)qw=const，由此可認為小的時間步長內(nèi)燒蝕表面的后退速度也是常數(shù)V-∞=const，即材料在某網(wǎng)格點是逐層燒蝕。

邊界條件也作相應(yīng)變換后得到

由式（5）可解得小行星內(nèi)溫度分布為

其中：a=λ/(ρcp)為熱擴散率；cp為定壓比熱。

將式（7）簡化后，得到在t時刻小行星表面燒蝕厚度為

由式（8）可得金屬材料表面熔融燒蝕速率

1.3 熔融剪切燒蝕模型與計算方法

鐵質(zhì)小行星材料在嚴酷熱環(huán)境下，一方面與氧發(fā)生燃燒，另一方面在氣動力/熱作用下，達到熔點熔化向后流動。在熔融剪切模型中，假設(shè)燒蝕表面存在一層薄的液體層，所有化學反應(yīng)都發(fā)生在液體層表面，且達到平衡。下面建立液體層表面能量守恒、質(zhì)量守恒和動量守恒關(guān)系式。通過求解這些關(guān)系式，即可獲得燒蝕特性。

在固–液交界面會發(fā)生如下熔化反應(yīng)

在液體層外表面，存在液態(tài)層的蒸發(fā)和燃燒反應(yīng)

考慮對流換熱、輻射傳熱、反應(yīng)熱與燒蝕帶走的熱量，則物面的能量平衡方程為

其中：haw為空氣在壁溫時的焓；hr為恢復焓；qr為輻射熱；?Q為反應(yīng)熱；qc為冷壁熱流；CT為熔化前材料的比熱；ρT為材料密度；V-∞為燒蝕速率；?H-∞為單位質(zhì)量熔化潛熱；ks為固體熱傳導系數(shù)；Tw為物面溫度。

空氣在壁溫時的焓haw表達式為

其中：h為靜焓；C為質(zhì)量分數(shù)。

反應(yīng)熱 ?Q用下式計算

其中：?QFe在計算中取正值。熔化潛熱?H-∞的表達式為

其中，?hm=269.55 kJ/kg。

式（15）中其它各量表達式為

對式（13）左側(cè)化簡得到

假定邊界層內(nèi)化學凍結(jié)，化學反應(yīng)在壁面進行，且假定普朗特數(shù)和Lewis數(shù)Le=1情況下，由物面的質(zhì)量守恒關(guān)系的邊界層擴散方程可得到

設(shè)式（22）中KFe,e=0，并除以ρeueCH有

由式（23）和（24）可得

在液體層表面有如下動量守恒關(guān)系

其中：ρT為固體材料密度；ρL為液體層密度；δ為液體層厚度；τw為剪切應(yīng)力；rb為物體徑向坐標。

式（13）、（25）、（27）～（28）是聯(lián)系4個未知量B-∞、TW、αFe、Mˉ的超越方程組，可用迭代法求解。其中B-∞可以從式（13）和B-∞的定義得到

聯(lián)立方程組求解得，然后給出燒蝕速率V-∞

2 鐵隕石燒蝕計算結(jié)果分析

鐵隕石燒蝕試驗是在CARDC超高速研究所20 MW電弧加熱器設(shè)備開展。為便于鐵隕石燒蝕理論研究與試驗測量結(jié)果對比分析，根據(jù)模擬的參數(shù)設(shè)計成球錐模型，頭部半徑R=20 mm，半錐角9°，高度70 mm（圖1）。鐵隕石熱物性參數(shù)為[19]：密度7 827.234 kg/m3，熔點1 846 K，比熱364 ± 19 J/(kg·K)，42.9 ± 15.5 W/(m·K)，輻射系數(shù)0.671，熔化潛熱269.55 kJ/kg。鐵隕石模型燒蝕試驗狀態(tài)參數(shù)有3個（表1），試驗時間4 s。

根據(jù)1.2節(jié)建立的鐵隕石熔融燒蝕計算模型，本文分別對表1的3個試驗狀態(tài)進行了計算。第Ⅰ狀態(tài)駐點平均線燒蝕計算速率1.55 mm/s，第Ⅱ狀態(tài)駐點平均線燒蝕計算速率2.14 mm/s，第Ⅲ狀態(tài)駐點平均線燒蝕計算速率2.50 mm/s。燒蝕開始的時間隨熱流的增加提前，得到的駐點燒蝕速率也隨駐點熱流的增加而增大，而地面試驗測得的駐點燒蝕速率隨駐點壓力的增加而增大（圖5，表2）。計算與試驗在燒蝕速率差異的原因是鐵隕石燒蝕計算模型為熔融模型，假定表面材料溫度達到熔點就發(fā)生質(zhì)量損失，沒有考慮鐵隕石試驗高溫熔融流動剪切作用。

根據(jù)1.3節(jié)建立的鐵隕石熔融剪切燒蝕計算模型，本文對鐵隕石球錐試件（圖1）試驗狀態(tài)開展計算分析，計算與試驗測量結(jié)果比較見表2。鐵隕石球錐試件燒蝕試驗的計算結(jié)果如圖6～9所示。圖6～7橫坐標為時間，圖8、圖9橫坐標為沿物面弧長（從駐點算起），圖6縱坐標為駐點溫度，圖7、圖9縱坐標為燒蝕速率，圖8縱坐標為液態(tài)層厚度。試驗狀態(tài)Ⅰ～Ⅲ的試驗條件是熱流遞增、駐點壓力遞減。熱流越高，導致溫度越高（圖6）、越早燒蝕，狀態(tài)Ⅰ～Ⅲ開始燒蝕的時間分別為3、1.5和1 s（圖7）。穩(wěn)定燒蝕后，駐點燒蝕速率隨試驗駐點壓力增加而增大（圖7）。液體層的黏性系數(shù)隨壁溫的增加而降低，因此熱流和溫度越高（從狀態(tài)Ⅰ～Ⅲ），液體層越?。▓D8）。從液體層厚度沿物面的分布看，在頭部區(qū)域附近較厚，在身部較薄，與試驗觀察的現(xiàn)象一致（圖3）。頭部高熱流區(qū)域，鐵隕石材料在高溫下熔融并不斷被氣流吹走，部分覆蓋在身部低熱流區(qū)，而模型身部大部分區(qū)域熱流低，隕石材料還未發(fā)生熔融燒蝕（圖9）。

圖6 駐點溫度隨時間變化Fig.6 Temperature of stagnant point changes with time

圖7 駐點燒蝕速率隨時間變化Fig.7 Ablation rate of stagnant point changes with time surface

圖9 燒蝕速率隨物面變化Fig.9 Ablation rate changes with object surface

計算與試驗對比如表2所示，表2中包含了熔融剪切燒蝕模型與不考慮剪切的熔融燒蝕模型計算結(jié)果。在熔融剪切燒蝕模型下，狀態(tài)Ⅰ～Ⅲ的駐點線燒蝕速率分別為1.91、1.86和 1.80 mm/s，隨試驗狀態(tài)的變化規(guī)律在定性上與試驗相符。而熔融模型得到的駐點燒蝕速率與試驗相反，且熱流越高，計算結(jié)果與試驗測量結(jié)果相差越大。這說明熔融層的剪切是燒蝕模型中必須考慮的因素。由表2還可看到，表面蒸發(fā)速率與液態(tài)層流失速率相比是小量，說明熔融層的剪切流失是燒蝕的主要機理。

3 結(jié) 論

本文針對鐵質(zhì)小行星以極高速進入地球大氣層，在嚴酷氣動加熱作用下發(fā)生高溫熔融燒蝕問題，分別建立了熔融燒蝕計算模型和熔融剪切燒蝕計算模型。通過對電弧加熱器試驗條件下鐵隕石球錐模型試驗狀態(tài)的計算分析，獲得了鐵隕石高溫熔融層剪切流失的燒蝕機理。下一步，將根據(jù)建立的鐵質(zhì)小行星進入大氣層高溫熔融剪切燒蝕模型，評估以極高速進入地球大氣層鐵質(zhì)小行星熔融燒蝕情況，為小行星撞擊地球防御提供參考。