小天体撞击地球大气层的空爆问题研究

耿淑娟,周炳红,韩鹏,郑伟,李明涛

(1.中国科学院国家空间科学中心,北京 100190;2.中国科学院大学天文与空间科学学院,北京 100049)

1 引言

小天体在撞击地球过程中,最先与地球大气层接触。在进入大气层的过程中,小天体在气动力、气动热的作用下,能量不断减少,产生烧蚀、破碎现象,作用剧烈时将在空中发生爆炸,该现象称为空爆。探究空爆背后的物理原因,探索小天体的解体机制,对小天体撞击地球的灾害评估具有重要意义。

由于小天体进入大气的过程是一个复杂的多因素耦合的过程,况且地球上存有观测数据的小天体撞击事件十分有限,加之小天体轨道数据不够完善,使得人们对小天体空爆的物理机制仍未形成较为深刻的认识。

国内,柳森等人首先开展了小天体撞击地球大气层的相关研究,对小天体撞击地球大气层的物理问题及国内外发展现状、模型方法进行了综述[1]。在国外,Melosh、Collns等人提出了一系列小天体进入大气模型,并对该过程中的其他效应进行了研究[2-4]。科学家们发展了几种不同类型的模型[5,6],对小天体撞击地球大气过程进行了数值模拟,在模拟空爆高度方面取得较大进展。但对于能量沉积过程、空爆机理的研究仍未取得重大突破。

模型主要包括煎饼 (Pancake)模型[7]、碎片云 (FCM)模型[8]以及高保真爆炸流体动力学模型[9,10]。

煎饼模型将小天体假设为一个均质的连续体,并将小天体破碎、在空中发生爆炸的过程假设为一个连续的变形过程。该模型认为小天体在进入大气过程中,当其触发特定的飞行条件时,开始发生“破碎”,破碎表现为小天体受到挤压发生变形。[7]此类模型有明确的空爆高度的概念,当小天体截面半径增大到一定程度时,认为小天体发生空爆。

碎片云模型将小天体假设为球体,将其破碎分为几次进行,每次破碎,小天体会分裂为几个主要碎片,主要碎片仍为球体,其结构强度有所增加,可再次发生破碎,而碎片云则不再继续分裂。[8]该模型继承了煎饼模型中的变形方程,分裂次数变为多次,理论上不存在空爆高度的概念,小天体在撞击大气层过程中的变化主要依赖能量沉积进行描述。

爆炸流体动力学直接数值模拟针对特定的撞击事件或进入大气过程中的一个特定的方面进行研究。[6]通过直接模拟小天体撞击大气层过程中的流动、传热以及烧蚀相变的能量沉积过程,发现小天体首先发生煎饼状变形,变形较大 (大于两倍直径)时开始发生瑞利—泰勒不稳定性和开尔文—霍姆赫兹不稳定性,最后破碎为碎片、融化液滴、蒸汽的混合物,初步揭示了空爆的物理机制。但由于该方法计算过程复杂,时间成本高,较少用于撞击参数的影响研究。

三者对比,煎饼模型计算较为简便,所需时间成本较低,同时可以捕捉小天体撞击大气层过程中能量沉积、空爆解体等显著特征。该模型模拟过程中的能量变化,对于分析空爆的物理机制有很大帮助。

能量沉积包主要由阻力作用和烧蚀作用引起,本文将首次讨论两种作用对能量沉积的贡献,并由此对空爆机制进行分析。

本工作采用小天体撞击地球大气层的煎饼模型,以车里雅宾斯克事件为对照,以能量沉积过程为基础,分析、探究了小天体的空爆机制,初步得出了空爆的物理原因。

2 物理问题及数学模型

将小天体视为均质球体,主要输入参数有：小天体撞击大气层的初始速度、撞击角度、小天体的直径、小天体的密度。其中,撞击角度指小天体的撞击速度与当地水平线之间的夹角,如图1所示。

图1 撞击角示意图Fig.1 The impact angle

煎饼模型为一类模型的统称,其核心思想是：将小天体视为连续的变形体,将小天体的破碎过程视为连续体受挤压发生变形的过程,当变形条件被触发时,小天体在气动压力下发生变形。在变形过程中,小天体迎风截面的半径逐渐增大,当增大到一定程度时,认为小天体发生空爆。[7]

尽管核心思想相同,但不同的煎饼模型之间仍可能存在较大差异,比如小天体形状假设的差异、控制方程的不同、各参数之间的差别、终止判据的不同等等。[2,7,11]

Christopher F.Chyba等人于1993年提出了较为经典的煎饼模型,用于模拟1908年的通古斯卡事件,并列举了几组不同参数的小天体撞击地球大气层的模拟结果。[7]该模型假设小天体为底面直径与高相等的均质圆柱体,考虑地球曲率的影响,忽略地球自转的影响。在小天体撞击地球大气层的过程中,当小天体的驻点压力大于小天体本身的结构强度时,小天体将受到挤压而发生形变,形状逐渐趋于扁平,迎风截面也逐渐增大。

