摘要：为探究循环冲击损伤后大理岩的静态断裂力学特征，基于有限差分(FDM)-离散元(DEM)耦合的建模技术构建了三维分离式霍普金森压杆(SHPB)数值模型，其中杆件系统和岩石试件分别采用FLAC3D和PFC3D程序建模。利用该模型对中心直切槽半圆盘(NSCB)试样进行了恒定子弹速度下的循环冲击，随后对受损试样进行静态三点弯曲断裂试验。通过编写Fish程序，提取试样断裂面数据，对断裂面进行重构并定量计算表面粗糙度。通过与相关室内试验结果的对比分析，验证了本文数值分析的合理性与可靠性。模拟结果表明，随着循环冲击次数的增加，试样内部微裂纹、破碎颗粒均增加。连接力场分布混乱，部分力链发生断裂。力链的变化是试样力学性能劣化的根本原因。在静态三点弯曲断裂试验中，冲击5次后的试样的静态断裂韧度较天然试样降低53.35%。试样在静载过程中会产生微裂纹和碎块的数量随着循环冲击次数的增加而增加。断裂面粗糙度随循环冲击次数的增加而增加。研究结论可为工程实践提供一定指导意义。
关键词：岩石力学;有限差分-离散元耦合;霍普金森;循环冲击;断裂韧度
随着我国经济建设的飞速发展，地下资源和空间的开发利用规模日益增大[1,2]，势必带来大量的地下岩体开挖掘进工作，而钻爆法在我国仍是目前最常用的施工方法[3]。炸药爆炸破碎岩体的同时，也会导致地下工程预保留围岩的动态扰动、损伤甚至破裂[4,5]。循环爆破后，围岩将受到不同程度的累积动态损伤[6,7]，导致其在承受地应力重分布带来的静态荷载时断裂韧度下降、承载能力大幅劣化，极易引发各种工程事故。因此，循环冲击损伤后岩石的静态断裂特征研究具有较强的工程意义。
学者们对脆性岩石在循环冲击作用下的疲劳特性进行了大量研究[8-10]，但现有成果主要集中在循环冲击下岩石损伤、破裂直至完全破坏全过程。如Li等[11]采用霍普金森压力杆系统(splitHopkinsonpressurebar,SHPB)对花岗岩进行了循环冲击加载试验，结果表明，当冲击气压在一定范围内时，每次冲击都会对岩石造成动态损伤但不至于完全破坏;林大能等[12]研究了循环冲击作用下围压、冲击气压和冲击次数对岩石动态损伤的影响;王彤等[13]利用动静组合加载装置进行不同轴压、不同冲击气压下的循环冲击试验。而在受损岩石的静态断裂特性研究方面，损伤诱因也多以环境因素为主，如左建平等[14]探究了不同温度影响后花岗岩的细观断裂机制;贺晶晶等[15]分析冻融损伤对砂岩断裂性能的劣化影响，并对试样断裂破坏面的形貌特征进行了扫描分析;杨建峰等[16]研究了不同程度水损伤作用对泥岩断裂力学特性的影响。
低能量密度的循环冲击能导致岩石的累积动态损伤乃至破裂，但是并不能致使其完全破坏，此时岩石仍有一定的承载能力[17]。但是，目前涉及循环冲击损伤后岩石静态承载能力的研究鲜有报道。虽然付安琪等[18]利用SHPB系统对中心直切槽半圆盘(notchedsemi-circularbend，NSCB)大理岩试样进行了循环冲击损伤处理，然后对其进行静态断裂试验，分析了动态损伤对大理岩断裂力学性能的劣化影响。但是岩石材料离散性较大，试验结论适用范围有限。此外，由于室内试验研究条件的局限，该论文未能对动态损伤累积过程及静态断裂力学行为劣化机制进行深入研究。
随着现代计算技术的不断发展，多种数值分析方法被用来构建SHPB模型并进行动态冲击试验，如有限元软件ABAQUS[19]、AUTODYN[20]，及离散元软件PFC[21,22]等。由于杆件由高强度钢材加工而成，而岩石试样是强度相对较低的脆性介质，在冲击试验过程中，杆件仅发生形变，而岩样则将发生破裂甚至破碎。因此，杆件系统利用连续介质模型模拟，同时岩石试样利用离散介质模拟是较为理想的数值建模方法。此外，仅采用三维离散元方法模拟SHPB试验时，只能通过提高杆件颗粒接触粘结强度来近似满足“应力均匀性假定”[23]，并且杆件模型颗粒间的孔隙率对模拟结果准确性也有一定影响。而连续-非连续耦合的建模方法则可以上述缺陷进行有效的改善，并能显著提高计算效率，但目前这种耦合技术很少被应用于SHPB试验的模拟。
鉴于此，本文利用Itasca公司开发的FLAC3D(基于有限差分方法FDM)和PFC3D(基于离散元方法DEM)系列程序构建了三维SHPB耦合模型。利用该模型对模拟NSCB试样沿厚度方向进行不同次数的循环冲击预损伤试验，随后对受损试样进行静态断裂韧度模拟试验，验证了相关室内试验研究的结论，并结合微裂纹场、力链场等机理信息对动态损伤累积过程及断裂力学性能劣化机制进行深入研究。
1模型描述
1.1循环冲击试验
三维FDM-DEM耦合SHPB模拟系统如图1所示。其中杆件利用FLAC3D中的线弹性模型建立，尺寸及材料参数与室内试验一致[18];NSCB岩样是基于PFC3D构建的离散颗粒粘结体，并选择了能够重现加载过程中的微裂缝萌生、聚结和宏观裂缝形成的平行黏结模型(parallelbondingmodel,PBM)，其尺寸构型与室内试验一致[18]。与大理岩三点弯曲室内试验结果反复对比反馈，通过“试错法”对模拟试样细观参数进行标定，最终获得一组能够反映脆性大理岩力学特性的细观参数，见表1。
应力波对比如图2所示，可以看出模拟波形与实验波形有很好的一致性。入射波和反射波之和近似等于透射波，并且在峰值后一定时间内仍保持相等。由此可见，本次模拟可实现试样的动态受力平衡。与室内试验一致[18]，所有试样分为6组，利用SHPB系统对各组试样沿厚度方向分别进行了0~5次冲击速度v均为4m/s的等能量冲击，以获得损伤程度不同的6组岩样。
1.2静态三点弯曲断裂试验
预损伤试验完成后，将子弹与杆件移除。在受损试样的底部加设两根支撑钢棒，并在试样上侧以0.06mm/min的恒定速率施加荷载[18]，如图3所示。图4给出了完整试样静态三点弯曲试验中荷载-位移曲线的室内试验[18]和数值模拟结果对比。可见模拟获得的试样峰值荷载及破坏位移均与试验结果一致，即表1所示细观参数较为合理。
2循环冲击结果分析
2.1应力波信号
图5给出了首次冲击时杆件中的应力波传播过程。t=0μs时，子弹撞击杆件，产生的压缩波由入射杆左端进入SHPB系统。t=409μs时，压缩波到达试样的左端，随即产生拉伸波(即反射波)向入射杆左端方向传播。剩余应力波到达试样右端时，再次发生反射-透射，透射波作为压缩波向右端传播。需要注意的是，在反射拉伸波之后，有一个压缩波沿着入射杆向左端传播(图5黑色线框内)。这是因为在相对较小的冲击速度下(4m/s)，试件在加载期间(入射波上升段)未被完全破坏。在卸载期间(入射波下降段)，试样中储存的部分应变能将沿着入射杆释放，即发生回弹现象，从而产生压缩波，这与Li[23]的结论一致。
图6为循环冲击过程中的应力波信号叠加图。从图中可以看出，入射波幅值在每个冲击周期中基本重合，说明模拟实现了等幅循环加载。随着循环次数的增加，透射波的幅值越来越小，而反射波则越来越大。主要原因为在循环冲击作用下，岩石内部损伤累积，岩石孔隙率增加，而波阻抗降低[24]。与应力传播过程图一致，在卸载段，入射杆上采集到的应力波信号有一段为负值，对应上述回弹现象。
2.2动态应力-应变曲线
图7为模拟试样在循环冲击过程中的动态应力-应变曲线，可划分为三个阶段：弹性变形阶段、裂纹扩展阶段和应力卸载阶段[24]。在弹性变形阶段，轴向应变几乎随动态应力的增加而线性增加。随着应力的进一步增大，试样进入微裂纹扩展阶段，微裂纹萌生、扩展并相互作用，导致曲线呈现非线性行为。应力卸载过程对应于曲线的峰后段。
图8为峰值应力试验值与模拟值变化对比。当冲击次数由1增长至5时，峰值应力模拟值从105.45MPa降低到75.76MPa，降幅为28.16%，试验值由103.94MPa降低到72.30MPa，降幅为30.44%。另一方面，随着冲击次数n的增大，峰值应变增大。随着冲击次数的增加，原生裂纹扩展，同时大量的新裂纹萌生，卸载阶段试样恢复的应变能降低。
2.3动态损伤演变
PFC3D程序中，试样颗粒之间的粘结断裂即视为产生了微裂纹，微裂纹扩展严重区域常伴随着颗粒脱落，形成碎块。图9给出了循环冲击过程中模拟试样裂纹场及碎块场的演变过程。从图9中可以直观了解到循环冲击过程中试样内部动态损伤演变过程：随着冲击次数的增多，试样内部微裂纹及碎块数目均有明显上升。图9b红色框内为撞击掉落颗粒，与室内试验结果(图9c)中观察到的白斑及边缘掉落颗粒一致，为岩石颗粒在动载作用下相互错动导致[25]。
为定量分析，图10给出了微裂纹数量Nc及碎块数量Nf随加载次数n的变化规律。可见试样损伤演变过程可以划分为急剧增长阶段，缓慢发展阶段和急剧增长阶段：当n=1时，首次冲击造成的损伤较大，Nc及Nf分别由0增长至3558和32;当n=2-4时，单次造成的损伤较小，冲击4次后的Nc及Nf分别为6402和58;当n=5时，单次造成的损伤较大，Nc及Nf分别增长至11331和198。
模拟试样颗粒之间具有粘结力，形成的力链在受到荷载作用时会改变，甚至断裂。因此，力链演变能够揭示受损试样的动态损伤累积机理[26]。图11给出了每次冲击后试样接触力场(剖面图)演变过程。
由图可见，在动态荷载冲击作用下，试样内部力链分布均匀性变差，部分力链发生断裂。当n=5时，试样力链出现区域性缺失(黑色虚线框内)。此外，试样部分区域出现明显的应力集中(红色虚线框内)，且均与室内试验中试样白斑及颗粒掉落位置一致。力链的断裂及缺失使试样在受到后续荷载作用时更易被破坏，这是受损试样力学特征劣化的根本原因。
为了验证上述微观结构变化，图12给出了受损试样内部结构偏光显微结果。穿晶裂纹和沿晶裂纹的数量可以在一定程度上反映试样结构的损伤程度。可以看出，天然状态下矿物颗粒完整，没有明显微裂纹。然而，经过动态冲击后，颗粒之间粘结强度降低。n=1时，可以观察到明显裂纹，说明岩石内部结构已被破坏。当n=2、3时，裂纹数目增多，且周围晶体颗粒的破碎程度明显增加。当n=4、5时，裂纹已发生合并，发展为宏观裂纹，周围的晶体颗粒几乎完全破碎，大理岩试样内部结构的动态损伤已经非常严重。
3三点弯曲结果分析
3.1荷载–位移曲线
图13给出了受损试样静态荷载 位移曲线的模拟结果及试验结果。模拟曲线包含：(1)线弹性阶段：荷载随位移线性增长;(2)脆性破坏阶段：试样发生脆性断裂，荷载瞬间跌落。随着循环冲击次数的增长，试样内部损伤程度加剧，裂纹数目增长，曲线压密段持续时间相应增长，而弹性段斜率随着冲击次数的增长而降低。模拟试样曲线整体趋势与室内试验结果一致。
3.3试样破坏特征
图15给出了模拟试样断裂破坏后的裂纹场及碎块场演变。可以看到，受损试样的破坏形式均为典型三点弯曲断裂破坏：从预制裂缝的顶端处起裂，扩展至加载点附近。
图16给出了静态加载过程中新增微裂纹及碎块数量随循环冲击次数n的变化规律。随冲击次数的增加，受损试样力链断裂部分越来越多，导致静态加载过程中更多的微裂纹萌生并扩展。n由0增长至5时，新增微裂纹及碎块数量分别由1429、11增长至1887、34。当外部荷载作用于试样，其内部所形成的集中应力大于试样损伤阈值时，接触断裂，微裂纹产生。而当外部荷载、参与试样变形的总接触数量一定时，接触断裂现象越明显，说明试样将越易发生宏观破裂。结合图13(a)所示，较天然试样，n=1~5时，试样刚度下降幅度分别为27.23%,36.05%,39.38%,55.53%,76.59%。因此，试样的抗变形能力随着冲击次数的增加而逐渐减小。
对破坏后的试样断裂面利用surfer软件进行后处理，得到断裂面重建图，如图17所示。为定量表征断裂面粗糙度，本次研究采用断裂面轮廓高度均值(Sa)，断裂面轮廓最大最小高度差值(Sz)对断裂面进行精确化描述。
由图18可知，断裂面粗糙度随冲击次数的增加而增加。较天然状态，循环冲击5次受损试样断裂面的Sa，Sz分别增长了32.06%，27.70%。在动态荷载作用下，微裂纹在试样内部随机产生，介质非均匀性、非连续性程度升高，致使试样破裂面粗糙度上升，与室内试验结论一致[18]。
4结论
本文基于有限差分法FDM及离散元方法DEM耦合思想，利用FLAC3D及PFC3D软件构建了三维SHPB模型。利用该模型对NSCB试样进行沿厚度方向的等能量循环冲击，随后对受损试样进行静态三点弯曲试验，探究了大理岩试样在循环冲击作用下动态损伤累积机理及受损试样的静态断裂力学特征劣化机制。主要结论如下：
(1)压缩作用下的应力波信号与试验结果较吻合，从而验证了FEM-DEM耦合SHPB系统用于动力加载的可行性。此外，微裂纹场、碎块场、接触力链结构、峰值强度、断裂韧度等特征均与试验结果吻合。
(2)在较低动态荷载作用下，应力-应变曲线会因试样在卸载段释放应变能而产生回弹。随着循环冲击次数的增加，试样虽未发生整体破碎，但内部微裂隙、破碎颗粒、断裂力链数量均增加，其动态峰值应力降低，动态损伤不断累积。
(3)受损试样的静态断裂韧度较天然试样明显劣化，破坏应变则呈上升趋势。随着冲击次数的增加，试样力链断裂现象愈发明显。试样介质非均匀性、非连续性程度在循环冲击过程中不断升高，致使试样破裂面粗糙度上升断裂面粗糙度增加。——论文作者：张涛1，蔚立元1，苏海健1，罗宁1，魏江波2
相关期刊推荐：《爆炸与冲击》杂志报道爆炸力学学科领域的国内外最新科技成果，反映学术前沿进展及水平，促进学术交流，创造本学科领域青年人才良好的成长环境，推进爆炸理论和应用、抗爆与爆炸安全技术的发展。主要读者对象为从事爆炸物理、岩土爆炸、工程和器材爆炸等方面的科技人、员工程设计人员和高等院校师生。
