廖祜明，龔自正，宋光明，樊宗岳，楊宏濤，黎波*

(1.北京航空航天大學(xué)能源與動力工程學(xué)院，北京 102206；2.北京衛(wèi)星環(huán)境工程研究所，北京 100094；3.北京大學(xué)工學(xué)院，北京 100871；4.云翼超算(北京)軟件科技有限公司，北京 100027)

1 引言

地球自誕生以來就不斷受到小行星(Asteroid)的撞擊，這是人類生存面臨的重大潛在威脅之一。如何減緩或避免小行星撞擊地球已經(jīng)成為學(xué)術(shù)界和國際社會關(guān)注的熱點問題[1]。在天文學(xué)上，定義軌道在距離太陽1.3AU(1AU=1.496×108km)范圍內(nèi)，距離地球軌道最小距離在0.3AU范圍內(nèi)的小行星為近地小行星(Near-Earth Asteroid, NEA)，截至2022年9月26日，已發(fā)現(xiàn)NEA 29901顆[2]。當(dāng)小行星與地球距離0.05AU時，有可能被地球引力俘獲，改變運行軌道與地球相撞，而直徑140m以上的小行星，撞擊地球的威力足以造成區(qū)域性災(zāi)難，毀掉一個中等大小的國家，因此把距離地球軌道最小距離在0.05AU范圍內(nèi)、直徑大于140m的小行星定義為具有潛在碰撞威脅的小行星(Potentially Hazardous Asteroid, PHA)[3]。目前已發(fā)現(xiàn)的PHA有2289顆[2]，是防御任務(wù)的重點目標(biāo)。

在能夠?qū)HA提前預(yù)警的前提下，根據(jù)預(yù)警時間長短以及目標(biāo)小行星尺寸的不同，將小行星分裂成碎片或者改變小行星軌道是避免其撞擊地球的兩種基本方式。具體來說：一是對于預(yù)警時間充裕且尺寸較小的PHA，可采用長期作用力來改變小行星軌道(包括太空拖船、引力拖車、質(zhì)量推進(jìn)器、激光燒蝕、表面噴漆、離子束等)；二是對于尺寸較小、預(yù)警時間較短或尺寸較大且預(yù)警時間較長的PHA，可采用動能撞擊技術(shù)直接改變小行星運行軌道；而對于預(yù)警時間很短且直徑大于500m的PHA，目前研究認(rèn)為摧毀或部分摧毀小行星的核爆炸(表面爆炸、對峙爆炸以及穿透爆炸)是進(jìn)行防御的有效手段[4-8]。

2022年9月27日，在太空經(jīng)歷了長達(dá)10個月的飛行后，執(zhí)行美國國家航空航天局(NASA)“雙小行星重定向測試”(DART)任務(wù)的衛(wèi)星成功撞擊了直徑約163m的Dimorphos小行星，這是全球首次行星防御技術(shù)的演示，也是NASA首次嘗試在太空中對小行星的運行進(jìn)行偏轉(zhuǎn)，作為NASA整體行星防御戰(zhàn)略的一部分。DART對小行星Dimorphos的撞擊展示了一種可行的防御技術(shù)，觀測DART任務(wù)撞擊效果的“小行星撞擊監(jiān)視器”(HERA)與“小行星撞擊偏轉(zhuǎn)評估計劃”(AIDA)[9]是正在實施的針對動能撞擊技術(shù)的在軌驗證試驗，將進(jìn)一步評估動能撞擊對于Dimorphos公轉(zhuǎn)速度和軌道的影響[4]。預(yù)先計算顯示此次撞擊將給Dimorphos的運行速度造成0.2mm/s的改變，在長時間的飛行過程中其軌跡將達(dá)到一個可觀的累積偏轉(zhuǎn)量。由于試驗研究存在周期長、成本高與測量技術(shù)受限等困難與挑戰(zhàn)，數(shù)值計算成為了輔助預(yù)測防御措施有效性的高效手段之一。

