馬天寶，粟鑫，郝莉

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081; 2.北京建筑大學(xué) 理學(xué)院，北京 100044)

在航天技術(shù)高速發(fā)展的今天，在軌航天器正時刻受到外層空間碎片的撞擊威脅.這些空間碎片分布在距地球表面2 000 km的低地軌道上，擁有極大的數(shù)量和極高的運(yùn)行速度，與航天器發(fā)生碰撞時其撞擊速度可達(dá)10 km/s以上，因此現(xiàn)今大多數(shù)的航天器多數(shù)采用在航天器艙壁外加裝Whipple屏[1]的防護(hù)形式來防護(hù)空間碎片.在對空間碎片防護(hù)技術(shù)的研究中，由于碰撞速度很高，受限于目前實(shí)驗(yàn)器材發(fā)展水平，實(shí)驗(yàn)研究存在很大的困難，因而數(shù)值模擬成為了研究航天器碎片防護(hù)技術(shù)中超高速碰撞問題的重要方法[2-3].

目前，超高速碰撞問題的數(shù)值模擬多采用光滑粒子流體動力學(xué)方法(smoothed particle hydrodynamics，SPH)及其相關(guān)的衍生方法進(jìn)行模擬.但SPH方法本身也存在一些問題，如存在拉力不穩(wěn)定性，施加位移邊界條件存在困難等.最優(yōu)輸運(yùn)無網(wǎng)格方法是一種顯式增量更新拉格朗日無網(wǎng)格計(jì)算方法，由M.Ortiz等[4]于2010年提出.該方法的理論體系由最優(yōu)輸運(yùn)理論、物質(zhì)點(diǎn)采樣法以及局部最大熵?zé)o網(wǎng)格近似[5]結(jié)合而成，與SPH方法相比，OTM方法具有以下優(yōu)勢：① 最優(yōu)輸運(yùn)理論為OTM方法中的顯式時間積分算法提供了可靠的理論基礎(chǔ)，保證了時間離散下系統(tǒng)的動量守恒與辛守恒，同時質(zhì)量守恒自動滿足無需額外計(jì)算，可以有效地克服現(xiàn)有顯式分析中能量不守恒的問題；② OTM在計(jì)算邊界處具有Kronecker-delta性質(zhì)，可以更好地引入位移邊界條件；③ OTM方法采用物質(zhì)點(diǎn)與節(jié)點(diǎn)結(jié)合的方式進(jìn)行空間離散，由于插值與積分在不同的離散點(diǎn)進(jìn)行，完全避免了SPH方法中存在的拉力不穩(wěn)定性.本文采用OTM方法模擬彈丸超高速碰撞單層及雙層薄靶板問題，通過與實(shí)驗(yàn)結(jié)果的對比，驗(yàn)證了OTM方法模擬超高速碰撞問題的有效性.

1 OTM方法

1.1 OTM方法理論

在參考構(gòu)型中，連續(xù)介質(zhì)運(yùn)動過程由變形映射φ:Ω0×[t0,t]→3進(jìn)行描述.假設(shè)材料的彈性響應(yīng)、塑性響應(yīng)以及比熱容是相互獨(dú)立的，因此系統(tǒng)的內(nèi)能密度為

W=We,vol+We,dev+Wp+Wh,

(1)

式中：We,vol為彈性響應(yīng)的體積部分;We,dev為彈性響應(yīng)的偏離項(xiàng);Wp為塑性響應(yīng);Wh為熱儲量.得到系統(tǒng)的內(nèi)能密度之后，可以推出以下變分結(jié)構(gòu),

(2)

式中：K為動能密度；B為單位質(zhì)量體積力；T為施加于Neumann邊界Γt上的牽引力.該變分結(jié)構(gòu)描述了材料系統(tǒng)的能量變化率.從變分原理出發(fā),對等效勢能Φ的變分取極值，再對該變分結(jié)構(gòu)進(jìn)行時間離散.在給定的時間增量[tn,tn+1]上，若tn時刻的狀態(tài)變量φn已知，則tn+1時刻下的狀態(tài)變量φn+1可以通過最小化增量能量近似獲得，即

(3)

根據(jù)最優(yōu)運(yùn)輸理論，時間間隔內(nèi)增量動能的精確最小值在[tn,tn+1]上可以表示為初始質(zhì)量密度和最終質(zhì)量密度之間的Wasserstein距離[6]，即

(4)

變分結(jié)構(gòu)(3)中除慣性項(xiàng)以外的其他項(xiàng)則采用后向差分格式進(jìn)行近似，可以得到以下半離散化(時間離散)的增量變分結(jié)構(gòu)，

(5)

