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

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

1 引言

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

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

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

2 數(shù)值仿真模型

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

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

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

(1)

(2)

(3)

(4)

2.2 狀態(tài)方程

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

(5)

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

(6)

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

(7)

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

(8)

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

2.3 強度模型

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

(9)

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

(10)

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

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

(11)

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

Y=(1-D)Yi+DYd

(12)

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

3.1 數(shù)值仿真算例

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

表1 材料參數(shù)

3.2 結(jié)果與討論

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

圖1 垂直撞擊成坑的濺射過程，箭頭方向為粒子速度方向，箭頭顏色為沿撞擊反向(x軸方向)的速度分量，其中黃色為向外濺射，藍色為向內(nèi)壓縮；圖中標(biāo)注了撞擊后20ms瞬時的撞擊坑，以及根據(jù)速度場預(yù)測的最終撞擊坑(直徑約16m，深度約6m)

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

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

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

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

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

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

4 結(jié)論

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