广泛存在于岩体中的节理是岩体在长期地质作用下形成的一种没有位移或微量位移的断裂构造（夏才初等，2002）。在外力作用下，节理会产生闭合、滑移等非线性力学行为，使岩体呈现出不连续性、非均质性、各向异性和非线性等特征（沈明荣等，2015），也会导致应力波波幅和能量的衰减、传播延迟及高频滤波。国内外学者从理论、试验和数值模拟3个方面对应力波在节理岩体中的传播规律进行了研究。

在理论研究方面，Miller（1977）和Schoenberg（1998）采用位移不连续理论来计算应力波在节理处的传播。Pyrak-Nolte（1988，1996） 根据位移不连续理论推导出应力波在单一线性变形节理处传播的透反射解析解。Crouch（1976）采用位移不连续理论对应力波在线性变形节理处的传播进行研究。Cai et al.（2000）和Zhao et al.（2001）采用非线性变形不连续模型与特征值法相结合的方法，对P波垂直入射到单一法向非线性变形干性节理的问题进行研究，获得了节理透反射系数的数值解。

在试验研究方面，多数学者采用分离式霍普金森压杆测试系统（SHPB）研究节理倾角、节理厚度、节理面几何性质和充填介质对应力波传播的影响。田振农等（2008）通过试验研究了应力波在一维节理岩体中的传播规律。Huang et al.（2016）采用SHPB系统，研究了大振幅应力波在充填节理岩体中的传播规律。Li et al.（2011）利用SHPB系统获得了单一节理岩石试样的透反射系数，试验结果与理论值相吻合。李娜娜等（2015）采用SHPB试验系统研究了岩体节理接触面积对应力波传播规律的影响。Li et al.（2019）采用SHPB系统研究了不同角度的单平面节理对应力波的影响，发现节理角度越大，波的衰减越明显，节理的动态强度和吸能率先减小后增大。颜鑫（2021）利用SHPB系统获得了岩石动态抗压强度、动态应力—应变关系、破碎破裂特征与能量耗散规律，揭示了充填节理介质特性对节理岩体应力波传播规律的影响。

与试验研究相比，数值模拟方法具有高效、经济且重复性好等优点，因此被广泛应用于节理岩体中应力波传播问题的研究。颗粒流离散元法（Particle Flow Code，PFC）从细观角度采用离散颗粒和不同黏结模型模拟岩石颗粒间的胶结状态，通过牛顿第二定律以及内部的时间步迭代积分的方式求解位移、接触力和速度。因而采用PFC研究应力波在节理处的传播问题时，不但能够直观显示应力波的传播过程，而且能够分析动载下岩体节理的动力响应过程，揭示应力波与节理相互作用的细观力学机理。Li et al.（2014）基于PFC建立了SHPB的数值测试系统，研究了应力波在霍普金森杆中的传播过程以及试样的破坏过程，从细观上揭示了SHPB试验中试样的动态响应。曾超等（2018）利用PFC研究了节理接触面积和节理厚度对应力波传播特性的影响。Resende et al.（2010）利用PFC模拟应力波在节理处的传播过程，分析透反射波和节理的几何变化情况。许鹏（2017）利用PFC程序模拟了爆炸应力波与不同角度节理相互作用的动态过程。

综上所述，国内外学者对应力波与节理相互作用的研究已取得诸多成果，但有关应力波与节理相互作用动力响应过程及细观力学机理方面的研究很少，需要进一步深入研究。因此，本文首先建立了非线性变形节理的颗粒流模型，然后采用PFC方法模拟分析了应力波在非线性变形节理处传播的透反射规律，最后通过分析节理颗粒运动、变形闭合、节理颗粒接触力链及接触数目的演化特征揭示应力波与节理相互作用的机理。

PFC中需要通过标定细观参数使得颗粒流模型反映出来的宏观参数符合实际情况，在进行模拟计算之前需要赋予岩石模型特定的接触模型，并通过不断试错将得到的模拟结果与室内试验结果进行对比验证，找到一组合适的细观参数用于本次的模拟计算。

