耿淑娟,周炳紅,韓鵬,鄭偉,李明濤

(1.中國(guó)科學(xué)院國(guó)家空間科學(xué)中心,北京 100190;2.中國(guó)科學(xué)院大學(xué)天文與空間科學(xué)學(xué)院,北京 100049)

1 引言

小天體在撞擊地球過(guò)程中,最先與地球大氣層接觸。在進(jìn)入大氣層的過(guò)程中,小天體在氣動(dòng)力、氣動(dòng)熱的作用下,能量不斷減少,產(chǎn)生燒蝕、破碎現(xiàn)象,作用劇烈時(shí)將在空中發(fā)生爆炸,該現(xiàn)象稱為空爆。探究空爆背后的物理原因,探索小天體的解體機(jī)制,對(duì)小天體撞擊地球的災(zāi)害評(píng)估具有重要意義。

由于小天體進(jìn)入大氣的過(guò)程是一個(gè)復(fù)雜的多因素耦合的過(guò)程,況且地球上存有觀測(cè)數(shù)據(jù)的小天體撞擊事件十分有限,加之小天體軌道數(shù)據(jù)不夠完善,使得人們對(duì)小天體空爆的物理機(jī)制仍未形成較為深刻的認(rèn)識(shí)。

國(guó)內(nèi),柳森等人首先開(kāi)展了小天體撞擊地球大氣層的相關(guān)研究,對(duì)小天體撞擊地球大氣層的物理問(wèn)題及國(guó)內(nèi)外發(fā)展現(xiàn)狀、模型方法進(jìn)行了綜述[1]。在國(guó)外,Melosh、Collns等人提出了一系列小天體進(jìn)入大氣模型,并對(duì)該過(guò)程中的其他效應(yīng)進(jìn)行了研究[2-4]。科學(xué)家們發(fā)展了幾種不同類型的模型[5,6],對(duì)小天體撞擊地球大氣過(guò)程進(jìn)行了數(shù)值模擬,在模擬空爆高度方面取得較大進(jìn)展。但對(duì)于能量沉積過(guò)程、空爆機(jī)理的研究仍未取得重大突破。

模型主要包括煎餅 (Pancake)模型[7]、碎片云 (FCM)模型[8]以及高保真爆炸流體動(dòng)力學(xué)模型[9,10]。

煎餅?zāi)Ｐ蛯⑿√祗w假設(shè)為一個(gè)均質(zhì)的連續(xù)體,并將小天體破碎、在空中發(fā)生爆炸的過(guò)程假設(shè)為一個(gè)連續(xù)的變形過(guò)程。該模型認(rèn)為小天體在進(jìn)入大氣過(guò)程中,當(dāng)其觸發(fā)特定的飛行條件時(shí),開(kāi)始發(fā)生“破碎”,破碎表現(xiàn)為小天體受到擠壓發(fā)生變形。[7]此類模型有明確的空爆高度的概念,當(dāng)小天體截面半徑增大到一定程度時(shí),認(rèn)為小天體發(fā)生空爆。

碎片云模型將小天體假設(shè)為球體,將其破碎分為幾次進(jìn)行,每次破碎,小天體會(huì)分裂為幾個(gè)主要碎片,主要碎片仍為球體,其結(jié)構(gòu)強(qiáng)度有所增加,可再次發(fā)生破碎,而碎片云則不再繼續(xù)分裂。[8]該模型繼承了煎餅?zāi)Ｐ椭械淖冃畏匠?分裂次數(shù)變?yōu)槎啻?理論上不存在空爆高度的概念,小天體在撞擊大氣層過(guò)程中的變化主要依賴能量沉積進(jìn)行描述。