對增量變分結(jié)構(gòu)的時間離散格式(5)求極值，將在時間段[tn,tn+1]內(nèi)產(chǎn)生動量守恒和能量守恒，使得系統(tǒng)的整個變形過程化為一系列離散的運(yùn)動過程.

在OTM方法中采用物質(zhì)點(diǎn)xp與節(jié)點(diǎn)xa相結(jié)合的方式對半離散增量變分公式進(jìn)行無網(wǎng)格空間離散.計(jì)算過程中，材料的物理信息存儲在物質(zhì)點(diǎn)上,而計(jì)算鄰域的動力學(xué)信息存儲在節(jié)點(diǎn)上，如圖1所示.

計(jì)算過程中，物質(zhì)點(diǎn)的運(yùn)動通過節(jié)點(diǎn)的動力學(xué)信息插值獲得，采用標(biāo)準(zhǔn)的Ritz-Galerkin法來近似位移場.假設(shè)N表示離散域內(nèi)的總節(jié)點(diǎn)數(shù)，則位移傳輸映射表示為

(6)

式中：Na(x)為節(jié)點(diǎn)xa的插值函數(shù)，插值函數(shù)Na(x)將采用局部最大熵?zé)o網(wǎng)格插值函數(shù).LME插值函數(shù)支持寬度自適應(yīng)改變，其導(dǎo)數(shù)滿足連續(xù)性要求，在計(jì)算邊界處具有Kronecker-delta性質(zhì)，LME的具體理論可詳見文獻(xiàn)[5].同樣地，基于Ritz-Galerkin法，物質(zhì)點(diǎn)在tn+1時刻的位置通過節(jié)點(diǎn)近似獲得，

(7)

式中：NH(xp,n)代表物質(zhì)點(diǎn)xp在tn時刻下的鄰域，其在每個時間步中將得到更新.

通過上述無網(wǎng)格近似方案對位移場進(jìn)行近似，增量變分勢能δΦn被完全地離散化，根據(jù)變分的極值條件，可得到系統(tǒng)的全離散格式力平衡方程，即

(8)

式中Ma,n+1代表節(jié)點(diǎn)xa在tn+1時刻的集中質(zhì)量，即

(9)

式中mp為物質(zhì)點(diǎn)xp所攜帶的質(zhì)量，節(jié)點(diǎn)xa在tn+1時刻的加速度則通過中心差分的形式近似獲得，因此節(jié)點(diǎn)坐標(biāo)的更新形式為

xa,n+1=xa,n+(tn+1-tn)×

(10)

式中l(wèi)a,n為節(jié)點(diǎn)xa在tn+1時刻的動量，即

(11)

節(jié)點(diǎn)xa在tn+1時刻的節(jié)點(diǎn)內(nèi)應(yīng)力和節(jié)點(diǎn)外力分別為

(12)

(13)

通過求解每個節(jié)點(diǎn)上的力平衡方程(8)，可以獲得整個系統(tǒng)的變形場.

1.2 數(shù)值迭代計(jì)算流程

基于OTM方法所編制程序的主要算法流程如圖2所示.

① 幾何模型的初始化.一般地，將d維的問題域(幾何模型)Ω0進(jìn)行網(wǎng)格劃分生成初始節(jié)點(diǎn)集{xa,0,a=1,2,…,N}，并取每個單元的形心作為一個物質(zhì)點(diǎn)的初始位置，組成物質(zhì)點(diǎn)集{xp,0,p=1,2,…,M}，物質(zhì)點(diǎn)的初始體積可設(shè)為該物質(zhì)點(diǎn)所在單元占據(jù)的空間大小.

② 物質(zhì)點(diǎn)和節(jié)點(diǎn)的初始化.在初始時刻n=0時，利用第1步離散得到的節(jié)點(diǎn)集合物質(zhì)點(diǎn)集，初始化節(jié)點(diǎn)坐標(biāo)xa,n、節(jié)點(diǎn)插值函數(shù)Na(xp,n)及其導(dǎo)數(shù)、節(jié)點(diǎn)力fa,n、節(jié)點(diǎn)動量la,n、節(jié)點(diǎn)質(zhì)量矩陣Ma,n+1，初始化物質(zhì)點(diǎn)坐標(biāo)xp,n、物質(zhì)點(diǎn)鄰域NH(xp,n)、物質(zhì)點(diǎn)體積vp,n、物質(zhì)點(diǎn)密度ρp,n以及物質(zhì)點(diǎn)質(zhì)量mp=ρp,nvp,n.

③ 節(jié)點(diǎn)坐標(biāo)更新.進(jìn)行tn→tn+1的瞬態(tài)分析，求解全離散格式的力平衡方程(8)，或節(jié)點(diǎn)坐標(biāo)更新公式(10)，顯示計(jì)算得到節(jié)點(diǎn)在tn+1時刻的坐標(biāo){xa,n+1,a=1,2,…,N}.

