無限域橫觀各向同性巖體介質(zhì)中折線裂紋邊界元方法

王炳軍，肖洪天，孫凌志，岳中琦

（1.山東科技大學(xué) 土木建筑學(xué)院 山東省土木工程防災(zāi)減災(zāi)重點實驗室，山東 青島 266590；2.香港大學(xué) 土木工程系，香港 999077）

1 引 言

眾所周知，巖體材料中往往含有斷層、節(jié)理、裂隙、孔隙或微裂紋等非連續(xù)面缺陷。在外部荷載作用下，應(yīng)力重分布使巖體中的非連續(xù)面缺陷逐漸增長發(fā)展，進(jìn)而產(chǎn)生宏觀斷裂并產(chǎn)生新的貫通滑移面，從而造成巖體發(fā)生斷裂破壞。在線彈性斷裂力學(xué)中，這些非連續(xù)面缺陷均被模擬為裂紋。應(yīng)力強(qiáng)度因子是用來表征裂紋尖端應(yīng)力場強(qiáng)度的重要物理量，也是判斷裂紋是否擴(kuò)展的一個重要參數(shù)，即線彈性斷裂力學(xué)中的應(yīng)力強(qiáng)度因子準(zhǔn)則。因此，對于裂紋問題的分析，首要任務(wù)是精確計算裂尖處的應(yīng)力強(qiáng)度因子。

對于巖體材料，不論裂紋賦存于巖體內(nèi)部還是局部表面，其本質(zhì)上是三維的，其解析解難以得到，需要借助于數(shù)值方法。三維裂紋問題一直倍受學(xué)術(shù)界和工程界的關(guān)注和重視，人們利用有限元方法[1]、邊界元方法[2]、基金資助：法[3]等研究了裂紋面光滑的三維裂紋問題的應(yīng)力強(qiáng)度因子。

巖體材料在力學(xué)特性上往往呈現(xiàn)為各向異性。由于巖體材料的各向異性，往往造成裂紋在擴(kuò)展過程中裂紋面背離原來的平面，即裂紋面往往不是沿直線擴(kuò)展，從而導(dǎo)致擴(kuò)展后的裂紋變?yōu)榱鸭y面非光滑的折線裂紋。因此，采用數(shù)值方法，精確計算折線裂紋的應(yīng)力強(qiáng)度因子，是巖體介質(zhì)中三維裂紋擴(kuò)展模擬的重要基礎(chǔ)。

目前，關(guān)于折線裂紋問題，國內(nèi)外開展了一些研究。很多學(xué)者采用Muskhelishvili復(fù)勢分析[4－5]、保角映射法[6]、積分變換法[7]、積分方程與有限元的混合法[8]、體積力法[9]、奇異積分方程法[10－12]和有限元法[13]等研究了折線裂紋問題。相對于其他數(shù)值方法，邊界元方法在處理裂紋問題方面有明顯優(yōu)勢，同時邊界元方法特別適合處理無限域和半無限域問題。顯然，巖體力學(xué)問題大多為無限域或半無限域問題。為此，本文在岳中琦等[14]發(fā)展的適用于光滑裂紋問題的雙材料對偶邊界元方法的基礎(chǔ)上，對該方法做了進(jìn)一步發(fā)展，引入了用于離散非光滑邊界的非連續(xù)單元和相應(yīng)的形函數(shù)，建立了用于分析折線裂紋問題的對偶邊界元方法，最后用該方法計算了橫觀各向同性巖體材料中折線裂紋的應(yīng)力強(qiáng)度因子。

2 無限域雙層橫觀各向同性材料的對偶邊界元方法

2.1 基于雙層橫觀各向同性材料基本解的對偶邊界積分方程

岳中琦[15]提出了雙層橫觀各向同性材料的基本解。岳中琦等[14]發(fā)展了基于該基本解的對偶邊界元法。該種類型的對偶邊界元法建立了兩套積分方程，即位移邊界積分方程和面力邊界積分方程。文獻(xiàn)[14]所發(fā)展的邊界元法適用于分析有限域雙層材料中的裂紋問題。當(dāng)然，建議的方法適應(yīng)于分析無限域雙層材料中的裂紋問題。此時，只需面力邊界積分方程，并且退化為

建議的方法適用于分析裂紋面光滑的三維裂紋問題，關(guān)于該方法的公式推導(dǎo)及數(shù)值計算過程參見文獻(xiàn)[2]。本文對該方法做了進(jìn)一步發(fā)展，發(fā)展后的方法適于分析裂紋面非光滑折線裂紋問題。下面僅對本文方法所做的發(fā)展給出說明。

2.2 邊界積分方程的數(shù)值方法

