表面裂紋疲勞擴展和壽命計算的高效高精度數值分析方法

柴國鐘,呂君,鮑雨梅,姜獻峰,丁浩

1.浙江工業大學 機械工程學院，杭州 310012 2.義烏工商職業技術學院 機電信息學院，義烏 322000

表面裂紋疲勞擴展和壽命計算的高效高精度數值分析方法

柴國鐘1，*,呂君1, 2,鮑雨梅1,姜獻峰1,丁浩1

1.浙江工業大學 機械工程學院，杭州 310012 2.義烏工商職業技術學院 機電信息學院，義烏 322000

基于彈性力學理論，建立了三維裂紋應力強度因子計算的混合邊界元法基本理論和數值求解技術；針對表面裂紋疲勞擴展過程中，需要計算每個裂紋擴展步下的應力強度因子，從而需要重復計算大型非對稱系數矩陣問題，提出了僅在初始裂紋狀態下一次計算主控矩陣，對于隨后的疲勞裂紋擴展，只需做非常小規模矩陣的計算，且以顯式形式給出應力強度因子解而無需求解大型線性代數方程組的方法，大大提高了計算效率；針對疲勞裂紋擴展過程中，單元需要不斷重新劃分問題，由于混合邊界元法中的主控矩陣與裂紋無關，故只需對裂紋表面單元進行重新劃分，對半橢圓表面裂紋，由于將其映射到單位半圓上劃分單元，而單位半圓上的單元在疲勞擴展過程中不變，從而通過映射關系自動重新劃分裂紋表面單元。最后，通過若干算例和試驗，考核了本文方法的精度和可靠性。本文的研究為工程結構表面裂紋疲勞擴展和壽命計算的高效高精度數值分析建立了理論基礎和實現方法。

表面裂紋；疲勞擴展；高效高精度分析；混合邊界元法；應力強度因子

表面裂紋是航天航空、機械、能源、海洋工程等裝備和結構部件中最常見的缺陷形式[1]，其疲勞斷裂也常常是這些裝備破壞的主要原因之一，因此，表面裂紋的疲勞擴展和壽命分析也是工程結構強度分析的重要內容之一。

有限元法的主要優點是通用性強，可以分析任意幾何形狀受任意載荷的裂紋體，且精度高。其主要缺點是需要對整個裂紋體進行離散，用于表面裂紋疲勞擴展分析，由于疲勞擴展過程中裂紋形狀和尺寸不斷變化，對于每個裂紋擴展步，裂紋前沿單元都要重新劃分和反復計算應力強度因子，非常繁雜，用于工程實際裝備和結構的表面裂紋疲勞擴展和壽命分析具有很大的困難。

本文基于彈性力學理論，建立一種表面裂紋應力強度因子計算和疲勞擴展分析的混合邊界元法基本原理和數值求解技術，該方法的主要優點是，只需在初始裂紋狀態下對整個結構進行一次完整的數值計算，對于隨后的每個疲勞裂紋擴展步，只需對裂紋表面單元(隨著疲勞擴展自動跟進劃分)做非常小規模矩陣的重復計算，且每一擴展步下計算的應力強度因子是“精確”(數值意義上)的。

若干算例和試驗的考核表明，本文方法具有高效和高精度的優點。本文的研究為工程結構表面裂紋疲勞擴展和壽命計算的高效高精度數值分析建立了理論基礎和實現方法。

1 基本理論與分析方法

1.1 混合邊界元法

對于圖1所示的三維表面裂紋體，混合邊界積分方程為[10]

(1)

(2)

即當奇異點p位于非裂紋邊界Г=ΓⅠ+ΓⅡ上時，采用利用單位力基本解(第一基本解)的第一類邊界積分方程，而當奇異點p位于裂紋面ГC上時，采用利用單位位移不連續基本解(第二基本解)的第二類邊界積分方程。

為了數值求解式(1)和式(2)，將圖1所示裂紋體的整個邊界離散劃分為n個單元，即

式中：

(3)