在動能撞擊數(shù)值仿真領(lǐng)域，張韻[10]等采用物質(zhì)點法(Material Point Method，MPM)進(jìn)行了鋁彈超高速撞擊碎石堆結(jié)構(gòu)(rubble-pile structure)和完整結(jié)構(gòu)(monolithic structure)兩種不同小行星材料模型的動能撞擊效果評估；Jutzi[11]等提出了一種改進(jìn)的SPH方法，模擬了小行星災(zāi)難性破壞的過程，結(jié)果展示了材料屈服極限、剪切強(qiáng)度和孔隙率在碰撞中的相對重要性；姜宇[12]等將小行星Bennu視作由球型顆粒聚成的碎石堆，采用軟球離散元耦合引力N體模型對其進(jìn)行了碰撞數(shù)值模擬；Raducan[13-15]等以DART任務(wù)作為研究案例，使用基于ALE方法的iSALE軟件，對類DART衛(wèi)星撞擊小行星作了大量的數(shù)值模擬，總結(jié)了衛(wèi)星幾何形狀、撞擊速度與角度和小行星材料屬性等參數(shù)對撞擊坑坑徑和撞擊動量轉(zhuǎn)換效率的影響；馬鑫[16]等采用High-explosive-burn材料模型和JWL狀態(tài)方程描述聚能爆炸成型彈丸(EFP)，采用Johnson-Cook材料模型和Gruneisen狀態(tài)方程描述混凝土球體靶板，利用有限元方法進(jìn)行深空撞擊載荷仿真，形成了撞擊速度、靶板強(qiáng)度、靶板密度、靶板體積不同參數(shù)條件下的撞擊坑坑徑變化規(guī)律。在核爆仿真領(lǐng)域，李毅[17]等在歐拉型沖擊動力學(xué)仿真軟件NTS中加入能量源，模擬了核爆裝置在不同深度爆炸對小行星產(chǎn)生的偏轉(zhuǎn)與破壞效應(yīng)；湯文輝[18]等采用有限元方法以及能夠描述氣化反沖的PUFF物態(tài)方程進(jìn)行了核爆炸偏轉(zhuǎn)小行星數(shù)值模擬；Patrick K[19]等使用Spheral++[20]軟件中的自適應(yīng)SPH算法和N體引力算法模擬了核爆過程和碎片重組過程，研究結(jié)果表明在實際操作中需要合理預(yù)估核爆沖擊能量，權(quán)衡風(fēng)險系數(shù)?？梢?，國內(nèi)外研究人員在小行星防御的數(shù)值仿真領(lǐng)域已開展相關(guān)研究工作并取得了一定的成果。

然而，由于小行星防御及與小行星撞擊地球涉及的極高速、極高溫以及強(qiáng)耦合的極端苛刻環(huán)境，如何高效高精度預(yù)測材料在高溫、高壓、高應(yīng)變率等極限熱-力-化學(xué)條件下的動態(tài)響應(yīng)機(jī)制一直是數(shù)值仿真模擬的難點之一。目前圍繞小行星防御的數(shù)值計算在基礎(chǔ)理論、計算方法、材料模型與高性能計算方面都還面臨許多困難與挑戰(zhàn)?；A(chǔ)理論方面，由于在核爆、動能撞擊、激光驅(qū)動等將小行星分裂成碎片或改變軌道的防御措施中，以及小行星進(jìn)入地球大氣與撞擊地球表面的過程中，大量存在著燒蝕、解體、爆炸、火球、撞擊成坑、反濺碎片云、地震等一系列復(fù)雜的物理化學(xué)和力學(xué)現(xiàn)象，需要理論模型能夠綜合考慮其中的熱流固耦合、熱力化學(xué)耦合等多物理場強(qiáng)耦合效應(yīng)；計算方法方面，目前開展的工作主要采用基于網(wǎng)格和無網(wǎng)格方法兩大類。基于網(wǎng)格的數(shù)值方法中，傳統(tǒng)拉氏有限元法受網(wǎng)格畸變困擾，難以處理超大變形問題；歐氏方法具有網(wǎng)格不變形的固有優(yōu)勢，但如何精確模擬與加載速率和加載路徑相關(guān)的材料響應(yīng)、準(zhǔn)確跟蹤物質(zhì)界面和固液氣動態(tài)相變?nèi)匀皇秦酱鉀Q的問題，同時在計算動態(tài)裂紋擴(kuò)展與多體碰撞等非線性行為時還需進(jìn)一步研究；任意拉格朗日-歐拉(Arbitrary Lagrangian Eulerian, ALE)[21-23]法在一定程度上緩解了網(wǎng)格畸變帶來的困難，但ALE需要復(fù)雜的網(wǎng)格管理和映射技術(shù)，同時也面臨超大變形時網(wǎng)格畸變的問題[24]。無網(wǎng)格法由于不需要進(jìn)行網(wǎng)格離散及采用高階插值函數(shù)，在大變形問題上應(yīng)用更為廣泛。常用的無網(wǎng)格方法有光滑粒子流體動力學(xué)法(Smoothed Particle Hydrodynamics, SPH)[25,26]、再生核粒子法(Reproducing Kernel Particle Method, RKPM)[27]、軟球離散元法(Soft-sphere Discrete Element)[28]以及物質(zhì)點法(Material Point Method, MPM)[29]等，這些方法在超高速動能毀傷和空間碎片防護(hù)研究中發(fā)揮了巨大作用，但在對多物理場強(qiáng)耦合的處理、穩(wěn)定性和效率等方面仍然存在不足[30]。材料模型方面，與航天飛行器經(jīng)過人類精心設(shè)計的材料和結(jié)構(gòu)不同，小行星的材料和結(jié)構(gòu)是自然形成的，不僅呈現(xiàn)出各向異性的特征，而且充滿孔隙和裂紋，因而對典型材料模型和高狀態(tài)相變數(shù)據(jù)都提出了更多的需求[1]。

