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

張曉天，賈光輝，黃 海

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

1 引 言

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

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

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

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

2 節點分離方法原理

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

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

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

圖2 節點復制Fig.2 Replication of nodes

3 畸變侵蝕方法

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

（1）畸變模式A 的判別

將一個六面體單元的8個節點分別向慣性坐標系投影，得到每個節點的坐標。畸變判別準則為

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

（2）畸變模式B的判別

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

式中：r2為體積擴大倍數閾值。

（3）畸變模式C的判別

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

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

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

圖4 網格自穿透的判別Fig.4Judgment of mesh self－piercing

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

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

4 節點分離方法基于LS－dyna的軟件實現

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

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

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

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

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

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

5 實驗驗證

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

使用節點分離的Lagrange有限元算法和SPH無網格算法對此實驗進行數值模擬，并與實驗結果以及手冊中提供Euler有限元算法結果［8］進行對比。材料模型均為Johnson－Cook材料模型，材料參數見表1；狀態方程為Gruneison狀態方程，材料參數分別為［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材料參數Table 1Johnson－Cook material parameters

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

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

圖8 計算結果對比Fig.8 Comparison of results

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

6 結 論

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

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

［1］ 樂莉，閆軍，鐘秋海.超高速撞擊仿真算法分析［J］.系統仿真學報，2004，16（9）：1941－1943.YUE Li，YAN Jun，ZHONG Qiu－hai.Simulations of debris impacts using three different algorithms［J］.Journal of System Simulation，2004，16（9）：1941－1943.

［2］ 賈光輝，黃海，胡震東.超高速撞擊數值仿真結果分析［J］.爆炸與沖擊，2005，25（1）：47－53.JIA Guang－hui，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］ 高重陽.內部爆炸載荷下柱殼結構破裂問題的研究［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.