赵永, 杨天鸿, 于庆磊, 张鹏海
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-03-17
基金项目：国家重点基础研究发展计划项目(2016YFC0801602)；国家自然科学基金资助项目(51574059，51574060)。
作者简介：赵永(1991-)，男，山东临沂人，东北大学博士研究生;
杨天鸿(1968-)，男，辽宁抚顺人，东北大学教授，博士生导师。

摘要：对预制裂纹花岗岩试件进行了单轴压缩试验, 得出了不同裂纹倾角下的岩石强度与弹性模量变化规律.建立了不同裂纹倾角下的岩石颗粒流计算模型, 并根据室内试验结果标定了计算模型的细观参数.通过编制FISH语言提取了裂纹扩展过程中的裂纹信息, 并进行了详细分析.研究结果表明:随着预制裂纹倾角的增大, 裂纹起裂应力、起裂角不断增大, 裂纹扩展由最初的翼型裂纹占主导向次生裂纹占主导转变; 通过定义的裂纹聚集系数分析得出随着裂纹倾角的增大, 裂纹聚集系数曲线由光滑变得曲折, 裂纹扩展路径由集中向分散过渡.
关键词：预制裂纹PFC2D单轴压缩裂纹扩展分形分析
Analysis of Pre-existing Crack Propagation Process Based on Discrete Element Method
ZHAO Yong, YANG Tian-hong, YU Qing-lei, ZHANG Peng-hai
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Yong, E-mail: zhaoyongrock@163.com
Abstract: Uniaxial compression experiments were carried out for granite specimen with pre-existing cracks. Variation laws of rock strength and elastic modulus under different crack angles were obtained, and the crack propagation process was recorded. The calculation model for rock particle flow under different flaw inclinations was established and the micro-parameters of calculation model were calibrated based on the laboratory testing results. Crack information during the crack propagation process was extracted by FISH language and was also analyzed in detail. The results show that with the increase of crack inclination, both the crack initiation stress and initiation angle increase continuously. The crack propagation types change from wing-crack-dominated to secondary-crack-dominated. Through the analysis of crack aggregation coefficient, it is found that with the increase of crack angle, the curve of crack aggregation coefficient changes from smooth to tortuous, and the crack propagation path is under the transition from concentrated to decentralized.
Key Words: pre-existing crackPFC2Duniaxial compressivecrack propagationfractal analysis
岩体中广泛存在着规模不同、产状各异的各类裂纹、节理等缺陷, 致使岩体力学性质具有显著的非连续性、各向异性以及非线性等特征.预制裂隙岩石裂纹扩展规律研究是从细观、宏观角度了解和掌握岩石破坏机理的一条重要途径.关于含裂隙岩体的变形、强度及其裂纹扩展方面的研究一直是岩土工程界中的前沿方向^[1].国内外许多学者对此进行了研究, 结论可总结为预制单裂纹的岩石类脆性材料在压剪作用过程中主要产生两类裂纹^[2-7]:翼裂纹(张拉裂纹)和次裂纹(剪切裂纹).Sagong等^[4]、Eberhardt等^[5]通过预制单节理单轴压缩试验得出岩石内部的裂隙分布特性对岩石材料的强度、变形以及破坏模式具有重要的影响.Horii等^[6]、Prudencio等^[7]研究了多裂纹的裂纹贯通机制, 指出宏观的破坏主要与裂纹分布、侧压及参与贯通的裂纹有关.上述研究对裂纹的起裂机制和扩展贯通模式有了深刻的认识.近年来, 利用数值方法研究节理岩体力学性质成为目前较为有效的方法.Tang等^[8]、Wong等^[9]通过离散元或有限元方法模拟了单轴压缩下预制裂纹的力学问题, 研究岩体强度、裂纹扩展等参数与预制裂纹的关系.前述研究多集中于预制裂纹对岩石强度、变形等规律分析, 而对预制裂纹试件在加载过程中裂纹扩展的动态过程、分形特征、裂纹类型的转变等方面并未进行系统研究.
岩石中的原生缺陷使得岩石的变形和破坏过程表现出强烈的复杂性、模糊性, 传统的力学方法难以较好地描述岩石的裂纹扩展规律.本文依靠室内试验和离散元方法, 研究裂纹扩展类型分类并借助自定义的裂纹聚集度参数研究裂纹扩展规律, 分析不同裂纹倾角下裂纹扩展的非线性变化特征.
1 不同裂纹倾角单轴压缩室内试验试验选取芝麻黑花岗岩, 预制裂纹采用高压水射流切割而成, 试件尺寸70 mm×50 mm×150 mm, 以试件中心为参考点, 预制对称的穿透型矩形裂纹, 裂纹尺寸30 mm×2 mm×50 mm.预制裂纹的角度分别为15°, 30°, 45°, 60°和75°, 加工而成的部分试件如图 1所示.单轴压缩试验采用位移控制加载, 加载速率为0.001 mm/s.
图 1(Fig. 1)

