刘京铄，曹 平，范金星，蒲成志，黄雪姣

(1．中南大学 资源与安全工程学院，湖南 长沙 410083；2.湖南水利水电职业技术学院，湖南 长沙 410131)



预制双裂隙类岩材料断裂特性的无网格法数值模拟

刘京铄1,2，曹 平1，范金星2，蒲成志1，黄雪姣1

(1．中南大学 资源与安全工程学院，湖南 长沙 410083；2.湖南水利水电职业技术学院，湖南 长沙 410131)

摘要:在伺服控制单轴加载试验机上对预制两裂隙类岩试件进行压缩试验，观测试件的破坏模式及宏观裂纹的扩展路径。在简要阐述无网格法基本原理的基础上，将基于加权移动最小二乘法的无网格法应用于含裂隙类岩材料断裂特性的数值模拟，对含预制双裂隙类岩试件的断裂特性及其裂纹扩展路径进行数值计算与分析，给出裂尖应力场及位移场的变化关系，并将数值计算所得的裂纹扩展路径与实验结果进行比较，验证该方法的正确性与精确性。

关键词：预制裂隙；类岩材料；单轴压缩；无网格法；数值分析；裂纹扩展

岩体的变形与强度特性不仅取决于岩石本身，而且很大程度上受限于岩体内部随机分布的节理、裂隙等不连续面。裂隙岩体内部往往由于裂纹尖端应力集中而导致裂纹的起裂、扩展、贯通，进而引起岩体的失稳破坏。由于现场原位试验难以操作，现阶段对裂隙岩体破坏过程的研究主要基于室内类岩材料模型试验展开[1]。关宝树等[2]通过预制贯通裂隙模型对含不同裂隙倾角试件的强度进行了试验；Wong等[3-4]和通过预制裂隙材料对预制裂隙贯通机制进行了研究；张平等[5]采用预制闭合裂隙试件研究了动载下裂隙的扩展与贯通过程。在试验和理论分析的基础上，许多学者结合数值分析方法，对裂纹扩展问题展开了研究[6-8]。但目前的数值分析方法主要基于有限单元法，对于裂纹扩展问题，由于裂纹的扩展方向不能事先确定，因此在计算过程中需要不断的重新划分网格以模拟裂纹的扩展过程，这将可能导致有限元网格的严重扭曲并且严重影响解的精度；而目前发展迅速的无网格方法由于不需要单元网格, 点与点之间只在进行插值构造时才临时建立联结关系, 可以很方便地修改节点分布, 故非常适用于研究裂纹扩展等断裂分析领域。因此近年来国内外许多学者，如Belytschko等[9-11]对无网格方法进行了大量的研究并将该方法应用于求解各类岩土工程问题。庞作会等[12]将无网格法应用于求解边坡开挖问题。ZHANG等[13]将无网格思想应用于求解节理岩体问题。寇晓东等[14]将无网格法应用于求解裂纹扩展问题。庞作会等[15]利用无网格方法求解岩土力学中的不连续面问题。许多学者基于无网格法优于其它数值方法而将其对静动态裂缝扩展进行数值模拟并取得了很好的效果；而本文将在基于实验研究的基础上对预制双裂隙类岩体断裂特性与宏观裂纹扩展用无网格法进行数值模拟，通过试验现象与计算结果的对比分析，验证了运用无网格法求解裂缝扩展及断裂问题的精确性和合理性。

1试件制作与试验

1.1试件制备

试件采用内部尺寸 (长×宽×高)为 200 mm×150 mm×30 mm的钢模，用C42.5白水泥与自来水经1.05 mm孔径筛分过的细砂按体积比2∶1∶1在室温下浇筑而成；采用在试件中预埋0.4 mm厚的铝合金薄片的方法预制贯通裂隙，预制裂隙长度为20 mm，将试件在标准条件下养护24 h后拔出铝合金片，并将试件置于静水中养护14 d，待养护结束后检测裂隙闭合情况。

试件中裂隙分布方式如图1所示：α为裂隙倾角;β为岩桥倾角;b为岩桥长度。本试验中α取45°，β取90°，b取40 mm。

图1　裂隙与岩桥分布图Fig.1 Distribution of cracks and rock-bridge

1.2试验仪器与加载条件

本次试验在高精度微机控制电液伺服万能试验机上进行(如图2所示)，并配合DCS-200加载控制系统，在200 N/S的加载速率下，观察并记录试件加载过程中裂尖宏观裂纹扩展路径及试件最终的宏观破坏模式。实验过程中，在试件上下受压端与机头钢块之间布置预先涂抹黄油的橡皮垫，以减弱端部效应的影响；在试件中部放置千分表，表征试件受压过程中的横向变形特性。

图2　电液伺服控制试验机Fig.2 Electro-hydraulic servo control test machine

