叶少锋，周 伟，常晓林，马 刚

(1. 武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2. 武汉大学 水工岩石力学教育部重点实验室，武汉 430072)

0 引 言

在众多坝型，土石坝因为其具有充分利用当地材料、对坝址地形地质条件要求较低、抗震性好等优势而在坝工设计中得到十分广泛的运用。其中心墙坝是土石坝常用的坝型之一。众所周知，心墙水力劈裂是大坝安全的一个重要隐患。在心墙坝的设计过程中，水力劈裂是一个必须重点研究和防范的关键问题。

许多学者对心墙水力劈裂问题进行了许多系统深入并且卓有成效的研究。黄文熙[1]在1982年指出水力劈裂是土体中因水压力抬高导致裂缝产生与扩展的一种物理现象。这一定义很好地反映了水力劈裂现象的本质。张丙印等[2]提出一种水压楔劈模型来解释水力劈裂裂缝的形成机理。朱俊高等[3]认为高水力梯度是形成水力劈裂的最根本因素。以上三种对水力劈裂形成机理的观点具有一定的相似之处。在心墙水力劈裂判定方面，殷宗泽等[4]提出将心墙内的总应力与心墙表面受到的水压力进行对比，以此分析是否会发生水力劈裂，其中总应力是由土体中有效应力与孔隙水压力叠加得到的。在总结并比较总应力法、综合法和有效应力法三种水力劈裂判别方法的基础上，陈五一等[5]提出一种新的安全系数来判定水力劈裂。在数值模拟方面，李全明等[6]通过综合比奥固结理论与弥散裂缝理论，提出了一种有限元计算方法来数值模拟水力劈裂现象。不过该方法需要预设裂缝位置，而且无法模拟裂缝贯穿整个心墙的过程。杨艳等[7]从细观角度出发，应用离散元方法研究心墙水力劈裂的过程，该研究表明裂缝形成的原因是土体受到的张拉应力大于抗拉强度。在水力劈裂室内物理试验方面，袁俊平等[8]自制了一套水力劈裂试验设备，在心墙试样中设置不同的初始裂缝长度，分析试验结果得出：试样的初始裂缝长度越大，则形成水力劈裂所需的水压力越小。

综合以上学者的研究可以看出，心墙水力劈裂是高水力梯度作用下导致裂缝发生与扩展的一种流固耦合的物理现象。现有的研究方法或理论主要关注于土料应力、孔隙水压力和裂缝等方面的其中一两点，不够全面。作为一种能够有效模拟裂缝萌生与扩展的数值模型，内聚力模型已经广泛运用于岩石等多种材料的开裂行为研究[9-11]，也为分析水力劈裂过程中裂缝的扩展行为提供了新手段。本文采用考虑孔隙水压力的内聚力模型对心墙水力劈裂现象进行了数值模拟，并且分析了心墙水力劈裂的发生机理。该方法不仅模拟了流固耦合作用，而且能够直观地跟踪裂缝萌生、扩展直至最终形成贯穿性裂缝的整个过程中孔隙水压力、应力变形等各个方面的状态。通过具体分析，该数值模拟所获得的成果符合现有的室内物理试验结果与理论，为今后心墙坝设计人员分析和解决水力劈裂问题提供了新手段。

1 考虑孔隙水压力的内聚力模型的基本原理

1.1 基本思路

