我国地质环境条件复杂，构造活动强烈，在此背景下的工程项目受到诸如动载、高应力等因素的影响^[1]，尤其随着工程不断向深部推进，岩体往往处于一种高围压、三向不等的应力状态（σ₁ ≠ σ₂ ≠ σ₃）^[2]。因此，研究在真三轴围压下岩石的变形失稳和损伤破坏机理对于科学准确地预测和评价地下岩石工程的稳定性，防止重大地质灾害事故的发生具有重要的工程意义^[3]。

分离式霍普金森压杆（split Hopkinson pressure bar, SHPB）已被广泛应用于研究岩石类材料在高应变率下的动态特性，使材料在10⁰～10² s⁻¹的应变率范围内变形，并实现可靠的动态力学响应监测^[4]。相关学者已对此开展了大量研究：Sun等^[5]通过对类岩石材料SHPB试验模拟，讨论了不同脉冲波形对试件动态响应特征的影响，得到了岩石的动强度与应变速率的关系，并提出了类岩石材料的动强度预测模型。刘晓辉等^[6-7]试验研究了煤岩在冲击载荷下的动态力学特性，考察了煤岩在动态破坏下的能量释放特征，指出弹性能的释放程度与煤岩破碎的块度呈反比。Li等^[8]也对裂隙大理岩在单轴应力条件下的力学特性和断裂行为进行了相关试验，而针对类岩石材料在复杂应力下的动力学性能研究，徐松林等^[9]研制了一套真三轴静载下的三维SHPB试验系统，此系统可利用压杆对立方体试样施加冲击载荷，并提供3个方向0～100 MPa的围压静载，实现试件在不同应力环境下的动态力学响应。基于此，也有不少学者进行了相关研究。Luo等^[10]从能量积累耗散的角度研究岩石的渐进破裂机制，综合考虑了地应力、冲击次数和应变率等因素，分析了岩石在动三轴加载下的累积损伤特征。袁良柱等^[11]结合真三轴冲击试验结果，开展数值模拟探讨了混凝土试样的应变率效应，发现应变率的升高会带来混凝土强度的提高，并与横向惯性约束具有较强的耦合作用。Xu等^[12]基于Drucker-Prager强度准则讨论了中间应力对动强度的影响，指出载荷路径和应变率具有明显的相关性，据此得出了随着应变率的提高，试样的摩擦角φ增大，黏聚力c逐渐减小。Chen等^[13]从中主应力效应角度出发，研究高应变率和不同应力状态下混凝土的应变率效应和中间主应力效应，并引入修正的Drucker-Prager模型，描述了中间主应力和加载速率间的关系，获得了较好的试验结果。

综上可知，针对类岩石材料在单轴和三轴SHPB动态试验方面，学者们已开展了大量的工作，而对真三轴SHPB条件下深部大理岩的动态力学特征研究尚显不足，且通过离散元和有限差分耦合方法分析大理岩动力学响应的成果鲜有报道。本文中，在已有研究的基础上，基于PFC (particle flow code)-FLAC (fast Lagrangian analysis of continua)耦合思路，采用平行黏结（parallel bond model, PBM）接触模型，从应力-应变曲线、裂纹与能量演化机理等方面对深部大理岩的动态响应进行深入分析，着重探讨真三轴压缩条件下大理岩的动力学特性和损伤演化规律，以期可为深部岩体工程的稳定性评价及灾害预测提供参考。

1.   模型建立与参数选取

1.1   数值耦合理论

PFC是一种基于离散元法（discrete element method, DEM）的计算软件，在模拟颗粒材料大变形时的破坏形态方面具有明显优势，作用在颗粒上的所有力都由力-位移决定时，一般采用牛顿运动方程的中心差分形式来描述每个颗粒的运动^[14]，如图1所示，x^(A)、x^(B)和x^(C)为节点横坐标；d为两颗粒圆心的距离；R^(A)和R^(B)分别为颗粒A和B的半径。颗粒与颗粒的接触力（法向和切向）与相对位移成线性关系：

