小行星進入與撞擊效應評估模型敏感性研究1)

黨雷寧 柳森 白智勇 石義雷

(中國空氣動力研究與發展中心超高速空氣動力研究所,四川綿陽 621000)

引言

近地小行星對地球的撞擊是人類生存和發展面臨的潛在威脅之一.學術界已經基本認同:6500 萬年前一顆直徑10 km 的小行星撞擊現墨西哥希克蘇魯伯地區,導致包括恐龍在內的全球物種大滅絕,史稱K-T 事件[1].1908 年6 月30 日,在俄羅斯西伯利亞的Tunguska 發生小行星撞擊事件,地面大范圍的樹木被破壞,釋放出約千萬噸TNT 當量的能量,受影響面積相當于美國華盛頓市[2].2013 年2 月15 日,在俄羅斯Chelyabinsk,一個尺寸近20 m 的流星在空中爆炸解體,釋放出約50 萬噸TNT 當量的能量,地面建筑物受損,數千名地面人員受到不同程度的傷害[3].

小行星撞擊地球存在復雜的物理化學和力學過程,引起多種類型的進入與撞擊效應及災害[4].小行星在進入地球大氣過程中發生空中爆炸,沖擊波向地面傳播,引起地面超壓損傷;爆炸產生的熱輻射傳至地面,造成熱輻射損傷.小行星穿透大氣后,對陸地的直接撞擊產生撞擊坑,對海洋的直接撞擊可能引發海嘯,直接撞擊還可能把部分物質拋向空中,產生反濺碎片云.對于千米量級的小行星,撞擊產生的反濺碎片云將會引發全球氣候變化.

計算分析小行星進入與撞擊效應,有數值仿真和工程計算兩種手段.數值仿真方法保真度高、耗時長,適用于具體撞擊事件或災害類型的精細分析.Boslough 和Crawford[5]以及Robertson[6]采用沖擊動力學數值方法,對小行星進入過程的解體、空中爆炸以及撞擊海洋海嘯開展了研究.Aftosmis 等[7]以及Collins 等[8]基于可壓縮空氣動力學數值方法,研究了空中爆炸產生爆炸波的傳播過程,評估了工程方法中常用的靜態點爆炸模型的準確性.Johnston 等[9]采用高溫氣體動力學數值模擬手段,耦合高溫燒蝕流場和輻射計算,研究了進入過程產生的熱輻射.在國內,張慶明[10]以及趙君堯[11]利用沖擊動力學數值軟件,初步研究了撞擊坑形成過程.李毅等[12]初步開展了空中爆炸與撞擊海面的數值仿真.

與數值仿真相比,工程方法采用快速預測的模型,能夠在短期內開展大量狀態的災害評估,并獲得可接受準度的計算結果.1993 年,Chyba 等[2]提出并采用Pancake 解體模型,研究了Tunguska 小行星的空中爆炸,獲得與數值仿真、地面試驗一致的空爆高度.同一時期,Hills 和Coda[13]研究了解體、空中爆炸、超壓、熱輻射、撞擊坑、海嘯等各種類型的進入與撞擊效應模型,特別是基于美國核試驗數據給出地面超壓半徑的擬合公式.2005 年,Collins 等[4]廣泛研究了各種類型的進入與撞擊效應工程模型,推導出相應的解析公式,研制了網頁版的小行星撞擊災害評估軟件,放置在互聯網上供開放使用.2017 年,Collins 等[8]又對該軟件進行了改進.2002 和2017 年,美國NASA的科學定義組[14-15](Science Definition Team)利用當時工程模型發展的成果,考慮小行星物理特性和進入參數的不確定性以及撞擊頻率,給出小行星撞擊年度受影響人口及其概率.美國NASA 還于2015—2017年,組織開發了小行星進入與撞擊風險分析評估系統PAIR(Probabilistic Asteroid Impact Risk)[16-17],設定小行星撞擊地球的場景,可定量分析進入和撞擊過程對地球人口和基礎設施的危害,并構成了美國政府、軍方、學術界小行星撞擊聯合應急演習[18-19]的主要內容.在國內,柳森等[20]梳理歸納了小行星撞擊地球的超高速問題,并采用文獻中已有的模型,建立了小行星進入與撞擊效應分析評估軟件AICA,利用該軟件對2019 年行星防御國際會議假想撞擊事件開展分析評估[21].周炳紅等[22]對小行星進入過程能量沉積開展了計算分析.

小行星進入與撞擊效應評估模型是工程計算的基礎與核心.隨著人們對小行星進入與撞擊過程認識的不斷深入,隨著與數值仿真最新成果的結合,這些工程模型不斷迭代和發展,計算結果的準度不斷提高,已經成功地應用于Tunguska,Chelyabinsk,Ko?ice,Bene?ov 等實際撞擊事件的計算分析[2,23-25],也在多次行星防御演習中獲得應用[18-19].盡管如此,當前的工程模型輸入參數不確定度大[15,17],給模型的使用帶來困擾,也較大地影響了危害評估結果.首先,天文觀測給出的潛在防御目標小行星的物理特性和運動參數不確定度大,即材料、尺寸、強度與進入參數通常在較大范圍.其次,由于基礎研究不充分,工程模型中的燒蝕系數、流明效率、擴散系數不確定度大.最后,工程模型特別是某些解體模型存在自由參數.