本文数值模拟是在连续—离散耦合分析方法(FDEM)的框架下实现的。FDEM方法最早是由Munjiza等提出，其基本思路是在心墙结构内划分有限元网格(本文采用三角形单元)，然后在所有相邻单元的公共边上插入无厚度的界面单元，通过界面单元的断裂来模拟裂缝的萌生与扩展，具体内容可以参考文献[12-17]。其中界面单元是该方法的关键所在，本文采用内聚力模型来定义界面单元的力学行为。由于内聚力模型只能模拟裂缝的形成，为了实现对水力劈裂流固耦合作用的模拟，必须在内聚力模型中考虑渗流作用，建立考虑孔隙水压力的内聚力模型，简称孔压内聚力模型。其模拟渗流的思路是，将流体看作是不可压缩的，并且流体前缘与裂缝尖端是重合的，计算时不考虑流体的滞后效应；有限元网格内部孔隙流体流动满足Darcy定律，界面单元裂缝内流体流动包括沿裂缝面的切向流动和垂直裂缝面的法向流动。在水压力的驱动下，流体的切向流动驱使裂缝向前扩展，法向流动使裂缝内流体渗透到周边介质中。在计算过程中，水压力随着渗流作用进行传递，改变模型应力状态，影响界面单元地张开与闭合状态；反之，模型的应力状态和界面单元地张开与闭合状态的变化必然引起孔隙与裂缝内流体渗流的变化。两者形成相互制约和相互作用的关系。这种关系就是流固耦合。

在FDEM框架下运用孔隙水压力内聚力模型模拟水力劈裂有两个优势。一是能够克服连续介质力学模型模拟裂缝的困难，有效模拟裂缝的萌生与扩展；二是既能考虑水力劈裂过程中孔隙流体和裂缝流体渗流作用的差异，又能反映裂缝流体与孔隙流体的流量交换。

1.2 孔隙流体渗流方程

假定心墙颗粒骨架和流体均具有不可压缩性，流体的连续性方程为：

(1)

式中：εv为体积应变；q为渗透流速矢量。

流体渗流的 Darcy 定律为：

q=k▽h

(2)

式中：k为渗透系数矢量；h为水头。

1.3 单元应力平衡方程

单元的应力平衡方程为：

σij,j+fi=0

(3)

式中：σij表示总应力；fi表示体力。

本文水力劈裂分析时均采用有效应力进行计算，其中应力状态以拉应力为正，孔隙水压力pw以压应力为正，有效应力σ′ij的表达式如下：

σ′ij=σij+δijpw

(4)

由式(3)和式(4)可得出用有效应力表示的平衡方程：

σ′ij-pw,j+fi=0

(5)

1.4 内聚力模型

内聚力模型理论起源于断裂力学。在20世纪60年代，Dugdale[18]和Barrenblatt[19]为了解决线弹性断裂力学中裂缝尖端处应力奇异性问题先后提出了内聚力模型(CZM)理论。该理论认为裂缝尖端处应力为裂缝开度的函数。本文通过采用内聚力模型定义界面单元受到荷载作用时应力与位移的响应关系来描述界面单元的力学行为。界面单元变形分为弹性变形、起裂、损伤演化和完全失效4个阶段；单元起裂之前，发生弹性变形；起裂后，进入损伤阶段，刚度逐渐降低；当刚度降低至0时，单元完全失去承载能力(如图1所示)。

图1 内聚力模型与界面单元Fig.1 Schematic representation of the cohesive zone model and cohesive interface element

在界面单元出现损伤之前，其线弹性的本构关系如下：

(6)

式中：tn、ts分别表示法向应力与切向应力；Kij为刚度矩阵；εn、εs分别表示法向应变与切向应变，本文计算时未考虑应力、应变的法向分量与切向分量之间的耦合关系。

法向应变、切向应变定义为：

(7)

式中：δn、δt分别表示法向位移和切向位移；T0为界面单元的本构厚度。因为本文使用无厚度的界面单元，所以取T0=1，界面单元的应变与相应方向上的位移大小相等。界面单元的线弹性本构方程可写成：

(8)

式中：kn、ks分别表示法向刚度和切向刚度。

为了模拟裂缝萌生过程，需要定义有效的起裂准则。本文假定材料破坏是由拉应力和剪应力共同作用下导致的，所以模拟时采用二次应力准则判断界面单元是否开裂。当各个方向的应力与其相应的临界应力的比值的平方和达到1时，开始形成裂缝，公式如下：

(9)

式中：t0n为界面单元法向应力的临界值，即为心墙土料的抗拉强度；t0s为界面单元切向应力的临界值，即为心墙土料的抗剪强度；〈〉为Macaulay括号：

(10)