在此实验条件下，对完整类岩试件进行测试，得到材料的物理力学参数见表1。

表 1完整类岩试件单轴压缩下物理力学参数

Table 1 Mechanical parameters of intact rock-like materials under uniaxial compression

密度/(g·cm-3)弹性模量/MPa单轴抗压强度/MPa单轴抗拉强度/MPa泊松比2.0192.272e323.132.750.2251

本次实验发现，当岩桥倾角为90°时，预制裂隙试件破坏形式主要为预制裂隙尖端微裂纹搭接、贯通、破坏而导致的张拉破坏。这是因为当岩桥走向近似平行于最大主应力方向时，岩桥恰好处于由裂隙面的相对滑动而引发的拉伸应力路径上。以拉应力为驱动力发育的微裂纹主要有翼形裂纹和纵向拉伸裂纹。当裂隙倾角为90°时，以翼裂纹为主。

典型裂尖微裂纹的搭接、贯通、破坏模式示意图如图3所示。试件加载过程中，裂纹首先在预制裂隙相邻两裂尖处开始扩展，使得应力强度因子上升，相邻两裂尖发育出两条翼形裂纹并很快搭接，导致裂纹贯通，岩桥断裂，但此过程比较缓慢；在预制裂隙相邻两裂尖处微裂纹搭接贯通之后，在后续施加荷载的持续作用下，裂纹继续从远离岩桥的两端裂尖处开始扩展，该处两端裂尖应力强度因子迅速增大，发育出翼裂纹，并最终导致试件整体破坏。

图3　预制裂隙尖端微裂纹的搭接、贯通裂纹破坏示意图Fig.3 Microcracks’ overlapping, coalescence and failure pattern at the tip of prefabricated cracks

2无网格法基本原理

2.1紧支试函数加权残量法

(1)

(2)

(3)

(4)

由此得到的系数矩阵为稀疏阵，提高了加权残量法的计算效率。选用不同的试函数，就得到不同的无网格法。

2.2加权移动最小二乘近似

(5)

(6)

将求解域Ω用N个节点离散，在每个节点xI(I=1,2,…,N)处定义一个权函数ωI(x)=ω(x-xI)，该权函数只在节点xI附近的一个有限区域ΩI中大于0，在该邻域外为0，即该函数是紧支的。区域ΩI称为ωI(x)的影响域，如图4所示。

图4　节点的影响域Fig.4 Domain of influence of x.

最小二乘法要求近似函数在各已知节点处误差的加权二阶范数最小，对局部近似函数构造加权二阶能量模，即

(7)

取最小值。由式(7)对ai(x)求导可得:

A(x)a(x)=B(x)u

(8)

(9)

(10)

由式(8)可得待定系数a(x)：

(11)

将式(11)代入式(6)，可得：

(12)

(13)

图5　局部近似和全局近似Fig.5 Local approximation and global approximation

3裂尖处裂纹扩展数值模拟

3.1裂纹扩展原则

本文选用最大周向拉应力准则作为裂纹的断裂扩展判据，即当最大周向拉应力强度因子达到临界值时，裂纹开始扩展。断裂准则如下：

(14)

式中，Kθc为周向拉应力强度因子临界值；KI与KII分别为I型与II型裂纹的应力强度因子；θ为裂纹扩展角度。

裂纹沿最大周向拉应力强度因子的方向开裂，开裂角度θ可根据下式确定：

KIsinθ+KII(3cosθ-1)=0

(15)

利用裂尖应力场可求出KI与KII，将计算出的KI与KII代入式(15)，可得开裂角度θ为：

(16)

3.2裂纹扩展过程数值计算

对试件进行准静态分析，判断裂尖应力是否达到最大周向拉应力测定值。图6～8分别给出了宏观裂纹产生之后裂纹尖端应力场的变化规律。从下图可以看出，裂纹尖端的应力已经远远大于材料本身的最大拉应力，说明材料此时已经破坏，宏观裂纹已经产生并继续扩展，τxy和σy发生振荡。

从裂尖位移变化关系图9～10可以看出，在裂纹扩展期间，裂纹扩展方向主要以y方向为主，这与实验所得裂纹路径的方向基本符合。

图6　裂纹尖端的σxFig.6 σxat the tip of the crack

图7　裂纹尖端的σyFig.7 σyat the tip of the crack

图8　裂纹尖端的τxyFig.8 τxyat the tip of the crack

图9　X方向的位移随裂尖的变化Fig.9 X direction’s displacement at the tip

图10　y方向的位移随裂尖的变化Fig.10 Y direction’s displacement at the tip

在利用无网格法进行裂纹扩展过程数值模拟时候，采用了裂尖能量释放法和小步长扩展技术。计算流程如图11所示。

