考慮殘余應(yīng)力的加筋板結(jié)構(gòu)疲勞裂紋應(yīng)力強度因子計算方法研究

王藝陶，馮國慶，任慧龍

（1.中國船舶科學(xué)研究中心，江蘇 無錫 214082；2.哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150010）

0 引 言

加筋板結(jié)構(gòu)是船舶與海洋工程結(jié)構(gòu)最重要、最基本的組成部分，在日常作業(yè)和運營過程中海浪的交變載荷作用下，容易產(chǎn)生疲勞缺陷進而導(dǎo)致疲勞破壞。同時，船體結(jié)構(gòu)在焊接過程中，不可避免地存在著相對復(fù)雜的焊接殘余應(yīng)力。近年來，關(guān)于殘余應(yīng)力對船舶結(jié)構(gòu)裂紋擴展影響的研究受到越來越多的關(guān)注[1-4]，關(guān)于殘余應(yīng)力計算方法的討論也不斷深入[5-7]。研究表明，殘余應(yīng)力場的存在會直接影響到結(jié)構(gòu)的疲勞裂紋擴展，使焊接結(jié)構(gòu)的疲勞裂紋擴展行為呈現(xiàn)出與非殘余應(yīng)力場中擴展不同的特性[8]。因此，研究含有初始裂紋加筋板結(jié)構(gòu)的疲勞裂紋擴展過程的數(shù)值計算方法，將有限元數(shù)值仿真方法與理論模型相結(jié)合，對加筋板結(jié)構(gòu)焊縫周圍殘余應(yīng)力分布進行數(shù)值仿真，并給出殘余應(yīng)力場中裂紋尖端應(yīng)力強度因子的計算結(jié)果具有十分重要的意義。

本文以有限元數(shù)值方法為基礎(chǔ)結(jié)合Green函數(shù)應(yīng)力強度因子計算理論，建立一套考慮焊接殘余應(yīng)力影響的加筋板結(jié)構(gòu)裂紋擴展的模擬方法。通過對含中心穿透裂紋平板結(jié)構(gòu)、含對接焊縫及中心穿透裂紋平板結(jié)構(gòu)等一系列結(jié)構(gòu)的討論，給出外部載荷作用下和焊接殘余應(yīng)力作用下有限元數(shù)值計算應(yīng)力強度因子的方法，并驗證其有效性；然后，結(jié)合裂紋擴展Paris公式，建立一套考慮焊接殘余應(yīng)力影響的裂紋擴展模擬流程，為同類結(jié)構(gòu)裂紋擴展規(guī)律研究以及裂紋擴展壽命安全性評估工作提供重要參考和借鑒。

1 含中心裂紋平板結(jié)構(gòu)裂紋應(yīng)力強度因子計算方法

1.1 應(yīng)力強度因子計算理論

采用J方法，對應(yīng)力強度因子進行計算。J積分描述的是由于裂紋的存在所吸收的能量，而應(yīng)變釋放率G描述的則是產(chǎn)生新裂紋面所需要的能量。對于線彈性材料而言，應(yīng)力強度因子K、路徑無關(guān)積分J積分和應(yīng)變能力釋放率G可以通過常數(shù)聯(lián)系起來，并且J積分和應(yīng)變能力釋放率是等價的。

式中：ui為位移矢量的分量；ds為積分路徑Γ上的微小增量；ω為應(yīng)變能力因子。

如圖1（b）所示，通過散度定理，用裂紋尖端附近一個有限區(qū)域來代替積分回路進行J積分的計算?；谶@樣的思路，（1）式可以轉(zhuǎn)化為：

如圖1（a）所示，考慮任意一個圍繞裂紋尖端的逆時針回路Γ，數(shù)學(xué)表達式如下：

針對二維問題可具體展開為：

1.2 應(yīng)力強度因子計算有限元建模方法