图 1 不同裂纹倾角的花岗岩试件Fig.1 Granite specimen with different crack angles

图 2列出了各角度花岗岩在单轴压缩下的抗压强度和弹性模量变化.从图 2中可以看出, 花岗岩在不同裂纹倾角下的抗压强度存在明显差异, 随着预制裂纹倾角的增大, 抗压强度增加, 并且在15°~45°之间抗压强度随倾角增加较为缓慢; 在45°~75°之间增加较快, 增长量显著增大.弹性模量的变化规律与抗压强度相同.随着裂纹倾角向轴向方向靠近, 裂纹开裂、扩展受到阻碍, 岩石强度增大.
图 2(Fig. 2)

图 2 岩石抗压强度与弹性模量随裂纹倾角的变化Fig.2 The compressive strength and elastic modulus vary with different crack angles (a)—抗压强度；(b)—弹性模量.

2 颗粒流计算模型建立与参数标定计算模型尺寸及预制裂纹尺寸与试验试件一致.加载速度为0.001 m/s, 通过记录模拟过程中墙体与颗粒间的接触力、位移求出单轴压缩条件的宏观力学行为, 通过预制测量点来获取横向应变, 进而验证泊松比等参数.以室内试验确定出的宏观力学参数作为标定目标进行“试错”标定^[10], 直到模型的宏观力学响应满足要求为止.表 1为最终确定的细观力学参数.标定试验的结果如图 3所示.
表 1(Table 1)

表 1 细观力学参数Table 1 Micro mechanical parameters 参数类型 细观参数 数值 颗粒基本参数 最小粒径/mm 0.4 粒径比 1.66 密度/(kg·m-3) 2 788 颗粒接触模量/GPa 12.0 颗粒接触刚度比 1.5 摩擦系数 0.78 平行黏结参数 平行黏结模量/GPa 12.0 平行黏结法向强度均值/MPa 55 平行黏结法向强度方差/(MPa)2 8.4 平行黏结切向强度均值/MPa 65 平行黏结切向强度方差/(MPa)2 8.4 平行黏结接触刚度比 1.5

表 1 细观力学参数 Table 1 Micro mechanical parameters

图 3(Fig. 3)

图 3 PFC单轴压缩试验模型及试验结果Fig.3 PFC model and result of uniaxial compression test

从图 3可以看出, 预制裂纹试件在单轴压缩模拟过程中可分为4个阶段:
1) 初始压缩阶段.该阶段轴向应力较小, 发生弹性变形, 未出现岩石微破裂.预制裂纹周边未出现微裂纹, 裂纹周边应力未达到起裂应力值.
2) 裂纹稳定阶段.随着轴向应力的增大, 裂纹尖端应力值达到起裂应力, 尖端裂纹萌生、缓慢扩展, 形成翼型拉裂纹, 裂纹扩展方向向加载方向偏转.
3) 裂纹加速扩展阶段.随着轴向压力的进一步增大, 岩石裂纹数增多, 成核发展, 翼型裂纹逐渐向加载端延伸.裂纹端部开始出现反向倾斜剪切裂纹.
4) 峰后阶段.当轴向应力达到峰值强度后, 裂纹贯通, 翼型裂纹与反向倾斜裂纹迅速扩展, 轴向应力迅速下降, 表现出较强的脆性.
3 预制裂纹扩展的离散元分析3.1 裂纹扩展过程分析图 4给出了基本的裂纹扩展类型, 其中θ为起裂角, d为起裂距离, β为裂纹倾角, 不同的裂纹扩展模式、路径与岩石材料性质有关.
图 4(Fig. 4)

