應用構型一致的顯式方法預測結構沖擊動響應

初東陽，戎宇飛，周章濤，汪俊，王海坤

（1.中國船舶科學研究中心，江蘇 無錫 214082；2.深海技術科學太湖實驗室，江蘇 無錫 214082）

沖擊碰撞及侵徹對海洋、航天航空裝備的結構安全造成嚴重威脅，而數值仿真是預測在該類載荷條件下裝備結構動態響應的重要技術手段。在沖擊碰撞及侵徹載荷下，裝備結構主要表現出局部大變形、損傷斷裂失效、彈塑性瞬態接觸等強非線性行為[1-2]，這為該類問題的有限元模擬帶來極大的困難。在結構有限元方面，固定網格的拉格朗日法由于跟蹤材料點便于處理固體材料本構及邊界條件，同時又具有較小的計算復雜度，因而得到普遍應用。目前，廣泛流行的沖擊動力學商業軟件如LS-DYNA[3]、ABAQUS[4]等，多使用固定網格的更新拉格朗日格式，其方便實現率相關本構模型，同時又與顯式中心差分方法相結合，實現應力波傳播捕捉，且極大地提高了計算效率[5-8]。然而，大量的實踐證明，傳統的拉格朗日格式在模擬高速碰撞沖擊等涉及到強幾何、材料及接觸非線性問題時往往遇到問題，極易產生網格畸變或時間增量步大大減小，最終使計算無法進行[9]。因此，針對該類問題的模擬，開始采用任意朗格朗日歐拉方法、無網格方法等以彌補固定網格有限元缺陷[10-13]。另一方面，Zhou 等[14]對經典顯式拉格朗日格式在求解強非線性問題產生網格畸變的原因進行了深入研究，認為在動量方程數值求解中，對未知構型變量使用中間構型量替代造成精度損失是誘發單元畸變的重要因素。為了提高求解的嚴謹性，任意參考構型（ARC）理論被提出[15-16]，其不再回避在數值計算中必須使用變量中間構型值的事實，將中間構型納入有限變形理論，引入了參考中間構型的應力應變度量，將理論與數值相統一，極大增強了固定網格拉格朗日格式的健壯性。

為了改善經典結構顯式動力學算法，結合更新拉格朗日格式、蛙跳式顯式中心差分法、罰函數法等模擬大變形碰撞及侵徹問題所面臨的網格畸變問題，本文基于任意參考構型理論，建立了構型一致的動量方程拉格朗日有限元格式，更準確地描述結構大變形。通過結合前增量位移的時間積分方法和共軛梯度算法，使得顯式增量步更新后，接觸非侵入條件仍能夠滿足，避免了實體間相互侵入。通過與商業軟件LS-DYNA 比較對經典Taylor 桿高速撞擊及侵徹問題的模擬結果，驗證了基于構型一致的顯式算法的可靠性。同時，其可有效抑制大變形碰撞過程的單元畸變，提高計算程序的健壯性和預測準確度，使該算法在模擬大變形碰撞問題方面具有優勢。

1 理論方程及數值算法

1.1 控制方程

將變形體在t時刻的構型表示為?t，其邊界包含力邊界與速度邊界，且二者滿足如圖1 所示。對該變形體，其大變形問題的動量方程可以當前構型為參考構型表達為：

圖1 沖擊碰撞問題示意圖Fig.1 Sketch of impact and collision problem

當考慮到該變形體與其他變形體間存在接觸碰撞，力邊界中還包含接觸邊界令g為兩變形體接觸面上點對間的相對間隙，其法向與切向分量分別為gn、gt，則g=gn+gt，gn·gt=0。沿接觸面的法向非侵入條件可表示為：

設sn為接觸面的法向力，則法向接觸 Kuhn-Tucker 條件[17]寫為：

1.2 構型一致的方程弱形式

動量方程的積分弱形式為：

