交叉裂紋各裂尖應(yīng)力強度因子同時快速求解的Williams單元

韓林君楊綠峰4

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

0 引言

許多結(jié)構(gòu)或結(jié)構(gòu)構(gòu)件，諸如航空航天工程、道路工程和房屋建筑工程等普遍存在裂紋缺陷，在外荷載作用下，裂紋極易擴展形成交叉裂紋。裂紋尖端應(yīng)力—應(yīng)變場的強弱程度通常采用應(yīng)力強度因子(Stress Intensity Factors，簡記為SIFs)來表征。傳統(tǒng)的有限元法計算裂尖SIFs時需在裂尖人為確定奇異區(qū)尺寸，然后利用線性外推獲得裂尖SIFs，然而多條裂紋交叉使得裂尖SIFs不再相互獨立，將受各裂尖奇異區(qū)尺寸和間距等因素影響，很難同時快速準確獲得各裂尖SIFs。因此，建立能夠同時快速準確求解裂尖具有相關(guān)性的交叉裂紋各裂尖SIFs的算法，對含交叉裂紋結(jié)構(gòu)安全評價具有重要意義。

外荷載作用下交叉裂紋受力較為復(fù)雜，目前主要采用數(shù)值算法進行模擬。章青等[1]提出了廣義擴展有限元法，并基于此方法研究了交叉裂紋裂尖SIFs及其擴展路徑；石路楊等[2]建立了求解多裂紋擴展的擴展有限元法，結(jié)合算例計算了含交叉裂紋體裂尖SIFs，并模擬了多裂紋擴展及交叉匯合過程；DUAX等[3]采用擴展有限元法研究了交叉裂紋裂尖SIFs，證明了該方法能夠有效地解決交叉裂紋問題。然而采用擴展有限元法解決交叉裂紋問題時，需要在位移場中引入由交叉點產(chǎn)生的附加位移項，使位移場形式更為復(fù)雜。文獻總結(jié)發(fā)現(xiàn)，簡化為十字交叉的裂紋也備受關(guān)注。BROCK等[4]將十字形臂視為位錯陣列，并將十字交叉裂紋受雙向拉伸載荷問題簡化為一對耦合的積分方程，進而研究其動態(tài)擴展；MOUSKOS等[5]采用解析函數(shù)理論和Busemann-Chaplygin技術(shù)，研究了平面應(yīng)變情況下十字交叉裂紋的擴展問題；郭懷民等[6-7]采用復(fù)變函數(shù)法研究了不對稱橢圓孔邊裂紋問題和星形裂紋問題，通過改變幾何參數(shù)，將問題蛻化為十字交叉裂紋問題，并給出了I型和II型裂尖SIFs的解析解；余敏等[8]采用彈性復(fù)勢法研究了壓電材料中非對稱十字交叉裂紋反平面問題，并推導(dǎo)出了裂尖III型廣義強度因子、能量釋放率和應(yīng)變能密度因子表達式；牛燚煒[9]基于廣義粒子動力學(xué)法對含預(yù)制十字裂紋的巖石材料進行了二維和三維的數(shù)值模擬，得到預(yù)制裂紋對巖石材料的破壞機理和力學(xué)性能的影響；CHEN[10]給出了十字交叉裂紋問題的數(shù)值解法，得到了非等長正交十字裂紋裂尖SIFs，并與解析法結(jié)果進行了對比；李浪花[11]采用含十字交叉初始裂紋的試件研究了單向受力初始裂紋的起裂和擴展，并通過ABAQUS軟件進行數(shù)值模擬以驗證實驗結(jié)果。

楊綠峰等[12-15]基于廣義參數(shù)有限元法，研究建立了廣義參數(shù)Williams單元(簡記為W單元)，目前已實現(xiàn)不同斷裂模式下的單條裂紋和相互獨立的多條裂紋裂尖SIFs的快速直接求解。本文結(jié)合廣義參數(shù)有限元思想，研究建立了交叉裂紋各裂尖SIFs同時快速求解的W單元計算格式，因各裂尖之間不再相互獨立，進一步分析了相關(guān)參數(shù)，如薄板尺寸、裂紋夾角與奇異區(qū)尺寸等對交叉裂紋各裂尖SIFs的影響。算例證明本文方法能同時求解交叉裂紋各裂尖SIFs，且具有較高的計算精度。

1 模型的建立及離散

