張曉天，賈光輝，黃 海

（北京航空航天大學(xué)宇航學(xué)院，北京100083）

1 引 言

超高速碰撞數(shù)值模擬技術(shù)是帶有破碎環(huán)節(jié)的沖擊動(dòng)力學(xué)數(shù)值模擬方法的重要應(yīng)用。超高速碰撞問(wèn)題的特點(diǎn)是沖擊速度快、材料變形大且破碎程度大，碎片往往以云狀出現(xiàn)，具有高度非線性。

沖擊動(dòng)力學(xué)數(shù)值模擬方法主要包括Lagrange有限元方法、SPH 無(wú)網(wǎng)格方法和Euler有限元方法。傳統(tǒng)Lagrange有限元方法在模擬大變形和破碎問(wèn)題中會(huì)產(chǎn)生嚴(yán)重的網(wǎng)格畸變，使得計(jì)算難以進(jìn)行。斷裂侵蝕機(jī)制用刪除大變形單元的方法解決了網(wǎng)格畸變的問(wèn)題。但是在超高速碰撞問(wèn)題中，由于變形區(qū)域較大，會(huì)有大量的單元被刪除。這妨礙了破碎數(shù)值模擬中碎片的形成，難以模擬碎片云的演化。另外，沖擊動(dòng)能以勢(shì)能的形式存入大變形單元中并進(jìn)一步傳播，刪除大量單元的同時(shí)使得系統(tǒng)總能量大幅損失，造成變形模擬的失真，如圖1。SPH 方法利用離散的粒子群模擬連續(xù)的物質(zhì)，使用無(wú)網(wǎng)格技術(shù)解決了網(wǎng)格畸變問(wèn)題，在碎片云模擬方面有突出的優(yōu)勢(shì)［1－2］，是目前使用最廣泛的超高速碰撞數(shù)值模擬方法。但是SPH 方法由于要進(jìn)行粒子搜索，計(jì)算效率低于有限元方法［3］。同時(shí)，該算法的受拉不穩(wěn)定性會(huì)影響結(jié)果的準(zhǔn)確度［4］。由于使用無(wú)網(wǎng)格技術(shù)，物質(zhì)邊界不明確是SPH 方法的又一不足。Euler有限元方法使用空間網(wǎng)格，網(wǎng)格不隨物質(zhì)變形而變形，因此不存在網(wǎng)格畸變問(wèn)題，沖擊模擬計(jì)算可以有效進(jìn)行，但是同樣難以刻畫(huà)物質(zhì)邊界。鑒于現(xiàn)有方法仍存在不足，不斷改進(jìn)現(xiàn)有方法和探索新的思路仍是沖擊破碎模擬算法研究的熱點(diǎn)。

本文中，對(duì)Lagrange有限元方法進(jìn)行改進(jìn)，用節(jié)點(diǎn)分離方法和畸變侵蝕方法分別解決斷裂機(jī)制問(wèn)題和網(wǎng)格畸變問(wèn)題，以替代斷裂侵蝕機(jī)制。

圖1 Lagrange算法算例［1］Fig.1 Asimulation example of the Lagrange method［1］

2 節(jié)點(diǎn)分離方法原理

節(jié)點(diǎn)分離方法是對(duì)傳統(tǒng)Lagrange單元的一種改進(jìn)，它將節(jié)點(diǎn)按參與構(gòu)成的單元的個(gè)數(shù)進(jìn)行復(fù)制，為這些單元各分配一個(gè)空間位置上重疊的節(jié)點(diǎn)，然后用斷裂準(zhǔn)則對(duì)所有空間位置重疊的節(jié)點(diǎn)族分別“捆綁”（約束相對(duì)自由度），當(dāng)某節(jié)點(diǎn)的力學(xué)狀態(tài)滿足斷裂準(zhǔn)則時(shí)，此節(jié)點(diǎn)所在節(jié)點(diǎn)族解離，兩兩節(jié)點(diǎn)間不再傳力，以此方式來(lái)描述材料的斷裂特性。

節(jié)點(diǎn)分離方法由有限元前處理中的模型節(jié)點(diǎn)復(fù)制、初始約束添加和計(jì)算中的約束解除等3個(gè)步驟組成。