具体实施过程中，在每次加载步末，判断裂纹是否扩展。若裂纹开始扩展，则保持荷载水平，进行能量释放。在能量释放过程中，若裂纹系统能回到一个新的平衡状态，则说明裂纹在当前荷载水平下是稳定的，可以继续施加荷载进行下一轮新的迭代过程；若系统无法回到新的平衡状态，则说明裂纹发生了失稳扩展，系统发生破坏，计算结束。

为了探讨裂纹的扩展路径，假定裂纹每次扩展长度为Δa=2.15 mm，整个路径共模拟15次，计算裂尖的坐标，给出裂纹扩展路径。由于裂隙是近似对称扩展存在的，表2给出了裂纹扩展过程中上部裂纹尖端位置坐标，图12给出了裂纹扩展路径与实验结果的对比。

图11　裂纹扩展过程计算流程图Fig.11 Flow chart of numerical calculation of crack growth

Table 2 Coordinates of the tip of the upper crack during the crack propagation process

步骤X/mmY/mm1002-2.322.153-4.414.304-6.136.455-7.458.606-8.4710.757-9.0212.908-9.6615.059-10.5017.2010-10.9119.3511-10.9521.5012-11.0623.6513-11.3025.8014-11.3827.9515-11.6430.10

(a)数值计算所得裂隙扩展路径;(b)裂纹扩展路径对比图12　实验与数值模拟所得裂纹扩展路径对比Fig.12 Comparison chart between the experiment and simulation results of the crack growth path

从图12可以看出，计算所得裂纹扩展路径与实验所得裂纹扩展路径基本吻合，说明利用无网格方法模拟裂纹扩展问题的正确性及解决不连续扩展问题的有效性。

4结论

1) 多裂隙岩体裂纹扩展过程受裂隙倾角和岩桥倾角的影响，当裂隙倾角为45°，岩桥倾角为90°时，即岩桥接近平行于最大主应力方向时，预制裂隙试件破坏形式以张拉破坏为主。

2) 预制裂纹首先在相邻两裂尖处开始扩展，裂纹贯通岩桥破坏之后，远离岩桥的两端裂尖处裂纹开始扩展并发育出翼裂纹。

3)无网格法基于加权移动最小二乘法，能同时处理连续与不连续问题，在提高数值解精度的同时也方便了试函数的建立，且前后处理过程相对于有限元更简单，适用于裂纹扩展的数值模拟；而裂纹扩展路径的模拟结果与实验结果基本吻合，说明无网格法能正确有效的解决不连续问题。

4) 通过数值模拟所得到的裂尖应力场和位移场的变化关系，可为进一步研究裂隙岩体的破坏机理提供依据。

参考文献：

[1] 张平, 李宁, 贺若兰. 含裂隙类岩石材料的局部化渐进破损模型研究[J]. 岩石力学与工程学报, 2006, 25(10): 2043-2050.

