DOI: 10.11817/j.issn.1672-7207.2015.04.035
基于颗粒离散元法的岩石节理面剪切破坏细观机理
王刚1,2，张学朋1，蒋宇静1,3，李博3，吴学震1，黄娜1
(1. 山东科技大学 矿山灾害预防控制省部共建国家重点实验室培育基地，山东 青岛，266590；
2. 山东科技大学 山东省土木工程防灾减灾重点实验室，山东 青岛，266590；
3. 长崎大学 工学部，日本 长崎，852-8521)
摘要：基于颗粒流理论，以含不同粗糙度系数的岩石节理面为研究对象，利用PFC2D程序在恒定法向荷载(CNL)和恒定法向刚度(CNS)边界条件下进行直接剪切试验模拟，分别从宏观和细观角度探讨节理在不同边界条件下剪切过程中的力学演化规律和破坏机制。模拟值与JRC-JCS模型的预测值具有良好的一致性。研究结果表明：岩石节理面的抗剪强度及其对应的剪切位移以及残余强度随法向应力的增加而增大。CNL边界条件下的抗剪强度和残余强度比CNS边界条件的低，而法向位移比CNS边界条件的高；在整个剪切过程中，粒间接触力的分布和残余阶段最大接触力的变化趋势在2种边界条件下存在不同。节理面凸起位置容易发生接触压力集中现象，细观裂纹也主要聚集在凸起的位置且张拉裂纹占绝对优势。不论在CNL边界条件下还是在CNS边界条件下，试件在剪切过程中主要经历弹性阶段、稳定破裂传播阶段、非稳定破裂传播阶段以及残余阶段的稳定破裂传播阶段。
关键词：岩石节理；粗糙度；直剪试验；颗粒流理论；恒定法向刚度
中图分类号：O319.56 文献标志码：A 文章编号：1672-7207(2015)04-1442-12
Meso-mechanism research on shear failure of rock joint based on particle discrete element method
WANG Gang1, 2, ZHANG Xuepeng1, JIANG Yujing1, 3, LI Bo3, WU Xuezhen1, HUANG Na1
(1. State Key Laboratory of Mining Disaster Prevention and Control Co-founded by Shandong Province and
Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266590, China;
2. Shandong Provincial Key Laboratory of Civil Engineering Disaster Prevention and Mitigation,
Shandong University of Science and Technology, Qingdao 266590, China;
3. Faculty of Engineering, Nagasaki University, Nagasaki 852-8521, Japan)
Abstract: The shear test of rock joint with different joint roughness coefficients was simulated using the particle flow code (PFC), under constant normal load (CNL) and constant normal stiffness (CNS) boundary conditions, in order to assess the macro- and meso-mechanical properties and failure mechanism of rock joints subjected to shear. The simulation method was verified by comparing simulation results with those estimated by Barton’s empirical equation. The results show that the peak and residual shear strength and shear displacement at peak shear stress increase with the increase of normal load. The peak and residual shear strength under CNL boundary condition are less than those under CNS boundary condition, while the normal displacement is greater. Due to the difference in normal stress under CNL and CNS boundary conditions during shear processes, different distributions of contact forces and the maximum contact forces can be observed at residual stage. Concentration of compression force on contacts is prone to occur near the major asperities, thereby micro-cracks are also concentrated near the corresponding asperities, where tensile cracks dominate the micro-cracks. There are four main stages during a shear process under both boundary conditions: elastic stage, stable fracture propagation stage, unstable fracture propagation stage and stable fracture propagation stage at residual stage.
Key words: rock joint; roughness; direct shear test; particle flow theory; constant normal stiffness
裂隙岩体宏观上是由节理或裂隙切割成一块块、相互排列与咬合着的岩块所组成的。岩体中往往具有明显的地质遗迹，如不整合、褶皱、断层、节理、劈理、解理、层理、片理等。由于节理的存在造成岩体介质的不连续性和各向异性，从而显著的改变岩体的各类力学性质(如抗压强度，弹性模量，抗剪强度)。近年来，岩石节理在地下洞室和边坡稳定性分析研究中的支配地位越来越受到重视，国内外已从理论及模型试验等方面对岩石节理面的剪切强度和变形特性进行了大量研究：自Patton[1]提出关于齿形结构面的双线性模型，又有许多抗剪强度准则相继提出[2-3]，最具代表性的是Barton等[3]提出的JRC-JCS模型；随后又在试验和理论模型方面解释节理剪切的峰后效应及凸起退化[4-7]。但是以上研究主要集中在恒定法向荷载的边界条件，此种边界条件只能模拟作用于节理面法向应力(自重)不变的如未支护的边坡等问题(见图1(a)，其中σ为法向应力；τ为抗剪强度)；由于深埋岩石节理变形，会引起围压的变化，同时围压变化会引起节理面的法向应力变化(见图1(b))，因此，恒定法向荷载的边界条件不再适合此类工程问题。从20世纪80年代起，开始对恒定法向刚度的条件进行研究：开始的剪切装置使用弹簧来表现围压影响及围岩的变形刚度[8-9]，此类实验方法不能反映剪切过程中节理面的变形对垂直刚度产生的影响，且不能自由设定垂直刚度系数。对此，Jiang等[10]提出用数控系统代替弹簧模拟节理面围岩的变形刚度特性，克服了上述缺点。随着计算机技术的发展，数值模拟已经与理论分析和模型试验一起成为科学研究的三大支柱。针对岩石节理面力学特性的研究，数值模拟已经成为一种有效的研究手段。目前国内外采用的方法主要是有限元[11-12]、边界元[13]和离散元[14-17]，但主要还是以恒定法向荷载的边界条件为主。针对目前剪切试验主要集中在单一的恒定法向荷载边界条件的研究现状，本文作者基于颗粒流理论，借助PFC2D程序的wall刚性边界和clump刚性边界分别实现岩石节理面直接剪切试验恒定法向荷载(CNL)和恒定法向刚度(CNS) 2种边界条件，并分别从宏观和细观角度探讨节理在不同边界条件下剪切过程中的力学演化规律和破坏机制。
图1 CNL和CNS剪切实验概略图
Fig. 1 Sketch map of direct shear test under CNL and CNS boundary conditions
1 岩石材料的颗粒流程序表达
1.1 颗粒流理论
PFC2D[18]利用刚性圆形颗粒集合体来表征材料，其中颗粒之间独立运动，只在接触节点上相互作用。PFC涉及的计算只需要简单的力-位移定律、牛顿运动定律和几个参数来控制颗粒、接触的相互作用。通过力-位移定律更新接触力，通过牛顿运动定律寻求颗粒与边界的位置，构成颗粒新的接触，其基本原理如图2所示。
1.2 本构模型
在PFC中，材料的本构关系是通过接触的本构模型来模拟的。颗粒之间的接触本构模型有接触刚度模型、滑动模型和黏结模型。接触刚度模型提供了接触力和相对位移的弹性关系；滑动模型则控制切向和法向接触力的关系，使得接触颗粒可以发生相对移动；而黏结模型是限制总的切向和法向力，使得颗粒在黏结强度范围内保持接触。
图2 颗粒流理论计算循环过程
Fig. 2 Calculation cycle in particle flow theory
黏结模型又分为接触黏结和平行黏结[19]。对于模拟岩石材料，需要设置黏结模型来表征颗粒之间胶结物的存在。接触黏结模型可以看作是1对作用在接触点的具有恒定法向和切向刚度的弹簧；平行黏结可以看作是一组均匀分布在接触面上以接触点为中心的具有恒定法向和切向刚度的弹簧(见图3，其中：Rmax和Rmin分别为两接触颗粒中较大颗粒半径和较小颗粒半径；Kn和Ks分别为颗粒法向和切向接触刚度；
和
分别为平行黏结的切向和法向刚度；Rbond为平行黏结半径；Fn和Fs分别为颗粒之间的法向和切向接触力)。在PFC中，颗粒可以在法向和切向自由运动，并且颗粒之间可以发生旋转。这一旋转会使颗粒之间产生矩，但是接触黏结不能抵抗矩(只能传递力)，而平行黏结既可以传递力又可以传递矩。
图3 平行黏结模型示意图
Fig. 3 Sketch map of parallel bond model
使用接触黏结模型时，由于黏结断裂后只要颗粒还保持接触，接触刚度仍然有效，因此，接触黏结的断裂不会对宏观刚度产生太大影响，这与岩石的破裂机制不符。使用平行黏结模型时，宏观刚度由接触刚度和黏结刚度组成，因此，黏结破裂立即导致宏观刚度的下降。在这个意义上来讲，平行黏结通过在拉伸或剪切断裂时刚度相应降低更逼真地模拟岩石类材料。因此本文中选用的是平行黏结模型。
2 岩石细观物理力学性质参数
由于PFC模拟采用细观物理力学性质参数来表征颗粒及黏结的力学性质且这些细观参数无法从室内试验直接获取。因此，在数值模型进行计算之前，需要对模型的细观参数进行标定。在此过程中，需要进行一系列与室内试验或现场条件类似的模型试验，并将模拟结果与室内试验或原位测试结果进行对比，采用“试错法”反复改变细观参数[18]，直到模型的宏观力学响应满足要求为止。
本文选用黄岛国家石油储备库地下水封石油洞库花岗岩的力学参数(见表1)，室内三轴压缩试验采用TAW-2000电液伺服岩石三轴仪进行。颗粒流程序中通过“试错法”反复调整对比，使得颗粒流模型的宏观力学参数(见表1)接近花岗岩室内试验力学参数，最终确定的细观物理力学性质参数见表2。围压在6 MPa (与室内三轴压缩试验对应) 下试件最终失效形式如图4所示。表1中内聚力和内摩擦角的试验值与模拟值相差较大，这是BMP模型选用圆形颗粒造成的。
表1 花岗岩物理力学参数与PFC模拟值
Table 1 Physico-mechanical parameters and PFC results
表2 花岗岩细观物理力学性质参数
Table 2 Mesoscopic physico-mechanical parameters of granite
图4 围压6 MPa下花岗岩试件受压破坏形态
Fig. 4 Failure modes of granite in compression at confining stress of 6 MPa
3 岩石节理面直接剪切模拟实验
3.1 直接剪切数值试验模型的建立
岩石试件的高×直径为100 mm×36 mm，孔隙率设置在0.16[20]。CNL边界条件下的剪切盒由8块刚性墙体组成(见图5(a))，其中1号、2号、3号和8号墙体组成下剪切盒，4号、5号、6号和7号墙体组成上剪切盒；CNS边界条件下的剪切盒中用clump墙体替代CNL边界条件下的6号墙体(见图5(b))。图5中：v为剪切方向施加的速度，σ为施加的法向应力。试样生成之后，初始化试件内部应力场至0.1 MPa。
图5 不同边界条件下岩石节理直接剪切试验数值模型
Fig. 5 Numerical models of joint rock direct shear test under different boundary conditions
通过剖面线的点坐标文件建立连续的墙体来还原节理面轮廓线，将与连续墙体接触的颗粒进行标记并定义为组，改变组的细观参数，从而实现节理面。将节理面的平行黏结强度设置为0。
3.2 边界条件的实现
在PFC中，模型边界有2种形式，一种是wall刚性边界，另一种是颗粒集合体边界。由于wall刚性边界不能参与运动方程的计算，所以，荷载只能以速度的方式施加，而不能直接施加力或力矩，但可以通过借助FISH函数的伺服机制[18]来实现力的施加。颗粒集合体边界的一种方式是将颗粒定义为组(group)，然后直接对其施加速度、力或力矩，但这样会存在边界变形问题，无法实现刚性边界；另一种是将边界颗粒定义为clump，既能对其直接施加速度、力或力矩，又可将其运动看作是单一刚体运动。
1) 恒定法向荷载边界条件。通过伺服机制调整墙体的速度来实现恒定荷载的施加，从而实现恒定法向荷载边界条件(CNL)。其基本原理如下。
首先给定法向荷载σrequired，然后通过监测6号墙体上不平衡力获得其上的实际应力σmeasured。然后通过下式换算该墙体上所需施加的速度
为
(1)
式中：G是伺服控制参数。
每一时步内由6号墙体移动所引起的墙体不平衡力的最大增量为
(2)
式中：Nc为6号墙体上的接触数；kn[w]为上述接触的平均刚度。因此，6号墙体上平均应力的改变量为
(3)
式中：A为6号墙体的面积。为了模型计算过程中的稳定性，Δσ[w]的绝对值须小于Δσ的绝对值。因此，引入松弛系数α(取值在0~1范围内)，可得
(4)
将式(1)和式(3)代入式(4)并整理可得伺服控制参数G：
(5)
再将式(5)代入式(1)中即可得到施加在6号墙体上的垂直加载速度，从而实现恒定荷载的施加。
2) 恒定法向刚度边界条件。由于恒定法向刚度边界条件(CNS)要求施加的荷载随法向变形量的变化时刻改变，因此荷载是随时间变化的量，在PFC中需要借助FISH循环语句来实现。若使用wall刚性边界，由于伺服机制也是借助FISH循环语句施加，从而增加计算的复杂性，且在2个循环相互嵌套的过程中无法预测对结果造成的影响，因此，选用clump刚性边界来实现恒定法向边界条件。其基本原理为
(6)
式中：σ(t+Δt)和σ(t)为施加的法向应力；K为节理岩体围岩法向刚度；Δy为节理面的法向位移。深埋洞室节理围岩刚度的计算公式如下[21]：
(7)
式中：E为岩石弹性模量；υ为岩石泊松比；r为受影响的围岩范围。
将clump刚性边界的不平衡力实时反馈至式(1)从而实时改变施加在clump上的荷载，从而实现恒定法向刚度边界条件。
3.3 直接剪切数值试验的实现
在剪切过程中，下剪切盒固定不动，4号墙体作为剪切加载墙，其水平方向位移量作为剪切位移，水平方向不平衡力除以节理面的水平投影面积作为剪切应力。在CNL边界条件下，6号墙体用伺服控制保证恒定法向应力，其垂直方向位移量作为法向位移。在CNS边界条件下，使用clump刚性边界替代6号刚性墙体来实现恒定法向刚度，其垂直方向位移量作为法向位移，垂直方向不平衡力除以节理面的水平投影面积作为法向应力。
在试验过程中，需要保持整个过程为拟静力加载状态，这样才能使系统能够及时对试件内部产生的不平衡力重新进行调整，保证结果的可靠性[17]。因此，经过不断调试，本文最终确定选用的水平方向加载速率为0.05 m/s，即相应于试验中的4.25×10-8 s/步，剪切1.0 mm需要加载493 300时步，这一加载速率已足够使整个剪切过程保持在拟静力的状态下。
3.4 直接剪切模拟实验值与经验公式对比
Barton等[22]根据大量的试验，在统计分析的基础上提出了CNL边界条件下的经验公式，同时考虑了正应力和不规则结构面表面特征对剪切强度的影响，是目前应用最广泛的强度公式。Barton结构面抗剪强度τ的经验公式如下：
(8)
式中：σn为作用于节理面上的正应力；JRC为节理面粗糙度系数；JCS为节理面壁的抗压强度；φb为基本内摩擦角。JCS与节理面的风化程度有关，未风化节理面的抗压强度可选为完整岩石的单轴抗压强度[23]。根据模拟结果得到岩石的单轴抗压强度为102 MPa，基本内摩擦角φb选取经验值37°[24]。节理面选取Barton提出的粗糙面曲线[22]中的5条作为研究对象，分别是粗糙度系数JRC为2~4，6~8，12~14，16~18和18~20的粗糙面曲线(见图6)。
CNL边界条件下数值模拟结果与JRC-JCS经验模型预测值的对比如图7所示。对比模拟值与经验公式预测值可以看出：在粗糙度系数JRC较小时，模拟值与预测值相差相对较大，这主要是PFC中颗粒模拟节理面时本身具有一定的粗糙性造成的，粗糙度系数JRC越小，颗粒模拟节理面本身的粗糙性表现的越明显。总体上岩石节理面的数值模拟结果与JRC-JCS经验公式预测结果具有良好的一致性，从而验证了模型的正确性。
图6 计算中采用的粗糙节理剖面图
Fig. 6 Rough joint profiles used in calculation
图7 Barton经验公式预测值与PFC模拟值对比
Fig. 7 Comparison of PFC simulation results with calculated values by Barton’s empirical equation
3.5 不同边界条件下直接剪切数值试验结果分析
参照表2中花岗岩的力学参数，选取受影响围岩的范围为2 m，根据式(6)求得法向刚度为12 GPa/m。对粗糙度系数JRC为2~4，6~8和12~14的节理面在不同边界条件、初始法向应力为10，12，15，18和20 MPa作用下进行直接剪切模拟研究，图8~10列出了部分模拟结果。
在CNL边界条件下，当法向应力σn为10，12，15，18和20 MPa时，岩石节理面的剪切特征具体表现为：1) JRC为2~4的岩石节理面的峰值剪切应力分别为13.11，14.13，16.64，18.07和18.87 MPa，相应的剪切位移分别为0.337 5，0.360 7，0.405 5，0.442 5和0.452 9 mm，残余强度分别为0.30，0.34，0.40，0.42和0.45 MPa；2) JRC为6~8的岩石节理面的峰值剪切应力分别为14.95，16.05，18.27，19.77和20.74 MPa，相应的剪切位移分别为0.367 6，0.392 9，0.447 1，0.478 3和0.491 5 mm，残余强度分别为0.30，0.35，0.44，0.50和0.60 MPa；3) JRC为12~14的岩石节理面的峰值剪切应力分别为15.72，16.93，18.98，21.12和22.25 MPa，相应的剪切位移分别为0.353 7，0.378 0，0.424 2，0.465 0和0.487 3 mm，残余强度分别为0.35，0.40，0.55，0.59和0.65 MPa。
对比上述数据可以看出：随着恒定法向应力增加，岩石节理面的抗剪强度及其对应的剪切位移以及残余强度都相应增加；抗剪强度随粗糙系数增加会明显增加。
分析法向位移-剪切位移曲线可以看出：在剪切初始阶段，法向位移变化较小，且为负值(取节理相对的2个面垂直移向对方为负)，表明节理面处于压密阶段；随着剪切过程的进行，法向位移开始由负变为正，且随剪切位移的增加逐渐增加，表明节理面发生剪胀。对比不同法向应力作用下的法向位移-剪切位移曲线可以看出：法向应力越大，剪切初始阶段的压密效应越明显；在剪胀阶段，法向应力越低，法向位移增加的速率越快，法向膨胀量越大。
在CNS边界条件下，当法向应力σn为10，12，15，18和20 MPa时，岩石节理面的剪切特征具体表现为：1) JRC为2~4的岩石节理面的峰值剪切应力分别为12.26，13.60，16.72，18.33和19.77 MPa，相应的剪切位移分别为0.278 1，0.302 0，0.371 3，0.408 2和0.430 2 mm，残余强度分别为0.35，0.40，0.43，0.47和0.50 MPa；2) JRC为6~8的岩石节理面的峰值剪切应力分别为15.71，17.07，19.02，21.39和22.41 MPa，相应的剪切位移分别为0.341 0，0.370 3，0.407 6，0.459 6和0.485 4 mm，残余强度分别为0.45，0.50，0.55，0.60和0.70 MPa；3) JRC为12~14的岩石节理面的峰值剪切应力分别为15.99，17.67，21.01，21.39和21.95 MPa，相应的剪切位移分别为0.315 9，0.378 9，0.425 9，0.473 8和0.440 6 mm，残余强度分别为0.65，0.70，0.73，0.77和0.80 MPa。
图8 JRC=2~4时CNL和CNS边界条件下不同粗糙节理面的剪切特性
Fig. 8 Shear behavior of joints with different roughnesses under CNL and CNS boundary conditions when JRC=2-4
图9 JRC=6~8时CNL和CNS边界条件下不同粗糙节理面的剪切特性
Fig. 9 Shear behavior of joints with different roughnesses under CNL and CNS boundary conditions when JRC=6-8
图10 JRC=12~14时CNL和CNS边界条件下不同粗糙节理面的剪切特性
Fig. 10 Shear behavior of joints with different roughness under CNL and CNS boundary conditions when JRC=12-14
从CNS边界条件下的试验结果可以看出：剪切应力、法向位移随剪切位移的变化趋势与CNL边界条件下的基本一致；但CNS边界条件下法向应力随法向位移增加而增大，其增加量会随初始法向应力增大而相应降低，因而法向位移的增加量会相应减小。
综合对比CNL和CNS 2种边界条件下的剪切试验结果可以看出：CNL边界条件下的峰值剪切应力要低于CNS边界条件下的峰值剪切应力；CNS边界条件下的残余强度较CNL边界条件下会有所提高。由法向应力、法向位移-剪切位移对比曲线可以看出：CNS边界条件下由于法向应力增加，法向位移与CNL边界条件下相比，会明显的降低。
图11所示为不同粗糙度下的峰值剪切应力随初始法向应力的变化趋势。由图11可以看出：节理面的摩擦角随粗糙度系数增加而增加，内聚力同样随粗糙度系数增加而增加。边界条件的变化只是影响两者数值大小，而对整体变化趋势影响不大。
图11 峰值剪切应力与法向应力关系曲线
Fig. 11 Relationship between peak shear stress and initial normal stress
4 试件细观力学响应及破坏机理
以粗糙度系数JRC为12~14的岩石节理面为例详细研究岩石节理面在不同边界条件下直接剪切过程中的细观力学响应及破坏机理。
4.1 粒间接触力分布及细观裂纹发育特征
试验过程中对试件内部的接触力分布和裂纹的发育分布进行了动态跟踪。图12和图13所示为在不同边界条件、初始法向应力为10，12，15，18和20 MPa作用下，剪切位移达到10.0 mm时试件内部的接触力分布和裂纹发育状况。
由图12和图13可以看出：接触力集中现象主要发生在节理面凸起位置，且以接触压力为主。正因为如此，细观裂纹也主要分布在节理面两侧，主要聚集在凸起的位置且张拉裂纹占绝对优势。这一点符合“压致张拉裂纹”的力学机制[20, 25]。节理面凸起位置的接触力集中程度随初始法向应力的增加而加剧且范围加大；裂纹数量也随初始法向应力的增加迅速增加，且分布范围向节理面的延伸方向和试件的纵深方向扩展。
对比同一初始法向应力、不同边界条件下试件内部粒间接触力及裂纹分布可以看出：CNS边界条件下由于法向应力随法向位移增加而增大，所以，接触力集中程度加剧，从而裂纹的数量也明显增加且纵深向扩展范围加大。在剪切过程中法向应力不断变化，其微观接触状态也不断变化，从而导致产生的裂纹纵深向扩展方向与CNL边界条件下的扩展方向有所不同。
图12 CNL边界条件时剪切位移达到10.0 mm时试件内部接触力及裂纹分布
Fig. 12 Contact force and micro-cracks distribution at shear displacement of 10.0 mm under constant normal load
图13 CNS边界条件时剪切位移达到10.0 mm时试件内部接触力及裂纹分布
Fig. 13 Contact force and micro-cracks distribution at shear displacement of 10.0 mm under constant normal stiffness
4.2 粒间接触力及细观裂纹演化特征
图14所示为不同边界条件、20 MPa初始法向应力作用下剪切应力、法向位移、裂纹数目和破断频数随剪切位移的演化曲线。为了研究剪切过程中试件内部粒间接触力和细观裂纹的演化特征，在图14中设置编号为A(a)~G(g)的7个监测点，分别对应于剪切位移0(0)，0.26(0.18)，0.35(0.32)，0.49(0.44)，1.00(0.44)，5.00(0.44)和10.00 mm(图中未标注)。图15和图16所示分别为监测点A(a)至F(f)对应的粒间接触力和细观裂纹分布情况。
由图14可知：剪切初始阶段，剪切应力曲线基本呈线弹性，且此阶段没有裂纹产生，破断频数为0。剪切应力达到B(b)点(宏观上的启裂强度)时，试件开始产生细观裂纹，在应力到达屈服点C(c)之前，破断频数一直比较稀疏，振幅随着裂纹的张开、扩展有所增大，表明试件进入稳定破裂传播阶段。当剪切应力从C(c)点发展到D(d)点时，剪切应力曲线呈非线弹性，并伴有小幅度突降，直到达到峰值剪切应力。此阶段裂纹数量明显增加，且破断频数呈现快速增长趋势，表明试件从稳定破裂传播阶段转向非稳定破裂传播阶段。剪切应力超过峰值D(d)点以后，试件表现出脆性发展特征，剪切应力大幅度下降，裂纹数量急剧增加，破断频数出现最大值。随后剪切应力曲线平稳发展，裂纹数量保持恒定的速率增长，破断频数降至较低水平，表明岩石又恢复了稳定破裂的发展态势。
图14 不同边界条件和20 MPa初始法向应力作用下剪切应力、法向位移、裂纹数量与破断频数演化曲线
Fig. 14 Evolutions of curves of shear stresses, normal displacements, micro-cracks and rupture frequencies at initial normal stress of 20 MPa under different boundary conditions
由图15和图16可以看出：剪切位移为0 mm，即点A(a)处，试件完成初始法向荷载加载，由于只有法向荷载作用，此时粒间接触力是压力，且沿竖直方向均匀分布。此阶段没有细观裂纹产生。当试件剪切至B(b)点时，节理面附近开始出现裂纹。在CNS边界条件下，由于试件上部clump加载板刚度大的原因，与加载板接触的颗粒之间产生少量的裂纹。由于受剪切荷载的影响，粒间接触力的方向发生明显偏转。4号墙体作为施力边界，造成上部试件左边接触力明显大于右边接触力，且接触力的方向从左往右逐渐由水平方向向节理面偏转，节理面上的接触力向荷载加载端偏转。由于试件下部固定，所以，其粒间接触力分布与上部恰好相反。在此阶段，节理面附近开始出现接触压力集中现象，但是不明显，CNL边界条件和CNS边界条件下最大接触力分别为6.986×104 N和7.048×104 N。
随着剪切的进行，试件内部接触力偏转现象更加明显。当剪切经过C(c)点和D(d)点，节理面附近接触压力集中程度加剧，最大接触力由C(c)点的1.260×105 N(CNL)和1.468×105 N(CNS)增加到D(d)点的1.694×105 N(CNL)和1.847×105 N(CNS)。受不断增大的接触力作用，节理面凸起发生破坏，产生大量的裂纹。当剪切至峰值D(d)点，裂纹基本贯穿整个节理面，此现象在CNS边界条件下更加明显。
图15 CNL边界条件和20 MPa初始法向应力作用下粒间接触力和细观裂纹分布
Fig. 15 Contact force and micro-cracks distribution at initial normal stress of 20 MPa under constant normal load
图16 CNS边界条件和20 MPa初始法向应力作用下粒间接触力和细观裂纹分布
Fig. 16 Contact force and micro-cracks distribution at initial normal stress of 20 MPa under constant normal stiffness
当剪切超过D(d)点以后，剪切进入峰后阶段。此时较大的凸起已被剪断，接触压力的集中位置逐渐转移至次一级的凸起。因此，在CNL边界条件下，在剪切峰后的D，E和F点及其以后更大剪切位移点的接触力会出现波动上升的趋势(如表3所示)。但是，在CNS边界条件下，由于法向应力的不断变化，峰后接触力不会出现波动上升的趋势，而是持续的增加(如表4所示)，这也说明了CNS边界条件下细观裂纹数较多的原因。
随着剪切位移的逐渐加大，节理面上凸起大量被剪断压碎，从而沿节理延伸方向产生大量裂纹，同时逐渐向试件的纵深方向扩展，最终形成1条较明显的宏观剪切破碎带。
表3 20 MPa恒定法向力作用下粒间最大接触力
Table 3 The maximum contact forces at normal stress of 20 MPa under CNL boundary condition
表4 20 MPa初始法向力作用下粒间最大接触力
Table 4 The maximum contact forces at initial normal stress of 20 MPa under CNS boundary condition
5 结论
1) 模拟值和JRC-JCS模型预测值具有良好的一致性，验证了模型的正确性。
2) 在CNL边界条件下，岩石节理的抗剪强度及其对应的剪切位移以及残余强度随恒定法向应力的增加而提高。在剪切过程中，试件在垂直于节理面方向上经历压密和膨胀2个阶段，且压缩量和膨胀量受法向应力的影响。
3) CNL边界条件的抗剪强度和残余强度比CNS边界条件的低，而法向位移比CNS边界条件的高。在CNS边界条件下，法向应力不断增加导致粒间接触力集中程度要明显比CNL边界条件的高，因而裂纹数量也明显增加且裂纹的扩展方向有所不同。在CNS边界条件下剪切的残余阶段，粒间最大接触力由于法向应力的不断变化，不会出现CNL边界条件下的波动上升趋势，而是持续增加。
4) 接触力集中现象主要发生在节理凸起位置，且以接触压力为主。正因为如此，细观裂纹也主要聚集在凸起的位置且张拉裂纹占绝对优势。接触力集中程度随初始法向应力的增大而加剧且范围加大。裂纹数量也随法向应力增加迅速增加，且分布范围向节理的延伸方向和试件的纵深方向扩展。
5) 剪切初始阶段试件内部不存在裂纹，在弹性阶段剪切到一定位移时，会首先在较大凸起位置产生裂纹。此后，裂纹的扩展主要经历平稳发展、急剧增加和平稳发展3个阶段，表明试件主要经历弹性阶段、稳定破裂传播阶段、非稳定破裂传播阶段以及残余阶段的稳定破裂传播阶段。
6) 受剪切应力和法向应力的影响，在剪切过程中粒间接触力发生重分布，且主要偏向施力边界。凸起位置接触压力集中，导致产生细观裂纹。剪切进一步发展，较大凸起被剪断，接触力逐渐向次一级凸起转移，从而产生沿节理延伸方向产生大量裂纹，同时逐渐向试件的纵深方向扩展，最终形成1条较明显的宏观剪切破碎带。
参考文献：
[1] Patton F D. Multiple modes of shear failure in rock [C]// Proceedings of the 1st International Congress on Rock Mechanic. Lisbon: International Society for Rock Mechanics, 1966: 509-513.
[2] Ladany B, Archambault G. Simulation of shear behavior of a jointed rock mass[C]//Proceedings of 11th US Symposium on Rock Mechanics, Berkeley: American Rock Mechanics Association, 1970: 105-123.
[3] Barton N, Choubey V. The shear strength of rock joints in theory and practice[J]. Rock Mechanics, 1977, 10: 1-54.
[4] Plesha M E. Constitutive models for rock discontinuties with dilatancy and surface degradation[J]. International Journal for Numerical and Analytical Methods in Geomachanics, 1987, 11: 345-362.
[5] 夏才初, 宋英龙, 唐志成, 等. 反复直剪试验节理强度与粗糙度变化的研究[J]. 中南大学学报(自然科学版), 2012, 43(9): 3589-3594.
XIA Caichu, SONG Yinglong, TANG Zhicheng, et al. Shear strength and morphology characteristic evolution of joint surface under cyclic loads[J]. Journal of Central South University (Science and Technology). 2012, 43(9): 3589-3594.
[6] 曹平, 范祥, 蒲成志, 等. 节理剪切试验及其表面形貌特征变化分析[J]. 岩石力学与工程学报, 2011, 30(3): 408-485.
CAO Ping, Fan Xiang, PU Chengzhi, et al. Shear test of joint and analysis of morphology characteristic evolution of joint surface[J]. Chinese Journal of Rock and Engineering, 2011, 30(3): 408-485.
[7] 肖卫国, 兑关锁, 陈铁林, 等. 剪胀和破坏耦合的节理岩体本构模型的研究[J]. 岩石力学与工程学报, 2009, 28(12): 2535-2543.
XIAO Weiguo, DUI Guansuo, CHEN Tielin, et al. A constitutive model coupling of dilatancy and degradation with rock joints[J]. Chinese Journal of Rock and Engineering, 2009, 28(12): 2535-2543.
[8] Brahim B, Gerard B. Laboratory of shear behavior of rock joints under constant normal stiffness conditions[C]// Khair A W. Rock Mechanics as a Guide for Efficient Utilization of Natural Resources. Rotterdam: Balkema, 1989: 899-906.
[9] Indraratna B, Haque A, Aziz N. Laboratory modeling of shear behavior of soft joints under constant normal stiffness conditions[J]. Geotechnical and Geological Engineering, 1998, 16: 17-44.
[10] Jiang Y, Xiao J, Tanabashi Y, et al. Development of an automated servo-controlled direct apparatus applying a constant normal stiffness condition[J]. International Journal of Rock Mechanics & Mining Sciences, 2001, 41: 275-286.
[11] Selvadural A P S, Yu Q. Mechanics of a discontinuity in a geomaterial[J]. Computers and Geotechnics, 2005, 32(2): 92-106.
[12] 王长军, 杨远志, 李勇会, 等. 接触面力学性能研究及数值分析[J]. 水利与建筑工程学报, 2008, 6(4): 74-76.
WANG Changjun, YANG Yuanzhi, LI Yonghui, et al. Research and numerical analysis for mechanical behavior of rough interface[J]. Journal of Water Resources and Architectural Engineering, 2008, 6(4): 74-76.
[13] Beer G, Poulsen B A. Efficient numerial modeling of faulted rock using the boundary element method[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1997, 31(5): 485-506.
[14] Park J W, Song J J. Numerical simulation of a direct shear test on a rock joint using a bonded-particle model[J]. International Journal of Rock Mechanics & Mining Sciences, 2009, 46: 1314-1328.
[15] 夏才初, 宋应龙, 唐志成, 等. 粗糙节理剪切性质的颗粒流数值模拟[J]. 岩石力学与工程学报, 2012, 31(8): 1545-1552.
XIA Caichu, SONG Yinglong, TANG Zhicheng, et al. Particle flow numerical simulation for shear behavior of rough joints[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(8): 1545-1552.
[16] 周喻, Misra A, 吴顺川, 等. 岩石节理直剪试验颗粒流宏细观分析[J]. 岩石力学与工程学报, 2012, 31(6): 1245-1256.
ZHOU Yu, Misra A, WU Shunchuan, et al. Macro-and meso-analyses of rock joint direct shear test using particle flow theory[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(6): 1245-1256.
[17] Asadi M S, Rasouli V, Barla G. A bonded particle model simulation of shear strength and asperity degradation for rough rock fracture[J]. Rock Mechanics and Rock Engineering, 2012, 45: 649-675.
[18] Itasca Consulting Group Inc. Manual of particle flow code in 2-dimension (Version 3.10)[M]. Minneapolis: Itasca Consulting Group Inc, 2005.
[19] Cho N, Martin C D, Sego D C. A clumped particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2007, 44: 997-1010.
[20] Potyondy D O, Cundall P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2004, 41: 1329-1367.
[21] Jiang Y, Tanabashi Y, Mizokami T. Shear behavior of joints under constant normal stiffness conditions[C]//Proceedings of the Second Asian Rock Mechanics Symposium(ISRM2001-2nd ARMS). Beijing, 2001: 247-250.
[22] Barton N, Choubey V. The shear strength of rock joints in theory and practice[J]. Rock Mechanics, 1977, 10(1/2): 1-54.
[23] 赵坚. 岩石节理剪切强度的JRC-JMC新模型[J]. 岩石力学与工程学报, 1998, 17(4): 349-357.
ZHAO Jian. A new JRC-JMC shear strength criterion for rock joint[J]. Chinese Journal of Rock Mechanics and Engineering, 1998, 17(4): 349-357.
[24] 叶金汉, 郗绮霞, 夏万仁, 等. 岩石力学参数手册[M]. 北京: 水利电力出版社, 1991: 411.
YE Jinhan, XI Qixia, XIA Wanren, et al. Rock mechanical parameter manual[M]. Beijing: Water Conservancy and Electric Power Press, 1991: 411.
[25] Cundall P A, Potyondy D O, Lee C A. Micromechanics-based models for fracture and breakout around the mine-by tunnel[C]//Proceedings of the Excavation Disturbed Zone Workshop, Designing the Excavation Disturbed Zone for a Nuclear Repository in Hard Rock. Toronto: Canadian Nuclear Society, 1996: 113-122.
(编辑 杨幼平)
收稿日期：2014-04-06；修回日期：2014-06-17
基金项目(Foundation item)：国家自然科学基金资助项目(51279097，51479108，51379117，51009086)；山东科技大学研究生创新基金资助项目(YC140357)；山东省优秀中青年科学家科研奖励基金资助项目(BS2010HZ015)(Projects (51279097, 51479108, 51379117, 51009086) supported by the National Natural Science Foundation of China; Project (YC1403517) supported by the Graduate Innovation Foundation of Shandong University of Science and Technology; Project (BS2010HZ015) supported by the Scientific Research Foundation for Outstanding Young Scientists of Shandong Province, China)
通信作者：王刚，博士，副教授，从事岩石力学与工程方面研究；E-mail：wanggang1110@gmail.com