對小行星防御場景的高效高精度仿真是一個對基礎(chǔ)理論、計算方法、材料模型不斷豐富完善的過程，其對數(shù)值計算方法的基礎(chǔ)理論框架提出了具有通用性和可擴(kuò)展性的需求，從而實現(xiàn)不斷融合新的物理化學(xué)現(xiàn)象、新的材料模型和更加真實的邊界條件，形成小行星防御專用仿真平臺，更加真實高效地仿真各種小行星防御場景。最優(yōu)運輸無網(wǎng)格方法(Optimal Transportation Meshfree, OTM)[31]是一種針對動態(tài)沖擊問題提出的顯式增量更新拉格朗日無網(wǎng)格方法。OTM采用基于變分原理的多物理場自主強(qiáng)耦合理論框架，為了求解熱流固耦合問題，在OTM基礎(chǔ)上通過引入能量守恒方程和熱流固耦合變分本構(gòu)，擴(kuò)展形成熱力強(qiáng)耦合最優(yōu)運輸無網(wǎng)格方法(Hot Optimal Transpor-tation Meshfree, HOTM)[32]，進(jìn)一步地通過結(jié)合分布式多進(jìn)程并行與共享內(nèi)存多線程并行策略，實現(xiàn)了大規(guī)模雙層混合并行最優(yōu)運輸無網(wǎng)格方法(massively parallel OTM method, pOTM)[33]，顯著提高了計算效率。在材料失效方面，在OTM系列方法的基礎(chǔ)上，采用基于能量準(zhǔn)則的裂紋擴(kuò)展算法(EigenFracture)[34]，以材料能量釋放率作為失效判據(jù)，克服傳統(tǒng)裂紋擴(kuò)展數(shù)值方法裂紋擴(kuò)展路徑網(wǎng)格相關(guān)、不收斂、無法清晰表征材料真實變形和失效物理機(jī)制的缺點。通過有機(jī)融合這些方法形成的ESCAAS高度非線性多物理場強(qiáng)耦合極限力學(xué)仿真平臺，如圖 1所示，實現(xiàn)在統(tǒng)一的框架下求解材料超大變形、熔化、氣化、沖擊爆炸、流固與熱流固耦合、自由表面、多體接觸、碎裂、層裂與碎片云等復(fù)雜物理現(xiàn)象的多物理場、多尺度強(qiáng)耦合問題，為小行星防御仿真預(yù)測提供了具有潛力的解決方案。

圖1 基于OTM的極限力學(xué)仿真理論框架