图  1  颗粒接触模型

Figure  1.  Particle contact model

下载: 全尺寸图片 幻灯片

式中：n_i为接触面单位法向量；U_n为颗粒法向重叠量；K_n和K_s分别为法向和切向的刚度；{F_i^s}为时间步长∆t开始时存在的剪切力矢量；∆U_s为在单位步长内产生的切向重叠增量。

FLAC是一种基于显式有限差分法（finite difference method, FDM）的软件，通过将连续体划分为四边形单元，将求解域替换为有限数量的网格节点^[15]。为了消除沙漏变形问题，FLAC将每个单元细分为2个重叠的等应变三角形子单元集，如图2所示，${u_i^{\left( {\mathrm{a}} \right)}}$、${u_i^{\left( {\mathrm{b}} \right)}}$为节点速度；∆t为边长；S⁽¹⁾、S⁽²⁾为三角形单元边；${n_i^{\left( 1 \right)}}$、${n_i^{\left(2 \right)}}$为对应边的单位法向向量。作用于单个节点处的力可由三角形单元力得到：

图  2  单元示意图

Figure  2.  Diagrams of elements

下载: 全尺寸图片 幻灯片

式中：$\Sigma$表示连接到该节点的所有子元素的总和，n_iS为三角形单元边的法向矢量力，σ_ij为单元节点力。

考虑到FLAC模型中作用在边界元上的节点力不能直接从PFC模型中得到，给出了一种将面合力转化为节点力的节点力分配方法^[16]，如图3所示，F_wx和F_wy分别为作用在边界墙面不同方向（x和y）上的合力；M_w为绕z轴的合成力矩，逆时针为正；f_1x和f_1y为节点1(x₁, y₁)在FLAC中的等效节点力，f_2x和f_2y为节点2(x₂, y₂)的节点力。

图  3  耦合边界节点力传递模型

Figure  3.  The force transfer model of coupling boundary nodes

下载: 全尺寸图片 幻灯片

节点力与边界墙面力的关系可用静力等效原理描述为：

1.2   细观参数标定与试验验证

本文数值工具含有限差分单元部分和离散元部分。其中，有限差分单元部分包括入射杆、透射杆，离散元部分对应大理岩试样。为保证研究可合理反映室内试验，建立与实际试验尺寸相符的杆件系统^[9]。入射杆和透射杆为方杆，横截面尺寸均为50 mm×50 mm，入射杆长L₁=2.5 m，透射杆长L₂=2.0 m。在有限差分网格设置上，入射杆轴向划分250份，透射杆轴向划分200份，入射杆、透射杆周方向均划分10份，最终总计有限单元45 000个。试样模型尺寸为50 mm×50 mm，其内部均匀性好且无裂隙。整个试样包含21 368个颗粒、73 187个平行黏结键和2 139个线性黏结键。数值建模如图4所示。

图  4  耦合数值模型

Figure  4.  Modelling for coupling simulation

下载: 全尺寸图片 幻灯片

计算分析所有杆件均采用高强度合金钢，其密度为7 900 kg/m³，弹性模量为210 GPa，泊松比为0.25。试样大理岩单轴抗压强度为215.30 MPa，弹性模量为107.74 GPa。

PFC软件中内嵌有多种模型，如CBM模型、PBM模型和SJM模型等。采用这些模型可以针对材料的不同物理力学特性，有选择性地进行模拟标定。其中，PBM模型能较好地模拟试样在受压时内部受到的拉伸和剪切破坏作用，且模拟效果较好，现已被广泛应用于研究岩石类材料的相关问题^[17]。因此，本文中选用PBM模型对大理岩的细观参数进行标定。

在模型初始细观参数中，试样宏观参数与细观参数间的关系较复杂。大理岩试样的弹性模量E主要受颗粒模量${E}_{\mathrm{p}}$和黏结模量$\overline{{E}_{\mathrm{p}}}$的影响，其泊松比ν与颗粒刚度${k}_{\mathrm{n}}$和${k}_{\mathrm{s}}$相关，抗压强度σ_c受黏结法向强度$\overline{{\sigma }_{\mathrm{c}}}$和黏结切向强度$\overline{{\tau }_{\mathrm{c}}}$的影响较大，而对其他参数的敏感度较低^[18-19]。宏细观参数与细观参数间的主要关系如下：