由于控制方程參考的是當前構型，式（4）中的微分、積分需在當前構型上進行，而實際計算時當前構型是未知的，故無法對式（4）的進行精確計算。現有經典顯式算法多對相關變量用中間構型的值進行近似，如LS-DYNA 中積分時采用的變形率值是其介于當前構型與上一時刻構型間的中間構型的值[3]。采用不一致構型的物理量計算會引入誤差，在單元大變形及接觸強非線性情況下，極易誘發單元畸變。為了提高式（4）的求解精度，根據任意參考構型（ARC）理論現將其映射到參考任意構型[15]，其形式可表達為：

1.3 本構模型

考慮到材料發生塑性變形，變形率D可分解為彈性和塑性部分：D=De+Dp。采用由Cauchy 應力的Truesdell 應力率表示次彈性本構：

式中：C?T為該應力率對應的材料剛度矩陣。這里塑性流動滿足J2 流動法則，采用Johnson-cook 屈服條件[18]：

材料失效采用考慮應力三軸度、應變率效應及熱軟化的Johnson-Cook 準則，臨界失效應變為[18]：

1.4 顯式算法

對該方程求解時需滿足tn+2時刻接觸約束條件，即：

以上計算方法在自主開發的結構有限元求解器FAT 中實現，其主要包括初始化、核心求解器及結果輸出部分。通過Hypermesh 等軟件輔助完成前處理生成固定格式的數據輸入文件，后處理則在Tecplot 軟件中進行。

2 數值算例

2.1 Taylor 桿撞擊

一長為100 mm、半徑為10 mm 的金屬圓桿以一定軸向初速度撞擊尺寸為89 mm×89 mm×80 mm、無約束的金屬塊狀靶體。使用本文方法與LS-DYNA軟件對該過程進行模擬，并比較計算結果。撞擊桿與靶體的材料均使用簡單的理想塑性模型，具體材料屬性見表1。采用八節點六面體減縮積分單元進行離散，桿與靶體網格特征長度分別取為2、2.8 mm。LS-DYNA 軟件中的接觸模擬采用*CONTACT_ERODING_SURFACE_TO_SURFACE 命令默認設置。以下為撞擊桿設置不同的初始速度進行數值計算及結果分析。

表1 撞擊桿及靶體材料參數Tab.1 Material parameters of impact bar and target

首先將撞擊桿的初速度設為400 m/s，通過本文方法和LS-DYNA 軟件計算的100 μs 時桿的變形狀態，結果如圖2 所示。桿的碰撞前端產生局部塑性大變形而出現明顯的橫向墩粗，二者計算結果完全吻合。通過撞擊桿上不同位置處質點的速度歷史對本文方法與LS-DYNA 的計算結果進行定量比較，如圖3所示，其中d表示質點到碰撞面的垂直距離。由圖3可知，二者對撞擊桿的動響應過程模擬結果完全一致，尤其是碰撞沖擊前期，一系列小的速度波動也均能得到合理捕捉。

圖2 撞擊桿初速度400 m/s、對桿在100 μs 時的變形狀態模擬結果Fig.2 Numerical results for deformation of bar at 100 μs in the case of initial velocity being 400 m/s:a) plastic strain;b) vertical displacement

將撞擊桿初速度設置為800 m/s，再次用本文方法與LS-DYNA 進行計算，二者的計算結果如圖4 和圖5 所示。可見，相較于初速度為400 m/s，初速度為800 m/s 時，桿碰撞前端局部塑性變形更加顯著。高速碰撞導致靶體變形產生明顯彈坑，桿體頭部沿彈坑邊界發生大的壓剪塑性變形，形成蘑菇蓋狀鈍化。大的塑性流動使桿體前端局部網格發生大的橫向伸長以及縱向縮短，嚴重影響網格質量。在不考慮損傷失效的情況下，撞擊桿單元被壓剪到一定程度時，極易產生單元畸變，使得計算無法進行。通過圖4 可知，本文方法能較好地捕捉到撞擊桿的沖擊鈍化過程，在大變形碰撞下，單元抗畸變能力顯然要高于 LSDYNA。LS-DYNA 使用罰函數接觸算法允許實體間相互侵入，同時顯式算法中物理量采用中間構型值引入了計算誤差，二者疊加，降低了計算精度，使其在計算該大變形碰撞問題早期便發生單元畸變，無法得到正確結果。通過該算例可以驗證構型一致的顯式計算方法的健壯性，也說明其在模擬大變形碰撞問題方面具有優勢。