本文將對OTM/HOTM極限力學(xué)仿真理論框架進(jìn)行介紹，并探討其在隕石超高速撞擊成坑、動能撞擊、激光燒蝕驅(qū)動、核爆攔截等偏移小行星軌道以及摧毀小行星等防御措施中的潛在應(yīng)用方向和相關(guān)案例，為小行星防御提供相關(guān)理論與技術(shù)支撐。

2 OTM/HOTM極限力學(xué)仿真理論框架

本節(jié)簡要表述極限力學(xué)問題的理論基礎(chǔ)與基于OTM/HOTM方法的數(shù)值求解框架。

2.1 拉格朗日描述控制方程

不失一般化，以一般帶耗散熱力耦合問題中Ω0?R3的連續(xù)物體進(jìn)行討論，其運動過程可表述為隨時間變化的變形映射φ:Ω0×[t0,t]→R3，其中Ω0代表參考構(gòu)型，[t0,t]為運動過程總時間。X∈Ω0表示材料質(zhì)點在參考構(gòu)形中的位置，x=φ(X,t)表示材料質(zhì)點在現(xiàn)時構(gòu)形Ωt=φ(Ω0,t)的位置，則拉格朗日描述下質(zhì)量守恒、動量守恒與能量守恒方程為：

ρ(x)J=ρ0

(1)

(2)

(3)

2.2 多物理場強(qiáng)耦合變分框架

結(jié)合熱力耦合變分原理[35]，建立與運動控制方程、初始條件與邊界條件等效的變分結(jié)構(gòu)。在連續(xù)介質(zhì)力學(xué)中，與速率問題相對應(yīng)的變分結(jié)構(gòu)可表示為：

(4)

(5)

式中：A表示Helmholtz自由能，φ*為粘性能量耗散，ψ*為塑性變形能量耗散，χ為與熱傳導(dǎo)相關(guān)的能量耗散。矢量G=-(?)0T)/T描述熱能狀態(tài)，?0表示在參考坐標(biāo)系下的偏導(dǎo)。Θ(F,T,Z)=?U/?N表示材料內(nèi)部溫度場由內(nèi)能推導(dǎo)得到，并在熱平衡狀態(tài)下滿足Θ=T?？梢姷刃菽苡上到y(tǒng)存儲的能量(Helmholtz自由能)和系統(tǒng)耗散的能量(熵、粘性耗散、塑性耗散、熱傳導(dǎo)耗散)所組成。其中，Helmholtz自由能可分解為：

(6)

2.3 最優(yōu)運輸無網(wǎng)格數(shù)值離散

在有限元法中，每個單元由若干節(jié)點組成特定形狀，各節(jié)點與單元有固定的從屬連接關(guān)系，每個單元攜帶了一定量的材料質(zhì)量并具有相應(yīng)的體積，材料的響應(yīng)在單元積分點上求解，插值則在單元節(jié)點上進(jìn)行。與有限元法類似，OTM方法中采用物質(zhì)點xp與節(jié)點xa相結(jié)合的方式進(jìn)行空間離散，如圖 2所示，其中下標(biāo)a表示節(jié)點索引，p表示物質(zhì)點索引，下同。OTM中的節(jié)點即為有限元單元中的節(jié)點，而物質(zhì)點則為有限元單元中的積分點，不同于有限元的是OTM中節(jié)點和物質(zhì)點不再有固定的連接關(guān)系，物質(zhì)點與節(jié)點的聯(lián)系將在每一步計算過程中動態(tài)構(gòu)建，即此“單元”沒有固定的形狀，其形狀將根據(jù)物質(zhì)的變形可以是任意的形狀，由于解放了網(wǎng)格的束縛，因此克服了拉氏方法在處理大變形問題時因網(wǎng)格畸變帶來的困難[36]。

圖2 空間離散示意圖[36]

材料的動力學(xué)信息存儲在節(jié)點上，包括位移、速度、加速度與溫度等。tk+1時刻的節(jié)點的位移、速度及當(dāng)前的溫度場更新公式為：

(7)

(8)

(9)

從tk時刻到tk+1時刻的位移傳輸(變形)映射φk→k+1為：

φk→k+1=∑a∈NH(xp,k)xa,k+1Na(xp,k)

(10)