界面单元起裂后，即进入裂缝扩展阶段，此时采用标量损伤因子来表示界面单元刚度的退化的程度。单元应力与位移的表达式为：

(11)

ts=(1-D)ksδs

式中：D表示标量损伤因子，当D=0时，表明界面单元未出现损伤；当D=1时，表明单元完全裂开，失去承载能力。

定义界面单元的等效应力和等效位移如下：

(12)

式中：teff为等效应力；δm为等效位移。

通过定义界面单元在损伤过程中耗散的能量，即可确定单元的损伤演化过程。耗散能量是指通常所说的断裂能，它等于应力-分离曲线下面所包含面积的大小。本文损伤演化准则采用基于断裂能的Benzeggagh-Kenane 准则[20]，其表达式如下：

GS=Gs,GT=GS+Gn

(13)

式中：GC为复合断裂能；GCn为Ⅰ型断裂能；GCs为Ⅱ型断裂能；本文中，η取2。

基于断裂能的损伤演化准则中界面单元的软化准则分为线性软化和指数软化。由于指数软化曲线更接近心墙料拉伸断裂的应力-应变曲线，本文采用指数软化准则。指数软化准则中，标量损伤因子的表达式如下：

(14)

式中：G0为界面单元开始出现损伤时对应的弹性能；δ0m、δfm分别指界面单元开始出现损伤时与完全失效时对应的等效位移。

1.5 水力裂缝模型

在考虑孔隙水压力的内聚力模型中，裂缝中流体的流动现象如图2所示。流体在裂缝中的流动可分为法向流动和切向流动，切向流动驱使裂缝向前扩展，法向流动反映了流体的滤失效应，即部分流体渗透到周边介质中。

图2 界面单元流体流动Fig.2 Flow within cohesive interface element

流体的切向流动遵守润滑方程：

qd=-kt▽p

(15)

式中：q为界面单元的体积流率密度矢量；d为界面单元地张开宽度；▽p为沿着界面单元的水压力梯度；kt为切向渗透系数，由雷诺方程求得：

(16)

式中：μ为流体的黏滞系数。

假设裂缝的上、下表面为可透水层，其渗透系数即为界面单元的滤失系数，在水压力的驱动下流体产生法向流动，法向流速的计算公式如下：

qt=ct(pi-pt)

qb=cb(pi-pb)

(17)

式中：qt和qb分别为界面单元上、下表面的流速；ct和cb分别为界面单元上、下表面的滤失系数；pt和pb分别为与界面单元上、下表面相邻的多孔介质的孔隙水压力；pi为界面单元内部的孔隙水压力。

2 基于内聚力模型的心墙水力劈裂数值模拟

2.1 内聚力模型力学参数选取

模型参数的选取是数值模拟的一个关键因素。笔者查阅了大量文献，尚未找到将内聚力模型运用于土石坝心墙材料研究的文献资料，因而缺乏心墙材料的内聚力模型力学参数的参考数据。为此，笔者在进行心墙水力劈裂的数值模拟前，先采用内聚力模型对文献[21]中的心墙土料试样的巴西圆盘劈裂试验进行数值模拟，初步确定心墙材料内聚力模型参数。由于堆石坝横河向的宽度一般远大于顺河向的宽度，结构变形特点更接近于平面应变，所以本文以平面应变的角度来数值模拟巴西圆盘劈裂和水力劈裂试验。

巴西圆盘劈裂物理试验采用直径为39.1 mm，高度为80 mm的圆盘试样。试验分别选取5组不同含水率的试样和5组不同密度的试样。本文对其中密度最大且比较接近工程实际情况的一组试样进行模拟，其密度为1.65 g/cm3，含水率为16%。数值模型为直径39.1 mm的圆形平面模型，划分为4 464 个界面单元和3 024 个三角形实体单元，具体模型如图3所示。数值模拟时通过上下两个刚性板以一定的速率向圆盘中心移动来进行加载，本文刚性板的加载速率为9×10-3m/s。随着刚性板的移动，圆盘逐渐形成上下贯通的裂缝，如图4所示。计算结束后，通过提取刚性板受到的反力与圆盘裂缝方向的径向位移，绘制出拉应力与径向位移的关系曲线并与文献[21]的物理试验结果进行比较。其中由反力计算拉应力的公式如下：

