Ⅰ-Ⅱ混合型應(yīng)力強(qiáng)度因子分析的三角形 Williams單元★

徐 華 劉文祥

(廣西大學(xué)土木建筑工程學(xué)院 工程防災(zāi)與結(jié)構(gòu)安全教育部重點實驗室，廣西南寧 530004)

工程結(jié)構(gòu)中不可避免地含有類似裂紋的缺陷存在，并且以Ⅰ-Ⅱ混合型裂紋最為常見。為了分析這些裂紋缺陷對工程結(jié)構(gòu)安全性的影響，通常將應(yīng)力強(qiáng)度因子作為重要的研究參量。因其與裂紋體幾何構(gòu)形、外荷載等因素密切相關(guān)，難以直接測得，往往采用數(shù)值方法獲得。

隨著大型計算機(jī)計算能力的增強(qiáng)，有限元方法得到了廣泛的應(yīng)用，并且在工程結(jié)構(gòu)分析中具有較強(qiáng)的適應(yīng)性。早在19世紀(jì)60年代，Swedllow和Kabayashi等利用常規(guī)有限元獲得裂尖應(yīng)力強(qiáng)度因子，但該方法需加密裂尖區(qū)離散網(wǎng)格，即需大量的自由度，才能使應(yīng)力強(qiáng)度因子達(dá)到一定的精度［1］，為了減少計算工作量，Levy，Tracey，Hellen 和 Akin 等［2］分別提出了幾類代表性的奇異元。由于這些奇異元存在缺乏剛性位移與常規(guī)元不易協(xié)調(diào)等不足之處，后來 Henshell［3］和 Barsoum［4］分別獨立提出了 1/4 邊中結(jié)點奇異等參元克服了這一缺陷，這種單元網(wǎng)格可以劃分的相當(dāng)粗，就能滿足工程精度要求，但計算結(jié)果受人為主觀因素影響較大。Williams［5］以裂尖為極點建立局部坐標(biāo)系，推導(dǎo)了以級數(shù)形式表示的位移場和應(yīng)力場，這是斷裂力學(xué)中的第一個局部解，基于此，Cheung［6］和 Leung［7］分別利用 Williams 級數(shù)建立了有限條模型和二層分型有限元模型進(jìn)行結(jié)構(gòu)斷裂分析，楊綠峰等［8，9］建立了斷裂分析的Williams單元。

文獻(xiàn)［8，9］中忽略了裂尖微三角的影響，且只對純Ⅰ型和Ⅱ型問題進(jìn)行了研究，本文在此基礎(chǔ)上，建立了Ⅰ-Ⅱ混合型裂紋裂尖應(yīng)力強(qiáng)度因子分析的三角形Williams單元，該單元很好的克服了裂尖奇異性問題，且單元的位移模型中含有與應(yīng)力強(qiáng)度因子直接相關(guān)的參數(shù)，可以直接確定應(yīng)力強(qiáng)度因子。

1 裂尖奇異區(qū)離散方式和奇異區(qū)位移場

對于含裂紋的彈性體，圍繞裂尖O確定應(yīng)力奇異域(如圖1所示)，其外圍是常規(guī)區(qū)域，利用普通有限元法進(jìn)行離散、分析。以裂尖O為中心用一組射線沿轉(zhuǎn)角方向?qū)⑵娈愑螂x散為多個扇形條元。任取其中的一個扇形條元OAB作為典型單元，說明奇異域單元分析過程:首先利用一組環(huán)繞裂尖O且相互平行的折線(記為Γ1，Γ2，…，Γn)，進(jìn)一步將扇形條元離散為一組幾何形狀相似的梯形單元和一個三角形微子域OCD，其中任意兩相鄰折線Γi+1，Γi到裂尖 O 的距離之比為常數(shù) α，α∈(0，1)，分別稱 n，α 為奇異區(qū)徑向離散數(shù)和徑向離散比例因子。

圖1 裂尖附近區(qū)域離散單元示意圖

利用Williams級數(shù)，建立Ⅰ-Ⅱ混合型裂紋裂尖區(qū)整體位移場:

其中，u0，v0為裂尖位移;ai，bi中 i=1，2，…;m 為待定系數(shù)，且a1，b1分別與應(yīng)力強(qiáng)度因子KⅠ和KⅡ有關(guān);對于平面應(yīng)變問題，式中系數(shù)κ=3－4μ;對于平面應(yīng)力問題，系數(shù)κ=(3－μ)/(1+μ)。

根據(jù)式(1)定義的位移場，可以在求得位移模型中的待定參數(shù)ai，bi后，利用系數(shù)a1，b1得到裂尖應(yīng)力強(qiáng)度因子:

裂尖奇異區(qū)典型單元OAB劃分的n+1個微單元包括:裂尖三角形微單元ODC、最外層的過渡微單元ABFE以及梯形CDEF中的n－1個梯形微單元，如圖1所示。隨著n的增大，直線CD逐漸逼近裂尖O，但不能到達(dá)。因此文中將考慮裂尖微三角ODC對總剛的貢獻(xiàn)。

梯形CDEF包含有n－1個自相似的梯形微單元，對于其中的任意微單元e，根據(jù)普通有限元理論建立單元位移場{u}e，或稱為局部位移場，且有:

由于局部位移場必須與總體位移場相協(xié)調(diào)，所以可以將微單元e的全部節(jié)點(本文采用八節(jié)點四邊形單元)坐標(biāo)依次代入單元OAB的總體位移場表達(dá)式(1)中，求得:

式(4)將微單元位移場的節(jié)點位移用總體位移場的待定參數(shù)表示出來，矩陣［T］e稱為位移場轉(zhuǎn)換矩陣，且有:

這里，(ri，θi)為微單元 e的節(jié)點 i(i=1，2，…，8)在坐標(biāo)系Oxy中的極坐標(biāo)。

將式(4)代入式(3)，可得:

式(6)表示的單元位移場中節(jié)點參數(shù){c}不局限于特定的物理意義。

2 奇異區(qū)三角形Williams單元剛度方程

1)梯形CDAB的剛度方程。梯形CDAB包括過渡單元ABFE和梯形條元集合EFCD，分別形成兩塊區(qū)域的控制方程后疊加，即可得到該區(qū)域的剛度方程，具體推導(dǎo)過程可仿照文獻(xiàn)［8，9］，這里不再贅述。現(xiàn)考慮裂紋面不受外力的情況，形成梯形CDAB的剛度方程為:

其中，分塊剛度矩陣的下標(biāo)r，s分別表示位于常規(guī)區(qū)和奇異區(qū)的節(jié)點;上劃線表示常規(guī)區(qū)單元剛度分塊矩陣。

2)微單元OCD的廣義參數(shù)剛度方程。

三角形OAB的第n+1層單元緊鄰裂尖O，建立六節(jié)點三角形廣義參數(shù)單元位移場:

其中，［T(n+1)］為三角形單元轉(zhuǎn)換矩陣，且有:

其中，(ri，θi)，i=1，2，…，6，表示三角形微單元六個節(jié)點的極坐標(biāo)。

所以，三角形單元的廣義參數(shù)剛度方程為:

其中:

集合梯形單元CDAB和裂尖三角形ODC的剛度方程，容易形成OAB的廣義參數(shù)剛度方程，可以分塊表示為:

式(12)即為三角形Williams單元OAB的剛度方程。集合求解域內(nèi)全部離散單元的剛度方程可以形成結(jié)構(gòu)總體剛度方程。在該方程中引入邊界條件并求解，即可根據(jù)式(2)直接求得應(yīng)力強(qiáng)度因子KⅠ和KⅡ。

3 算例分析

帶有一條中心斜裂紋的彈性有限寬板，一端固定，其余端自由，其中裂紋長度 2a=0.10 m，板寬2w=1.0 m，板高 h=1.25 m，裂紋傾斜角為γ，在y軸方向受到均布拉力σ=30.0 kN/m。彈性模量 E=2 ×105MPa，泊松比 μ =0.167，板厚 t=0.01 m，如圖 2 所示。利用本文建立的混合型裂紋分析的廣義參數(shù)三角形Williams單元求解裂紋尖端處KⅠ和KⅡ，并與文獻(xiàn)［2］給出的解析解作比較。

