曲線裂紋裂尖SIFs 等效分析的 廣義參數Williams 單元確定方法

徐 華，鄧 鵬，藍淞耀，劉祖容，楊綠峰,2

(1.廣西大學土木建筑工程學院，工程防災與結構安全教育部重點實驗室，廣西防災減災與工程安全重點實驗室，廣西，南寧 530004； 2.廣西壯族自治區住房和城鄉建設廳，廣西，南寧 530028)

當前線彈性平面斷裂力學研究領域中，主要圍繞直線裂紋的開裂模式、擴展和相關斷裂參數等開展研究。然而，裂紋在擴展過程中裂尖合力方向可能發生變化，導致擴展角變化，即使初始裂紋為直線型，絕大部分經擴展后的裂紋最終都以曲線的型式存在，因而曲線裂紋更符合實際工程。曲線裂紋具有非線性裂紋邊界條件，與直線裂紋的線性裂紋邊界條件相比，其求解難度更高。因而，研究曲線裂紋的相關斷裂參數及開裂機理意義大、難度高。

SIFs 作為表征裂尖附近應力-應變場強弱程度的重要參量，是目前裂紋問題的重點研究對象。近年來，國內外多個課題組對曲線裂紋裂尖SIFs 開展了研究，其分析方法主要歸納為半解析法、解析法和數值法三大類。半解析法在曲線裂紋的分析中有著較廣泛的應用，通常將問題轉化為奇異積分方程進行求解：Nik Long 等[1]提出將復變函數多重彎曲裂紋問題轉化為超奇異積分方程，并通過雙圓弧裂紋裂尖SIFs 求解的算例證明了該方法的正確性；Aridi 等[2]根據該方法求解了直裂紋與曲裂紋相互作用情況下的裂尖SIFs；Elfakhakhre 等[3]采用修正的復勢自由牽引邊界條件，將半平面彈性彎曲裂紋問題轉化為奇異積分方程對裂尖SIFs 進行求解；Monfared 等[4]則基于位錯技術，導出了非均勻平面上多條彎曲裂紋的I 型和II 型SIFs 求解的Cauchy奇異積分方程。在解析法方面，一些研究者結合復變函數和保角映射理論，分析了多種典型的冪函數曲線裂紋，其中：胡元太等[5]利用Stroh 法及映射法研究了沿拋物線分布的各向異性曲線裂紋問題，但只適用于較平坦的拋物線；魏雪霞等[6]依據復變函數和保角映射推導新的保角變換公式，得到了對稱拋物線曲線裂紋問題的I 型SIFs 解析解；郭懷民 等[7]將魏雪霞的方法推廣至3 次～6 次冪函數類對稱曲線裂紋的裂尖SIFs 求解。數值法則是目前曲線裂紋問題分析的主流方法，其優點在于能夠模擬各種復雜邊界，適用范圍廣，當前最常用的是邊界元法和有限元法：王銀邦等[8]采用邊界元法研究了含曲線裂紋圓柱的扭轉問題，同課題組陸孜子等[9]將研究對象擴展至含曲線裂紋的任意柱體；Elangovan等[10]利用有限元法給出了非加勁彎曲板的裂尖SIFs；Judt 等[11]基于有限單元法改進了路徑無關積分法并精確計算出曲線裂紋尖端SIFs；Choi 等[12]提出了曲線裂紋中SIFs 分析的等幾何法，與標準有限元方法相比，相互作用積分域具有更高的應力-應變場連續性。綜上所述，解析法和半解析法均適用范圍較小且理論推導復雜；而數值法具有廣泛的適用性，其中邊界元法和有限元法在曲線裂紋問題的分析中雖然最常見，但二者均需通過復雜的后處理獲取SIFs，由此導致二次誤差。楊綠峰等[13]提出的W 單元以SIFs 為基本未知量進行求解，既繼承有限元法的通用性，又避免后處理引入的人為誤差，且精度高，在斷裂分析中優勢明顯；同課題組徐華等[14]在前期研究成果上提出了W 單元與矩陣壓縮相結合的SIFs 快速解法，該方法能夠提高復雜荷載作用下帶裂紋薄板裂尖SIFs 的求解效率，且有較高的計算精度，但W 單元應用于裂紋問題分析時，裂紋面須滿足σθ=0、τρθ=0(θ=±π)，而曲線裂紋在直接利用其求解時不能滿足該邊界條件。