爆炸流體動(dòng)力學(xué)直接數(shù)值模擬針對(duì)特定的撞擊事件或進(jìn)入大氣過(guò)程中的一個(gè)特定的方面進(jìn)行研究。[6]通過(guò)直接模擬小天體撞擊大氣層過(guò)程中的流動(dòng)、傳熱以及燒蝕相變的能量沉積過(guò)程,發(fā)現(xiàn)小天體首先發(fā)生煎餅狀變形,變形較大 (大于兩倍直徑)時(shí)開(kāi)始發(fā)生瑞利—泰勒不穩(wěn)定性和開(kāi)爾文—霍姆赫茲不穩(wěn)定性,最后破碎為碎片、融化液滴、蒸汽的混合物,初步揭示了空爆的物理機(jī)制。但由于該方法計(jì)算過(guò)程復(fù)雜,時(shí)間成本高,較少用于撞擊參數(shù)的影響研究。

三者對(duì)比,煎餅?zāi)Ｐ陀?jì)算較為簡(jiǎn)便,所需時(shí)間成本較低,同時(shí)可以捕捉小天體撞擊大氣層過(guò)程中能量沉積、空爆解體等顯著特征。該模型模擬過(guò)程中的能量變化,對(duì)于分析空爆的物理機(jī)制有很大幫助。

能量沉積包主要由阻力作用和燒蝕作用引起,本文將首次討論兩種作用對(duì)能量沉積的貢獻(xiàn),并由此對(duì)空爆機(jī)制進(jìn)行分析。

本工作采用小天體撞擊地球大氣層的煎餅?zāi)Ｐ?以車?yán)镅刨e斯克事件為對(duì)照,以能量沉積過(guò)程為基礎(chǔ),分析、探究了小天體的空爆機(jī)制,初步得出了空爆的物理原因。

2 物理問(wèn)題及數(shù)學(xué)模型

將小天體視為均質(zhì)球體,主要輸入?yún)?shù)有：小天體撞擊大氣層的初始速度、撞擊角度、小天體的直徑、小天體的密度。其中,撞擊角度指小天體的撞擊速度與當(dāng)?shù)厮骄€之間的夾角,如圖1所示。

圖1 撞擊角示意圖Fig.1 The impact angle

煎餅?zāi)Ｐ蜑橐活惸Ｐ偷慕y(tǒng)稱,其核心思想是：將小天體視為連續(xù)的變形體,將小天體的破碎過(guò)程視為連續(xù)體受擠壓發(fā)生變形的過(guò)程,當(dāng)變形條件被觸發(fā)時(shí),小天體在氣動(dòng)壓力下發(fā)生變形。在變形過(guò)程中,小天體迎風(fēng)截面的半徑逐漸增大,當(dāng)增大到一定程度時(shí),認(rèn)為小天體發(fā)生空爆。[7]

盡管核心思想相同,但不同的煎餅?zāi)Ｐ椭g仍可能存在較大差異,比如小天體形狀假設(shè)的差異、控制方程的不同、各參數(shù)之間的差別、終止判據(jù)的不同等等。[2,7,11]

Christopher F.Chyba等人于1993年提出了較為經(jīng)典的煎餅?zāi)Ｐ?用于模擬1908年的通古斯卡事件,并列舉了幾組不同參數(shù)的小天體撞擊地球大氣層的模擬結(jié)果。[7]該模型假設(shè)小天體為底面直徑與高相等的均質(zhì)圓柱體,考慮地球曲率的影響,忽略地球自轉(zhuǎn)的影響。在小天體撞擊地球大氣層的過(guò)程中,當(dāng)小天體的駐點(diǎn)壓力大于小天體本身的結(jié)構(gòu)強(qiáng)度時(shí),小天體將受到擠壓而發(fā)生形變,形狀逐漸趨于扁平,迎風(fēng)截面也逐漸增大。

本文所用模型,在小天體飛行的基本控制方程、小天體截面半徑變化方程、發(fā)生空爆的判斷方式等方面,與Chyba等人提出的煎餅?zāi)Ｐ突疽恢?同時(shí)也考慮了地球曲率的影響,忽略地球自轉(zhuǎn)。不同的是,本文假設(shè)小天體為均質(zhì)球體,并在一些參數(shù)上進(jìn)行了細(xì)微的調(diào)整。小天體的速度V、撞擊角γ、軌道半徑r、路程中心角ψ、質(zhì)量m隨時(shí)間t的變化由公式 (1a)至 (1e)給出。

