小天體撞擊地球大氣層的空爆問題研究

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

煎餅?zāi)Ｐ蜑橐活惸Ｐ偷慕y(tǒng)稱,其核心思想是：將小天體視為連續(xù)的變形體,將小天體的破碎過程視為連續(xù)體受擠壓發(fā)生變形的過程,當(dāng)變形條件被觸發(fā)時,小天體在氣動壓力下發(fā)生變形。在變形過程中,小天體迎風(fēng)截面的半徑逐漸增大,當(dāng)增大到一定程度時,認為小天體發(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)的影響。在小天體撞擊地球大氣層的過程中,當(dāng)小天體的駐點壓力大于小天體本身的結(jié)構(gòu)強度時,小天體將受到擠壓而發(fā)生形變,形狀逐漸趨于扁平,迎風(fēng)截面也逐漸增大。

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

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

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

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

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

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

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

2 模型合理性驗證及空爆機制研究

模型通過對車里雅賓斯克事件進行模擬,以驗證其合理性。

通過對真實視頻材料等觀測數(shù)據(jù)的分析[12-14],車里雅賓斯克事件的主要物理特征已經(jīng)較為明確,為模型驗證提供了豐富的對比材料。

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

通過對錄像的分析與提取,得到了車里雅賓斯克事件的小天體視幅輻照度光曲線,如圖3所示。[14]

圖2 車里雅賓斯克事件小天體破碎過程Fig.2 Meteoroid fragmentation stages in Chelyabinsk event

圖3 車里雅賓斯克事件光曲線Fig.3 Fireball visual magnitude irradiance light curve in Chelyabinsk event

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

車里雅賓斯克事件主要參數(shù)如表1所示。

表1 車里雅賓斯克事件參數(shù)表Tab.1 Parameters of Chelyabinsk event

2.1 能量沉積計算方法

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

能量沉積曲線采取如下計算方式：

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

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

2.2 能量沉積曲線對比

通過仿真,得到了車里雅賓斯克事件發(fā)生空爆及之前的能量沉積曲線,如圖4所示。

圖4 車里雅賓斯克事件能量沉積曲線Fig.4 Energy deposition of Chelyabinsk event

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

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

3 空爆機制分析

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

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

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

圖5 阻力項能量沉積曲線Fig.5 Energy deposition caused by resistance

圖6 燒蝕項能量沉積曲線Fig.6 Energy deposition caused by ablation

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

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

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

圖8和圖9分別給出了阻力項與燒蝕項對能量沉積貢獻的占比。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 總結(jié)與討論

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

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

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

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

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

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