廖祜明，黎 波，樊 江，焦立新，于帥超，林健宇，裴曉陽

（1. 北京航空航天大學(xué)能源與動力工程學(xué)院，北京 102206；2. 云翼超算(北京)軟件科技有限公司，北京 100027；3. 北京大學(xué)工學(xué)院，北京 100871；4. 中國工程物理研究院流體物理研究所，四川 綿陽 621900）

地球軌道上散布著大量的空間碎片和微流星體，若它們與在軌運行航天器相撞，相對撞擊速度可達(dá)十幾千米每秒，對航天器造成巨大威脅。大型碎片撞擊可直接撞毀航天器；中型碎片撞擊會造成航天器部件功能失效；而小型碎片撞擊可在航天器表面造成累積損傷，導(dǎo)致部件功能下降[1]。因此，空間碎片超高速撞擊防護是航天器結(jié)構(gòu)設(shè)計須重點考慮的問題。厘清空間碎片超高速撞擊的力學(xué)行為是防撞設(shè)計的重要基礎(chǔ)。碎片超高速撞擊的物理過程非常復(fù)雜，撞擊過程中激波在材料中傳播，同時激活各種相互耦合與競爭的能量耗散機制，引起材料的非彈性大變形、相變（液化、氣化、凝固等）、裂紋萌生及擴展和層裂、碎裂等復(fù)雜的物理過程。采用數(shù)值方法模擬超高速撞擊面臨如何解決速度過大導(dǎo)致的網(wǎng)格畸變，以及如何準(zhǔn)確描述超高速撞擊過程中各種能量耗散模式之間的競爭與耦合關(guān)系等巨大挑戰(zhàn)。

目前分析超高速撞擊問題的數(shù)值計算方法主要包括基于網(wǎng)格的方法和無網(wǎng)格方法兩大類?；诰W(wǎng)格的數(shù)值方法中，傳統(tǒng)拉氏有限元法受網(wǎng)格畸變困擾，難以處理超大變形問題；歐氏方法具有網(wǎng)格不變形的固有優(yōu)勢，但如何跟蹤歷史變量和物質(zhì)界面仍然是亟待解決的問題，同時在模擬動態(tài)裂紋擴展、多體碰撞與應(yīng)變率相關(guān)的材料熱力學(xué)行為時還需進一步研究；任意拉格朗日-歐拉（arbitrary Lagrangian-Eulerian, ALE）[2-3]法在一定程度上緩解了網(wǎng)格畸變帶來的困難，但ALE 繼承了拉氏有限元的大部分缺點，同樣面臨大變形時網(wǎng)格畸變的問題[4]。無網(wǎng)格法由于不需要進行網(wǎng)格離散及采用高階插值函數(shù)，有利于解決大變形問題，在超高速撞擊問題上應(yīng)用更為廣泛。常用的無網(wǎng)格方法有光滑粒子流體動力學(xué)法（smoothed particle hydrodynamics, SPH）[5]、再生核粒子法（reproducing kernel particle method, RKPM）[6]以及物質(zhì)點法（material point method, MPM）[7]等。有大量的研究應(yīng)用以上方法[8-9]及其改進方法[10]對超高速撞擊問題進行仿真分析并取得了豐富的研究成果，但在算法的精度、穩(wěn)定性、健壯性、效率等方面仍然存在許多不足[11]。

最優(yōu)運輸無網(wǎng)格方法（optimal transportation meshfree, OTM）[12]是一種專門針對動態(tài)沖擊問題提出的顯式增量更新拉格朗日無網(wǎng)格方法。OTM 方法采用局部最大熵?zé)o網(wǎng)格插值函數(shù)（local maximum entropy,LME）[13]，在邊界上滿足Kronecker-Delta 屬性，解決了傳統(tǒng)無網(wǎng)格法在處理位移邊界條件時的缺陷；同時在節(jié)點的基礎(chǔ)上增加物質(zhì)點，以物質(zhì)點作為積分點有效避免了計算結(jié)果在拉伸載荷下的不穩(wěn)定性；在時間離散上采用了最優(yōu)運輸理論（optimal transport theory）[14]，保證了其時間離散形式在理論上嚴(yán)格滿足動量守恒，具有嚴(yán)格的收斂性數(shù)學(xué)證明；在并行計算上結(jié)合了共享內(nèi)存多線程并行與分布式多進程并行實現(xiàn)了大規(guī)模多層混合并行[15]，顯著提高了計算效率；在材料失效方面，采用本征侵蝕裂紋擴展算法（Eigenerosion），以材料固有屬性能量釋放率作為失效判據(jù)[16-18]，克服了傳統(tǒng)裂紋擴展數(shù)值方法無法清晰表征材料真實變形和失效物理機制的缺點，理論證明收斂到Griffith 理論，與網(wǎng)格無關(guān)。OTM 方法的這些特性為超高速撞擊模擬預(yù)測提供了高效精準(zhǔn)的解決方案[19]。本文采用基于OTM 方法的極限力學(xué)仿真軟件ESCAAS，對空間碎片超高速撞擊問題進行數(shù)值模擬，并與實驗結(jié)果進行對比驗證。