含中心穿透裂紋平板結(jié)構(gòu)是研究應(yīng)力強度因子計算過程中最典型也是最基礎(chǔ)的結(jié)構(gòu)之一，由于結(jié)構(gòu)簡單且應(yīng)力強度因子手冊已經(jīng)給出了其應(yīng)力強度因子的解析解，因此，常常作為應(yīng)力強度因子計算的首選目標。在研究基于有限元方法計算應(yīng)力強度因子的過程中，就首先以含中心穿透平板結(jié)構(gòu)為研究對象，建立有限元模型，采用不同尺寸網(wǎng)格對結(jié)構(gòu)進行劃分，并采用J積分方法計算裂紋尖端應(yīng)力強度因子，通過將有限元數(shù)值計算結(jié)果與理論計算結(jié)果的比較，驗證有限元數(shù)值方法的正確性，并討論網(wǎng)格大小對應(yīng)力強度因子計算結(jié)果的影響，同時，給出應(yīng)力強度因子在板厚方向不同位置計算結(jié)果的變化情況，最終，確定合適的有限元網(wǎng)格尺寸及計算點位置[9]。

含中心穿透裂紋的平板長l=2 000 mm，寬b=560 mm，厚度t=8 mm，中心裂紋長2a，平板兩端受均勻拉應(yīng)力σ=134 MPa。模型的材料參數(shù)E=2.1×105MPa，泊松比μ=0.3。

在進行有限元模型建模時，采用分區(qū)域控制網(wǎng)格疏密程度的劃分方法，將模型分為三大部分（如圖2所示），即接近裂紋尖端及裂紋延長線范圍內(nèi)為區(qū)域一，網(wǎng)格密度最大；遠離裂紋及其延長線范圍為區(qū)域三，網(wǎng)格密度最??；中間部分為過渡區(qū)域，這種方法可以兼顧裂紋尖端應(yīng)力強度因子的精度與計算效率。

圖2有限元模型網(wǎng)格劃分區(qū)域Fig.2 Zoning of finite element model

圖3有限元模型邊界條件及加載Fig.3 Boundary constraint and loading

其中：a為中心裂紋半寬；b為平板寬度；σ為平板兩端均布法向拉力。

模型邊界條件和加載情況如圖3所示。由于模型的對稱性，僅對結(jié)構(gòu)的1/2進行建模，在中心裂紋兩側(cè)，采用對稱邊界條件進行約束，平板中心邊緣位置處A、B兩點采用全約束，模型兩端加載均布拉力σ。有限元模型及網(wǎng)格劃分情況，如圖4所示。

對含中心穿透裂紋的平板結(jié)構(gòu)進行有限元數(shù)值計算，如圖5所示為有限元模型計算結(jié)果的應(yīng)力云圖及裂紋尖端附近Y方向應(yīng)力結(jié)果，這里可以清楚顯示出裂紋張開情況及裂紋尖端高應(yīng)力區(qū)的分布范圍。

含中心穿透裂紋平板結(jié)構(gòu)裂紋尖端應(yīng)力強度因子的理論計算式如下式所示：

圖4有限元模型邊界條件及加載 Fig.4 Constraint and loading condition

圖5 Y方向應(yīng)力云圖Fig.5 Y component stress contour plot

在計算過程中，裂紋尖端區(qū)域細化網(wǎng)格尺寸分別取為10 mm，5 mm，2 mm，1 mm，0.5 mm和0.3 mm，其中選取有代表性的10 mm和0.5 mm的細化模型如圖6所示。

圖6裂紋區(qū)網(wǎng)格細化情況Fig.6 Mesh refinement

裂紋尖端附近網(wǎng)格尺度對應(yīng)力強度計算精度有一定的影響。于是，將不同網(wǎng)格尺度有限元模型計算結(jié)果與理論公式計算結(jié)果進行對比，以裂紋長度為a=50 mm，σ=134 MPa的計算結(jié)果為例，如表1所示。由結(jié)果可知，有限元模型計算結(jié)果與理論公式計算結(jié)果基本一致，證明有限元數(shù)值計算結(jié)果是正確有效的。

表1應(yīng)力強度因子有限元結(jié)果與理論結(jié)果對比Tab.1 Comparison of finite element results and theoretical results

