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

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

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

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

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

1 修正的RTCL失效準(zhǔn)則

1.1 RTCL失效準(zhǔn)則

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

(1)

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

(2)

其中，

(3)

式中：D為累計(jì)損傷，當(dāng)一單元在其厚度方向上的所有積分點(diǎn)都滿足D≥1時(shí)，該單元?jiǎng)h除；εf為失效應(yīng)變，是RTCL準(zhǔn)則唯一輸入?yún)?shù)；f(η)為關(guān)于應(yīng)力三軸度的函數(shù)，在不同應(yīng)力三軸度區(qū)間內(nèi)累計(jì)損傷函數(shù)不同。

1.2 RTCL失效準(zhǔn)則修正

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

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

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

以平板試件5 mm模型及光滑圓棒1 mm模型為例，多次嘗試設(shè)定失效應(yīng)變直至仿真與試驗(yàn)的載荷-位移曲線吻合，仿真失效前的應(yīng)力云圖與試驗(yàn)同樣出現(xiàn)頸縮現(xiàn)象，如圖2所示，此時(shí)得到對(duì)應(yīng)網(wǎng)格尺寸下的失效應(yīng)變數(shù)值分別為0.62與1.05，同理獲得多組數(shù)據(jù)并將結(jié)果匯總?cè)绫?所示。

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

表2 不同網(wǎng)格尺寸下有限元模型對(duì)應(yīng)失效應(yīng)變Tab.2 Failure strain of FEMs with different mesh sizes

(4)

圖3 失效應(yīng)變與網(wǎng)格尺寸關(guān)系Fig.3 Relationship between failure strain and element size

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

(5)

(6)

此時(shí)，I-RTCL失效準(zhǔn)則的累計(jì)損傷表達(dá)式由式(2)轉(zhuǎn)變?yōu)?/p>

(7)

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

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

(1) 編寫的VUMAT子程序可以從ABAQUS主程序中獲得材料參數(shù)、狀態(tài)變量等參數(shù)信息用于材料的失效判定。在獲得參數(shù)后，首先假定應(yīng)變?cè)隽繛閺椥詰?yīng)變，使用廣義胡克定律計(jì)算試探應(yīng)力后計(jì)算其等效應(yīng)力，將等效應(yīng)力與屈服應(yīng)力對(duì)比判斷此增量步下該積分點(diǎn)是否處于屈服狀態(tài)。

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

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

圖4 VUMAT子程序計(jì)算流程圖Fig.4 VUMAT subroutine calculation flow chart

2 I-RTCL失效準(zhǔn)則試驗(yàn)驗(yàn)證及分析

2.1 拉伸試件設(shè)計(jì)

為模擬不同應(yīng)力狀態(tài)下的金屬失效，按圖5尺寸制作多種單軸拉伸試件，按表3劃分為不同網(wǎng)格尺寸與類型的有限元模型，試件與部分有限元模型如圖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在被拉伸時(shí)連接段(即試件最先失效處)的受力狀態(tài)接近于純剪切，而SST-2的連接段相較于SST-1存在2.5 mm的錯(cuò)位，因此連接段除受剪力外還會(huì)受到拉力。通過控制PNT-Rn與SNB-Rn的缺口半徑來控制連接段受力狀態(tài)，總體上，隨著缺口半徑n的增大，應(yīng)力三軸度η不斷增大。

表3 有限元模型網(wǎng)格尺寸與類型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 拉伸試驗(yàn)與仿真概況

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

2.3 結(jié)果分析

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

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

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

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

2.3.2 載荷-位移曲線

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

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

圖8 試驗(yàn)與仿真對(duì)應(yīng)載荷-位移曲線Fig.8 Force-displacement relation in test and simulation

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

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

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

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

表5 (續(xù))

圖9 0.4 mm計(jì)算工況失效單元應(yīng)力三軸度Fig.9 Stress triaxiality of failure element in 0.4 mm calculation condition

3 結(jié) 論

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

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