為克服上述缺陷，本文在裂尖附近截取斜率呈單調變化的曲線微段，以該曲線微段兩端點的距離為半徑，并以裂尖為圓心建立等效區。在截取裂紋微段兩端引切線并交叉成折線段，以該折線段近似代替曲線微段，則該折線與裂尖相連的一段處于θ=±π 上。經此改進，裂尖局部滿足裂紋面自由的邊界條件，即σθ=0、τρθ=0(θ=±π)。本文重點研究裂尖等效區的選取，并通過算例分析證明曲線裂紋尖端局部等效原則的合理性。

1 曲線裂紋尖端等效區計算格式

1.1 裂尖曲線段等效處理

在平面任意位置建立整體坐標系XOY，曲線裂紋可用曲線函數Y=H(X)表示。如圖1 所示，在裂尖處建立局部直角坐標系xoy和極坐標系ρθ：直角坐標系以o為原點，沿o點切線裂尖指向方向為x軸正方向，其逆時針旋轉90°為y軸；極坐標系以o為極點，ox為極軸ρ，逆時針旋轉方向為正方向。局部坐標系x軸與整體坐標系X軸所成角度為γ。設裂尖o在整體坐標系下的坐標為(X0,Y0)，整體坐標系與局部坐標系的坐標轉換關系如下：

圖1 曲線裂紋整體坐標系和裂尖局部坐標系 Fig.1 The global coordinate of curved crack and local coordinate of crack tip

圖2 曲線裂紋裂尖局部等效區 Fig.2 The local equivalent region of curved crack tip

如圖2 所示，在裂尖附近微小區域截取斜率呈單調變化的裂紋曲線微段ΓA，其截取處裂紋上下表面分別用點A1和A2表示。分別在A1、A2處作切線，并與x軸分別交于點B1、B2，令∠A1B1o=∠A2B2o=ψ。取點A1、A2的局部極坐標為(ρA,θA)，線段oB1、oB2的長度為ρB。以ψmin為臨界值，當ψ≥ψmin時，ΓA的曲率足夠小，此時可對ΓA進行等效替換并建立等效區：

1) 以o點為圓心，ρB為半徑建立圓形區域ΩB，稱為等效奇異區，線段oB1、oB2分別為裂紋上下表面。根據角度區間(-π, π)將奇異區均等分為p個相同的扇形條元，如圖2 所示，p=8；

2) 以o點為圓心，ρB為內徑，ρA為外徑，建立圓環區域ΩA，稱為等效常規區，該區域以線段A1B1和A2B2近似等效原裂紋的上下表面，并同樣離散為p個常規單元，除裂紋上下表面連接的2 個單元外，其余p-2 個單元形狀大小均相同；

3) ΩA+ΩB合稱為曲線等效區。

由于等效區尺寸很小，在等效區內對曲線進行等效處理時，因線段A1B1和A2B2在點A1、A2處與裂紋曲線相切，線段oB1、oB2是在點o處與裂紋曲線相切，故用折線A1B1o和A2B2o擬合ΓA可保證較高的擬合度，并能保證等效奇異區裂紋面自由的邊界條件，即σθ=0、τρθ=0(θ=±π)。

根據等效常規區單元劃分規則可知：∠A1oB1應小于扇形條元弧心角，即∠A1oB1＜2π/p。因曲線微段ΓA足夠小，可近似看作圓弧，則有∠oA1B1=∠A1oB1=｜θA-π｜。根據幾何關系，綜合等效區建立規則，得到如下關系式：

