自然界中岩体处于复杂应力状态，在进行深部洞室开挖过程中，临空面会释放岩体中积累的应变能，造成应力重分布，导致岩体内部产生新裂纹，并且裂纹与裂纹之间扩展贯通，最终会造成岩体的失稳破坏^[1-3]。以往的常规三轴压缩试验忽略了中间主应力(σ₂)的影响，传统的莫尔-库仑和霍克-布朗准则对真三轴应力状态下岩石力学表征难以适用，而中间主应力对岩石在真三轴应力状态下的力学特性影响显著^[4]。近年来，许多学者对深埋大理岩进行真三轴压缩试验，以探究其在高地应力下的破坏机制。Zheng等^[4]探讨了锦屏大理岩在真三轴压缩下的残余强度特征，发现对于较高σ₃(最小主应力)，随着σ₂的增大残余强度先减小后增大，对于较低σ₃，则随着σ₂的增大而逐渐减小；Feng等^[5]对花岗岩、大理岩和砂岩等进行了真三轴循环加卸载试验，指出岩石黏聚力对σ₃和σ₂不敏感，内摩擦角最小值不受σ₃和σ₂的影响，但其最大值受两者的影响；Gao等^[6]对含节理的锦屏大理岩进行真三轴压缩试验，发现节理大理岩强度对σ₃的敏感性要大于对σ₂的敏感性。

然而，在实验室进行真三轴试验需要先进且昂贵的实验设备，试件的制备也费时耗力，并且试验中很难观察到细观的破坏过程。近年来，数值模拟方法已成为传统研究方法的有力补充。其中，离散单元法作为一个有效的数值分析工具，广泛应用于岩石力学研究^[7-8]。Yang等^[9]研究了射孔角度和射孔布置对同步水力压裂裂缝扩展机制和岩石破坏模式的影响，指出大主应力差下较小的射孔间距更有利于中心岩体的破裂。Lee等^[10]对非平行双裂隙花岗岩试样进行单轴压缩试验，并且开展了颗粒流程序模拟，发现岩桥裂纹的贯通主要通过拉伸裂纹的扩展聚集。Zhang等^[11]对含两个平行裂隙的矩形类岩石试件进行单轴压缩数值模拟，总结出岩桥贯通两种新模式。由于实际工程问题复杂，简化为二维问题后，其计算结果存在一定差距，故开展三维数值研究十分必要。Duan等^[12]对砂岩展开真三轴压缩离散元模拟，重点分析σ₃、σ₂对微裂纹取向、黏结力取向、组构和配位数演化的影响，指出随着σ₃水平的增加，σ₂对颗粒尺度响应的影响逐渐减弱。Zhang等^[13]采用平行黏结模型先标定大理岩的力学性能，然后建立真三轴压缩数值模型，在数值岩样上预制裂隙，研究裂隙大理岩在真三轴应力状态下的破坏特征。Zhang等^[14]采用黏结粒子模型研究花岗岩在真三轴应力状态下的变形和强度特征，消除了端部效应，隔离了中间主应力的影响。

综上，中间主应力对岩石力学特性影响较大，且目前对真三轴应力环境下的岩石内部裂隙扩展演化过程难以做到实时观察。为此，拟采用三维离散元颗粒流程序(PFC^3D)对大理岩试件开展真三轴应力状态下的微观裂纹扩展机制数值研究，并基于全应力-应变曲线特征，深入探析裂纹扩展过程与受压变形过程当中的内在联系，同时与实验观测结果进行比较，力求得出具有参考价值的结论。

1 数值模型的建立 1.1 数值模型

PFC5.0内嵌有十种接触模型^[8]，其中，平行黏结模型在模拟岩石类材料受压破坏时，能够沿着法向或者切向破坏，可较好地反映此类材料的破坏形式，已广泛应用于岩石类材料相关问题分析中^[15]，故本文选择该模型进行后继研究。大理岩试样数值模型尺寸为50 mm×50 mm×100 mm，假定无内部缺陷。采用与文献[16]常规三轴压缩数值试验相同的平行黏结模型细观参数，大理岩颗粒的半径在0.9~1.356 mm，避免颗粒尺寸对力学性能的影响。其中，真三轴数值模型包含26 018个颗粒，112 093个平行键黏结，4 504个线性黏结，如图 1所示。

1.2 真三轴压缩模拟加载方式

数值模拟中的加载应力路径是仿照岩石真三轴压缩试验的加载方式^[17-18]，如图 2所示，模拟的加载路径由以下3个阶段组成：

