劉赫崴，劉 昆，王秀飛，王正耀

(江蘇科技大學 船舶與海洋工程學院，江蘇 鎮(zhèn)江 212003)

隨著計算機硬件及有限元算法的不斷發(fā)展，有限元數值仿真方法已成為評估工程結構安全性的重要手段。對于船舶碰撞擱淺過程中的結構破壞問題，選擇合適的失效準則是保證仿真計算結果合理準確的關鍵。

國內外學者開展了相關研究工作，其中使用較多的臨界等效塑性應變準則認為當等效塑性應變達到臨界值時，材料出現失效破壞。該準則形式簡單，輸入參數少，但在應力狀態(tài)對失效的影響方面考慮不足。Alsos等[1]結合Hill[2]的局部頸縮分析準則和Bressan等[3]的剪應力準則而提出Bressan & Wlliams-Hill(BWH)失效準則，并將局部頸縮的出現作為斷裂產生的判定依據，避免了局部頸縮區(qū)內過大變形，因而具有網格尺度敏感度低的優(yōu)點，但僅適用于平面應力單元且忽略了彎曲對失效的影響。Johnson等[4]提出的通用斷裂應變模型同時考慮應力三軸度、應變率和溫度效應對失效的影響，具有物理意義清晰、參數易標定等優(yōu)點，但其不能準確描述較高應變率下的斷裂行為。T?rnqvist[5]結合Rice-Tracey準則[6]和Cockcroft-Latham[7]準則提出RTCL準則，相較于上述準則，其定義較為簡單，僅需要失效應變εf作為唯一參數，且考慮全局應力三軸度，對應力狀態(tài)復雜的船體結構失效問題[8-9]模擬效果較好。然而相關研究[10-11]表明，在使用RTCL準則進行斷裂失效模擬時，不同的網格尺寸下仿真結果存在較大差異，需要改變失效參數εf的設定，才能較好地模擬網格尺寸差異較大的船體碰撞仿真。

綜上，本文考慮網格單元尺寸的影響，對現有RTCL失效準則進行修正，建立了修正的RTCL失效準則(I-RTCL)，并通過二次開發(fā)嵌入大型非線性有限元程序中。在此基礎上，充分考慮船舶與海洋工程碰撞問題特點，設計開展了不同應力狀態(tài)下的單軸拉伸試驗與仿真分析，驗證了提出的I-RTCL失效準則在不同應力狀態(tài)、不同網格尺寸下的適用性。

1 修正的RTCL失效準則

1.1 RTCL失效準則

一般認為，以靜水應力與等效應力的比值構造出的無量綱量參數應力三軸度η可以較好地表征塑性變形過程中應力狀態(tài)對韌性斷裂的影響，并經常被用于失效準則的構建[12]。Bao等[13]發(fā)現在負應力三軸度條件下，斷裂由剪切機制掌控；在高應力三軸度范圍內(η≥1/3))，斷裂由孔長大機制所主導；而在介于二者之間的低應力三軸度條件下，斷裂由兩種機制共同作用產生。

(1)

此外，Cockcroft-Latham準則能夠準確預測韌性剪切斷裂(η≤1/3)，Rice-Tracey準則可以正確模擬孔洞生長導致的斷裂(η≥1/3)，T?rnqvist將二者組合并結合Bao等[14]有關斷裂截止值的研究(η≤-1/3時不再發(fā)生斷裂)，得到適用于全局應力三軸度的RTCL準則表示為

(2)

其中，

(3)

式中：D為累計損傷，當一單元在其厚度方向上的所有積分點都滿足D≥1時，該單元刪除；εf為失效應變，是RTCL準則唯一輸入參數；f(η)為關于應力三軸度的函數，在不同應力三軸度區(qū)間內累計損傷函數不同。

1.2 RTCL失效準則修正

由式(2)可知，RTCL準則定義的累計損傷與失效應變密切相關，而在有限元軟件中網格尺寸L又影響失效應變的設定，因此為降低RTCL準則在仿真中因L產生的誤差，增設L為影響失效的自變量并以此構建I-RTCL準則。參照文獻[15]制造平板試件與光滑圓棒試件，開展單軸拉伸試驗及仿真，按表1建立不同網格尺寸、類型的有限元模型，尺寸、試件及部分有限元模型如圖1所示。

