基于pOTM方法的超高速碰撞并行數(shù)值模擬研究

張?zhí)忑?，馬天寶，郝 莉

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

自20世紀(jì)60年代以來，隨著航天事業(yè)的快速發(fā)展，太空碎片也隨之增多。而這些大部分位于近地軌道的太空碎片相對(duì)于航天器的速度可達(dá)10 km/s以上，其極大的動(dòng)能可以對(duì)在軌航天器造成非常嚴(yán)重的威脅。由于太空碎片體積小、數(shù)量大，大部分航天器采用加裝Whipple[1-2]防護(hù)屏的方式來防護(hù)太空碎片。對(duì)于超高速碰撞問題，用理論分析方法建立一個(gè)普適的動(dòng)力學(xué)理論模型需要涉及到多種領(lǐng)域，目前很難實(shí)現(xiàn)；現(xiàn)有的實(shí)驗(yàn)技術(shù)也很難達(dá)到所需要的環(huán)境；因此采用數(shù)值模擬方法研究超高速碰撞問題成為當(dāng)下一種重要的方法。

對(duì)于超高速碰撞問題的數(shù)值模擬研究，無網(wǎng)格法因其自身的優(yōu)勢(shì)成為解決問題的主要方法。由Lucy等[3-4]提出的光滑粒子流體動(dòng)力學(xué)方法，是目前最廣泛使用的方法。然而SPH方法在超高速碰撞問題中，仍具有一定的缺點(diǎn)和局限性。其主要問題在于插值函數(shù)不滿足Kronecker-delta屬性，無法直接施加唯一邊界條件，缺乏有效的數(shù)值積分方式，具有嚴(yán)重的零能模式以及拉力不穩(wěn)定性問題等。

近十年來發(fā)展了一種新的數(shù)值方法，最早由M.Ortiz等[5-6]提出的一種顯式增量更新拉格朗日無網(wǎng)格計(jì)算方法，最優(yōu)輸運(yùn)無網(wǎng)格方法(Optimal Transportation Meshfree method，OTM)。相比于SPH方法[7]，OTM方法集合了拉格朗日方法和無網(wǎng)格法的所有優(yōu)點(diǎn)：質(zhì)量守恒自動(dòng)滿足；邊界處滿足Kronecker-delta屬性；插值和積分在不同的離散點(diǎn)上進(jìn)行，避免了拉力不穩(wěn)定性問題。本文采用pOTM方法模擬彈丸超高速正撞擊單層靶板過程，通過與文獻(xiàn)[8]結(jié)果進(jìn)行對(duì)比，驗(yàn)證pOTM方法模擬超高速碰撞問題的適用性，同時(shí)分析了不同量級(jí)規(guī)模模型以及不同核心數(shù)對(duì)CPU并行計(jì)算耗費(fèi)時(shí)間的影響。

1 OTM方法基本理論

幾何模型Ω0∝R3所經(jīng)歷的一般力學(xué)過程由變形映射進(jìn)行描述。假設(shè)變形梯度為

F=FeFp

(1)

式中：Fe表示彈性變形；Fp表示塑性變形。因此系統(tǒng)的內(nèi)能密度可表示為

W=We,vol(J)+We,dev(FFp-1)+Wp(Fp)+Wh

(2)

式中：J=detF為變形梯度的Jacobian，We,vol與We,dev分別為彈性響應(yīng)的體積部分和偏離項(xiàng)，Wp是塑性響應(yīng)，Wh是熱儲(chǔ)量。從而可以推出以下變分結(jié)構(gòu)：

(3)

式中：K是動(dòng)能密度；B是單位質(zhì)量體積力；T是施加于邊界Γt上的牽引力。通過對(duì)等效勢(shì)能Φ取變分可以獲得線性動(dòng)量守恒方程和能量守恒方程，即：

δΦ=0

(4)

1) 時(shí)間離散。從變分原理出發(fā),對(duì)變分結(jié)構(gòu)Φ取極值，再進(jìn)行時(shí)間離散。對(duì)于給定的[tn,tn+1]，若tn時(shí)刻的變量{φn}已知，則tn+1時(shí)刻的變量{φn+1}可以通過最小化增量近似獲得，即：

(5)

采用后向差分格式進(jìn)行近似，可以得到以下增量變分結(jié)構(gòu)：

