花岗岩真三轴加、卸载破坏的细观损伤机制研究

刘婕，丛宇,2，张黎明,2，王在泉,2

(1.青岛理工大学理学院，山东青岛，266033；2.青岛理工大学蓝色经济区工程建设与安全协同创新中心，山东青岛，266033)

随着地下工程埋深的增加，岩体所赋存的地质环境更为复杂。在深部地下工程开挖过程中，围岩体的初始三向受力状态发生改变，1个或者几个方向的应力被卸载，经常会导致岩爆、塌方、冒顶等灾害[1−3]。这些灾害不仅影响工程进度，而且威胁施工人员的生命安全，因此，研究真三轴应力状态下岩体破坏机制具有重要意义。

早期岩石加、卸载试验以圆柱试样的常规三轴试验为主，首先将岩样加载到一定应力，再以一定的速率卸载，研究常规三轴卸载路径下岩样石力学特性。DAI 等[4−5]通过研究发现，常规三轴试验只能对圆柱试样施加相同的侧向应力，卸载时2个方向的侧向应力同时被卸载，这与实际开挖过程中巷道围岩应力状态的变化不符。为此，研究者开展了一系列室内真三轴试验。何满潮等[6−7]采用真三轴试验模拟岩爆发生，得到了花岗岩岩爆发生的过程；SU等[8−9]通过真三轴卸载试验分析了岩爆弹射破坏过程的特征与规律，并通过开展不同切向应力加载速率下应变型岩爆试验研究发现，随着加载速率增大，岩样的破坏模式逐渐由无弹射的脆性破坏(劈裂破坏)转化为动力岩爆弹射破坏；许文松等[10]采用真三轴卸载扰动岩石测试系统对大理岩进行不同主应力加卸载试验，得到了黏聚力和内摩擦角的变化规律；LI等[11]通过真三轴卸载试验发现，试样高宽比和中间主应力对硬岩的破坏模式、峰值强度和严重程度均有影响；LI等[12]通过真三轴压缩试验，研究了砂岩在不同压力下的渐进损伤变形特征和体积应变变化规律，建立了含微裂纹缺陷的三维岩石理论模型；SI等[13−14]对立方细粒花岗岩试件进行了三轴卸载压缩试验，发现卸载对细粒花岗岩具有明显的强度减弱效应。

由于试验条件受限，室内试验难以观察试样内部裂纹发展过程，研究者利用三维颗粒流软件PFC3D 进行室内试验的数值模拟。PFC3D 是基于细观离散元理论的相关程序，可以更好地再现岩石类材料的基本特性，特别是颗粒物质的动力响应、岩石类介质的破裂及裂纹发展[15−17]。SUN等[18]基于人工神经网络等方法对PFC3D 中岩石三轴压缩模拟试验的微观参数进行了预测，并获得了可以准确反映硬岩性质的微观参数；马春驰等[19]利用PFC3D 进行常规三轴卸载颗粒流数值模拟，对不同围压下的岩爆效应进行分析发现，低围压卸载与高围压卸载下岩样的破坏均以张拉破坏为主导；丛怡等[20]进行了常规三轴加、卸试验的数值模拟，提出了通过围压−轴向应变关系来获取数值模拟试验卸载速率与室内试验卸载速率之间的对应关系，解决了PFC3D 中卸载速率的定量控制问题；张杨[21]通过PFC3D 模拟真三轴加载试验条件下岩石的损伤破坏过程发现，中主应力对岩石破坏起先保护后加剧损伤的作用，岩石的破坏主要是黏结拉破坏。吴顺川等[22]进行了真三轴卸载试验的模拟研究，发现卸载岩爆试验进程可分为平静期、局部颗粒弹射期、发展期及最终爆发期4种状态。

综上可知，虽然研究者在室内真三轴加、卸载试验方面进行了大量研究，但对岩石真三轴加、卸载破坏过程中能量演化特征研究尚不充分。三维颗粒流模拟刚起步，主要通过PFC3D 对岩石试验的破坏过程进行了模拟研究，对破坏过程中细观能量的演化以及裂纹扩展方面的研究较少。为此，本文在花岗岩真三轴加、卸载试验的基础上，利用颗粒流离散元程序PFC3D 探索真三轴加、卸载应力路径下花岗岩的破坏机制。