(4)

圖1 表面裂紋
Fig.1 Surface crack

將式(4)代入式(1)和式(2)，得到離散后的邊界積分方程為

cij(p)uj(p)=

(5)

ti(p)=

(6)

將整個裂紋體邊界分為圖1所示3類，并采用相應的3類單元：第一類邊界(通常邊界ΓⅠ)——Ⅰ類單元；第二類邊界(裂紋嘴上下邊界(含裂紋擴展)ΓⅡ)——Ⅱ類單元；第三類邊界(裂紋表面ГC)——裂紋單元。

第一類邊界(通常邊界ΓΙ)：采用通常的4-8節點等參元，單元形式如圖2所示，其插值函數為Ni(ξ1,ξ2)(i=1,2,…,8)[11]。

因此單元內任一點的整體坐標、位移和面力可分別表示為

(7)

第二類邊界(裂紋嘴上下邊界(含裂紋擴展)ΓⅡ)：由于裂紋嘴上下表面位移不連續，因此不能采用與第一類邊界相同的單元插值函數。本文對該類邊界提出如下的非協調等參元。

圖2 4-8節點等參元
Fig.2 Four-eight nodal iso-parameter element

圖3 邊界的非協調等參元
Fig.3 Incompatible iso-parameter element on boundary

邊界的單元劃分如圖3(a)所示。對于單元1，將圖2所示的通常的8節點單元的2、5節點向單元內移到ξ2=-1/2處；對于單元2，將節點1、5、2向單元內移到ξ2=-1/2處；對于單元3，將節點1、5向單元內移到ξ2=-1/2處，如圖3(b)～圖3(d)所示。

對于單元1有

(8a)

對于單元2有

(8b)

對于單元3有

(8c)

式中：Nl(l=1,2,…,8)為圖2所示通常的4-8節點等參元插值函數。

(9)

對于圖5所示橢圓裂紋，其相對位移基本函數為

(10)

橢圓及部分橢圓裂紋是工程結構中最常見的三維裂紋，結構中的深埋裂紋通常可以簡化成深埋橢圓裂紋，平板表面裂紋及圓筒中的軸向表面裂紋通常可以簡化為半橢圓表面裂紋，孔邊角裂紋通常可以簡化為1/4橢圓表面裂紋等，而這些裂紋模型都可看做是圖5所示橢圓裂紋中截取的一部分,因此它們的相對位移基本函數都可以用式(10)表示，只是對于不同的裂紋，它們的相對位移權函數不同。

然后利用式(9)將式(5)和式(6)中的最后一個積分項分別改寫為

(11)

圖4 平片裂紋
Fig.4 Planar crack

圖5 橢圓裂紋
Fig.5 Elliptical crack

(12)

(13)

最后將式(13)代入式(12)，再返回式(5)和式(6)，式(5)和式(6)可進一步寫為

圖6 半線性權函數單元
Fig.6 Semi-linear weight function element

(14)

(15)

1.2 數值求解技術

由式(14)和式(15)可知，在離散形式的邊界積分方程形成過程中，需計算大量的如式(16)～式(18)所示的奇異積分。

(16)

(17)

Sk∈ΓC

(18)

對于1/R的奇異積分式(16)和1/R2的奇異積分式(17)采用通常邊界元中的方法解決[11]。

對于1/R3的超奇異性積分，首先將式(18)改寫為

(19)

不失一般性，設單元Sk位于Ox1x2平面，n(p)=n(q)=(0,0,1),如圖7所示。令

Fl(q)=Nl(q)f(q)

(20)

并將其在奇異點p附近展開為Taylor級數,即

(21)

其中:

圖7 裂紋單元超奇異積分計算
Fig.7Hypersingular integral calculation of crack element

(22)

(23)

將式(21)代入式(19)，將式(19)表示為兩個積分之和

Il=I0+I

(24)

其中:

(25)

(26)

因此，式(19)的積分等價于式(25)和式(26)表示的兩個積分。

