應(yīng)用構(gòu)型一致的顯式方法預(yù)測(cè)結(jié)構(gòu)沖擊動(dòng)響應(yīng)

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

（1.中國(guó)船舶科學(xué)研究中心，江蘇 無錫 214082；2.深海技術(shù)科學(xué)太湖實(shí)驗(yàn)室，江蘇 無錫 214082）

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

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

1 理論方程及數(shù)值算法

1.1 控制方程

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

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

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

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

1.2 構(gòu)型一致的方程弱形式

動(dòng)量方程的積分弱形式為：

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

1.3 本構(gòu)模型

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

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

材料失效采用考慮應(yīng)力三軸度、應(yīng)變率效應(yīng)及熱軟化的Johnson-Cook 準(zhǔn)則，臨界失效應(yīng)變?yōu)閇18]：

1.4 顯式算法

對(duì)該方程求解時(shí)需滿足tn+2時(shí)刻接觸約束條件，即：

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

2 數(shù)值算例

2.1 Taylor 桿撞擊

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

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

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

圖2 撞擊桿初速度400 m/s、對(duì)桿在100 μs 時(shí)的變形狀態(tài)模擬結(jié)果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

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

圖4 撞擊桿初速度800m/s、通過本文方法獲得的80 μs 時(shí)刻塑性應(yīng)變分布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，對(duì)不同時(shí)刻局部變形的模擬結(jié)果Fig.5 Numerical results for local deformation at different moments in the case of initial velocity being 800 m/s

2.2 高速侵徹加筋板

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

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

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

圖6 對(duì)250 μs 時(shí)刻子彈及加筋板塑性變形的模擬結(jié)果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 結(jié)論

1）構(gòu)型一致顯式計(jì)算方法基于任意參考構(gòu)型理論，建立參考已知構(gòu)型的動(dòng)量方程拉格朗日有限元格式，更準(zhǔn)確地描述與計(jì)算結(jié)構(gòu)大變形。通過結(jié)合前增量位移的時(shí)間積分方法和接觸算法，避免了網(wǎng)格間侵入現(xiàn)象。

2）通過對(duì)大變形碰撞及侵徹問題的數(shù)值計(jì)算，并與LS-DYNA 軟件結(jié)果對(duì)比，驗(yàn)證了構(gòu)型一致顯式計(jì)算方法的可靠性。

3）驗(yàn)證算例同時(shí)顯示，該算法可有效抑制大變形碰撞過程的單元畸變，提高計(jì)算程序的健壯性和對(duì)動(dòng)態(tài)響應(yīng)預(yù)測(cè)的準(zhǔn)確度，使其在模擬大變形碰撞問題方面具有優(yōu)勢(shì)。