圖4 撞擊桿初速度800m/s、通過本文方法獲得的80 μs 時刻塑性應變分布Fig.4 Plastic strain distribution at 80 μs obtained from our method in the case of initial velocity being 800 m/s:a) deformation of impact bar and the target;b) deformation of single bar;c) profile plot of deformed bar

圖5 撞擊桿初速度800 m/s，對不同時刻局部變形的模擬結果Fig.5 Numerical results for local deformation at different moments in the case of initial velocity being 800 m/s

2.2 高速侵徹加筋板

使用本文方法與LS-DYNA 軟件對子彈高速垂直侵徹加筋板的過程進行模擬，并比較計算結果。初始時刻，子彈的初速度為800 m/s，其右邊緣與加筋板的軸板對齊。加筋板與子彈材料均為高強度鋼，強度及失效模型為Johnson-Cook 模型，材料參數見表2。采用八節點六面體減縮積分單元進行離散。該算例涉及接觸非線性、幾何大變形及從三軸壓到三軸拉多種應力狀態下的材料損傷失效，可考察計算方法對復雜瞬態沖擊非線性問題的仿真能力。

表2 子彈及加筋板材料參數Tab.2 Material parameters of bullet and reinforced plate

本文方法與LS-DYNA 軟件對250 μs 時刻侵徹狀態的模擬結果如圖6 所示。在單核計算的情況下，二者的計算耗時分別為1.8 h 與0.45 h，本文方法計算成本有所增加。一方面與更精細的大變形及接觸計算有關，其次也與自研程序代碼仍有待優化有關。通過圖6 可知，二者對加筋板以及子彈動響應的模擬結果基本一致。此時子彈已穿透外面板及加筋面板，加筋軸板受到子彈的剮蹭也發生了一定的塑性變形。此外，使用本文方法模擬的臨近失效位置，材料的塑性應變要稍高于LS-DYNA 的結果。圖7a 中比較了兩者模擬所得的子彈運動速度變化歷史，可見其捕捉的子彈降速過程基本一致。顯式計算的準確度可以通過能量平衡進行考察。圖7b 輸出了系統的總能量，其由系統的動能與彈塑性變形能之和得到。由圖7 可知，本文方法所得的總能量更接近于系統的初始總能量，由于在該過程中無外力功，故使用本文方法的計算過程更加接近能量平衡過程。其中，未平衡的部分主要由沙漏控制做功、人工黏性能及單元刪除處理等引起。這也說明了使用本文方法所計算的應力與變形具有更高的可信度。

圖6 對250 μs 時刻子彈及加筋板塑性變形的模擬結果Fig.6 Numerical results for plastic deformation of bullet and reinforced plate at 250 μs:a) bullet;b) reinforced plate

圖7 彈速及總能量比較Fig.7 Comparison of bullet velocity and total energy:a) bullet velocity;b) total energy

3 結論

1）構型一致顯式計算方法基于任意參考構型理論，建立參考已知構型的動量方程拉格朗日有限元格式，更準確地描述與計算結構大變形。通過結合前增量位移的時間積分方法和接觸算法，避免了網格間侵入現象。

2）通過對大變形碰撞及侵徹問題的數值計算，并與LS-DYNA 軟件結果對比，驗證了構型一致顯式計算方法的可靠性。

3）驗證算例同時顯示，該算法可有效抑制大變形碰撞過程的單元畸變，提高計算程序的健壯性和對動態響應預測的準確度，使其在模擬大變形碰撞問題方面具有優勢。