對于積分I0，由式(23)可知，由于ΔFlp,q=OR3而使被積函數的奇異性消除，因此可直接利用Gauss求積公式計算。

I=2μA·

(27)

由式(27)可見，為了計算I值，只需研究有限部分積分

(28)

其中k1和k2的取值為：k1=k2=0;k1=k2=1;k1=2,k2=0;k1=0,k2=2。

如圖7所示，令

x1-ξ1=rcosθ

x2-ξ2=rsinθ

則

(29)

這樣就將二維的有限部分積分化為一維的有限部分積分。將式(22)代入式(29)，對r積分并取有限部分，得

(30)

其中:

(31)

式(31)中的被積函數已無奇異性，可用數值積分進行計算。為了便于利用Gauss求積分式計算，將積分變量再作如下變換。

設單元Sk的周邊由M條簡單曲線段所組成，如圖7所示，將式(30)對θ的積分變換為對描述單元周邊各曲線段的廣義坐標γ的積分。對于單元周邊第m條曲線段，有

(32)

Jm可由式(31)和式(32)導出,即

(33)

將式(32)和式(33)代入(30)，得到

(34)

對于(部分)橢圓，將其映射到(部分)單位圓上劃分單元，映射關系為

(35)

橢圓裂紋及單位圓上的單元劃分如圖8所示。

由式(35)，可以方便地導出式(33)的Jm。

(36)

J1

(37)

(38)

(39)

圖8 橢圓裂紋及單位圓上的單元劃分
Fig.8Elliptical crack and element division of unit circle

由于式(36)和式(38)分別表示以φ和ρ為參數的雙曲線族和橢圓族，故將圖8(a)所示單元稱為雙曲線-橢圓單元。

求解式(14)和式(15)，得到基本線性代數方程組:

(40)

由此得到裂紋前緣周邊任意D點(見圖4)所在單元的權函數Wj(D)j=z,n,t,進而確定D點的應力強度因子為

(41)

式中：E為彈性模量。對于(部分)橢圓裂紋(見圖5)。

1.3 疲勞裂紋擴展和壽命分析方法

工程中計算疲勞裂紋擴展的一般形式的方程為

(42)

式中：αi(i=1,2,…,m)為材料常數和環境及其他力學量。其中最常用的是Paris公式(見式(43),對于表面裂紋的表面點，一般認為其閉合效應與裂紋最深點不同，通常采用ΔKeff=0.9ΔK代替式(43)中的ΔK，本文也采用該方法)。

(43)

或

(44)

式中：a0和ac分別為初始和終止裂紋深度。

由式(41)、式(43)和式(44)可知，要高效和精確計算表面裂紋疲勞擴展，需要解決兩個關鍵問題：① 每一個裂紋擴展步下，由于需要計算當前的應力強度因子(幅)，因此對應于每一個擴展步，單元網格的自動生成;② 由于式(40)左邊為大型、滿秩、非對稱矩陣，計算矩陣系數和求解式(40)需要耗費大量的時間，計算效率很低。

本文通過以下技術解決以上兩個關鍵問題，從而實現表面疲勞裂紋擴展的高效高精度計算。

1) 通過適當的排列，使得形成的線性代數方程組具有以下形式：

(45)

式中：A0為主控矩陣，其只與邊界Γ的幾何形狀和邊界條件有關(即式(14)右邊的前2個積分)，而與裂紋擴展無關，故只需在初始裂紋狀態下計算一次，在隨后的疲勞裂紋擴展過程中，無需再重復計算，而只需計算矩陣B1、B2和C，而這3個矩陣占整個矩陣的計算量為

式中：NC和N0分別為裂紋單元節點數和非裂紋單元節點數，對于一般的表面裂紋問題，NC<100(本文研究中取NC=88)，而對于一般的實際工程構件，N0一般為幾千到幾萬，故每一個疲勞裂紋擴展步下的矩陣計算量約比初始裂紋狀態下的矩陣計算量小1～2個量級。