材料的物理信息存儲在物質(zhì)點上，包括質(zhì)量、體積、密度、變形、應(yīng)力與材料內(nèi)部參數(shù)等。物質(zhì)點的運動通過節(jié)點的動力學(xué)信息插值獲得。在計算過程中，物質(zhì)點在tk+1的位置和溫度更新，通過對鄰域內(nèi)的節(jié)點進(jìn)行插值獲得：

xp,k+1=∑a∈NH(xp,k)xa,k+1Na(xp,k)

(11)

(12)

式中：Na(xp,k)為局部最大熵?zé)o網(wǎng)格插值函數(shù)(Local Maximum Entropy, LME)[37,38]，NH(xp,k)代表物質(zhì)點xp的鄰域。LME插值函數(shù)滿足嚴(yán)格的非負(fù)性，以及0階和1階連續(xù)性要求，同時在邊界滿足Kronecker-Delta屬性。此外，在LME插值函數(shù)中，通過調(diào)整γ值，插值函數(shù)作用范圍可從局部有限元無縫過渡到全局無網(wǎng)格，如圖 3所示，這種衰減特性建立了與高斯徑向基函數(shù)的聯(lián)系，其重要作用是只有少量的節(jié)點對待求函數(shù)有明顯的貢獻(xiàn)，從而大大降低計算成本。

圖3 局部最大熵插值函數(shù)[38]

通過將上述離散的參數(shù)代入半離散變分公式中，并實施駐點條件，可以獲得完全離散的應(yīng)力平衡方程和熱平衡方程：

(13)

(14)

2.4 基于物理的裂紋擴(kuò)展算法

材料在極限條件下的裂紋擴(kuò)展過程包含了多種能量傳播與耗散方式的結(jié)合與競爭，準(zhǔn)確地描述這些特征需要采用基于物理的裂紋擴(kuò)展算法。OTM方法中采用基于本征侵蝕的EigenFracture裂紋擴(kuò)展方法，在該方法中以能量的方式來描述整個材料的響應(yīng)過程，每個物質(zhì)點代表一小塊物質(zhì)，允許物質(zhì)點失效，物質(zhì)點失效代表該小塊物質(zhì)內(nèi)部產(chǎn)生裂紋并形成自由表面，而自由表面的產(chǎn)生將伴隨能量耗散，裂紋路徑則由各種能量耗散的競爭與耦合共同決定。在EigenFracture中物質(zhì)點失效的判斷準(zhǔn)則為：

(15)

圖4 黑色點為一組失效物質(zhì)點組成的裂紋，圓圈內(nèi)的紅點為位于裂紋尖端的物質(zhì)點的鄰域η

由于EigenFracture算法采用材料固有參數(shù)能量釋放Gc為失效判據(jù)，不需要對裂紋的方向、大小進(jìn)行顯式描述，具備簡便的三維幾何與拓?fù)浣Y(jié)構(gòu)的處理方式，通過平均化能量克服了由于網(wǎng)格分布形式而帶來的收斂性問題和裂紋路徑網(wǎng)格相關(guān)的問題，嚴(yán)格的數(shù)學(xué)證明收斂于Griffith解[39]。此外，由于EigenFracture方法從材料失效物理機(jī)制出發(fā)，考慮了材料內(nèi)部能量耗散的各種機(jī)制，包括塑性變形、裂紋擴(kuò)展和相變，使得精確預(yù)測脆性或者韌性材料在不同溫度以及載荷下的裂紋擴(kuò)展成為可能。

3 應(yīng)用場景與結(jié)果分析

本節(jié)主要探討OTM/HOTM極限力學(xué)仿真方法在小行星防御技術(shù)中的應(yīng)用方向。

3.1 動能撞擊