其中,Cd為大氣阻力系數(shù),本程序中取為1;ρ為即時(shí)大氣密度,隨高度發(fā)生變化;A為小天體的即時(shí)迎風(fēng)截面積,由式A=πR2算得,其中R為小天體的即時(shí)迎風(fēng)截面半徑;μ則為地球引力常數(shù),取3.986×105km3/s2;Cl為升力系數(shù),由于小天體自身并無(wú)驅(qū)動(dòng)力,故Cl取0;Ch為傳熱系數(shù),取為0.1;Q為小天體的燒蝕熱,與小天體的材料、結(jié)構(gòu)等因素有關(guān),由于本文假設(shè)小天體為均質(zhì)球體,因此,在本文中,Q主要由小天體的材質(zhì)決定。大氣密度ρ采用冪指數(shù)模型計(jì)算：

其中,ρ0=1.2256 kg/m3為地平線處的大氣密度,h為高度,在程序計(jì)算中為小天體軌道高度r與地球半徑Re之差,單位為m;H為大氣標(biāo)高,為7254 m。

小天體進(jìn)入大氣層后,當(dāng)駐點(diǎn)壓力P大于小天體結(jié)構(gòu)強(qiáng)度S0時(shí),將開(kāi)始發(fā)生變形,對(duì)應(yīng)實(shí)際過(guò)程中的破碎過(guò)程。駐點(diǎn)壓力由式 (3)求得：

在駐點(diǎn)壓力的作用下,小天體受到擠壓發(fā)生變形,呈現(xiàn)出越來(lái)越扁平的趨勢(shì),其迎風(fēng)截面半徑也逐漸增大,截面半徑R的變化率隨時(shí)間t的變化由式 (4)給出：

其中,ρm為小天體的密度。

小天體迎風(fēng)截面直徑D=2R,當(dāng)D變?yōu)槌跏冀孛嬷睆紻0的N倍時(shí),則認(rèn)為小天體發(fā)生了空爆。其中,N為經(jīng)驗(yàn)參數(shù),通常在 5～10之間[2,7],本文在程序合理性驗(yàn)證中取為6,該取值下,程序模擬結(jié)果與該事件真實(shí)情況最為符合。在后續(xù)的主要撞擊參數(shù)掃描中,N的取值保持6,不再變化。

2 模型合理性驗(yàn)證及空爆機(jī)制研究

模型通過(guò)對(duì)車?yán)镅刨e斯克事件進(jìn)行模擬,以驗(yàn)證其合理性。

通過(guò)對(duì)真實(shí)視頻材料等觀測(cè)數(shù)據(jù)的分析[12-14],車?yán)镅刨e斯克事件的主要物理特征已經(jīng)較為明確,為模型驗(yàn)證提供了豐富的對(duì)比材料。

據(jù)觀察,車?yán)镅刨e斯克事件中,小天體于97km高空處進(jìn)入大氣,90km處開(kāi)始形成沖擊波,83km處初步產(chǎn)生塵埃及部分碎片,輻射峰值發(fā)生在29.7km左右[14]。圖2[14]給出了該事件小天體的破碎過(guò)程。

通過(guò)對(duì)錄像的分析與提取,得到了車?yán)镅刨e斯克事件的小天體視幅輻照度光曲線,如圖3所示。[14]

圖2 車?yán)镅刨e斯克事件小天體破碎過(guò)程Fig.2 Meteoroid fragmentation stages in Chelyabinsk event

圖3 車?yán)镅刨e斯克事件光曲線Fig.3 Fireball visual magnitude irradiance light curve in Chelyabinsk event

對(duì)比圖3中視星等和能量沉積兩條曲線,能夠看到,20～50km高度的范圍內(nèi),輻照度與能量沉積曲線的趨勢(shì)相一致,且均在30km高度左右達(dá)到峰值。

