焦藝菲，程彬 *，寶音賀西, 2*

(1.清華大學(xué)航天航空學(xué)院，北京100084；2.內(nèi)蒙古工業(yè)大學(xué)，呼和浩特010051)

1 引言

太陽(yáng)系內(nèi)已發(fā)現(xiàn)的小天體超過(guò)數(shù)百萬(wàn)顆，其中近地小天體(包括近地小行星、彗星等)超過(guò)3萬(wàn)顆，據(jù)估計(jì)，直徑50m以上的近地小天體的真實(shí)數(shù)量可能超過(guò)10萬(wàn)顆，但目前發(fā)現(xiàn)的僅有2%，由于軌道與地球接近，這些小天體可能與地球相撞并引發(fā)區(qū)域級(jí)的災(zāi)難[1-3]。為了應(yīng)對(duì)近地小行星的潛在威脅，國(guó)內(nèi)外學(xué)者提出了眾多防御手段，包括動(dòng)能撞擊、核爆炸、長(zhǎng)期作用力偏轉(zhuǎn)等，其中動(dòng)能撞擊方案旨在利用航天器的高速撞擊來(lái)改偏小天體的軌道，進(jìn)而避免與地球相撞，是目前最有效可行的防御手段之一[4-10]。

在2022年中國(guó)航天日大會(huì)上，國(guó)家航天局表示，我國(guó)將著手組建近地小行星防御系統(tǒng)，并計(jì)劃在2025年左右實(shí)施一次對(duì)某顆有威脅的小行星既進(jìn)行抵近觀(guān)測(cè)，又實(shí)施就近撞擊改變其軌道的技術(shù)實(shí)驗(yàn)[11]。小行星撞擊防御任務(wù)需要對(duì)高速撞擊過(guò)程進(jìn)行盡可能精確的觀(guān)測(cè)，地基望遠(yuǎn)鏡和繞地的觀(guān)測(cè)衛(wèi)星難以實(shí)現(xiàn)，而使用探測(cè)器接近觀(guān)測(cè)能夠?qū)ψ矒魹R射過(guò)程和撞擊的直接結(jié)果進(jìn)行記錄，但是高速撞擊濺射的碎片可能對(duì)近距離的探測(cè)器造成威脅。此前國(guó)外航天機(jī)構(gòu)也開(kāi)展過(guò)多次深空的高速撞擊和接近觀(guān)測(cè)任務(wù)，如2005年的“深度撞擊”(Deep Impact)探測(cè)器高速撞擊彗星9P/Tempel，采用近距離飛越觀(guān)測(cè)的方式，最小飛越距離為500km，并成功觀(guān)測(cè)到撞擊產(chǎn)生的大量發(fā)光濺射塵埃，但由于觀(guān)測(cè)時(shí)間窗口較短且觀(guān)測(cè)距離較遠(yuǎn)，成像質(zhì)量較差且無(wú)法觀(guān)測(cè)到最終的撞擊坑[12]；2009年的“月坑觀(guān)測(cè)與感知衛(wèi)星”(LCROSS)月球撞擊任務(wù)使用火箭上面級(jí)撞擊月球表面，在撞擊后4min內(nèi)，使用飛船對(duì)噴射羽流進(jìn)行近距離觀(guān)測(cè)并傳輸數(shù)據(jù)，隨后飛船再次撞擊月球[13]；2019年“隼鳥(niǎo)”2(Hayabusa 2)探測(cè)器在小行星“龍宮”(Ryugu)表面進(jìn)行了小型人工撞擊試驗(yàn)，在撞擊側(cè)向約1km位置放置小型相機(jī)進(jìn)行觀(guān)測(cè)，探測(cè)器則逃離到小行星背面約20km處以免被濺射碎片撞擊[14]；2022年9月，“雙小行星重定向測(cè)試”(DART)探測(cè)器成功撞擊Didymos雙小行星系統(tǒng)的子星Dimorphos，攜帶的立方星LICIACube在撞擊前分離，以55km的距離飛越觀(guān)測(cè)撞擊濺射情況，但目前發(fā)布的撞擊圖像無(wú)法用于精確評(píng)估軌道偏轉(zhuǎn)效果[15-17]；此外，計(jì)劃于2024年發(fā)射的Hera任務(wù)將再次造訪(fǎng)Didymos雙星并對(duì)DART的撞擊結(jié)果進(jìn)行全面評(píng)估[18]。