2) 通過矩陣運算，求得式(45)左邊系數矩陣的逆矩陣為

(46)

其中:

由此可以求得式(45)的解為

(47)

由于在疲勞裂紋擴展過程的每個步中，只需要計算當前裂紋的應力強度因子(式(47)中對應的未知量U)，故式(47)還可以進一步簡化為

(48)

式中：A0只需要在裂紋初始狀態下計算一次，隨后的每次疲勞裂紋擴展，只需計算B1、B2和C，又由于式(48)給出了U的顯式解，無需求解大型線性代數方程組，大大提高了計算效率。

3) 對于工程中的疲勞裂紋，例如最常見的半橢圓表面裂紋、拐角1/4橢圓表面裂紋等，將它們均映射到(部分)單位圓上劃分單元(見圖8)。疲勞裂紋擴展中雖然裂紋尺寸a、b發生變化，但單位圓上的單元始終不變，因此疲勞擴展過程中，通過映射關系式(35)，實際物理平面上的單元自動調整(自動重新劃分)。

由于針對每一個新的裂紋擴展步，重新精確計算(數值意義上)當前應力強度因子，所以是一種高精度的疲勞裂紋擴展計算方法；而對于每一個新的裂紋擴展步，計算工作量只是第一次計算的1/10～1/100左右，效率高，且整個計算只需一次輸入初始參數，隨后全部自動完成。

2 算例與試驗考核

2.1 算 例

算例主要考核混合邊界元法計算應力強度因子的精度。圖9(a)所示遠場受均勻拉伸應力σ作用的板內半橢圓表面裂紋，幾何尺寸為：a/w=a/h=0.2。裂紋面上單元劃分為Nφ×Nρ=11×4=44，b/a=0.2與b/a=1.0的計算結果與Raju和Newman, Jr[2]利用三維有限元法用了近萬個自由度在大型計算機上求得結果進行比較，如圖9(a)和圖9(b)所示，圖中Φ為第二類橢圓積分。由圖9可知，兩者結果吻合很好。

算例表明，本文所建立的混合邊界元法是可靠的，用于計算表面裂紋應力強度因子與有限元等方法具有同等的精度。

圖9平板半橢圓表面裂紋受均勻拉伸應力時的 應力強度因子
Fig.9 Stress intensity factor of semi-elliptical surface crack in plate subjected to uniform tensile stress

2.2 表面裂紋疲勞擴展試驗考核

為了說明混合邊界元法計算疲勞裂紋擴展的高效和高精度，進行試驗考核。

試樣為圖9(a)所示遠場受拉伸應力σ的平板半橢圓表面裂紋試樣，材料為16Mn鋼，幾何尺寸為w×t×h=38 mm×10 mm×120 mm、初始裂紋尺寸和應力幅如表1所示。初始裂紋通過厚度為0.13 mm的砂輪加工而成，共3個試樣。

試驗在Instron試驗系統上進行，試驗過程中“降載留痕”，試驗結束后，打開斷口，進行高溫氧化，通過“降載留痕+氧化著色”可以清晰地看到斷口上的“Beachmark”，由此測得疲勞裂紋擴展量與循環次數的關系。

將3個試樣的試驗結果用Paris公式(見式(43))進行疲勞裂紋擴展速率計算(回歸分析)，得到16Mn鋼Paris公式材料常數為：C=4.41×10-8，n=3.61。數值分析中，利用對稱性，取1/4幾何模型，模型表面單元劃分和邊界約束如圖10所示，表面節點數和單元數分別為1 185和372(通常邊界ΓⅠ上的單元有357個，邊界ΓⅡ上非協調單元有15個)，裂紋面上單元數為Nφ×Nρ=11×4=44。

表1 試樣幾何尺寸、初始裂紋尺寸與壓力幅

圖10 數值分析模型
Fig.10 Numerical analysis model

疲勞裂紋擴展及應力強度因子的計算結果與試驗結果如圖11所示。由圖11可見，計算結果與試驗結果非常一致。