④ 物質(zhì)點(diǎn)的變形更新.根據(jù)節(jié)點(diǎn)坐標(biāo)，此時節(jié)點(diǎn)就已產(chǎn)生一個位移，同時物質(zhì)點(diǎn)將產(chǎn)生變形.采用Jaumann應(yīng)力率模型更新物質(zhì)點(diǎn)的應(yīng)力偏張量.

⑤ 根據(jù)第3步得到的節(jié)點(diǎn)坐標(biāo)，由式(7)更新物質(zhì)點(diǎn)的坐標(biāo)xp,n+1，更新物質(zhì)點(diǎn)的體積vp,n+1=det(Jp,n→n+1)vp,n與密度ρp,n+1=mp/vp,n+1，其中mp為物質(zhì)點(diǎn)xp的質(zhì)量，在計(jì)算過程中保持不變，Jp,n→n+1為tn→tn+1時物質(zhì)點(diǎn)xp的Jacobian矩陣.

⑥ 結(jié)合前兩步得到的數(shù)據(jù)，通過本構(gòu)關(guān)系獲得物質(zhì)點(diǎn)的應(yīng)力σp,n+1，其中應(yīng)力張量的偏張量部分由Jaumann應(yīng)力率模型給出，球張量部分由物質(zhì)點(diǎn)的密度ρp,n+1通過狀態(tài)方程給出；更新tn+1時刻由于物質(zhì)點(diǎn)變形引起的節(jié)點(diǎn)力fa,n+1，更新節(jié)點(diǎn)動量la,n及節(jié)點(diǎn)質(zhì)量矩陣Ma,n+1.

⑦ 鄰域和插值函數(shù)的更新.根據(jù)第5步得到的物質(zhì)點(diǎn)坐標(biāo)xp,n+1、體積vp,n+1與密度ρp,n+1，更新tn+1時刻下的鄰域NH(xp,n+1)；同時更新節(jié)點(diǎn)的插值函數(shù)Na(xp,n+1)及其導(dǎo)數(shù).

⑧ 在完成第7步后，即代表完成了物質(zhì)點(diǎn)和節(jié)點(diǎn)在一個時間步內(nèi)的力學(xué)動態(tài)分析，判斷當(dāng)前時間步n，如果已計(jì)算至最后一個時間步則結(jié)束計(jì)算，否則轉(zhuǎn)入第3步.

2 彈丸超高速碰撞薄靶問題的數(shù)值模擬

2.1 彈丸和靶板的本構(gòu)模型及狀態(tài)方程

為描述超高速碰撞中高溫高壓高應(yīng)變率對靶板材料屈服強(qiáng)度的影響，采用Johnson-Cook損傷模型[7]，即

σY=(1-D)(A+Bεn)×

(14)

D=∑(Δε/εf)，

(15)

損傷變量D將在每一個時間步上得到累積，其中εf是斷裂塑性應(yīng)變，與材料的應(yīng)力三軸度、應(yīng)變率和溫度有關(guān)，其描述如下，

εf=[D1+D2exp(D3σ*)]×

(16)

式中：D1～D5為材料參數(shù);σ*為應(yīng)力三軸度;σm為平均正應(yīng)力;σeq為Mises等效應(yīng)力.

超高速碰撞問題中常用的狀態(tài)方程為Mie-Grüneisen狀態(tài)方程，表達(dá)式為

p=pH+Γ0ρ0(e-eH)，

(17)

式中：p和e分別為靜水壓力和比內(nèi)能;ρ0為材料的初始密度；Γ0為材料的Grüneisen系數(shù).Mie-Grüneisen狀態(tài)方程的shock形式定義如下，

(18)

式中：μ=ρ/ρ0-1，pH和eH分別為Hugoniot曲線上靜水壓力和比內(nèi)能的參考值；S為沖擊波速度U與波后質(zhì)點(diǎn)速度up之間線性關(guān)系的斜率，通常取U=C0+Sup；C0為體積聲速.

2.2 彈丸超高速碰撞單層靶板

針對彈丸超高速撞擊單層薄靶的過程，采用OTM方法進(jìn)行模擬，數(shù)值工況參照Hiermaier等[8]的實(shí)驗(yàn)過程，彈丸和靶板皆采用LY12硬鋁合金材料，材料模型采用Johnson-Cook損傷模型和Mie-Grüneisen狀態(tài)方程，具體材料參數(shù)見表1、表2.

表1 本構(gòu)模型參數(shù)

表2 狀態(tài)方程參數(shù)