折線裂紋幾何模型如圖1所示。對于其中離開折線EF 的光滑裂紋面，采用原方法的3種離散單元，具體參見文獻(xiàn)[2]。對于靠近折線EF 邊界處的非光滑裂紋面，考慮到裂紋面位于折線上各點的外法線方向不惟一，為滿足被積函數(shù)的連續(xù)性，本文在原方法基礎(chǔ)上，新增加了如圖2所示的3種非連續(xù)單元。3種單元相應(yīng)的形函數(shù)參見文獻(xiàn)[16]。

單元內(nèi)任意點的坐標(biāo)表示為式中：Nα(α=1，9)為單元的插值函數(shù)；(ξ1，ξ2)為單元局部坐標(biāo)系。

如圖2(a)所示，Ⅳ型單元內(nèi)任意點的間斷位移Δui與節(jié)點間位移關(guān)系式為

如圖2(b)、2(c)所示，Ⅴ和Ⅵ型單元任意點的間斷位移Δui與節(jié)點位移關(guān)系式為

圖2 3種類型9節(jié)點非連續(xù)單元Fig.2 Three types of nine-nodes discontinuous element

2.3 應(yīng)力強(qiáng)度因子的計算

在裂尖處節(jié)點建立局部三軸正交坐標(biāo)系（x1，x2，x3），其中x1軸垂直于裂紋面，x2軸為垂直于裂紋面與切向于裂紋面的平面的交線，x3軸與裂尖長度方向相切。節(jié)點間斷位移定義為

根據(jù)文獻(xiàn)[17]，在距離裂尖r 處的節(jié)點間斷位移與裂尖應(yīng)力強(qiáng)度因子的關(guān)系式為

式中：Δ u=(Δ u1，Δ u2，Δu3)T；k=(KI，KII，KIII)T為3種模式的應(yīng)力強(qiáng)度因子；L為Barnett-Lothe張量，依賴于局部坐標(biāo)系下裂尖材料的各向異性特點。

為便于比較，將計算出的應(yīng)力強(qiáng)度因子KI、KII、KIII分別除以含裂紋長為2c 的無限大平板受均勻拉伸荷載p 時的應(yīng)力強(qiáng)度因子，進(jìn)行無量綱化，即

3 方法準(zhǔn)確性驗證

如圖1所示，當(dāng)裂紋彎折角度θ=0°時，折線裂紋即變?yōu)楣饣鸭y，也就是說，光滑裂紋是折線裂紋的特例。本文首先利用修改后的方法計算了無限域橫觀各向同性巖體中矩形光滑裂紋（含折線）的應(yīng)力強(qiáng)度因子，然后與原矩形光滑裂紋（不含折線）的數(shù)值結(jié)果進(jìn)行對比，以驗證修改后的該方法的準(zhǔn)確性。其中，本文采用的橫觀各向同性巖體彈性系數(shù)見表1，矩形裂紋的AB 邊長為2c，BC 邊長為4c，裂紋面作用法向均布壓力q，其中c=1 m，q=p=1 GPa。裂紋面單元網(wǎng)格如圖3所示。應(yīng)用修改后的方法計算得到的最大值SIF=（KI/p (π c)1/2）=0.904 73，原矩形光滑裂紋的最大值SIF=0.904 74[2]，表明了本文方法的正確性。

表1 橫觀各向同性巖體彈性系數(shù)Table 1 Elastic coefficients of transversely isotropic rock material

圖3 折線裂紋面單元網(wǎng)格Fig.3 Meshes of rectangular kinked crack

4 計算結(jié)果分析

下面分析無限域橫觀各向同性巖體（彈性系數(shù)見表1）中的折線裂紋，其幾何模型如圖1所示。其中ABEF 裂紋面平行于橫觀各向同性體的各向同性面，ECDF 裂紋面傾斜于橫觀各向同性體的各向同性面。圖4為設(shè)定的作用于折線裂紋面的3種荷載分布方案。分別為整個裂紋面作用法向均布力、僅ABEF 裂紋面作用法向均布力和僅ECDF 裂紋面作用法向均布力。

圖5為無限域橫觀各向同性介質(zhì)中折線裂紋在方案1即整個裂紋面作用法向均布力時的裂尖應(yīng)力強(qiáng)度因子隨裂紋面彎折角度的變化。從圖中可以看出，在橫觀各向同性介質(zhì)中，裂紋面ABEF 平行于介質(zhì)各向同性面，而裂紋面ECDF 與介質(zhì)各向同性面呈傾斜角度θ，由此彎折型矩形裂紋裂尖應(yīng)力強(qiáng)度因子不再沿折線EF 對稱，同時表明了各向異性對裂尖應(yīng)力強(qiáng)度因子的影響。