图 4 单轴压缩条件下预制裂纹扩展模式Fig.4 Pre-existing crack propagation patterns under uniaxial compression

通过室内试验过程中的图像记录结合数值模拟过程中的裂纹扩展图像记录, 按照不同性质裂纹的生成顺序, 将裂纹扩展分为3个阶段.表 2中的各阶段的图片为该阶段结束时的图像.
表 2(Table 2)

表 2 裂纹扩展形态描述Table 2 Description of crack growth pattern

表 2 裂纹扩展形态描述 Table 2 Description of crack growth pattern

15°时(峰值应力78.92 MPa):第一阶段翼型拉裂纹在47.10 MPa时开始起裂, 起裂角θ=79°, 起裂距离d=6.94 mm.第二阶段, 预制裂纹两侧出现向下扩展的倾斜次生裂纹, 翼型裂纹稳步扩展.第三阶段, 翼型裂纹扩展缓慢, 右侧倾斜剪切裂纹快速扩展直至试件破坏.
30°时(峰值应力88.40 MPa):第一阶段翼型裂纹在48.40 MPa时开始起裂, 起裂角θ=88°, 起裂距离d=2.57 mm, 裂纹扩展形式同15°, 但裂纹周边出现较多零星微裂纹.第二阶段裂纹右侧尖端出现倾斜次生剪切裂纹, 同时左侧开始萌生次生倾斜裂纹.第三阶段右侧翼型裂纹、左右倾斜次生剪切裂纹同时扩展, 剪切裂纹扩展程度较高.
45°时(峰值应力97.83 MPa):第一阶段翼型裂纹在50.20 MPa时起裂, 起裂角θ=110°, 起裂距离d=2.27 mm.由于起裂应力较高, 在翼型裂纹生成前, 岩石内部弱基质先产生微破裂.进入第二阶段, 裂纹尖端产生次生倾斜裂纹, 左侧倾斜裂纹发展迅速.第三阶段次生倾斜裂纹迅速扩展, 左侧翼型裂纹缓慢扩展, 而右侧翼型裂纹几乎不再扩展.
60°时(峰值应力120.10 MPa):第一阶段预制裂纹尖端起裂时, 岩石其他区域微裂纹数增多, 翼型裂纹起裂应力增大, 为65.20 MPa, 起裂角θ=124°, 起裂距离d=0.44 mm.第二阶段出现准共面次生裂纹.第三阶段左侧出现倾斜次生裂纹迅速扩展, 同时右侧的准共面次生裂纹也迅速扩展.
75°时(峰值应力135.23 MPa):该角度下裂纹扩展较复杂, 第一阶段裂纹起裂应力值高, 为84.5 MPa, 起裂角θ=182°, 起裂距离d=0.在预制裂纹尖端起裂时, 岩石内部出现较多的微破裂, 且倾斜次生裂纹与准共面次生裂纹同时萌生.第二阶段次生倾斜裂纹与准共面次生裂纹共同扩展, 时间较短, 随后进入第三阶段, 准共面次生裂纹迅速发展.
3.2 裂纹扩展路径的非线性分析岩石脆性破裂具有明显的非线性的特点, 从宏观到微观岩石内部都具有非均值的特点, 内部岩性丰富, 孔隙裂隙分布复杂等, 这便造成了岩石在破裂过程中裂纹的扩展路径、扩展规律复杂的特点.
为了深入地了解裂纹扩展规律, 本文定义了裂纹聚集系数C_s, 原理如图 5所示.裂纹越聚集, C_s值越小, 裂纹越分散, C_s值越大.其计算公式为
图 5(Fig. 5)