完整花岗岩试样单轴压缩试验中所用试样是直径为50 mm，高度为100 mm的圆柱体，岩石密度为2 700 kg/m³。运用PFC^2D建立完整试样的二维数值模型，模型尺寸为50 mm×100 mm（宽×高），保证在单轴压缩时模型颗粒数在适当范围内能加快计算速度并保证一定的精度，花岗岩典型的颗粒尺寸大于0.5 mm，因此选择平均半径稍大于0.25 mm的颗粒尺寸建立模型，能够很好地表现花岗岩试样在受载荷过程中的基本特征。将数值模型的粒径范围设置为0.25~0.42 mm，颗粒半径服从均匀分布，孔隙率设置为15%，颗粒密度则设为3 150 kg/m³。黏结模型采用平行黏结模型（PB模型），反复修改并对PB模型中的线性组有效模量E^*、平行黏结组有效模量E˙*、法向切向刚度比k^*、黏结法向切向刚度比k˙*、颗粒尺寸、颗粒摩擦系数μ和平行黏结摩擦角φ˙等细观参数进行单轴压缩模拟试验。根据加载速度、模拟速度和实际模拟情况对参数进行微调，最终结果见表1。

根据标定参数模拟出的试样破坏状态以及应力—应变曲线与室内试验结果的对比如图1和图2所示。由图可知，模拟试验的应力—应变曲线与室内试验试样的应力—应变曲线吻合较好，破坏形式基本相似。

以往学者对节理非线性特性进行了诸多研究，通过建立不同的函数来描述节理法向闭合变形特性，如双曲线函数模型（Bandis et al.，1983）、指数函数模型（Malama et al.，2003）、改进的指数函数模型（Malama et al.，2003）、半对数函数模型（夏才初等，2002）和幂函数模型（Swan，1983）。其中，双曲线函数模型能够较好地描述节理闭合特性，其方程如下：

式中：σ为法向应力；d为节理闭合量；α和β均为常数。当σ →∞时，α/β= d_m，d_m为节理最大闭合量。

根据节理刚度K的定义有：

当σ = 0时，联立式（1）和式（2）可得节理初始刚度K₀为

则在任意应力下节理刚度表示为

在现有的PFC节理模型中，光滑节理模型表现出的节理闭合特性是线性的，不能描述节理非线性闭合情况，因此非线性变形节理闭合模型需采用PFC软件中的FISH语言进行二次开发实现。

在PFC模型中，颗粒属于刚体，在力的作用下颗粒会产生重叠，使得颗粒间接触力大小和接触位移的大小被改变。在细观层面假设节理间接触位移δ为细观节理闭合量，两颗粒中半径较小的颗粒半径R为细观节理最大闭合量，则细观节理法向刚度可表示为

式中：δ为细观节理闭合量；R为细观节理最大闭合量。

利用FISH语言实现节理非线性变形特性的程序流程如图3所示。首先判断颗粒与颗粒之间的接触类型是否为光滑节理模型，若不是光滑节理模型则保持其原有的状态不做修改，若接触模型为光滑节理模型，则再判定接触间距是否小于0，若接触间距小于0，则找到接触颗粒中较小颗粒的半径R，节理间接触位移δ为法向位移，自定义一个变量为perc=(1-δR)，再确定一个细观初始刚度K˙0，通过式（5）计算获得细观节理法向刚度K˙。

选择室内物理试验中分形维数最小的花岗岩节理试样，即节理粗糙度系数最小的试样，JRC值为0~2，使用CAD绘图软件将节理面描绘出来，用PFC命令流导入模型中，建立的岩石节理颗粒流模型如图4所示。

设置细观初始刚度为30 GPa/m，采用改进设置后的非线性变形节理闭合模型进行单轴压缩数值试验，根据标定参数模拟出的节理闭合变形曲线如图5所示。

由图5可知，改进设置后的节理闭合颗粒流模型能够较好地描述节理法向闭合曲线的非线性特性。

建立带有非线性变形特性节理的岩石模型，模型的尺寸为1 000 mm×50 mm（长×宽），颗粒密度为3 150 kg/m³，孔隙率为15%，颗粒半径在0.50~0.84 mm之间，颗粒服从均匀分布。将距离模型左右两边界宽度为3 mm范围内的颗粒设置为边界颗粒，目的在于方便模型动态载荷及边界条件的施加（石崇等，2018）。

本研究采用接触黏结模型，原因是该模型具有更简单的个体响应，从而产生更便于预测的宏观模型（Resende et al.，2010）。经过初步模拟对比分析，相比平行黏结模型，接触黏结模型能够更好地反映应力波传播的过程，且更易于标定参数。由于应力波加载过程属于动载荷加载，其中的接触黏结参数根据颗粒半径、模型尺寸、孔隙率和颗粒数目进行多次试错模拟和调整，对系统有效模量E^*进行标定（Zhou et al.，2017；王卫华等，2021），得到有效模量值为62 GPa。