(6)

2) 空間離散。OTM方法采用兩套粒子(物質(zhì)點(diǎn)與節(jié)點(diǎn))對(duì)幾何模型進(jìn)行空間離散。其離散方案如圖1所示。

圖1 物質(zhì)點(diǎn)與節(jié)點(diǎn)空間離散方案示意圖

在tn→tn+1過程中，物質(zhì)點(diǎn)的更新通過鄰域節(jié)點(diǎn)信息進(jìn)行插值獲得。則位移傳輸映射φn→n+1可以表示為

(7)

式中：N表示離散域內(nèi)的總節(jié)點(diǎn)數(shù);φa表示構(gòu)形中節(jié)點(diǎn)的位置；Na(x)表示節(jié)點(diǎn)xa的插值函數(shù)[7]。同樣地，物質(zhì)點(diǎn)xp在tn+1時(shí)刻的位置xp,n+1及其變形梯度Fp,n+1可近似為

Fp,n+1=▽?duì)課→n+1(xp,n)Fp,n

(8)

通過上述無網(wǎng)格離散方案，可得到全局系統(tǒng)的力平衡方程，即：

(9)

因此，節(jié)點(diǎn)xa在tn+1時(shí)刻坐標(biāo)的更新形式為

(10)

式中：Ia,n為節(jié)點(diǎn)xa在tn時(shí)刻的線性動(dòng)量，即：

(11)

通過求解每個(gè)節(jié)點(diǎn)的力平衡方程式(9)，即可獲得系統(tǒng)的整個(gè)位移變形場(chǎng)。

2 并行計(jì)算

在OTM方法求解過程中，其計(jì)算負(fù)荷主要集中在節(jié)點(diǎn)力的求解以及物質(zhì)點(diǎn)鄰域的更新及其相應(yīng)的節(jié)點(diǎn)插值函數(shù)與導(dǎo)數(shù)的計(jì)算。由于其計(jì)算主要集中在物質(zhì)點(diǎn)上進(jìn)行，各個(gè)物質(zhì)點(diǎn)之間沒有相應(yīng)的數(shù)據(jù)交換，所以此方法非常適合進(jìn)行并行化計(jì)算。

pOTM是一種基于多核心并行化架構(gòu)的大規(guī)模并行計(jì)算方法，利用Fortran語言通過MPI[9]編寫了一套基于共享節(jié)點(diǎn)的劃分方案，將求解模型進(jìn)行分布式計(jì)算分解。

2.1 離散域分區(qū)

在pOTM方法中，首先將物質(zhì)點(diǎn)劃分成不同的子集分配到各個(gè)子核心中進(jìn)行分布式計(jì)算，具體劃分如圖2所示。通過計(jì)算tn時(shí)刻PI子核心中各物質(zhì)點(diǎn)的鄰域，得到各個(gè)子核心所需要的鄰域范圍，把這個(gè)稱為Range Box[10]。 Range Box中包含PI子核心中所有物質(zhì)點(diǎn)的鄰域范圍內(nèi)的節(jié)點(diǎn)。不同Range Box之間重疊部分被定義為Shadow Box，處于Shadow Box中的節(jié)點(diǎn)被定義為共享節(jié)點(diǎn)，其數(shù)據(jù)信息存儲(chǔ)于共享數(shù)據(jù)交換表中。

圖2 各核心物質(zhì)點(diǎn)與節(jié)點(diǎn)劃分示意圖

2.2 數(shù)值迭代流程

(12)

根據(jù)更新后物質(zhì)點(diǎn)鄰域Np,k+1，更新各子核心PI其tn+1時(shí)刻的Range Box，然后判斷時(shí)間步，若tn+1=tk，代表計(jì)算結(jié)束，得到t0→tk時(shí)間段內(nèi)物質(zhì)點(diǎn)和節(jié)點(diǎn)的動(dòng)力學(xué)信息；若tn+1≠tk，則進(jìn)行下一步迭代，直至tn+1=tk為止。

pOTM方法的主要求解過程可以描述為：

1) 將d維幾何模型Ω0進(jìn)行單元?jiǎng)澐郑瑢⑵潆x散為一組初始的節(jié)點(diǎn)集{xa,0,a=1,2,…,N}和一組初始的物質(zhì)點(diǎn)集{xp,0,p=1,2,…,M}；

