节理是岩体在长期地质运动中形成的一种结构面。节理的存在影响着岩体的结构和力学特性。在岩土工程中广泛应用的爆破技术以及地震、岩爆等地质灾害均会产生应力波，而应力波在节理岩体中的传播与在连续介质中的传播完全不同。节理会导致应力波在传播过程中发生衰减，加剧能量的耗散。因此，研究节理对应力波传播过程的影响对指导工程实践、优化爆破方案和防灾抗震等具有重要作用。

为此，许多学者从理论和试验方面对应力波在节理裂隙中的传播问题进行了一系列的研究。Zhao et al.（2001）和俞缙等（2007，2008）在线弹性位移不连续理论（Schoenberg，1980；Pyrak-nolte et al.，1987）的基础上发展改进了节理非线性位移不连续理论（王卫华等，2006），并对节理的相关参数进行了研究。节理面的粗糙程度与节理岩体的动态力学特性及应力波在节理面的传播规律密切相关。为了模拟实际节理面，学者们进行了大量室内试验研究。Barton et al.（1977）通过倾斜和拉压等试验提出用粗糙度系数（JRC）描述节理面的起伏程度，并给出10种典型的节理剖面，对应JRC取值范围为0~20；谢和平（1992）提出用分形几何的概念描述节理的粗糙度，量化了节理的粗糙程度；Li et al.（2011）考虑分形几何描述节理粗糙程度，进一步建立了分形损伤节理模型，利用激光表面仪和SHPB试验系统研究了节理面分形维数与应力波透射、反射系数的关系；Li et al.（2017）采用SHPB试验系统研究了节理吻合系数（JMC）（赵坚，1997）对岩石节理动态特性及应力波传播的影响。在模拟试验方面，Resende et al.（2010）利用PFC软件模拟应力波与节理岩体的相互作用，研究了应力波在岩体中传播的微观过程及过程中节理几何形状的变化；曾超等（2018）应用PFC软件研究了节理接触面积和节理厚度对应力波传播特性的影响。一些学者通过在PFC软件中建立SHPB系统数值模型，对应力波传播的微观过程进行了研究。Li et al.（2014）建立了基于PFC的SHPB系统数值模型，对完整试样进行了不同速度下的冲击试验，研究了应力波在霍普金森杆中的传播过程以及试样的破坏过程，从微观上揭示了SHPB试验中试样的动态响应； Zhou et al.（2017）运用SHPB试验系统颗粒流模型，研究了岩石在冲击载荷下的破坏过程及内在机理；Qiu et al.（2017）运用PFC软件，对不同倾角的层状岩石进行了SHPB数值冲击试验，研究了层状岩石的动态破裂特性。以上研究没有进一步分析应力波在不同粗糙度节理中传播的微观机理，也没有对比不同粗糙度的节理试样微观破坏特征。

本文利用PFC软件建立了SHPB系统的颗粒流模型，对不同粗糙度系数的节理岩样进行数值冲击试验，研究了应力波在节理岩石中的传播过程，分析了试样两端的应力平衡关系及节理粗糙度系数对应力波透射系数的影响规律，从微观角度比较了不同粗糙度节理的裂纹扩展和破坏机理。通过研究应力波在粗糙节理试样中的传播规律及其对试样的微观破坏机理，进而为岩土工程设计提供参考。

根据中南大学SHPB试验装置各主要部件的几何尺寸（李夕兵，2014），在PFC^2D中构建SHPB数值模型，对SHPB系统中包括入射杆、透射杆、纺锤形冲头和试样在内的各主要部件建模。为保证结果的准确性，入射杆的长度需确保应力波可在试样中反射9~10次，且该时间内，冲击端产生的弯曲波不传到试样中。同时，入射杆和透射杆的长度需要保证初始的反射波和透射波在相应端面产生的反射波不会影响试样的应力平衡。据此，模型中入射杆和透射杆长度分别设定为2.0 m和1.5 m，直径均为50 mm。纺锤形冲头尺寸与实验室中的一致。数值模拟中赋予组成冲头的所有颗粒相同的速度，使冲头作为整体运动。综合考虑数值模型的精度和计算机性能，选择颗粒半径范围为0.9~3.0 mm，颗粒间孔隙率为0.12，颗粒半径服从均匀分布。实验室中SHPB系统部件的密度为7 800 kg/m³，根据质量守恒定律计算，数值模型的颗粒密度为7 894.7 kg/m³。将模拟过程设置为在透射杆末端产生的拉伸波到达试样之前停止数值计算，从而省去了吸收装置的建模过程。为了在数值试验中监测应力波的传播过程，分别在入射杆上设置4个测量圆（A、B、C、D），在透射杆上设置2个测量圆（E、F）。