1) 在静水条件下(σ₁=σ₂=σ₃)，通过伺服控制函数控制6个墙体以一定位移速率单调加载试样，直到达到预设的σ₃水平。

2) 通过伺服控制函数控制x方向上两个平行墙体保持σ₃恒定，同时控制y、z方向两对平行墙体以恒定位移速率继续加载，保证σ₁和σ₂以相同的速率单调上升，直到达到预定的σ₂。

3) 通过伺服控制函数控制x、y方向上两对平行墙体保持σ₂、σ₃恒定，控制z方向平行墙体保证σ₁单调升高，直至峰值(σ_{1, peak})，然后继续加载至σ₁后下降到σ_{1, peak}的70%时终止。

1.3 真三轴模拟应力路径及破坏角的测量

对大理岩数值试样分别进行7组不同应力控制的真三轴压缩破坏试验，即将最小主应力σ₃分为5、20、50、80、100、120和150 MPa 7个应力水平，每个应力水平下，σ₂的大小根据模拟过程中峰值强度确定，涵盖了从σ₂=σ₃到σ₂=σ₁的整个范围(表 1)。

图 3展示出破坏面角度θ的测量方法，左边为颗粒破坏后的fragment显示(同种颜色的颗粒表示同种破碎块体)，中间为微裂纹的空间分布(深色为剪切裂纹，浅色为拉伸裂纹)，右边为颗粒位移矢量空间分布。当σ₂=σ₃时没有明显的破坏角；当σ₂>σ₃时，平行黏结键断裂，无法传递力，故无法形成力链，则在平行黏结键断裂处产生微裂纹，微裂纹的积累导致破裂面的产生。而且，破坏面平行于中间主应力方向，各个应力路径下破坏角的范围为60°~65°(见图 3(b)、表 1)，对于恒定σ₃，破坏角随着σ₂的增大而增大。

2 真三轴压缩数值结果分析 2.1 中间主应力对真三轴破坏应力的影响

表 1列出了破坏时的所有真三轴压缩应力条件，对于每个最小主应力σ₃，σ_{1, peak}是σ₂的函数，如图 4所示。可以看出，σ_{1, peak}随σ₂的增加而增大，且当σ₂达到某一固定平台值时，超过此值σ_{1, peak}逐渐下降，这种变化趋势不是左右对称的，相关硬岩真三轴压缩试验结果一致^[19-20]。

由上可知，σ_{1, peak}既是σ₃的函数，又是σ₂的函数。σ_{1, peak}这种先升后降的趋势符合二阶多项式方程特征。曲线图 4中的虚线组成常规三轴压缩(σ₂=σ₃)和常规三轴拉伸(σ₂=σ₁)两种极限情况，σ_{1, peak}与σ₃的变化可以用下式来线性拟合，这实为传统的三轴破坏判据，即Mohr-Coulomb准则。

对于常规三轴压缩(σ₂=σ₃)

$ {\sigma _{1, {\rm{ peak }}}} = 150.84 + 3.68{\sigma _3}, {R^2} = 0.999 $

(1)

对于常规三轴拉伸(σ₁=σ₂)

$ {\sigma _{1, {\rm{ peak }}}} = 169.28 + 5.40{\sigma _3}, {R^2} = 0.999 $

(2)

对于给定的σ₃，σ_{1, peak}与σ₂的关系曲线表现出强度上的差异性，这揭示了传统的三轴压缩试验以及经典莫尔-库仑和霍克-布朗等破坏准则存在不足。

取峰值应力50%处的割线模量计算岩样的弹性模量E，图 5展示E在σ₃应力水平下随σ₂的变化关系，可以看出，在每个最小主应力水平下，弹性模量随着中间主应力的增加先急速增加，再逐渐趋近于平稳，并且最小主应力越大，平稳增加段所占比例就越大，表明最小主应力对弹性模量的增加有约束作用。采用Logistic函数可直观表示E和σ₂的变化关系，拟合方程列于图 5中，拟合度较高。

2.2 基于八面体理论的真三轴破坏应力分析

八面体理论认为，当岩石内部八面体剪应力达到某一临界值时，岩石将进入破坏阶段。为了研究大理岩试样的空间破坏特性与中间主应力之间的关系，并且表 1和图 4真三轴破坏时的应力状态都可以由两个主应力不变量之间的单一关系表示，即破坏时的八面体剪应力(τ_{oct, f})和破坏时的八面体正应力(σ_{oct, f})：