圖4 折線裂紋面作用的法向均布力Fig.4 Normal uniform loads acted on surface of kinked crack

圖5 方案1均布荷載作用下折線裂紋應(yīng)力強(qiáng)度因子Fig.5 SIF K1of kinked crack against Barnett-Lothe tensor L for case 1

同時還可知，在整個裂紋面作用法向均布力時，折線矩形裂紋裂尖應(yīng)力強(qiáng)度因子小于光滑的非彎折矩形裂紋。這是因為，對于折線型矩形裂紋，折線兩側(cè)裂紋面外法線方向不同，兩側(cè)裂紋面上的法向力方向自然也不同。因此，在法向均布力作用下，折線兩側(cè)裂紋面ABEF 和ECDF 在張開時存在相互抑制效應(yīng)，即相互限制了裂紋面的張開位移。且這種相互抑制效應(yīng)隨彎折角度θ 的增大而增強(qiáng)。從圖5可知，隨彎折角度θ 的增大，裂尖應(yīng)力強(qiáng)度因子逐漸減小。同時還可發(fā)現(xiàn)，相對于光滑的矩形裂紋，折點E 點處的應(yīng)力強(qiáng)度因子的減小幅度最大。即這種抑制效應(yīng)對裂尖上的折點E 點的影響最大。

圖6 方案2均布荷載作用下折線裂紋應(yīng)力強(qiáng)度因子Fig.6 SIF K1of kinked crack against Barnett-Lothe tensor L for case 2

圖6為無限域橫觀各向同性介質(zhì)中折線裂紋在方案2即僅ABEF 裂紋面作用法向均布力時裂尖應(yīng)力強(qiáng)度因子隨裂紋面彎折角度的變化。從圖可以看出，當(dāng)法向均布力僅作用于平行于各向同性面的裂紋面ABEF 時，隨裂紋面彎折角度的變化，平行于EF 折線的AB 裂尖應(yīng)力強(qiáng)度因子不變，表明裂紋面的彎折對AB 裂尖的斷裂力學(xué)特性無影響；垂直于EF 折線的AF 和BE 裂尖，裂紋面的彎折僅對裂紋面折點E 和F 處裂尖應(yīng)力強(qiáng)度因子產(chǎn)生影響；而對垂直于EF 折線的EC 和DF 裂尖，隨裂紋面彎折角度的增大，裂尖應(yīng)力強(qiáng)度因子減小，當(dāng)裂紋面彎折角度θ=90°時，其應(yīng)力強(qiáng)度因子減小為0。

圖7 方案3均布荷載作用下折線裂紋應(yīng)力強(qiáng)度因子Fig.7 SIF K1of kinked crack against Barnett-Lothe tensor L for case 3

圖7為無限域橫觀各向同性介質(zhì)中折線裂紋在方案3即僅ECDF 裂紋面作用法向均布力時裂尖應(yīng)力強(qiáng)度因子隨裂紋面彎折角度的變化。從圖中可以看出，對于ABEF 裂紋面，裂尖應(yīng)力強(qiáng)度因子隨彎折角度的增大而減小，這與在方案1和2中所體現(xiàn)出的彎折角度對裂尖應(yīng)力強(qiáng)度因子的影響規(guī)律相同；對ECDF 裂紋面，隨彎折角度的增大，相鄰裂尖的應(yīng)力強(qiáng)度因子變化規(guī)律相反。具體為，對于EC和DF 裂尖，隨彎折角度的增大，其應(yīng)力強(qiáng)度因子增大，CD 裂尖則正好相反。這種影響本質(zhì)上是由于裂紋體的各向異性引起的，因為彎折角度即為裂紋面ECDF 傾斜于橫觀各向同性體的各向同性面的角度，裂紋面ECDF 傾斜于各向同性面角度的變化引起裂尖應(yīng)力強(qiáng)度因子的變化。

5 結(jié) 論

（1）在適用于光滑裂紋問題的對偶邊界元方法基礎(chǔ)上，通過引入新的非連續(xù)單元和相應(yīng)的形函數(shù)，建立了用于分析折線裂紋問題的對偶邊界元方法，用數(shù)值算例驗證了本文所建立方法的正確性。

（2）在整個矩形折線裂紋面作用法向均布力時，計算結(jié)果表明，折線兩側(cè)裂紋面在張開時存在較強(qiáng)的相互抑制效應(yīng)，這種抑制效應(yīng)導(dǎo)致折線裂紋的裂尖應(yīng)力強(qiáng)度因子隨裂紋面彎折角度的增大而逐漸減小，且折點處減小最為明顯。

