随着经济建设的飞速发展，地下资源和空间的开发利用规模日益扩大^[1-2]，大量地下岩体被开挖掘进，而钻爆法是目前最常用的施工方法^[3]。炸药在爆炸破碎岩体的同时，也会导致地下工程预保留围岩的动态扰动、损伤甚至破裂^[4-5]。循环爆破后，围岩将受到不同程度的累积动态损伤^[6-7]，导致其在承受地应力重分布带来的静态荷载时断裂韧度降低、承载能力下降，引发各种工程事故。因此，研究循环冲击损伤后岩石的静态断裂特征具有较强的工程意义。

对脆性岩石在循环冲击作用下的疲劳特性已有大量研究^[8-10]，但现有成果主要集中在循环冲击下岩石损伤、破裂直至完全破坏的全过程。如Li等^[11]采用分离式霍普金森压杆（split Hopkinson pressure bar，SHPB）对花岗岩进行了循环冲击加载实验，结果表明，当冲击气压在一定范围内时，每次冲击都会对岩石造成动态损伤但不至于完全破坏；林大能等^[12]研究了循环冲击作用下围压、冲击气压和冲击次数对岩石动态损伤的影响；王彤等^[13]利用动静组合加载装置进行不同轴压、不同冲击气压下的循环冲击实验。而在受损岩石的静态断裂特性研究方面，损伤诱因也多以环境因素为主，如左建平等^[14]探究了不同温度影响后花岗岩的细观断裂机制；贺晶晶等^[15]分析了冻融损伤对砂岩断裂性能的劣化影响，并对试样断裂破坏面的形貌特征进行了扫描分析；杨健锋等^[16]研究了不同程度水损伤作用对泥岩断裂力学特性的影响。

低能量密度的循环冲击能导致岩石的累积动态损伤乃至破裂，但是并不能使其完全破坏，此时岩石仍有一定的承载能力^[17]，但目前涉及循环冲击损伤后岩石静态承载能力的研究鲜有报道。虽然付安琪等^[18]利用SHPB系统对中心直切槽半圆盘（notched semi-circular bend，NSCB）大理岩试样进行了循环冲击损伤处理，然后对其进行静态断裂实验，分析了动态损伤对大理岩断裂力学性能的劣化影响，但是岩石材料离散性较大，实验结论适用范围有限。此外，由于室内实验条件的局限，该文献中未对动态损伤累积过程及静态断裂力学行为的劣化机制进行深入研究。

随着现代计算技术的不断发展，多种数值分析方法被用来构建SHPB模型并进行动态冲击实验，如有限元软件ABAQUS^[19]、AUTODYN^[20]及离散元软件PFC^[21-22]等。由于杆件由高强度钢材加工而成，而岩石试样是强度相对较低的脆性介质，在冲击实验过程中，杆件仅发生形变，而岩样则将发生破裂甚至破碎。因此，利用连续介质模型模拟杆件系统、离散介质模拟岩石试样是较为理想的数值建模方法。此外，仅采用三维离散元方法模拟SHPB实验时，只能通过提高杆件颗粒接触黏结强度来近似满足“应力均匀性假定”^[23]，并且杆件模型颗粒间的孔隙率对模拟结果的准确性也有一定影响。而连续-非连续耦合的建模方法则可以有效弥补这些缺陷，并能显著提高计算效率，但目前这种耦合技术很少被应用于SHPB实验的模拟。

鉴于此，本文中利用Itasca公司开发的FLAC^3D（基于有限差分方法FDM）和PFC^3D（基于离散元方法DEM）系列程序构建三维SHPB耦合模型。利用该模型对模拟NSCB试样沿厚度方向进行不同次数的循环冲击预损伤实验，随后对受损试样进行静态断裂韧度模拟实验，验证相关室内实验研究的结论，并结合微裂纹场、力链场等机理信息对动态损伤累积过程及断裂力学性能劣化机制进行深入研究。

1.   模　型

1.1   循环冲击实验

三维FDM-DEM耦合SHPB模拟系统如图1所示。其中杆件利用FLAC^3D中的线弹性模型建立，尺寸及材料参数与室内实验一致^[18]；NSCB岩样是基于PFC^3D构建的离散颗粒黏结体，并选择了能够重现加载过程中的微裂缝萌生、聚结和宏观裂缝形成的平行黏结模型（parallel bonding model，PBM），其尺寸构型与室内实验一致^[18]。与大理岩三点弯曲室内实验结果反复对比反馈，通过“试错法”对模拟试样细观参数进行标定，最终获得一组能够反映脆性大理岩力学特性的细观参数。