相比于飛越觀(guān)測(cè)的方式，伴飛、繞飛的方式可以對(duì)小天體的撞擊過(guò)程進(jìn)行長(zhǎng)時(shí)間、近距離的觀(guān)測(cè)[19,20]，因此能夠有效觀(guān)測(cè)到高速撞擊濺射、撞擊坑形成、低速碎片在不規(guī)則的弱引力場(chǎng)中演化的整個(gè)過(guò)程，既是對(duì)撞擊防御結(jié)果的直接評(píng)估，也對(duì)研究小天體的材料性質(zhì)、內(nèi)部結(jié)構(gòu)和理解小天體的形成與演化過(guò)程有重要價(jià)值。針對(duì)接近觀(guān)測(cè)任務(wù)的安全性需求，本文使用完全自主知識(shí)產(chǎn)權(quán)的小天體高速撞擊數(shù)值模擬軟件THU-SPHSOL，通過(guò)數(shù)值仿真研究了均質(zhì)巖石小天體的高速撞擊濺射過(guò)程，以及不同撞擊角度對(duì)濺射結(jié)果的影響，分析了高速濺射碎片的空間分布情況和觀(guān)測(cè)相機(jī)的安全區(qū)域，研究結(jié)果將為我國(guó)即將開(kāi)展的小行星撞擊防御任務(wù)提供重要參考。

2 數(shù)值仿真模型

小天體的高速撞擊和濺射過(guò)程涉及強(qiáng)非線(xiàn)性的力-熱耦合、流-固耦合，材料在高壓高應(yīng)變率的極端條件下發(fā)生大變形、斷裂和破碎，傳統(tǒng)的基于網(wǎng)格的數(shù)值方法難以求解，因此，基于無(wú)網(wǎng)格的光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics, SPH)方法[21]，我們開(kāi)發(fā)了完全自主知識(shí)產(chǎn)權(quán)的THU-SPHSOL軟件[22]，軟件面向航天動(dòng)力學(xué)和行星科學(xué)等領(lǐng)域的高速撞擊問(wèn)題，能夠?qū)崿F(xiàn)千萬(wàn)粒子量級(jí)的大規(guī)模高效計(jì)算，本節(jié)主要對(duì)軟件采用的物理模型進(jìn)行介紹。

2.1 拉格朗日描述的控制方程

由于高速撞擊產(chǎn)生的沖擊壓力遠(yuǎn)遠(yuǎn)大于材料的強(qiáng)度，材料的行為特性接近無(wú)粘的可壓縮流體，使用拉格朗日描述的連續(xù)介質(zhì)守恒方程作為控制方程[23]，包括連續(xù)性方程、動(dòng)量方程、能量方程和運(yùn)動(dòng)方程：

(1)

(2)

(3)

(4)

2.2 狀態(tài)方程