1 OTM 最優(yōu)運輸無網(wǎng)格方法

1.1 時間與空間離散

假設(shè)參考構(gòu)型中的連續(xù)介質(zhì) ?0R3運動過程由變形映射（位移） φ:?0×[t0;t]!R3描述，其中 [t0;t] 為運動過程的總時間。X2 ?0表示材料質(zhì)點在參考構(gòu)形中的位置，x=φ(X;t) 表示材料質(zhì)點在現(xiàn)時構(gòu)形?t=φ(?0;t)的位置，則拉格朗日描述下質(zhì)量守恒、動量守恒與能量守恒方程為：

式中：下標(biāo)k代表時間步長索引。式(5)為時間離散的變分公式，為了便于計算機求解，要對此半離散的變分公式進行空間離散得到時間-空間全離散格式。

在有限元法中，每個單元有特定的形狀且各單元間具有固定的連接關(guān)系，每個單元攜帶了材料一定量的質(zhì)量并具有相應(yīng)的體積，材料的響應(yīng)在單元的積分點上進行求解。與有限元法類似，OTM 方法中采用物質(zhì)點xp與節(jié)點xa相結(jié)合的方式對半離散增量變分公式進行無網(wǎng)格空間離散，其中下標(biāo)a表示節(jié)點索引，p表示物質(zhì)點索引。如圖1 所示，所有的場變量信息分別由節(jié)點xa和物質(zhì)點xp攜帶。物質(zhì)點可以理解為有限元中的積分點，而其攜帶的體積可以看作是有限元中的單元，而不同于有限元的是節(jié)點和物質(zhì)點不再有固定的連接關(guān)系，即此“單元”沒有固定的形狀，其形狀將根據(jù)物質(zhì)的變形可以是任意的形狀，由于解放了網(wǎng)格的束縛，因此對于大變形問題不存在網(wǎng)格畸變帶來的問題[21]。

材料的動力學(xué)信息存儲在節(jié)點上，包括位移、速度、加速度與溫度等。在計算過程中tk+1時刻節(jié)點位置xa;k+1和節(jié)點速度 φ˙a;k+1的更新公式為：

從tk時刻到tk+1時刻的位移傳輸（變形）映射 φk!k+1為：

