移除單元法中單元失效應變的確定

張美蘭，陳 震，趙 晟

上海交通大學船舶海洋與建筑工程學院，上海200240

0 引 言

近年來，數值計算方法被廣泛應用于船舶碰撞、擱淺、穿甲和爆炸損傷等工程結構分析中，較好地模擬了船體結構在復雜外力條件下的破壞問題［1-4］。其中，基于失效應變的移除單元法就是模擬大型結構破壞過程的重要方法之一，該方法通過定義單元失效應變來描述局部結構因變形過大而逐漸失去承載能力并最終導致破壞的行為。利用移除單元法，可以準確模擬船體結構的失效路徑和模式、觀察結構的變形破壞特征和獲得結構的極限承載能力，為優化結構設計提供可靠依據。

移除單元法的核心是單元失效應變的確定，合理地選取材料失效應變是獲得可靠結果的基礎。本文將以裂紋板拉伸斷裂試驗為基礎，研究單元失效應變的選取方法，分析網格尺寸對失效應變的影響，并討論誤差范圍，以為該方法的進一步應用提供依據。

1 問題的提出

當結構承受過大載荷時，材料會產生塑性應變，同時微觀組織發生分離、開裂，并最終發展成宏觀裂紋。該過程經歷了從微觀到宏觀的發展，因此，從材料的微觀特性出發，同時考慮宏觀的力學行為，建立材料破壞斷裂準則，是研究結構破壞問題的理想做法。然而，盡管國內外學者已對此進行過廣泛而深入的研究，但目前還是難以建立統一的破壞準則以應用于大型結構物的分析中［5-7］。在開展船舶與海洋結構物的破損強度分析時，通常是采用簡化的移除單元法模擬局部結構的破壞失效行為，從宏觀上關注裂紋的力學行為與結構整體響應的對應，斷裂準則不直接用于描述材料的微觀力學。當定義平均等效塑性應變達到某臨界值時，該單元會因失去承載能力而破壞，該臨界值稱為單元失效應變。

單元失效應變與材料拉伸斷裂應變不同，前者是為了模擬結構破壞失效而引入的單元失效判據，后者是材料自身的力學特性參量。失效單元的幾何尺寸通常比拉伸斷裂處的尺寸大很多。單元失效應變不能簡單地等于材料拉伸斷裂應變值；等效塑性應變反映的是單元面積范圍內的平均應變水平，單元尺寸不同，失效應變也有所差異。Simonsen 等［8］通過分析比較板拉伸試驗與數值計算結果，指出失效應變與單元尺寸之間關系密切，在應用失效應變作為單元失效準則時應考慮網格大小。Lehmann 等［9］在分析試驗結果的基礎上，認為塑性變形發生在一個有限的局部帶狀寬度內，當帶狀區域寬度達到某延展率時，板會發生斷裂。由于單元尺寸大于局部帶狀寬度，所以基于單元平均失效應變的斷裂準則依賴于網格尺寸大小。并且，通過測量裂紋周圍區域的厚度，可以得到失效應變與單元尺寸之間的關系。

單元失效應變的合理選取是采用移除單元法進行船體結構破壞損傷數值仿真的基礎。本文將通過開展系列裂紋板拉伸斷裂試驗，獲得不同裂紋尺寸下表征裂紋板宏觀力學特性的位移載荷曲線。同時，采用有限元軟件MSC.Dytran 進行數值仿真計算，分別以載荷極值、最大伸長量和能量為依據得到與網格尺寸相關的單元失效應變。

2 裂紋板拉伸斷裂試驗

為了獲得含破損結構的板的宏觀力學特性，本文開展了裂紋板拉伸斷裂試驗，得到了裂紋板的整體位移載荷曲線、載荷極值和最大伸長量等參數。試驗在某大學的工程力學實驗中心進行，共包含3 種裂紋長度的裂紋板，且均為中心貫穿型裂紋，尺寸如圖1 所示。試件材料為Q235 普通冷軋鋼，在拉伸試驗機上測得材料的力學參數為：彈性模量E = 157.6 GPa，極限拉伸強度σT= 299.5 MPa，屈服應力σy= 172 MPa，斷裂工程應變εf=0.33。圖2 所示為裂紋板試件的應力—應變曲線。

圖1 裂紋板試件尺寸圖Fig.1 The dimension of cracked plate specimen

圖2 裂紋板試件工程應力—應變曲線Fig.2 Normal stress-strain curve of the plate material

試驗加載裝置為MTS-880 材料疲勞試驗機，采用位移控制方式對試件的一端施加均布位移，同步輸出位移載荷數值。為保證試驗處于準靜態條件下，位移加載速率恒定為0.2 mm/min。在試驗過程中，隨著夾持端位移的增大，裂紋板沿對稱中心對稱變形，隨后，裂紋尖端出現明顯的撕裂狀，并沿水平方向穩定擴展，直至斷裂。圖3 所示為SPVT_040C 裂紋板試驗時裂紋開始擴展和即將斷裂的情況。