圖11 表面裂紋疲勞擴展數值分析與試驗結果的比較
Fig.11Comparison of numerical analysis and test results of fatigue growth of surface crack

3 結 論

1) 建立的混合邊界元法基本方程和數值求解技術基于精確的彈性力學理論和數學方法，因此是一種嚴格可靠的數值分析方法，可用于計算任意三維裂紋應力強度因子和基于線彈性斷裂力學的疲勞裂紋擴展分析，數值算例表明，混合邊界元法與有限元法具有同等級別的精度。

2) 混合邊界元法用于表面裂紋疲勞擴展計算，對于每個裂紋擴展步均能精確(數值意義)計算應力強度因子，因此這是一種高精度計算表面裂紋疲勞擴展的方法。

3) 除了在初始裂紋狀態下需要計算大型非對稱系數矩陣外，對于隨后的疲勞裂紋擴展，只需做小規模矩陣的計算(結構越大，這種優勢越明顯)，且以顯式形式給出應力強度因子解而無需求解大型線性代數方程組，大大提高了計算效率。

4) 疲勞裂紋擴展過程中，單元需要不斷重新劃分，由于混合邊界元法中的主控矩陣與裂紋無關，故只需對裂紋表面單元進行重新劃分，對半橢圓表面裂紋，由于單位半圓上單元在疲勞擴展過程中不變，從而通過映射關系自動重新劃分裂紋表面單元。

綜上所述，混合邊界元法用于表面裂紋疲勞擴展計算，只需一次輸入初始數據即可進行全疲勞過程的高效和高精度自動分析。本文的研究為工程結構表面裂紋疲勞擴展和壽命計算的高效高精度數值分析建立了理論基礎和實現方法。

[1] 胡博, 于潤橋, 徐偉津. 人工槽模擬GH4169渦輪盤表面裂紋缺陷的微磁檢測[J]. 航空學報, 2015, 36(10): 3450-3456.