本文所用模型,在小天体飞行的基本控制方程、小天体截面半径变化方程、发生空爆的判断方式等方面,与Chyba等人提出的煎饼模型基本一致,同时也考虑了地球曲率的影响,忽略地球自转。不同的是,本文假设小天体为均质球体,并在一些参数上进行了细微的调整。小天体的速度V、撞击角γ、轨道半径r、路程中心角ψ、质量m随时间t的变化由公式 (1a)至 (1e)给出。

其中,Cd为大气阻力系数,本程序中取为1;ρ为即时大气密度,随高度发生变化;A为小天体的即时迎风截面积,由式A=πR2算得,其中R为小天体的即时迎风截面半径;μ则为地球引力常数,取3.986×105km3/s2;Cl为升力系数,由于小天体自身并无驱动力,故Cl取0;Ch为传热系数,取为0.1;Q为小天体的烧蚀热,与小天体的材料、结构等因素有关,由于本文假设小天体为均质球体,因此,在本文中,Q主要由小天体的材质决定。大气密度ρ采用幂指数模型计算：

其中,ρ0=1.2256 kg/m3为地平线处的大气密度,h为高度,在程序计算中为小天体轨道高度r与地球半径Re之差,单位为m;H为大气标高,为7254 m。

小天体进入大气层后,当驻点压力P大于小天体结构强度S0时,将开始发生变形,对应实际过程中的破碎过程。驻点压力由式 (3)求得：

在驻点压力的作用下,小天体受到挤压发生变形,呈现出越来越扁平的趋势,其迎风截面半径也逐渐增大,截面半径R的变化率随时间t的变化由式 (4)给出：

其中,ρm为小天体的密度。

小天体迎风截面直径D=2R,当D变为初始截面直径D0的N倍时,则认为小天体发生了空爆。其中,N为经验参数,通常在 5～10之间[2,7],本文在程序合理性验证中取为6,该取值下,程序模拟结果与该事件真实情况最为符合。在后续的主要撞击参数扫描中,N的取值保持6,不再变化。

2 模型合理性验证及空爆机制研究

模型通过对车里雅宾斯克事件进行模拟,以验证其合理性。

通过对真实视频材料等观测数据的分析[12-14],车里雅宾斯克事件的主要物理特征已经较为明确,为模型验证提供了丰富的对比材料。

据观察,车里雅宾斯克事件中,小天体于97km高空处进入大气,90km处开始形成冲击波,83km处初步产生尘埃及部分碎片,辐射峰值发生在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高度的范围内,辐照度与能量沉积曲线的趋势相一致,且均在30km高度左右达到峰值。

车里雅宾斯克事件主要参数如表1所示。

表1 车里雅宾斯克事件参数表Tab.1 Parameters of Chelyabinsk event

2.1 能量沉积计算方法

小天体在进入地球大气层时,部分动能通过大气阻力的作用、热烧蚀作用转移到大气层中,该过程为能量沉积。而能量沉积曲线用于描述不同高度处小天体在单位高度上损失的动能,反映小天体在撞击大气层过程中的能量变化及其随高度的变化情况,是探究空爆过程及其物理机制的重要依据。

能量沉积曲线采取如下计算方式：

式中,dE为两步计算之间小天体的动能减少量,dh为对应两步计算之间小天体的高度变化量,同样地,dm、dv也为对应两步计算之间质量、速度的减少量。各变化量取为原变化量的绝对值,以保证结果为正。

式 (5)也给出了能量沉积各项分量,通过简化,保留一阶项,略去二阶、三阶项,则能量沉积的主要贡献分为两项,如公式 (6)所示。该式主要用于2.3中拐点分析部分,不用于能量沉积的直接计算。

2.2 能量沉积曲线对比

通过仿真,得到了车里雅宾斯克事件发生空爆及之前的能量沉积曲线,如图4所示。

图4 车里雅宾斯克事件能量沉积曲线Fig.4 Energy deposition of Chelyabinsk event

图4中,黑色实线代表车里雅宾斯克事件真实的能量沉积曲线,基于对该事件的直接观测数据分析得出;绿色及蓝色的线条皆为碎片云模型模拟所得的能量沉积曲线[8];而红色粗体线条则为本程序模拟而得的能量沉积曲线。

从四条曲线的对比可以看出,本文得到的能量沉积曲线与真实曲线的整体趋势、重要节点基本一致,尤其在40km后的部分,二者几乎完全重合。在发生空爆的30km高度处,能量沉积值也与真实情况相一致。因此,可以认为,该方法能够初步模拟小天体撞击大气层的过程,其关于空爆高度及能量沉积的计算是可信的。

3 空爆机制分析

探究空爆机制,并不能仅着眼于发生空爆的短暂瞬间,而应当从小天体进入大气至发生空爆的整个过程出发。能量沉积曲线是记录该过程的重要依据,从该曲线出发,分析能量沉积过程中小天体的变化情况,对于揭露空爆机制有重要意义。

从图4可以看到,三类能量沉积曲线均存在一个较为明显的拐点,说明在该拐点处,小天体发生了某种变化,改变了能量沉积随高度变化的速率。