式中：E为试样的弹性模量，$\nu$为试样泊松比，σ_c为试样抗压强度，ϕ_E、$\phi_\nu$和ϕ_c分别为弹性模量、泊松比和抗压强度的相关函数，$\mu$为颗粒摩擦因数。

根据平行黏结模型中宏观参数与细观参数的关系，采用“试错法”对数值模拟结果进行校准微调，得到一组大理岩细观参数，如表1所示，使用该组参数所得计算结果与室内试验结果匹配度较高。

表  1  大理岩细观参数

Table  1.  Microscopic parameters of marble

细观参数	含义	标定值
Rmin/mm	最小颗粒半径	0.9
Rmax/Rmin	最大、最小颗粒半径比	1.4
Ec/GPa	颗粒接触模量	30
kn/ks	颗粒刚度比	1.5
${ \overline{{E}_{\mathrm{c}}}}$/GPa	平行黏结接触模量	10
${\overline{{k}_{\mathrm{n}}}/\overline{{k}_{\mathrm{s}}}}$	平行黏结刚度比	1.5
μ	颗粒摩擦因数	0.5
${ \overline{{\sigma }_{\mathrm{c}}}}$/MPa	黏结法向强度	70
${ \overline{{\tau }_{\mathrm{c}}} }$/MPa	黏结切向强度	44
λ	黏结半径乘子	1

下载: 导出CSV 

| 显示表格

应力-应变曲线的数值模拟结果与室内试验结果的对比如图5所示，可以看出，模拟所得应力-应变曲线与试验曲线吻合较好。以单轴压缩为例，数值模拟所得试样单轴抗压强度为218.75 MPa，弹性模量为101.65 GPa，与室内试验结果的差距较小，在误差允许范围内。值得注意的是，室内单轴试验的应力-应变曲线在峰前阶段存在明显的上凹段，且峰后屈服平台较明显。这是因为，真实的大理岩内部存在许多细小的孔隙，在试验过程中受到压缩作用，这些孔隙发生闭合，且在较高应变率下大理岩材料会由脆性逐渐向延性过渡，与Zhao等^[20]的研究报道类似；而数值模拟中假定试样均匀、内部无缺陷，应力-应变曲线在峰前阶段未呈现上述现象。

图  5  不同压缩条件下试件的应力-应变曲线对比

Figure  5.  Comparison of stress-strain curves of specimens under different compression conditions

下载: 全尺寸图片 幻灯片

在单轴、三轴压缩条件下，大理岩试样的实际破碎情况与模拟结果的对比如图6所示。从实际破碎情况分析，在冲击荷载作用下，试样的主要破碎区在受击端面，特别是单轴试验下的破碎情况尤为明显。在三轴试验中，由于围压的存在，试样迎冲面受损较轻，仅试样的边角处略微破碎。观察可见，数值结果能够较好地表征大理岩试样在单、三轴压缩条件下的实际裂纹分布。综上可知，FDM-DEM理论方法适用于真三轴SHPB试验分析，可沿用该套PFC-FLAC耦合理论进行后续深层次的研究。

图  6  不同压缩条件下试件的破碎形态对比

Figure  6.  Comparison of failure modes of specimens under different compression conditions

下载: 全尺寸图片 幻灯片

2.   系统验证与方案选取

2.1   一维应力传播与动态应力平衡

SHPB测试系统一般要满足以下2个假设：杆中一维应力波假设和试样内部应力均匀化假设^[21]。通过数值模拟对3个主应力方向预压力状态为（5 MPa，10 MPa，15 MPa）的岩样进行动态真三轴试验，选取入射弹性波应力峰值为200 MPa。以入射杆左端500 mm处为起点，每间隔150 mm选取1、2、3等3个测点，用来监测经过该测点的应力波信号。