HU B, YU R Q, XU W J. Micro-magnetic NDT for surface crack defect in a GH4169 turbine disc simulated by artificial groove[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3450-3456 (in Chinese).

[2] RAJU I S, NEWMAN J C, Jr. Stress-intensity factors for a wide range of semi-elliptical surface cracks in a finite-thickness plate[J]. Engineering Fracture Mechanics, 1979, 11(4): 817-829.

[3] 國家標準化管理委員會. 在用含缺陷壓力容器安全評定: GB/T 19624—2004[S]. 北京: 中國標準出版社, 2004.

Standardization Administration of China. Safety assessment for in-service pressure vessels containing defects: GB/T 19624—2004[S]. Beijing: China Standard Press, 2004 (in Chinese).

[4] LIN X B, SMITH R A. Finite element modelling of fatigue crack growth of surface cracked plates: Part I: The numerical technique[J]. Engineering Fracture Mechanics, 1999, 63(5): 503-522.

[5] LIN X B, SMITH R A. Finite element modelling of fatigue crack growth of surface cracked plates: Part II: Crack shape change[J]. Engineering Fracture Mechanics, 1999, 63(5): 523-540.

[6] LIU C H, CHU S K. Prediction of shape change of corner crack by fatigue crack growth circles[J]. International Journal of Fatigue, 2015, 75: 80-88.

[7] YU P S, GUO W L. An equivalent thickness conception for prediction of surface fatigue crack growth life and shape evolution[J]. Engineering Fracture Mechanics, 2012, 93: 65-74.

[8] GUCHINSKY R, PETINOV S. Numerical modeling of the surface fatigue crack propagation including the closure effect[J]. International Journal for Computational Methods in Engineering Science and Mechanics, 2016, 17(1): 1-6.

[9] 李有堂, 張洋洋. 壓力容器中表面裂紋在高周疲勞下的擴展規律[J]. 蘭州理工大學學報, 2015, 41(6): 168-172.

LI Y T, ZHANG Y Y. Surface crack growth pattern of pressure vessel with high-cycle fatigue[J]. Journal of Lanzhou University of Technology, 2015, 41(6): 168-172 (in Chinese).

[10] CHAI G Z, FANG Z M, JIANG X F, et al. Hybrid boundary element analysis for surface cracks[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(23-26): 2069-2086.

[11] WATSON J O. Advanced implementation of the boundary element method for two and three-dimensional elastostatics[J]. Developments in Boundary Element Methods, 1979, 61: 31-63.

[12] IOAKIMIDIS N I. Application of finite-part integrals to the singular integral equations of crack problems in plane and three-dimensional elasticity[J]. Acta Mechanica, 1982, 45(1-2): 31-47.

Ahighlyefficientandaccuratenumericalanalysismethodforfatiguepropagationofsurfacecrackandlifeprediction

CHAIGuozhong1,*,LYUJun1, 2,BAOYumei1,JIANGXianfeng1,DINGHao1

1.CollegeofMechanicalEngineering,ZhejiangUniversityofTechnology,Hangzhou310012,China2.SchoolofMechatronics&IT,YiwuIndustrial&CommercialCollege,Yiwu322000,China

Thebasictheoryandnumericalsolvingtechniqueofthehybridboundaryelementmethodforthecalculationofthestressintensityfactorsofthree-dimensionalcrackisestablishedbasedonthetheoryofelasticmechanics.Inanalyzingthefatiguepropagationofthesurfacecrack,thestressintensityfactorateachcrackpropagatingstepneedstobecalculated,andaccordinglythelargenon-symmetriccoefficientmatrixshouldbecomputedrepeatedly.Amethodisproposedthatthemastermatrixiscalculatedonlyonceintheinitialcrackstate,andthenaverysmall-scalematrixiscalculatedduringthesubsequentfatiguecrackpropagation.Thesolutionforthestressintensityfactorisalsogiveninanexplicitformwithoutsolvinglarge-scalelinearalgebraicequationssothatthecalculationefficiencyisimprovedgreatly.Toaddresstheproblemofcontinuousdivisionandremeshingofelementsduringthefatiguecrackpropagating,thehybridboundaryelementmethodisappliedasthemastermatrixisindependentofthecrack,andthereforeonlyremeshingofelementsonthecracksurfaceisrequired.Forthesemiellipticalsurfacecrack,theelementsonthecracksurfaceareremeshedaccordingtothemappingrelationshipasthecrackismappedintothesemicircleinmeshingandtheelementsintheunitsemicirclehavenochangeduringthefatiguepropagation.Theaccuracyandreliabilityoftheproposedmethodareverifiedbyseveralexamplesandexperiments.Theresearcheffortsmayprovidethetheoreticalfoundationandtherealizationmethodforhighlyefficientandaccuratenumericalanalysisofthesurfacecrackfatiguepropagationandthelifepredictionofengineeringstructures.

surfacecrack；fatiguepropagation；highlyefficientandaccurateanalysis；hybridboundaryelementmethod；stressintensityfactor

2017-03-29；

2017-04-06；

2017-04-10；Publishedonline2017-04-261801

URL：http://hkxb.buaa.edu.cn/CN/html/20171217.html

NationalNaturalScienceFoundationofChina(51275471)

.E-mailchaigz@zjut.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.221291

2017-03-29；退修日期2017-04-06；錄用日期2017-04-10；網絡出版時間2017-04-261801

http://hkxb.buaa.edu.cn/CN/html/20171217.html

國家自然科學基金(51275471)

.E-mailchaigz@zjut.edu.cn

柴國鐘,呂君,鮑雨梅，等. 表面裂紋疲勞擴展和壽命計算的高效高精度數值分析方法J. 航空學報，2017,38(12):221291.CAIGZ,LYUJ,BAOYM,etal.AhighlyefficientandaccuratenumericalanalysismethodforfatiguepropagationofsurfacecrackandlifepredictionJ.ActaAeronauticaetAstronauticaSinica,2017,38(12):221291.

V215.5+5

A

1000-6893(2017)12-221291-12

徐曉)