無限板孔邊裂紋問題的高精度解析權函數解

趙曉辰，吳學仁, *，童第華，徐武，陳勃，胡本潤

1. 中國航發 北京航空材料研究院，北京 100095 2. 上海交通大學 航空航天學院，上海 200240

無限板中的圓孔邊徑向裂紋是工程實際中常見的一類裂紋幾何，該裂紋問題的求解方法在航空航天、機械運載、土木建筑、油氣開采等工業領域中有非常廣泛的需求?？走吜鸭y在各種復雜載荷下的高精度應力強度因子(SIF)K和裂紋面張開位移(COD)等關鍵斷裂力學參量，是計算和預測含孔邊裂紋結構的疲勞壽命和剩余強度的重要前提[1]。雖然有限元法(Finite Element Method，FEM)和邊界元法(Boundary Element Method，BEM)等數值方法能求解各種復雜載荷情況下的孔邊裂紋問題，但數值方法對計算資源和使用者的經驗均有頗高的要求。特別是在需要求解裂紋長度變化不斷和(或)多種復雜載荷情況的裂紋問題時，數值方法的計算和校驗工作量都很大，經濟性較差。因此，尋求高效準確的孔邊裂紋問題求解方法有重要的實際意義。

權函數法(WFM)是求解裂紋問題的一種高效手段[2-6]。它把影響應力強度因子的變量——裂紋體的幾何(包括載荷/位移邊界劃分)和所承受的載荷作了分離。權函數僅與裂紋幾何及載荷/位移邊界(ST/SU)的劃分相關，一旦確定，則僅需通過對權函數和假想裂紋處應力分布的乘積作簡單積分就能夠得到作為裂紋長度的函數的K和COD。因此在處理復雜載荷情況或同一裂紋幾何多種載荷情況的裂紋問題，或對許多裂紋長度(例如疲勞裂紋擴展)求解時，權函數法可以在極大地提高計算效率的同時，保證結果的高精度。

權函數法總體上可以分為解析法和數值法2種。在解析權函數法方面，目前主要有：①Wu-Carlsson基于一種參考載荷下的K和CMOD參考解的權函數[4]；②Fett-Munz[5]、Glinka-Shen[6]分別發展的基于多種(一般為2種)參考載荷下的K求得的權函數。Wu-Carlsson解析權函數法能夠直接得到作為裂紋長度α連續函數形式的權函數，不但能方便地求解連續變化α的K和COD，而且所得權函數受參考載荷形式及參考解精度的影響很小[7-8]，在現有解析權函數法中，其求解精度和魯棒性均具有明顯優勢。在數值權函數方面，除了采用傳統有限元法求解權函數外，Wagner-Millwater[9]近年來提出了一種基于復變函數泰勒級數展開的權函數(Weight function Complex Taylor Series Expansion，WCTSE)法。這種方法基于復變函數的泰勒級數展開理論，結合復變有限元計算，能夠得到任意裂紋體幾何的寬α范圍離散裂紋長度的數值形式權函數。Jing和Wu[10]通過對WCTSE法的深入研究，進一步提高了其求解精度。在對復變有限元及裂紋問題建模方面積累了豐富經驗的前提下，WCTSE法所求得的數值權函數具有很高的精度[7-8, 11]。但由于WCTSE法本質上是一種數值方法，對每個裂紋尺寸都需要進行一次復變有限元計算，因此其求解效率遠低于解析權函數法。而WCTSE法的高精度特點，則使得它在對通過其他方法得到的權函數的精度評價驗證方面具有獨特的優勢，因此可以作為評價解析權函數精度的“基準”。

本文利用Wu-Carlsson解析權函數法求得了無限板孔邊單裂紋(N=1)和對稱雙裂紋(N=2)的權函數；并通過與權函數直接對應的格林函數(Green’s Function，GF)，分別與Shivakumar-Forman[12]和Newman[13]的高精度擬合表達式以及WCTSE數值權函數結果作了全面對比驗證。此外，還指出了近期文獻中的孔邊裂紋權函數[14-15]的明顯錯誤。在此基礎上用得到的解析權函數公式計算了孔邊裂紋在多種典型載荷作用下的應力強度因子，為復雜受載的孔邊裂紋問題提供了形式統一、便捷高效的高精度求解手段。

1 無限板孔邊裂紋的解析權函數

考慮無限板圓孔邊的徑向裂紋問題，即單裂紋(N=1)和對稱雙裂紋(N=2)，如圖1所示。