图7为3个测点处记录下来的应力波信号。从图7可以得知，3个测点处记录的应力波波形、幅值大小完全一致；入射波从测点3传播到测点1过程中没有发生应力波的能量损耗，且应力波通过试样部分变为反射波又传回入射杆，反射波从测点1到3同样未发生衰减，表明波在弹性杆件中传播时，其应力衰减可忽略，因此该数值模型可以满足一维应力波假设。图8为动态压缩下试样的动态应力平衡曲线，其中，I表示入射波，R表示反射波，T表示透射波。根据曲线数据处理，入、反射波叠加曲线I+R与透射波T重合度较高，且在峰前部分几乎完全重合，表明在试验全过程中，试样两侧的应力处于平衡状态，测试结果有效。

图  7  3个测点应力信号记录

Figure  7.  Recorded stress signals at three measuring points

下载: 全尺寸图片 幻灯片

图  8  动态应力平衡图

Figure  8.  Diagram of dynamic stress balance

下载: 全尺寸图片 幻灯片

2.2   模拟方案

为了深入研究复杂应力条件下深部大理岩受动荷载影响的力学特性，通过控制入射应力和围压来模拟不同工况下的试验加载。其中，对入射应力的模拟采用直接施加半正弦应力波的方式，采用3种幅值（σ_in 为150、175和200 MPa）来模拟不同入射应力下的冲击试验；沿着冲击方向的预静态压力（下简称“轴压”，σ₁）分别设定为0、5、10、15和20 MPa。模拟加载由以下3个阶段组成：

(1)在静水条件下（σ₁ = σ₂ = σ₃），通过控制6个伺服墙使试样处于预定围压下；

(2)通过wall-ball耦合代码，建立有限元与离散元间的连接，使耦合界面处单元与颗粒间存在力与速度的传递通道；

(3)在入射杆的端部施加预定大小的半正弦应力波，并设置模型动态求解时间开始运行，直至系统完成计算。

3.   数值结果分析

3.1   动态应力-应变曲线

图9为在3种不同入射应力σ_in下，通过轴压σ₁得到的试样应力-应变曲线。由图9(a)可知，在3种入射应力下试样的峰值强度依次为115.2、131.3和138.7 MPa。其中，当入射应力σ_in为150 MPa时，峰后段应变未继续增大，而是随着应力的降低而逐渐减小，产生了明显的“回弹”。其原因是：在试验加载阶段应变能存贮在试样中，由于入射应力较低且未达到试样的动态破坏强度，在峰后段存贮的应变能逐渐释放。随着入射应力的提高，试样应力-应变曲线的回弹现象逐渐消失，当入射应力达到200 MPa时，表现较明显的应变软化特征。

图  9  不同轴压试样在不同入射应力下的应力-应变曲线

Figure  9.  Stress-strain curves of the specimens in different axial compression states (σ₁, 5 MPa, 10 MPa) under different incident pressures

下载: 全尺寸图片 幻灯片

当轴压σ₁从0逐渐增大至20 MPa，试样在3种入射应力下的峰值强度呈逐渐下降趋势。以入射应力为200 MPa时为例，试样在3种轴压σ₁下的峰值应力分别为138.7、126.4和117.2 MPa，随着轴压的升高而降低。此外，观察图9中的曲线可知，轴压相同时，3组入射应力下试样应力-应变曲线的峰前段几乎重合；在加载阶段，试样的弹性模量未随着入射应力的变化而发生明显改变。

3.2   峰值应力、应变与预压力和入射应力的关系