图  1  SHPB模拟系统

Figure  1.  Simulation of the SHPB system

下载: 全尺寸图片 幻灯片

模拟应力波如图2所示，可以看出，模拟波形的入射波和反射波之和近似等于透射波，并且在峰值后一定时间内仍保持相等。可见，本次模拟可实现试样的动态受力平衡。与室内实验一致^[18]，所有试样分为6组，利用SHPB系统对各组试样沿厚度方向分别进行了0～5次冲击速度恒定的等能量冲击，以获得损伤程度不同的6组岩样。

图  2  模拟应力波

Figure  2.  Stress waves obtained from the numerical simulation

下载: 全尺寸图片 幻灯片

1.2   静态三点弯曲断裂实验

预损伤实验完成后，将子弹与杆件移除。在受损试样的底部加设两根支撑钢棒，并在试样上侧以恒定速率施加荷载^[18]，如图3所示。图4中给出了完整试样静态三点弯曲实验中荷载-位移曲线的室内实验^[18]与数值模拟结果对比。可见模拟获得的试样峰值荷载及破坏位移均与实验结果一致。

图  3  静态三点弯曲模拟实验

Figure  3.  Simulation of the static three-point bending test

下载: 全尺寸图片 幻灯片

图  4  天然大理岩试样实验与模拟结果比较^[18]

Figure  4.  Comparison of the experimental and numerical results of the static three-point bending test of a natural marble sample^[18]

下载: 全尺寸图片 幻灯片

2.   循环冲击结果分析

2.1   应力波信号

图5中给出了首次冲击时杆件中的应力波传播过程。当t=0 μs时，子弹撞击杆件，产生的压缩波由入射杆左端进入SHPB系统。当t=409 μs时，压缩波到达试样的左端，随即产生拉伸波（即反射波）向入射杆左端方向传播。剩余应力波到达试样右端时，再次发生反射-透射，透射波作为压缩波向右端传播。需要注意的是，在反射拉伸波之后，有一个压缩波沿着入射杆向左端传播（图5中黑色线圈内）。这是因为在相对较低的冲击速度下，试件在加载期间（入射波上升段）未被完全破坏。在卸载期间（入射波下降段），试样中储存的部分应变能将沿着入射杆释放，即发生回弹现象，从而产生压缩波，这与Li等^[23]的结论一致。

图  5  应力波在杆件中的传播过程

Figure  5.  Stress wave propagation in the bars

下载: 全尺寸图片 幻灯片

图6为循环冲击过程中的应力波信号叠加曲线。从图6中可以看出，入射波幅值在每个冲击周期中基本重合，说明模拟实现了等幅循环加载。随着循环冲击次数n的增加，透射波的幅值越来越小，而反射波的幅值则越来越大。主要原因为，在循环冲击作用下，岩石内部损伤累积，岩石孔隙率增加，而波阻抗降低^[24]。与应力传播过程一致，在卸载段，入射杆上采集到的应力波信号有一段为负值，对应上述回弹现象。

图  6  循环冲击加载过程中应力波信号变化

Figure  6.  Variation of the stress wave signals during the cyclic impact loading

下载: 全尺寸图片 幻灯片

2.2   动态应力-应变曲线

图7为模拟试样在循环冲击过程中的动态应力-应变曲线，可划分为3个阶段：弹性变形阶段、裂纹扩展阶段和应力卸载阶段^[24]。在弹性变形阶段，轴向应变几乎随动态应力的增加而线性增加。随着应力的进一步增大，试样进入微裂纹扩展阶段，微裂纹萌生、扩展并相互作用，导致曲线呈现非线性行为。应力卸载过程对应于曲线的峰后段。

图  7  模拟动态应力-应变曲线

Figure  7.  Dynamic stress-strain curves obtained from numerical simulations

下载: 全尺寸图片 幻灯片

图8为峰值应力的实验值与模拟值变化情况。当冲击次数由1增加至5时，峰值应力模拟值从75.1 MPa降低到64.9 MPa，降幅13.5%，远小于实验值的30.4%，这主要是因为天然试样内部存在大量缺陷，这些缺陷会在动态冲击过程中进一步发育，加剧试样荷载传递能力的劣化。另一方面，随着冲击次数n的增加，峰值应变增大，裂纹扩展，同时大量新裂纹萌生，卸载阶段试样恢复的应变能降低。

图  8  峰值应力实验值与模拟值比较^[18]

Figure  8.  Comparison of the dynamic peak stress results obtained from experiments and numerical simulations^[18]

下载: 全尺寸图片 幻灯片

2.3   动态损伤演变

