无限域横观各向同性岩体介质中折线裂纹边界元方法

王炳军，肖洪天，孙凌志，岳中琦

（1.山东科技大学 土木建筑学院 山东省土木工程防灾减灾重点实验室，山东 青岛 266590；2.香港大学 土木工程系，香港 999077）

1 引 言

众所周知，岩体材料中往往含有断层、节理、裂隙、孔隙或微裂纹等非连续面缺陷。在外部荷载作用下，应力重分布使岩体中的非连续面缺陷逐渐增长发展，进而产生宏观断裂并产生新的贯通滑移面，从而造成岩体发生断裂破坏。在线弹性断裂力学中，这些非连续面缺陷均被模拟为裂纹。应力强度因子是用来表征裂纹尖端应力场强度的重要物理量，也是判断裂纹是否扩展的一个重要参数，即线弹性断裂力学中的应力强度因子准则。因此，对于裂纹问题的分析，首要任务是精确计算裂尖处的应力强度因子。

对于岩体材料，不论裂纹赋存于岩体内部还是局部表面，其本质上是三维的，其解析解难以得到，需要借助于数值方法。三维裂纹问题一直倍受学术界和工程界的关注和重视，人们利用有限元方法[1]、边界元方法[2]、基金资助：法[3]等研究了裂纹面光滑的三维裂纹问题的应力强度因子。

岩体材料在力学特性上往往呈现为各向异性。由于岩体材料的各向异性，往往造成裂纹在扩展过程中裂纹面背离原来的平面，即裂纹面往往不是沿直线扩展，从而导致扩展后的裂纹变为裂纹面非光滑的折线裂纹。因此，采用数值方法，精确计算折线裂纹的应力强度因子，是岩体介质中三维裂纹扩展模拟的重要基础。

目前，关于折线裂纹问题，国内外开展了一些研究。很多学者采用Muskhelishvili复势分析[4－5]、保角映射法[6]、积分变换法[7]、积分方程与有限元的混合法[8]、体积力法[9]、奇异积分方程法[10－12]和有限元法[13]等研究了折线裂纹问题。相对于其他数值方法，边界元方法在处理裂纹问题方面有明显优势，同时边界元方法特别适合处理无限域和半无限域问题。显然，岩体力学问题大多为无限域或半无限域问题。为此，本文在岳中琦等[14]发展的适用于光滑裂纹问题的双材料对偶边界元方法的基础上，对该方法做了进一步发展，引入了用于离散非光滑边界的非连续单元和相应的形函数，建立了用于分析折线裂纹问题的对偶边界元方法，最后用该方法计算了横观各向同性岩体材料中折线裂纹的应力强度因子。

2 无限域双层横观各向同性材料的对偶边界元方法

2.1 基于双层横观各向同性材料基本解的对偶边界积分方程

岳中琦[15]提出了双层横观各向同性材料的基本解。岳中琦等[14]发展了基于该基本解的对偶边界元法。该种类型的对偶边界元法建立了两套积分方程，即位移边界积分方程和面力边界积分方程。文献[14]所发展的边界元法适用于分析有限域双层材料中的裂纹问题。当然，建议的方法适应于分析无限域双层材料中的裂纹问题。此时，只需面力边界积分方程，并且退化为

建议的方法适用于分析裂纹面光滑的三维裂纹问题，关于该方法的公式推导及数值计算过程参见文献[2]。本文对该方法做了进一步发展，发展后的方法适于分析裂纹面非光滑折线裂纹问题。下面仅对本文方法所做的发展给出说明。

2.2 边界积分方程的数值方法

折线裂纹几何模型如图1所示。对于其中离开折线EF 的光滑裂纹面，采用原方法的3种离散单元，具体参见文献[2]。对于靠近折线EF 边界处的非光滑裂纹面，考虑到裂纹面位于折线上各点的外法线方向不惟一，为满足被积函数的连续性，本文在原方法基础上，新增加了如图2所示的3种非连续单元。3种单元相应的形函数参见文献[16]。