2) 在初始時(shí)刻n=0時(shí)，初始化物質(zhì)點(diǎn)坐標(biāo)xp,n和節(jié)點(diǎn)坐標(biāo)xa,n，將物質(zhì)點(diǎn)集劃分成不同的子集并分配到各個(gè)子核心中，初始化各子核心的Range Box以及相應(yīng)的節(jié)點(diǎn)插值函數(shù)Na(xp,n)以及節(jié)點(diǎn)插值函數(shù)導(dǎo)數(shù)(Na(xp,n)；

4) 根據(jù)節(jié)點(diǎn)坐標(biāo)xn+1，此時(shí)節(jié)點(diǎn)產(chǎn)生位移，同時(shí)物質(zhì)點(diǎn)將產(chǎn)生變形。利用局部最大熵插值函數(shù)計(jì)算物質(zhì)點(diǎn)的局部增量變形梯度Fp,n→n+1=xa,n+1?▽Na(xp,n)和變形梯度Fp,n+1=Fp,n→n+1°Fp,n來計(jì)算物質(zhì)點(diǎn)的變形；

5) 根據(jù)第4步得到物質(zhì)點(diǎn)從tn→tn+1時(shí)刻的運(yùn)動(dòng)和變形，更新物質(zhì)點(diǎn)坐標(biāo)xp,n+1、體積νp,n+1=det (Fp,n→n+1)νp,n以及密度ρp,n+1=mp/vp,n+1，其中mp為物質(zhì)點(diǎn)質(zhì)量；

6) 根據(jù)第4步和第5步得到的數(shù)據(jù)，通過本構(gòu)關(guān)系獲得物質(zhì)點(diǎn)的應(yīng)力σp,n+1，其中應(yīng)力偏張量由Jaumann應(yīng)力率模型給出，應(yīng)力球張量由物質(zhì)點(diǎn)的密度ρp,n+1通過狀態(tài)方程給出；