動能撞擊防御技術(shù)是指撞擊器以一定的速度和角度撞擊小行星，使其自旋狀態(tài)和軌道發(fā)生改變。撞擊器可以選用航天器、火箭甚至可操控的小行星。這一技術(shù)的關(guān)鍵在于掌握動能撞擊過程的動態(tài)響應(yīng)和動能撞擊的能量傳遞規(guī)律。動能撞擊防御方法技術(shù)簡單、啟動迅速、靈活性強(qiáng)、作用效果明顯，是一種實際可行的成熟技術(shù)。由于動能撞擊過程中涉及高溫、高壓、高應(yīng)變率等極限條件下材料的高度非線性動力學(xué)響應(yīng)，在發(fā)生超高速動能撞擊的局部區(qū)域，小行星表面材料在瞬間經(jīng)歷高溫、高壓、高應(yīng)變率極端狀態(tài)后，將發(fā)生大變形、斷裂、破碎及熔化、氣化乃至等離子體化等復(fù)雜的力學(xué)、物理過程及其耦合作用。溫度從300K增到104K，壓強(qiáng)從0.1MPa增至1TPa，應(yīng)變率高達(dá)107s-1。其中，熔化和氣化是超高速撞擊下材料動態(tài)響應(yīng)的重要特征，是超高速撞擊過程中能量轉(zhuǎn)化的主要表現(xiàn)形式[3]。對這些現(xiàn)象的描述直接影響著動能撞擊過程中小行星動態(tài)響應(yīng)和動能撞擊的能量傳遞規(guī)律的準(zhǔn)確性。基于OTM通過有機(jī)結(jié)合物質(zhì)點空間離散、最優(yōu)運輸理論時間離散、局部最大熵?zé)o網(wǎng)格近似，以及基于物理的裂紋擴(kuò)展算法和固液氣全域材料模型，為動能撞擊仿真提供了解決方案，如圖5所示。

圖5 基于OTM的超高速撞擊仿真模擬

在OTM超高速撞擊仿真中，采用J2粘塑性模型(J2-power law viscoplasticity)描述材料應(yīng)變硬化、應(yīng)變率硬化和高溫軟化等效應(yīng)。高溫高壓下材料變形與溫度、壓力的關(guān)系可采用SESAME狀態(tài)數(shù)據(jù)庫或Grüneisen狀態(tài)方程。通過采用Lindemann修正模型[40]描述熔化溫度的體積依賴性。目前該方案已成功應(yīng)用于金屬[41-43]、尼龍[44]、陶瓷[45]、冰[46]等材料沖擊碰撞仿真，如圖 6所示。圖 7為利用OTM方法初步模擬動能撞擊過程中的彈性波傳播、塑性變形、裂紋擴(kuò)展，以及由于塑性與狀態(tài)方程引起的溫度升高而導(dǎo)致的材料相變之間的相互耦合與競爭關(guān)系的過程。

圖6 OTM超高速撞擊仿真應(yīng)用

圖7 小行星動能撞擊仿真

為了詳細(xì)描述小行星內(nèi)部結(jié)構(gòu)在超高速撞擊效應(yīng)下的動態(tài)熱力學(xué)響應(yīng)，可對復(fù)雜的小行星體建立分層復(fù)合材料模型，基于OTM的多場耦合變分理論框架，將對小行星體撞擊過程中的彈性波傳播、塑性變形、裂紋擴(kuò)展以及由于塑性與狀態(tài)方程引起的溫度升高而導(dǎo)致的材料相變之間的相互耦合與競爭關(guān)系進(jìn)行精準(zhǔn)的描述，量化狀態(tài)方程、材料強(qiáng)度以及熱效應(yīng)對超高速撞擊下行星體動態(tài)熱力學(xué)響應(yīng)的影響。

3.2 激光燒蝕

激光燒蝕驅(qū)動技術(shù)是指使用強(qiáng)激光照射小行星表面，利用表面燒蝕產(chǎn)生的等離子體噴射所帶來的反作用力，使小行星自旋狀態(tài)和軌道發(fā)生改變。大量研究認(rèn)為，激光燒蝕驅(qū)動技術(shù)是一種高效的空間碎片清除技術(shù)[47-50]，基于同樣的作用機(jī)理可用于小行星防御。激光驅(qū)動改變小行星軌道是一種非接觸式防御方法，其關(guān)鍵問題在于燒蝕驅(qū)動小行星機(jī)理、對驅(qū)動效果的影響因素、驅(qū)動的動力學(xué)模型。