因此,需要研究模型輸入參數變化對輸出參數的影響,即開展進入與撞擊評估模型敏感性分析.Hills 和Coda[13]研究了不同小行星尺寸、進入速度和材料對超壓損傷范圍的影響,屬于早期(1993 年)的工作,采用的工程模型與目前有所不同.Mathias[26]研究了輸入參數對傷亡人口的影響,重點在于受影響人口,空爆模型也與目前不同.Wheeler 等[23]開展了FCM 解體模型的敏感性研究.這些研究還不夠充分.著重分析輸入參數對能量沉積的影響,缺乏對空爆高度和地面損傷范圍的研究.

本文利用課題組建立的小行星進入與撞擊效應分析評估軟件AICA[20],研究輸入參數的變化對能量沉積、空爆高度和地面損傷范圍等輸出參數的影響,期望得出一些有意義的結論,指導軟件的使用和下一步的研究工作,同時提高對小行星進入與撞擊效應的認識.

1 小行星進入與撞擊效應分析評估軟件AICA 介紹

1.1 AICA 軟件中的工程模型

目前,AICA 軟件包含了進入彈道方程、質量損失方程、解體判據、解體模型、空中爆炸模型,地面超壓損傷范圍模型、地面熱輻射損傷范圍模型等小行星進入與撞擊效應工程模型,暫未包含撞擊坑、海嘯和全球氣候變化模型.下面逐一闡述這些模型.

1.1.1 彈道方程

在AICA 軟件中,考慮到后續研究精確求解彈道的需要,所用彈道方程包含圓球自轉地球模型以及經緯度變量.若僅考慮飛行中的阻力,則單個飛行體的彈道方程[29]為

其中,h為飛行高度;θ 和Φ 分別是經度、緯度;V是飛行速度大小;γ 是飛行航跡角或進入角,定義為速度矢量與當地水平面的夾角,速度矢量在當地水平面上方為正;ψ 是航向角,定義為速度矢量在當地水平面的投影與當地緯線方向的夾角;Cd是阻力系數;A是橫截面積;rt是小行星到地心的距離;ω 是地球自轉角速度;采用美國1976 年標準大氣計算大氣密度ρair.

在本文的進入與撞擊效應工程模型中,有兩種類型的單個飛行體:以剛體運動的小行星和解體后的碎片,以變形體運動的解體后的碎片云(見1.1.4節).為了在后續分析中有所區別,本文分別記這兩種類型單個飛行體阻力系數為Cd,f和Cd,c.

對于小行星和碎片,阻力系數與外形和飛行姿態相關,橫截面積與外形相關.小行星和碎片是非規則外形,由于天文觀測方面的基礎研究不充分,大多數進入過程研究都假定圓球外形[4,8,13-18,23-28].Wheeler 等[24]在圓球外形的假設下利用改進的FCM 解體模型,獲得了與多個小行星進入過程能量沉積觀測數據一致的計算結果.Ceplecha 和Revelle[30]、Borovi?ka 和Kalenda[31]、Borovi?ka 等[32]在運動方程中利用參數ΓA′考慮外形和姿態對阻力系數的影響,成功地開展了多個小行星進入過程觀測數據分析及彈道重構.其中,Γ 是Ceplecha 定義的阻力系數,等于0.5Cd,f;A′是形狀因子,圓球的形狀因子為1.21.小行星的氣動升力與外形相關,受姿態影響較大.一方面,小行星在進入前往往存在自旋[33],進入大氣后自旋有可能繼續存在[34];另一方面非規則外形的姿態運動受氣動力矩、燒蝕引射和解體爆炸影響,變化規律非常復雜,迄今為止還沒有相關研究.當前的研究不考慮升力影響,從研究結果上看與觀測數據基本一致[2,24,30,32].

對于以變形體運動的碎片云,根據Pancake 模型[2,13],其外形是壓扁的“薄煎餅”形狀,橫截面積A=πr2(r為半徑).由于碎片云不是剛體,在運動中不斷變形,僅需要考慮阻力.

Wheeler 等[16-17]在PAIR 軟件中以及Collins等[4,8]在其危害評估軟件中對這兩種類型的單個飛行體,都采用1.0 的阻力系數.白智勇等[35]采用高超聲速化學非平衡流動數值模擬手段,獲得了直徑1 m的圓球小行星在速度10～20 km/s 的阻力系數,范圍為0.87～0.92.Ceplecha 和Borovi?ka[30-32]在多個進入過程觀測數據分析及彈道重構中,采用的小行星及碎片ΓA′范圍為0.7～1.14,若以圓球形狀因子換算,則Cd,f(等于2Γ)的范圍為1.15～1.88.

基于以上分析,本文假定小行星及碎片是圓球外形,阻力系數范圍為0.8～2.0.通過考察阻力系數對進入與撞擊效應的影響,間接研究外形和姿態的影響.對于以變形體飛行的碎片云,由于缺乏研究支撐,無法獲得阻力系數取值范圍,阻力系數取Wheeler 和Collins 研究中的1.0.