綜上所述，曲線微段截斷點A1(A2)可按下列原則確定：

1) 曲線微段ΓA的斜率呈單調變化；

2)θA滿足式(3)的要求。

確定點A1(A2)后，即可通過點A1(A2)與圓心o建立裂尖曲線等效區。根據式(4)，通過點A1(A2)的局部極坐標(ρA,θA)，便可確定等效奇異區ΩB的半徑ρB，進而確定等效常規區和等效奇異區的范圍。

1.2 等效奇異區Williams 級數位移場與SIFs

按W 單元離散規則對ΩB區域進行網格離散，如圖3 所示。以裂尖o為圓心，α(0＜α＜1)為比例分別建立半徑為αρ B,α2ρB, ···,α nρB的n個同心圓，將各條元分為裂尖扇形微單元和n層相似子單元，即W 單元，其最外層單元中的一條邊與常規區相連，稱之過渡單元。

忽略裂尖扇形微單元的剛度貢獻，在W 單元內根據平面斷裂分析的Williams 級數[15]位移場表達式，對其位移場進行改進并截取級數的前m+1項，如式(5)所示：

圖3 裂尖等效奇異區網格離散 Fig.3 Discretization of equivalent singular region at crack tip

式中：u、v分別表示x、y坐標方向的位移分量；ai、bi分別表示待定系數，因無特定物理意義，稱為廣義參數，其值由外荷載與邊界條件確定；

fi,11(θ) ,fi,12(θ) ,fi,21(θ) ,fi,22(θ)為三角函數，詳見文獻[13]。

待定系數a1、b1分別對應位移表達式中的ρ1/2項和應力表達式中的ρ-1/2項，同I、II 型SIFs 有如下關系：

整理式(5)，將其改為矩陣形式如下：

式中：δ表示整體位移場；φ表示廣義參數，且有：

1.3 等效奇異區整體剛度方程集成

根據有限元中8 節點等參數單元理論，利用變分原理建立奇異區條元內任意子單元剛度方程：

根據式(5)可得：

式中：()kT為第k層子單元轉換矩陣，且：

將式(11)代入式(10)，并在等式兩邊左乘T(k)T

可得：

整理后得廣義剛度方程：

式中：K(k)表示第k層子單元廣義剛度矩陣；f(k)表示第k層子單元廣義荷載列陣，因奇異區內不存在外荷載，故f(k)=0。

根據式(13)可知，第2 層到第n層子單元剛度方程存在相同乘子φ，故可疊加如下：

扇形條元最外層的過渡單元，有3 個節點位于等效常規區與等效奇異區的邊界上，剩余5 個節點位于等效奇異區內，根據節點所處區域的不同，將過渡單元剛度方程進行分塊并與式(14)疊加，即可得到整個條元的剛度方程如下：

根據式(15)集成所有條元剛度方程，即等效奇異區剛度方程，分塊表示如下：

2 總剛度方程的集成與SIFs 求解

將所有常規單元剛度集成并區分出等效奇異區邊界上的2p+1 個節點，將剛度方程分塊如下：

疊加式(16)與式(17)得到總剛度方程如下所示：

3 算例分析

例1. 含中心圓弧裂紋薄板，以板中心為原點建立整體坐標系XOY，并在各裂尖分別建立局部坐標系，如圖4 所示。取H和W遠大于圓弧R，即為無限大板，板厚t=1 cm。圓弧半徑為R，弧心角為2β，關于Y軸對稱，彈性模量20.0 GPa，泊松比μ=0.167；該板受雙向拉伸荷載作用， 荷載大小為σ=1 kN/cm2。結合本文等效處理原則和W 單元求解該算例中裂尖SIFs，并與ANSYS 中1/4 奇異單元計算結果和應力強度因子手冊[16]的解析解作比較。W 單元模型的等效奇異區網格離散如圖3 所示，取條元劃分個數p=8；1/4 奇異單元模型的奇異區取周向離散單元數為8，其單元大小和數目與W 單元相同；兩種計算模型的常規區均通過ANSYS 進行網格自由離散，單元數均控制在500 個左右。