根據權函數理論[2-6]，對權函數m(α,γ)和無裂紋體在假想裂紋處的應力分布σ(ξ)的乘積沿裂紋長度作積分，即可得到任意載荷條件下的應力強度因子為

(1a)

(1b)

式中：α=a/W,ξ=x/W,γ=ξ/α,a和x分別為裂紋的長度和坐標；W為裂紋體的特征尺寸；σ(ξ)為無裂紋體在假想裂紋處的應力分布；σ0為σ(ξ)的歸一因子。對于無限板孔邊裂紋，取孔的半徑作為特征尺寸，即W=R，相關參量定義見圖1。這里采用Wu-Carlsson解析權函數法對孔邊裂紋問題進行分析。作為一種邊緣裂紋，孔邊裂紋的Wu-Carlsson權函數形式為[4]

(2)

(3)

式(2)表明，βi(α)是確定權函數m(α,γ)的前提，而βi(α)又由函數Fi(α)確定。所以Fi(α)的推導成為Wu-Carlsson權函數的關鍵環節。根據文獻[4]，確定邊緣裂紋Fi(α)的條件一般可包括：① 裂紋尖端的COD-K關系；② 參考載荷情況的K的自洽條件；③ 參考載荷作用下的裂紋嘴張開位移(CMOD)。Wu-Carlsson利用條件①和條件②在其權函數專著中給出了孔邊裂紋的解析權函數[4]，并進而求解得到各種復雜載荷情況下的孔邊裂紋應力強度因子[4, 16]。多位學者基于這2個條件，成功應用Wu-Carlsson權函數方法求解了各種孔邊裂紋問題。近年來的典型例子有：有限寬板圓孔邊單/雙裂紋分析以及對參考應力強度因子擬合方法的改進[17]、無限板孔邊不等長雙裂紋的混合型問題的KI和KII以及裂紋張開位移[18]、周期性孔邊共線裂紋的K和COD[19]。

本文研究發現，增加條件③能明顯提高α較小時的求解精度。值得注意的是，裂紋面位移曲線在裂紋嘴處的曲率為0同樣適用于邊緣裂紋，但對多連域裂紋問題，此條件對精度的影響很小，且會增加求解的復雜性。這里采用條件①、條件②和條件③推導孔邊裂紋的解析權函數。推導Fi(α)所需前提為：同一種參考載荷下的無量綱應力強度因子fr(α)和裂紋嘴張開位移Vr(α)的解。

以裂紋面均布應力(σ(ξ)/σ0=1)作為參考載荷，則有

F1(α)=4fr(α)

(4a)

(4b)

(4c)

式中：

(5)

(6)

(7)

為避免影響后續積分計算，將Φ(α)中的變量α替換為s；ρn和τn分別為fr和Vr的擬合多項式系數?？紤]到當無量綱裂紋長度α>2.0時，孔對裂尖K的影響很小，此時可把孔視為裂紋長度的一部分作為無限板中心裂紋問題處理。因此本文僅考慮α≤2.0范圍的解析權函數。α=0時，邊緣裂紋的無量綱應力強度因子和裂紋嘴張開位移有精確解，分別為[20]：fr(0)=1.121 5，Vr(0)=2.908 6。用有限元法計算在裂紋面均布載荷作用下，α=0.1～3.0范圍內多個裂紋長度的K和CMOD，結合α=0的精確解進行擬合，得到式(6)中fr(α)和式(7)中Vr(α)的多項式系數(如表1所示)，且多項式的最大擬合誤差約為0.1%。

由式(6)可計算出式(5)，然后將式(5)～式(7)代入式(4)求得Fi(α)，進而由式(3)確定系數βi(α)，再由式(2)便能求得無限板孔邊裂紋(N=1，2)的解析權函數m(α,γ)。表2給出了部分裂紋長度的βi(α)值。

為方便工程應用，對由式(3)計算得到的βi(α)值用式(8)作擬合：

(8)

式中:多項式系數bin在表3中給出。

以上多項式的擬合誤差，除βi→0的局部區域外，均不超過1.0%；用式(8)計算的格林函數與式(3)相比，偏差不超過0.3%。

表1 無限板孔邊裂紋無量綱應力強度因子和裂紋嘴張開位移的多項式系數Table 1 Fitted coefficients of dimensionless SIFs and CMODs for radial crack(s) at circular hole in infinite plate

表3 無限板孔邊裂紋解析權函數系數βi(α)的擬合多項式系數Table 3 Fitted polynomial coefficients for βi(α) of radial crack(s) at circular hole in infinite plate

2 無限板孔邊裂紋解析權函數精度檢驗