PFC^3D程序中，试样颗粒之间的黏结断裂即视为产生了微裂纹，微裂纹扩展严重区域常伴随着颗粒脱落，形成碎块。图9中给出了循环冲击过程中模拟试样裂纹场及碎块场的演变过程。从图9中可以直观了解到循环冲击过程中试样内部动态损伤演变过程：随着冲击次数的增加，试样内部微裂纹及碎块数目均有明显上升。图9（b）红色圈内为撞击掉落颗粒，与室内实验结果（图9（c））中观察到的白斑及边缘掉落颗粒一致，为岩石颗粒在动载作用下相互错动导致^[25]。

图  9  动态损伤演变过程

Figure  9.  Evolution process of the dynamic damage

下载: 全尺寸图片 幻灯片

为定量分析，图10中给出了微裂纹数量N_c及碎块数量N_f随加载次数n的变化规律。可见试样损伤演变过程可以划分为急剧增长阶段、缓慢发展阶段和急剧增长阶段：当n = 1时，首次冲击造成的损伤较大，N_c及N_f分别由0增长至3558和32；当n = 2～4时，单次冲击造成的损伤较小，冲击4次后，N_c及N_f分别为6402和58；当n = 5时，单次冲击造成的损伤较大，N_c及N_f分别增长至11331和198。

图  10  微裂纹及碎块数量随冲击次数的变化规律

Figure  10.  Variation of the microcrack number and fragment number with the cyclic impact number

下载: 全尺寸图片 幻灯片

模拟试样颗粒之间具有黏结力，形成的力链在受到荷载作用时会改变，甚至断裂。因此，力链演变能够揭示受损试样的动态损伤累积机理^[26]。图11中给出了每次冲击后试样接触力场（剖面图）的演变过程。

图  11  循环冲击作用下试样接触力场演变过程

Figure  11.  Evolution of the contact force field of a sample under cyclic impact loading

下载: 全尺寸图片 幻灯片

由图11可知，在动态荷载冲击作用下，试样内部力链分布均匀性变差，部分力链发生断裂。当n = 5时，试样力链出现区域性缺失（黑色虚线圈内）。此外，试样部分区域出现明显的应力集中（红色虚线圈内）现象，且均与室内实验中试样白斑及颗粒掉落位置一致。力链的断裂及缺失使试样在受到后续荷载作用时更易被破坏，这是受损试样力学特征劣化的根本原因。

穿晶裂纹和沿晶裂纹的数量可以在一定程度上反映试样结构的损伤程度。为了验证上述微观结构变化，图12中给出了受损试样内部结构偏光显微结果。可以看出，天然状态下矿物颗粒完整，没有明显微裂纹。然而，经过动态冲击后，颗粒之间黏结强度降低。当n = 1时，可以观察到明显的裂纹，说明岩石内部结构已被破坏；当n = 2～3时，裂纹数目增多，且周围晶体颗粒的破碎程度明显增加；当n = 4～5时，裂纹已发生合并，发展为宏观裂纹，周围的晶体颗粒几乎完全破碎，大理岩试样内部结构的动态损伤已经非常严重。

图  12  受损大理岩的微观图像

Figure  12.  Photomicrographs of the damaged marbles

下载: 全尺寸图片 幻灯片

3.   三点弯曲结果分析

3.1   荷载-位移曲线

图13中给出了受损试样静态荷载-位移曲线的模拟结果及实验结果。模拟曲线包含：（1）线弹性阶段，荷载随位移线性增长；（2）脆性破坏阶段，试样发生脆性断裂，荷载瞬间跌落。随着循环冲击次数的增加，试样内部损伤加剧，弹性段斜率随着冲击次数的增加而降低。模拟试样曲线整体趋势与室内实验结果一致。

图  13  模拟静态荷载-位移曲线对比^[18]

Figure  13.  Comparison of the static load-displacement curves obtained from experiments and numerical simulations^[18]

下载: 全尺寸图片 幻灯片

3.2   宏观力学特性

岩石静态Ⅰ型断裂韧度K_ⅠC的计算方法如下^[27]：

式中：p_max为试样破坏时的峰值荷载，R为试样半径，B为试样厚度，Y为无量纲应力强度因子，与预制裂缝长度a及试样支撑间距S有关，计算公式如下：

式中：α为无量纲预制裂缝长度，α = a/R；S/2R是无量纲支撑间距。本次模拟中α = 0.2，S/2R = 0.7。受损试样静态断裂韧度及破坏位移随循环冲击次数n的变化规律如图14所示。

图  14  断裂参数实验值和模拟值比较^[18]