参考Li et al.（2019）研究中试样的尺寸，由PFC生成完整试样模型（φ50 mm×60 mm），模型包含的颗粒总数为5 855个，颗粒半径范围为0.30~0.45 mm，孔隙率为0.1。在完整试样的基础上，按照不同的JRC值生成不同的节理面，建立不同JRC值的节理试样模型。为了使应力波通过接触面时接触面应力均匀，将透射杆、入射杆、冲头以及试样的边界均设置为一列排列规则、半径为0.6 mm的颗粒，细观参数与对应部件内部颗粒一致。包含完整试样和节理试样的SHPB系统颗粒流模型见图1。

采用试错法对SHPB数值模型的细观参数进行标定。SHPB系统各部分采用接触黏结（CB）模型（Cundall et al.，2004）。实际上，系统在试验过程中几乎不损坏，因此，设定黏结强度为1×10¹⁰⁰ MPa。弹性杆和冲头的阻尼系数设置为零，该参数设置保证了应力波在弹性杆中有效传播而不产生明显衰减。采用文献（Zhou et al.， 2017）中对同一套SHPB系统部件的细观参数标定结果，得到法向刚度和剪切刚度分别为6.68×10¹¹ N/m和2.45×10¹¹ N/m。SHPB试验装置的颗粒流模型的细观参数标定结果见表1。当冲击速度为10.45 m/s时，物理试验测得的入射波与数值模拟测得的入射波的比较如图2所示。从图中可以看出，物理试验测得的入射波与数值模拟测得的入射波曲线形状几乎相同，说明数值模型可以代替实际试验装置对岩石的动态力学性质进行研究，即验证了颗粒流模型模拟SHPB系统的可行性。

试验中采用红砂岩材料，将其加工成φ50 mm×60 mm的圆柱体试样，试样端部打磨光滑，平整度小于0.02 mm，数量若干。将部分试样切割成等长的2个部分，作为JRC值为0~2的节理岩石试样。将红砂岩试样赋予平行黏结（PB）模型（石崇等，2018），对完整试样模型细观参数进行标定。首先依据红砂岩的杨氏模量和泊松比等物理力学性质，通过单轴压缩和单轴拉伸数值试验，初步标定试样的微观参数。其次，通过在SHPB颗粒流模型中对试样反复进行数值冲击试验，不断调整试样微观参数，直至模拟波形与物理试验波形一致。最终校准的完整试样细观参数见表2，波形结果对比见图3。

对含有节理的试样模型细观参数进行标定时，将光滑节理（SJ）模型（Mehranpour et al.，2017）赋予节理面两侧颗粒，非节理面颗粒仍然赋予PB模型。然后根据其他学者的研究（Mehranpour et al.，2017；胡万瑞，2017）初步选定节理面参数，反复对节理试样进行模拟冲击试验，不断调整模型参数，直至数值模拟得到的应力波波形与物理试验所得波形一致。最终校准的节理试样细观参数见表3，波形结果对比见图4。

通过数值模拟的方法研究节理粗糙度对应力波传播的影响，建立具有不同JRC值的节理试样数值模型如图5所示。为保证试样在冲击过程中不发生破坏，即在试样处于弹性变形阶段研究应力波的传播规律，将冲头速度设置为3 m/s，分别冲击完整试样和节理试样。

在SHPB试验中，保证试样两端应力平衡是试验数据真实可靠的前提条件。因此，在模拟试验中，通过FISH语言自定义函数分别实时监测试样两端最外侧一列规则颗粒的应力变化，计算式如下：