应变率效应是指材料的性质随其形变量的变化而发生改变的现象，其定义为应变变化的时间效应^[22]。图10为预压力状态为(0 MPa, 5 MPa, 10 MPa)时在3种入射应力下试样的应变时程曲线。由图10可知，试样峰值应变受冲击荷载作用表现出了明显的应变率效应。随着入射应力由低到高，试样的峰值应变依次为6.45×10⁻³、7.96×10⁻³和11.8×10⁻³，呈现逐渐增大的趋势，且增大幅度也随入射应力的提高而变大。此外，与入射应力为175和200 MPa时不同，入射应力为150 MPa时，试样的应变随着时间的推移逐渐减小为0。这是因为，较低的入射应力未能达到试样的破坏强度，试样在加载后阶段仍然保持弹性特性，产生的形变随着时间的推移而逐渐恢复。

图  10  三轴围压(0 MPa, 5 MPa, 10 MPa)试样在3种入射应力下的应变时程曲线

Figure  10.  Strain-time histories of specimens in triaxial compression (0 MPa, 5 MPa, 10 MPa) under three different incident pressures

下载: 全尺寸图片 幻灯片

图11为3种入射应力下岩样的峰值应力和峰值应变随轴压σ₁的变化曲线。峰值应变是指峰值应力对应的应变。此组数据通过改变轴向压力σ₁（0、5、10、15、20 MPa）来模拟大理岩试样在三向围压(σ₁, 5 MPa, 10 MPa)状态下的动力学变化特性。图11显示，在峰值应力（或动态压缩强度）方面，以σ₁ = 5 MPa为例，随着入射应力的提高，大理岩试样的动态压缩强度逐渐升高，呈现明显的应变率增强效应。入射应力从150 MPa升高到175 MPa，大理岩试样的峰值应力升高幅值为16.4 MPa，显著大于入射应力从175MPa升高到200 MPa时试样峰值应力的增幅7.1 MPa。这是因为，高入射应力下，试样裂纹的萌生和扩展要滞后于载荷的增大^[6]；应变能存贮在试样内，表现出来就是动态压缩强度随加载速率的提高而提高，而增幅随应变率的升高而逐渐变小。此外，试样的动态压缩强度也随轴压σ₁的升高而呈逐渐降低的趋势，其大致呈负相关。在峰值应变方面，同等条件下由于应力的升高必然导致试样产生更大的变形^[23]，因此其变化规律与峰值应力有相关性。

图  11  不同入射压力下试样的峰值应力和峰值应变随轴压的变化

Figure  11.  Vairations of peak stresses and peak strains of specimens under different incident pressures with axial pressure

下载: 全尺寸图片 幻灯片

3.3   相同入射应力下预压力对强度的影响

以上分析了仅改变轴压${\sigma }_{1}$的情况下，大理岩试样在真三轴条件下的动态力学特性。本节中将综合探讨在入射应力一定的情况下（200 MPa），3个方向围压（${\sigma }_{1}, {\sigma }_{2}$,${\sigma }_{3}$）对大理岩试样力学特性的影响，其中以试样的动态峰值破坏强度为主要研究对象。

在三轴静载试验中，围压增大会导致试样的强度有所提高，而随着中间主应力${\sigma }_{2}$从${\sigma }_{2}$=${\sigma }_{3}$增大到${\sigma }_{2}$=${\sigma }_{1}$时，试样强度呈现先升高后降低的变化规律^[24]。在三轴冲击试验中，试样强度受围压的影响与其方向相关。图12为不同约束条件下围压对大理岩试样动态压缩强度的影响，初始预压力选取为5、10、15、20 MPa。从图12得知，改变冲击方向预压力${\sigma }_{1}$（即轴压）时，试样动态压缩强度的变化规律与改变横向预压力${\sigma }_{2}$、纵向预压力${\sigma }_{3}$时的不一致。改变轴向压力时，试样动态压缩强度随轴向压力${\sigma }_{1}$的增大呈线性降低，当${\sigma }_{1}$从5 MPa增大至20 MPa时，试样峰值应力从130.5 MPa降低到117.8 MPa，降低了9.7%。这主要是因为，在冲击方向上施加的预压增大，诱导该方向上裂纹萌生与扩展，导致试样的抗压强度降低。然而，横向压力${\sigma }_{2}$和纵向压力${\sigma }_{3}$的影响趋势与轴向压力${\sigma }_{1}$的相反。以纵向压力为例，${\sigma }_{3}$的增大导致了大理岩动态强度的提高，当${\sigma }_{3}$从5 MPa增大至20 MPa时，试样动态压缩强度从113.7 MPa升到124.6 MPa，提升了9.6%。从机理上分析，纵向预压${\sigma }_{3}$提供了垂直于冲击方向的约束，抑制了岩样内裂纹的萌发、扩展和汇聚，一定程度上提高了抗压强度。由于上图模拟的工况偏少，动态抗压强度与侧向围压间的变化规律有待进一步探讨。