Figure  14.  Comparison of the fracture parameters obtained from experiments and numerical simulations^[18]

下载: 全尺寸图片 幻灯片

当循环冲击次数n从0增加至5时，断裂韧度由1.28 MPa·m^0.5减小至 1.07 MPa·m^0.5，破坏位移由 0.21 mm增长至0.27 mm。随着n的增加，断裂韧度模拟值及实验值均呈下降趋势。破坏位移模拟值及实验值变化过程基本一致。值得注意的是，模拟试样参量的变化程度远小于室内实验中的天然试样，这主要是由于模拟试样内部不存在天然缺陷，而这些缺陷会在动态冲击中进一步扩展，汇聚，进而对后续静力学参数变化产生影响。

3.3   试样破坏特征

图15中给出了模拟试样断裂破坏后的裂纹场及碎块场演变。从图5中可以看到，受损试样的破坏形式均为典型的三点弯曲断裂破坏：从预制裂缝的顶端处起裂，扩展至加载点附近。

图  15  模拟试样断裂破坏形态

Figure  15.  Fracture and failure patterns of the simulated sample

下载: 全尺寸图片 幻灯片

图16中给出了静态加载过程中新增微裂纹及碎块数量随循环冲击次数n的变化规律。随着冲击次数的增加，受损试样力链断裂部分越来越多，导致静态加载过程中更多的微裂纹萌生并扩展。n由0增加至5时，新增微裂纹及碎块数量分别由1429、11增长至1887、34。当外部荷载作用于试样、内部形成的集中应力大于试样损伤阈值时，接触断裂，微裂纹产生。而当外部荷载、参与试样变形的总接触数量一定时，接触断裂现象越明显，试样越易发生宏观破裂。结合图13(a)所示，与天然试样相比，当n = 1～5时，试样刚度下降幅度分别为27.23%、36.05%、39.38%、55.53%、76.59%。因此，试样的抗变形能力随着冲击次数的增加而逐渐减小。

图  16  新增微裂纹及碎块数量随冲击次数的变化情况

Figure  16.  Variation of the new microcrack and fragment number with the cyclic impact number

下载: 全尺寸图片 幻灯片

利用surfer软件对破坏后的试样断裂面进行后处理，得到断裂面重建图，如图17所示。为定量表征断裂面粗糙度，采用断裂面轮廓高度均值（S_a）、断裂面轮廓最大最小高度差值（S_z）对断裂面进行精确描述。

图  17  模拟试样断裂面形貌

Figure  17.  Topography of the fracture surface of the numerical sample

下载: 全尺寸图片 幻灯片

断裂面轮廓高度均值：

断裂面轮廓最大最小高度差：

式中：N为裂隙面测点数；z_i为裂纹面内i点处的高度；z_max、z_min分别表示裂纹面内最高点、最低点处的高度。

由图18可知，断裂面粗糙度随冲击次数的增加而增加。与天然状态相比，循环冲击5次受损试样断裂面的S_a、S_z分别增长了32.06%、27.70%。在动态荷载作用下，微裂纹在试样内部随机产生，介质非均匀性、非连续性程度升高，致使试样破裂面粗糙度上升，与室内实验结论一致^[18]。

图  18  断裂面粗糙度变化

Figure  18.  Variation of the fracture surface roughness

下载: 全尺寸图片 幻灯片

4.   结　论

基于有限差分法FDM及离散元方法DEM耦合思想，利用FLAC^3D及PFC^3D软件构建了三维SHPB模型。利用该模型对NSCB试样进行沿厚度方向的等能量循环冲击，对受损试样进行静态三点弯曲实验，探究了大理岩试样在循环冲击作用下动态损伤累积机理及受损试样的静态断裂力学特征劣化机制，主要结论如下。

（1）压缩作用下的应力波信号与实验结果较吻合，从而验证了FEM-DEM耦合SHPB系统用于动力加载的可行性。此外，微裂纹场、碎块场、接触力链结构、峰值强度、断裂韧度等特征均与实验结果吻合。

（2）在较低动态荷载作用下， 应力应变曲线会因试样在卸载段释放应变能而产生回弹。随着循环冲击次数的增加，试样虽未发生整体破碎，但内部微裂隙、破碎颗粒、断裂力链数量均增加，其动态峰值应力降低，动态损伤不断累积。

（3）受损试样的静态断裂韧度较天然试样明显劣化，破坏应变则呈上升趋势。随着冲击次数的增加，试样力链断裂现象更加明显。试样介质非均匀性、非连续性程度在循环冲击过程中不断升高，致使试样破裂面粗糙度上升、断裂面粗糙度增加。