圖3 SPVT_040C 板裂紋擴展及斷裂圖Fig.3 The pictures of crack propagation and fracture of the SPVT_040C plate

表1 所示為裂紋板試件尺寸及由拉伸試驗所得的載荷極值和最大伸長量。圖4 所示為裂紋板的位移載荷曲線，由圖可看出，兩次由SPVT_030C裂紋板拉伸試驗所得的位移載荷曲線形式、載荷極值和最大伸長量數據基本一致。

表1 裂紋板試件尺寸及試驗結果Tab.1 Dimension and experimental results of cracked plates

圖4 裂紋板試驗位移載荷曲線Fig.4 Experimental displacement-load curves of cracked plates

3 數值模擬

本文采用MSC.Dytran 模擬裂紋板的拉伸斷裂過程，模型采用四邊形板單元建立，尺寸與試驗試件相同，在模型的中心位置，刪除了與裂紋長度相同的一組單元以模擬初始裂紋。網格尺度用裂紋尖端單元邊長s 與板厚t 的比值表示，取s/t =0.625～10。在模型一端施加恒定速度的位移，同時輸出另一端的約束反力。圖5 所示為網格尺度s/t =1.25 時的SPVT_040C 裂紋板計算模型及其裂紋開始擴展和即將斷裂的位移云圖。

圖5 SPVT_040C 裂紋板模型和位移云圖（s/t =1.25）Fig.5 SPVT_040C model and displacement contours（s/t =1.25）

計算模型材料參數由裂紋板試件試驗結果得到，采用硬化指數形式表示，即其中，σT為真實應力，εT為真實應變，C = 565 MPa，n=0.28。圖6 所示為擬合所得模型材料的真實應力—應變曲線。

圖6 模型材料的真實應力—應變曲線Fig.6 True stress-strain curve of material for simulation

材料失效模型采用最大塑性應變失效，單元等效塑性應變失效準則采用V-M 失效準則

式中：ε1，ε2，ε3分別為3 向主應變；εep為單元等效塑性應變；εf為單元失效應變。

4 單元失效應變的確定

針對各網格尺寸計算模型，分別定義不同的單元失效應變，計算得出該失效應變下的位移載荷曲線。圖7 所示為裂紋長度為20 mm，網格尺度s/t = 1.25 下各失效應變的計算結果。由圖可見，試驗與有限元計算結果均包含有載荷快速增加段、緩慢上升至極值段以及卸載段。單元失效應變不同，卸載段出現的位置也明顯不同。單元失效應變越小，初始裂紋尖端處的單元便會越早達到失效應變值而失去承載能力，裂紋發生擴展。模型若整體較早出現了卸載，載荷極值和裂紋板完全斷裂時的伸長量便也會較小。

圖7 不同的εf 計算的位移載荷曲線（s/t =1.25）Fig.7 Simulated displacement-load curves by different failure strains（s/t =1.25）

根據裂紋板拉伸斷裂位移載荷曲線的這些特點，將以裂紋長度20 mm 為例，分別以載荷極值、最大伸長量和能量為準則研究單元失效應變的選取問題。

4.1 載荷極值準則

以裂紋板拉伸斷裂過程的載荷極值為依據，結合有限元計算和試驗值確定單元失效應變值。圖8 所示為各網格尺寸下由不同失效應變計算所得的載荷極值繪制的載荷—應變曲線。分別用試驗載荷極值與該曲線相交，可得出對應網格尺寸下以載荷極值為準則的單元失效應變（εf_L）。

圖8 載荷極值準則確定εf_LFig.8 Failure strain on maximum load criterion

圖9 所示為不同網格尺寸下，依據載荷極值準則的失效應變計算所得的位移載荷曲線。由圖可見，盡管計算所得的載荷極值大小與試驗接近，但載荷極值出現的位置和裂紋板斷裂時的最大伸長量均小于試驗值，位移載荷曲線整體上與試驗偏差較大，并且偏差是隨網格尺寸的增加而增大。例如，當網格尺度s/t = 0.625 時，計算所得的最大伸長量約為試驗值的70%，而當s/t = 5 時則約為50%。用以載荷極值為準則的失效應變計算裂紋板的拉伸斷裂難以全面反映其力學特性。

圖9 εf_L 計算的位移載荷曲線Fig.9 Simulated displacement-load curves by εf_L

4.2 最大伸長量準則

以裂紋板斷裂時的最大伸長量為依據，用與載荷極值準則失效應變相同的方法確定不同網格尺寸下的單元失效應變（εf_D），如圖10 所示。