如圖1所示，為含多條裂紋且至少任意兩條裂紋相交的交叉裂紋示意圖，將模型的幾何中心作為坐標原點o，水平方向為x軸，逆時針旋轉(zhuǎn)90°為y軸，建立整體坐標系xoy；以裂尖i為例說明局部坐標系建立規(guī)則：將裂尖作為坐標原點oi，裂紋擴展方向為xi軸，逆時針旋轉(zhuǎn)90°為yi軸，建立局部坐標系xioiyi，圖1所示計算模型中含有z條裂紋，不考慮分支裂紋的情況，共有2z個裂尖，因此需建立2z個局部坐標系。

因裂尖應(yīng)力具有奇異性，故將整個計算模型離散為裂尖奇異區(qū)和外圍常規(guī)區(qū)。外圍常規(guī)區(qū)可利用ANSYS有限元軟件中8節(jié)點等參單元自動離散，減少前處理工作量。裂尖奇異區(qū)離散過程如下：選取以裂尖oi為中心點，邊長為l的正方形區(qū)域為奇異區(qū)，將奇異區(qū)邊界進行8等分離散，并分別與中心點oi連接，如圖2所示。選取條元oiL1L2作為典型單元進行說明，利用一組環(huán)繞裂尖oi且相互平行的折線將條元劃分為n層梯形微單元和紅色裂尖三角形單元，n稱為W單元內(nèi)的徑向離散單元數(shù)。任意相鄰折線Γn、Γn-1到裂尖oi的距離之比為常數(shù)，稱為奇異區(qū)徑向離散比例因子α。W單元最外層微單元L1L2L3L4的邊L1L2為常規(guī)區(qū)與奇異區(qū)的交界，其上3個黑色實心圓形表示常規(guī)單元節(jié)點，另外的5個空心圓形位于奇異區(qū)內(nèi)，為虛節(jié)點，故稱微單元L1L2L3L4為過渡單元。根據(jù)上述離散方式，可將任一裂尖奇異區(qū)等分為8個W單元。

圖1 交叉裂紋示意圖
Fig.1 Diagram of cross cracks

圖2 裂尖奇異區(qū)離散示意圖
Fig.2 Discretization element aroundthe crack tip in singular region

2 交叉裂紋各裂尖SIFs同時快速求解的W單元

2.1 整體位移場與局部位移場

以圖1所示交叉裂紋任意裂尖i為例，裂尖奇異區(qū)整體位移場可用改進的Williams級數(shù)表示，并截取前m+1項，即：

(1)

(2)

式中：Ki,I、Ki,II分別為i裂尖的I型SIFs與II型SIFs。

將式(1)改寫為矩陣相乘形式：

[ω]i=[F(r,θ)]iφi，

(3)

由于φi中的參數(shù)ai,j、bi,j是廣義的，無具體物理意義，故稱之為廣義參數(shù)。裂尖奇異區(qū)的整體位移場由改進的Williams級數(shù)建立，由于裂尖三角形微單元較小，可忽略其對條元oiL1L2的貢獻。故W單元為條元內(nèi)n層互相似的梯形微單元的集合。

以奇異區(qū)內(nèi)的任一梯形微單元為典型單元，根據(jù)有限單元法可得其梯形微單元局部位移場：

(4)

由總體位移場控制局部位移場求出各個結(jié)點的位移表達式，對于圖2所示裂尖i周圍的梯形單元有：

(5)

(6)

由式(4)和式(5)，可得梯形微單元廣義參數(shù)位移場：

(7)

2.2 建立單元剛度方程

由“2.1節(jié)”可知：整個模型離散為2種單元類型，即：奇異區(qū)W單元與外圍常規(guī)區(qū)8節(jié)點等參單元，下面分別介紹這兩種單元剛度形成過程。

2.2.1 奇異區(qū)W單元剛度方程

對于奇異區(qū)內(nèi)任意第k層梯形微單元(2≤k≤n)，根據(jù)變分原理，建立其常規(guī)單元剛度方程：

(8)

則第k層梯形單元的廣義剛度方程可以表示為：

(9)

將條元內(nèi)第2層到第n層單元的剛度進行疊加，得到條元奇異區(qū)剛度方程：

(10)

對于奇異區(qū)最外層的過渡單元，如圖2所示，建立其常規(guī)單元剛度矩陣并分塊表示為：

(11)

進而得到過渡單元廣義剛度方程：

(12)

將式(10)和式(12)疊加，得到條元o1L1L2，即1個W單元的剛度方程：

(13)

根據(jù)式(13)集成單個裂尖8個W單元剛度方程：

(14)