7) 根據(jù)tn+1時(shí)刻下物質(zhì)點(diǎn)xp,n+1的鄰域以及節(jié)點(diǎn)坐標(biāo)xa,n+1，更新各個(gè)子核心的Range Box，以及其相對(duì)應(yīng)的節(jié)點(diǎn)插值函數(shù)Na(xp,n+1)及其導(dǎo)數(shù)(Na(xp,n+1)；

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

3 彈丸超高速碰撞數(shù)值模擬

3.1 計(jì)算模型及材料模型參數(shù)

采用pOTM方法模擬球形彈丸超高速正撞擊單層靶板的過程。模型結(jié)構(gòu)中球形彈丸的直徑為1 cm，撞擊速度為6.18 km/s，單層靶板的尺寸為5 cm×5 cm×0.4 cm，球形彈丸和單層靶板的材料均為鋁，其密度取為2.78 g/cm3。具體離散方案如圖3所示。

左側(cè)為局部放大圖，空心為物質(zhì)點(diǎn)，實(shí)心為節(jié)點(diǎn)

本構(gòu)模型采用Johnson-Cook材料損傷模型[11]，它可以描述在高溫、高壓、高應(yīng)變率情況下金屬材料的強(qiáng)度極限以及失效過程，即：

(13)

而Johnson-Cook材料損傷模型中的斷裂由以下累積損傷定律導(dǎo)出：

(14)

對(duì)于狀態(tài)方程，本文選用了較為常用的Mie-Grüneisen狀態(tài)方程，其表達(dá)式為

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

(15)

其本構(gòu)模型及狀態(tài)方程的具體參數(shù)如表1和表2。數(shù)值工況參照Hiermaier[8]的實(shí)驗(yàn)數(shù)據(jù)。

3.2 網(wǎng)格收斂性驗(yàn)證

采用pOTM方法模擬球形彈丸超高速正撞擊單層靶板的過程，分析不同網(wǎng)格密度模型規(guī)模下碎片云的形態(tài)。圖4為20 μs時(shí)，不同網(wǎng)格密度模型規(guī)模模擬結(jié)果。

表1 Johnson-Cook損傷模型參數(shù)

表2 Mie-Grüneisen狀態(tài)方程參數(shù)

圖4 不同網(wǎng)格密度模擬結(jié)果

從圖4可以看出： 當(dāng)模型規(guī)模偏少時(shí)，其碎片云形狀比較不規(guī)則；隨著模型規(guī)模的增大，碎片云形狀開始變得更加圓潤，且與實(shí)際效果吻合得更好。

假定布點(diǎn)密度為網(wǎng)格密度模型規(guī)模之比的三次方，且10萬量級(jí)模型規(guī)模的布點(diǎn)密度為1，則可以得出相應(yīng)的布點(diǎn)密度對(duì)碎片云具體參數(shù)的影響，如表3所示。

表3中，l為碎片云長度，w為碎片云寬度，Δ為相對(duì)誤差。由表3可得，對(duì)于布點(diǎn)密度大的模型，碎片云的基本參數(shù)變化不大，誤差在1.5%左右。因此可以認(rèn)為在保證足夠規(guī)模的網(wǎng)格密度下，計(jì)算結(jié)果隨布點(diǎn)密度的增加趨于穩(wěn)定。同時(shí)也可以在保證計(jì)算精度的前提下，適當(dāng)減少布點(diǎn)密度來降低計(jì)算壓力和加快計(jì)算速度。

表3 碎片云的具體參數(shù)

3.3 結(jié)果及分析

通過高性能計(jì)算機(jī)集群進(jìn)行并行模擬計(jì)算，采用24個(gè)核心數(shù)，計(jì)算模型被離散為45 634個(gè)物質(zhì)點(diǎn)和58 916個(gè)節(jié)點(diǎn)，其數(shù)量級(jí)約為10萬量級(jí)，其中彈丸包含5 704個(gè)物質(zhì)點(diǎn)和7 366個(gè)節(jié)點(diǎn)，單層靶板包含39 930個(gè)物質(zhì)點(diǎn)和51 550個(gè)節(jié)點(diǎn)，時(shí)間步長由CFL條件決定。圖5是球形彈丸超高速正撞擊單層靶板過程中，5 μs、10 μs、15 μs以及20 μs時(shí)刻碎片云的演化結(jié)果。圖6是20 μs時(shí)刻碎片云的速度云圖(撞擊方向)和密度云圖。

圖5 pOTM方法模擬碎片云的演化過程

圖6 20 μs時(shí)刻碎片云的速度云圖和密度云圖

取其20 μs時(shí)刻的演化結(jié)果與文獻(xiàn)[8]結(jié)果進(jìn)行對(duì)比(圖7)，本文所模擬得到的結(jié)果，靶板前端有明顯的碎片飛濺，靶板后方的大部分碎片云質(zhì)量集中在碎片云前端中軸線附近，得到的結(jié)果和文獻(xiàn)[8]的結(jié)果吻合較好。

表4為采用pOTM方法模擬彈丸超高速正撞擊單層靶板以及文獻(xiàn)[8]結(jié)果的碎片云具體參數(shù)對(duì)比。

表4中，dc為靶板彈坑直徑，Δ為l/w與實(shí)驗(yàn)結(jié)果的相對(duì)誤差。綜合表3、表4可以得出，pOTM方法模擬結(jié)果碎片云形態(tài)比較圓潤，與文獻(xiàn)[8]相比其相對(duì)誤差比較小，且與實(shí)際結(jié)果吻合較好[12]，比較適合對(duì)超高速碰撞問題進(jìn)行深入研究。

圖7 不同數(shù)值方法模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比

表4 碎片云的模擬結(jié)果

3.4 并行效率分析

采用pOTM方法在高性能計(jì)算機(jī)集群上對(duì)彈丸超高速正撞擊單層靶板進(jìn)行并行模擬研究，對(duì)比了不同核心數(shù)、不同物質(zhì)點(diǎn)數(shù)和節(jié)點(diǎn)數(shù)下的多次并行計(jì)算結(jié)果，分析迭代 1 000 個(gè)循環(huán)所需要耗費(fèi)的CPU時(shí)間，得到運(yùn)行時(shí)間隨核心數(shù)變化情況。圖8為58 916個(gè)節(jié)點(diǎn)以及45 634個(gè)物質(zhì)點(diǎn)迭代1 000個(gè)時(shí)間步耗費(fèi)的時(shí)間隨核心數(shù)的變化情況，其數(shù)量級(jí)約為十萬量級(jí)。圖9為638 196個(gè)節(jié)點(diǎn)以及50 3762個(gè)物質(zhì)點(diǎn)，其數(shù)量級(jí)為約百萬量級(jí)。

隨著核心數(shù)的增加，計(jì)算所耗費(fèi)的時(shí)間將逐漸減小，但是核心過多會(huì)導(dǎo)致MPI通信時(shí)間增加，從而導(dǎo)致整體計(jì)算時(shí)間增加。因此存在一個(gè)核心數(shù)使得并行計(jì)算耗費(fèi)的時(shí)間最少。

對(duì)于十萬量級(jí)的計(jì)算，其串行迭代1 000個(gè)時(shí)間步需要耗費(fèi)7 975 s，使用24個(gè)核心進(jìn)行并行計(jì)算，迭代1 000個(gè)時(shí)間步需要耗費(fèi)624 s，其加速比[13]約為7 975/6 240=12.79，效率約為12.79/24=53.3%。

圖8 十萬量級(jí)計(jì)算時(shí)間隨核心數(shù)變化圖

圖9 百萬量級(jí)計(jì)算時(shí)間隨核心數(shù)變化圖

對(duì)于百萬量級(jí)的計(jì)算，使用24個(gè)核心進(jìn)行并行計(jì)算，迭代1 000個(gè)時(shí)間步需要耗費(fèi)7 336 s，其加速比約為11.5，使而用28個(gè)核心需要耗費(fèi)7 280 s，其加速比約為11.6。因此，雖然繼續(xù)增加核心數(shù)，耗費(fèi)的時(shí)間可以繼續(xù)減少，但是減少的幅度變得很小，沒有了實(shí)際意義。

因此對(duì)于十萬量級(jí)的來說，采用24個(gè)核心數(shù)計(jì)算求解時(shí)，其加速比可以達(dá)到12.79，且隨著核心數(shù)的增加，加速比不會(huì)產(chǎn)生明顯的提升。對(duì)于百萬量級(jí)，采用28個(gè)核心時(shí)，其加速比為11.6，其并行所耗費(fèi)的時(shí)間可以降到最低。

從上述2個(gè)結(jié)果分析可以得出，如果模型規(guī)模在百萬量級(jí)之下時(shí)，可以提供最多為13的加速比，即并行計(jì)算所耗費(fèi)的時(shí)間為串行計(jì)算所耗費(fèi)時(shí)間的1/13。因此本并行程序可以很大程度上節(jié)省CPU計(jì)算時(shí)間。

4 結(jié)論

1) 在彈丸超高速正撞擊單層板把的數(shù)值模擬中，采用pOTM方法模擬所得到的結(jié)果與文獻(xiàn)[8]中的實(shí)驗(yàn)和模擬結(jié)果吻合較好。結(jié)果顯示碎片在單層靶板前端有明顯的飛濺，靶板后方的大部分碎片云質(zhì)量集中在碎片云前端中軸線附近，且靶板后方外圍的碎片云質(zhì)量分布比較均勻。彈丸超高速正撞擊單層靶板所形成的彈坑直徑以及碎片云的具體參數(shù)均與文獻(xiàn)[8]中的結(jié)果相差無幾，說明pOTM方法很適合來進(jìn)行求解超高速碰撞相關(guān)的數(shù)值模擬問題；

2) 針對(duì)不同量級(jí)模型規(guī)模以及不同核心數(shù)對(duì)并行效率的影響，對(duì)十萬量級(jí)以及百萬量級(jí)的模型規(guī)模進(jìn)行了分析，取其迭代1 000個(gè)時(shí)間步所耗費(fèi)的CPU時(shí)間進(jìn)行比較，結(jié)果顯示，計(jì)算所耗費(fèi)的時(shí)間會(huì)隨著核心數(shù)的增加逐漸減小，但是核心過多也會(huì)導(dǎo)致通信時(shí)間的增加，因此存在最優(yōu)核心數(shù)，對(duì)于百萬量級(jí)28核心數(shù)為其最優(yōu)核心數(shù)，其加速比可以達(dá)到12，可以很大程度上節(jié)省CPU計(jì)算時(shí)間。