圖11 所示為不同網格尺寸下，依據最大伸長量準則確定的失效應變計算所得的位移載荷曲線。由圖可見，在相同網格尺寸下，以最大伸長量為準則的單元失效應變均大于以載荷極值為準則的失效應變，載荷極值的出現位置與試驗值較接近，其大小略高于試驗結果。與載荷極值準則失效應變相比，根據最大伸長量準則的失效應變計算所得的載荷極值的誤差較小，當網格尺度s/t =0.625 時，載荷極值誤差約為6.31%；當s/t = 5 時，誤差約為16.67%。

圖10 最大伸長量準則確定εf_DFig.10 Failure strain on maximum displacement criterion

圖11 εf_D 計算的位移載荷曲線Fig.11 Simulated displacement-load curves by εf_D

載荷極值的誤差源于移除單元法無法細致描述裂紋尖端復雜的彈塑性應力分布情況，而只能用單元平均的應力—應變反映該單元整體尺寸范圍內的應力—應變狀態。并且，由材料本構關系曲線可知，材料進入塑性段后，隨著塑性應變的增加，應力變化極度平緩。單元失效應變的改變對結構承載能力的影響不大，而對伸長量的影響則更為明顯。

4.3 能量準則

在裂紋板拉伸斷裂過程中，外力持續做功，所做功的大小即為位移載荷曲線下所圍的面積，是裂紋板整體抵抗外力能力的表征，如圖12 所示。

式中：P 為裂紋板拉伸載荷；δ 為裂紋板拉伸位移；E 為外力功。

圖12 外力做功圖Fig.12 Work done by external force

以外力功為依據，結合有限元計算和試驗值確定單元失效應變。根據式（2），將各網格尺寸下不同失效應變計算所得的位移載荷曲線進行積分，得出外力所作的功繪制的能量—失效應變曲線。用試驗外力功（P020C= 670.335 kN·mm）與該曲線相交，得出對應網格尺寸下能量準則的單元失效應變（εf_E），如圖13 所示。

圖13 能量準則確定εf_EFig.13 Failure strain on maximum energy criterion

圖14 所示為不同網格尺寸下，依據能量準則的失效應變計算所得的位移載荷曲線。由圖可見，載荷極值和最大伸長量與試驗值均有一定的偏差，載荷極值略高于試驗值，最大伸長量略小于試驗值。當網格尺度s/t = 0.625 時，載荷極值誤差約為6.31%，最大伸長量約為7.96%；當s/t = 5時，載荷極值誤差約為14.75%，最大伸長量約為18.02%。該準則綜合考慮了裂紋板拉伸過程中的載荷與位移因素，反映了裂紋板整體承受外載荷的能力，位移載荷曲線能較好地反映裂紋板的力學特性。

4.4 不同準則失效應變誤差分析

針對裂紋長度分別為30 mm 和40 mm 的裂紋板拉伸試驗，分別開展以上工作，得出了不同裂紋長度下的單元失效應變。圖15 所示為以最大伸長量為準則的3 種裂紋長度下的單元失效應變及其平均值曲線。由圖可見，在不同裂紋長度情況下，各網格尺寸下的單元失效應變很接近，取其平均值可以作為該網格尺寸的失效應變。圖16 所示為分別以載荷極值、最大伸長量和能量準則確定的3 種裂紋長度失效應變平均值與網格尺寸的關系曲線。

圖14 εf_E 計算得位移載荷曲線Fig.14 Simulated displacement-load curves by εf_E

圖15 最大伸長量準則的失效應變曲線Fig.15 Failure strain curves on maximum displacement criterion

以裂紋長度20 mm 為例，依據3 種準則的單元失效應變計算得出了載荷極值、最大伸長量和能量，并分別將計算結果與試驗值進行了比較，誤差如表2 所示。

圖16 3 種準則的失效應變曲線Fig.16 Failure strain curves of three criterions

式中：R計算為有限元計算結果；R試驗為試驗值；ω為誤差，%。

由表可見，計算結果與試驗值的誤差是隨網格尺度的增加而增大。由載荷極值準則的失效應變計算所得的最大伸長量和能量均小于試驗值，且誤差較大。以最大伸長量和能量為準則的失效應變計算所得的結果與試驗值較接近，載荷極值均略高于試驗值。

5 算例驗證

本文分別用載荷極值、最大伸長量和能量準則對文獻［10］中的裂紋板拉伸斷裂試驗進行了數值計算，以驗證單元失效應變選取的合理性。裂紋板的長、寬為500 mm，厚度1.6 mm，初始裂紋長度分別為60，30 和15 mm。材料參數為：彈性模量E = 198.3 GPa，屈服應力σy= 296.1 MPa，極限拉伸強度σT= 362.1 MPa，工程斷裂應變εf= 0.369。根據以上參數，建立有限元計算模型。因文獻［10］中試驗材料的斷裂應變與本文的試樣材料很接近，故以圖16 的平均網格尺度—失效應變曲線得出的不同網格尺度所對應的失效應變作為數值計算中的單元失效應變。

