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

（中國(guó)空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所，綿陽(yáng) 621000）

引 言

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

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

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

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

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

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

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

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

圖2 試驗(yàn)和測(cè)試方法示意圖Fig.2 Diagram of experimental set-up

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

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

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

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

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

1.1 氣動(dòng)熱計(jì)算方法

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

1.2 熔融燒蝕模型與計(jì)算方法

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

一維燒蝕過(guò)程的導(dǎo)熱微分方程為[18]

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

初始條件為

邊界條件為

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

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

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

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

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

將式（7）簡(jiǎn)化后，得到在t時(shí)刻小行星表面燒蝕厚度為

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

1.3 熔融剪切燒蝕模型與計(jì)算方法

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

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

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

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

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

空氣在壁溫時(shí)的焓haw表達(dá)式為

其中：h為靜焓；C為質(zhì)量分?jǐn)?shù)。

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

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

其中，?hm=269.55 kJ/kg。

式（15）中其它各量表達(dá)式為

對(duì)式（13）左側(cè)化簡(jiǎn)得到

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

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

由式（23）和（24）可得

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

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

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

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

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

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

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

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

圖6 駐點(diǎn)溫度隨時(shí)間變化Fig.6 Temperature of stagnant point changes with time

圖7 駐點(diǎn)燒蝕速率隨時(shí)間變化Fig.7 Ablation rate of stagnant point changes with time surface

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

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

3 結(jié) 論

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