单元内任意点的坐标表示为式中：Nα(α=1，9)为单元的插值函数；(ξ1，ξ2)为单元局部坐标系。

如图2(a)所示，Ⅳ型单元内任意点的间断位移Δui与节点间位移关系式为

如图2(b)、2(c)所示，Ⅴ和Ⅵ型单元任意点的间断位移Δui与节点位移关系式为

图2 3种类型9节点非连续单元Fig.2 Three types of nine-nodes discontinuous element

2.3 应力强度因子的计算

在裂尖处节点建立局部三轴正交坐标系（x1，x2，x3），其中x1轴垂直于裂纹面，x2轴为垂直于裂纹面与切向于裂纹面的平面的交线，x3轴与裂尖长度方向相切。节点间断位移定义为

根据文献[17]，在距离裂尖r 处的节点间断位移与裂尖应力强度因子的关系式为

式中：Δ u=(Δ u1，Δ u2，Δu3)T；k=(KI，KII，KIII)T为3种模式的应力强度因子；L为Barnett-Lothe张量，依赖于局部坐标系下裂尖材料的各向异性特点。

为便于比较，将计算出的应力强度因子KI、KII、KIII分别除以含裂纹长为2c 的无限大平板受均匀拉伸荷载p 时的应力强度因子，进行无量纲化，即

3 方法准确性验证

如图1所示，当裂纹弯折角度θ=0°时，折线裂纹即变为光滑裂纹，也就是说，光滑裂纹是折线裂纹的特例。本文首先利用修改后的方法计算了无限域横观各向同性岩体中矩形光滑裂纹（含折线）的应力强度因子，然后与原矩形光滑裂纹（不含折线）的数值结果进行对比，以验证修改后的该方法的准确性。其中，本文采用的横观各向同性岩体弹性系数见表1，矩形裂纹的AB 边长为2c，BC 边长为4c，裂纹面作用法向均布压力q，其中c=1 m，q=p=1 GPa。裂纹面单元网格如图3所示。应用修改后的方法计算得到的最大值SIF=（KI/p (π c)1/2）=0.904 73，原矩形光滑裂纹的最大值SIF=0.904 74[2]，表明了本文方法的正确性。

表1 横观各向同性岩体弹性系数Table 1 Elastic coefficients of transversely isotropic rock material

图3 折线裂纹面单元网格Fig.3 Meshes of rectangular kinked crack

4 计算结果分析

下面分析无限域横观各向同性岩体（弹性系数见表1）中的折线裂纹，其几何模型如图1所示。其中ABEF 裂纹面平行于横观各向同性体的各向同性面，ECDF 裂纹面倾斜于横观各向同性体的各向同性面。图4为设定的作用于折线裂纹面的3种荷载分布方案。分别为整个裂纹面作用法向均布力、仅ABEF 裂纹面作用法向均布力和仅ECDF 裂纹面作用法向均布力。

图5为无限域横观各向同性介质中折线裂纹在方案1即整个裂纹面作用法向均布力时的裂尖应力强度因子随裂纹面弯折角度的变化。从图中可以看出，在横观各向同性介质中，裂纹面ABEF 平行于介质各向同性面，而裂纹面ECDF 与介质各向同性面呈倾斜角度θ，由此弯折型矩形裂纹裂尖应力强度因子不再沿折线EF 对称，同时表明了各向异性对裂尖应力强度因子的影响。

图4 折线裂纹面作用的法向均布力Fig.4 Normal uniform loads acted on surface of kinked crack

图5 方案1均布荷载作用下折线裂纹应力强度因子Fig.5 SIF K1of kinked crack against Barnett-Lothe tensor L for case 1

同时还可知，在整个裂纹面作用法向均布力时，折线矩形裂纹裂尖应力强度因子小于光滑的非弯折矩形裂纹。这是因为，对于折线型矩形裂纹，折线两侧裂纹面外法线方向不同，两侧裂纹面上的法向力方向自然也不同。因此，在法向均布力作用下，折线两侧裂纹面ABEF 和ECDF 在张开时存在相互抑制效应，即相互限制了裂纹面的张开位移。且这种相互抑制效应随弯折角度θ 的增大而增强。从图5可知，随弯折角度θ 的增大，裂尖应力强度因子逐渐减小。同时还可发现，相对于光滑的矩形裂纹，折点E 点处的应力强度因子的减小幅度最大。即这种抑制效应对裂尖上的折点E 点的影响最大。

