雙向應力場中表面裂紋應力強度因子的權函數法

袁奎霖，周忠華，趙 峰，洪 明

（大連理工大學 船舶工程學院 工業裝備結構分析國家重點實驗室，遼寧 大連116024）

0 引 言

船舶等大型焊接結構在交變外載荷作用下，在應力集中較嚴重的焊趾部位極易萌生多處微裂紋,隨著微裂紋生長與合體形成近似于半橢圓形狀的表面裂紋[1]。表面裂紋沿板厚和板寬方向不斷擴展，直至貫穿鋼板形成穿透型裂紋甚至引發斷裂破壞。為了保證結構服役安全性，表面裂紋尖端應力強度因子的計算是損傷容限設計中極其重要的一部分。

對于半橢圓表面裂紋應力強度因子的計算，已有國內外學者對此進行了大量研究，其中有限元法是最常用的方法之一。Newman和Raju[2-3]通過建立三維有限元模型對遠端拉伸載荷和彎矩作用下的有限厚度平板表面裂紋應力強度因子進行了分析，并且總結了相應的經驗公式。Bowness等[4]針對T型接頭中的表面裂紋，結合三維有限元分析，給出了考慮焊趾處應力集中影響的應力強度因子放大系數Mk的計算公式。Shiratori等[5]基于線性疊加原理將遠場載荷與結構局部應力集中的影響轉化為裂紋面上的近場應力載荷，分別計算了在0～3次多項式表達的沿板厚方向變化應力場作用下平板表面裂紋應力強度因子。然而，實際焊接結構由于受形狀復雜性和焊接殘余應力的影響，焊趾處應力載荷往往沿板寬和板厚方向同時變化（如圖1所示），并且需要高次多項式進行擬合[6]，這對應力強度因子的有限元計算造成了很大困難。

與有限元法相比，權函數法是求解復雜應力場作用下裂紋應力強度因子的高效計算方法。應力強度因子的權函數法最早由Bueckner[7]和Rice[8]提出，并證明了權函數是裂紋體的幾何特性，一旦確定便可用來計算任意載荷分布下的應力強度因子，所需計算僅是一個與裂紋面應力乘積的積分。Glinka等[9-10]基于權函數與裂紋開口位移之間的關系，提出了適用于穿透型裂紋和半橢圓表面裂紋的一維權函數統一形式，并被美國石油協會規范API-579[11]所采用。Wang等[12-13]在已有統一權函數基礎上，通過有限元計算和Shiratori[5]的數據，分別發展了小半長比（a/c=0.05～1.0）和大半長比（a/c=0.6～2.0）的平板表面裂紋權函數。然而，上述權函數僅適用于裂紋面上應力載荷單向變化（即沿板厚方向）的情況，對于實際工程中遇到的雙向變化應力情況則難以直接應用。針對這一問題，Wang和Glinka[14]基于無限體中裂紋面上集中力載荷誘導的應力強度因子解析解的特性，提出了一種任意載荷作用下橢圓形埋藏裂紋的二維權函數統一形式。Jin和Wang[15]考慮了有限板厚的邊界效應，進一步提出了平板表面裂紋的二維統一權函數，但遺憾的是仍只驗證了沿板厚方向的最高三次分布應力。此外，現有權函數[11-15]的裂紋深度比適用范圍主要集中在a/T=0.2～0.8,而對于焊接接頭的初始表面裂紋尺寸（例如裂紋初始深度a0=0.5 mm，板厚T=20 mm[16]，a/T=0.025）仍需要改進。

圖1 焊趾表面裂紋附近應力場Fig.1 2-D stress field near weld toe

圖2 有限厚度平板半橢圓表面裂紋形狀與坐標系Fig.2 Geometry and coordinate system for semi-elliptical surface crack in finite thickness plate