图 5 裂纹聚集系数示意图Fig.5 Schematic diagram of crack aggregation coefficient

(1)

式中:N为划分的网格数目; M为单网格中微裂纹的数目; d为微裂纹中心到宏观裂纹的垂直距离.
式(1)的物理意义为一定网格数目下裂纹聚集程度的平均值.该公式需作如下假定:单网格中宏观裂纹呈直线分布, 整体试件的宏观裂纹由多个单网格中的宏观裂纹组成; 单网格中单条微裂纹无聚集意义, 因此出现两个微裂纹才开始计算; 除此之外, 需要强调的是每个试件的模拟计算需进行两次, 因为加载初期所产生的微裂纹数目少, 拟合得到的宏观裂纹与实际偏差较大, 随着加载的进行后续微裂纹数目的增多所拟合的宏观裂纹愈加趋近实际宏观裂纹.第一次模拟计算求得模拟结束时所有网格中的宏观裂纹方程, 第二次模拟提取对应网格的宏观裂纹方程, 对整个加载过程中产生的微裂纹求解裂纹聚集系数.
按照上述求解方法对各角度试件加载过程中的裂纹聚集系数进行求解, 结果如图 6所示.
图 6(Fig. 6)

图 6 不同裂纹倾角下的裂纹聚集系数Fig.6 Crack aggregation coefficient under different crack angles (a)—15°；(b)—30°；(c)—45°；(d)—60°；(e)—75°.

从图 6中可以看出, 15°时, 裂纹聚集系数在第一阶段曲线平滑, 裂纹扩展路径集中, 线性扩展.进入第二阶段岩石进入不稳定扩展阶段, 裂纹聚集系数增加.进入第三阶段, 结合表 2中15°试件的破坏图以及图 6a中的裂纹聚集系数变化曲线及拟合曲线的斜率可以看出, 岩石达到损伤应力后, 裂纹聚集、贯通形成宏观裂纹, 裂纹聚集系数波动增大, 急剧增加.裂纹倾角为30°~60°时, 第一阶段裂纹聚集系数波动增强, 扩展路径越加分散, 进入第二、第三阶段后, 裂纹扩展路径分散程度更高.当裂纹倾角为75°时, 裂纹在整个过程扩展时都较波动, 第一阶段较多的裂纹是由于岩石基质强度较低的区域发生破坏, 而不是全部集中于预制裂纹周边, 所以裂纹聚集系数较大; 进入第二阶段后裂纹开始在预制裂纹周边扩展, 扩展路径开始集中, 裂纹聚集系数较小; 进入第三阶段后, 产生大量的微裂纹, 宏观裂纹形成, 聚集系数波动上升.
从图 6a可以看出, 整个聚集系数符合指数分布, 通过对15°~75°的试件按照C_s=A+Be^－ε/t进行拟合, 表 3列出了各角度拟合后的拟合公式中的常数值以及拟合相关系数(75°除外).可以看出，随着角度的增大拟合相关系数降低, 裂纹扩展越不稳定, 裂纹扩展路径越复杂.裂纹开裂的初始聚集系数A值增大, 裂纹扩展越复杂、越分散, 而75°(图 6e)曲线波动过大, 无法利用指数分布进行拟合.从图 6中可以看出, 初始聚集系数与岩石破坏后的聚集系数随着预制裂纹倾角的增大而增大.
表 3(Table 3)

表 3 裂纹聚集系数的拟合常数Table 3 Fitting constants of the crack aggregation coefficient at different crack angles 裂纹倾角/(°) A B t R2 15 2.944E-4 1.094E-10 -2.889E-4 0.980 30 4.624E-4 1.812E-10 -3.328E-4 0.965 45 4.765E-4 1.164E-8 -4.861E-4 0.837 60 5.185E-4 2.690E-17 -1.824E-4 0.788