2.2.2 外圍常規(guī)區(qū)8節(jié)點等參單元剛度方程

根據(jù)有限單元法，常規(guī)單元剛度矩陣可由下式得到：

(15)

進行整體坐標與局部坐標轉(zhuǎn)換時，有面積微元公式：

dxdy=|J|dξdη。

(16)

將單元剛度矩陣進行高斯計算可得到：

(17)

式中：Hi、Hj為高斯積分點權(quán)函數(shù)；N為高斯積分點數(shù)；|J|為雅可比行列式。

根據(jù)變分原理，建立上述交叉裂紋模型的常規(guī)區(qū)所有單元的剛度方程并進行集成，可分塊表示為：

(18)

式中：下標r表示常規(guī)區(qū)內(nèi)所有節(jié)點。

2.3 建立整體控制方程

根據(jù)2.2.1節(jié)與2.2.2節(jié)中建立的單個奇異區(qū)所有W單元和外圍常規(guī)區(qū)所有8節(jié)點等參單元剛度方程。將2z個裂尖奇異區(qū)所有W單元剛度方程和常規(guī)區(qū)所有單元剛度方程進行集成，得到整個求解域內(nèi)的整體控制方程：

(19)

引入位移和應(yīng)力邊界條件，對整體控制方程進行求解，即可求得各裂尖奇異區(qū)的廣義參數(shù)列陣φ1,φ2，…，φ2z，進而從各裂尖廣義參數(shù)列陣中提取ai,1、bi,1，根據(jù)式(2)計算得到交叉裂紋各裂尖SIFs。

3 算例分析

圖3 含中心交叉裂紋彈性薄板Fig.3 Elastic thin plate with cross cracks

本算例將圖1所示的多條交叉裂紋模型簡化為如圖3所示的中心十字交叉裂紋矩形薄板。板寬為2W，高為2H，板厚為單位厚度，以十字交叉裂紋交點為坐標原點o，水平方向為x軸，逆時針旋轉(zhuǎn)90°為y軸，建立整體坐標系xoy；將裂尖作為坐標原點oi，裂紋擴展方向為xi軸，逆時針旋轉(zhuǎn)90°為yi軸，建立局部坐標系xioiyi(其中i=A,B,C,D，詳見圖3)。沿x軸水平裂紋長度為2a，過中心的斜裂紋長度為(b+c)，其中b=e+d，c=e-d，e和d均為斜裂紋長度參數(shù)，裂紋夾角為γ；材料參數(shù)泊松比為μ，彈性模量為E；板邊受平行于y軸和x軸的均布拉應(yīng)力分別為σ和λσ作用，其中λ為荷載參數(shù)，表示板邊平行于y軸和x軸所承受的拉應(yīng)力比值。

采用本文方法計算十字交叉裂紋各裂尖SIFs時，裂尖奇異區(qū)采用W單元分析，外圍常規(guī)區(qū)利用ANSYS有限元軟件中8節(jié)點等參單元自動離散，并生成節(jié)點和單元信息文本作為本文方法FORTRAN95編程計算的前處理文件。由圖3可知，根據(jù)十字交叉裂紋中兩條裂紋的位置關(guān)系，可將其分為正交十字裂紋和斜十字裂紋。下面將以這兩種情況進行研究：

3.1 正交十字裂紋

根據(jù)Ⅱ型裂紋的定義可知：因計算模型存在對稱性，當(dāng)承受對稱荷載時，裂紋上下表面未發(fā)生相對滑移，故在本節(jié)研究中所有裂尖僅存在I型SIFs。

3.1.1 等長正交十字裂紋

取裂紋半長a=b=c=10 cm，裂紋夾角γ=90°，板邊受平行于y軸和x軸的均布拉應(yīng)力分別為σ和λσ作用，其中σ=1.0 kN/cm。由于模型與外部載荷均具有對稱性，則有KI,A=KI,C，KI,B=KI,D。將本文方法計算結(jié)果與應(yīng)力強度因子手冊[16]解(見應(yīng)力強度因子手冊增訂版95頁)：

(20)

作比較。