表2 SPVT_020C 板計算結果誤差Tab.2 The deviation of simulated results（SPVT_020C）

圖17 所 示為NP16-60，NP16-30 和NP16-15裂紋板計算位移載荷曲線與試驗結果的比較。由圖可見，計算載荷極值均略高于試驗結果，且用載荷極值準則的失效應變計算的位移載荷曲線與試驗偏差較大，當NP16-30 板的網格尺度s/t = 6.25時，最大伸長量偏差達44.06%。能量和最大伸長量準則的失效應變計算結果與試驗值較接近，例如，NP16-15 板的網格尺度s/t = 3.125 時，計算載荷極值以最大伸長量為準則所得的誤差約為6.89%，以能量為準則的誤差為6.41%。

圖17 3 種裂紋板計算位移載荷曲線Fig.17 Simulated displacement-load curves of three kinds of cracked plates

6 結 論

本文通過移除單元法數值模擬了裂紋板的拉伸斷裂過程，確定了3 種準則的單元失效應變，并分析討論了數值計算結果，結論如下：

1）采用能量準則確定單元失效應變時，綜合考慮了裂紋板拉伸過程中的載荷和位移因素，且計算所得的載荷極值和最大伸長量與試驗值間的偏差均較小，能較好地反映裂紋板的力學特性。

2）依據最大伸長量準則的失效應變計算所得的位移載荷曲線與試驗值吻合較好，載荷極值略高于試驗結果，適用于優先考慮裂紋板伸長位移的數值計算。而以載荷極值為準則的失效應變計算結果則偏差較大，難以全面反映裂紋板拉伸斷裂的力學特性。

3）移除單元法無法細致地描述裂紋尖端復雜的彈塑性應力分布情況，導致計算結果與試驗值有一定的誤差，但對模擬大型結構物的破壞過程仍具適用性。

［1］姚熊亮，楊樹濤，張阿漫. 爆炸載荷作用下艦船板架的變形與斷裂研究綜述［J］. 中國艦船研究，2009，4（1）：1-7，12.YAO Xiongliang，YANG Shutao，ZHANG Aman. Re?view for deformation and fracture behavior of ship stiff?ened plate subjected to blast loading［J］.Chinese Jour?nal of Ship Research，2009，4（1）：1-7，12.

［2］王善，蓋京波，楊世全. 艦船結構在爆炸載荷作用下的破壞研究綜述［J］. 艦船科學技術，2006，28（3）：8-11.WANG Shan，GAI Jingbo，YANG Shiquan. Study sum?mary of ship structure under blast load［J］. Ship Sci?ence and Technology，2006，28（3）：8-11.

［3］祈恩榮，崔維成.船舶碰撞和擱淺研究綜述［J］.船舶力學，2001，5（4）：67-80.QI Enrong，CUI Weicheng. A state-of-the-art review on ship collision and grounding［J］.Journal of Ship Me?chanics，2001，5（4）：67-80.

［4］梅志遠，朱錫，劉潤泉. 船用加筋板架爆炸載荷下動態響應數值分析［J］. 爆炸與沖擊，2004，24（1）：80-84.MEI Zhiyuan，ZHU Xi，LIU Runquan. Dynamic re?sponse researches of ship’s stiffened plate structure un?der explosive load［J］. Explosion and Shock Waves，2004，24（1）：80-84.

［5］趙效東，張樂山.海洋工程結構物碰撞失效準則探討［J］.船海工程，2012，41（2）：144-148.ZHAO Xiaodong，ZHANG Leshan. Research of failure criterion of offshore structural collision［J］. Ship and Ocean Engineering，2012，41（2）：144-148.

［6］陳永念，譚家華.數值仿真中單元密度對材料失效應變的影響［J］.船海工程，2007，36（6）：1-4.CHEN Yongnian，TAN Jiahua. Mesh sensitivity of ele?ment failure strain in numerical simulation［J］. Ship and Ocean Engineering，2007，36（6）：1-4.

［7］SIMONSEN B C，T?RNQVIST R. Experimental and numerical modelling of ductile crack propagation in large-scale shell structures［J］. Marine Structures，2004，17（1）：1-27.

［8］SIMONSEN B C，LAURIDSEN L P. Energy absorption and ductile failure in metal sheets under lateral inden?tation by a sphere［J］. International Journal of Impact Engineering，2000，24（10）：1017-1039.

［9］LEHMANN E，PESCHMANN J. Energy absorption by the steel structure of ships in the event of collisions［J］.Marine Structures，2002，15（4/5）：429-441.

［10］PAIK J K，THAYAMBALLI A K. Ultimate strength of ageing ships［J］.Journal of Engineering for the Mar?itime Environment，2002，216（1）：57-77.