(18)

式中：σt为拉应力；F为刚性板受到的反力；w和l分别表示试样直径和厚度。其中，对于二维数值模拟，试样的厚度取为1。

图3 巴西圆盘劈裂试样网格划分Fig.3 Mesh configuration of Beazilian samples

图4 巴西圆盘劈裂数值试样裂缝图Fig.4 Cracks of numerical Beazilian samples

运用内聚力模型所需的力学参数如表1所示，包括密度ρ、弹性模量E、泊松比υ、抗拉强度tn、抗剪强度ts、法向刚度kn、切向刚度ks、Ⅰ型断裂能GCn、Ⅱ型断裂能GCs等。本文数值模型中，没有考虑到土料的含水量这一因素，但是含水量不同，试样的力学性能也会不同，含水量对试样的影响会通过力学参数体现出来。数值模拟时，不断调整各个参数进行试算，使试样最大拉应力与其对应的径向位移与物理试样结果相吻合。最终各个参数的取值如表1所示，数值试验与物理试验的拉应力和径向位移曲线如图5所示。从图5可以看出，数值试验曲线在拉应力达到最大值后，急剧下降至0，而物理试验曲线则没有这一过程。这主要是由于物理试验最后压裂试样的过程很短暂，仪器来不及记录数据，试样就已经完全失去承载能力，而数值试验则没有这个问题，因而数值试验记录了最后压裂的过程曲线。除了上面这一差异，该数值试验的应力变形曲线基本与物理实验结果相符合，因此初步确定的内聚力模型参数能够模拟心墙试样的劈裂力学行为。

表1 内聚力模型力学参数Tab.1 Mechanical parameters of cohesive zone model

图5 物理试验与数值试验的拉伸应力-位移曲线对比[21]Fig.5 Experimental tensile stress-displacement response comparised to numerical predictions[21]

2.2 心墙水力劈裂数值模拟

2.2.1 建立计算模型

本文选取文献[2]中的室内水力劈裂试验数据资料进行水力劈裂二维数值模拟。该实验物理模型如图6所示，心墙试样为边长20 cm的立方体，试样中部有一层长10 cm、厚0.5 cm的砂层，表示心墙的初始裂缝或者缺陷。进行试验时，先在模型的左侧和上侧通过气囊施加初始应力σ3和σ1使试样固结，然后在左侧施加水压力，通过观察窗口和测量渗流量的变化判断是否形成水力劈裂。

图6 物理模型示意图[2]Fig.6 Sketch of physical model[2]

数值模型取20 cm×20 cm的正方形，模型中部有一道长为10 cm的初始裂缝或缺陷。模型由3 628 个实体单元和5 362 个界面单元组成，如图7所示。由于水力劈裂数值试验是采用考虑孔隙水压力作用的内聚力模型，所以计算时所用参数除了表1中力学参数外，还有流体计算参数，包括渗透系数k、黏滞系数μ、重度γw、孔隙比e和滤失系数ct与cb等，参数的选取参照文献[22]，具体数值见表2。

表2 流体计算参数Tab.2 Computational parameters for fluid

图7 数值模型示意图Fig.7 Sketch of numerical model

2.2.2 数值模拟结果分析

本文模拟心墙水力劈裂的过程如下：首先对模型上表面和上游表面分别施加σ1和σ3的面荷载进行固结；然后在模型上游施加水压力；一开始取一较小的水压力值，如果计算结果没有形成裂缝，则设置比较大的水压力值重新进行计算；如此重复上述步骤，直到计算结果出现贯穿上下游的裂缝为止(如图8所示)，此时施加的水压力就是劈裂水压力值。

图8 水力劈裂后的贯穿性裂缝Fig.8 Penetrating cracks after hydraulic fracturing