1.1.2 質量損失方程

小行星由于燒蝕而不斷損失質量,可以直接使用文獻中廣泛應用的質量損失方程[2,4,23-28]

其中,m是小行星瞬時質量,σ 是燒蝕系數,CH是熱流系數,Q是燒蝕熱.CH是小行星尺寸、飛行速度、高度的函數,Q與小行星的材料、結構相關.基于當前學術界的認識[2,17,36-37],燒蝕系數σ 取值范圍為:3.5×10?10～7.0×10?8kg/J.

1.1.3 解體判據

在天文學中,廣泛使用如下解體判據[2,4,17,36]

其中,P是小行星駐點壓力,S是解體強度.即當駐點壓力超過強度,小行星解體.

小行星與隕石的強度不同.隕石經歷了嚴酷的氣動力、熱而降落至地面,代表了小行星結構中強度最大的部分.而小行星通常是包含裂紋、空穴的碎石堆結構,強度要小得多.根據Popova 等[38](2011)對大量觀測結果的總結,0.1～106kg 質量范圍的石質小行星的初次解體強度為0.1～10 MPa,而材料的平均拉伸強度為30 MPa,壓縮和屈服強度更大.

Weibull 和Sweden[39]統計分析了地球上的巖石強度,認為強度和質量服從如下統計規律

其中,Sc是巖石質量mc的強度,Sp是同一巖石材料質量mp的強度;α 是強度指數,0 ＜α ＜1,不同材料有不同的取值.Cotto-Figueroa 等[40]通過對不同尺寸隕石開展強度試驗,研究了兩種隕石材料強度的統計散布,并與小行星進入觀測得到的解體強度對比,初步認為小行星材料也服從Weibull 定律.

一些解體模型應用Weibull 定律計算解體后的碎片強度[8,23],此時下標“c”和“p”分別表示解體后的子碎片以及解體前的父碎片.Wheeler 等[23]在FCM解體模型的敏感性分析中,認為α 在0.1～0.3 范圍內取值可以得到合理的、能夠代表小行星結構變化對解體影響的能量沉積結果,在其利用PAIR 軟件開展撞擊風險分析中也取這個范圍[17].NASA 科學定義組在2017 年的報告中[15],對含水和不含水的石質小行星分別采用0.1 及0.2 的取值,并認為鐵質小行星不解體.

基于此,本文在敏感性分析中,認為強度指數α的取值范圍為0.1～0.3.

1.1.4 解體模型

解體后,需要應用解體模型確定小行星的存在形式.目前存在多種解體模型,包括連續的Pancake模型[2,13]和Park 模型[41]、離散的“collective wake”模型[42]和“independent wake”模型[43],以及把離散和連續模型綜合起來的耦合模型[44]及其改進后的FCM模型[23-24].

FCM 模型認為,小行星每次解體為一系列大尺度的碎片和一個由小碎片、液滴組成的碎片云.所有的碎片都是獨立的飛行體,其運動服從方程組(1)和(2).碎片云的運動按照Pancake 模型計算,其整體運動在服從方程組(1)和(2)的同時,橫截面積還在前后壓力差的作用下逐漸增大.Hills 和Coda[13]根據能量守恒關系,推導出碎片云變形中橫向擴散速度公式

其中,Vdisp是碎片云的擴散速度,V是飛行速度,ρair和ρ 分別是大氣及小行星密度.Cdisp是擴散系數,Hills 和Coda[13]認為其量值為3.5,Passey 和Melosh略有不同,認為其量值為3.0.Wheeler 等[23]在FCM模型的敏感性分析中,認為Cdisp的變化對能量沉積影響很小.因此在本文中,不研究Cdisp的變化對模型輸出的影響.

1.1.5 空中爆炸模型

在進入大氣過程中,小行星可發生多次空中爆炸.目前的空中爆炸模型尚不能處理多次空爆情形,而認為空中爆炸是發生在某高度的點爆炸[2,4,7,13-17],稱為靜態點爆炸模型.點爆炸的高度稱為空爆高度,點爆炸的能量稱為空爆能量.

Chyba 等[2]以及Collins 等[4]基于Pancake 模型,認為當碎片云半徑擴散到小行星初始半徑的某倍數時,發生空中爆炸.NASA Ames 在建立PAIR 軟件之初,認為空爆高度是瞬時動能減小到初始進入動能一半時的高度[26].NASA 科學定義組2003 和2017 年的報告[14-15]以及在PAIR 軟件目前的版本[17]中,都認為空中爆炸發生在能量沉積峰值處.

Hills 和Coda[13]在1993 年認為空爆能量是全部動能,1998 年又認為空爆能量是全部動能的一部分[36],Toon 等[45]認為這個比例是50%.Boslough和Crawford[5]發現沖擊動力學數值模擬預測的地面超壓大于點爆炸模型預測結果,而Collins 等[8]卻認為在距離較遠處兩種模型預測的地面超壓基本一致.Wheeler 等[17]基于上述Boslough 的結果,認為把空爆能量定義為全部動能可以更準確地預測地面超壓,并將其應用在PAIR 軟件中.