本文基于Wang等[14-15]提出的二維權函數統一形式，通過創建裂紋半長比a/c=0.05～1.0，裂紋深度比a/T=0.01～0.8的平板表面裂紋三維有限元模型，分別計算了裂紋最深點和表面點的應力強度因子并將其作為參考解，提出一組裂紋形狀適用范圍更廣的表面裂紋二維權函數。在此基礎上，通過在裂紋面上施加兩組不同的最高階次為六次的復雜應力載荷進行計算，與有限元結果對比，驗證了本文所提出的權函數的準確性。

1 二維權函數統一形式

基于疊加原理，理論上采用權函數法可以計算任意載荷條件下的應力強度因子。對于一維貫穿型裂紋，基于權函數法的裂紋應力強度因子計算公式如下：

式中：x為裂紋面坐標，a為裂紋長度，σ（x）為無裂紋體假想裂紋處應力分布，m（x,a）為權函數。m（x,a）表示作用于x點處的單位集中力載荷誘導的裂紋尖端應力強度因子。Glinka和Shen[9-10]指出，對于一維和二維裂紋，權函數可由（2）式的統一形式表示：

式中：M1、M2和M3為權函數系數。需指出，對于半橢圓表面裂紋，（2）式僅適用于應力分布沿板厚單向變化的情況。

對于雙向變化應力場中的二維裂紋，應力強度因子可由權函數m（x,y;P′）和應力分布σ（x,y）在裂紋面全域S上的雙重積分計算得到：

可知，m（x,y;P′）表示（x,y）點處的單位集中力載荷在裂紋前緣P′點處誘導的應力強度因子。Wang等[14-15]對于圖2所示平板半橢圓表面裂紋提出了二維權函數統一形式：

Ghajar等[17]研究表明，當n=1時即如（5）式所示的權函數統一形式對于不同形狀的半橢圓表面裂紋都具有良好的計算精度：

圖3 半橢圓表面裂紋各參數示意圖Fig.3 Parameters for semi-elliptical surface crack

圖4 表面裂紋有限元模型Fig.4 Typical FE model for surface crack

式中：s是載荷作用點P到裂紋前緣的最短距離，ρ是載荷作用點P到應力強度因子計算點P′的距離，θ代表P′在裂紋前緣的位置，r與φ為載荷作用點P的極坐標（見圖3）。求解權函數中的未知系數M，一般需要已知某種載荷下的一系列高精度應力強度因子作為參考解，可采用文獻中已知數據或進行有限元計算得到。

2 平板表面裂紋有限元分析

由于現有文獻中關于淺表面裂紋（a/T＜0.2）的應力強度因子數據鮮有報道，本文自行創建了裂紋半長比a/c=0.05、0.1、0.2、0.4、0.6、0.8和1.0，裂紋深度比a/T=0.01、0.05、0.1、0.2、0.4、0.6和0.8，共49個平板表面裂紋三維有限元模型。如圖4所示，采用對稱邊界條件創建了1/4裂紋模型，載荷加載方式為在裂紋面上直接施加應力分布載荷，材料模型為線彈性體，楊式模量E=2.06E5 MPa，泊松比υ=0.3。有限元軟件為ANSYS 16.0，采用的單元為20節點的SOLID 186高階體單元，其中裂尖單元采用奇異單元處理，應力強度因子計算方法為位移外插法。

為了驗證有限元模型的準確性，在裂紋面上施加如下四種應力場：

式中：σ0為名義應力，a表示裂紋深度，取n=0,1,2,3。求得的應力強度因子按（7）式進行歸一化處理：

式中：F為邊界效應修正因子，Q為第二類橢圓積分的近似解。參考Shiratori[5]的結果，對最深點和表面點的邊界修正因子F的計算結果進行對比，所有結果的誤差均在5%以內，部分整理結果如表1所示。驗證結果表明本文創建的表面裂紋有限元模型具有良好的計算精度，并認為淺表面裂紋（a/T＜0.2）的應力強度因子有限元結果可作為后續的參考解。