$ \begin{array}{*{20}{c}} {{\tau _{{\rm{oct }}, {\rm{f}}}} = \left[ {{{\left( {{\sigma _{1, {\rm{ peak }}}} - {\sigma _2}} \right)}^2} + {{\left( {{\sigma _2} - {\sigma _3}} \right)}^2} + } \right.}\\ {{{\left. {{{\left( {{\sigma _3} - {\sigma _{1, {\rm{ peak }}}}} \right)}^2}} \right]}^{1/2}}/3} \end{array} $

(3)

$ {\sigma _{{\rm{oct }}, {\rm{f}}}} = \left( {{\sigma _{1, {\rm{ peak }}}} + {\sigma _2} + {\sigma _3}} \right)/3 $

(4)

图 6(a)是在τ_{oct, f}-σ_{oct, f}区域中重新绘制破坏时的应力条件，可以看出，在大理岩数值试样当中，τ_{oct, f}随着σ_{oct, f}的升高而持续上升，尽管上升速率有所下降。因为所建数值模型没有明显缺陷，破坏时的主要原因是剪切而不是压实^[21]，与2.4节的结论相同。同时，为了更好地确定在恒定σ₃应力条件下σ₂在σ₁和σ₃之间的相对位置，引入应力比b^[17]，对于确定的破坏应力状态，可以获得良好的数据基础。应力比b定义为

$ b = \left( {{\sigma _2} - {\sigma _3}} \right)/\left( {{\sigma _1} - {\sigma _3}} \right) $

(5)

对图 6(a)中所示的数据点进行二阶多项式方程拟合，得到式(6)的τ_{oct, f}-σ_{oct, f}变化方程。

$ {\tau _{{\rm{oct}}, {\rm{f}}}} = 32.47 + 0.62{\sigma _{{\rm{oct}}, {\rm{f}}}} - 0.00017\sigma _{{\rm{oct}}, {\rm{f}}}^2, {R^2} = 0.96 $

(6)

由图 6(b)可以看出，在τ_{oct, f}-σ_{oct, f}曲线区域中数据点表现出一定的分散性，但是分散中包含潜在规律，即对于每个级别的τ_{oct, f}，最小的σ_{oct, f}一致地出现在σ₂=σ₃时，最大的σ_{oct, f}一致地匹配σ₂=σ₁应力状态。

在图 6(c)中，重新绘制了与图 6(a)和(b)相同的应力状态数据点，从式(5)可看出, σ₂=σ₃和σ₂= σ₁时，应力比b=0和b=1。与此同时，图 6(c)为每个数据点分配了与其表示的b值相对应的颜色(侧面显示的颜色条)，易见低σ₃下个体b的离散度非常低，σ₂=σ₃和σ₂=σ₁的拟合曲线(图中虚线)明确地描绘b=0和b=1时的应力状态，即分别对应岩石处于广义三轴压缩状态(σ₁>σ₂=σ₃)和广义三轴拉伸状态(σ₁=σ₂>σ₃)。对比图 6(a)和(c)可以看出，对于0到1的任何其他b值，都可以得到一条拟合良好的曲线，并且对于恒定的b值，τ_{oct, f}随着σ_{oct, f}线性增加，形成线性的破坏包络线。

2.3 变形特征与裂纹扩展的联系

图 7给出了σ₃=150 MPa下的偏应力与3个主应变分量的关系。限于篇幅，仅讨论一个最小主应力水平下的关系图。易见随σ₂/σ₁的增大，应变ε_y由拉应变转为压应变，并且压应变也随之增大。增加σ₂/σ₁值也会导致在最小应力方向上位移膨胀量的提升。在第1个拐点之前，3个主应力、主应变大小相同，压缩方向一致且为正；在第1个拐点与第2个拐点之间，最大主应力和中间主应力大小相同，ε_z和ε_y由于两个方向的压缩而一致且为正，ε_x由于在该方向上开始扩张而为负；第2个拐点之后，ε_y与ε_z由重合变为分离，两个方向上的应变变化较为明显。