表1 有限元模型網格尺寸與種類Tab.1 Mesh sizes and types of FEMs

圖1 尺寸、試件及有限元模型(mm)Fig.1 Dimensions, specimens and FEMs(mm)

以平板試件5 mm模型及光滑圓棒1 mm模型為例，多次嘗試設定失效應變直至仿真與試驗的載荷-位移曲線吻合，仿真失效前的應力云圖與試驗同樣出現頸縮現象，如圖2所示，此時得到對應網格尺寸下的失效應變數值分別為0.62與1.05，同理獲得多組數據并將結果匯總如表2所示。

圖2 試驗與仿真載荷-位移曲線Fig.2 Load-displacement curve in experiment and simulation

表2 不同網格尺寸下有限元模型對應失效應變Tab.2 Failure strain of FEMs with different mesh sizes

(4)

圖3 失效應變與網格尺寸關系Fig.3 Relationship between failure strain and element size

對于平板試件，a=0.717 2，b=-1.219，c=0.534 7；對于光滑圓棒，a=0.738 7，b=-0.506，c=0.342 2。

(5)

(6)

此時，I-RTCL失效準則的累計損傷表達式由式(2)轉變?yōu)?/p>

(7)

1.3 I-RTCL失效準則的VUMAT開發(fā)

利用Fortran語言編寫用戶子程序VUMAT以自定義I-RTCL失效準則，圖4展示了本文編寫子程序判定單元是否失效的計算流程，其在每個增量步內完成如下工作：

(1) 編寫的VUMAT子程序可以從ABAQUS主程序中獲得材料參數、狀態(tài)變量等參數信息用于材料的失效判定。在獲得參數后，首先假定應變增量為彈性應變，使用廣義胡克定律計算試探應力后計算其等效應力，將等效應力與屈服應力對比判斷此增量步下該積分點是否處于屈服狀態(tài)。

(2) 如等效應力小于屈服應力則未屈服，積分點不產生累積損傷，所得試探應力即為此增量步結束時的應力，之后更新能量和狀態(tài)變量，子程序結束，進入下一增量步；如等效應力大于屈服應力則發(fā)生屈服，表示材料產生塑性變形，此時需要將彈性試探應力返回到屈服面上來計算出真實的應力。由所得真實應力計算應力三軸度，根據此時的應力狀態(tài)利用I-RTCL準則計算累計損傷D，判斷積分點是否失效。

(3) 如D<1則未發(fā)生失效，子程序更新能量和狀態(tài)變量后，進入下一增量步；如D≥1則發(fā)生失效，失效積分點刪除，當單元厚度方向上所有積分點均被刪除時，控制單元失效的狀態(tài)變量賦值為0，單元刪除，進入下一增量步。

圖4 VUMAT子程序計算流程圖Fig.4 VUMAT subroutine calculation flow chart

2 I-RTCL失效準則試驗驗證及分析

2.1 拉伸試件設計

為模擬不同應力狀態(tài)下的金屬失效，按圖5尺寸制作多種單軸拉伸試件，按表3劃分為不同網格尺寸與類型的有限元模型，試件與部分有限元模型如圖6所示，包含平板剪切試件：SST-1，SST-2；缺口平板試件：PNT-Rn，n=6，n=10，n=14(n代表缺口半徑，下同)；缺口圓棒試件：SNB-Rn，n=6，n=9，n=18。

由圖5可知，SST-1在被拉伸時連接段(即試件最先失效處)的受力狀態(tài)接近于純剪切，而SST-2的連接段相較于SST-1存在2.5 mm的錯位，因此連接段除受剪力外還會受到拉力。通過控制PNT-Rn與SNB-Rn的缺口半徑來控制連接段受力狀態(tài)，總體上，隨著缺口半徑n的增大，應力三軸度η不斷增大。

表3 有限元模型網格尺寸與類型Tab.3 Mesh sizes and types of FEMs

圖5 拉伸試件尺寸(mm)Fig.5 Dimension of the specimens (mm)