在HOTM方法中，熱流固耦合整體求解框架結(jié)合固液氣相敏全域材料模型為小行星的激光燒蝕仿真提供了理想解決方案，而施加在小行星表面的激光可通過熱通量邊界條件進(jìn)行模擬，頂部平坦的超高斯分布熱流可以較好地近似物理激光束[51]，該激光模型可以被描述為：

(16)

式中：A為材料的激光吸收率，P為激光功率，r為激光半徑，xc(t)、yc(t)為激光中心坐標(biāo)，可用任意隨時間的變化的函數(shù)描述。通過這種模型可容易地控制激光軌跡、速度等參數(shù)。圖 8(a)表示n取不同值時的超高斯熱流的分布，其中橫坐標(biāo)為該點到激光圓心距離與激光半徑之比，縱坐標(biāo)為取值。在HOTM方法中使用光線追蹤算法來動態(tài)獲取材料表面施加熱流條件的節(jié)點集。計算域內(nèi)其他部分的加熱則通過熱傳導(dǎo)來實現(xiàn)。除了激光熱通量之外，還可以施加熱對流邊界條件和熱輻射。

在HOTM方法中，小行星表面物質(zhì)氣化形成的等離子體對小行星形成的反沖壓力可由施加在表面結(jié)點的力邊界條件實現(xiàn)，具體為采用Anisimov[52]等提出的隨溫度變化的模型來模擬反沖壓力，表達(dá)式如下：

(17)

圖8 (a)高斯階數(shù)n取值不同時的高斯熱流分布；(b)金屬鋁的反沖壓力曲線

圖9 (a)金屬粉末床激光選區(qū)熔融工藝(SLM)仿真；(b)超高速激光金屬熔覆工藝(EHLA)仿真

激光燒蝕技術(shù)使用大功率激光照射到小行星表面使其表面物質(zhì)溫度急劇上升，熔化及氣化并產(chǎn)生等離子體。劇烈氣化的物質(zhì)向外噴射產(chǎn)生很大的反沖壓力推動小行星運動，使其速度發(fā)生變化并改變軌跡，如圖 10所示。相比于通過爆炸改變PHA軌跡的方法，激光實際部署地點與PHA一般相距較遠(yuǎn)，實際照射到PHA表面的激光能量密度較低，所以不存在PHA表面碎片飛出的風(fēng)險，因而采用這種技術(shù)更加安全。

圖10 小行星激光燒蝕驅(qū)動仿真

3.3 核爆攔截

核爆是應(yīng)對短預(yù)警時間、大尺寸小行星撞擊的技術(shù)。核爆炸防御方式有兩種：一是利用核爆裝置直接炸毀PHA；二是利用核爆直接炸毀小行星或改變小行星的軌道以避免其與地球相撞，是近地小行星防御最主要的手段之一。根據(jù)PHA尺寸、材質(zhì)、結(jié)構(gòu)的不同，可選擇表面爆炸、對峙爆炸以及穿透爆炸三種方式[4]。①表面爆炸：針對小體積的PHA，可以采用作用能量較大的表面爆炸或淺地下爆炸的方式，使PHA分裂成數(shù)塊碎片；②對峙爆炸：對峙爆炸是在距離PHA表面一定距離時引爆核裝置，爆炸產(chǎn)生的熱中子、X射線以及γ射線輻射PHA表面，產(chǎn)生高溫，引發(fā)PHA表面物質(zhì)的噴射，噴射時產(chǎn)生的推力使PHA發(fā)生偏轉(zhuǎn)。此外，爆炸產(chǎn)生的部分碎片與PHA發(fā)生撞擊，傳遞動能。兩種作用效果疊加，實現(xiàn)防御目的。對峙爆炸是規(guī)避爆炸碎片威脅的有效方法之一，適用于防御體積較大的PHA；③穿透爆炸：穿透核爆炸是指核裝置穿入PHA內(nèi)部一定深度處發(fā)生爆炸。該方法的優(yōu)勢在于，除了核爆產(chǎn)生的爆炸能量外，爆炸引起的表面沖擊波能夠擴(kuò)大作用威力，穿透深度很淺的爆炸就足以改變PHA的運行軌道。