車?yán)镅刨e斯克事件主要參數(shù)如表1所示。

表1 車?yán)镅刨e斯克事件參數(shù)表Tab.1 Parameters of Chelyabinsk event

2.1 能量沉積計(jì)算方法

小天體在進(jìn)入地球大氣層時(shí),部分動(dòng)能通過(guò)大氣阻力的作用、熱燒蝕作用轉(zhuǎn)移到大氣層中,該過(guò)程為能量沉積。而能量沉積曲線用于描述不同高度處小天體在單位高度上損失的動(dòng)能,反映小天體在撞擊大氣層過(guò)程中的能量變化及其隨高度的變化情況,是探究空爆過(guò)程及其物理機(jī)制的重要依據(jù)。

能量沉積曲線采取如下計(jì)算方式：

式中,dE為兩步計(jì)算之間小天體的動(dòng)能減少量,dh為對(duì)應(yīng)兩步計(jì)算之間小天體的高度變化量,同樣地,dm、dv也為對(duì)應(yīng)兩步計(jì)算之間質(zhì)量、速度的減少量。各變化量取為原變化量的絕對(duì)值,以保證結(jié)果為正。

式 (5)也給出了能量沉積各項(xiàng)分量,通過(guò)簡(jiǎn)化,保留一階項(xiàng),略去二階、三階項(xiàng),則能量沉積的主要貢獻(xiàn)分為兩項(xiàng),如公式 (6)所示。該式主要用于2.3中拐點(diǎn)分析部分,不用于能量沉積的直接計(jì)算。

2.2 能量沉積曲線對(duì)比

通過(guò)仿真,得到了車?yán)镅刨e斯克事件發(fā)生空爆及之前的能量沉積曲線,如圖4所示。

圖4 車?yán)镅刨e斯克事件能量沉積曲線Fig.4 Energy deposition of Chelyabinsk event

圖4中,黑色實(shí)線代表車?yán)镅刨e斯克事件真實(shí)的能量沉積曲線,基于對(duì)該事件的直接觀測(cè)數(shù)據(jù)分析得出;綠色及藍(lán)色的線條皆為碎片云模型模擬所得的能量沉積曲線[8];而紅色粗體線條則為本程序模擬而得的能量沉積曲線。

從四條曲線的對(duì)比可以看出,本文得到的能量沉積曲線與真實(shí)曲線的整體趨勢(shì)、重要節(jié)點(diǎn)基本一致,尤其在40km后的部分,二者幾乎完全重合。在發(fā)生空爆的30km高度處,能量沉積值也與真實(shí)情況相一致。因此,可以認(rèn)為,該方法能夠初步模擬小天體撞擊大氣層的過(guò)程,其關(guān)于空爆高度及能量沉積的計(jì)算是可信的。

3 空爆機(jī)制分析

探究空爆機(jī)制,并不能僅著眼于發(fā)生空爆的短暫瞬間,而應(yīng)當(dāng)從小天體進(jìn)入大氣至發(fā)生空爆的整個(gè)過(guò)程出發(fā)。能量沉積曲線是記錄該過(guò)程的重要依據(jù),從該曲線出發(fā),分析能量沉積過(guò)程中小天體的變化情況,對(duì)于揭露空爆機(jī)制有重要意義。

從圖4可以看到,三類能量沉積曲線均存在一個(gè)較為明顯的拐點(diǎn),說(shuō)明在該拐點(diǎn)處,小天體發(fā)生了某種變化,改變了能量沉積隨高度變化的速率。

為探究該拐點(diǎn)所對(duì)應(yīng)的物理過(guò)程,即導(dǎo)致該拐點(diǎn)產(chǎn)生的物理原因,根據(jù)公式 (6),分別將能量沉積的主要兩項(xiàng)進(jìn)行展示,第一項(xiàng)為對(duì)應(yīng)速度變化,主要作用因素為大氣阻力的影響,不妨稱之為阻力項(xiàng),其能量沉積貢獻(xiàn)如圖5所示;第二項(xiàng)為,對(duì)應(yīng)質(zhì)量變化,主要作用因素為燒蝕,此處稱之為燒蝕項(xiàng),其能量沉積貢獻(xiàn)如圖6所示。圖7表示能量沉積總量。