为了让杆件能够完整地演示应力波传播过程不会因为应力波幅值过大使得模型被冲散而受到破坏，将切向和法向黏结强度数值设置为1×10¹⁰⁰ Pa，而后再添加节理面。首先在CAD程序中将节理面描绘出来，然后采用PFC命令流将其导入模型中，最后将节理面上下相邻的颗粒作为节理面上的颗粒，节理采用光滑节理模型（Mehranpour et al.，2017）。本文所建立的含非线性变形节理岩杆的颗粒流模型如图6所示。

节理面颗粒半径、颗粒密度、孔隙率和摩擦系数等基本参数与岩体颗粒性质保持一致，光滑节理模型中的黏聚力和节理摩擦角采用默认值，接触刚度将在后文进行选取分析。

应力波选用半正弦波，通过在左端施加动应力的方法加载入射应力波（石崇等，2018），右端设置为吸收边界条件（汤云坤，2021），其他各面设置为自由面。本次模拟研究采用控制变量法进行，通过分析节理颗粒间接触数目、节理接触力链分布及变化情况来研究刚度、波幅值和波频率对应力波透反射波幅值和透反射系数的影响。

将透射波峰值与入射波峰值的比值定义为透射系数模拟值，将反射波峰值与入射波峰值的比值定义为反射系数模拟值（Li et al.，2010），即：

式中：T、R表示透射系数和反射系数；V_I、V_T和V_R分别表示入射波、透射波和反射波质点速度。

非线性变形节理处应力波传播的透射系数T₀和反射系数R₀的理论解为式（8）和式（9）（王卫华等，2006）。

式中：K_e为节理的等效刚度；ω为波的圆频率；Z为介质的波阻抗，Z=ρc。

将幅值为5 m/s、频率为10 kHz的半正弦波垂直入射加载到含有相同节理的岩杆中，在距离左边界0.1 m和0.9 m处用测量圆监测入反射波和透射波，岩石密度为2 700 kg/m³，波速为4 191 m/s，后续模拟均采用该波速和岩石密度。

在节理处应力波传播的颗粒流模拟中，颗粒模型的细观初始刚度是最重要的细观参数之一。通过模拟分析可得到不同细观初始刚度下透反射系数与模拟等效刚度值，根据透反射系数理论，由式（8）和式（9）可反算出理论等效刚度，得到的模拟等效刚度和理论等效刚度如表2所示。

由表2可知，模拟等效刚度与理论等效刚度十分接近，证明本文第2节中所建立非线性变形节理颗粒流模型的正确性。

图7~图9分别是细观初始刚度为100，200，300 GPa/m时的速度波形图和应力波形图。从图中可以看出，反射和透射应力波波形均与入射波波形相似。随着细观初始刚度的增加，透射波的幅值呈现逐渐递增的趋势，反射波的峰值呈现逐渐递减的趋势。

不同等效刚度下透反射系数模拟结果与理论值对比如图10所示。由图10可知，透反射系数模拟结果与理论值具有很好的一致性，验证了本文模拟计算的正确性。由图10可以看出，透射系数随节理等效刚度的增加而增大，反射系数则与之相反，随节理等效刚度的增加而减小。节理等效刚度在0~800 GPa/m时，透射系数增长迅速；当节理等效刚度大于800 GPa/m之后，透射系数增长缓慢并趋于平缓。

不同幅值应力波在节理处传播时的透反射波速度波形图和应力波形图如图11和图12所示。

当入射波幅值增大时，节理闭合量随之增大，闭合量越大则节理刚度越大，相应的透射系数也越大。从图11和图12中可以看出，振幅为1 m/s和5 m/s的应力波波形与入射波波形相似，均保持良好的半正弦波形态。另外，观察速度波形图可以发现，振幅越大的入射应力波在节理处的透射波幅值越大，相应的透射系数也越大。这是因为当入射波幅值增大时，节理闭合量随之增大，闭合量越大则节理刚度越大，因而透射系数也越大。

图13所示为透反射系数随应力波幅值的变化情况。由图13可以看出，当应力波幅值在1~10 m/s范围内时，随着应力波幅值的增加，透射系数增大，反射系数减小，透射系数随振幅的增大呈近似线性增加。