1 花岗岩真三轴破坏试验

1.1 试样制备

花岗岩岩样取自某地下洞库，中细粒结构，块状构造，主要矿物成分为石英、钾长石、斜长石、黑云母。选取同批次花岗岩开展室内力学试验研究，对现场收集的新鲜岩块进行钻芯、打磨，按照GBT 50266—1999“工程岩体试验方法标准”制成长×宽×高为100 mm×50 mm×50 mm 的标准试样。

1.2 真三轴试验

试验在中科院武汉岩土力学研究所研制的真三轴测试系统上进行。该测试系统由全刚性力学加载系统和伺服控制与数据采集系统2个子系统构成，可真实模拟岩石在深部地层三维应力下的环境。试验过程中，在垂直方向施加最大主应力σ1，水平方向分别施加中间主应力σ2和最小主应力σ3。

真三轴加载试验过程如下：同步增加σ2和σ3到设定围压(σ2=40 MPa，σ3=30 MPa)；保持σ2和σ3不变，以0.1 mm/min 的速度加载竖向荷载σ1至岩样破坏，利用数据采集系统记录岩样变形全过程的应力、应变等参数。

真三轴卸载试验过程如下：同步增加σ2和σ3到其目标值(σ2=40 MPa，σ3=30 MPa)；以0.1 mm/min速度加载σ1至真三轴加载试验峰值强度σc的80%；保持3 个方向应力不变2 min 后，通过位移控制以不同的速率(0.3，0.45和0.55 mm/min)单面卸载σ3，利用数据采集系统记录岩样变形全过程的应力、应变等，卸载应力路径如图1所示。

图1 卸载路径示意图Fig.1 Schematic diagram of unloading path

2 花岗岩真三轴破坏试验结果分析

2.1 花岗岩真三轴试验应力−应变曲线分析

真三轴加、卸载试验应力(σ1)−应变(ε)曲线分别如图2和图3所示。根据裂纹体积应变曲线和总体积应变曲线[23]，加载试验应力−应变曲线可以分为初始压密阶段OA、弹性变形阶段AB、裂纹稳定扩展阶段BC、裂纹不稳定扩展阶段CD和峰后破坏阶段DE。卸载试验应力−应变曲线则可以分为初始压密阶段OA、弹性变形阶段AB、裂纹稳定扩展阶段BC和卸荷破坏阶段CD。

图2 真三轴加载试验应力−应变曲线Fig.2 Stress−strain curves of true triaxial loading test

图3 真三轴卸载试验应力−应变曲线Fig.3 Stress−strain curve of true triaxial unloading test

2.2 花岗岩真三轴试验宏观破坏特征

图4所示为花岗岩真三轴加、卸载试验岩样的破坏形式，其中蓝色表示张拉裂纹，红色表示剪切裂纹。花岗岩试件的破坏模式为张剪复合破坏，岩样表面分布有多条剪切裂纹和张拉裂纹。岩样侧向扩容主要沿卸荷方向进行，宏观裂纹产生并主要沿着方向σ2扩展，最后贯通形成了贯穿整个岩样的主破裂面。对比加、卸载应力路径下岩样的裂纹分布发现，在2种应力路径下，岩样在沿σ2方向的裂纹分布位置存在差异；在卸载路径下，岩样在沿σ2方向的面上裂纹的分布近似呈“V”型，裂纹贯通形成了2个主破裂面，而加载路径下裂纹贯通仅形成了1个主破裂面。

图4 花岗岩真三轴加载与卸载试验宏观破坏特征Fig.4 Macroscopic failure characteristics of granite in true triaxial loading and unloading tests

2.3 能量演化特征分析

依据能量守恒理论，忽略试验过程中的能量损失，岩样吸收的总能量等于外力所做的功。在真三轴试验中，外力做的功即为试验机对岩样施加的轴压与围压，分别导致岩样产生轴向与侧向变形。因此，真三轴应力状态下岩样实际吸收的总能量Utotal为[24-27]：

式中：σ1,i和ε1,i分别为应力−应变曲线计算点i所对应的最大主应力和最大主应变；σ2,i和ε2,i分别为计算点i的中间主应力和中间主应变；σ3,i和ε3,i分别为计算点i的最小主应力和最小主应变。