通過設(shè)置不同的W/a分析板的尺寸效應(yīng)對裂尖SIFs的影響規(guī)律。由式(20)可知：裂尖SIFs與裂紋半長a、荷載參數(shù)λ及外荷載σ相關(guān)，故同時討論了無限大板界限與荷載參數(shù)λ之間的關(guān)系。由圖4可知，無論荷載參數(shù)λ如何取值，均在滿足W/a≥13時，KI,A與KI,B皆穩(wěn)定收斂于應(yīng)力強度因子手冊解，可以忽略板的尺寸對裂尖SIFs的影響，即可視為無限大板。當(dāng)λ=0.5時，隨著W/a的逐漸增大，KI,A呈現(xiàn)減小趨勢，KI,B則呈現(xiàn)增大趨勢；當(dāng)λ=1.0時，板承受四邊均布拉應(yīng)力，即σ=λσ，所以KI,A=KI,B。隨著W/a的逐漸增大，KI,A和KI,B均呈現(xiàn)減小趨勢；當(dāng)λ=2.0時，隨著W/a的逐漸增大，KI,A呈現(xiàn)增大趨勢，KI,B則呈現(xiàn)減小趨勢。

同時以W/a=13為例，研究荷載參數(shù)λ與SIFs的關(guān)系，將本文方法計算結(jié)果與應(yīng)力強度因子手冊解進行對比分析，計算結(jié)果高度吻合，相對誤差極小，如圖5所示。隨著荷載參數(shù)λ的逐漸增大，KI,A逐漸減小，KI,B逐漸增大，且呈線性規(guī)律變化，與應(yīng)力強度因子手冊解變化一致，驗證了本文方法的正確性及高精性。

圖4W/a與λ的關(guān)系
Fig.4 Relationship betweenW/aandλ

圖5KI與λ的關(guān)系
Fig.5 Relationship betweenKIandλ

3.1.2 非等長正交十字裂紋

取裂紋半長a=15 cm，且有a≠b，d=0，通過改變裂紋長度參數(shù)e來改變裂紋半長b和c，進而改變式(21)中參數(shù)t的取值。選取板尺寸W=200 cm，H=13b，此時可將模型看作無限大板，取板厚為單位厚度，裂紋夾角γ=90°，板四周承受均布拉應(yīng)力σ=λσ=1.0 kN/cm。將本文方法計算結(jié)果與應(yīng)力強度因子手冊解進行對比。其中應(yīng)力強度因子手冊解(見應(yīng)力強度因子手冊增訂版95頁)為：

圖6 Origin擬合曲線Fig.6 Origin fitting curve

(21)

式中：t=a/(a+b)；F(t)由應(yīng)力強度因子手冊查得。

將本文方法計算非等長正交十字裂紋裂尖SIFs的結(jié)果采用Origin軟件進行多因素擬合分析，擬合曲線如圖6所示，擬合曲線參數(shù)見圖6中的表格，擬合得到計算表達式如下：

F(t)=5.595 22t4-12.026 71t3+
6.941 41t2+0.488 14t+0.011 78，

(22)

式中：t=a/(a+b)。

3.2 斜十字裂紋

斜十字裂紋作為十字交叉裂紋更為普遍的情況，研究其裂紋夾角、奇異區(qū)尺寸等對各裂尖SIFs的影響具有重要意義。以圖3為基本計算模型，取裂紋長度b=e+d，c=e-d，其中e=a，裂紋夾角為γ，受雙向均布拉應(yīng)力σ=λσ=1.0 kN/cm。采用本文方法計算斜十字裂紋各裂尖SIFs，并將計算結(jié)果與應(yīng)力強度因子手冊解或ANSYS解進行對比分析。

根據(jù)小節(jié)3.1中對含垂直交叉裂紋薄板無限大板界限的判定，本節(jié)取板尺寸W=H=350 cm，板厚為單位厚度，水平裂紋半長a=20 cm，裂紋夾角γ=45°。分別取不同d值，討論d/a對各裂尖SIFs的影響，對應(yīng)的應(yīng)力強度因子手冊解為[16]：

(23)

式中：F由應(yīng)力強度因子由手冊查得。

本文方法計算結(jié)果如圖7所示，應(yīng)力強度因子手冊解可由式(23)計算得到，將本文方法計算結(jié)果與應(yīng)力強度因子手冊解進行對比，對于I型和II型SIFs，相對誤差均在5 %以內(nèi)，需要說明：該應(yīng)力強度因子手冊解的修正系數(shù)由應(yīng)力強度因子手冊中曲線圖進行取值，即本算例中應(yīng)力強度因子手冊解存在著人為誤差。由圖7可知，隨著d/a的增大，KI,B、KI,C、KII,A和KII,B呈現(xiàn)逐漸增大的趨勢，KI,A、KI,D和KII,C呈現(xiàn)逐漸減小的趨勢，而KII,D則呈現(xiàn)先減小后增大的趨勢。