狀態(tài)方程是高速撞擊過(guò)程中控制材料行為特性的關(guān)鍵部分，描述了材料密度ρ、壓力p和比內(nèi)能e之間的關(guān)系，在小天體高速撞擊場(chǎng)景中，Tillotson狀態(tài)方程是最常使用的狀態(tài)方程之一，在低壓/低能量時(shí)滿(mǎn)足Hugoniot關(guān)系，在高壓(高能量)時(shí)趨于Thomas-Fermi極限[24]。當(dāng)材料處于壓縮(ρ>ρ0)或冷膨脹(ρ<ρ0且e

(5)

當(dāng)材料處于熱膨脹狀態(tài)時(shí)，

(6)

對(duì)于中間狀態(tài)，則根據(jù)比內(nèi)能e進(jìn)行線(xiàn)性插值，

(7)

狀態(tài)方程還決定了材料的聲速或縱波波速c

(8)

式中：G為剪切模量，對(duì)于使用Tillotson狀態(tài)方程描述的材料，根據(jù)式(5)～(8)計(jì)算材料聲速[25]。

2.3 強(qiáng)度模型

高速撞擊的沖擊壓力隨著時(shí)間和距離發(fā)生衰減，此時(shí)壓力項(xiàng)和偏應(yīng)力項(xiàng)將同時(shí)主導(dǎo)固體材料的響應(yīng)特性，因此屈服失效和斷裂等與強(qiáng)度相關(guān)的特性必須加以考慮，通常按照理想彈塑性模型進(jìn)行處理。在彈性階段，假設(shè)偏應(yīng)力率和偏應(yīng)變率成比例，即滿(mǎn)足(增量型)胡克定律，并考慮剛體轉(zhuǎn)動(dòng)的影響，此時(shí)偏應(yīng)力率的更新形式為：

(9)

材料是否進(jìn)入塑性階段則由屈服條件決定，對(duì)于小行星等地質(zhì)材料，采用Lund模型，考慮了內(nèi)聚力和內(nèi)摩擦力對(duì)屈服強(qiáng)度的影響[26,27]，

(10)

材料的損傷程度描述了材料抗剪切和抗拉伸能力的降低，通常使用一維標(biāo)量D作為損傷的度量，且0≤D≤1，其中D=0對(duì)應(yīng)完好無(wú)損的材料，D=1對(duì)應(yīng)完全失效的材料。對(duì)于小行星等脆性地質(zhì)材料，損傷的形式主要是內(nèi)部缺陷導(dǎo)致的拉伸損傷，因此采用基于Weibull形式缺陷分布的損傷累積模型[28,29]。根據(jù)Grady Kipp模型[30]，材料損傷對(duì)抗拉伸能力的削弱通過(guò)壓力項(xiàng)的修正實(shí)現(xiàn)，

p′=(1-D)p,p<0

(11)

損傷對(duì)抗剪切能力的削弱則通過(guò)修正屈服強(qiáng)度來(lái)實(shí)現(xiàn)，在Lund模型中，材料完全失效時(shí)屈服強(qiáng)度與壓力成正比，即Yd=μdp，此時(shí)損傷修正的Lund強(qiáng)度為

Y=(1-D)Yi+DYd

(12)

3 仿真結(jié)果與討論

3.1 數(shù)值仿真算例

面向我國(guó)即將開(kāi)展的小行星撞擊防御任務(wù)，典型的撞擊場(chǎng)景為質(zhì)量500kg的航天器以5000m/s的相對(duì)速度垂直撞擊目標(biāo)小行星，但考慮到航天器的撞擊位置偏差以及小行星的表面地形，真實(shí)撞擊場(chǎng)景下撞擊方向可能偏離撞擊表面的法向，撞擊濺射分布也可能完全不同，因此分別設(shè)置撞擊方向偏離表面法向0°、15°和30°三組工況，對(duì)應(yīng)算例(1)～(3)并假設(shè)撞擊局部的表面平坦，研究不同撞擊角度對(duì)結(jié)果的影響。航天器采用密度2700kg/m3的鋁材料，使用Von Mises模型和Tillotson狀態(tài)方程描述，小行星則采用密度2700kg/m3的玄武巖材料，不考慮孔隙度，使用Lund強(qiáng)度模型、損傷累積模型和Tillotson狀態(tài)方程描述，材料參數(shù)見(jiàn)表1。為提高離散精度并減小計(jì)算量，將計(jì)算域限制為撞擊點(diǎn)周?chē)睆?0m的范圍，SPH粒子初始間距取0.132m，此時(shí)撞擊器和小行星分別表示為81個(gè)和300萬(wàn)個(gè)SPH粒子，每個(gè)粒子的質(zhì)量為6.21kg，撞擊仿真總時(shí)長(zhǎng)取20ms。

表1 材料參數(shù)

3.2 結(jié)果與討論

以垂直撞擊為例，算例1的數(shù)值仿真結(jié)果如圖1所示，顯然，20ms的瞬時(shí)撞擊坑并不是最終形態(tài)，整個(gè)成坑過(guò)程可能持續(xù)幾秒到幾十秒，但20ms對(duì)于研究撞擊濺射的較高速碎片是足夠的，此時(shí)圖1中大部分區(qū)域都已經(jīng)完全失效，沖擊壓力也相比撞擊前期大幅衰減，此時(shí)粒子的運(yùn)動(dòng)速度對(duì)后續(xù)的演化過(guò)程起主導(dǎo)作用，具有向外速度分量的粒子(圖中黃色箭頭)將繼續(xù)濺射，而具有向內(nèi)速度分量的粒子(圖中藍(lán)色箭頭)在沖擊作用下向內(nèi)壓縮，因此不會(huì)發(fā)生濺射，這也與Z-model的模型預(yù)測(cè)相符[32,33]。