图 8展示出σ₃在80 MPa下微裂纹数量、每10时步所产生的微裂纹数和偏应力与轴向应变的关系。对比常规三轴压缩偏应力-应变曲线，真三轴大理岩破坏微裂纹的演化过程也可分为4个阶段：第1个阶段为线弹性阶段(OB段)，其中3个主应力伺服至σ₃时到达A点，数值试样不断被压缩，包含初始压密区，几乎没有微裂纹和声发射事件产生，B点所对应的偏应力在40%σ_{1, peak}附近；第2个阶段为裂纹稳定扩展阶段(BC段)，试样声发射事件呈稳定增长态势，岩石损伤呈稳定速率上升，C点所对应的偏应力在80%σ_{1, peak}附近；第3个阶段为裂纹非稳定扩展阶段(CD段)，D点对应峰值强度，从C点新生裂纹张开度明显增大，岩样出现不可恢复的塑性变形，并且随着σ₂增加该阶段占比增大；第4个阶段为峰后破坏阶段(DE段)，D点之后，不同应力条件下微裂纹的扩展速率不尽相同，但是随着宏观裂隙的形成，该阶段后期的微裂纹扩展速率逐渐降低。

从图 8可以看出，随着σ₂的增大，试样峰后表现出由塑性破坏转向脆性破坏的趋势。当σ₃保持不变时，随着σ₂增大CD阶段所占比重增加，并且剪切裂纹数目慢慢大于拉伸裂纹数目，破坏方式从拉伸破坏转变为拉剪混合破坏。对于恒定的σ₃，不同应力路径下σ₂的大小影响微裂纹的产生速度，即损伤的变化速率，且σ₂越大则峰后损伤速率下降得越快。

2.4 微裂纹形成对宏观破裂面的影响

图 9为不同应力路径下微观裂纹倾向倾角与其相对数量关系的赤平极射投影。对比图 9(a)与(d)可看出，初始围压对微观裂纹的倾向倾角有较大影响：初始围压较低时微观裂纹倾向分布较均匀，裂纹主要平行于最大主应力方向，类似单轴压缩时裂纹的赤平极射投影^[22]；σ₃=σ₂=80 MPa属于初始围压较高的状态，裂纹的倾向开始往中间主应力方向(90°)和最小主应力方向(0°)上集中，并且裂纹的倾角慢慢变大。从图 9(b)和(c)可以看出，因σ₂>σ₃时出现明显的破坏角，裂纹倾向以及数量分布从中间主应力方向逐渐向最小主应力方向扩展，并且随着中间主应力的增大，裂纹倾角也慢慢偏离最大主应力方向，最终扩展的角度形成实际的宏观破坏角(位于60°~65°)。

从图 10(a)可以看出，对于低σ₃应力水平，拉伸裂纹占比较大；当σ₃>50 MPa之后，随着σ₂增大，剪切裂纹占比开始增大，破坏方式从拉伸破坏向拉剪混合破坏转变。在图 10(b)和(c)中，对于恒定的σ₃，微裂纹数目随σ₂先下降然后再上升，岩样损伤发展呈现出“对勾”型变化趋势，即开始下降快，而上升速度滞缓。当σ₃=σ₂时，随着σ₃的增大，微裂纹数目快速上升；而当σ₁=σ₂时，随着σ₃的增大，微裂纹数目上升速率较为缓慢，后期微裂纹数目几乎没有变化。从这两图中的虚线可以看出高中间主应力抑制损伤的上升趋势。

图 11对比了不同真三轴应力路径下大理岩数值模拟破坏形态与试验后的宏观破坏形态^[4]。主裂纹为剪切破坏，次生裂纹为拉伸破坏，二者吻合较好，且与文献[5-6]中大理岩的破坏形态基本相同。此外，图 9与10的结果也体现了大理岩真三轴宏观破坏方式，故本文所建模型是合理的。

3 结论

1) 基于平行黏结模型的三维PFC细观数值模拟，能较好地从微观角度观察大理岩的破坏模式，中间主应力σ₂对岩样的力学响应影响显著，且微裂纹的赤平极射投影可动态展示微裂纹的倾向倾角和数量分布规律。

2) 八面体剪切应力与八面体正应力可准确地拟合出此大理岩真三轴压缩条件下的破坏应力，随着σ₂的增加，大理岩试样破坏模式由拉伸破坏向拉剪破坏转变，并且应力-应变曲线的峰后段由塑性状态逐渐向脆性状态过渡。

3) 随着σ₂的增加，偏应力与3个主应变曲线将出现明显拐点，中间主应变由拉伸向压缩转变；根据应力-应变曲线形状，可将微裂纹演化过程分为4个阶段，岩石损伤演化受σ₂影响显著，呈现出“√”型的变化趋势。