基于目前的研究現狀,可以認為空中爆炸是發生在能量沉積峰值處以全部初始動能為能量的點爆炸.

1.1.6 地面超壓損傷模型

1993 年,Hills 和Coda[13]基于空中爆炸的靜態點爆炸模型以及美國核試驗數據,擬合出4 psi 地面超壓損傷半徑的經驗公式,Wheeler 等[17]將其用于PAIR 軟件中.該公式為

其中,rp為地面超壓超過4 psi 的半徑,單位是km;h是空爆高度,單位是km;E是空爆能量,單位為MT.MT 是能量單位,指的是106噸TNT 當量的能量,1 MT=4.184×1015J.4 psi 的超壓可以使木框架房屋倒塌[4].

2005 年,Collins 等[4]對美國核試驗數據進行擬合,推導出另一套地面超壓的經驗公式,并于2017 年進行了改進[8].Collins 的經驗公式把自空爆點正下方(文獻中稱為ground zero) 沖擊波的影響范圍分為正則反射與馬赫反射兩個區域,并分別進行計算.在正則反射區域中

其中,p是地面超壓,單位是Pa;r1和h1分別是1 KT爆炸當量時距空爆點正下方的水平距離及空爆高度,可通過縮尺率計算

式(9)中的1000 是把km 轉換為m、MT 轉換為kT所需的因子.在馬赫反射區域中

正則反射區域與馬赫反射區域的分界點為

Aftosmis 等[7]通過數值模擬,考察了在推導式(7)～式(11) 中使用的爆炸當量與距離的縮尺率關系,認為浮力破壞了縮尺率關系,在爆炸能量小于數十MT 時,可以使用縮尺率,而當爆炸能量大于200 MT 時,必須考慮浮力的影響.NASA 科學定義組在2017 年的報告中[15],把美國核試驗數據與Aftosmis 的數值模擬結果綜合起來,評估空中爆炸引起的地面超壓.

在AICA 軟件中,暫未包含數值模擬結果,主要使用Hills 和Coda 的經驗公式(式(7))計算地面超壓損傷范圍;而在爆炸能量較低、超壓小于1 psi 時,使用Collins 的經驗公式(式(8)～式(11)) 計算超壓損傷半徑.圖1 針對1KT 空爆能量下4 psi 超壓損傷半徑與空爆高度的關系,對比了Hills 和Coda 經驗公式、Collins 經驗公式的計算結果與美國核試驗數據(圖中標注為Glasstone&Dolan).可以看到,兩種經驗公式都能獲得與美國核試驗數據基本一致的結果.

圖1 Hills 和Coda 的經驗公式與Collins 經驗公式計算結果的對比(1 kT 空爆能量,4 Psi 超壓)Fig.1 Comparison of Hill&Coda’s and Collins’empirical formula for 1 kT airburst energy and 4 psi overpressure

1.1.7 地面熱輻射損傷模型

Collins 等[4]基于美國核試驗數據和理論分析,推導出在地面空爆時火球熱輻射損傷范圍的經驗公式,Wheeler 等[17]將其修改到適用于在某高度的空爆情形

其中,rr是地面熱輻射損傷半徑,單位km;h是空爆高度,單位km;E是空爆能量,單位為MT.η 是流明效率,是輻射能占初始動能的比例.Collins[4]認為η不確定度較大,取值范圍為:1.0×10?4～1.0×10?2.Φi是熱暴露,即單位面積的傳熱.Zi是1MT TNT 當量的熱暴露,對于3 級燒傷,Zi=0.42 MJ/(m2MT1/6).3級燒傷所需能量大于把草點燃的能量[4].AICA 軟件也包含了1 級、2 級燒傷的輸出.本文中若無特別說明,熱輻射損傷半徑指3 級燒傷的損傷半徑.

1.2 AICA 軟件的計算流程

AICA 軟件的計算流程如圖2 所示.該軟件用彈道計算把各工程模型串聯起來.通過數值求解彈道方程(1)和質量損失方程(2)獲得單體運動參數和能量沉積.在彈道計算中通過式(4) 判斷是否解體,若解體則通過FCM 模型得到解體后碎片及碎片云的數量和質量分數,通過式(5)計算解體后子碎片的強度,通過式(6) 計算碎片云的擴散.當全部碎片和碎片云都完全燒蝕或落地,則計算整個進入過程的能量沉積,并根據空中爆炸模型獲得空爆高度和空爆能量.最后基于空爆高度和空爆能量,按照式(7)～式(13)計算地面超壓損傷半徑和熱輻射損傷半徑.

1.3 AICA 軟件驗證

作者2018 年的論文[20],通過Chelyabinsk 小行星進入和Tunguska 小行星進入兩個算例,初步驗證了AICA 軟件.本文在此基礎上補充危害評估結果,進一步對AICA 軟件進行驗證.