圖2 帶中心斜裂紋彈性板

其中，F(xiàn)Ⅰ和FⅡ是關(guān)于裂紋傾角及a/w的函數(shù)，可以通過邊界配置法求得或者直接通過文獻(xiàn)［2］第348～349頁查得。這里三角形 Williams單元的重要參數(shù)取值［8，9］分別為:m=10，α =0.9，n=300。

表1列出了傾斜角 γ 分別取30°，45°，60°和75°時應(yīng)力強(qiáng)度因子KⅠ和KⅡ的計算值，并與文獻(xiàn)［2］的結(jié)果進(jìn)行比較。

表1 應(yīng)力強(qiáng)度因子與裂紋傾斜角的關(guān)系

表1列出了本文方法和解析法計算應(yīng)力強(qiáng)度因子KⅠ和KⅡ的結(jié)果，可以看出，兩種方法的結(jié)果基本吻合，最大誤差不超過0.8%，證明三角形Williams單元具有很高的計算精度。另外，圖3給出了三角形Williams單元所得應(yīng)力強(qiáng)度因子與傾斜角γ的關(guān)系曲線，由圖可見，裂紋傾斜角γ在第一象限內(nèi)時，KⅠ和KⅡ隨γ的增加而呈現(xiàn)不同的變化規(guī)律，前者逐漸減小，后者表現(xiàn)出先增后減的規(guī)律。

圖3 KⅠ和KⅡ隨著γ變化

4 結(jié)語

本文利用Williams級數(shù)，研究建立了Ⅰ-Ⅱ混合型裂紋應(yīng)力強(qiáng)度因子分析的廣義參數(shù)三角形William單元。算例分析表明:對于含中心裂紋單向拉伸板，三角形Williams單元計算結(jié)果與應(yīng)力強(qiáng)度因子手冊中提供的解析解吻合較好，隨著傾斜角γ的增加，KⅠ逐漸減小，KⅡ表現(xiàn)為先增后減，與實際相符。三角形Willams單元原理簡單，便于構(gòu)造和應(yīng)用，具有較高的計算精度和計算效率。

［1］黃觀鴻，杜家吉.混合型應(yīng)力強(qiáng)度因子KⅠ、KⅡ的有限元計算［J］.天津大學(xué)學(xué)報，1982(1):25-37.

［2］中國航空研究院.應(yīng)力強(qiáng)度因子手冊(增訂版)［M］.北京:科學(xué)出版社，1993.

［3］R.D.Henshell，K.G.Shaw.Crack tip finite elements are unnecessary［J］.International Journal for Numerical Methods in Engineering，1975(9):495-507.

［4］R.S.Barsoum.On the use of isoparametric finite elements in linear fracture mechanics［J］.International Journal for Numerical Methods in Engineering，1976(10):25-37.

［5］M.L.Williams.On the stress distribution at the base of a stationary crack［J］.Journal of Applied Mechanics，1957(24):109-114.

［6］C.W.Woo，Y.H.Wang，Y.K.Cheung.The mixed mode problems for the cracks emanating from a circular hole in a finite plate［J］.Engineering Fracture Mechanics，1989，32(2):279-288.

［7］A.Y.T.Leung，R.K.L.Su.Mixed-mode two-dimensional crack problem by fractal two level finite element method［J］.Engineering Fracture Mechanics，1995，51(6):889-895.

［8］楊綠峰，徐 華，李 冉，等.廣義參數(shù)有限元法計算應(yīng)力強(qiáng)度因子［J］.工程力學(xué)，2009，26(3):48-54.

［9］楊綠峰，徐 華，彭 俚，等.斷裂問題分析的Williams廣義參數(shù)單元［J］.計算力學(xué)學(xué)報，2009，26(1):33-39.