式中：σ_SI和σ_ST分别为试样两端的实时入射端应力和透射端应力；F_Ij和F_Ti分别为试样中与入射杆和透射杆端面第j、i个接触颗粒的轴向接触力；N_I和N_T分别为试样与入射杆、透射杆端面接触的颗粒总数目；r和δ分别为试样的半径和厚度。

采用应力平衡因子η（Zhou et al.，2017）来刻画试样两端应力平衡情况，计算公式如下：

图6为完整试样、JRC值为4~6的节理试样和JRC值为16~18的节理试样入射端、透射端应力和应力平衡因子随时间的变化曲线。从图中可以看出，应力波在各试样内传播的过程中，应力平衡因子经历了相似的不同变化阶段，以t₁到t₅时刻为特征点，将试样应力状态划分为5个变化阶段，描述如下：①在t₁~t₂时间段，应力波开始到达试样入射端，并向透射端传播，应力波传播过程中，入射端应力不断增大，透射端应力为0 MPa，因此，应力平衡因子维持2不变；②在t₂~t₃时间段，应力波开始到达试样透射端，并在透射端发生透反射，反射应力波向入射端传播，在此时间内应力波在试样中往复传播，在节理面、入射端和透射端不断发生透反射，应力平衡因子逐渐减小，并趋于0；③在t₃~t₄时间段，从应力平衡因子首次达到零值，到试样两端应力达到峰值，整个过程中应力平衡因子波动起伏较缓，说明试样到达稳定的应力平衡状态；④在t₄~t₅时间段，应力峰值过后，应力平衡因子仍能维持在零值上下，但波动明显强于上一时间段，在应力接近0 MPa时，虽然两端应力差值逐渐减小，但比值相对较大，显现出应力平衡因子的剧烈波动，这一时段应力平衡因子的变化说明在峰后阶段试样未受到破坏，仍能够作为一个整体承受应力载荷。⑤在t₅之后的时间段，应力平衡因子为2，说明试样透射端已与透射杆分离，透射端应力为0，入射端仍有部分颗粒与入射杆接触。

图6（a）与图6（b）和图6（c）对比，节理试样t₁~t₂阶段持续时间为28 μs，完整试样t₁~t₂阶段持续时间为25 μs。在节理试样中，应力波首次传播到节理面时，存在使节理面发生闭合的过程，使得t₁~t₂阶段持续时间比完整试样更长。在t₂~t₃阶段，从完整试样到JRC值为4~6的节理试样再到JRC值为16~18的节理试样，入射端应力曲线越来越平缓，是因为节理面反射的应力波比透射端的反射波更早到达入射端，与入射应力波叠加，使节理试样入射端应力的增长比完整试样更缓慢，且节理粗糙度越大，节理面的反射波越强，应力增长放缓的效应越明显。

通过求解应力波穿过节理试样的透射系数来研究节理对应力波传播的影响，透射系数T的定义是：

式中：T为透射系数；max|ε_T|和max|ε_I|分别为最大透射脉冲和最大入射脉冲。

图7为相同冲击速度下，通过透射杆及入射杆上的测量圆所获得的完整试样和不同JRC节理试样中应力随时间的变化过程，表4为完整试样和不同JRC节理试样在3 m/s冲击速度下的透射系数。从图7和表4中可以看出，由于节理面对应力波的衰减作用，节理试样比完整试样的透射系数更小；节理JRC值越大，节理面越粗糙，透射系数越小。这与物理试验得到的节理对应力波传播的影响规律相同。