图  12  不同方向预压力对岩样峰值强度的影响

Figure  12.  Effect of pre-pressures in different directions on the peak strength of specimens

下载: 全尺寸图片 幻灯片

3.4   破碎形态分析

PFC软件能用于监测岩石真三轴试验过程中的声发射现象，软件内部的离散裂隙网络（discrete fracture network, DFN）模块可模拟岩石在试验过程中的破裂损伤状态，引用fracture代码文件可在胶结模型的两端生成一个不参与力学计算的微裂纹标识，代表该胶结处达到受力极限状态破坏并产生裂隙^[25]。在本文研究中，将单位时间内试样内部黏结破坏所形成的微裂纹当成一次声发射事件。

图13为在2种入射应力（175、200 MPa）下监测到的大理岩试样内部声发射事件数量和应力-应变曲线关系。可以发现，在峰值应力段之前声发射事件处于一个缓慢积攒期，在这个时间段由于试样内部应力逐渐升高，裂纹开始萌生。此后，声发射事件开始增多，达到峰值强度后，微裂纹相互贯通、扩展，声发射事件大量发生并达到峰值，称为峰后爆发期。对比图13(a)和13(b)发现，当入射应力较高时，峰前段单位时间内的声发射事件普遍较多。此外，峰后段声发射试件数量峰值也更大：175 MPa入射应力下声发射事件数量峰值为50左右，而200 MPa入射应力下声发射事件数量峰值接近70，明显大于前者。

图  13  不同入射应力下大理岩试样内部声发射事件数量和应力-应变曲线关系

Figure  13.  Numbers of acoustic emission events in specimens and stress-strain curves under different incident stresses

下载: 全尺寸图片 幻灯片

为揭示真三轴条件对受压试样破坏模式的影响，选取围压(5 MPa, 10 MPa, 15 MPa)试验条件（入射应力200 MPa）展示三轴承压岩石试样的渐进式破坏过程，如图14～15所示。不同破坏形式以不同的颜色加以区分，红色和黄色裂纹分别表示拉伸和剪切破坏。由图14可知，随着入射应力波的加载，试样冲击加载面出现裂纹并逐渐向内扩展。图15为试样的正视图，试样的左侧为冲击加载面，荷载从这个面将应力波传递至试样中。应力波透过试样后，左侧最先受压并出现裂纹，在520 μs时刻，应力逐渐向中心区域集中，在试样中部出现微裂纹。在541 μs时刻，应力集中区进一步变宽，聚集在冲击荷载面和试样内部的一条斜裂纹分布带。随着冲击应力持续施加，微裂纹进一步扩展和汇聚，形成了宏观层面上的破坏形态。

图  14  不同时刻的岩样破坏形貌

Figure  14.  Failure morphologies of rock specimens at different times

下载: 全尺寸图片 幻灯片

图  15  岩样破坏正视图

Figure  15.  Front views of failed specimens

下载: 全尺寸图片 幻灯片

图16为试样在三轴围压（5 MPa，10 MPa，15 MPa）下、采用200 MPa入射应力加载过程中应力和裂纹数目演化时程曲线。可以看出，在真三轴围压条件下试样的破坏模式主要以张拉作用形成的拉伸破坏为主，剪切裂纹相对较少。图16(a)显示，在应力峰后的阶段，岩样处于裂纹爆发期，裂纹从萌生、扩展到最终稳定经历时间很短。记录得到的拉伸裂纹有26700条，占总裂纹数的83.6%；剪切裂纹为5250条，占总裂纹数的16.4%。由图16(b)可以看出，试样首先由张拉破坏形成拉伸裂纹，550 μs后才开始有压剪裂纹出现。