如圖3所示，鋁制球形彈丸半徑為5 mm，鋁靶板的尺寸為50 mm×50 mm×4 mm，鋁球的撞擊速度設(shè)為與實(shí)驗(yàn)[8]一致的6.18 km/s.計(jì)算過程中，時間步長由CFL條件給出，每個時間步長的量級約為0.1 μs，室溫設(shè)為300 K.

圖4為20 μs時OTM方法計(jì)算得到的碎片分布，將其與Hiermaier的實(shí)驗(yàn)[8]及Hiermaier數(shù)值模擬結(jié)果[8]進(jìn)行對比，可以發(fā)現(xiàn)OTM得到的碎片云，整個碎片云的大部分質(zhì)量中位于外泡碎片云的前端，氣泡的大致形貌較為圓潤，在形貌上和分布上都與實(shí)驗(yàn)結(jié)果吻合較好.

表3給出了OTM方法計(jì)算得到的碎片云形貌的具體特征參數(shù)，并與傳統(tǒng)SPH方法[8]、ASPH方法[9]的模擬結(jié)果進(jìn)行對比，其中：d1和d2分別為不包含和包含彈坑邊緣粒子的彈坑直徑；l為碎片云的膨脹長度；w為碎片云寬度；Δ為l/w與實(shí)驗(yàn)結(jié)果的相對誤差.從表中可以看出，OTM方法得到的碎片云形貌與實(shí)驗(yàn)結(jié)果更加吻合，證明OTM方法非常適合于模擬超高速碰撞問題.

表3 不同數(shù)值方法的碎片云模擬結(jié)果對比

Tab.3 Comparison of different algorithms’ simulation result for debris cloud

數(shù)值方法D1D2l/mmw/mml/wΔ/%OTM方法模擬結(jié)果29.635.6102.478.21.315.8Hiermaier實(shí)驗(yàn)結(jié)果[2]27.534.51.39Hiermaier模擬結(jié)果[2]35.01.1120.1Liu模擬結(jié)果[7]28.9105.186.11.2212.2

2.3 彈丸超高速碰撞雙層靶板

針對鋁彈丸超高速碰撞雙層鋁板問題，采用OTM方法進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)結(jié)果[10]相對比.實(shí)驗(yàn)中采用了9 mm直徑的鋁球以5 km/s的速度正碰撞雙層鋁靶板，第1層靶板厚1 mm，第2層靶板厚2 mm，兩板的尺寸都設(shè)為50 mm×50 mm，兩板相距17 mm.彈丸和兩層靶板皆采用LY12硬鋁合金材料，材料模型與材料參數(shù)的選取與2.2節(jié)相同.

圖5給出了彈丸超高速碰撞雙層靶板的計(jì)算結(jié)果，并與實(shí)驗(yàn)結(jié)果相比較.從圖5中的結(jié)果對比可以看出：在彈丸撞擊第1層板4 μs后，實(shí)驗(yàn)中的一次碎片云與后板開始接觸，這一點(diǎn)與模擬結(jié)果相同；在一次碎片云撞擊后板6.8 μs后，模擬得到的二次碎片云長度為15.8 mm，而實(shí)驗(yàn)結(jié)果顯示為15.4 mm，相差2.6%，表明模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好.

3 結(jié) 論

本文采用最優(yōu)輸運(yùn)無網(wǎng)格方法(OTM)，對鋁制球形彈丸超高速撞擊鋁薄板形成碎片云的過程進(jìn)行了數(shù)值模擬，得到以下結(jié)論.

① 在超高速撞擊單層靶板的算例中，OTM方法計(jì)算得到的碎片云形態(tài)及內(nèi)核碎片云的形狀和分布均與實(shí)驗(yàn)吻合較好.整個碎片云的大部分質(zhì)量中位于外泡碎片云的前端，氣泡的大致形貌較為圓潤，碎片云外圍粒子分布較為平均.超高速撞擊的彈坑直徑、碎片云膨脹距離和寬度，與傳統(tǒng)SPH方法、ASPH方法的模擬結(jié)果進(jìn)行比較，OTM方法在碎片云的形貌上與實(shí)驗(yàn)結(jié)果更加接近，說明OTM方法適合模擬超高速碰撞等動態(tài)大變形問題.

② 針對球形彈丸超高速碰撞雙層靶板的情況，采用OTM方法進(jìn)行模擬，給出了幾個關(guān)鍵時刻超高速碰撞過程的模擬結(jié)果.模擬結(jié)果顯示，一次、二次碎片云的形貌與實(shí)驗(yàn)較為吻合，模擬給出的二次碎片云長度僅與實(shí)驗(yàn)結(jié)果相差2.6%，模擬結(jié)果良好.