本文根据室内物理试验σ1、σ3的3组初始应力组合方案(见表3)分别进行了水力劈裂数值模拟，数值模拟与物理试验的劈裂水压力值对比如表3与图9所示。根据图9所呈现的规律可以得出，随着初始竖向应力的增大，数值模拟与物理试验的劈裂水压力值均逐渐增大，劈裂水压力值与竖向应力大致呈线性相关。另外，数值模拟的劈裂水压力值要大于室内试验的劈裂水压力。笔者认为这可能源于3个原因。第一，室内心墙材料的物理特性受含水率、密度等因素的影响很大，两个试验并非采用同一心墙土料进行实验，水力劈裂试样的力学参数与巴西圆盘劈裂试样可能有所差异，本文直接采用巴西圆盘劈裂试样的参数来模拟水力劈裂，可能导致结果有一定误差。第二，室内物理试验在达到水力劈裂破坏前，经过了多次水的渗透作用，心墙土在水的反复作用下发生软化作用，强度有所下降，而数值模拟没能考虑到这一因素的影响。第三，数值模拟在考虑实际情况的基础上必须做一定合理假定和简化，这也会使数值试验与物理试验结果之间产生一些无法避免的误差。

表3 数值模拟结果[2]Tab.3 Numerical results[2]

图9 数值模拟结果与室内试验结果对比[2]Fig.9 Numerical and laboratory results[2]

为了从水压力角度探究心墙水力劈裂的发生机理，从第二组应力组合的数值试验分别提取水平裂缝内和垂直裂缝路径上的水压力值。图10为裂缝扩展至距上游面0.15 m处时刻的水平裂缝内水压力分布图。从图10 可以看出，水压力分布曲线可分为4段。第1段曲线横坐标区间为距离上游面0～0.10 m，此处为模型的初始裂缝或者缺陷，水压力等于上游面的水压力值0.33 MPa。第2段曲线横坐标区间为距离上游面0.10～0.14 m，此处界面单元损伤因子已经到达1，界面单元完全裂开形成张开裂缝，随着水流进入裂缝，水压力由初始裂缝向新生裂缝处传播，由于水的黏滞作用，水压力在向前传播的过程中逐渐减弱。第3段曲线横坐标区间为0.14～0.15 m，水压力急剧下降，呈现很大的水力梯度。此处为裂缝的最前端，在高水力梯度的作用下，界面单元出现损伤，损伤因子介于0～1，裂缝由闭合状态逐渐张开。第4段横坐标区间为距离上游面0.15～0.20 m，水压力值接近于0。此处界面单元完好未出现损伤，没有裂缝出现，水压力以Darcy定律在介质中传播，水压力值接近于下游面的初始水压力值。图11为垂直穿过距离上游面0.12 m处水平裂缝的路径上的水压力值分布图，提取水压力时刻与图10一样，其中横坐标等于0表示裂缝的位置，横坐标大于0表示位于裂缝的下方，横坐标小于0表示位于裂缝的上方。从图11中可以看出，裂缝中的水压力要远大于裂缝上下方介质中的孔隙水压力值。这是由于裂缝形成后，初始裂缝内的水压力迅速向裂缝传播，而裂缝上下方介质的孔隙水压力则是上游面初始水压力满足Darcy定律向下游传播的结果；在心墙材料较低渗透性与裂缝内水压力未能及时向周边介质消散的条件下，裂缝内水压力扩散速率要大于孔隙水压传播速度，因而裂缝内水压力大于上下方介质中的孔隙水压力。在裂缝内与上下方介质中的水压力差的作用下，上下方的介质分别向上与向下移动，裂缝逐渐张开，直到裂缝内水压力向周边介质消散达到水压平衡为止。综上所述，水力劈裂过程中，高水力梯度和裂缝与周边介质的水压差是劈裂裂缝前端的两个重要特征。结合文献[2]描述的“水压楔劈效应”可以看出，高水力梯度和裂缝与周边介质的水压差是裂缝萌生与扩展的两个重要原因。

图10 裂缝内的水压力Fig.10 Pressure in the fractures

图11 垂直水平裂缝路径上的水压力Fig.11 Pore pressure in a vertical path