為保證在任意載荷作用下K和COD的求解精度，必須對第1節得到的孔邊裂紋權函數的精度作檢驗確認。文獻中通常的檢驗方法是：用所求得的權函數計算在某些載荷情況下的K并與文獻中的高精度解比較。盡管這種比較K的方法可以作為評價權函數的一種間接手段，但實際上存在著嚴重的弊端。因為僅對某些載荷情況的K進行比較，并不能真實和全面地反映權函數自身的精度。其主要原因是：① 由于在權函數的推導中對參考載荷情況使用了K自洽條件，故當新的載荷情況與參考載荷沒有明顯差別時，用權函數計算得到的K必然與參考解K的精度很接近(某種形式的“遺傳性”)。但當新的載荷情況與參考載荷有明顯差別時，用權函數計算得到的K就可能出現較大甚至非常大的誤差；②K是由權函數m(α,γ)和裂紋面應力σ(ξ)的乘積沿整個裂紋長度α積分得到的，權函數的正/負誤差可能會由于積分的平均效應相互抵消，從而降低了K的誤差。所以用某些載荷情況下K的良好符合性來證明權函數本身的高精度，這種把多個因素交織在一起的做法，很可能造成假象并引起誤導。

解決權函數精度評價問題的唯一可靠途徑，是對與之相對應的格林函數進行比較[4,7-8,10]。格林函數表示在裂紋面受一對單位集中力(點載荷)P作用下的應力強度因子。利用格林函數評價權函數是對整個裂紋長度α范圍(即0≤γ≤1.0)逐點進行比較，所以能排除權函數對參考載荷情況的偏向性和基于K的積分平均效應。這種基于格林函數的評價方法由于采用逐點檢驗的方式，因此能夠準確反映權函數在裂紋面任意位置的精度。格林函數G(α,γ)與權函數m(α,γ)之間的關系非常簡單[4]：

(9)

2.1 無限板孔邊單裂紋的Shivakumar-Forman格林函數[12]

Shivakumar和Forman采用Muskhelishvili的復變函數理論[21]，得到了無限板孔邊單裂紋(N=1)的格林函數值，并擬合給出了包含30個系數的二重多項式[12]。采用本文的變量定義，其擬合式為

(10)

式中：

多項式C(δ,γ)中的系數Cm,n如表4所示。

表4 式(10)中的系數Cm,nTable 4 Coefficients Cm,n in Eq. (10)

2.2 無限板孔邊雙裂紋的Newman格林函數[13]

Newman用邊界配位法計算得到了無限板孔邊雙裂紋(N=2)封閉形式的格林函數擬合表達式[13]。采用本文的變量定義，其擬合式為

(11)

式中：

2.3 基于復變函數泰勒級數展開理論(WCTSE)的孔邊裂紋權函數

Wagner和Millwater在2012年提出了一種基于泰勒級數展開理論和復變函數有限元分析的高精度數值權函數(WCTSE)法[9]。Jing和Wu對這種方法作了深入研究，進一步提高了其求解精度[10]。對共線中心裂紋和圓盤邊緣裂紋等具有解析解的裂紋幾何的驗證表明，WCTSE法有很高的精度和魯棒性[7-8, 10-11]。

WCTSE法給出的權函數是數值形式的，需要通過擬合才能得到封閉形式的權函數表達式。Jing和Wu[10]的研究表明，對邊緣裂紋使用式(12)形式的擬合多項式能使權函數獲得最佳精度：

M2(1-γ)2+…+Mn(1-γ)n]

(12)

本文利用WCTSE法，對孔邊裂紋(N=1, 2)進行了復變有限元計算和權函數擬合。表5給出了按式(12)擬合的系數Mi。

2.4 無限板孔邊裂紋Wu-Carlsson解析權函數的精度驗證

將式(2)和式(12)分別代入式(9)可得到Wu-Carlsson解析權函數和WCTSE數值權函數的格林函數，并與Shivakumar-Forman(式(10))和Newman(式(11))的結果進行比較。

表5無限板孔邊裂紋(N=1，2)WCTSE法權函數(式(12))的擬合系數Mi

Table5WCTSEcoefficientsMi(Eq.(12))ofradialcrack(s) (N=1,2)atcircularholeininfiniteplate