圖2給出了一個(gè)前處理階段進(jìn)行模型的節(jié)點(diǎn)復(fù)制的示意圖。首先，建立一般的Lagrange有限元網(wǎng)格，然后將節(jié)點(diǎn)N 復(fù)制為N1、N2、N3、N4，分別分配給共用節(jié)點(diǎn)N 的4個(gè)實(shí)體單元。然后，為新創(chuàng)建的4個(gè)節(jié)點(diǎn)添加約束，使之擁有一致的自由度。4個(gè)節(jié)點(diǎn)真實(shí)的空間位置是重合的，圖2中人為將4個(gè)節(jié)點(diǎn)分離開(kāi)來(lái)，只為方便表達(dá)。

圖2 節(jié)點(diǎn)復(fù)制Fig.2 Replication of nodes

3 畸變侵蝕方法

節(jié)點(diǎn)分離概念解決了網(wǎng)格斷裂的問(wèn)題，這是斷裂分析的基本——斷裂機(jī)制，但是Lagrange有限元網(wǎng)格在大變形下發(fā)生畸變的弊端并未解決。傳統(tǒng)的Lagrange有限元斷裂算法使用斷裂侵蝕機(jī)制，將滿足斷裂條件的所有單元?jiǎng)h除。這種做法一方面實(shí)現(xiàn)了網(wǎng)格斷裂，另一方面處理了畸變單元。節(jié)點(diǎn)分離概念將滿足斷裂條件的節(jié)點(diǎn)組分離，并不刪除單元，但是網(wǎng)格畸變并沒(méi)有解決?；兦治g方法根據(jù)單元當(dāng)前的幾何形狀判斷是否發(fā)生畸變，并將滿足畸變判別準(zhǔn)則的單元?jiǎng)h除。由于六面體單元的計(jì)算精度高于四面體，本文的討論僅限于六面體單元。單元畸變分為3種模式：?jiǎn)蜗虺叨犬惓Ｔ龃螅Ｊ紸）、體積異常增大（模式B）和自穿透（模式C），見(jiàn)圖3。

（1）畸變模式A 的判別

將一個(gè)六面體單元的8個(gè)節(jié)點(diǎn)分別向慣性坐標(biāo)系投影，得到每個(gè)節(jié)點(diǎn)的坐標(biāo)。畸變判別準(zhǔn)則為

式中：Lmax為單向尺度閾值，Lch是劃分網(wǎng)格時(shí)給定的單元特征長(zhǎng)度，r1是單向增大閾值系數(shù)。如果沿某個(gè)慣性系坐標(biāo)軸方向投影尺度超過(guò)此閾值，則認(rèn)為該單元發(fā)生畸變。

（2）畸變模式B的判別

求解某時(shí)刻某1個(gè)單元在慣性系3個(gè)坐標(biāo)軸上的投影尺寸，進(jìn)而將3個(gè)投影尺寸的乘積作為當(dāng)前單元體積的簡(jiǎn)化度量。如果滿足下式，則認(rèn)為該單元發(fā)生畸變

式中：r2為體積擴(kuò)大倍數(shù)閾值。

（3）畸變模式C的判別

對(duì)于網(wǎng)格自穿透問(wèn)題，需要判斷某1個(gè)單元在任意1個(gè)節(jié)點(diǎn)處是否發(fā)生自穿透，如果有1個(gè)節(jié)點(diǎn)處發(fā)生，則認(rèn)為整個(gè)單元發(fā)生畸變。如圖4所示（圖中節(jié)點(diǎn)對(duì)應(yīng)關(guān)系見(jiàn)圖3）。