为了初步探讨水力劈裂过程中界面单元的应力状态，从第二组数值试验选取开裂路径上任意一个界面单元，提取其切向应力、法向应力以及损伤因子的历时曲线，如图12所示。图中，法向应力正值表示拉伸应力，切向应力正值表示逆时针方向，0～1 s的时间段为施加初始应力的阶段，1～3 s为施加水压力阶段。从图中可以看出，法向应力在施加初始应力阶段由0逐渐增加到0.3 MPa的压应力；在施加水压力阶段的前期，还未产生裂缝，孔隙水压力抵消了部分初始应力，压应力逐渐变小到0.2 MPa左右；在2.2 s时刻界面单元开始出现损伤，上游水压力传递到裂缝中，法向应力由压应力减小至0并转变成拉应力；当损伤因子达到1时，界面单元完成裂开失去承载能力，拉应力全部释放，法向应力降低至0。界面单元在施加水压力的前期受渗流作用存在部分切向应力，在界面单元裂开过程中，切向应力逐渐降低至0。从图中的曲线可以发现，界面单元在开裂过程中，法向应力的变化幅度以及起裂时的应力值要远高于切向应力。从以上的分析中可以看出，在二维水力劈裂数值模拟中，水力劈裂的破坏机理是以拉伸破坏为主。值得注意的是，在三维情况下，心墙结构的应力状态要复杂得多，其水力劈裂的破坏机理还有待进一步研究。

图12 界面单元切向应力、法向应力和损伤因子历时曲线Fig.12 Time-history curves of normal stress, shear stress and the damage variable of cohesive elements

为了进一步验证数值模拟与物理试验具有相同的一般性规律，本文针对第二组初始应力组合分别建立4组具有不同初始裂缝长度的模型进行水力劈裂数值模拟，其中将初始裂缝长度相对模型边长的百分比定义为裂缝相对长度。除了裂缝相对长度不同外，4组数值试验所施加的初始应力和约束条件等均相同。最终计算结果如图13所示，从图13中可以看出，劈裂水压力值随着初始裂缝相对长度的增加逐渐降低，此规律与文献[8]的物理试验结果相同。该规律表明在受到相同外界荷载和边界条件下，如果有较长的初始裂缝或者更大的缺陷，那么心墙更容易形成水力劈裂。因为初始裂缝越长，裂缝尖端处孔隙水压力集中程度越剧烈，导致形成更大的水力梯度，促进了裂缝的萌生与扩展。在心墙坝建设过程中，控制心墙上游面附近初始裂缝或缺陷的大小，将有效降低水力劈裂发生的几率。

图13 初始裂缝相对长度与劈裂水压力的关系Fig.13 Relation between initial crack length and fracturing pressure

综上所述，基于考虑孔隙水压力内聚力模型的心墙水力劈裂数值模拟，虽然计算出的劈裂水压力值偏大，但是，能够反映心墙水力劈裂的一般规律，并且合理地揭示了水力劈裂的发生机理。

3 结 语

(1)本文运用内聚力模型模拟了心墙试样巴西圆盘劈裂试验，计算得出的拉应力—径向位移曲线与物理实验基本相符，初步确定了心墙材料的内聚力模型力学参数。

(2)将考虑孔隙水压力的内聚力模型运用于心墙水力劈裂数值模拟，虽然由于一些不可避免的原因，与室内物理试验相比，结果中的劈裂水压力值大于物理试验；但是，劈裂水压力值与竖向初始应力的关系与物理试验结果相同，即随着竖向初始应力的增大，劈裂水压力值线性增大。

(3)通过探究数值试验中水压力在裂缝及其周边介质的分布，可以得出：高水力梯度和裂缝与周边介质的水压差是水力劈裂裂缝萌生与扩展的两个重要原因。

(4)提取开裂路径上一个界面单元的法向应力、切向应力以及损伤因子，并绘制出历时曲线。分析曲线的规律可以发现，在二维水力劈裂数值模拟中，水力劈裂的破坏机理是以拉伸破坏为主。