本算例中，W 單元的3 個重要參數α，m和n取參考值[13]：α=0.95，m=20，n=300。

由應力強度因子手冊可知，含有中心圓弧裂紋無限大板的SIFs 解析解為：

圖4 含中心圓弧裂紋無限大板 Fig.4 Infiniteplate with central arc crack

1) 在R與β取不同值的情況下，通過改變ψ選取不同大小的等效奇異區，計算A裂尖SIFs，并與解析解作比較，分析ψ對SIFs 的影響。

圖5 和圖6 所示，當R與β取不同值時，本文方法的SIFs 隨ψ變化的計算結果，所有計算結果與解析解的誤差均在5%以內，說明ψ的取值在(13π/18, 35π/36)區間均能得到正確結果。在(13π/18, 5π/6)區間內，隨著ψ增大，KI呈遞減趨勢，KII呈遞增趨勢，且二者均趨近于解析解；在(5π/6, 8π/9)區間內，KI、KII最接近解析解，誤差小于1%；當ψ＞8π/9 時，等效區的尺寸過小，W 單元在奇異區尺 寸過小的情況下誤差會增大。綜上所述，給出建議值ψmin=5π/6，且奇異區尺寸不宜過小，建議等效區截取曲線長度大于原曲線長度的1/5。

圖5 不同β 情況SIFs 隨ψ 變化的計算結果(R=10 cm) Fig.5 Calculation of SIFs with ψ changed under different β conditions (R=10 cm)

圖6 不同R 情況SIFs 隨ψ 變化的計算結果(β=π/2) Fig.6 Calculation of SIFs with ψ changed under different R conditions(β=π/2)

2) 根據建議，取ψ=5π/6，分別計算不同R、β情況下A裂尖的SIFs，并與ANSYS 計算結果和解析解進行比較，計算結果如圖7 和圖8 所示。

本文方法結果與解析解非常吻合，最大誤差為1.32%，而ANSYS 的1/4 奇異單元計算結果最大誤差為2.77%。由圖7 可知，當R=10 cm，β在(π/18, 4π/9)區間遞增時，KI先增大后減小，在β=π/4 時最大，KII呈遞增趨勢；由圖8 可知，當β=π/6 時，隨著R的增大，KI、KII都呈遞增趨勢。

圖7 SIFs 隨β 變化的計算結果(R=10 cm) Fig.7 Calculation of SIFs changing with β (R=10 cm)

圖8 SIFs 隨R 變化的計算結果(β=π/6) Fig.8 Calculation of SIFs changing with R(β=π/6)

圖9 含圓弧分岔裂紋無限大板 Fig.9 Infiniteplate with bifurcate arccrack

例2. 含圓弧分岔裂紋薄板，以板中心為原點，建立整體坐標系XOY，并在各裂尖分別建立局部坐標系，如圖9 所示。取H和W遠大于圓弧R，即為無 限大板，板厚t=1 cm。圓弧半徑為R=2S，S為直線段裂紋半長，弧心角為2β，關于X軸對稱，彈性模量20.0 GPa，泊松比μ=0.167；該板受雙向拉伸荷載作用，荷載大小為σ=1 kN/cm2。結合本文等效處理原則和W 單元求解該算例中裂尖 SIFs，并與ANSYS中1/4奇異單元計算結果和應力強度因子手冊[16]的解析解作比較。W 單元模型的等效奇異區網格離散如圖3 所示，取條元劃分個數p=8；1/4 奇異單元模型的奇異區取周向離散單元數為8，其單元大小和數目與W 單元相同；兩種計算模型的常規區均通過ANSYS 進行網格自由離散，單元數均控制在500 個左右。

本算例中，等效區根據建議取值如下：β在(0, π/2)區間時，ψ=π-β/5；β在[π/2, π]區間時，ψ=π-β/10。W 單元的3 個重要參數α，m和n取參考值[13]：α=0.95，m=20，n=300。