圖5 阻力項(xiàng)能量沉積曲線Fig.5 Energy deposition caused by resistance

圖6 燒蝕項(xiàng)能量沉積曲線Fig.6 Energy deposition caused by ablation

圖7 能量沉積總量Fig.7 Energy deposition

從圖中可以看到,第一項(xiàng)與第二項(xiàng)以及總量均在39km處存在拐點(diǎn),并呈現(xiàn)出相似的變化趨勢(shì)。

就整體變化趨勢(shì)而言,第二項(xiàng)與能量沉積總量的變化趨勢(shì)更為接近,兩條曲線幾乎保持平行。

圖8和圖9分別給出了阻力項(xiàng)與燒蝕項(xiàng)對(duì)能量沉積貢獻(xiàn)的占比。

圖8 阻力項(xiàng)能量沉積占比Fig.8 The ratio of energy deposition caused by resistance

圖9 燒蝕項(xiàng)能量沉積占比Fig.9 The ratio of energy deposition caused by ablation

從圖8、圖9可知,燒蝕項(xiàng)所占的比例始終高于77%,遠(yuǎn)大于阻力項(xiàng)能量沉積貢獻(xiàn)的占比,為能量沉積的主要貢獻(xiàn)。

結(jié)合圖6、圖7,可以判斷,在39km左右,小天體燒蝕作用加劇,使得質(zhì)量減少更加迅速,從而使該高度下小天體單位高度的動(dòng)能沉積量變大。

推測(cè)燒蝕加劇的結(jié)果與小天體的變形有關(guān)。

圖10為小天體的迎風(fēng)截面半徑與小天體高度之間的關(guān)系。

圖10 小天體迎風(fēng)截面半徑隨高度的變化Fig.10 Windward section radius variation with altitude

從圖中能夠看出,在小天體的高度到達(dá)39km之前,其截面半徑并未發(fā)生改變。此時(shí)小天體也存在燒蝕作用,但由于駐點(diǎn)壓強(qiáng)小于小天體屈服強(qiáng)度,因此燒蝕造成的質(zhì)量損失主要表現(xiàn)為小天體縱向長(zhǎng)度的減少,并未改變小天體的迎風(fēng)截面。而在39km處,小天體的迎風(fēng)截面半徑開(kāi)始發(fā)生變化,燒蝕作用也隨之加劇,能量沉積進(jìn)一步增加,最終導(dǎo)致了拐點(diǎn)的產(chǎn)生。

小天體開(kāi)始變形的條件為駐點(diǎn)壓強(qiáng)大于小天體本身的結(jié)構(gòu)強(qiáng)度,因此,為了對(duì)拐點(diǎn)產(chǎn)生原因進(jìn)行進(jìn)一步驗(yàn)證,本文改變了小天體的結(jié)構(gòu)強(qiáng)度,并呈現(xiàn)了總能量沉積曲線中拐點(diǎn)的位置變化情況,如圖11、圖12所示。

圖11 S=2S0時(shí)能量沉積曲線Fig.11 Energy deposition with doubled yield strength

圖12 S=0.5S0時(shí)能量沉積曲線Fig.12 Energy deposition with halved yield strength

圖11為S=2S0時(shí)得到的能量沉積曲線,從圖中可以明顯看到,拐點(diǎn)發(fā)生在35km左右;圖5正常情況下,拐點(diǎn)在39km左右;圖12中曲線為時(shí)的能量沉積曲線,其拐點(diǎn)發(fā)生在約42km高度處。其中,S表示小天體結(jié)構(gòu)強(qiáng)度,S0表示表1中基準(zhǔn)結(jié)構(gòu)強(qiáng)度。

