基于FDM-DEM 耦合的冲击损伤大理岩静态断裂力学特征研究*

2022-02-11蔚立元苏海健魏江波

爆炸与冲击 2022年1期

张涛，蔚立元，苏海健，罗宁，魏江波

（1. 中国矿业大学深部岩土力学与地下工程国家重点实验室，江苏徐州 221116；2. 西安科技大学地质与环境学院，西安陕西 710054）

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

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

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

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

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

1 模型

1.1 循环冲击实验

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

图1 SHPB 模拟系统Fig. 1 Simulation of the SHPB system

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

图2 模拟应力波Fig. 2 Stress waves obtained from the numerical simulation

1.2 静态三点弯曲断裂实验

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

图3 静态三点弯曲模拟实验Fig. 3 Simulation of the static three-point bending test

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

图5 应力波在杆件中的传播过程Fig. 5 Stress wave propagation in the bars

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

图6 循环冲击加载过程中应力波信号变化Fig. 6 Variation of the stress wave signals during the cyclic impact loading

2.2 动态应力-应变曲线

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

图7 模拟动态应力-应变曲线Fig. 7 Dynamic stress-strain curves obtained from numerical simulations

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

图8 峰值应力实验值与模拟值比较[18]Fig. 8 Comparison of the dynamic peak stress results obtained from experiments and numerical simulations[18]

2.3 动态损伤演变

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

图9 动态损伤演变过程Fig. 9 Evolution process of the dynamic damage

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

图10 微裂纹及碎块数量随冲击次数的变化规律Fig. 10 Variation of the microcrack number and fragment number with the cyclic impact number

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

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

图11 循环冲击作用下试样接触力场演变过程Fig. 11 Evolution of the contact force field of a sample under cyclic impact loading

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

图12 受损大理岩的微观图像Fig. 12 Photomicrographs of the damaged marbles

3 三点弯曲结果分析

3.1 荷载-位移曲线

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

图13 模拟静态荷载-位移曲线对比[18]Fig. 13 Comparison of the static load-displacement curves obtained from experiments and numerical simulations[18]

3.2 宏观力学特性

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

图14 断裂参数实验值和模拟值比较[18]Fig. 14 Comparison of the fracture parameters obtained from experiments and numerical simulations[18]

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

3.3 试样破坏特征

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

图15 模拟试样断裂破坏形态Fig. 15 Fracture and failure patterns of the simulated sample

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