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

徐 華，鄧 鵬，藍(lán)淞耀，劉祖容，楊綠峰,2

(1.廣西大學(xué)土木建筑工程學(xué)院，工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室，廣西，南寧 530004； 2.廣西壯族自治區(qū)住房和城鄉(xiāng)建設(shè)廳，廣西，南寧 530028)

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

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

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

1 曲線裂紋尖端等效區(qū)計算格式

1.1 裂尖曲線段等效處理

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 等效奇異區(qū)Williams 級數(shù)位移場與SIFs

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

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

圖3 裂尖等效奇異區(qū)網(wǎng)格離散 Fig.3 Discretization of equivalent singular region at crack tip

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

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

待定系數(shù)a1、b1分別對應(yīng)位移表達(dá)式中的ρ1/2項和應(yīng)力表達(dá)式中的ρ-1/2項，同I、II 型SIFs 有如下關(guān)系：

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

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

1.3 等效奇異區(qū)整體剛度方程集成

根據(jù)有限元中8 節(jié)點(diǎn)等參數(shù)單元理論，利用變分原理建立奇異區(qū)條元內(nèi)任意子單元剛度方程：

根據(jù)式(5)可得：

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

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

可得：

整理后得廣義剛度方程：

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

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

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

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

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

將所有常規(guī)單元剛度集成并區(qū)分出等效奇異區(qū)邊界上的2p+1 個節(jié)點(diǎn)，將剛度方程分塊如下：

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

3 算例分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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