(5)进一步模拟设置不同初始裂缝长度的水力劈裂试验，得到如下规律：初始裂缝长度越大，心墙越容易形成水力劈裂，控制心墙上游面附近初始裂缝或缺陷的大小，将有效降低水力劈裂发生的几率。

□

[1] 黄文熙. 对土石坝科研工作的几点看法[J]. 水利水电技术, 1982,(4):23-28.

[2] 张丙印, 李 娜, 李全明. 土石坝水力劈裂发生机理及模型试验研究[J]. 岩土工程学报, 2005,27(11):1 277-1 281.

[3] 朱俊高, 王俊杰, 张 辉. 土石坝心墙水力劈裂机制研究[J]. 岩土力学, 2007,28(3):487-492.

[4] 殷宗泽，朱俊高，袁俊平，等. 心墙堆石坝的水力劈裂分析[J]. 水利学报，2006，37(11)：1 348-1 353.

[5] 陈五一, 赵颜辉. 土石坝心墙水力劈裂计算方法研究[J]. 岩石力学与工程学报, 2008, 27(7):1 380-1 386.

[6] 李全明，张丙印，于玉贞，等. 土石坝水力劈裂发生过程的有限元数值模拟[J]. 岩土工程学报，2007，29(2)：212-217.

[7] 杨 艳, 周 伟, 常晓林, 等. 高心墙堆石坝心墙水力劈裂的颗粒流模拟[J]. 岩土力学, 2012,33(8):2 513-2 520.

[8] 袁俊平, 王启贵. 水力劈裂楔劈效应试验研究[J]. 岩土力学, 2015,36(12):3 360-3 364.

[9] 王永刚, 胡剑东, 王礼立. 金属材料层裂破坏的内聚力模型[J]. 固体力学学报, 2012,33(5):465-470.

[10] 胡 超, 周 伟, 常晓林, 等. 基于内聚力模型的高心墙堆石坝坝顶裂缝模拟及其成因分析[J]. 中南大学学报:自然科学版, 2014,45(7):2 303-2 310.

[11] 马 刚, 周创兵, 常晓林, 等. 岩石破坏全过程的连续-离散耦合分析方法[J]. 岩石力学与工程学报, 2011,30(12):2 444-2 455.

[12] Munjiza A. The combined finite-discrete element method[M]. Cichester：John Wiley and Son，Inc.，2004：233-243.

[13] Munjiza A, Bangash T, John N W M. The combined finite-discrete element method for structural failure and collapse[J]. Engineering Fracture Mechanics, 2004,71(46)：469-483.

[14] Munjiza A, Owen D R J, Bicanic N. A combined finite-discrete element method in transient dynamics of fracturing solids[J]. Engineering Computations, 1995,12(2)：145-174.

[15] Owen D R J, Feng Y T. The modelling of multi-fracturing solids and particulate media[J]. International Journal for Numerical Methods in Engineering, 2004, 60(1)：317-339.

[16] Zhou Wei, Yuan Wei, Ma Gang, et al. Combined finite-discrete element method modeling of rockslides[J], Engineering Computations, 2016,33(5):1-27.

[17] MA Gang, Zhou Wei, Chang Xiaolin, et al. Combined FEM/DEM modeling of triaxial compression tests for rockfills with polyhedral particles[J]. International Journal of Geomechanics, 2014,14(4):04014014.

[18] Dugdale D S. Yielding of steel sheets containing slits[J]. Journal of the Mechanics and Physics of Solids, 1960,8(2):100-104.

[19] Barenblatt G I. The mathematical theory of equilibrium cracks in brittle fracture[J]. Advances in applied mechanics, 1962,7(1):55-129.

[20] Benzeggagh M L, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[J]. Composites science and technology, 1996,56(4):439-449.

[21] 王惠敏. 黏性土抗拉特性试验研究[D]. 南京: 南京大学, 2011.

[22] 龚迪光, 曲占庆, 李建雄, 等. 基于ABAQUS平台的水力裂缝扩展有限元模拟研究[J]. 岩土力学, 2016,37(5):1 510-1 520.