對(duì)于節(jié)點(diǎn)B 來(lái)說(shuō)，初始時(shí)相連的3條邊為BA、BC、BF。令3個(gè)坐標(biāo)軸與3條邊固連，可以建立右手笛卡爾坐標(biāo)系B－ACF，3個(gè)基矢量記為｛e1，e2，e3｝。在沖擊過(guò)程中某時(shí)刻，該坐標(biāo)系的3個(gè)基矢量在慣性系O－xyz 中的坐標(biāo)形式可以由節(jié)點(diǎn)坐標(biāo)定量計(jì)算得出，記為。該坐標(biāo)系一般來(lái)說(shuō)是笛卡爾斜角坐標(biāo)系。由于慣性系是單位正交直角坐標(biāo)系，因此這個(gè)斜角坐標(biāo)系相對(duì)于慣性系的坐標(biāo)變換矩陣Jacobian就是由單位正交基到｛e＇1，e＇2，e＇3｝的變換矩陣。而Jacobian的行列式值則反映了斜角坐標(biāo)系B′－ACF 各基矢量之間耦合的程度。耦合程度越大，說(shuō)明角變形越大。如果Jacobian的行列式值由正變負(fù)，則說(shuō)明該坐標(biāo)系由右手系變?yōu)樽笫窒担▓D4中即是），這是由于節(jié)點(diǎn)F 穿透了面ABCD 造成的，即發(fā)生了自穿透。按此方法遍歷該單元的所有節(jié)點(diǎn)，即可判斷整個(gè)單元是否發(fā)生了自穿透。

圖3 六面體單元的3種畸變模式Fig.3 Three distortion modes of hexahedron elements

圖4 網(wǎng)格自穿透的判別Fig.4Judgment of mesh self－piercing

將3種畸變模式綜合，某單元只要發(fā)生了其中1種，即認(rèn)為該單元發(fā)生畸變，予以刪除。由于單元畸變的變形程度要大于單元斷裂時(shí)的變形程度，所以使用節(jié)點(diǎn)分離方法加畸變侵蝕方法雖然與斷裂侵蝕方法一樣要?jiǎng)h除單元，但是刪除的數(shù)目要少很多，保留了那些大變形但未畸變的單元，保留了系統(tǒng)能量，在常規(guī)網(wǎng)格密度下可以提高計(jì)算結(jié)果的真實(shí)性。

節(jié)點(diǎn)分離只是一種有限元斷裂數(shù)值模擬的概念，關(guān)于這種概念應(yīng)用的文獻(xiàn)較少，實(shí)現(xiàn)這個(gè)概念的具體技術(shù)途徑也不盡相同，對(duì)應(yīng)了不同的具體方法。有將節(jié)點(diǎn)分離概念用于模擬低速撞擊、膨脹破裂和侵徹問(wèn)題的［5－7］，但還沒(méi)有將節(jié)點(diǎn)分離方法用于超高速碰撞數(shù)值模擬的。

4 節(jié)點(diǎn)分離方法基于LS－dyna的軟件實(shí)現(xiàn)

本文中使用C＋＋編程調(diào)用LS－dyna求解器的技術(shù)途徑實(shí)現(xiàn)節(jié)點(diǎn)分離算法。首先，使用通用有限元軟件建立原始有限元網(wǎng)格。然后，通過(guò)自編程對(duì)原始網(wǎng)格進(jìn)行節(jié)點(diǎn)復(fù)制、節(jié)點(diǎn)集建立和約束添加，并通過(guò)添加材料、接觸、時(shí)間步、輸出等LS－dyna的K 文件配置關(guān)鍵字生成節(jié)點(diǎn)分離方法的計(jì)算K 文件。之后，流程控制由自編程序來(lái)完成。在LS－dyna中，節(jié)點(diǎn)集的創(chuàng)建使用＊SET＿NODE＿LIST 關(guān)鍵字實(shí)現(xiàn)，而通過(guò)＊CONSTRAINED＿TIED＿NODES＿FAILURE 關(guān)鍵字實(shí)現(xiàn)節(jié)點(diǎn)集的約束和斷裂應(yīng)變閾值達(dá)到時(shí)節(jié)點(diǎn)集的自動(dòng)解離。