圖3 是Chelyabinsk 事件的計算結果,其中進入彈道(圖3(a)) 和能量沉積結果(圖3(b)) 來自作者2018 年論文.計算條件見作者2018 年的論文[20],其中小行星物理特性及進入條件來自觀測結果,小行星強度、強度指數及碎片云質量分數來自Wheeler 等[17]的研究結果.可以看到進入軌跡和能量沉積曲線都與觀測符合較好;能量沉積曲線能夠復現 Chelyabinsk 小行星兩次空爆的特征.在Chelyabinsk 小行星進入后,當地收到建筑物玻璃破壞的報告,圖3(c) 列出玻璃無損、受損、嚴重受損以及收集到隕石的經緯度位置[3],根據Popova 等[3]的研究:500 Pa 超壓的損傷半徑可用來衡量玻璃破壞的最大范圍,1000 Pa 超壓的損傷半徑可用來衡量玻璃嚴重破壞的范圍.本文500 Pa 和1000 Pa 的超壓損傷范圍計算結果與觀測結果基本一致.

圖2 AICA 軟件流程圖Fig.2 Flow chart of AICA code

圖3 Chelyabinsk 流星事件計算結果Fig.3 Computational results of Chelyabinsk event

圖4 給出Tunguska 大爆炸事件的計算結果,其中能量沉積結果(圖4(a))來自作者2018 年論文.計算條件見作者2018 年的論文[20],其中小行星種類和初始進入參數采用Chyba 等的研究成果,解體強度、燒蝕系數和流明效率在石質小行星的范圍內取值,強度指數、每次解體后生成的子碎片(云)數及其質量分數通過假設給定.由圖4(a)可見,空爆高度為9.2 km,與文獻中的8.5 km[2]比較接近.4 psi 的超壓可以使木框架房屋倒塌,Wheeler 等[25]也把4 psi 作為研究Tunguska 森林倒伏的超壓閾值.由圖4(b)可見,4 psi 超壓損傷范圍能夠反映森林倒伏范圍,比森林倒伏范圍略小.當森林樹木受到熱輻射著火后,火勢會蔓延,導致觀測到的森林著火范圍(圖中標注為forest fire,observed)比最初的熱輻射燒傷范圍(圖中標注為radiation burn,observed)大一些[9].2 級燒傷相當于點燃闊葉林[4],其損傷半徑的計算結果與觀測的熱輻射燒傷范圍相當.本文的計算與觀測結果還存在一定偏差,原因如下:一是Tunguska 當地森林樹木倒伏與超壓的定量關系、當地森林樹木著火與熱輻射的定量關系還在繼續研究中;二是地形變化可能對地面損傷有影響.

圖4 Tunguska 大爆炸事件計算結果Fig.4 Computational results of Tunguska event

本節以Chelyabinsk 與Tunguska 兩次小行星進入事件為算例,說明在給定小行星物理特性、進入參數以及仔細選擇解體模型參數的情況下,AICA 軟件的計算結果是可靠的.

2 工程模型的敏感性分析

2.1 輸入參數范圍

輸入參數的范圍和基準值如表1 所示.本文在基準值的基礎上,通過改變一個或多個輸入參數,研究進入與撞擊效應工程模型的輸出對輸入參數的敏感性.其中重點研究的輸出參數為空爆高度hb,4 psi超壓損傷半徑rp和3 級燒傷半徑rr,也有對彈道和能量沉積的研究.

如表1 所示,小行星直徑d0變化范圍下限為60 m,旨在包含Tunguska 小行星直徑(目前學術界認為是50～80 m[25]);而20 m 直徑的Chelyabinsk 小行星不足以產生嚴重的地面超壓和熱輻射損傷,在表1 中沒有包含.直徑上限取300 m,不僅大于潛在威脅小行星的最小直徑140 m,而且使得地面災害以超壓和熱輻射為主,不用考慮海嘯和全球氣候變化等災害[15].表1 中進入速度V0的范圍取自Greenstreet等[46](2012 年)的數值模擬結果,忽略概率低的高于40 km/s 的值,基準值20 km/s 是行星防御研究的典型速度.進入角γ0的基準值?45?是進入角范圍中概率最大的取值[17].燒蝕系數σ、初次解體強度S0、強度因子α、流明效率η 的取值范圍在1.1 節中已經敘述,基準值取自范圍中間值或對數中間值.小行星及碎片的阻力系數Cd,f的取值范圍在1.1 節中已經敘述,基準值為Wheeler 和Collins 研究中的取值1.0.

近地小行星絕大多數都是石質小行星,密度不僅與具體組成(含水、不含水)、隕落概率以及孔隙率有關,也和強度相關[15,17].本文不研究密度對進入與撞擊評估模型輸出的影響,而根據Mathias 等[17]的結果把基準值取為石質小行星中密度概率最大的值.在NASA 科學定義組2017 年的報告中,直接讓碎片云質量分數等于50%來評估地面災害.Wheeler 等在用PAIR 軟件計算Chelyabinsk 流星事件的能量沉積中采用了大于80%的碎片云質量分數,在研究Tunguska 小行星可能物理特性時采用了80%的取值.在本文中,取80%作為碎片云質量分數的基準值,2.7 節中將會再次討論這一取值的合理性.

表1 進入與撞擊效應模型敏感性研究輸入參數基準值及變化范圍Table 1 Baseline value and variation range of AICA input parameters for sensitivity study

2.2 空爆高度、能量對地面損傷半徑的影響