图6 方案2均布荷载作用下折线裂纹应力强度因子Fig.6 SIF K1of kinked crack against Barnett-Lothe tensor L for case 2

图6为无限域横观各向同性介质中折线裂纹在方案2即仅ABEF 裂纹面作用法向均布力时裂尖应力强度因子随裂纹面弯折角度的变化。从图可以看出，当法向均布力仅作用于平行于各向同性面的裂纹面ABEF 时，随裂纹面弯折角度的变化，平行于EF 折线的AB 裂尖应力强度因子不变，表明裂纹面的弯折对AB 裂尖的断裂力学特性无影响；垂直于EF 折线的AF 和BE 裂尖，裂纹面的弯折仅对裂纹面折点E 和F 处裂尖应力强度因子产生影响；而对垂直于EF 折线的EC 和DF 裂尖，随裂纹面弯折角度的增大，裂尖应力强度因子减小，当裂纹面弯折角度θ=90°时，其应力强度因子减小为0。

图7 方案3均布荷载作用下折线裂纹应力强度因子Fig.7 SIF K1of kinked crack against Barnett-Lothe tensor L for case 3

图7为无限域横观各向同性介质中折线裂纹在方案3即仅ECDF 裂纹面作用法向均布力时裂尖应力强度因子随裂纹面弯折角度的变化。从图中可以看出，对于ABEF 裂纹面，裂尖应力强度因子随弯折角度的增大而减小，这与在方案1和2中所体现出的弯折角度对裂尖应力强度因子的影响规律相同；对ECDF 裂纹面，随弯折角度的增大，相邻裂尖的应力强度因子变化规律相反。具体为，对于EC和DF 裂尖，随弯折角度的增大，其应力强度因子增大，CD 裂尖则正好相反。这种影响本质上是由于裂纹体的各向异性引起的，因为弯折角度即为裂纹面ECDF 倾斜于横观各向同性体的各向同性面的角度，裂纹面ECDF 倾斜于各向同性面角度的变化引起裂尖应力强度因子的变化。

5 结 论

（1）在适用于光滑裂纹问题的对偶边界元方法基础上，通过引入新的非连续单元和相应的形函数，建立了用于分析折线裂纹问题的对偶边界元方法，用数值算例验证了本文所建立方法的正确性。

（2）在整个矩形折线裂纹面作用法向均布力时，计算结果表明，折线两侧裂纹面在张开时存在较强的相互抑制效应，这种抑制效应导致折线裂纹的裂尖应力强度因子随裂纹面弯折角度的增大而逐渐减小，且折点处减小最为明显。

（3）当仅在矩形折线裂纹的折线一侧的裂纹面作用法向均布力时，计算结果表明，折线两侧裂纹面在张开时的抑制效应较弱。

（4）计算结果还表明，当裂纹面作用法向均布力时，横观各向同性岩体中的矩形裂纹在发生弯折时，不仅折线两侧裂纹面在张开时存在抑制效应，同时裂纹面在弯折时，其断裂特性还受到岩体各向异性的影响，且各向异性对倾斜于各向同性面的裂纹面的影响最为明显。

[1]卿海，杨卫，吕坚.多条表面裂纹相互作用的应力分析[J].工程力学，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]王炳军，肖洪天，岳中琦.半无限横观各向同性介质中多裂纹相互作用分析[J].岩土力学，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]肖洪天，王炳军，岳中琦.非均匀荷载作用下层状岩体矩形裂隙力学特性分析[J].工程力学，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]张润梓，管大椿.拐折裂纹应力强度因子的体积力法研究[J].北京工业大学学报，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]管大椿.三点弯曲梁中折裂纹应力强度因子的有限元分析[J].北京工业大学学报，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.