圖1 垂直撞擊成坑的濺射過(guò)程，箭頭方向?yàn)榱Ｗ铀俣确较?，箭頭顏色為沿撞擊反向(x軸方向)的速度分量，其中黃色為向外濺射，藍(lán)色為向內(nèi)壓縮；圖中標(biāo)注了撞擊后20ms瞬時(shí)的撞擊坑，以及根據(jù)速度場(chǎng)預(yù)測(cè)的最終撞擊坑(直徑約16m，深度約6m)

盡管撞擊濺射的粒子仍然受到小天體、太陽(yáng)和其他行星的引力作用，太陽(yáng)光壓、太陽(yáng)磁場(chǎng)產(chǎn)生的電磁力甚至粒子之間的接觸碰撞也都將對(duì)粒子的運(yùn)動(dòng)產(chǎn)生影響[34,35]，但對(duì)于撞擊過(guò)程的短期(約幾十到數(shù)百秒)接近觀(guān)測(cè)，仍然可以假設(shè)所有濺射粒子都保持撞擊仿真結(jié)束時(shí)的速度做勻速直線(xiàn)運(yùn)動(dòng)。我們篩選所有速度大小超過(guò)某個(gè)閾值且具有向外的速度分量的粒子作為濺射粒子集合，分別取1m/s和10m/s作為速度閾值(遠(yuǎn)大于百米量級(jí)小天體的表面逃逸速度)，為了對(duì)比不同撞擊角度的濺射結(jié)果，根據(jù)濺射速度在yz平面的投影，進(jìn)一步篩選速度投影方向與y軸夾角在±15°以?xún)?nèi)的粒子進(jìn)行分析，如圖2所示，對(duì)于斜撞擊情況，篩選區(qū)域粒子的濺射速度具有明顯的不對(duì)稱(chēng)性，并對(duì)應(yīng)了不同撞擊角度的濺射分布特征。根據(jù)不同的篩選條件，表2統(tǒng)計(jì)了濺射粒子的數(shù)量分布，可以看出，垂直撞擊的濺射粒子最多，隨著撞擊角度的偏離，濺射粒子數(shù)量減少，撞擊坑也會(huì)更小；速度10m/s以上的粒子數(shù)量約為速度1m/s以上粒子數(shù)量的6%～7%，因此絕大部分濺射粒子都是低速的，但對(duì)于接近觀(guān)測(cè)的探測(cè)器，10m/s以上的濺射粒子具有較大的動(dòng)能且覆蓋范圍更遠(yuǎn)，因此相比于低速粒子更具威脅性。

表2 濺射粒子統(tǒng)計(jì)

圖2 撞擊坐標(biāo)系的定義和濺射粒子的篩選區(qū)域

對(duì)于濺射速度方向在所篩選區(qū)間的粒子，分別考慮1m/s和10m/s兩個(gè)速度閾值，不同撞擊角度(算例1～3)的濺射方向和空間分布如圖3所示?？梢钥闯?，對(duì)于垂直于小行星表面的撞擊，速度1m/s以上的濺射粒子與撞擊反向的夾角主要分布在20°～70°之間，而速度10m/s以上的濺射粒子角度主要分布在30°～60°，其中分布在40°～50°的粒子數(shù)量最多，此前的文獻(xiàn)研究和撞擊實(shí)驗(yàn)[36,37]同樣指出，垂直撞擊成坑過(guò)程中，不同位置發(fā)射的濺射物速度與水平方向角度分布在40°～55°之間，但整體的濺射角度接近45°，并使用標(biāo)度律來(lái)描述濺射位置、速度大小和累積質(zhì)量之間的關(guān)系，當(dāng)然斜撞擊情況下濺射角度還會(huì)隨著方位角發(fā)生變化[38,39]；如果只考慮對(duì)接近觀(guān)測(cè)造成威脅較大的、速度10m/s以上的濺射粒子，并假設(shè)這些粒子在短期的演化過(guò)程中做勻速直線(xiàn)運(yùn)動(dòng)，在隨后的1～100s時(shí)間內(nèi)，濺射粒子的空間位置也基本分布在以撞擊點(diǎn)為中心的錐形范圍內(nèi)，錐角區(qū)間與速度角度區(qū)間基本相同，都是30°～60°，且濺射分布關(guān)于撞擊面法向(對(duì)于垂直撞擊也是撞擊方向)軸對(duì)稱(chēng)。隨著撞擊方向的偏離，濺射粒子在撞擊點(diǎn)兩側(cè)的分布呈現(xiàn)出明顯的不對(duì)稱(chēng)性，下坡側(cè)的濺射粒子數(shù)量遠(yuǎn)多于上坡側(cè)，因此在隨后的1～100s時(shí)間內(nèi)，濺射粒子的空間位置關(guān)于撞擊方向不再旋轉(zhuǎn)對(duì)稱(chēng)，但基本上仍是關(guān)于撞擊面法向的軸對(duì)稱(chēng)分布，即使兩側(cè)的濺射粒子數(shù)量不同。

