节理的广泛存在会导致岩体的强度和稳定性显著降低，众多的岩体工程失稳问题与节理剪切滑移相关。因此，节理剪切特性一直是矿山、隧道、边坡以及水利水电等工程中研究的热点问题[1]。目前，国内外学者已在节理剪切强度、表面破坏特征、抗剪强度与表面形貌之间的关系以及岩石类型、边界条件和加载方式等方面开展了诸多研究[2-7]，对恒剪切速率和恒剪切应力条件下的节理剪切特性有了深入的认识。然而，在地震作用和开挖扰动等情况下节理会承受循环载荷[8]。在循环载荷作用下，节理会不断的滑移和闭合，从而导致其承载性能不断劣化[9]，其剪切行为也变的更为复杂[10]。可见，研究节理的循环剪切特性对于准确评价岩体工程的长期稳定性等具有重要意义。
目前，已有学者针对节理循环剪切特性开展了研究。LEE等[11] 进行了不同粗糙度节理的循环剪切试验，基于试验结果建立了考虑二阶微凸体退化的弹塑性本构模型。LI等[10]提出一种循环剪切作用下考虑节理粗糙退化和岩屑充填的节理剪切本构模型，并进行了试验验证。刘博等[12]开展了不同起伏角、不同试件强度和不同法向应力条件下的节理循环剪切试验，并分析了这些因素对剪胀角和剪切刚度的影响。NIKTABAR等[13]考虑法向刚度的影响，进行了恒定法向刚度和恒定法向应力边界条件下的循环剪切试验。MIRZAGHORBANALI等[14-15]也开展了一系列恒定法向刚度边界条件下的循环剪切试验，探究初始法向应力和剪切速率的影响。受限于现有试验设备和数值模拟方法，以上研究大都是采用位移控制剪切循环，即保持恒定剪切速率在一定剪切位移范围内往复剪切。然而，通过分析地震加速度波可发现，节理剪切速率是随时间不断变化的。为更好地还原节理在循环载荷作用下的受力状态，一些学者选择对节理施加循环剪切应力。DANG等[16-17]采用其研发的新型剪切试验装置，对光滑节理试件进行了法向循环应力条件下的剪切实验，并结合数值模拟对光滑节理的摩擦响应进行分析。杨小彬等[18]设计了恒定围压条件下的循环加载双剪试验，并利用CCD相机和数字散斑方法对试件滑移变形进行监测。但上述研究并未考虑节理粗糙度的影响。FATHI等[19]开展了峰前循环剪切试验，并通过自编数值模拟方法追踪了剪切过程中的节理表面接触特征，进一步阐明了循环剪切荷载对节理剪切机制的影响。刘新荣等[20-23]对含二阶起伏体的均质砂岩进行峰前循环加载剪切试验，并进行了数值模拟研究，从宏细观角度探究循环加载次数与加载幅值对节理剪切特性的影响规律。但上述研究主要针对弱地震作用情况，且未考虑加载频率的影响。综上可见，目前对于循环剪切载荷作用下岩石节理的宏细观剪切特性尚缺乏系统深入的研究。
鉴于以上认识，笔者首先开展节理剪切试验以标定模拟参数和确定循环载荷特征参数，然后进行了不同加载幅值、不同加载频率等循环剪切载荷作用下的粗糙节理剪切数值模拟，研究了岩石节理的宏细观剪切特性以及循环载荷特征参数对节理剪切特性的影响。研究成果可为节理岩体的长期稳定性控制设计等提供一定参考。
1 循环剪切的数值实现及模拟方案
1.1 节理剪切室内试验
为标定模拟参数和确定循环载荷特征参数，首先开展了节理剪切室内试验。试验中采用类岩石材料制备节理试件，该材料的质量配比为高强石膏∶水∶缓凝剂=1∶0.2∶0.005，单轴抗压强度和弹性模量分别为 38.29 MPa 和46.55 GPa。试件的尺寸为200 mm(长)×100 mm(宽)×100 mm(高)，节理的表面形貌通过3D激光扫描仪获得，其粗糙度系数为9.86，如图1所示。试件制作和节理面粗糙度计算过程详见文献[24-25]。
图1 类岩石材料制备节理试件
Fig.1 Rock-like material preparation of jointed specimens
试验采用如图2所示的日本长崎大学的MIS-233-1-55-03型法向应力伺服控制剪切试验系统开展[26]。试验中采用CNL边界条件，由于地下巷道、隧道等开挖后节理围岩的应力会向深部转移，支护范围内的围岩应力远低于原始地应力，根据相关经验选取法向应力为2 MPa。试验过程为：首先将节理面试件放置于特制的剪切盒内并将剪切盒安装至剪切试验机上，之后以0.5 MPa/min的加载速度在剪切盒上方施加初始法向应力至给定值2 MPa，然后以0.5 mm/min的速度在下部剪切盒上施加剪切力推动下半部分实现剪切错动，最大剪切位移设置为3 mm。剪切试验系统的指标参数及试验方法详见文献[26]。
图2 岩石节理剪切试验系统
Fig.2 Rock joint shear test system
1.2 数值模拟试验方案
1.2.1 数值模型的建立
在PFC中建立如图3所示粗糙节理数值模型，其尺寸为200 mm(长)×100 mm(高)，节理粗糙度与室内试验模型的一致。建模过程如下：首先建立4面墙体模拟剪切盒，并在剪切盒范围内按标定的细观参数随机生成颗粒；然后将粗糙曲线导入模型形成粗糙节理面，并对节理面上下部分颗粒进行分组，最后按照结构面坐标重新生成6面墙体。该建模方法能够保证节理面粗糙体不会在建模过程中因颗粒运动而错位。颗粒间的接触模型参照文献[20]选用平行黏结模型，节理处的接触选用光滑节理模型。
图3 粗糙节理数值模型
Fig.3 Numerical simulation model diagram
数值模拟中施加的模型边界条件与室内试验相同，具体如下：限制上部试件沿水平方向和下部试件沿法向方向的移动，即限制1号、2号、3号墙体的水平位移和6号墙体的法向位移，通过伺服函数控制3号墙体的法向速度在试件顶部施加恒定法向应力，对4号、5号和6号墙体设定相同的水平速度实现剪切载荷的施加。
1.2.2 物理力学参数标定
室内试验中无法直接获得颗粒接触模量、黏聚力、抗拉强度等细观参数，建立颗粒细观参数与岩石材料宏观力学参数之间的匹配关系至关重要。因此，结合室内单轴压缩试验和粗糙节理剪切试验的结果，采用 “试错法[27]”对试件细观参数进行了标定，模拟结果如图4所示，确定的参数见表1、2。由图4可知，数值模拟试验与室内试验结果具有较好的一致性，其中因天然岩体存在孔隙、裂隙导致最初的压密阶段有一定差别。
图4 室内试验与数值模拟曲线对比
Fig.4 Comparison of indoor test and numerical simulation test curves
表1 数值模型中颗粒和胶结的细观参数
Table 1 Microscopic parameters of the spherical particles and bonds in numerical model
圆形颗粒胶结参数数值参数数值弹性模量/GPa5.12弹性模量/GPa5.12法向切向刚度比1.20法向切向刚度比1.20摩擦因数0.45拉伸强度/MPa3.00最大半径/mm0.50剪切强度/MPa1.50最大最小半径比1.5摩擦角/(°)20密度/(kg·m-3)2 780
表2 数值模型中光滑节理的细观参数
Table 2 Microscopic parameters of smoothrock joints in numerical model
法向刚度/(GPa·m-1)剪切刚度/(GPa·m-1)摩擦因数抗拉强度/MPa黏聚力/MPa2002000.7500
1.2.3 循环剪切加载实现
循环剪切载荷的施加是数值模拟的关键，本文借助FISH语言将循环剪切应力路径嵌入伺服函数，主要包括3个方面：
其一为目标循环剪切应力的表达，使用FISH语言将循环剪切应力路径通过数学公式表达，本文采用的余弦波表达式为
τ=A+Bcos(2πft)
(1)
式中，τ为施加的循环剪切应力；A为初始剪切应力；B为加载幅值；f为加载频率；t为时间。
其二为节理面剪切应力的监测，主要通过颗粒与墙之间的水平载荷计算得出，具体为：分别监测试件左侧4号和右侧5号墙体与颗粒间的水平荷载，然后取平均值得出节理面上的剪切载荷，再除以节理面长度得到剪切应力。
其三为伺服控制，首先通过计算得出伺服调整系数G(计算过程详见PFC手册[28])，然后根据施加的循环剪切应力与节理实际所受剪切应力计算应力差Δσ，根据伺服调整系数和应力差计算得到4号、5号和6号墙体的速度，为
(2)
式中，
为3面墙体的速度；σm为监测到的剪切应力。
1.2.4 数值模拟方案
由第1.1节的剪切试验结果可知，本文所选用节理的峰值剪切应力为2.35 MPa，为避免节理直接发生失稳破坏，则模拟中循环剪切应力的加载幅值不能高于该值。地震波的频率一般在10 Hz以下，相关研究[16,29]中加载频率一般取0.1～2.0 Hz。因此，为探究加载幅值和加载频率2个特征参数对节理剪切特性的影响，本文依照室内剪切试验结果以及相关文献中的取值，选取了表3的10组循环剪切载荷特征参数在法向应力为2 MPa条件下进行数值模拟试验。
表3 循环剪切数值模拟方案
Table 3 Numerical simulation scenarios for cyclic shear
试验编号频率/Hz加载幅值/MPa1～50.1、0.5、1.0、1.5、2.01.76～1011.5、1.7、1.9、2.1、2.3
2 循环载荷下节理宏细观剪切特性
限于篇幅，仅以法向应力2 MPa时循环剪切载荷加载幅值1.7 MPa、加载频率1的工况为例，对岩石节理宏细观剪切力学特性进行分析。
2.1 剪切位移-剪切应力特征
图5为恒定剪切速率和循环剪切载荷试验的剪切位移-剪切应力曲线。由图5可以看出，随着剪切载荷的循环加卸，节理剪切应力-剪切位移曲线出现明显的滞回效应，形成多个滞回环，这是由于循环剪切载荷的施加导致节理面损伤，产生了不可恢复的塑性应变。由图5还可知，在不同时期，剪切应力曲线表现出的特征不同，大致可分为4个阶段：第Ⅰ阶段在直剪曲线峰值前，处于弹性区，为第1个循环，该阶段的剪切位移较大，滞回环较小；第Ⅱ阶段剪切位移较小，在直剪曲线峰值附近，与第Ⅰ阶段相同，该阶段滞回环较小，但该阶段曲线密集；第Ⅲ阶段处于直剪曲线的软化阶段，该阶段内形成了明显的滞回环，曲线较为疏松，且略低于目标剪切应力，这是由于该阶段节理已发生失稳破坏，节理抗剪强度下降，再加之循环剪切应力的频率为1，因此节理实际所受剪切应力难以达到目标应力；第Ⅳ阶段处于直剪曲线的残余阶段，该阶段内的滞回环较大，曲线非常疏松，同样略低于目标剪切应力。并且该阶段内部分剪切应力曲线高于直接剪切时的剪切应力曲线，这是主要是由于节理粗糙度退化以及循环剪切应力和恒定剪切速度2种加载方式剪切速率不同造成的，当进行循环剪切试验时采用的是伺服控制，会不断调整剪切速度来尽力达到目标剪切应力，而直接剪切试验则是采用恒定剪切速率进行剪切。
图5 循环剪切过程中剪切应力-剪切位移曲线
Fig.5 Shear stress-shear displacement curves during cyclic shear
2.2 剪切位移和法向位移特征
为探究节理表面粗糙体对节理位移场的影响，以图6所示的剪切初期(第3秒)时和剪切结束时的颗粒位移分布云图为例进行分析。由剪切初期的图6(a)可见，在节理表面粗糙体的作用下，节理的位移并非均匀分布，大致呈以下规律：① 上部节理试件水平方向固定，整体位移几乎为0，仅节理面附近的颗粒发生较小位移，最大值为0.076 mm；② 下部节理试件在剪切力作用下，整体上发生较大位移，但在节理表面粗糙体的作用下，节理面附近颗粒位移相对较小，最小0.11 mm，并且整个位移云图大致以节理中部较大凸起处为中心呈“V”型分布。由剪切结束时的图6(b)可见，此时节理试件的位移差别主要体现在粗糙体附近。这是由于，剪切结束时大部分粗糙体已被破坏，不会阻碍下部试件的剪切移动，即被破坏部分不再传递力和位移，因此除已被破坏的粗糙体外，试件其他部分的位移分布更均匀。
图6 节理剪切初期及结束时位移分布云图
Fig.6 Displacement distribution clouds at the beginning and end of joint shear
图7为剪切过程中剪切应力、剪切位移、法向位移和剪切速度随时间的变化规律。可以看出，随着剪切应力的循环加卸，剪切位移和法向位移曲线都随之波动，且不同阶段曲线波动幅度不同。在第Ⅰ阶段，剪切应力首次作用在节理上，节理产生较大的剪切位移；在第Ⅱ阶段，曲线波动较小，尤其是法向位移曲线，几乎看不出变化，每个循环内，剪应力峰值与位移峰值对应，这是由于该阶段基本处于弹性阶段，并未产生大量的塑性变形。在第Ⅲ阶段，曲线波动逐渐增大，每个循环内，位移峰值偏向于剪切应力峰值右侧，即在部分卸载阶段，节理剪切位移和法向位移也在增大，这是由于结构面开始产生较多的微裂纹，产生了较多损伤造成的，详见第2.3节。到了第Ⅳ阶段，位移曲线斜率进一步增大，且在每个循环内位移曲线明显不对称，即使在卸载阶段时，剪切位移和法向位移也会增大。这是由于在该阶段，节理损伤程度更大，节理的抗剪强度显著下降，因此每个循环内的剪切变形量更大。此外还可看出，在伺服函数的控制下，为达到目标循环剪切应力，剪切速度在不断变化，其波形图同真实地震速率波相似。在不同阶段内剪切速度的峰值不同，在第Ⅲ和第Ⅳ阶段内剪切速率较大，最大值为1.36 mm/min，即节理面损伤大时剪切速度波动大。
图7 剪切应力、剪切位移、法向位移和剪切速度随时间变化规律
Fig.7 Variation of shear stress-shear displacement-normal displacement-shear velocity with time
为更直观的分析循环剪切载荷对节理位移的影响，统计了多个循环结束时的剪切位移和法向位移，如图8所示，由图8可知，整个剪切过程可分为个4阶段，分别对应前文中剪切应力曲线的4个阶段：第Ⅰ阶段，初始阶段，为第1个循环，剪切位移由0增大到0.095 7 mm，该阶段曲线斜率为0.095 7 mm/次；第Ⅱ阶段，稳定阶段，由第2个循环到第25个循环，在该阶段随着循环次数的增加，节理剪切位移变化并不大，由0.095 7 mm到0.288 7 mm，累计增长了0.193 0 mm，曲线平均斜率为0.008 4 mm/次；第Ⅲ阶段，加速阶段，由第25个循环到第35个循环，在这个阶段，随着循环次数的增加，剪切位移增长明显，由0.288 7 mm增大到1.282 1 mm，增大了0.993 6 mm，曲线平均斜率为0.198 7 mm/次；第Ⅳ阶段，失稳阶段，为第35循环以后，这是由于在该阶段，节理面已发生失稳破坏，其表面损伤较大，抗剪强度迅速下降，因此，对循环剪切载荷的响应更剧烈，剪切位移变化更大，为1.717 8 mm，曲线斜率为0.277 1 mm/次。法向方向的变形较小，法向位移在4个阶段内分别增长了0.008 4、0.032 6、0.215 0、0.211 0 mm，曲线平均斜率分别为0.008 4、0.000 2、0.014 3、0.023 4 mm/次，其增长规律与剪切方向相同。
图8 循环次数与剪切位移关系曲线
Fig.8 Number of cycles versus shear displacement curves
2.3 裂纹及力链分布演化特征
模拟中，循环剪切应力的频率为1，即1 s一个循环，前0.5 s为加载阶段，后0.5 s为卸载阶段，如3.0～3.5 s为加载阶段，剪切应力由最小值0.5 MPa到最大值1.7 MPa。为探究循环剪切应力作用过程中裂纹的发育规律，记录了多个循环内加载过程中的裂纹分布情况，如图9所示，其中红色曲线为节理面。整体来看，裂纹主要分布在节理面附近，且裂纹数量随着循环剪切应力的施加不断增多。在30 s以前，裂纹增长较为缓慢，数量较少。半个周期内，即一次剪切应力加载过程中几乎未见裂纹增长；在30 s以后，可明显的看到节理面上分布着大量裂纹，同样的在半个周期内几乎未见裂纹明显增多，但是在36 s以后，可看到在半个周期内，节理面上裂纹明显增加，仅在3个循环后，裂纹遍布在节理面上。
由图9还可以发现,裂纹的分布位置与节理表面粗糙体相关，在36 s以前裂纹大都集中在较小的二阶凸起上，直到第36秒时，大致位于节理面中部的一阶凸起底部出现了裂纹，之后该凸起上产生大量裂纹，最后裂纹分布在整个节理面上。
图9 剪切过程中节理裂纹分布
Fig.9 Distribution of nodal cracks during shearing
图10为不同时刻节理试件的接触力分布，其中接触的方向代表受力方向，接触的粗细反映了受力的大小，黑线表示节理面。可以看出，随着循环剪切的进行，节理受力逐渐向节理面上集中。在前27 s内，节理试件接触力主要是从下部试件左侧到上部试件右侧，分布范围较广，但其最大值仅为25 kN左右；27 s 以后，节理试件接触力开始集中于节理表面2个粗糙体周围，且试件内部的接触力有所减小，此时的最大接触力超过45 kN。分析节理面接触力分布与节理表面粗糙体的关系可知，接触力主要分布在节理表面坡角[30-31]较大的凸起上，随着节理面凸起的不断破坏，接触力开始集中于节理面上一些抗剪能力较强的凸起上。此外还可知，在半个周期内节理试件接触力分布变化较大，在整数秒即剪切应力最小时，节理试件接触力以竖直方向为主，仅在小范围内有少许切向方向的接触力，随着剪切应力增大到峰值，在半数秒即剪切应力最大时时，节理切向接触力分布的范围和大小都显著增加。
图10 剪切过程中颗粒间的受力分布
Fig.10 Distribution of forces between particles during shearing
结合节理裂纹分布和力链分布演化规律可得到节理的失稳机理：对节理施加剪切应力后，节理表面坡角较大的凸起开始发挥抗剪作用，这些凸起有大有小，其中较小的二阶凸起因抗剪强度较弱先被破坏；随着这些凸起被破坏，节理表面的接触发生改变，接触力也随即发生改变，逐渐向较大的一阶凸起上集中；当接触力超过一阶凸起的极限强度后，节理整体的抗剪强度下降，发生失稳破坏。
2.4 裂纹发展与剪切应力的关系
图11为剪切过程中节理剪切位移、累计裂纹数量和裂纹率。由图11可知，节理剪切位移的增长趋势与累计裂纹数量的增长趋势总体一致，这表明节理的损伤程度直接体现在剪切位移的变化上。在前25个循环内，节理剪切位移和累计裂纹数量都随着循环剪切应力的施加呈线性增长；在25～30循环内，曲线增长速率明显增大；在30～36循环内，曲线增长速率又开始减小；在36～39循环内，曲线增速迅速增大；最后在39循环以后，曲线增长速率随着循环次数增加又有所减小。这是由于节理上不断产生裂纹导致节理损伤程度加重所致，具体如下：在前25个循环，裂纹率分布较稀松且基本为1，损伤较小，剪切位移增长较小；在25～30循环内，裂纹率有所增大且大都为1，但其分布开始变的密集，该段内大量二阶凸起被破坏，节理损伤较大，故剪切位移变化较大；在30～36循环内，裂纹率分布密度较大，但裂纹率较小，基本为1,极个别为2，该阶段为一阶凸起开始发生轻微破坏，节理损伤进一步累计，剪切位移进一步增大，但增大速率有所减缓；在36～39循环内，裂纹率密集且数值逐渐变大，此时一阶凸起破坏程度加深，节理发生失稳破坏，剪切位移以极大的增速增大，是节理失稳的前兆；最后在39～42循环内，裂纹率密集但最大数值逐渐下降，此时节理已发生失稳破坏，节理的抗剪强度显著下降，每个循环内节理发生较大的剪切变形，但每个循环内产生的裂纹数量逐渐减小。
图11 剪切过程中颗粒间的受力分布
Fig.11 Distribution of forces between particles during shearing
为进一步分析不同阶段内剪切应力与产生裂纹的关系，分别统计了0～5、15～20、25～30、39～44 s的裂纹率和循环剪切应力，如图12所示。在0～5 s内，裂纹率大都为1，且裂纹大都在剪切应力增大到接近峰值时产生，即在该时间段内，主要是在剪切应力加载到峰值附近时产生损伤；在15～20 s内，裂纹率都为1，此时裂纹仍主要是在剪切应力接近峰值时产生，在其他部分，也有少量裂纹产生；在25～30 s内，可明显看出裂纹出现的频率以及数值均大于前2个时间段，此时不仅在剪切应力增大阶段产生裂纹，在剪切应力卸载时也会产生大量裂纹，即表明该时间段内损伤加剧，不仅在剪切应力加载阶段，剪切应力卸载时也同样产生损伤；在35～40 s内，最大裂纹率增大到5，主要分布在剪切应力峰值处附近，其较大的裂纹率也都主要出现在剪切应力大于1.25 MPa部分，这是由于节理残余剪切应力大致为1.25 MPa，在循环剪切应力超过该值后，节理发生较大的剪切错动，因此产生较多裂纹。
图12 剪切过程中裂纹率以及累计裂纹数量演化规律
Fig.12 Evolution of crack rate and cumulative number of cracks during shearing
3 循环特征参数的影响
3.1 加载幅值的影响
为探究加载幅值对节理剪切特性的影响，分别监测了加载幅值为1.5、1.7、1.9、2.1和2.3 MPa时的剪切位移、剪切应力和裂纹数量，得到不同振幅条件下的剪切位移-剪切应力曲线、循环次数-剪切位移曲线、时间-累计裂纹数曲线，如图13所示。
由图13(a)可知，随着加载幅值的增大，曲线的密集度降低，且剪切位移-剪切应力曲线的4个阶段越来越不明显，在加载幅值为1.5 MPa时，仅有前2个阶段，在加载幅值为2.3 MPa时，第2和第3阶段几乎没有。这是由于单个循环内，加载幅值越小，节理越不易产生损伤破坏，剪切位移小，进而造成加载幅值越小，曲线前期越密集。
由图13(b)可知，在加载幅值为1.5 MPa时，经上百个循环后，节理剪切位移仍未达到目标值，仅为0.13 mm，这是由于该幅值过小，难以对节理造成较大的损伤，节理仅产生较小的变形，即表明在有限的循环次数内，扰动应力幅值存在一个扰动阈值，超过该值后才能诱使节理发生破坏，当然该阈值会随着节理损伤的累积而下降，即便是在扰动应力较小值时也能使节理产生损伤，如图12所示。此外可见，随着加载幅值增大，达到相同剪切位移需要的循环次数减小。且随着加载幅值的增大，循环次数-累计剪切位移曲线的稳定阶段更短，甚至消失，如幅值为2.2 MPa时，在初始阶段后，几乎直接进入加速阶段。图12右上角为统计的不同加载幅值条件下，达到目标剪切位移时扰动载荷的循环次数，在振幅为1.7、1.8、2.1和2.3 MPa条件下，节理试件分别经历了41.7、22.3、14.5和8.7个循环后达到目标位移，循环次数随加载幅值增大非线性减小。
图13 不同振幅条件下节理剪切应力、循环加载次数、裂纹发展规律
Fig.13 Graph of the number of cyclic loads,crack development patterns and energy changes under different amplitude conditions
由图13(c)可知，从剪切开始至剪切结束，加载幅值越大，相同时间内产生的裂纹越多。在振幅为1.5、1.7、1.9、2.1和2.3 MPa条件下，最终累计裂纹数量分别为18、1 793、1 513、1 759和1 770。随着振幅的增加累计裂纹数量增长趋势不同，在振幅为1.5 MPa时，曲线几乎水平，不再增长；在振幅为1.7时，曲线呈5段增长，具体见第2.3节分析；在振幅为1.9 MPa时，曲线大致呈4段增长；在振幅为2.1和2.3 MPa时，曲线基本呈3段增长。这是由于振幅的增加，影响了节理表面粗糙的破坏，在振幅较小时，节理表面的小粗糙体先被破坏，在经历足够的循环后较大的粗糙体才会被破坏，而振幅较大时，即作用在节理面上的力较大时，节理表面上的小粗糙体和大粗糙体同时被破坏，单位时间内即会产生较多的裂纹，因此节理表面粗糙体的渐进破坏规律不明显。
3.2 加载频率的影响
为探究加载频率对节理剪切特性的影响，分别监测了加载频率为0.1、0.5、1.0、1.5和2.0 Hz时的剪切位移、剪切应力和裂纹数量，得到不同频率条件下的剪切位移-剪切应力曲线、循环次数-剪切位移曲线、时间-累计裂纹数曲线，如图14所示。
由图14(a)可知，在5个不同频率条件下，节理均达到目标剪切位移，但曲线的形态相差较大。由第1阶段曲线可知，剪切位移相同时，加载频率越低剪切应力越小，由第4阶段曲线可知，加载频率越低，每个循环内的剪切变形越大。以上规律均表明，加载频率越低，一个循环内节理位移越大，即节理产成的损伤越多。
由图14(b)可知，随着加载频率由0.1增加到2.0 Hz，达到目标剪切位移所经历的循环次数分别为26.7、30.5、41.7、49.4、73.6，呈线性增加。频率越高稳定阶段和加速阶段越长，在频率较低，如0.1 Hz时，加速阶段很短暂，直接由稳定阶段到失稳阶段，即表明节理发生失稳破坏更突然。这是由于频率越低周期越长，一个循环内剪切应力作用的时间越长，节理积累的损伤也就越多，进而加速了节理的失稳破坏。
由图14(c)可知，加载频率不同，节理第1次产生裂纹的时间不同，频率越高时间越短。由第3.1节分析可知，循环剪切应力加载幅值超过1.5 MPa后才能对节理造成明显损伤，因此，在剪切初期，仅在扰动应力大于1.5 MPa部分才能造成节理损伤。由于频率高时剪切应力加载到阈值时间短，第1次产生裂纹的时间短。此外可知，加载频率越高累计裂纹曲线越光滑，频率为0.1和0.5 Hz时的曲线有明显的阶梯，这主要是由于频率小时，剪切应力增长速率慢造成的，当剪切应力处于阈值下，裂纹几乎不增多，当剪切应力超过阈值后，裂纹数量迅速增大，因此具有明显的阶梯，而频率大时，剪切应力变化较快，每个循环内损伤较小，产生的裂纹较少，固没有明显的阶梯状曲线。由图14还可知，频率不同最终裂纹数量也不同，频率由0.1～2.0 Hz，最终裂纹数量分别为1 446、1 700、1 793、1 789、1 551。
图14 不同加载频率条件下节理剪切应力、循环加载次数、裂纹发展规律
Fig.14 Graph of the number of cyclic loads,crack development patterns and energy changes under different loading frequencies
4 结 论
(1) 随着剪切载荷的循环加卸，节理剪切应力-剪切位移曲线出现滞回效应，形成多个滞回环，这是由于节理表面损伤所致。节理循环剪切过程可分为4个阶段，其中剪切初期滞回环较小且分布密集，剪切后期滞回环较大且分布疏松。
(2)剪切过程中节理的裂纹分布和接触力分布与节理表面粗糙体相关，在剪切初期裂纹和高接触力大都集中在二阶凸起上，随着二阶凸起破坏，剪切力逐渐作用到较大的一阶凸起上，在一阶凸起破坏后，节理发生失稳破坏。
(3)节理剪切位移的增长趋势与累计裂纹数量的增长趋势总体一致，主要包括缓慢增长和加速增长阶段。剪切初期仅有少量裂纹在应力峰值附近产生，而后期裂纹产生时对应的应力范围和裂纹数量都显著增加。
(4)加载幅值越大，作用在节理上的剪切力越大，每个循环内产生的节理损伤越多，节理剪切位移也越大，达到相同目标位移所需的循环次数与加载幅值呈反比。加载幅值存在有阈值，该阈值会随着节理损伤的累积有所降低。
(5)加载频率越大，一个循环内剪切应力作用的时间越短，每个循环内节理产生的损伤越少，节理剪切位移也越小，达到目标剪切位移所经历的循环次数与加载频率呈正比。加载频率越低，加速阶段越短，节理发生失稳破坏越突然。
[1] 蒋宇静,张孙豪,栾恒杰,等.剪切载荷作用下岩体结构面动态接触特征数值模拟[J].煤炭学报,2022,47(1):233-245.
JIANG Yujing,ZHANG Sunhao,LUAN Hengjie,et al.Numerical simulation of dynamic contact characteristics of rock discontinuities under shear loading[J].Journal of China Coal Society,2022,47(1):233-245.
[2] YIN Q,WU J J,ZHU C,et al.Shear mechanical responses of sandstone exposed to high temperature under constant normal stiffness boundary conditions[J].Geomechanics and Geophysics for Geo-Energy and Geo-Resources,2021,7(2):1-17.
[3] XIA Caichu,TANG Zhicheng,XIAO Weimin,et al.New peak shear strength criterion of rock joints based on quantified surface description[J].Rock Mechanics and Rock Engineering,2014,47(2):387-400.
[4] ASADI M S ,RASOULI V ,BAELA G.A laboratory shear cell used for simulation of shear strength and asperity degradation of rough rock fractures[J].Rock Mechanics and Rock Engineering,2013,46(4):683-699.
[5] ZHAO Zhihong,PENG Huan,WU Wei,et al.Characteristics of sh-ear-induced asperity degradation of rock fractures and implications for solute retardation[J].International Journal of Rock Mechanics and Mining Sciences,2018,105:53-61.
[6] 尹乾,靖洪文,孟波,等.CNL 和 CNS 边界条件下砂岩宏细观剪切力学特性[J].采矿与安全工程学报,2021,38(3):615-624.
YIN Qian,JING Hongwen,MENG Bo,et al.Macro meso shear mechanical properties of sandstone under CNL and CNS boundary conditions[J].Journal of Mining and Safety Engineering,2021,38(3):615-624.
[7] MENG Fanzhen,ZHOU Hui,LI Shaojun,et al.Shear behaviour and acoustic emission characteristics of different joints under various stress levels[J].Rock Mechanics and Rock Engineering,2016,49(12):4919-4928.
[8] 周友华,胡奉湘,饶云阶,等.地震应力产生传递与集中机制的研究[J].国际地震动态,2004(11):8-14.
ZHOU Youhua,HU Fengxiang,RAO Yunjie,et al.Researches on the mechanisms for the transform and concentration of seismic stresses[J].Recent Developments in Woeld Seismology,2004(11):8-14.
[9] HAN Guansheng,JING Hongwen,JIANG Yujing,et al.Effect of cyclic
loading on the shear behaviours of both unfilled and infilled rough rock joints under constant normal stiffness conditions[J].Rock Mechanics and Rock Engineering,2020,53(1):31-57.
[10] LI Yingchun,WU Wei,TANG Chun’an,et al.Predicting the shear
characteristics of rock joints with asperity degradation and debris backfilling under cyclic loading conditions[J].International Journal of Rock Mechanics and Mining Sciences,2019,120:108-118.
[11] LEE H S,PARK Y J,CHO T F,et al.Influence of asperity degradation on the mechanical behavior of rough rock joints under cyclic shear loading[J].International Journal of Rock Mechanics and Mining Sciences (Oxford,England:1997),2001,38(7):967-980.
[12] 刘博,李海波,刘亚群.循环剪切荷载作用下岩石节理变形特性试验研究[J].岩土力学,2013，34(9):2475-2481,2488.
LIU Bo,LI Haibo,LIU Yaqun.Experimental study of deformation behavior of rock joints under cyclic shear loading[J].Rock and Soil Mechanics,2013，34(9):2475-2481,2488.
[13] NIKTABAR S M Mahdi,RAO K Seshagiri,SHRIVASTAVA Amit Kumar.Automatic static and cyclic shear testing machine under constant normal stiffness boundary conditions[J].Geotechnical Testing Journal,2018,41(3):20170083.
[14] MIRZAGHORBANALI Ali,NEMCIK Jan,AZIZ Naj.Effects of sh-ear rate on cyclic loading shear behaviour of rock joints under constant normal stiffness conditions[J].Rock Mechanics and Rock Engineering,2014,47(5):1931-1938.
[15] MIRZAGHORBANALI Ali,NEMCIK Jan,AZIZ Naj.Effects of cyc-lic loading on the shear behaviour of infilled rock joints under constant normal stiffness conditions[J].Rock Mechanics and Rock Engineering,2014,47(4):1373-1391.
[16] DANG Wengang,KONIETZKY Heinz,CHANG Lifu,et al.Velocity-frequency-amplitude-dependent frictional resistance of planar joints under dynamic normal load (DNL) conditions[J].Tunnelling and Underground Space Technology,2018,79:27-34.
[17] DANG Wengang,KONIETZKY Heinz,FRÜHWIRT Thomas,et al.Cyclic frictional responses of planar joints under cyclic normal load conditions:Laboratory tests and numerical simulations[J].Rock Mechanics and Rock Engineering,2020,53(1):337-364.
[18] 杨小彬,周杰,宋义敏,等.循环加载岩石界面滑移位移演化特征试验研究[J].煤炭学报,2019,44(10):3041-3048.
YANG Xiaobin,ZHOU Jie,SONG Yimin,et al.Evolution characteristics of sliding displacement ofrock interface under cyclic loading[J].Journal of China Coal Society,2019,44(10):3041-3048.
[19] FATHI Ali,MORADIAN Zabihallah,RIVARD Patrice,et al.Sh-ear mechanism of rock joints under pre-peak cyclic loading condition[J].International Journal of Rock Mechanics and Mining Sciences,2016,83:197-210.
[20] 刘新荣,邓志云,刘永权,等.岩石节理峰前循环直剪试验颗粒流宏细观分析[J].煤炭学报,2019,44(7):2103-2115.
LIU Xinrong,DENG Zhiyun,LIU Yongquan,et al.Macroscopic and microscopic analysis of particle flow in pre-peak cyclic direct shear test of rock joint[J].Journal of China Coal Society,2019,44(7):2103-2115.
[21] 刘新荣,邓志云,刘永权,等.峰前循环剪切作用下岩石节理损伤特征与剪切特性试验研究[J].岩石力学与工程学报,2018,37(12):2664-2675.
LIU Xinrong,DENG Zhiyun,LIU Yongquan,et al.Experimental study on damage characteristics and shear characteristics of rock joints under cyclic shear before peak[J].Journal of Rock Mechanics and Geotechnical Engineering,2018,37(12):2664-2675.
[22] KOU Miaomiao,LIU Xinrong,TANG Shangding,et al.Experimental study of the prepeak cyclic shear mechanical behaviors of artificial rock joints with multiscale asperities[J].Soil Dynamics and Earthquake Engineering,2019,120:58-74.
[23] LIU X R ,KOU M M ,LU Y M,et al.An experimental investigation on the shear mechanism of fatigue damage in rock joints under pre-peak cyclic loading condition[J].International Journal of Fatigue,2018,106:175-184.
[24] 蒋宇静,张孙豪,栾恒杰,等.恒定法向刚度边界条件下锚固节理岩体剪切特性试验研究[J].岩石力学与工程学报,2021,40(4):663-675.
JIANG Yujing,ZHANG Sunhao,LUAN Hengjie,et al.Experimental study on shear characteristics of anchored jointed rock mass under constant normal stiffness boundary conditions[J].Journal of Rock Mechanics and Geotechnical Engineering,2021,40(4):663-675.
[25] 尹乾,靖洪文,孟波,等.恒定法向刚度条件下三维粗糙裂隙面剪切力学特性[J].岩石力学与工程学报,2020,39(11):2213-2225.
YIN Qian,JING Hongwen,MENG Bo,et al.Shear mechanical properties of 3D roughrock fracture surfaces under constant normal stiffness conditions[J].Chinese Journal of Rock Mechanics and Engineering,2020,39(11):2213-2225.
[26] JIANG Y,XIAO J,TANABASHI Y,et al.Development of an automated servo-controlled direct shear apparatus applying a constant normal stiffness condition[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(2):275-286.
[27] BAHAADDINI M,HAGAN P C,MITRA R,et al.Experimental and
numerical study of asperity degradation in the direct shear test[J].Engineering Geology,2016,204:41-52.
[28] Itasca Consulting Group Inc (ICG).PFC2D Version 4.0:Theory and
Background[Z].Itasca Consulting Group Inc,Minneapolis,2008.
[29] BAI Qingsheng,KONIETZKY Heinz,DANG Wengang.Microscopic modeling of frictional response of smooth joint under normal cyclic loading[J].Rock Mechanics and Rock Engineering,2022,55:169-186.
[30] GRASSELLI.Shear Strength of rock joints based on quantified surface description[D].Switzerland:Ecole Polytechnique Federale de
Lausanne,2001.
[31] DANG Wengang,WU Wei,KONIETZKY Heinz,et al.Effect of she-ar-induced aperture evolution on fluid flow in rock fractures[J].Computers and Geotechnics,2019,114:103152.