ZHANG Ping, LI Ning, HE Ruolan. Research on localized progressive damage model for fractured rocklike materials[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(10): 2043-2050.

[2] 关宝树, 熊火耀，翁汉民．裂隙岩体强度的试验研究[J]．西南交通大学学报, 1982(1)：12-23.

GUAN Baoshu ，XIONG Huoyao，WENG Hanmin．The experimental study of the strength of jointed rock mass[J]．Journal of Southwest Jiaotong University，1982 (1)：12-23．

[3] Wong R H C, Chau K T. Crack coalescence in a rock-like material containing two cracks[J]. International Journal of Rock Mechanics and Mining Sciences, 1998, 35(2):147-164．

[4] Bobet A, Einsmln H H. Fracture coalescence in rock-type materials under uniaxial and biaxial compression[J]. International Journal of Rock Mechanics and Mining Sciences, 1998, 35(7):863-888．

[5] 张平, 李宁, 贺若兰,等. 动载下两条断续预制裂隙贯通机制研究[J]. 岩石力学与工程学报，2006, 25(6):1210-1217.

ZHANG Ping, LI Ning, HE Ruolan, et al. Mechnism of fracture coalescence between two pre-existing flaws under dynamic loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(6):1210-1217.

[6] 车法星, 黎立云, 刘大安．类岩材料多裂纹体断裂破坏试验及有限元分析[J]．岩石力学与工程学报，2000，19(3):295-298．

CHE Faxing，LI Liyun，LIU Daan．Fracture experiments and finite element method analysis for multi-cracks body of rock-like material[J]. Chinese Journal of Rock Mechanics and Engineering, 2000，19(3):295-298．

[7] 王元汉, 徐钺, 谭国焕，等．岩体断裂的破坏机理与计算模拟[J]. 岩石力学与工程学报，2000, 19(4)：449-452．

WANG Yuanhan, XU Yue, TAN Guohuan, et al．Fracture mechanism and calculation of rock fracture[J]．Chinese Journal of Rock Mechanics and Engineering，2000, 19(4):449-452．

[8] 黄明利，唐春安，朱万成. 岩石破裂过程的数值模拟研究[J]．岩石力学与工程学报，2000, 19(4):468-471．

HUANG Mingli, TANG Chunan, ZHU Wancheng. Numerical simulation on failure process of rock[J]．Chinese Journal of Rock Mechanics and Engineering, 2000, 19(4):468-471．

[9] Belytschko T, Krongauz Y, Organ D, et al. Meshless methods: An overview and recent developments. comput[J]. Methods Appl Mech Engrg, 1996(139):3-47.

[10] Li S, Liu W K. Meshfree and particle methods and their applications[J]. Appl Mech Rev, 2002, 55(1):1-34.

[11] 张雄, 宋康祖, 陆明万. 无网格法研究进展及其应用[J]. 计算力学学报, 2003, 20(6):57-63.

ZHANG Xiong, SONGR Kangzu, LU Mingwan. Research progress and application of meshless method[J]. Chinese Journal of Computational Mechanics, 2003, 20(6):57-63.

[12] 庞作会,葛修润,王水林,等.无网格伽辽金法(EFGM)在边坡开挖问题中的应用[J].岩土力学,1999, 20(1):61-64.

PANG Zuohui, GE Xiurun, WANG Shuilin, et al. Application of element free galerkin method(EFGM) to the slope excavation[J]. Rock and Soil Mechanics, 1999, 20(1):61-64.

[13] ZHANG Xiong, LIU Mingwan. A 2-D meshless model for jointed rock structures[J]. Int J Numer Methods Engrg, 2000, 47(10):1649-1661.

[14] 寇晓东, 周维垣. 应用无单元法近似计算拱坝开裂[J]. 水利学报, 2000, 10:28-35.

KOU Xiaodong, ZHOU Weiyuan. The application of element-free menthod to approximate calculation of arch dam crack propagation[J]. Journal of Hydraulic Engineering, 2000, 10:28-35.

[15] 庞作会, 葛修润. 无网格伽辽金法(EFGM)模拟不连续面[J]. 工程地质学报， 2000, 8(3):364-368.

PANG Zuohui, GE Xiurun. Simulation discontinuty with element free galerkin method(EFGM)[J]. Journal of Engineering Geology, 2000, 8(3):364-368.

[16] 徐次达. 固体力学加权残量法[M]. 上海: 同济大学出版社, 1987.

XU Cida. Weighted residual mehtod in solid mechanics[M]. Shanghai: Tongji University Press, 1987.

[17] 张雄, 宋康祖, 陆明万. 紧支试函数加权残量法[J]. 力学学报, 2003, 35(1):43-49.

ZHANG Xiong, SONG Kangzu, LU Mingwan. Weighted residual method based on compactly supported trial functions[J]. Chinese Journal of Theoreyical and Applied Mechanics, 2003, 35(1):43-49.

(编辑蒋学东)

Application of element-free method to fracture analysis of rock-like materials with double prefabricated cracks

LIU Jingshuo1,2，CAO Ping1, FAN Jinxing2， PU Chengzhi1， HUANG Xuejiao1

(1.School of Resource and Safety Engineering, Central South University, Changsha 410083, China；2. Hunan Technical College of Water Resources and Hydropower, Changsha 410131, China)

Abstract：In order to investigate the failure mode and the crack expansion route, the compression test on rock-like specimens with two prefabricated cracks was carried out by utilizing the servo control uniaxial loading test machine. Afterward, the fundamental theory of the element-free method is briefly illustrated. Then, the proposed method, which is based on the weighted moving least square approximation, is extended to numerical simulation of fracture behavior of rock-like specimens with multiple cracks. Specifically, the fracture behavior and the crack propagation path of rock-like specimens are numerically calculated and analyzed. At the same time, the curves of the stress and displacement fields at the tip of the fissures are obtained. Finally, the numerical solutions are compared with the available experimental results to verify the advantages, accuracy and availability of element-free method to fracture and discontinuous deformation problems. It is demonstrate that this method is superior to other methods in terms of fracture analysis.

Key words：prefabricated cracks；rock-like material；uniaxial compression；element-free method；numerical analysis；crack propagation

中图分类号：TU45

文献标志码：A

文章编号：1672-7029(2016)03-0447-07

通讯作者：刘京铄(1979-)，男，湖南邵东人，讲师，博士研究生，从事岩土工程方面的研究；E-mail：ljs7327@163.com

基金项目：国家自然科学基金重点资助项目(51174228，51274249)；湖南省教育厅科技项目(14C0746)

收稿日期：2015-07-22