孔邊單裂紋(N=1)Mii=1i=2i=3i=4α=0.1 0.417 4 0.198 8α=0.2 0.347 8-0.308 9 0.865 0-0.462 7α=0.5 0.048 2-0.093 0 0.442 6-0.247 5α=1.0-0.174 5 0.148 2-0.036 8α=1.5-0.199 0 0.019 4 0.012 4α=2.0-0.230 2-0.030 7 0.024 4孔邊雙裂紋(N=2)Mii=1i=2i=3α=0.10.331 90.394 6-0.117 9α=0.20.292 20.197 5α=0.50.167 40.149 5α=1.00.138 20.117 3α=1.50.144 30.109 3α=2.00.157 20.105 9

2.4.1 孔邊單裂紋

對于孔邊單裂紋(N=1)，本文得到的Wu-Carlsson解析權函數與Shivakumar-Forman[12](式(10))和WCTSE法的結果均符合良好(見圖2(a))。所對應的格林函數與式(10)的差別基本在±2.0%以內(見圖2(b))；與WCTSE法最大差別約為2.0%(見圖2(c))。

這里需要指出，本文的解析權函數(N=1)和WCTSE法在γ=0時的格林函數符合很好(見圖2(c))，而Shivakumar-Forman公式[12]結果比其他2種方法偏高約3%，見表6。由此可以推斷，圖2(b)中在γ=0附近稍大的差別很可能是文獻[12]的精度偏低造成的。

αWu-CarlssonWCTSEDifference/%Shivakumar-Forman[12]Difference/%0.51.629 51.626 8-0.171.676 12.91.01.325 11.325 0-0.011.367 03.21.51.174 91.177 8 0.241.207 32.8

2.4.2 孔邊雙裂紋

對于孔邊雙裂紋(N=2)，本文得到的Wu-Carlsson解析權函數與Newman[13]和WCTSE法的結果均符合良好(見圖3(a))，所對應的格林函數與Newman的式(11)以及WCTSE格林函數的差別都不超過2%(見圖3(b)和圖3(c))。

以上比較表明，本文求得的孔邊單/雙裂紋Wu-Carlsson解析權函數不但具有很高的精度，而且是以無量綱裂紋長度α為變量的封閉解，能夠方便確定任意裂紋長度(0≤α≤2.0)的權函數解，這對疲勞裂紋擴展和壽命預測計算是十分有利的。采用這種統一的權函數解法，僅需變換K和CMOD參考解的多項式擬合系數即可分別求解孔邊單/雙裂紋問題。

最近Jin等[14-15]采用Glinka-Shen的基于2種參考載荷應力強度因子解的方法[6]得到了孔邊裂紋(N=1, 2)數值形式的權函數，進而通過擬合給出了封閉形式的孔邊裂紋權函數表達式。Jin等利用這些權函數計算的幾種其他載荷情況的K與文獻中的解符合良好，進而據此認為其權函數的精度得到了驗證。如前所述，實際上這種基于K的“驗證”很容易誤導。本文作者通過對權函數所對應的格林函數比較，對文獻[14-15]給出的權函數(表格數值和擬合公式)逐點進行了精度檢驗，發現其孔邊單/雙裂紋的權函數均存在顯著誤差：孔邊單裂紋結果與Shivakumar-Forman[12]及本文的解析權函數最大差別達±200%，如圖4所示；孔邊雙裂紋與Newman解[13]的差別雖然明顯小于單裂紋情況，但也達-8%～14%，如圖5所示。圖6則進一步給出了孔壁楔形加載時的應力強度因子差別，可見文獻[15]的孔邊單裂紋結果與其他3種權函數的結果差別極大。本文作者已對Jin等孔邊裂紋權函數[14-15]的精度問題提出質疑[22]。在此提請讀者注意。

3 應用高精度解析權函數求解孔邊裂紋問題

利用本文得到的高精度解析權函數可以高效求解各種復雜載荷情況下的孔邊裂紋問題，下文給出代表性的3種載荷情況的K計算結果。

3.1 裂紋嘴楔形載荷

考慮孔邊裂紋承受楔形載荷(即作用在裂紋嘴處的一對集中力P)，如圖6所示。該集中力可以表示為σ(ξ)/σ0=Pδ(0)，δ(0)是Dirac函數。此時的應力強度因子實際上就是孔邊裂紋在γ=0點的格林函數。用式(2)和式(3)確定的解析權函數的計算結果如圖6所示。Wu-Carlsson權函數結果與Shivakumar-Forman[12](N=1)和Newman[13](N=2)的解，以及與WCTSE結果之間的差別均不超過2%。而Jin權函數結果則明顯偏離，其單裂紋的K在α=0～0.16范圍甚至為負值(見圖6(a))，這顯然有悖于物理規律。