同時，隨著網(wǎng)格密度增加，網(wǎng)格尺度的逐漸減小，應(yīng)力強度因子的有限元模型計算結(jié)果逐漸增加，與理論公式計算結(jié)果間的誤差逐漸減小，當網(wǎng)格尺寸為1 mm時，誤差僅為-0.31%，但當網(wǎng)格尺度繼續(xù)減小，誤差又逐漸升高。可見，有限元網(wǎng)格的尺度并不是越小越好，權(quán)衡計算精度與計算效率的關(guān)系，選取裂紋尖端附近網(wǎng)格尺度為0.5 mm或1 mm進行計算，應(yīng)力強度因子計算結(jié)果均是相對準確的。兩種網(wǎng)格尺寸下應(yīng)力強度因子計算結(jié)果與理論值對比情況如圖7所示。

另外，由于平面應(yīng)力/平面應(yīng)變效應(yīng)的不同，同一裂紋長度下應(yīng)力強度因子在厚度方向存在一定的變化，如圖8所示，本文所涉及的應(yīng)力強度因子均取板厚中心處最大值。

圖7不同網(wǎng)格尺寸計算結(jié)果對比Fig.7 SIF in different mesh refinement

圖8應(yīng)力強度因子厚度方向變化Fig.8 SIF in different thickness

2 焊接殘余應(yīng)力下應(yīng)力強度因子的模擬方法

2.1 焊接殘余應(yīng)力對裂紋擴展的影響

在計算焊接殘余應(yīng)力分布下應(yīng)力強度因子時，需要首先確定結(jié)構(gòu)中焊接殘余應(yīng)力的分布。本文選取Faulkner模型對焊接殘余應(yīng)力分布進行模擬，具體分布形式如圖9所示，其中最大殘余拉應(yīng)力值σRESmax取為材料的屈服極限，殘余拉應(yīng)力區(qū)域范圍占板寬度的η倍，即焊接殘余拉應(yīng)力寬度為η×tplate，這里η通常取為3.5或4，而焊接殘余壓應(yīng)力的大小及范圍則由平衡條件得到。Dexter在相關(guān)研究中已經(jīng)驗證[10-11]，采用Faulkner殘余應(yīng)力分布模型模擬船舶加筋板結(jié)構(gòu)焊接殘余應(yīng)力形式合理可行。

圖9 Faulkner殘余應(yīng)力分布模型Fig.9 Faulkner’s residual stresses model

圖10裂紋區(qū)網(wǎng)格細化情況Fig.10 Mesh refinement

在計算殘余應(yīng)力分布下應(yīng)力強度因子的理論模型可基于Green函數(shù)建立，應(yīng)力強度因子KRES的計算公式表示為：

計算原理的示意圖如圖10所示。由圖中可見，一對單獨作用在裂紋面上載荷F被轉(zhuǎn)化為作用在微元dx上的分布應(yīng)力σ。

盡管焊接殘余應(yīng)力分布中既存在殘余拉應(yīng)力又存在殘余壓應(yīng)力，但對裂紋擴展速率產(chǎn)生影響的部分僅是高壓應(yīng)力區(qū)域這部分應(yīng)力，也就是說，焊接殘余應(yīng)力只會抑制疲勞裂紋的擴展，而不會使裂紋擴展速率增加。這是因為疲勞裂紋擴展的驅(qū)動力是應(yīng)力強度因子范圍ΔK，而真正對裂紋擴展起作用的是有效應(yīng)力強度因子ΔKeff，即只有使裂紋張開的那部分應(yīng)力強度因子才對裂紋擴展起作用。因此，當應(yīng)力比R≥0時，在考慮殘余應(yīng)力引起的應(yīng)力強度因子KRES與外部其他載荷引起的應(yīng)力強度因子Kapp進行疊加的情況下，焊接殘余應(yīng)力對裂紋擴展影響關(guān)系的判定關(guān)系可表示為：

2.2 焊接殘余應(yīng)力下應(yīng)力強度因子的計算

取含中心穿透裂紋及對接焊縫的平板結(jié)構(gòu)，其中長度L=2 m、寬度W=1 m、厚度B=8 mm，計算其在焊接殘余應(yīng)力區(qū)內(nèi)一系列裂紋長度下的應(yīng)力強度因子值，其中，中心穿透裂紋長度取為a=0.01~0.2 m。由Faulkner殘余應(yīng)力分布模型可知對接焊縫兩側(cè)殘余應(yīng)力的分布，如圖11所示為平板結(jié)構(gòu)尺度及其殘余應(yīng)力分布情況。