式中：Na為LME 局部最大熵插值函數(shù)，NH(xp;k代表物質(zhì)點xp在tk時刻下的鄰域。

材料的物理信息存儲在物質(zhì)點上，包括質(zhì)量、體積、密度、變形、應(yīng)力與材料內(nèi)部參數(shù)等。材料的響應(yīng)在物質(zhì)點上求解，物質(zhì)點攜帶的質(zhì)量在計算過程中保持不變，而體積和密度在不斷的變化，每個時間步都得到更新。物質(zhì)點的運動通過節(jié)點的動力學(xué)信息插值獲得[22]。在計算過程中，物質(zhì)點在tk+1的位置更新，通過對鄰域內(nèi)的節(jié)點進行插值獲得：

物質(zhì)點鄰域范圍NHxp;k在每個時間步進行更新，同時更新各節(jié)點上插值函數(shù)值。LME 插值函數(shù)Na滿足嚴(yán)格的非負(fù)性，以及零階和一階連續(xù)性要求，同時在邊界滿足Kronecker-Delta 屬性。此外，在LME 插值函數(shù)中，通過調(diào)整 γ 值，插值函數(shù)作用范圍可從局部有限元無縫過渡到全局無網(wǎng)格，如圖2 所示，這種衰減特性建立了與高斯徑向基函數(shù)的聯(lián)系，其重要作用是只有少量的節(jié)點對待求函數(shù)有明顯的貢獻(xiàn)，大大降低計算成本。

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

1.2 捕捉式多體動態(tài)接觸算法

OTM 方法采用物質(zhì)點和節(jié)點進行離散，并且動態(tài)更新局部鄰域，使其具備自動捕獲接觸（seizing contact）的能力。具體地，物體 ?A上的節(jié)點同時也可以作為物體 ?B中物質(zhì)點的鄰域，因此對兩物體的接觸可以不加特殊處理，如圖3 所示，其中Ra為節(jié)點xa的插值函數(shù)的作用范圍，Rb同理，vp;k代表物質(zhì)點速度，vb;k代表節(jié)點速度。

在捕捉式多體動態(tài)接觸算法中，當(dāng)碰撞節(jié)點進入物質(zhì)點鄰域后，其線性動量會被自然地抵消，如碰撞節(jié)點xa的線性動量la;k為：

式中：M為節(jié)點鄰域內(nèi)的物質(zhì)點數(shù)量，為mp為物質(zhì)點質(zhì)量。節(jié)點的動量為其鄰域物質(zhì)點動量的平均值，結(jié)合運算中采用的集中質(zhì)量法，節(jié)點的速度本身代表了局部物質(zhì)點速度的質(zhì)量加權(quán)平均值，平均了速度的所有分量，接觸時物體之間的法向速度和切向速度都被用于相互抵消，從而自動確定了碰撞體之間的接觸條件。因此，這種捕捉式接觸時的摩擦系數(shù)可視為無限大，適用于模擬碰撞終端拋射物與目標(biāo)物混合及伴隨的絕熱加熱過程。

1.3 本征侵蝕裂紋擴展算法

材料在極限條件下的裂紋擴展過程包含了多種能量傳播與耗散方式的結(jié)合與競爭，準(zhǔn)確地描述這些特征需要采用基于物理的裂紋擴展算法。OTM 方法中采用能量的方式來描述整個材料的響應(yīng)過程，提出了本征侵蝕（Eigenerosion）裂紋擴展算法。在該方法中，每個物質(zhì)點代表一小塊物質(zhì)，通過移除失效物質(zhì)點的方式來模擬材料中裂紋擴展和碎裂的形成，物質(zhì)點失效代表該小塊物質(zhì)內(nèi)部產(chǎn)生裂紋并形成自由表面，而自由表面的產(chǎn)生則必然伴隨著能量耗散，裂紋路徑是由彈性能量（彈性能量作為能量釋放機制促進斷裂）、材料特定的斷裂能量（與裂紋面積成正比的斷裂比能）、單調(diào)性和接觸約束之間的競爭引起的結(jié)果。物質(zhì)點失效的判斷準(zhǔn)則為：

式中：Gp表示物質(zhì)點xp的等效能量釋放率，下標(biāo)p代表計算域內(nèi)物質(zhì)點的索引；Cη表示局部鄰域， η 決定了其范圍，kCηk 表示局部鄰域總體積；λ為幾何校正因子；Vq為局部鄰域內(nèi)物質(zhì)點xq體積，下標(biāo)q代表鄰域內(nèi)物質(zhì)點的索引。如圖4 所示，黑色點為一組失效物質(zhì)點組成的裂紋，圓圈內(nèi)的紅點為位于裂紋尖端的物質(zhì)點的鄰域 η ，通過平均該物質(zhì)點鄰域內(nèi)的彈性能We得到物質(zhì)點的等效能量釋放率Gp，若Gp大于材料本身的臨界能量釋放率Gc時 ，該物質(zhì)點將失效不再參與計算。

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

2 銅飛片超高速撞擊鋁合金靶板數(shù)的值模擬與驗證

2.1 數(shù)值模擬模型

文獻(xiàn)[25]開展了柱形OFHC 銅飛片超高速撞擊Al6061-T6 鋁合金靶板實驗，并運用高速攝像系統(tǒng)拍攝了碎片云形貌。本節(jié)將采用基于OTM 方法開發(fā)的ESCAAS 軟件，建立同樣的超高速撞擊模型（如圖5 所示），對該超高速撞擊過程進行模擬并開展與實驗結(jié)果的對比驗證。

工況1：撞擊角度為5.4°，質(zhì)量為3 g 的柱形OFHC 銅飛片，以5.55 km/s 的速度撞擊Al6061-T6 鋁合金靶板，幾何模型如圖5(a)所示，靶板外側(cè)施加固定邊界條件，取銅飛片和靶板的1/2 對稱模型進行無網(wǎng)格離散，物質(zhì)點最小特征尺寸為0.14 mm，共計31.3 萬個物質(zhì)點，6.26 萬個節(jié)點；

工況2：撞擊角度為11.7°，質(zhì)量為10 g 的柱形OFHC 銅飛片，以5.12 km/s 的速度撞擊Al6061-T6 鋁合金靶板，幾何模型如圖5(b)所示，靶板外側(cè)施加固定邊界條件；取銅飛片和靶板的1/2 對稱模型進行無網(wǎng)格離散，物質(zhì)點最小特征尺寸為0.6 mm，共計40.6 萬個物質(zhì)點，22.9 萬個節(jié)點，如圖5(c)所示。

2.2 材料模型

超高速撞擊過程涉及材料應(yīng)變硬化、應(yīng)變率硬化和高溫軟化等效應(yīng)，J2 黏塑性模型（J2-power law viscoplasticity）能較好地模擬銅和鋁合金材料塑性響應(yīng)，其等效屈服應(yīng)力為：

表1 材料物性參數(shù)Table 1 Material parameters

表2 材料本構(gòu)模型參數(shù)Table 2 Parameters of material constitutive model

2.3 結(jié)果分析

工況1：在Linux 平臺下采用ESCAAS 模擬 3 g 柱形銅飛片，以5.4°的撞擊角度、5.55 km/s 的速度撞擊鋁合金靶板的過程。整個模擬采用16 核并行計算，10 μs 撞擊過程模擬用時約為1.33 h。圖6 為銅飛片穿透緩沖墻形成碎片的過程，模擬結(jié)果很好地呈現(xiàn)出碎片云的形態(tài)，包括碎片云呈現(xiàn)出頭部橢圓形、內(nèi)核碎片云的位置、外泡碎片云的形態(tài)以及反濺碎片云的形態(tài)等特征信息。

除了碎片云輪廓形貌的定性對比外，碎片云特征尺寸定量對比方面：鋁合金靶板前端碎片云輪廓軸向最大位置點Point 1 距靶板的軸向距離的OTM 模擬結(jié)果為(50.74±0.66) mm，實驗測量結(jié)果為51.3 mm，與實驗的相對誤差為(1.09±1.29)%；銅飛片前端碎片云輪廓軸向最大位置點Point 2 距靶板的軸向距離的OTM 模擬結(jié)果為(31.9±1.4) mm，實驗測量結(jié)果為33.5 mm，與實驗的相對誤差為(4.8±4.18)%；鋁合金靶板反濺碎片云的高度測量點Point 6 距靶板的軸向距離的OTM 模擬結(jié)果為(7.87±0.74) mm，實驗測量結(jié)果為7.76 mm，與實驗對比相對誤差為(1.42±9.54)%；銅飛片和靶板材料碎片徑向相交距離（測量點Point 4 和Point 5 的徑向位置）的OTM 模擬結(jié)果為(27.6±1.3) mm，實驗測量結(jié)果為27.9 mm，與實驗對比相對誤差為(1.07±4.66)%。

在超高速撞擊過程中，彈丸厚度/板厚(D/H)、打擊速度與角度的結(jié)合決定了沖擊波的分段分區(qū)效應(yīng)。本例中D/H=1.2，屬于中厚板情形，其主要失效模式為塑性變形、裂紋擴展、熔化三種能量轉(zhuǎn)換與耗散方式相結(jié)合。如圖8 所示，OTM 方法模擬了飛片撞擊下，靶板發(fā)生塑性變形引起靶板翻邊撕裂的過程。測量得到OTM 模擬的穿孔直徑為(20.9±0.4) mm，實驗測試值約為20.1 mm，相對誤差為(3.98±1.99)%。從上述分析可見仿真量化結(jié)果與實驗吻合良好。

工況2：模擬了10 g 柱形OFHC 銅飛片，以11.7°的撞擊角度，5.12 km/s 的速度撞擊Al6061-T6 鋁合金靶板的過程。整個模擬采用16 核并行計算，10 μs 撞擊過程模擬時間約為2.4 h。圖9 所示為不同時刻下碎片云的輪廓，由于銅飛片撞擊角度11.7°較大，銅飛片碎片云與靶板碎片云前端呈現(xiàn)出顯著偏轉(zhuǎn)。

圖10 為7.6 μs 時刻下，OTM 模擬與實驗測量的碎片云形貌對比情況。同樣地，軸向位置點Point 1 為靶板碎片云前端；軸向位置點Point 2 為銅飛片碎片云前端；軸向點Point 3 為撞擊之前靶板的前端；徑向位置點Point 4 和點Point 5 為銅飛片和靶板碎片徑向相交距離測量點；徑向位置點Point 6 為靶板后端反濺碎片云高度。此外，由于撞擊角度11.7°較大，銅飛片碎片云與靶板碎片云前端呈現(xiàn)出顯著偏轉(zhuǎn)，增加了碎片云偏轉(zhuǎn)角度與實驗測量結(jié)果的對比情況，包括碎片云前端偏離中線的角度 α ，靶板上側(cè)反濺碎片云偏離角 θ ，靶板下側(cè)反濺碎片云偏離角 β 。

銅飛片速度為5.12 km/s 時，鋁合金靶板前端碎片云輪廓軸向最大位置點Point 1 距靶板軸向距離的OTM 模擬結(jié)果為(47.2±1.2) mm，實驗測量結(jié)果為47.5 mm，與實驗的相對誤差為(0.63±2.53)%；銅飛片前端碎片云輪廓軸向最大位置點Point 2 距靶板軸向距離的OTM 模擬結(jié)果為(27.1±2.1) mm，實驗測量結(jié)果為26.0 mm，與實驗的相對誤差為(4.23±8.08)%；鋁合金靶板反濺碎片云的高度測量點Point 6 距靶板軸向距離的OTM 模擬結(jié)果為(32.0±2.6) mm，實驗測量結(jié)果為30.9 mm，與實驗對比相對誤差為(3.56±8.41)%；碎片云前端偏離中線的角度 α的OTM 模擬結(jié)果為(10.55±0.65)°，實驗測量結(jié)果為11.5°，與實驗對比相對誤差約為(8.26±5.65)%；靶板上側(cè)反濺碎片云偏離角 θ 的OTM 模擬結(jié)果為(109.5±3.5)°，實驗測量結(jié)果為100.8°，與實驗對比相對誤差約為(8.63±3.47)%；靶板下側(cè)反濺碎片云偏離角 β 的OTM 模擬結(jié)果為(109.75±3.25)°，實驗測量結(jié)果為108.0°，與實驗對比相對誤差約為(1.62±3.00)%，模擬結(jié)果與實驗結(jié)果吻合理想。

圖11 所示為工況2 銅飛片撞擊鋁合金靶板，本例中D/H=1.17 屬于中厚板，同樣地，在撞擊過程中可以觀察到鋁合金靶板具有顯著的塑性變形和翻邊效果。OTM 模擬的靶板穿孔直徑為32.125 mm，實驗測試值約為31.6 mm，相對誤差大約1.66%，穿孔孔徑與實驗高度吻合。

3 結(jié) 論

超高速撞擊涉及材料局部超大變形、裂紋擴展、層裂與碎裂、多體動態(tài)接觸、固液氣相變與多相混合等復(fù)雜的材料動態(tài)響應(yīng)過程，對傳統(tǒng)基于網(wǎng)格的數(shù)值方法提出了巨大挑戰(zhàn)。本文引入最優(yōu)運輸無網(wǎng)格方法（OTM）對空間碎片超高速撞擊過程進行數(shù)值分析，指出了OTM 在處理此類問題中有如下優(yōu)勢：

(1) OTM 方法基于變分原理，在統(tǒng)一的拉格朗日框架下求解物質(zhì)運動方程，實現(xiàn)了對不同材料在高溫、高壓、高應(yīng)變率等極限條件下多種能量耗散模式結(jié)合與競爭的現(xiàn)象進行了完善的考慮，在理論上保證了不同形式能量耗散的自主耦合分配；

(2) OTM 方法采用物質(zhì)點和節(jié)點空間離散，在時間與空間上的離散被嚴(yán)格地限制在變分框架內(nèi)，保證了質(zhì)量，能量及動量的守恒；結(jié)合局部最大熵?zé)o網(wǎng)格插值函數(shù)，克服了一般無網(wǎng)格法在精準(zhǔn)施加位移邊界條件上的不足，穩(wěn)定性和收斂性有完備的數(shù)學(xué)證明；

(3) OTM 方法采用了真實的裂紋擴展算法，考慮了材料內(nèi)部能量耗散的各種機制，包括塑性變形、裂紋擴展和相變，采用材料的固有屬性能量釋放率作為失效判據(jù)，而非簡單采用單元變形量、應(yīng)變大小，或是基于實驗的失效模型來作為材料的失效判據(jù)，克服了傳統(tǒng)裂紋擴展數(shù)值方法無法清晰表征材料真實變形和失效物理機制，以及不收斂和網(wǎng)格相關(guān)的缺點；

(4) OTM 方法及ESCAAS 軟件在計算彈孔直徑、碎片云形態(tài)方面與實驗結(jié)果吻合理想，驗證了其在超高速撞擊數(shù)值模擬方面的適用性，顯示出OTM 方法及ESCAAS 軟件可以作為超高速撞擊的有力數(shù)值分析手段。