圖3 不同撞擊場(chǎng)景的濺射角度分布，其中(a)～(c)為速度超過(guò)1m/s的部分，(d)～(f)為速度超過(guò)10m/s的部分，極坐標(biāo)系的橫向坐標(biāo)為濺射速度與撞擊反向(x軸方向)的夾角，徑向坐標(biāo)為每個(gè)角度區(qū)間的粒子數(shù)量，紅色箭頭表示撞擊方向，綠色線(xiàn)表示撞擊表面；不同撞擊場(chǎng)景的濺射空間分布，包括(g)～(i)，只考慮速度超過(guò)10m/s的濺射粒子，假設(shè)粒子從20ms(撞擊仿真的結(jié)束時(shí)刻)開(kāi)始做勻速直線(xiàn)運(yùn)動(dòng)，極坐標(biāo)系的橫向坐標(biāo)為粒子的空間位置與x軸反向的夾角，徑向坐標(biāo)為濺射粒子遠(yuǎn)離撞擊點(diǎn)的距離

在真實(shí)的小行星撞擊防御任務(wù)中，航天器的撞擊角度不僅取決于撞擊點(diǎn)的選取和小行星的形狀，還會(huì)受到導(dǎo)航和制導(dǎo)精度的影響[40,41]，因此真實(shí)撞擊角度可能分布在某個(gè)范圍內(nèi)。考慮直徑50m的球形小行星，假設(shè)標(biāo)稱(chēng)撞擊點(diǎn)指向球心，那么最終撞擊點(diǎn)可能是以標(biāo)稱(chēng)撞擊點(diǎn)為中心、某個(gè)距離范圍內(nèi)的任意位置：如果撞擊位置誤差在半徑6.5m范圍內(nèi)，對(duì)應(yīng)的撞擊角度總是小于15°，此時(shí)絕大部分撞擊濺射粒子分布在與撞擊反向(x軸方向)夾角20°～80°的錐形范圍內(nèi)，并且將在100s的時(shí)間內(nèi)到達(dá)1～10km外的區(qū)域，觀(guān)測(cè)器應(yīng)避開(kāi)該角度區(qū)間，或者在更遠(yuǎn)處觀(guān)測(cè)一段時(shí)間后重新尋找安全區(qū)域；如果撞擊位置誤差在半徑12.5m范圍內(nèi)，對(duì)應(yīng)的撞擊角度不超過(guò)30°，此時(shí)濺射粒子可能分布在與撞擊反向夾角0°～90°的半空間范圍內(nèi)，但如果觀(guān)測(cè)相機(jī)躲避到撞擊的背面，將難以拍攝到撞擊濺射過(guò)程，此時(shí)建議觀(guān)測(cè)相機(jī)放置在約10km外的區(qū)域，并在約100s后離開(kāi)該區(qū)域。

4 結(jié)論

本文結(jié)合我國(guó)即將開(kāi)展的小行星撞擊防御任務(wù)和接近觀(guān)測(cè)的安全性需求，使用小天體高速撞擊數(shù)值模擬軟件THU-SPHSOL對(duì)撞擊過(guò)程進(jìn)行數(shù)值仿真研究，介紹了軟件主要采用的物理模型和材料參數(shù)，數(shù)值仿真研究了均質(zhì)巖石小天體的高速撞擊濺射過(guò)程，以及不同撞擊角度對(duì)濺射結(jié)果的影響。研究表明，垂直撞擊的高速濺射粒子主要分布在與撞擊反向夾角為30°～60°的錐形范圍內(nèi)，但撞擊方向的偏離會(huì)引起濺射粒子數(shù)量和角度分布的不對(duì)稱(chēng)性，在真實(shí)撞擊場(chǎng)景下，觀(guān)測(cè)相機(jī)的安全區(qū)域范圍還將取決于撞擊位置的誤差范圍。本文的研究結(jié)果將為我國(guó)即將開(kāi)展的小行星撞擊防御任務(wù)提供重要參考。