基于節(jié)點分離Lagrange有限元方法的超高速碰撞碎片云數(shù)值模擬＊

張曉天，賈光輝，黃 海

（北京航空航天大學宇航學院，北京100083）

1 引 言

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

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

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

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

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

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

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

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

圖2 節(jié)點復制Fig.2 Replication of nodes

3 畸變侵蝕方法

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

（1）畸變模式A 的判別

將一個六面體單元的8個節(jié)點分別向慣性坐標系投影，得到每個節(jié)點的坐標?；兣袆e準則為

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

（2）畸變模式B的判別

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

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

（3）畸變模式C的判別

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

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

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

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

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

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

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

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

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

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

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

圖5 計算流程圖Fig.5 Computation flow chart

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

5 實驗驗證

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

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

圖7 實驗工況Fig.7 Experiment configuration

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

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

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

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

表2 穿孔直徑對比Table 2 Comparison of perforation diameters

6 結(jié) 論

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

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

［1］ 樂莉，閆軍，鐘秋海.超高速撞擊仿真算法分析［J］.系統(tǒng)仿真學報，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］ 王吉，王肖鈞，卞梁.光滑粒子法與有限元的耦合算法及其在沖擊動力學中的應用［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］ 高重陽.內(nèi)部爆炸載荷下柱殼結(jié)構(gòu)破裂問題的研究［D］.北京：清華大學，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］ 張慶明，黃風雷.超高速碰撞動力學引論［M］.北京：科學出版社，2000.