采用梯度溫度法模擬焊接殘余應(yīng)力分布時，需要在不含裂紋的完整結(jié)構(gòu)模型上施加溫度梯度，使結(jié)構(gòu)在溫度應(yīng)力作用下收縮或膨脹，從而產(chǎn)生相應(yīng)的殘余應(yīng)力分布。具體方法是：以焊縫所在位置為中心，以等間距向兩側(cè)施加3到5個溫度應(yīng)力梯度，范圍總寬度根據(jù)Faulkner模型控制在η×tplate（這里η取3.5）范圍內(nèi)。計算過程中通過多次嘗試反復(fù)迭代，選擇適合的溫度載荷值使其滿足Faulkner模型的殘余應(yīng)力分布形式。對于含中心穿透裂紋及對接焊縫平板結(jié)構(gòu)，這里模擬的殘余應(yīng)力分布為焊縫中心處最大殘余拉應(yīng)力取材料屈服極限350 MPa，分布情況如圖12所示。

圖11含中心穿透裂紋及對接焊縫平板結(jié)構(gòu)Fig.11 Center cracked plate with butt weld

圖12殘余應(yīng)力分布Fig.12 Residual stress distribution

由于平板結(jié)構(gòu)的對稱性，這里僅對平板結(jié)構(gòu)的1/2部分進行有限元建模。在確定模擬殘余應(yīng)力的溫度梯度后，在有限元模型中引入中心穿透裂紋，具體方法是將裂紋面處相應(yīng)位置的節(jié)點進行邊界釋放，從而形成中心穿透裂紋。在有限元分析過程中，為保持整體結(jié)構(gòu)的應(yīng)力平衡，原始的殘余應(yīng)力將進行重新分布以達到新的平衡。

圖13裂紋區(qū)網(wǎng)格細化情況及溫度梯度Fig.13 Mesh refinement and temperature loads

有限元模型網(wǎng)格劃分、邊界條件、溫度梯度及相應(yīng)的應(yīng)力云圖如圖13所示。取一系列裂紋長度計算殘余應(yīng)力分布下的應(yīng)力強度因子，并與理論近似計算值進行對比，結(jié)果如圖14所示。由圖中結(jié)果可知，有限元數(shù)值模擬計算結(jié)果（圖中FEM）與理論近似計算值（圖中Analysis）吻合良好?？梢?，以有限元數(shù)值模擬方法為基礎(chǔ)，采用梯度溫度載荷模擬殘余應(yīng)力分布，對應(yīng)力強度因子進行計算的結(jié)果是準確有效的。

圖14有限元模型與理論模型KRES計算結(jié)果對比Fig.14 Comparison of KRESfrom FEM and analysis

3 含中心穿透裂紋加筋板裂紋應(yīng)力強度因子計算方法

在驗證了有限元方法模擬焊接殘余應(yīng)力分布及應(yīng)力強度因子計算在含中心穿透裂紋及對接焊縫的平板結(jié)構(gòu)中的有效性后，將其進一步推廣到加筋板結(jié)構(gòu)，并將計算結(jié)果與已知文獻計算結(jié)果進行對比，證明其在含中心穿透裂紋加筋板結(jié)構(gòu)中計算的可靠性。

選取與文獻[10] 尺度相同的船體結(jié)構(gòu)中典型的加筋板結(jié)構(gòu)，詳細尺寸見圖15所示。加筋板中殘余應(yīng)力的分布情況如圖16所示，這里是加筋板結(jié)構(gòu)中殘余應(yīng)力測試結(jié)果與Faulkner模型模擬結(jié)果的對比情況，可見Faulkner模型對殘余應(yīng)力的模擬基本可以反映殘余應(yīng)力的真實分布。因此，在計算過程中仍然采用Faulkner模型殘余應(yīng)力分布對有限元模型進行加載。