按照岩石变形过程中能量性质的不同划分，可将岩样吸收的总能量分为弹性应变能Ue与耗散能Ud，即

真三轴试验中弹性应变能的计算公式为

式中：E为岩石的初始弹性模量；μ为岩石的泊松比。为方便计算，E和μ取真三轴加载试验中弹性变形段进行计算。

真三轴加载试验过程中能量−应力−应变曲线如图5(a)所示。可见：在初始压密阶段(OA段)，岩样吸收的能量较少，外力功主要以弹性应变能形式储存，总能量曲线与弹性应变能曲线基本重合；在线弹性阶段(AB段)；在该阶段前期，总能量、弹性应变能曲线稳定增长，耗散能基本无变化，后期由于岩样内部会产生少量微裂纹，耗散能继续缓慢增加，外力持续做功，总能量与弹性应变能继续增加，但是二者增加速率不同，并逐渐分离；在裂纹稳定扩展阶段(BC段)，岩样内部微裂纹持续萌生、扩展消耗的能量增多，耗散能稳定增加，弹性能增速变缓；在裂纹非稳定扩展阶段(CD段)，岩样内部微裂纹快速扩展、贯通，需要消耗大量的能量，耗散能增长速率迅速增大，弹性应变能增长速率明显减小，表明该阶段耗散能在能量分配中所占的比例进一步增大；在峰后破坏阶段(DE段)，宏观裂纹贯通形成主破裂面，导致岩样内部积聚的一部分弹性应变能以耗散能的形式释放，耗散能以较大速率持续增长，其量值甚至超过弹性能，表明岩样在宏观破裂面形成之前需要消耗大量的能量[25−26]。

真三轴卸载试验的能量过程曲线如图5(b)所示。可见，各能量−应变曲线在压密阶段、线弹性阶段、裂纹稳定扩展阶段的变化规律与加载试验的基本一致。卸载开始后，卸荷点C点处，轴向应力−应变曲线大幅度下降，耗散能大幅度增加且增速加快，弹性应变能快速释放，裂纹迅速扩展；随着最小主应力的减小，宏观裂纹贯通，岩样出现宏观破坏并沿卸荷方向迅速扩容，此时，耗散能曲线大幅度增长，弹性应变能进一步释放，岩样完全破坏。

对比图5(a)和(b)可知：在真三轴加、卸载2种应力路径下，花岗岩在压密阶段、线弹性阶段、裂纹稳定扩展阶段的宏观能量演化规律基本一致，在试验后期即破坏阶段的宏观能量演化规律有所差别；与加载应力路径相比，卸载应力路径下岩样破坏时各能量的变化更剧烈。在真三轴卸载试验中，轴向应力在卸载瞬间发生跌落，弹性应变能曲线迅速下降，且相较于加载路径曲线下降的速率更快、变化更陡；在卸载时，耗散能曲线迅速增加，增值为1.187 MJ/m3，而在加载路径下，耗散能在岩样破坏时的突增值为0.379 MJ/m3，表明卸载路径下岩样发生破坏的过程迅速。此外，卸载路径下岩样发生破坏时的耗散能量为2.089 MJ/m3，明显小于加载路径的耗散能3.111 MJ/m3，表明卸载路径下岩样发生破坏时消耗的能量更少，相较于加载路径更危险。

图5 真三轴加卸载试验能量曲线Fig.5 Energy curve of true triaxial loading and unloading test

3 花岗岩真三轴试验数值仿真模型

3.1 颗粒流仿真计算

数值模拟试验通过颗粒流离散元软件PFC3D完成。模型尺寸与室内试样尺寸一致。通过“墙”模拟模型边界，颗粒直径为2.50～3.75 mm，在“墙”组成的空间内生成规定数目的颗粒，颗粒半径在最大与最小半径范围内随机分布。利用伺服调节法使模型内部应力达到平衡状态，颗粒分布均匀并满足设定的半径以及孔隙率等，得到仿真计算模型。模型的建立过程如图6所示。

图6 数值模型建立过程Fig.6 Numerical model building processes

3.2 颗粒流仿真参数