在本文的模型中,地面損傷范圍與空爆高度、能量直接相關.為此,本文先討論空爆高度、能量與地面損傷半徑的關系,為下面的分析打下基礎.

利用Hills 和Coda 的經驗公式(式(7)),對表1 得到的空爆能量范圍,計算4 psi 超壓損傷半徑隨空爆高度的變化曲線,并如圖5 所示,圖中各條曲線表示不同的進入動能.正如Hills 和Coda[13]的研究并分析圖5,超壓損傷半徑隨空爆高度呈非線性變化;對于某空爆能量,存在使超壓損傷范圍最大的“最優高度”,小于此高度,空爆高度隨超壓損傷半徑遞增,大于此高度,則遞減;若空爆高度不變,則超壓損傷半徑隨空爆能量遞增.

圖5 不同空爆能量下4 psi 超壓損傷半徑隨空爆高度變化曲線Fig.5 Radius of 4 psi overpressure as a function of airburst altitude for burst energy levels computed from Table 1

熱輻射損傷半徑隨空爆高度、能量的變化則簡單得多.分析式(12)～式(13),空爆高度的減小、空爆能量的增大,都使得熱輻射損傷半徑增大.

2.3 直徑、進入速度和進入角的影響

圖6 給出在小行星直徑和進入速度都變化的情況下,空爆高度、4 psi 超壓損傷范圍、3 級燒傷范圍的分布云圖,圖中的黑色等值線表示空爆能量,單位為MT.在圖6 的計算中,除直徑和進入速度外的其它輸入參數取表1 中的基準值.根據圖6(a),大尺寸、小進入速度的小行星具有較低的空爆高度.能量沉積的主要影響機理是解體爆炸.在其他參數都相同的情況下,小行星尺寸越大,解體后的碎片尺寸就越大,就越有能力在較低的高度存在飛行且繼續解體的碎片,并對低高度的能量沉積產生貢獻.進入速度越小,就需要在較低的高度獲得足夠的動壓以解體,并在低高度沉積能量.

如前所述,地面損傷范圍是空爆高度和能量的函數.在圖6(a) 中,空爆高度總體上變化不大(0～20 km).而進入動能隨直徑的立方和速度的平方變化,對超壓損傷半徑的影響更大.因此,超壓損傷半徑總體上隨小行星直徑和進入速度,也即隨著空爆能量的增大而增大(圖6(b)).同理,3 級燒傷半徑主要受空爆能量影響,總體上隨著直徑和進入速度的增大而增大(圖6(c)).

圖6 小行星直徑和進入速度對空爆高度、4 psi 超壓損傷半徑以及3級燒傷半徑的影響(圖中黑色直線等值線為空爆能量,單位MT)Fig.6 Effects of pre-entry diameter and entry velocity on airburst altitude,radius of 4psi overpressure,and radius of 3rd degree burn.The black contour lines indicate burst energy levels with a unit of magaton(MT)

圖7 小行星進入角和進入速度對空爆高度、4 psi 超壓損傷半徑以及3 級燒傷半徑的影響Fig.7 Effects of entry angle and velocity on airburst altitude,radius of 4 psi overpressure and 3rd degree burn

圖7 給出不同進入速度下進入角對空爆高度和地面損傷半徑的影響,圖中其他輸入變量取表1 中的基準值.在小進入角下(大和小指絕對值,后文相同),小行星的水平速度可能大于第一宇宙速度,將會飛離地球而不會進入.圖7(a) 中的灰色區域給出速度12～40 km/s、其他參數取基準值時飛離地球的進入角范圍,約為?10?～0?.隨著進入角的增大,小行星穿透大氣層的時間更短,更有可能在較低的高度沉積動能,導致空爆高度降低.進入角對地面損傷半徑的影響可從空爆高度的變化解釋.在較大的進入角下,空爆高度較低且變化小(圖7(a)),對損傷半徑的影響小(圖7(b)).當進入角較小時,空爆高度可能增大到“最優高度”以上,在進入角向飛離地球的臨界值減小的過程中,空爆高度增大,使得地面超壓損傷半徑減小,也使得熱輻射損傷半徑減小.從圖7(b)還可看出,相同計算條件下4 psi 地面超壓半徑要大于3 級燒傷半徑,這與Wheeler 等[17]的結果相同.

2.4 燒蝕系數、初次解體強度和強度指數的影響

通過同時改變燒蝕系數、初次解體強度和強度指數,其他輸入參數取表1 中基準值,計算得到如圖8 所示的對空爆高度、4 psi 超壓損傷半徑、3 級燒傷損傷半徑的影響曲線,圖8(a) 和圖8(b) 中每條曲線代表不同的燒蝕系數,圖8(c) 中每條曲線代表燒蝕系數與地面損傷類型的組合.

圖8 燒蝕系數、初次解體強度和強度指數對空爆高度、4 psi 超壓損傷半徑以及3 級燒傷半徑的影響Fig.8 Effects of ablation coefficient,strength at 1st breakup and strength scaling exponent on airburst altitude,radius of 4 psi overpressure,and radius of 3rd degree burn