圖7 各裂尖SIFs與d/a的關(guān)系Fig.7 Relationship between SIFs and d/a at each crack tip

需要注意的是：對比本文解與應(yīng)力強度因子手冊解可知，本文方法的計算結(jié)果與應(yīng)力強度因子手冊(增訂版，97～98頁)中對于C、D裂尖的II型SIFs系數(shù)的判定相悖，將手冊解中C、D裂尖的II型SIFs系數(shù)相互調(diào)換可能為正確解。為此，本文采用ANSYS解進行佐證可得：應(yīng)力強度因子手冊(增訂版，97～98頁)中C、D裂尖的II型SIFs系數(shù)相互調(diào)換方為正確結(jié)論。

此外，本文分析了裂紋長度參數(shù)d=0時，不同裂紋夾角γ對裂尖SIFs的影響。使用本文方法分別計算裂紋夾角γ取15°～90°的裂尖SIFs，并與ANSYS計算結(jié)果進行對比，如圖8所示，最大相對誤差不超過3 %，計算結(jié)果非常吻合，進一步驗證了本文方法具有較高的計算精度。由結(jié)果分析得，對于不同的裂紋夾角γ，KI,A、KI,B大小相等，符號相同，KII,A、KII,B大小相等，符號相反。由于裂尖奇異區(qū)相對位置隨著裂紋夾角γ的變化而變化，取不同裂紋夾角γ與裂尖奇異區(qū)尺寸，研究其對裂尖SIFs的影響，本文方法計算結(jié)果如圖9所示，裂尖SIFs隨著裂紋夾角γ的變化而變化，但奇異區(qū)尺寸對裂尖SIFs幾乎無影響，即采用本文方法計算斜十字裂紋各裂尖SIFs時，計算結(jié)果對奇異區(qū)尺寸不敏感。

圖8 各裂尖SIFs與γ的關(guān)系
Fig.8 Relationship between SIFs andγ

圖9 各裂尖SIFs與l/2的關(guān)系
Fig.9 Relationship between SIFs andl/2

4 結(jié)論

本文建立了交叉裂紋各裂尖SIFs同時快速求解的W單元，以十字交叉裂紋為例，根據(jù)裂紋的相對位置，分別對正交十字裂紋或斜十字裂紋裂尖SIFs相關(guān)參數(shù)進行了敏感性分析。雖然交叉裂紋各裂尖SIFs之間存在相互影響，不再相互獨立，但采用本文方法分析交叉裂紋各裂尖SIFs發(fā)現(xiàn)其計算結(jié)果對裂尖奇異區(qū)尺寸等相關(guān)參數(shù)并不敏感，裂尖仍可按相互獨立的關(guān)系來分析，證明了W單元分析交叉裂紋時仍具有高精高效性。算例分析表明：

① 對于等長正交十字裂紋問題，當(dāng)板寬與裂紋長度滿足W/a≥13時，可忽略板的尺寸效應(yīng)，將模型當(dāng)作無限大板進行計算，且荷載參數(shù)λ對無限大板界限沒有影響；荷載參數(shù)λ和裂尖SIFs成線性變化規(guī)律，隨著荷載參數(shù)λ的逐漸增大，KI,A逐漸減小，KI,B逐漸增大。

② 對于非等長正交十字裂紋問題，根據(jù)應(yīng)力強度因子手冊求解各裂尖SIFs時，相應(yīng)修正系數(shù)可根據(jù)式(22)取值，從而簡化計算并得到較高精度的SIFs解。

③ 對于斜十字裂紋問題，當(dāng)裂紋夾角γ=45°時，若水平裂紋長度a保持不變，隨著斜裂紋長度參數(shù)d的增加，KI,A、KII,C和KI,D逐漸減小，而KII,A、KI,B、KII,B和KI,C逐漸增大，KII,D則出現(xiàn)先減后增的趨勢；應(yīng)力強度因子手冊中，對該模型C、D裂尖的II型SIFs系數(shù)錯誤，二者相互調(diào)換后可得到正確結(jié)果。當(dāng)斜裂紋長度參數(shù)d=0時，取不同裂紋夾角γ，KI,A、KI,B大小相等，符號相同，KII,A、KII,B大小相等，符號相反。隨著裂紋夾角γ的變化，各裂尖奇異區(qū)相對位置也隨之改變，當(dāng)裂尖奇異區(qū)尺寸改變時，對各裂尖SIFs幾乎無影響，表明采用本文方法分析交叉裂紋問題時，各裂尖SIFs對奇異區(qū)尺寸不敏感。