為模擬加筋板中的殘余應(yīng)力分布，采用梯度溫度法，在焊縫兩側(cè)3.5tplate寬度內(nèi)施加一系列溫度載荷，計算過程中經(jīng)過多次嘗試反復(fù)迭代，確定溫度梯度，使其滿足Faulkner模型的殘余應(yīng)力分布形式。由于加筋板結(jié)構(gòu)的對稱性，僅對1/4部分的加筋板結(jié)構(gòu)進行有限元建模，邊界條件設(shè)置及有限元模型如圖17-18所示。

圖15加筋板結(jié)構(gòu)示意圖Fig.15 Stiffened panel sketch

圖16 Faulkner殘余應(yīng)力模型與試驗值比較Fig.16 Faulkner residual stress vs tests

圖17有限元模型邊界條件及加載情況Fig.17 Boundary constrains and loading

圖18加筋板有限元模型Fig.18 Finite element model of stiffened panel

取一系列裂紋長度，計算殘余應(yīng)力分布影響下的應(yīng)力強度因子，并與理論近似計算值和文獻計算結(jié)果進行對比，結(jié)果如圖19所示?？梢?，理論近似計算值（圖中Analysis）與文獻計算結(jié)果（圖中Reference）吻合情況較好，而有限元數(shù)值模擬計算結(jié)果（圖中FEM）僅在焊縫附近與其他兩種計算方法存在一定差距，其余位置吻合情況較好。由于鄰近焊縫周圍，焊接殘余應(yīng)力表現(xiàn)為拉應(yīng)力，而拉伸殘余應(yīng)力對裂紋擴展速率沒有影響，所以這里的差異并不會影響最終計算結(jié)果。因此，以有限元數(shù)值模擬方法為基礎(chǔ)，采用梯度溫度載荷模擬殘余應(yīng)力分布，并對應(yīng)力強度因子進行計算的方法是可行且可靠的。

圖19 KRES計算結(jié)果對比Fig.19 Comparison of KRES

4 考慮殘余應(yīng)力的裂紋擴展數(shù)值模擬流程

應(yīng)力強度因子計算是裂紋擴展模擬的關(guān)鍵環(huán)節(jié)。得到了符合載荷特性與結(jié)構(gòu)特征的應(yīng)力強度因子計算值，便可進一步開展考慮殘余應(yīng)力裂紋擴展的數(shù)值模擬工作。

根據(jù)Paris公式的定義，疲勞裂紋擴展速率如下式所示：

其中：N為載荷循環(huán)的次數(shù)；C和m為疲勞裂紋擴展材料參數(shù)。將Paris公式整理并轉(zhuǎn)化為積分形式：

在對裂紋擴展過程進行模擬過程中，首先給定初始裂紋長度a0，根據(jù)有限元數(shù)值計算結(jié)果可得到對應(yīng)的應(yīng)力強度因子范圍ΔK0，選擇合適的載荷循環(huán)次數(shù)ΔN，則可根據(jù)（9）式計算得到這一時刻的裂紋長度增量Δa0，將初始裂紋長度a0與裂紋長度增量Δa求和，如（10）式所示，則可得到新一時刻的初始裂紋尺寸ai。

以此類推，經(jīng)過逐步迭代計算，則可獲得相應(yīng)載荷循環(huán)次數(shù)所對應(yīng)的裂紋長度。這里可將整個裂紋擴展計算過程采用APDL參數(shù)化設(shè)計語言進行編程處理，形成計算模塊嵌入有限元計算軟件中。裂紋擴展數(shù)值計算模塊主要分析流程如圖20所示。

于是，裂紋長度增量Δa則可表示為：

圖20模塊計算分析流程圖Fig.20 Flow chart of APDL module analysis process

5 數(shù)值計算方法的試驗驗證

5.1 加筋板裂紋擴展試驗概況

加筋板長1 000 mm，帶板寬560 mm，帶板中心位置為加強筋結(jié)構(gòu)，加強筋面板寬90 mm，厚8 mm，腹板高136 mm，厚6 mm，加筋板兩端圓孔為裝夾用螺栓孔。加筋板模型中間部分，加強筋與帶板相交處，采用線切割方法分別在加筋板腹板和帶板上加工長度相同的預(yù)制裂紋，位置及尺寸如圖21所示。