由應力強度因子手冊可知，含圓弧分岔裂紋無限大板的SIFs 解析解為：

取S=10 cm，β=π/12～11π/12 計算各裂尖SIFs隨β改變的變化趨勢，并將W 單元計算結果與ANSYS 計算結果和解析解進行對比。

由于板受雙向拉伸時，A裂尖裂紋上下表面無水平方向的相對位移，故KII始終為0，不進行圖示，其余計算結果如圖10 所示。對比KI,A和KI,B發現：W 單元計算結果和ANSYS 計算結果均與解析解非常吻合，排除值太小導致誤差值偏大的情況，其余值中W 單元最大誤差為2.34%，ANSYS 最大誤差為3.32%；對比KII,B發現：W 單元計算結果比ANSYS 計算結果更吻合解析解，其中W 單元最大誤差為1.89%，ANSYS 最大誤差為3.94%，且在β=π/3～3π/4 區間ANSYS 誤差明顯大于W 單元；隨著β的增大，KI,A先減小后增加，KI,B和KII,B絕對值先增加后減小，KI,A在β=2π/3 時最小，KI,B在β=5π/12 時最大，KII,B絕對值在β=7π/12 時最大。

圖10 SIFs 隨β 變化的計算結果(S=10 cm) Fig.10 Calculation of SIFs changing with β (S=10 cm)

4 結論

本文對曲線裂紋尖端局部區域進行等效處理并建立等效區，從而使曲線裂紋在裂尖等效奇異區內滿足σθ=0、τρθ=0(θ=±π)的邊界條件，結合等效處理原則與W 單元計算了含中心圓弧裂紋無限大板與含圓弧分岔裂紋無限大板的裂尖SIFs，證明了本文方法的合理性和正確性，得到結論如下：

(1)ψ的取值在(13π/18, 35π/36)區間與解析解的誤差均在5%以內，證明了本文對裂尖進行等效處理的合理性。在(13π/18, 5π/6)區間內，隨著ψ增大，KI呈遞減趨勢，KII呈遞增趨勢，且二者均趨近于解析解；在(5π/6, 8π/9)區間內，KI、KII誤差小于1%；當ψ＞8π/9 時，等效區的尺寸過小，W 單元在奇異區尺寸過小的情況下誤差會增大。綜上所述，給出建議值ψmin=5π/6，且奇異區尺寸不宜過小，建議等效區截取曲線長度大于原曲線長度的1/5。

(2) 當奇異區尺寸取建議值時，本文方法在含中心圓弧裂紋無限大板裂尖SIFs 的計算中，計算結果與解析解的最大誤差為1.32%，精度較ANSYS 計算結果更高；在含圓弧分岔裂紋無限大板裂尖SIFs 的計算結果中：W 單元和ANSYS 的KI,A、KI,B計算結果均與解析解非常吻合，除去值太小導致誤差值偏大的情況，其余值中W 單元最大誤差為2.34%，ANSYS 最大誤差為3.32%；W 單元的KII,B計算結果比ANSYS 更吻合解析解，其中W 單元最大誤差為1.89%，ANSYS 最大誤差為3.94%，且在β=π/3～3π/4 區間ANSYS 誤差明顯大于W 單元。算例結果證明，本文方法具有高精性和廣泛的適用性。

(3) 含中心圓弧裂紋無限大板受雙向拉伸時，R不變，隨著β在(π/18, 4π/9)區間遞增，KI先增大后減小，在β=π/4 時最大，KII呈遞增趨勢；當β不變時，隨著R的增大，KI、KII都呈遞增趨勢。

(4) 含圓弧分岔裂紋無限大板受雙向拉伸時，C不變，隨著β在(π/12, 11π/12)區間遞增，KI,A先減小后增加，在β=2π/3 時最小；KI,B和KII,B絕對值先增加后減小，KI,B在β=5π/12 時最大，KII,B絕對值在β=7π/12 時最大。