由圖8(a) 和圖8(b) 可見,燒蝕系數和初次解體強度對空爆高度影響較大,強度指數對空爆高度的影響較小.初次解體強度越大,小行星就越趨向于在較低的高度獲得足夠的動壓、初次解體并沉積能量,使得空爆高度減小.燒蝕系數越大,在下降相同高度的過程中碎片和碎片云就能沉積更多動能,使得在以后的飛行中需要沉積的動能減小,導致能量沉積峰值高度即空爆高度上移.

在圖8 的計算條件下,空爆高度總體上小于“最優高度”,與地面超壓損傷半徑呈遞增關系(圖5),與熱輻射損傷半徑呈遞減關系,使得隨著初次解體強度的增大以及燒蝕系數的減小,4 psi 超壓損傷半徑減小,3 級燒傷半徑增大(圖8(c)).總體上看,在圖8的計算條件下,地面損傷范圍變化較小.

2.5 碎片云質量分數和解體后子碎片數的影響

碎片云質量分數是影響能量沉積、空爆高度和地面損傷范圍的關鍵參數之一.

圖9(a)給出不同碎片云質量分數下進入過程的能量沉積情況.可以看到,能量沉積曲線不僅在整體上呈現出大的“鼓包”,而且在局部有小峰值.曲線的這種形態主要由小行星解體行為決定.若僅有一次解體,如Pancake 模型計算結果[2,13](相當于Ccloud=100%),或者當碎片云質量分數較大時(如70%,80%),峰值由初次解體決定,只存在一個大的“鼓包”.若碎片云質量分數較小(如30%～50%),則大質量碎片的再次解體在能量沉積中占主導作用,會呈現數個局部峰值.

由于本文把空爆高度取為能量沉積峰值所在的高度,碎片云質量分數較低時的局部峰值現象,將使空爆高度結果呈現較大波動,并影響地面損傷范圍.圖9(b)、圖9(c) 分別展示了空爆高度、地面損傷范圍隨碎片云質量分數的變化,圖中也給出解體后子碎片數n取2,4,6 的對比.可以看出,空爆高度總體上隨著碎片云質量分數的增大而減小.當碎片云質量分數小于50%,空爆高度和地面損傷范圍曲線存在明顯波動,且空爆高度量值較大;而當碎片云質量分數大于50%時,空爆高度和地面損傷范圍曲線則平滑變化,空爆高度量值小得多.

圖9 碎片云質量分數和解體后子碎片數對能量沉積、空爆高度和地面損傷范圍的影響(4 Psi 超壓半徑,3 級燒傷半徑)Fig.9 Effects of mass fraction of debris cloud and number of son fragments at each breakup on airburst altitude,radius of 4-psi overpressure,and radius of 3rd degree burn

當解體后的子碎片數n從2 增大到6,圖9(b)中的各條空爆高度曲線逐漸靠攏,n取4 和6 時的曲線已經非常接近,說明本文采用n=4 的基準值可以得到n取更大值時計算結果.還可看到,當碎片云質量分數大于70%時,解體后的子碎片數對空爆高度和地面損傷范圍影響較小.這是因為,在碎片云質量分數較大時,碎片云的減速和燒蝕主導了能量沉積過程.

2.6 流明效率和阻力系數的影響

在本文的進入與撞擊效應評估模型中,流明效率η 僅用于計算地面熱輻射損傷半徑(式(12)),因此僅影響熱輻射損傷范圍,而不影響超壓損傷范圍.

在前文的分析中,發現相同的計算條件下4 psi超壓損傷半徑總大于3 級燒傷半徑.那么在流明效率具有表1 所示的不確定度情況下,這個結論是否仍然成立？為此,計算了流明效率取表1 中最小、最大值時的3 級燒傷半徑,并與相同條件下的4 psi 超壓損傷半徑進行比較,如圖10 所示.圖10 的橫軸為進入速度,旨在通過速度變化說明不同空爆能量的影響.可見,3 級燒傷半徑在大多數情況下要小于4 psi超壓損傷半徑,而在進入速度或進入能量較大時,由流明效率最大值得到的3 級燒傷半徑略大.Aftosmis等[7]指出,當進入能量大于數百MT 時,公式(7)預測的地面超壓范圍小于精確的數值模擬結果.因此,可以認為在考慮流明效率不確定度的條件下,4 psi超壓損傷半徑大于3 級燒傷半徑.Tunguska 流星事件中樹木倒伏范圍大于著火范圍(見圖4(b) 或文獻[47]),也是旁證.

圖10 流明效率對3 級燒傷半徑的影響Fig.10 Effect of luminous efficiency on radius of 3rd degree burn

圖11 給出不同小行星及碎片阻力系數下的能量沉積曲線,圖中其他輸入參數來自表1 中的基準值.在高度100 km 到初次解體位置43 km 之間,小行星的氣動阻力是產生能量沉積的重要原因,阻力系數越大,能量沉積越大,但影響總體較小.小行星解體后,在當前較大的碎片云質量分數基準值下,對能量沉積起主要作用是碎片云的減速和燒蝕,碎片的阻力系數對能量沉積影響較小,也對空爆高度和地面損傷范圍影響較小.當然,在小的碎片云質量分數情況下,碎片的阻力系數可能對能量沉積有較大影響,但是正如2.7 節的討論,在本文研究的小行星直徑范圍內,80%的碎片云質量分數基準值有較大的合理性.此外,在圖11 的計算條件下,碎片云的減速對能量沉積貢獻較大.由于缺乏碎片云阻力系數基礎研究支撐,無法獲得合理的阻力系數變化范圍,本文沒有研究工程模型輸出對碎片云阻力系數的敏感性.