（3）當(dāng)僅在矩形折線裂紋的折線一側(cè)的裂紋面作用法向均布力時，計算結(jié)果表明，折線兩側(cè)裂紋面在張開時的抑制效應(yīng)較弱。

（4）計算結(jié)果還表明，當(dāng)裂紋面作用法向均布力時，橫觀各向同性巖體中的矩形裂紋在發(fā)生彎折時，不僅折線兩側(cè)裂紋面在張開時存在抑制效應(yīng)，同時裂紋面在彎折時，其斷裂特性還受到巖體各向異性的影響，且各向異性對傾斜于各向同性面的裂紋面的影響最為明顯。

[1]卿海，楊衛(wèi)，呂堅.多條表面裂紋相互作用的應(yīng)力分析[J].工程力學(xué)，2009，26(1):1－7.QING Hai，YANG Wei，Lü Jian.Stress analysis of multiple surface cracks with interactions[J].Engineering Mechanics，2009，26(1):1－7.

[2]王炳軍，肖洪天，岳中琦.半無限橫觀各向同性介質(zhì)中多裂紋相互作用分析[J].巖土力學(xué)，2012，33(8):2527－2535.WANG Bing-jun，XIAO Hong-tian，YUE Zhong-qi.Interaction between two rectangular cracks in a transversely isotropic medium of semi-infinite extent[J].Rock and Soil Mechanics，2012，33(8):2527－2535.

[3]肖洪天，王炳軍，岳中琦.非均勻荷載作用下層狀巖體矩形裂隙力學(xué)特性分析[J].工程力學(xué)，2012，29(12):108－113.XIAO Hong-tian，WANG Bing-jun，YUE Zhong-qi.Analysis of rectangular cracks subjected to non-uniform loads in layered rocks[J].Engineering Mechanics，2012，29(12):108－113.

[4]LO K K.Analysis of branched cracks[J].Journal of Applied Mechanics，1978，45(4):797－802.

[5]ISHIHAMA N.Kinked crack in a bending trapezoidal plate[J].International Journal of Fracture，2001，108(2):133－142.

[6]KITAGAW H，YUUKI R，OHIRA T.Crack-morphological aspects in fracture mechanics[J]，Engineering Fracture Mechanics，1975，7(3):515－520.

[7]THEOCARIS P S，MAKRAKIS G N.The kinked crack solved by mellin transform[J].Journal of Elasticity，1986，16(4):393－411.

[8]THEOCARIS P S，TSAMASPHYROS G，THEOTOKOGLOU E E.An alternating coupling of finite elements and singular integral equations for the solution of branched cracks in finite sheets[J].Engineering Fracture Mechanics，1984，20(4):583－589.

[9]張潤梓，管大椿.拐折裂紋應(yīng)力強(qiáng)度因子的體積力法研究[J].北京工業(yè)大學(xué)學(xué)報，1998，24(2):39－44.ZHANG Run-zi，GUAN Da-chun.A study of the stress intensity factors of the kinked cracks by the body force method(BFM)[J].Journal of Beijing Polytechnic University，1998，24(2):39－44.

[10]SUR U，ALTIERO N J.An alternative integral equation approach for curved and kinked cracks[J].International Journal of Fracture，1988，38(1):25－41.

[11]MELIN S.On singular integral equations for kinked cracks[J].International Journal of Fracture，1986，30(1):57－65.

[12]BLANCO C，MARTLNEZ-ESNAOLA J M，ATKINSON C.Kinked cracks in anisotropic elastic materials[J].International Journal of Fracture，1998，93(3):387－407.

[13]管大椿.三點彎曲梁中折裂紋應(yīng)力強(qiáng)度因子的有限元分析[J].北京工業(yè)大學(xué)學(xué)報，1994，20(4):19－23.GUAN Da-chun.The finite element analysis of the SIFs of off-centered kinked cracks in three-point bending beams[J].Journal of Beijing Polytechnic University，1994，20(4):19－23.

[14]YUE Z Q，XIAO H T，PAN E.Stress intensity factors of square crack inclined to interface of transversely isotropic bi-material[J].Engineering Analysis with Boundary Elements，2007，31(1):50－56.

[15]YUE Z Q.Elastic fields in two joined transversely isotropic solids due to concentrated forces[J].International Journal of Engineering Science，1995，33(3):351－369.

[16]PAN E，YUAN F G.Boundary element analysis of three dimensional crack in anisotropic solids[J].International Journal for Numerical Methods in Engineering，2000，48(2):211－237.

[17]TING T C T.Anisotropic elasticity:Theory and applications[M].New York:Oxford University Press，1996.