利用OTM方法可從核彈侵徹開始模擬侵徹過程因高速沖擊壓力、塑性變形升溫、摩擦生熱等共同作用引起炮彈內(nèi)部炸藥點火起爆、爆轟波在小行星內(nèi)部傳播、爆炸產(chǎn)生的高能氣體膨脹流固耦合結(jié)構(gòu)破壞的全過程。其中，采用Lee-Tarver三項點火增長模型[54]描述炸藥點火起爆，在此模型中，點火、生長和完成分為三個過程。當(dāng)炸藥完成點火增長起爆后，爆炸產(chǎn)生的爆轟波氣體以JWL狀態(tài)方程[55]的形式膨脹對結(jié)構(gòu)發(fā)生破壞作用，JWL狀態(tài)方程由Jones、Wikins及Lee提出的一種不顯含化學(xué)反應(yīng)、由實驗方法確定參數(shù)的半經(jīng)驗狀態(tài)方程，能較精確地描述爆轟產(chǎn)物膨脹驅(qū)動做功過程。JWL狀態(tài)方程及其等熵方程由三項組成，第一項主要在高壓區(qū)起作用，第二項在中壓區(qū)起作用，第三項在低壓區(qū)起作用。爆炸產(chǎn)生的高溫氣體先沖擊結(jié)構(gòu)表面，然后沖擊波在結(jié)構(gòu)內(nèi)部傳播，可引發(fā)結(jié)構(gòu)的大范圍碎裂效應(yīng)。通過結(jié)合EigenFracture裂紋擴(kuò)展算法，可描述爆轟波在結(jié)構(gòu)內(nèi)部的流固耦合破壞作用。圖 11初步展示了核彈穿入PHA內(nèi)部發(fā)生爆炸引起的表面沖擊波和核爆產(chǎn)生的爆炸能量共同作用使PHA分裂成許多碎片的過程。

圖11 小行星核爆攔截仿真

4 總結(jié)與展望

針對小行星防御數(shù)值仿真在基礎(chǔ)理論、計算方法、材料模型與高性能計算方面存在的困難與挑戰(zhàn)，本文介紹了OTM/HOTM極限力學(xué)仿真理論。OTM/HOTM有機(jī)融合了多物理場變分原理、基于物理的本構(gòu)理論和基于能量的斷裂力學(xué)理論，為小行星防御中綜合考慮熱流固耦合、熱力化學(xué)耦合等多物理場效應(yīng)強(qiáng)耦合復(fù)雜的物理化學(xué)和力學(xué)現(xiàn)象提供了理論基礎(chǔ)；在數(shù)值離散方面，OTM采用了最優(yōu)運輸時間積分、物質(zhì)點空間離散、局部最大熵?zé)o網(wǎng)格近似技術(shù)，并采用多線程多進(jìn)行混合并行策略實現(xiàn)大規(guī)模并行計算，克服了基于網(wǎng)格的數(shù)值方法和無網(wǎng)格方法在精度、穩(wěn)定性、健壯性、效率等方面的不足；基于這些特性所發(fā)展的ESCAAS極限力學(xué)仿真平臺實現(xiàn)了在統(tǒng)一的框架下求解材料超大變形、熔化、氣化、沖擊爆炸、流固與熱流固耦合、自由表面、多體接觸、碎裂、層裂與碎片云等復(fù)雜物理現(xiàn)象的多物理場、多尺度強(qiáng)耦合問題；最后，展示了OTM/HOTM方法及ESCAAS平臺在隕石超高速撞擊成坑、動能撞擊、激光燒蝕驅(qū)動、核爆攔截等小行星防御措施仿真中的有效性。由于OTM理論框架具有良好的可擴(kuò)展性，后續(xù)可在小行星材料模型、動能撞擊噴出的巖石羽流、破碎的石塊和塵埃、激光燒蝕驅(qū)動小行星機(jī)理、對驅(qū)動效果的影響因素、驅(qū)動的動力學(xué)模型、高能量密度的核爆炸藥材料模型、核爆機(jī)理及對小行星的破壞效果等一系列融合新的物理化學(xué)現(xiàn)象、新的材料模型和更加真實的邊界條件方面進(jìn)一步開展研究工作，形成小行星防御專用仿真平臺，更加真實高效地仿真各種小行星防御場景。