3.2 裂紋面多項式分布載荷

考慮裂紋面承受如圖7所示的冪函數分布應力。該載荷形式可以表示為

σ(ξ)/σ0=ξn

(13)

將式(13)和式(2)分別代入式(1)即可得到孔邊裂紋無量綱應力強度因子的Wu-Carlsson解析權函數解fn。經適當化簡，fn可以寫成包含系數βi(α)的解析表達式[4, 7-8]為

(14)

將式(3)的βi(α)代入式(14)得到孔邊裂紋在冪函數分布應力作用下的應力強度因子fn。為方便使用，表7列出了部分結果。

在很多情況下，連續分布的裂紋面應力可用無量綱坐標ξ的多項式表示。根據疊加原理，若孔邊裂紋裂紋面承受多項式分布載荷為

σ(ξ)/σ0=∑Cnξn

(15)

則對應的應力強度因子可由冪函數應力引起的無量綱應力強度因子fn計算得到：

f=∑Cnfn

(16)

需要注意的是，式(14)中的積分上下限分別為0和α，因此表7中給出的fn僅適用于應力分布在整個裂紋面上的情況(0≤ξ≤α)。

表7 孔邊裂紋在裂紋面冪函數分布載荷(σ(ξ)/σ0=ξn)下的無量綱應力強度因子fnTable 7 Dimensionless SIFs fn of radial crack(s) at circular hole under power stress (σ(ξ)/σ0=ξn) on crack surfaces

3.3 冷擠壓殘余應力場K的Wu-Carlsson權函數法求解

為延長帶圓孔飛機結構的疲勞壽命，常采用過盈配合的銷釘冷擠壓技術(或干涉配合)引入殘余應力。由于其自平衡特性，殘余應力場在近孔壁環形區域內為壓縮，遠離孔邊則為拉伸。文獻[23]針對彈性理想塑性材料，基于塑性不可壓縮和平面應變條件給出了正則化的冷擠壓殘余應力分布：

[0.5+ln(1+ρ)](1+ξ)-2}/[2ln(1+ρ)]

ξ≤ρ

(17a)

ξ>ρ

(17b)

式中：ρ為塑性區圓環的無量綱寬度(以孔半徑R作無量綱處理)，ρ=0.3，0.6，1.0時的裂紋面殘余應力分布如圖8所示。

當裂紋長度α≤ρ時，直接將式(17a)和式(2)代入式(1)即可求出孔邊裂紋在擠壓殘余應力作用下的應力強度因子；當α>ρ時，需要將式(17a)和式(17b)分別代入式(1)，并對應力和權函數的乘積進行分段積分(也可以采用分區段線性化求和的方法[4])：

(18)

在式(17)的擠壓殘余應力下，用式(18)計算得到的孔邊單裂紋的無量綱殘余應力強度因子(k=f·(πα)1/2)如圖9所示。值得注意的是，k(α)曲線隨著裂紋長度的增加趨于零，這是殘余應力場的自平衡性決定的k曲線固有特征[4]。

以上用解析權函數求解冷擠壓孔邊殘余應力強度因子的方法，同樣也可用來高效準確地求解圓孔由于過盈配合和噴丸引起的殘余應力場中的孔邊裂紋問題。

4 結 論

利用3個求解條件推導得到了無限板孔邊單/雙裂紋的Wu-Carlsson解析權函數，并與文獻中采用其他方法以及基于泰勒級數展開和復變函數有限元法(WCTSE)求得的權函數作了廣泛對比驗證。得到以下主要結論：

1) 本文求得的孔邊單/雙裂紋Wu-Carlsson解析權函數具有高精度。其對應的格林函數分別與基于Muskhelishvili復變函數理論的Shivakumar-Forman和邊界配位法的Newman擬合公式，以及WCTSE法的結果高度符合，最大差別在2.0%以內。

2) 裂紋嘴張開位移CMOD條件能夠進一步提高孔邊單/雙裂紋Wu-Carlsson解析權函數的精度，特別是當α值較小時。

3) 利用權函數計算某些載荷情況的K值并據此對權函數精度進行驗證的做法，不僅不能準確全面反映權函數的真實精度，而且會得到誤導的結論?；诟窳趾瘮抵档闹瘘c對比是評價權函數精度的唯一正確途徑。

4) 對裂紋嘴楔形集中力、冪函數分布應力和冷擠壓殘余應力場孔邊裂紋問題的求解表明，本文的孔邊裂紋Wu-Carlsson解析權函數為工程結構的孔邊裂紋問題提供了高效高精度的求解手段。