图14所示为节理面颗粒间接触力链随时间的变化情况。其中，灰色颗粒为节理面两侧的颗粒，绿色颗粒为节理面颗粒；颗粒间的红色短线表示力链，力链的粗细表示力的大小，力链越粗表示接触力越大。在0.13 ms时刻，应力波未到达节理面，节理面颗粒间的接触保持原状，几乎没有法向位移，接触力接近于0；0.16 ms时刻应力波开始穿过节理面，节理面颗粒间的接触数目增多，接触力增大，法向位移增大，力链在整个节理面都有分布，部分应力波开始产生透反射现象；0.20 ms时刻应力波波峰穿过节理面，颗粒接触数目达到最大，节理面部分颗粒间的法向位移达到最大，整体节理面颗粒间法向位移增加，力链在整个节理面保持均匀分布，使应力波持续均匀穿过节理面；0.23 ms时刻颗粒接触数目明显下降，法向位移变小，力链整体变细，节理面间传递的接触力变小；随着时间继续增加，颗粒间接触数目和法向位移进一步减小，力链极细，节理面两侧颗粒逐渐分离，0.26 ms时刻应力波全部穿过节理面。

图15所示为应力波穿过节理的过程中节理面颗粒接触数目的变化情况。以波幅值为5 m/s为例进行分析，在0.13 ms时应力波波头开始到达节理面，此时颗粒接触数目约为49个，在0.16 ms时刻应力波到达节理面，随后节理面两边的颗粒在应力波的作用下发生运动产生法向位移，节理面颗粒接触数目也逐渐增多，随着波峰穿过界面，0.23 ms以后颗粒接触数目逐渐下降，在0.26 ms时刻，应力波已离开节理面，此时节理面接触数目与初始时刻的接触数目接近，之后由于节理面速度不为0，节理面持续运动，节理面的法向位移继续增大，在某一时刻导致大量节理颗粒分离，接触数目骤减，接触数目曲线呈断崖式下滑。从振幅为1，3，5，7 m/s时不同接触数目曲线可以看出，颗粒接触数目的最大值随着振幅的增大而增大，但当振幅达到某一个临界值后，接触数目的最大值不再增加，这是由节理面模型生成时的颗粒数目所决定的。从细观角度来看，节理面颗粒接触数目越多则颗粒间接触力链形成越紧密，接触数目最大值持续时间也越长，应力波越容易穿过节理面传播。

不同频率应力波在节理处传播时的透反射波速度波形图和应力波形图如图16~图18所示。

由图16~图18可知，透射波幅值随着频率的增加而减小，透射系数减小，反射波幅值随着频率的增加而增大，反射系数增大。当频率为20 kHz时，透反射波均出现振荡，这种振荡随着波频率的增加而表现得越剧烈。

图19所示为透射系数随着应力波频率的变化曲线。由图19可知，整体上应力波透射系数随入射波频率的增大而减小，这也印证了节理的高频滤波作用。

图20所示为不同频率应力波穿过节理时节理面颗粒接触数目的变化曲线。

由图20可知，在不同频率应力波作用下，节理面颗粒接触数目的最大值是相同的，这是由节理面建模时颗粒数目决定的。另外，随着频率的增加，颗粒接触数目达到最大值的时间越早，持续时间也越短，产生断崖式下跌的时间会更早到来，这意味着节理面颗粒间分离得越快，应力波穿越节理传播越困难。

利用PFC软件中FISH语言对光滑节理模型进行二次开发，建立了非线性变形节理岩体颗粒流模型，并模拟了应力波在非线性变形节理处的传播，分析了节理刚度、波幅值和波频率对透反射系数的影响，还从接触力链和节理间颗粒接触数目的细观角度分析应力波在非线性变形节理处的传播过程。通过研究，得到以下3点认识：

（1）随着节理刚度的增加，透射系数增大，透射系数越大，反射系数则越小，当节理刚度达到某一阈值后透射系数增加不明显。

（2）振幅对应力波在非线性变形节理处传播有很大影响，入射波的振幅越大，透射系数越大，反射系数则越小，细观层面上表现为节理面颗粒接触持续时间变长，导致应力波更容易穿过节理面。

（3）应力波透射系数随入射波频率的增大而减小。低频率范围内，随着频率的增加透射系数急剧下降，高频率范围内，随着频率的增加透射系数缓慢降低，节理表现出高频滤波的特性，高频率波比低频率波透射系数更小。