在SHPB系统颗粒流模型中，将冲头速度设置为11 m/s，对不同JRC值的节理试样进行模拟冲击试验，使试样发生破坏，研究裂纹扩展和试样破坏的微观机理。表5为具有不同JRC值的节理试样裂纹数量统计。图8为节理试样内部应力及裂纹数量随时间的变化情况。图9以JRC值为4~6的节理试样和JRC值为16~18的节理试样为代表，记录了应力波在节理试样传播过程中不同时刻的微观裂纹分布及微观颗粒接触力链分布。t₁时刻，两节理试样透射端接触界面颗粒开始出现接触力，入射端颗粒接触力明显大于透射端，试样内部颗粒间接触力不均匀，试样尚未产生裂纹；t₂时刻，节理试样颗粒内部接触力基本均匀，试样开始产生裂纹；t₃时刻，节理试样处于峰值应力状态，试样出现少量的微裂隙；t₄时刻，节理试样少量内部颗粒间接触力不均匀，微裂隙数量明显增多；t₅时刻，节理试样内部颗粒间接触力不均匀，可以看到明显的宏观裂缝；t₆时刻，节理试样内部颗粒间仅有极少接触力，试样完全破坏，裂纹数目几乎不再增长。

从图8和图9中可以看出：

（1）入射应力相同，节理形貌不同，节理试样的力学特性、裂纹扩展和破坏形态各不相同。

（2）在冲击荷载下，由于节理的动态强度较低，裂纹先在节理面处形成，随着荷载的继续作用，裂纹逐渐从节理面向试样完整部分尤其是试样端面扩展，裂纹的合并和贯通形成宏观裂隙，导致试样最终的破坏，这与实际试验中观察的现象一致。

（3）试样的力链分布情况可以反映试样两端及内部的应力情况，在应力波刚作用于入射端时，应力波还未传播至透射端，入射端的应力明显大于透射端。随着应力波在试样中的透反射，试样两端的接触力趋于相等，试样两端逐渐达到应力平衡，且一直维持到峰后阶段。

（4）由表5可知，在到达峰值应力前，节理试样内部出现少量裂纹，在峰后阶段，张拉裂纹和剪切裂纹数量快速增长；在整个加载过程中，裂纹主要在峰后阶段形成，张拉裂纹数量远多于剪切裂纹，且张拉裂纹占主导地位，说明节理试样主要以拉伸劈裂破坏为主，与实际试验中节理试样的破坏形态一致。

（5）由表5可知，当入射应力相同时，节理越粗糙，节理的动态强度越低，节理试样越容易破坏，产生的裂纹数目越多。节理试样JRC值从4~6增加至16~18时，节理试样的强度随之降低，最终产生的裂纹数目则由2 050个增长至2 377个。

（1）通过分析试样在数值冲击试验中两端的应力分布和平衡关系，发现当应力波传到试样并作用一段时间后，试样两端的应力达到平衡状态，并且平衡状态一直维持到峰后阶段，进一步证明采用PFC^2D模拟SHPB试验的可行性以及模拟计算结果的准确性。

（2）对比分析节理试样和完整试样在数值冲击试验中两端应力变化的典型阶段，发现应力波首次从入射端向透射端传播过程中，节理面发生闭合，延迟了波的传播；相较于完整试样，由于节理面的存在，节理试样比完整试样存在更多的透反射，初次到达入射端的反射波更早，与入射波叠加，使入射端应力增长更缓慢，且JRC值越大，节理面的反射波越强，应力增长放缓的效应越明显。

（3）对具有不同JRC值的节理试样进行相同速度的数值冲击试验，发现完整试样的透射系数大于节理试样的透射系数，节理的粗糙度对应力波的传播有影响，且JRC值越大，节理面越粗糙，透过节理的应力波越弱，透射系数越小。

（4）通过分析具有不同形貌的节理试样在冲击荷载下的裂纹扩展过程，发现裂纹的演化控制着节理试样的破坏过程和最终破坏形态。在到达峰值应力前，试样内部出现少许张拉裂纹和剪切裂纹，到达峰值应力后，试样内部张拉裂纹和剪切裂纹快速增长。在整个加载过程中，裂纹主要在峰后阶段形成，张拉裂纹数量远多于剪切裂纹，且张拉裂纹占主导地位，说明节理试样主要以拉伸劈裂破坏为主，与实际试验中节理试样的破坏形态一致。裂纹起始于节理面处，并逐渐扩展、合并和贯通形成宏观裂隙，导致节理试样最终破坏。节理面越粗糙，动态强度越低，试样越容易破坏，产生的裂纹数越多。