颗粒间的接触设置为平行黏结模型。平行黏结是为颗粒间提供一定截面形状和厚度的黏结材料，将颗粒黏结在一起，能够同时传递力和力矩，可以更好地反映岩石材料的塑性特征[21]。由于PFC中并没有给出宏细观参数之间明确的对应关系，为保证数值模拟结果能逼近实际岩石的宏观力学响应，先依据宏观临界强度初步选定细观参数，然后通过试错法进行细观参数标定，即调整各个细观参数以确保数值模型所得应力−应变关系、峰值强度、峰值应变室内试验相吻合。当数值模拟与室内试验得到的应力−应变曲线基本一致时，认为该组细观参数为满足计算所用细观参数。

根据丛宇等[28]得到的岩石类材料宏观力学特征与细观参数之间的定量关系，主要通过调整颗粒−颗粒接触模量Ec使数值模型的宏观弹性模量与室内试验相近，调整与得到合适的峰值轴向应力，调整与得到与室内试验一致的破坏形式以及裂纹分布。通过试错法最终标定的细观参数如表1所示。根据表1，模拟真三轴加载试验的应力−应变曲线及破坏形式与室内试验的对比如图7所示，数值模型中红色部分表示剪切裂纹，蓝色部分表示张拉裂纹。从图7可见：数值模拟的应力−应变曲线与试验曲线基本重合，破坏方式、裂纹的分布也基本一致；在初始加载阶段，由于实际岩样微裂隙的压密，室内试验曲线为下凹曲线，而数值模拟在颗粒生成过程中，球体已在其重力加速度作用下进行压密，因此，数值试验没有初始压密阶段。以上结果表明：该组细观参数能准确描述花岗岩真三轴破坏的主要力学特征。

表1 PFC数值模型细观参数Table 1 PFC numerical model meso-structure parameters

图7 室内试验与数值模拟应力−应变曲线与破坏形式对比Fig.7 Comparison of stress−strain curve and failure form between indoor test and numerical simulation

3.3 颗粒流仿真实现

颗粒流数值模拟的应力路径与室内试验的一致。PFC内嵌Fish语言，可通过伺服调节法控制模型“墙体”运动，进而实现加载与卸载的过程。与室内试验以mm/min 为卸载速率的单位不同，PFC3D程序中利用时间步来记录加、卸载的过程，卸载速率单位为mm/步。参考丛怡等[20]建立二者对应关系的方法，可以实现室内试验PFC3D 模拟试验之间的单位转换，即认为相同应力路径下，若数值模拟与室内试验最小主应力−轴向应变曲线斜率一致，则认为两者卸荷速率相同，得到的对应关系为室内试验卸载速率0.30 mm/min 相当于PFC3D程序卸载速率0.01 mm/步。图8所示为卸荷速率为0.3 mm/min 的室内试验与模拟试验中的围压−轴向应变曲线。

图8 室内与模拟试验中的最小主应力−轴向应变曲线Fig.8 Minimum principal stress−axial strain curve in laboratory and simulation test

真三轴加载试验的模拟过程主要可以分为2个过程。首先，利用PFC3D 中的伺服控制给四周“墙体”施加围压，使σ2和σ3达到设定值；然后保持σ2和σ3恒定，通过伺服控制使顶面和底面的“墙”以恒定的速度相向运动，实现模拟室内试验中轴向加载过程，加载至岩样破坏。

真三轴卸载试验的数值模型是在加载试验模型的基础上进行调节得到的，主要分以下3 个过程：生成六面“墙”，利用伺服控制给四周“墙体”施加围压，使σ2和σ3达到设定值；保持σ2和σ3不变，利用伺服控制使顶面及底面的“墙”以恒定的速度相向运动模拟轴向加载；加载轴压至真三轴压缩试验峰值应力80%后，利用伺服控制保持轴压σ1恒定，移动σ3方向的一面“墙”以实现单面卸载最小主应力，直至试样发生破坏。

4 数值模拟结果分析

4.1 破坏细观破坏特征

图9所示为数值模拟试验中岩样的破坏形式，其中红色表示剪切裂纹，蓝色表示张拉裂纹。从图9可见：颗粒流模拟试验中岩样的破坏模式也为张剪复合破坏，张拉裂纹与剪切裂纹主要从岩样两端向中间扩展，形成贯通的宏观主破裂面，在与宏观主破裂面近似垂直的方向形成次要破坏面。对比图9与图4可知：数值模拟试验与室内试验的岩样沿σ2方向上裂纹的分布基本一致，沿σ3方向的裂纹分布有所不同。由于室内试验的岩样内部存在原始裂隙，岩样破坏时在沿σ3方向产生了有少量的剪切裂纹，而数值模拟试验的模型中没有原始裂隙，因此，岩样模型在沿σ3方向无明显贯通的裂纹。