为探究该拐点所对应的物理过程,即导致该拐点产生的物理原因,根据公式 (6),分别将能量沉积的主要两项进行展示,第一项为对应速度变化,主要作用因素为大气阻力的影响,不妨称之为阻力项,其能量沉积贡献如图5所示;第二项为,对应质量变化,主要作用因素为烧蚀,此处称之为烧蚀项,其能量沉积贡献如图6所示。图7表示能量沉积总量。

图5 阻力项能量沉积曲线Fig.5 Energy deposition caused by resistance

图6 烧蚀项能量沉积曲线Fig.6 Energy deposition caused by ablation

图7 能量沉积总量Fig.7 Energy deposition

从图中可以看到,第一项与第二项以及总量均在39km处存在拐点,并呈现出相似的变化趋势。

就整体变化趋势而言,第二项与能量沉积总量的变化趋势更为接近,两条曲线几乎保持平行。

图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%,远大于阻力项能量沉积贡献的占比,为能量沉积的主要贡献。

结合图6、图7,可以判断,在39km左右,小天体烧蚀作用加剧,使得质量减少更加迅速,从而使该高度下小天体单位高度的动能沉积量变大。

推测烧蚀加剧的结果与小天体的变形有关。

图10为小天体的迎风截面半径与小天体高度之间的关系。

图10 小天体迎风截面半径随高度的变化Fig.10 Windward section radius variation with altitude

从图中能够看出,在小天体的高度到达39km之前,其截面半径并未发生改变。此时小天体也存在烧蚀作用,但由于驻点压强小于小天体屈服强度,因此烧蚀造成的质量损失主要表现为小天体纵向长度的减少,并未改变小天体的迎风截面。而在39km处,小天体的迎风截面半径开始发生变化,烧蚀作用也随之加剧,能量沉积进一步增加,最终导致了拐点的产生。

小天体开始变形的条件为驻点压强大于小天体本身的结构强度,因此,为了对拐点产生原因进行进一步验证,本文改变了小天体的结构强度,并呈现了总能量沉积曲线中拐点的位置变化情况,如图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时得到的能量沉积曲线,从图中可以明显看到,拐点发生在35km左右;图5正常情况下,拐点在39km左右;图12中曲线为时的能量沉积曲线,其拐点发生在约42km高度处。其中,S表示小天体结构强度,S0表示表1中基准结构强度。

小天体屈服强度高时,达到变形条件对应的高度较低,与拐点产生的高度规律相符。因此,小天体能量沉积曲线的拐点是小天体受到挤压开始发生变形所致,而变形引起的主要能量损失则源于烧蚀。

小天体在进入大气的过程中,首先在气动力、气动热的作用下,发生纵向烧蚀,在该阶段,由于大气密度较低,小天体的速度减少几乎可以忽略;随着大气密度的增加,小天体的驻点压强超过其自身的结构强度,小天体开始发生横向的变形,在该过程中,烧蚀作用加剧,小天体的能量以热、光等形式向外辐射;小天体的质量损失过大,其烧蚀将减弱。烧蚀作用最为剧烈时,对应小天体的外在表现即为空爆。空爆由剧烈的烧蚀作用引起。

空爆发生后,小天体在压力作用下发生破碎,破碎后的小天体碎片面质比增大,单个碎片在烧蚀作用下剥落、蒸发,形成图2中I、J所示的燃烧状烟雾。

4 总结与讨论

小天体撞击地球首先接触的便是地球大气层,空爆为小天体撞击地球大气层的重要物理特征,研究空爆对小天体防御有着重要意义。

本文基于Chyba等人于1993年提出的煎饼模型,对模型进行了适当的调整,对小天体撞击地球大气层中的空爆现象进行了分析,并研究了小天体撞击地球大气层的能量沉积过程,探究了空爆背后的物理机制。

首先以车里雅宾斯克事件为基础,对模型的合理性进行了验证,证明该煎饼模型程序能够用于小天体撞击地球大气层中空爆高度及能量沉积的计算,同时对能量沉积曲线中的拐点,对能量沉积过程进行了分析,并进一步提出了空爆对应的物理过程,发现空爆主要由剧烈的烧蚀作用引起。

尽管模型不同,但该结论与碎片云模型并不矛盾。碎片云模型中,虽没有明确的空爆概念,但仍存在几次小天体破碎的过程,使小天体的总表面积增大,烧蚀作用也会因此加剧。

此外,在煎饼模型中,发生空爆的判断标准为小天体的迎风截面半径变为初始半径的N倍,N为经验参数,一般取值在5～10之间,本文基于对车里雅宾斯克事件的模拟,将N取为6。

小天体撞击地球大气层的过程是十分复杂的,包括了声、光、热、冲击波等多种效应[2,17,18],本文对空爆这一重要现象进行了探讨,分析了小天体发生空爆的物理原因。若需对小天体撞击地球进行整体评估,则应当在此基础上逐步完善各种物理效应,形成一个整体的评估方案,以对行星防御提供更为有利的支持。