程序調(diào)用LS－dyna求解器，進(jìn)行顯式積分計(jì)算至T 時(shí)刻終止，在求解結(jié)束后調(diào)用LS－prepost后處理軟件讀取d3plot計(jì)算結(jié)果，并將節(jié)點(diǎn)單元信息輸出為文本文件。C＋＋程序讀取文本文件中的計(jì)算結(jié)果數(shù)據(jù)，針對(duì)數(shù)據(jù)進(jìn)行畸變侵蝕分析，確定要?jiǎng)h除的嚴(yán)重畸變單元。程序創(chuàng)建LS－dyna重啟動(dòng)文件，將要?jiǎng)h除的單元，按＊DELETE＿ELEMENT＿SOLID 關(guān)鍵字的格式寫(xiě)入LS－dyna重啟動(dòng)文件中，并將2T 作為下一個(gè)計(jì)算終止時(shí)刻寫(xiě)入重啟動(dòng)文件。之后再次調(diào)用LS－dyna求解器進(jìn)行重啟動(dòng)分析，求解至2T 時(shí)刻終止，以此類(lèi)推循環(huán)進(jìn)行。流程如圖5所示。

在這個(gè)流程中，節(jié)點(diǎn)分離模型的創(chuàng)建由通用前處理軟件配合自編代碼再輔助以人工完成，而使用自編程序調(diào)用LS－dyna頻繁進(jìn)行重啟動(dòng)分析的目的在于每隔T 時(shí)刻，停止計(jì)算并進(jìn)行畸變侵蝕分析，刪除嚴(yán)重畸變單元。而T 即為畸變侵蝕分析步長(zhǎng)，這個(gè)步長(zhǎng)并不需要等于LS－dyna求解器的顯式積分步長(zhǎng)，只要保證適時(shí)刪除畸變單元保證計(jì)算正常有效進(jìn)行即可。

圖6給出了節(jié)點(diǎn)分離加畸變侵蝕方法對(duì)球形鋁彈丸超高速撞擊鋁板的數(shù)值模擬算例，可見(jiàn)本方法具有碎片云模擬能力。

圖5 計(jì)算流程圖Fig.5 Computation flow chart

圖6 節(jié)點(diǎn)分離方法碎片云Fig.6 Debris simulated by the node－separation method

5 實(shí)驗(yàn)驗(yàn)證

圖7為文獻(xiàn)［8］中的3層板防護(hù)結(jié)構(gòu)超高速撞擊實(shí)驗(yàn)的示意圖。彈丸為Al 2024，防護(hù)屏均為Al 6061。彈丸直徑7.9mm，撞擊速度為5.5km／s，結(jié)構(gòu)布局參數(shù)如圖7所示。

使用節(jié)點(diǎn)分離的Lagrange有限元算法和SPH無(wú)網(wǎng)格算法對(duì)此實(shí)驗(yàn)進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)結(jié)果以及手冊(cè)中提供Euler有限元算法結(jié)果［8］進(jìn)行對(duì)比。材料模型均為Johnson－Cook材料模型，材料參數(shù)見(jiàn)表1；狀態(tài)方程為Gruneison狀態(tài)方程，材料參數(shù)分別為［11］：γ＝1.97，c1＝5.386km／s，s1＝1.339，T0＝300K，c＝884J／（kg·K）。材料的斷裂應(yīng)變：Al 2024取0.8，Al 6061取1.0。

圖7 實(shí)驗(yàn)工況Fig.7 Experiment configuration

表1 Johnson－Cook材料參數(shù)Table 1Johnson－Cook material parameters

圖8為撞擊后25μs時(shí)3種算法計(jì)算結(jié)果對(duì)比。其中，圖8（a）為文獻(xiàn)［8］中提供的Ouranos軟件計(jì)算的結(jié)果，該軟件基于Euler 2D 算法；圖8（b）為本文中提出的節(jié)點(diǎn)分離加畸變侵蝕方法的計(jì)算結(jié)果；圖8（c）～（d）為本文中使用LS－dyna的SPH 3D 算法的計(jì)算結(jié)果。在后3個(gè)算例中，均使用雙CPU 并行計(jì)算。

表2給出了3層板穿孔情況、計(jì)算規(guī)模和計(jì)算時(shí)間。d1、d2為前、中板穿孔直徑，ε1、ε2為誤差。由穿孔直徑可知，本文算法具有較高精度，明顯優(yōu)于Ouranos方法結(jié)果，略優(yōu)于SPH－2結(jié)果。在同等計(jì)算規(guī)模（節(jié)點(diǎn)數(shù)目）下，本文方法和SPH 方法計(jì)算花費(fèi)時(shí)間相當(dāng)，但是本文方法結(jié)果明顯優(yōu)于SPH 方法。另外，SPH 方法的結(jié)果強(qiáng)依賴于粒子密度，SPH－2算例結(jié)果明顯優(yōu)于SPH－1結(jié)果且與節(jié)點(diǎn)分離方法精度相當(dāng)；然而，SPH－2計(jì)算花費(fèi)時(shí)間是節(jié)點(diǎn)分離方法的10倍以上。