图9 数值模拟试样破坏形式Fig.9 Numerical simulation of specimen failure form

4.2 细观能量分析

PFC3D 中的能量分析是从细观角度考虑岩样内部能量变化，细观能量主要包括边界能、黏结能、摩擦能、动能和应变能。其中，边界能是由与模型边界接触的球能量累计而成的，是边界作用力与位移的乘积；黏结能是储存在模型中所有颗粒间平行黏结模型储存的应变能，是克服颗粒间黏结力做的功；摩擦能是颗粒间滑动消耗的能量，是摩擦力与位移的乘积；动能是系统颗粒运动消耗的能量；应变能是颗粒与颗粒间累计储存的能量[14]。

图10所示为加、卸载应力路径下数值模拟的细观能量演化规律。从图10可见：在加载前期，模型处于被压缩的状态，“墙体”边界一直对模型做正功，边界能量曲线基本呈增长趋势，增长速率逐渐增大，轴向应变达到一定值(加载时为0.050 1，卸载时为0.037 1)，模型开始发生明显的侧向变形，增长速率逐渐减缓，整个过程中边界能明显高于其他能，表明试验过程中模型吸收的能量主要来源于边界能。黏结能变化特征与应力−应变曲线形式基本对应，这是因为随着微裂隙的产生、扩展，颗粒间的黏结发生张拉断裂，存储在平行黏结模型中的应变能释放。黏结能仅小于边界能，是整个试验消耗能量的主体。当轴向应变达到一定值时，模型发生破坏，黏结能与应变能出现峰值，由于微裂纹进一步扩展、贯通，颗粒间的黏结断裂并发生剪切滑动，此时，摩擦能的增长速率增大；与其他细观能量相比，动能变化不明显，表明试验过程中颗粒运动不剧烈。

对比图10(a)和图10(b)可知：在轴向应力上升阶段，在加、卸载2种应力路径下，各细观能量的变化趋势一致；在轴向应力下降阶段，2种路径下各细观能量曲线的变化则有所不同；在加载应力路径下，当轴向应力达到峰后97.4%σc时，应变能与黏结能达到最大值，摩擦能曲线的增长速率增大；而在卸载应力路径下，随着σ3的卸载，轴向应力逐渐减小并在最后发生陡降，应变能与黏结能也是在最后破坏时出现陡减，同时，摩擦能曲线近似呈垂直上升趋势，由0.0459 MJ/m3迅速增到0.113 MJ/m3，各能量曲线变化率明显大于加载路径时的变化率，表明在卸载路径下花岗岩模型发生破坏更剧烈。

图10 加、卸载数值模拟试验的细观能量转化Fig.10 Meso-level energy conversion of loading and unloading numerical simulation test

4.3 裂纹扩展过程分析

图11所示为真三轴加载数值模拟试验中轴向应力分别加载到起裂应力、损伤应力、峰值应力以及峰后70%时细观裂纹的分布情况，其中红色表示剪切裂纹，蓝色表示张拉裂纹。从图11可见：当轴向应力加载到起裂应力时，裂纹开始萌发，张拉裂纹与损伤裂纹数量均比较少；加载到损伤应力时，张拉裂纹零散遍布于整个岩样，剪切裂纹开始萌发，零星分布；达到峰值强度时张拉裂纹布满岩样，剪切裂纹分布零散；加载到峰后70%时，张拉裂纹密集布满岩样，剪切裂纹主要由两端向中间发展，贯通形成宏观的剪切破坏面带。

图11 真三轴加载试验破坏的细观损裂纹演化过程Fig.11 Micro-damage crack evolution processes of true triaxial loading test failure