小天體屈服強(qiáng)度高時(shí),達(dá)到變形條件對(duì)應(yīng)的高度較低,與拐點(diǎn)產(chǎn)生的高度規(guī)律相符。因此,小天體能量沉積曲線的拐點(diǎn)是小天體受到擠壓開(kāi)始發(fā)生變形所致,而變形引起的主要能量損失則源于燒蝕。

小天體在進(jìn)入大氣的過(guò)程中,首先在氣動(dòng)力、氣動(dòng)熱的作用下,發(fā)生縱向燒蝕,在該階段,由于大氣密度較低,小天體的速度減少幾乎可以忽略;隨著大氣密度的增加,小天體的駐點(diǎn)壓強(qiáng)超過(guò)其自身的結(jié)構(gòu)強(qiáng)度,小天體開(kāi)始發(fā)生橫向的變形,在該過(guò)程中,燒蝕作用加劇,小天體的能量以熱、光等形式向外輻射;小天體的質(zhì)量損失過(guò)大,其燒蝕將減弱。燒蝕作用最為劇烈時(shí),對(duì)應(yīng)小天體的外在表現(xiàn)即為空爆。空爆由劇烈的燒蝕作用引起。

空爆發(fā)生后,小天體在壓力作用下發(fā)生破碎,破碎后的小天體碎片面質(zhì)比增大,單個(gè)碎片在燒蝕作用下剝落、蒸發(fā),形成圖2中I、J所示的燃燒狀煙霧。

4 總結(jié)與討論

小天體撞擊地球首先接觸的便是地球大氣層,空爆為小天體撞擊地球大氣層的重要物理特征,研究空爆對(duì)小天體防御有著重要意義。

本文基于Chyba等人于1993年提出的煎餅?zāi)Ｐ?對(duì)模型進(jìn)行了適當(dāng)?shù)恼{(diào)整,對(duì)小天體撞擊地球大氣層中的空爆現(xiàn)象進(jìn)行了分析,并研究了小天體撞擊地球大氣層的能量沉積過(guò)程,探究了空爆背后的物理機(jī)制。

首先以車?yán)镅刨e斯克事件為基礎(chǔ),對(duì)模型的合理性進(jìn)行了驗(yàn)證,證明該煎餅?zāi)Ｐ统绦蚰軌蛴糜谛√祗w撞擊地球大氣層中空爆高度及能量沉積的計(jì)算,同時(shí)對(duì)能量沉積曲線中的拐點(diǎn),對(duì)能量沉積過(guò)程進(jìn)行了分析,并進(jìn)一步提出了空爆對(duì)應(yīng)的物理過(guò)程,發(fā)現(xiàn)空爆主要由劇烈的燒蝕作用引起。

盡管模型不同,但該結(jié)論與碎片云模型并不矛盾。碎片云模型中,雖沒(méi)有明確的空爆概念,但仍存在幾次小天體破碎的過(guò)程,使小天體的總表面積增大,燒蝕作用也會(huì)因此加劇。

此外,在煎餅?zāi)Ｐ椭?發(fā)生空爆的判斷標(biāo)準(zhǔn)為小天體的迎風(fēng)截面半徑變?yōu)槌跏及霃降腘倍,N為經(jīng)驗(yàn)參數(shù),一般取值在5～10之間,本文基于對(duì)車?yán)镅刨e斯克事件的模擬,將N取為6。

小天體撞擊地球大氣層的過(guò)程是十分復(fù)雜的,包括了聲、光、熱、沖擊波等多種效應(yīng)[2,17,18],本文對(duì)空爆這一重要現(xiàn)象進(jìn)行了探討,分析了小天體發(fā)生空爆的物理原因。若需對(duì)小天體撞擊地球進(jìn)行整體評(píng)估,則應(yīng)當(dāng)在此基礎(chǔ)上逐步完善各種物理效應(yīng),形成一個(gè)整體的評(píng)估方案,以對(duì)行星防御提供更為有利的支持。