圖6 拉伸試件及有限元模型(mm)Fig.6 Tensile specimen and FEMs (mm)

2.2 拉伸試驗與仿真概況

對相同網格尺寸的模型分別使用I-RTCL準則與RTCL準則開展拉伸仿真，對比修正前后仿真斷裂時刻拉伸位移與試驗值的誤差以驗證I-RTCL優(yōu)勢?？紤]到平板與圓棒試件結構上差異較大，在仿真時分別選取用式(5)與式(6)修正的失效準則。

2.3 結果分析

2.3.1 失效單元應力狀態(tài)

以NRB-Rn為例，其失效單元應力三軸度變化曲線如圖7所示，在加載初期曲線存在振蕩情況，隨后曲線穩(wěn)定在均值附近，失效單元應力狀態(tài)較為穩(wěn)定，說明利用單軸拉伸試件可有效模擬某一應力狀態(tài)下的金屬斷裂失效行為。

圖7 失效單元應力三軸度變化曲線Fig.7 Stress triaxiality of failure element

表4 不同試件失效單元平均應力三軸度Tab.4 Mean stress traxiality of failure element in different specimens

2.3.2 載荷-位移曲線

選取SST-1，SST-2，PNT-R6與NRB-R6拉伸試驗與仿真的載荷-位移曲線如圖8所示，可以看到：

(1) 修正前后的仿真與試驗測量載荷-位移曲線吻合均較好，但使用I-RTCL失效準則的仿真斷裂時刻拉伸位移分布更為集中，說明在多種應力狀態(tài)、多種網格尺寸下，相比原RTCL準則本文提出的I-RTCL失效準則精度更高。

圖8 試驗與仿真對應載荷-位移曲線Fig.8 Force-displacement relation in test and simulation

(2) 純剪切試件SST-1的拉伸曲線如圖8(a)、圖8(b)所示，其仿真曲線后半段載荷明顯高于試驗結果，結合0.4 mm工況下失效單元應力三軸度隨拉伸產生的變化(圖9)不難看到：在仿真后期失效單元的三軸度仍不斷上升，單元呈現拉伸狀態(tài)，因此載荷仍呈上升趨勢。由于載荷隨拉伸位移增大持續(xù)上升，極限強度也表現出與斷裂位移相似的尺寸效應，斷裂拉伸位移更集中的I-RTCL準則對應仿真的強度極限也更加接近試驗值。

(3) 表5為RTCL準則與I-RTCL準則在不同網格尺寸下仿真斷裂位移及其誤差，可以看到對于中低應力三軸度狀態(tài)(SST-1，SST-2，PNT-Rn)，修正后仿真的誤差明顯降低；而對于高應力狀態(tài)，存在修正效果一般的情況(NRB-R9-0.5 mm，NRB-R18-1 mm)；若修正前誤差較小，經過修正誤差仍保持在較低水平(SST-1-0.8 mm，PNT-R6-2 mm，NRB-R9-1 mm)。

(4) 對于缺口試件，缺口半徑越小，斷裂發(fā)生越快，分析產生此現象的原因主要如下：一是缺口半徑越小，試件失效區(qū)域應力三軸度越高，越容易發(fā)生斷裂；二是缺口半徑越大，試件發(fā)生變形的連接段越長，其產生相同應變需要拉伸的位移也越長。

表5 仿真斷裂位移及誤差Tab.5 Fracture displacements and errors in simulations

表5 (續(xù))

圖9 0.4 mm計算工況失效單元應力三軸度Fig.9 Stress triaxiality of failure element in 0.4 mm calculation condition

3 結 論

(1) 相同模型，采用不同的網格尺寸進行仿真計算得到的斷裂位移相差較大，網格尺寸對失效分析準確性影響較大，失效應變隨網格尺寸的增大而降低，并逐漸趨于定值。

(3) 使用I-RTCL準則進行多應力狀態(tài)、多網格尺寸下的仿真，其計算精度高于原準則，其中對于中低應力三軸度，提高顯著；對于高應力三軸度，精度有所提升但存在修正效果一般的情況；此外，若修正前誤差較小，經過修正誤差仍保持在較低水平。