图12所示为真三轴卸载数值模拟试验轴向应力分别加载到起裂应力、峰前80%(即卸荷点)以及σ3卸载到15 MPa，10 MPa和破坏时细观裂纹的分布情况。从图12可见：当轴向应力加载到起裂应力时，裂纹开始萌发，张拉裂纹与损伤裂纹数量均比较少；当轴向应力加载到峰前80%后开始卸载σ3，张拉裂纹零散分布，剪切裂纹数量极少；当σ3卸载到15 MPa和10 MPa时，张拉裂纹增多并零散分布，剪切裂纹增加不明显，零星分布；在破坏时，张拉裂纹与剪切裂纹迅速增多，张拉裂纹与剪切裂纹主要由两端向中间发展，最后贯通形成1条主剪切带以及2条次剪切带。

图12 真三轴卸载试验破坏的细观损裂纹演化过程Fig.12 Micro-damage crack evolution processes of true triaxial unloading test

图13所示为竖向荷载σ1和裂纹数量关系。从图13可见：在加载应力路径下，张拉裂纹与剪切裂纹增长趋势基本一致，初始阶段裂纹数量增加很少，裂纹曲线接近于水平直线，随着荷载增加，裂纹增多，曲线斜率逐渐增大；当达到峰值应力之后，裂纹曲线呈直线上升状态；在卸载应力路径下，试样内部裂纹在卸载后开始逐渐扩展，并随着最小主应力的减小，数量逐渐增多；在岩样破坏时，裂纹数量急剧增加，且剪切裂纹的增长速率大于张拉裂纹的增长速率。

图13 裂纹扩展数量Fig.13 Crack growth number

图14所示为加、卸载应力路径下不同应力水平裂纹数量。从图14可见：在加、卸载路径下，岩样内部张拉裂纹数量均明显多于剪切裂纹数量；在加载应力路径下，随着荷载增加，剪切裂纹数在总裂纹数中的占比由0 增加到28.8%，在卸载应力路径下，随着最小主应力的卸载，剪切裂纹数在总裂纹数中的占比由0.1%增加到21.0%。

图14 不同应力水平裂纹数量Fig.14 Number of cracks at different stress levels

由图13和14可知：真三轴卸载路径下花岗岩模型首先发生的是张拉破坏，且达到失稳破坏阶段后，张拉裂纹数量仍明显多于剪切裂纹数量，张剪裂纹比均大于2，说明颗粒间的黏结破坏以张拉破坏为主。此外，进入破坏阶段时，模型的张剪裂纹比均有所下降，在加载路径下，张剪裂纹比由5.77 降到2.47，在卸载路径下，由38.5 降到3.76，表明模型发生宏观破坏时张拉裂纹与剪切裂纹同时发育，因此，模型表现为张剪复合破坏形式。

5 结论

1)在加、卸载应力路径下，花岗岩宏观能量演化规律的不同主要体现在后期接近破坏阶段，在卸载路径下，弹性应变能曲线迅速下降，相较于加载路径，曲线下降的速率更快、变化更大；耗散能曲线在卸载瞬间迅速增加，增长的幅度与速率明显比加载路径下的大，这表明卸载路径下岩样发生破坏需要的时间更少、更迅速。

2)在不同应力路径下，克服颗粒间黏结力做功的黏结能仅小于边界能，是消耗能量的主体，当轴向应变达到一定值时，黏结能与应变能出现峰值，颗粒间的黏结发生剪切断裂，摩擦能成为消耗能量的主体，岩样发生明显的变形破坏。

3)在加、卸载应力路径下，各细观能量曲线演化特征的不同主要表现在破坏阶段，在卸载路径下，随着σ3减小，应变能和黏结能先逐渐增大并在最后破坏时出现陡降，同时，耗散能曲线近似呈垂直上升趋势，各能量曲线变化率要明显大于加载路径下的变化率，表明在卸载路径下模型发生破坏更剧烈。

4)在相同围压下，岩样内部张拉裂纹数量均多于剪切裂纹数量，且卸载应力路径下的裂纹总数少于加载路径下的裂纹总数，表明卸载路径下岩样达到破坏时的内部损伤较少；在卸载试验中，岩样破坏时，裂纹数量大量增加，表明在卸载应力路径下，损伤累积的过程较短，岩样快速破坏。

5)在加、卸载应力路径下，花岗岩的破坏形式均为张剪复合破坏；花岗岩模型发生破坏时的颗粒黏结破裂机制以张拉破裂为主，剪切破裂为辅。