圖11 小行星和碎片的阻力系數對能量沉積的影響Fig.11 Effect of drag coefficient of asteroid and fragment on energy deposition

2.7 討論

本章通過同時改變1～3 個表1 中的輸入參數,研究進入與撞擊效應評估模型對輸入參數敏感性.從研究方法上看,沒有覆蓋更多變量的同時變化.然而,正如前文對表1 的說明,輸入參數的基準值代表了取值范圍中概率最大或文獻研究的推薦值,因此分析結果仍然具有一定的意義.

由2.3～2.6 節的分析可見,除碎片云質量分數之外的其他輸入參數,對空爆高度和地面損傷范圍的影響比較“平滑”,沒有出現“跳躍”的變化.可以簡單推論,在這些輸入參數基準值下僅同時改變少數變量得到的定性分析結論,與更多輸入變量同時變化的情形相比,不會有質的變化.

在碎片云質量分數小于50%時,空爆高度和地面損傷范圍隨碎片云質量分數的變化趨勢出現了波動.本文沒有研究碎片云質量分數小于50%時,小行星直徑、進入速度、進入角、初次解體強度和強度指數的影響.Wheeler 等[24]在利用改進的FCM 模型研究Ko?ice,Bene?ov,Tagish Lake 等小能量進入事件能量沉積的過程中,碎片云質量分數總體大于50%;在進入能量更大的Chelyabinsk 事件能量沉積計算中,碎片云質量分數大于80%[23].而對于近代以來進入能量最大的Tunguska 事件,Chyba 等[2]在計算時采用了Pancake 模型,意味著碎片云質量分數為100%.Wheeler 等和Chyba 等的能量沉積計算結果都與觀測或數值模擬符合較好.因此,在本文研究的小行星尺寸條件下(下限是Tunguska 小行星),有理由認為采用80%碎片云質量分數基準值獲得的結果,能夠代表大多數進入情況.

3 結束語

小行星進入與撞擊效應評估模型輸入參數不確定度大,給模型的使用帶來困擾,也較大地影響了危害評估結果.本文首先闡述了模型的由來和發展現狀.然后分析了輸入參數的取值范圍,并以范圍中概率最大或文獻的推薦值作為敏感性研究的基準值.最后利用小行星進入與撞擊效應分析評估軟件AICA,在輸入參數基準值的基礎上,通過改變一個或多個輸入參數,研究能量沉積、空爆高度、4 psi 超壓損傷半徑和3 級燒傷半徑等模型輸出參數對輸入參數的敏感性.研究的模型輸入參數包括:小行星直徑,進入速度,進入角,燒蝕系數,初次解體強度,強度指數,碎片云質量分數,解體后的子碎片數,流明效率,小行星和碎片阻力系數.

對于直徑60～300 m、進入速度12～40 km/s 的石質小行星,本文的研究指出:

(1)大尺寸、小進入速度的小行星具有較低的空爆高度.超壓損傷半徑和熱輻射損傷半徑總體上隨著小行星尺寸和進入速度的增大而增大.

(2)在小進入角下,小行星的水平速度可能大于第一宇宙速度,將會飛離地球而不會進入.

(3) 隨著初次解體強度的增大或燒蝕系數的減小,空爆高度減小.小行星及碎片阻力系數、強度指數對空爆高度和地面損傷范圍影響較小.

(4)碎片云質量分數是影響能量沉積、空爆高度和地面損傷半徑的關鍵參數之一.若碎片云質量分數小于50%,能量沉積曲線出現數個局部峰值,使得空爆高度和地面損傷半徑評估結果出現波動.空爆高度總體上隨著碎片云質量分數的增大而減小.通過分析Wheeler,Chyba 等對若干小行星進入事件能量沉積的計算結果,可以認為在本文研究的小行星直徑范圍內,較大的碎片云質量分數基準值(如80%)能夠代表大多數進入情況.

(5)流明效率僅影響熱輻射損傷范圍,不影響超壓損傷范圍.在本文的計算條件下,3 級燒傷半徑小于4 psi 超壓損傷半徑,小行星進入與撞擊危害以空中爆炸產生的地面超壓為主.

小行星撞擊地球危害評估是行星防御研究的重要方向之一,當前國內的研究處于起步階段,國外的研究中還有大量問題沒有解決.通過本文的敏感性分析,可以認為在后續的工作中,一方面需要在天文和進入過程觀測中,建立小行星尺寸、材料、進入速度、進入角、強度等物理和力學特性的高保真度模型;另一方面,需要在進入過程研究中,著重研究燒蝕、空中爆炸、沖擊波效應和熱輻射問題,持續改進工程模型并減小相關參數的不確定度.