圖8 計(jì)算結(jié)果對(duì)比Fig.8 Comparison of results

表2 穿孔直徑對(duì)比Table 2 Comparison of perforation diameters

6 結(jié) 論

提出了節(jié)點(diǎn)分離加畸變侵蝕的Lagrange有限元數(shù)值模擬方法，給出了基于LS－dyna軟件的一種實(shí)現(xiàn)途徑。使用該方法對(duì)超高速碰撞問(wèn)題進(jìn)行了數(shù)值模擬，與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果、數(shù)值模擬結(jié)果和SPH方法模擬結(jié)果做了對(duì)比，表明了該方法的可行性和有效性。

使用節(jié)點(diǎn)分離的Lagrange有限元方法模擬超高速碰撞問(wèn)題，有效利用了Lagrange有限元方法的各種優(yōu)勢(shì)，如計(jì)算速度較快、穩(wěn)定性好、物質(zhì)邊界明確等。因此，這種方法在超高速碰撞數(shù)值模擬方面，可以成為SPH 無(wú)網(wǎng)格方法的有效補(bǔ)充和替代。

［1］ 樂(lè)莉，閆軍，鐘秋海.超高速撞擊仿真算法分析［J］.系統(tǒng)仿真學(xué)報(bào)，2004，16（9）：1941－1943.YUE Li，YAN Jun，ZHONG Qiu－h(huán)ai.Simulations of debris impacts using three different algorithms［J］.Journal of System Simulation，2004，16（9）：1941－1943.

［2］ 賈光輝，黃海，胡震東.超高速撞擊數(shù)值仿真結(jié)果分析［J］.爆炸與沖擊，2005，25（1）：47－53.JIA Guang－h(huán)ui，HUANG Hai，HU Zhen－dong.Simulation analyse of hypervelocity impact perforation［J］.Explosion and Shock Waves，2005，25（1）：47－53.

［3］ 王吉，王肖鈞，卞梁.光滑粒子法與有限元的耦合算法及其在沖擊動(dòng)力學(xué)中的應(yīng)用［J］.爆炸與沖擊，2007，27（6）：522－528.WANG Ji，WANG Xiao－jun，BIAN Liang.Linking of smoothed particle hydrodynamics method to standard finite element method and its application in impact dynamics［J］.Explosion and Shock Waves，2007，27（6）：522－528.

［4］ Liu G R，Liu M B.Smoothed particle hydrodynamics：A meshfree particle method［M］.Singapore：World Scientific Publishing Company，2003：30－32.

［5］ Norman F K Jr.Comparison of two modeling approaches for thin－plate penetration simulation［C］∥6th International LS－DYNA Users Conference.Detroit，MI，USA，2000.

［6］ Lee M，Kim E Y，Yoo Y H.Simulation of high speed impact into ceramic composite systems using cohesive－law fracture model［J］.International Journal of Impact Engineering，2008，35（12）：1636－1641.

［7］ 高重陽(yáng).內(nèi)部爆炸載荷下柱殼結(jié)構(gòu)破裂問(wèn)題的研究［D］.北京：清華大學(xué)，2000.

［8］ IADC WG3members.Protection manual（AI 20.1）［EB／OL］.Inter－Agency Space Debris Coordination Committee，2010－07.http：／／www.iadc－online.org／index.cgi？item＝docs＿pub.

［9］ Johnson G R，Cook W H.Selected Huguenots：EOS［C］∥7th International Ballistics.The Hague，Netherlands，1983.

［10］ Fish J.AL 6061－T6－elastomer impact simulations［EB／OL］.The Scientific Computation Research Center.http：／／www.scorec.rpi.edu／REPORT／2005－11.pdf.

［11］ 張慶明，黃風(fēng)雷.超高速碰撞動(dòng)力學(xué)引論［M］.北京：科學(xué)出版社，2000.