图  16  不同类型裂纹演化曲线

Figure  16.  Evolution curves of different types of cracks

下载: 全尺寸图片 幻灯片

3.5   能量演化分析

荷载作用下岩样内部能量跟踪对分析岩石破坏机理、推导基于能量的本构参数和准则具有重要意义^[26]。在SHPB试验中，由于试样内部裂纹形成、扩展而损耗的能量，可看作入射波的弹性波能减去透射能和反射能之和。真三轴SHPB系统中试样耗散的能量包括相对滑动的摩擦能、碎片的动能、产生新破裂面累积的破裂能以及其他形式的能量，如声音和热量等，实际试验当中这些能量难以直接测量和划分。在动态载荷作用下，冲击过程中记录了5种能量：(1)杆在入射方向上施加的总能量，即输入能W；(2)储存在相邻颗粒的弹性键和黏结接触中的应变能E_s；(3)试样中颗粒因接触间摩擦和阻尼作用而损耗的耗散能E_d；(4)颗粒运动而产生的动能E_k；(5)打破颗粒间黏结接触键所需的胶结破坏能E_ce。其中，耗散能E_d可进一步划分为摩擦能E_df和阻尼能E_dd，应变能E_s可进一步划分为颗粒应变能E_sp和胶结应变能E_sce。各能量有如下关系：

图17为冲击载荷200 MPa下应力和能量的时程曲线。依据图17(a)，在峰前阶段，随着动载荷的增大，大部分输入功以应变能的形式存储在试样中。当达到峰值应力点后，试样破碎，应变能逐渐释放。同时，耗散能、动能和胶结破坏能开始上升，结合图16分析，这是由于峰后阶段试样产生大量裂纹，形成各种碎块需要消耗较多能量，且试样变形导致动能逐渐增大。值得注意的是，动应力峰值时刻应变能也几乎达到峰值，随后应力下降伴随着应变能减少，这表明存储的应变能部分沿压杆方向释放。图17(b)展示了应变能和耗散能的能量细分组成。在应变能方面，由于试样接触模型中以平行黏结接触占据主导。因此，在峰前加载阶段，能量积蓄部分以胶结应变能居多，占总应变能的62.4%；在耗散能方面，其主要组成为颗粒位置移动引起的接触摩擦能，岩样内部阻尼能受颗粒间重叠量变化的影响。

图  17  应力和能量时程曲线

Figure  17.  Curves of stress and energy time histories

下载: 全尺寸图片 幻灯片

4.   结　论

聚焦深部大理岩的动力学特性，采用离散元与有限差分耦合法，对真三轴SHPB试验中不同弹速和不同预压条件下岩样的峰值强度、裂纹数、声发射事件以及能量演化规律等进行了深入研究，得到的主要结论如下。

(1)基于PFC-FLAC耦合法模拟得到的大理岩应力-应变曲线、裂纹分布特征与室内真三轴SHPB试验结果吻合度较高，且其结果满足一维应力波假设和应力均匀化假设，该耦合计算方法可适用于复杂应力状态下大理岩的真三轴动态特性数值研究。

(2)大理岩的峰值强度和峰值应变随轴向静压力σ₁的提高而降低；当σ₁值一定时，峰值强度又随入射应力的升高呈现上升趋势；相同入射应力下，轴向压力σ₁对大理岩的强度有削弱作用，而侧向围压σ₂和σ₃能在一定程度上提高试样的峰值强度。

(3)试样声发射爆发期稍滞后于应力峰值，试样中裂纹以拉伸裂纹为主，占总裂纹数的83.6%，并可见明显斜剪切带；应力峰值时刻试样达到应变储能极限，以占比62.4%的胶结应变能为主，随着峰后应力的降低，应变能释放转化为耗散能或动能等。