表 3 裂纹聚集系数的拟合常数 Table 3 Fitting constants of the crack aggregation coefficient at different crack angles

这说明裂纹聚集系数与岩石的强度和模量随预制裂纹角度的变化呈正相关关系.预制裂纹的角度改变了岩石的裂纹扩展, 影响了裂纹的起裂.起裂应力随着预制裂纹倾角的增大而增大, 较大的起裂应力状态下, 岩石整体强度、模量越大, 微裂纹的产生则聚集于岩石中较弱的结构中, 微裂纹越分散聚集系数则越大.
4 结论1) 通过花岗岩不同倾角预制裂纹的单轴压缩试验, 得出:随着裂纹倾角的增大, 花岗岩的抗压强度与弹性模量增加, 裂纹倾角对岩石的宏观力学性质影响显著.
2) 通过对不同预制裂纹倾角下的模型进行数值模拟得出:随着裂纹倾角的增加, 裂纹起裂应力增大、起裂角增大, 翼型裂纹起裂点不断向裂纹尖端移动; 裂纹倾角越大, 裂纹扩展越复杂, 翼型拉破坏裂纹减弱, 次生倾斜剪切破坏增强.
3) 分析加载过程中定义的裂纹聚集系数C_s值得出:裂纹倾角在15°~60°的C_s曲线可由指数函数进行拟合, 拟合的相关系数逐渐降低, C_s曲线由光滑变得曲折, 75°时曲折度最大, 裂纹扩展由低角度时的扩展路径集中到高角度时的路径分散过渡.
参考文献

[1]	赵洪宝, 胡桂林, 李伟, 等. 预制裂隙岩石裂纹扩展规律的研究进展与思考[J].地下空间与工程学报, 2016, 12(sup2): 899–906. ( Zhao Hong-bao, Hu Gui-lin, Li Wei, et al. Research progress and thinking on the crack propagation law of pre-fractured rock[J].Chinese Journal of Underground Space and Engineering, 2016, 12(sup2): 899–906.)
[2]	Yang S Q. Discrete element modeling on fracture coalescence behavior of red sandstone containing two unparallel fissures[J].Engineering Geology, 2014, 178(6): 28–48.
[3]	Cao R H, Cao P, Lin H, et al. Mechanical behavior of brittle rock-like specimens with pre-existing fissures under uniaxial loading:experimental studies and particle mechanics approach[J].Rock Mechanics & Rock Engineering, 2016, 49(3): 763–783.
[4]	Sagong M, Bobet A. Coalescence of multiple flaws in a rock-model material in uniaxial compression[J].International Journal of Rock Mechanics & Mining Sciences, 2002, 39(2): 229–241.
[5]	Eberhardt E, Stead D, Stimpson B, et al. The effect of neighbouring cracks on elliptical crack initiation and propagation in uniaxial and triaxial stress fields[J].Engineering Fracture Mechanics, 1998, 59(2): 103–115.DOI:10.1016/S0013-7944(97)00095-7
[6]	Horii H, Nemat-Nasser S. Brittle failure in compression:splitting, faulting and brittle-ductile transition[J].Philosophical Transactions of the Royal Society of London, 1986, 319(1549): 337–374.DOI:10.1098/rsta.1986.0101
[7]	Prudencio M, Jan M V S. Strength and failure modes of rock mass models with non-persistent joints[J].International Journal of Rock Mechanics & Mining Sciences, 2007, 44(6): 890–902.
[8]	Tang C A, Kou S Q. Crack propagation and coalescence in brittle materials under compression[J].Engineering Fracture Mechanics, 1998, 61(3/4): 311–324.
[9]	Wong L N Y, Li H Q. Numerical study on coalescence of two pre-existing coplanar flaws in rock[J].International Journal of Solids & Structures, 2013, 50(22/23): 3685–3706.
[10]	Itasca Consulting Group. PFC2D users' manual (version 3. 1)[M]. Minnesota: Itasca Consulting Group Inc, 2004.