圖21加筋板試樣幾何尺寸及預(yù)制裂紋Fig.21 Detailed geometry of specimen and wire-electrode cutting pre-cracks

表2試樣加載信息Tab.2 Test specimen matrix

試驗采用單側(cè)軸向加載，應(yīng)力比為0.1，載荷形式為正弦循環(huán)載荷，共進行5組試驗，各組試驗的加載信息如表2所示。試驗過程中，每隔相等時間間隔記錄裂紋長度，試樣端部疲勞裂紋遠場應(yīng)力監(jiān)測持續(xù)試驗整個過程。圖22所示為試樣安裝就位后的情況。

5.2 試驗結(jié)果和數(shù)值計算結(jié)果的比較

將拉伸載荷作用下，加筋板帶板裂紋擴展有限元數(shù)值計算結(jié)果與試驗結(jié)果對比，結(jié)果如圖23所示。由殘余應(yīng)力作用下應(yīng)力強度因子的計算結(jié)果可知，拉伸載荷作用下加筋板帶板上殘余應(yīng)力壓應(yīng)力值較小，對裂紋擴展速率影響不大，因而在裂紋擴展過程模擬中是否計及殘余應(yīng)力對計算結(jié)果影響不大，總體上計算結(jié)果與試驗結(jié)果吻合情況良好。

圖22加筋板試樣安裝就位Fig.22 Specimen in place

圖23帶板裂紋數(shù)值計算結(jié)果與試驗結(jié)果對比Fig.23 Comparison of simulation results and test results for plate

圖24腹板裂紋數(shù)值計算結(jié)果與試驗結(jié)果對比Fig.24 Comparison of simulation results and test results for web

加筋板腹板裂紋擴展有限元數(shù)值計算結(jié)果與試驗結(jié)果對比情況如圖24所示。由殘余應(yīng)力作用下應(yīng)力強度因子的計算結(jié)果可知，加筋板腹板殘余壓應(yīng)力較大，且由于殘余應(yīng)力對裂紋擴展速率的影響是一個累積的效果，因而在裂紋擴展的中后期，殘余壓應(yīng)力的效果開始顯現(xiàn)，從結(jié)果上看，計及殘余應(yīng)力的計算結(jié)果與不考慮殘余應(yīng)力的結(jié)果開始產(chǎn)生差距，且計及殘余應(yīng)力的計算結(jié)果明顯更加符合試驗中的實測結(jié)果。同時可知，本研究中加筋板腹板裂紋受殘余應(yīng)力影響較大，且對裂紋擴展有一定的抑制作用。

6 結(jié) 論

本文開展了焊接殘余應(yīng)力影響下加筋板結(jié)構(gòu)裂紋應(yīng)力強度因子計算方法研究，分析了焊接殘余應(yīng)力對裂紋擴展的影響，給出了常規(guī)外載荷以及殘余應(yīng)力下應(yīng)力強度因子的數(shù)值模擬方法，建立了考慮殘余應(yīng)力的裂紋擴展數(shù)值模擬流程，并對數(shù)值計算方法進行了試驗驗證，為后續(xù)含復(fù)雜組合裂紋加筋板結(jié)構(gòu)裂紋擴展模擬奠定基礎(chǔ)。通過本文研究，可得如下結(jié)論：

（1）通過對不同網(wǎng)格劃分規(guī)則下，含中心裂紋平板結(jié)構(gòu)應(yīng)力強度因子理論計算值與有限元數(shù)值模擬值的比較結(jié)果可知，網(wǎng)格尺度為1 mm時，計算精度最好，與理論值誤差可控制在0.5%以下；

（2）給出了焊接殘余應(yīng)力對裂紋擴展影響的判斷關(guān)系表達式。分析結(jié)果表明，焊接殘余應(yīng)力不會加速裂紋擴展速率，但是在特定條件下會對裂紋擴展速率起到抑制作用；

（3）提出了基于溫度梯度法模擬焊接殘余應(yīng)力的有限元方法，分別對含中心裂紋平板、含初始裂紋的加筋板結(jié)構(gòu)殘余應(yīng)力下應(yīng)力強度因子進行了計算，計算結(jié)果與理論值和試驗值吻合情況良好。