表1 有限元計算得到的邊界修正因子F與Shiratori參考解[5]對比，a/c=0.6Tab.1 Comparison of boundary correction factors from present FEM calculation and Shiratori et al[5],a/c=0.6

續表1（b）

3 權函數求解

圖5 單元子分法示意圖Fig.5 Schematic illustration of element subdivision technique

采用上述方法分別對裂紋半寬比a/c=0.05～1.0，裂紋深度比a/T=0.01～0.8的不同形狀表面裂紋的表面點和最深點處對應的權函數參數M進行了求解，結果匯總見表2。為了方便工程應用，采用了最小二乘法對表2中的M值進行多項式擬合得到參數公式（詳見附錄），擬合優度R2均大于0.99。

4 權函數的有效性驗證

為了驗證權函數的有效性，將由一階參考載荷得到的權函數擴展應用到高階載荷的情況，權函數計算得到的高階載荷下應力強度因子與有限元結果進行對比。

4.1 沿板厚單向變化的應力分布

式中a表示裂紋深度。將（10）式所表達的沿板厚方向變化的線性分布應力和非線性分布應力載荷施加在裂紋面上，計算所得的應力強度因子經歸一化后得到邊界修正因子F。如圖6所示，將有限元法與權函數法的計算結果進行對比，并統計了兩者之間的誤差（見圖7）。除a/T=0.01，a/c=0.8和1.0的高階應力m=4～6情況，其他裂紋形狀的計算誤差均在10%以內，滿足實際工程需要。需要說明的是，相對誤差超過10%的邊界修正因子F絕對值很小且均小于0.1，因此并不影響本權函數的適用性。

圖7 權函數結果與有限元結果之間的誤差分析Fig.7 Error analysis of boundary correction factors from weight function method and FEM

4.2 沿板厚和板寬雙向變化的應力分布

在船體結構中縱骨端部（如圖1所示的焊接接頭）是疲勞強度校核的關鍵節點，焊趾附近的應力分布沿板厚和板寬方向變化。為了模擬上述情形，采用（11）式中所表達的雙向變化非線性應力載荷施加在裂紋面上，計算所得的應力強度因子經歸一化得到邊界修正因子F。如圖8所示，將有限元法與權函數法的計算結果進行對比，并統計了兩者之間的誤差（見圖9）。相對誤差超過10%的情況出現在a/T=0.01，a/c=1.0的高階應力m=4～5時，其他裂紋形狀的計算誤差均在10%以內，可滿足實際工程需要。

式中，a表示裂紋深度，c表示裂紋半寬。

圖8 權函數結果與有限元結果對比Fig.8 Comparison of boundary correction factors from weight function and FEM

圖9 權函數結果與有限元結果之間的誤差分析Fig.9 Error analysis of boundary correction factors from weight function method and FEM

5 結 論

本文在Wang等提出的二維權函數統一形式的基礎上，采用三維有限元計算得到裂紋半長比a/c=0.05、0.1、0.2、0.4、0.6、0.8和1.0，裂紋深度比a/T=0.01、0.05、0.1、0.2、0.4、0.6和0.8的平板表面裂紋最深點和表面點的應力強度因子作為參考解，得到了裂紋形狀適用范圍更廣的權函數解，并得到以下結論：

（1）由低階參考載荷得到的權函數可擴展應用到高階應力分布，本文所提出的權函數可高效準確地求解最高六階的沿板厚方向和板寬方向變化的復雜應力場作用下的應力強度因子。

（2）對于表面點和最深點，半長比a/c=0.05～1.0，a/T=0.05～0.8之間，權函數結果與有限元結果相比，誤差均在10%以內，滿足實際工程應用需要。

（3）本文所提出的權函數可應用于考慮結構應力集中、焊接殘余應力和焊后表面處理（如噴丸、超聲沖擊）等影響的實際焊接結構中表面裂紋擴展的數值預報。

裂紋表面點（θ=0）：