冲击作用下煤岩动态破坏机理的FDEM模拟研究

程树范， 高 睿， 曾亚武， 张嘉凡， 陈世官

(1. 武汉大学 土木建筑工程学院， 武汉 430072；

2. 西安科技大学理学院， 西安 710054；

3. 西安科技大学 建筑与土木工程学院， 西安 710054)

煤炭作为重要能源矿产对我国经济发展起到了举足轻重的作用[1]。在矿井建设和煤矿开采过程中不可避免地存在大型机械掘进、常规爆破等冲击扰动，随着开采深度的增加，相应的开采强度随之增加，由冲击作用诱发的矿井灾害愈发频繁。与准静态作用下的渐进破坏不同，冲击作用下煤岩破坏过程具有显著的突发性，难以有效的预防，因此有必要对煤岩冲击破坏规律进行系统研究。

分离式霍普金森压杆(Split Hopkinson pressure for, SHPB)装置，可有效地避免冲击载荷作用下对待测试件进行直接测量，是材料高应变率试验的重要手段[2-3]。近年来，利用大直径SHPB装置开展的各类岩石动态破坏规律研究已取得众多阶段性成果，如纪杰杰等[4]开展了冲击荷载作用下花岗岩和砂岩碎块的分形特征研究，认为分形维数可实现对岩石在冲击荷载作用下的破碎特性、力学特性和破碎耗能特性的定量研究；
苗磊刚等[5]利用SHPB装置研究不同应变率下岩-煤-岩组合体破碎特征和碎块分形变化规律；
李成杰等[6]则进行了煤岩组合体的冲击压缩试验，分析组合体的破碎特征和能量耗散规律，发现煤岩复合体发生动力灾害的能量阈值更低，实际工程中的危险性更大。刘少虹等[7]针对冲击地压防治，开展了一维的动静组合SHPB试验，研究表明煤岩动态强度和碎片分形维数随加载应力波能量的增大而增大，随静载的增大则呈现先增大后减小的趋势，说明煤层卸压有利于提高其对高应变率荷载的承载能力；
王春等[8]则采用该装置开展了一维静载与频繁扰动共同作用下的动力学试验。平琦等[9]进行了被动围压条件下的煤矿砂岩SHPB试验，试验结果表明被动围压可以大大改善砂岩的动态力学性能，降低岩爆风险；
解北京等[10]利用SHPB试验比较了煤岩与砂岩的动力学特性，认为层理煤岩的弹性模量小，冲击荷载作用下的塑性变形能力强。受限于测试技术，对岩石试件冲击破坏过程的实时观察还无法实现，数值模拟技术在一定程度上解决了这一困难。研究表明，根据试验数据进行标定的数值模型一般具有较高可靠性，如凌天龙等[11]通过室内试验对砂岩的HJC(Holmquist-Jchnson-Cook)模型主要进行标定，并进行了SHPB冲击压缩的数值模拟，模拟得到的岩石动态应力-应变曲线与试验高度一致；
针对煤岩，李成武等[12-13]则基于HJC模型和单元失效法，模拟了煤岩的动态破坏，研究了煤岩冲击破坏规律和被动围压的强化效应。就数值模拟方法而言，凌天长等将岩石动态破坏规律描述为材料的非线性问题，而忽略裂纹扩展引起的几何非线性，虽然宏观上与试验结果一致，但就岩石材料的脆性特征，其机制值得商榷；
李成武等采用的侵蚀算法，即单元失效模型，虽然在一定程度上可以反映煤岩的破碎过程，但在较高应变率加载条件下，煤岩试件可能会出现过大的质量损失，进而导致计算失真。耦合有限-离散元方法(finite discrete element method, FDEM)，通过零厚度节理单元的失效来模拟岩石裂纹扩展过程，有效解决了质量不守恒的问题。Gui等[14-15]在离散元的基础上，通过Cohesive单元的嵌入，在UDEC(universal distinct element code)中模拟了岩石的多种动态破坏过程，认为该模型可以模拟岩石材料的应变率效应；
Wu等[16]采用类似的方法模拟了滚刀破岩，分析了围压对刀具破岩效果的影响；
Wang等[17]和Wang等[18]则采用FDEM方法分别模拟了锤头破岩和岩石切削，取得了较好效果。吴志军等[19]比较了FDEM和侵蚀算法的模拟效果，认为FDEM避免了单元删除带来的能量异常，模拟效果更真实。

本文通过SHPB冲击压缩试验，研究了煤岩强度和弹性模量的应变率相关性；
并采用零厚度节理单元模拟煤岩内节理结构，建立了煤岩FDEM模型，通过室内试验数据对内聚力单元参数进行了标定；
最后讨论了FDEM方法在模拟煤岩动态破坏时的可靠性，采用本文标定的参数，模拟了煤岩的冲击破坏过程，并分析煤岩破坏规律。

1.1 试验方案

试验所用的煤岩样品采自陕北某矿区，经取芯、切割和打磨等工序，将新鲜的大块样品分别加工成Φ50 mm×100 mm、Φ50 mm×50 mm和Φ48 mm×25 mm三种尺寸试件。前两规格的试件主要用于准静态的单轴压缩和巴西劈裂试验，第三种规格的则进行动态的SHPB试验，得到的煤岩基础力学参数如表1所示。

表1 煤岩基础力学参数Tab.1 basic mechanic parameters of coal

SHPB装置，是进行材料中高应变率试验的推荐设备，包括压杆试验系统、能量回收系统、信号采集系统和数据处理系统四个部分，如图1所示。

图1 SHPB试验装置原理图Fig.1 Schematic illustration of SHPB experimental device

试验中通过压杆的应力波信号对试件动态应力-应变进行间接测量，对于脆性材料需采用“三波法”进行数据处理，其表达式为[20]

(1)

(2)

(3)

本次试验使用的SHPB试验装置直径为50 mm，配套的超动态应变仪采集频率为105Hz。试件安装前对试件的端面平行度和垂直度进行了复查，并涂抹凡士林，减小压杆和试件间的摩擦作用。为减少弥散效应所产生的P-C震荡的影响，在入射杆前端粘贴了紫铜滤波片。子弹的冲击冲击速度由气炮压力控制，经测试气炮压力大于0.25 MPa后(子弹速度大于2 m/s)，试件将出现不同程度的破坏。据此，本次SHPB试验共分五组，每组三个试件，各组的预设气压值分别为0.25 MPa、0.30 MPa、0.35 MPa、0.40 MPa和0.45 MPa。剔除无效数据后各分组的子弹速度平均值分别为1.94 m/s、2.42 m/s、3.08 m/s、3.50 m/s和4.15 m/s。

1.2 试验结果分析

根据采集的应变信号，利用式(1)～式(3)的“三波法”对各试验分组的应变信号进行处理，消去时间变量后得到不同应变率条件下的煤岩动态应力-应变曲线如图2所示。

图2 煤岩动态应力-应变曲线Fig.2 Dynamic stress-strain curves of coal

SHPB试验通过子弹和入射杆的弹性撞击加载，其加载强度(入射波峰值)可通过撞击速度控制，而加载时间(入射波长度)则主要由子弹长度决定[21]，本次试验使用的子弹长度为20 cm，入射波的加载时长约为220 μs，由于加载时长有限，应变率较低时(79.78 s-1和84.36 s-1)的加载总应变不足以使试件完全破碎，因此动态应力-应变曲线包含短暂的回弹段，相应的峰值应力和变形模量并不是该应变率条件下的动态强度和弹性模量。根据试验结果，应变率98.05 s-1、119.22 s-1和135.85 s-1的煤岩动态强度因子DIF(dynamic strength factor)(动态抗压强度与准静态抗压强度之比)分别为1.92、2.08和2.23，煤岩抗压强度的应变率相关性比较明显。根据已有的经验模型[22]，在较高应变率条件下，DIF与加载应变率间存在如下函数关系

(4)

根据本文的试验结果，模型参数β的拟合值为0.13 2，对应的确定系数R2=0.953，可见式(4)的经验模型也适用于本文煤岩。

煤岩动态弹性模量大于准静态弹性模型，而动态弹性模量(切线模量)受应变率影响不明显，表现为较弱的正相关性，随着应变率由98.05 s-1增加到135.85 s-1煤岩动态模量仅由6.30 GPa增长为6.75 GPa，增长幅度为7.1%。

图3 煤岩破坏形态Fig.3 Failure modes of coal samples

由图3可知，随着冲击速度的增加，煤岩的破碎程度和破坏形态有明显的区别。当冲击速度在2 m/s时，加载完成后的煤岩试件基本完整，仅在入射面出现了少量非贯穿裂纹和角部脱落，因此具备较强的卸载回弹能力；
当冲击速度达到2.5 m/s时，试件发生断裂，并破碎为若干碎块，其中较大的碎块约占试件总体积的1/2，可以完成部分卸载，因此应力-应变曲线包含有较短的回弹过程。随着冲击速度进一步增加，煤岩的破碎更加明显，试件在加载阶段已破碎为数块，不再具备卸载回弹能力，煤岩试件在入射面一侧破碎较透射面严重，沿加载方向的裂纹较多，并伴随有片状脱落。

虽然通过破坏形态和应力-应力曲线，可以简单地分析煤岩的破坏规律，但由于冲击加载过程十分短暂，仅为0.1 ms数量级，且煤岩内部无法直接观测，关于煤岩破坏过程的研究仍需要借助于数值模拟来实现。

2.1 内聚力模型理论

在宏观尺度上，岩石材料具有连续介质的特点，而在细观尺度上，又具有离散的颗粒聚合结构，采用何种理论来分析岩石的动态破坏规律，便成为建立岩石数值模型的关键所在。内聚力模型属于耦合有限-离散元方法(CFDEM，combined FDEM)范畴，通过内聚力单元(cohesive element)与弹性实体单元的共节点连接将离散颗粒连接成连续整体，兼顾了有限元和离散元方法的优势，是模拟岩石材料破坏的理想方法。内聚力模型的生成方法，如图4(a)所示。其中插入的内聚力单元没有厚度，但在几何上是二维的，其实质是横观各向同性的实体单元，具有两个位移分量和两个应力分量，即法向的位移δn和应力tn；
切向位移δs和对应方向上的应力ts。在法向和切向受力均可以产生如图4(b)所示的相对位移。

(a) 单元的生成

(b) 单元破坏模式图4 内聚力模型Fig.4 Cohesive model

相较于传统有限元法中借助实体单元失效来模拟岩石的破碎，内聚力模型可以通过零厚度单元的失效，更加真实地模拟了这一过程，避免了质量不守恒和能量异常损失。

在切向，内聚力单元的应力与位移满足式(5)的本构方程。

(5)

式中：ts为切向应力，δs为切向位移，ks为切向刚度，Ds为切向的损伤变量；
δs0和δsf分别为切向的损伤起始位移和失效位移。

在法向上，内聚力单元具有与切向相似的应力-位移关系，但只考虑拉应力产生的损伤和失效，本构方程如式(6)所示。

(6)

式中：kn为法向刚度，δn为法向相对位移，Dn为法向的损伤变量；
δn0和δnf分别为法向的的损伤起始位移和失效位移。

由式(5)和式(6)可知，当δn>δn0或δs>δs0时内聚力模型将产生损伤，对于损伤变量Dn和Ds，本文均采用线性损伤演化方程进行描述，即：

(7)

(8)

式中：δn max和δs max分别为内聚力单元在加载过程中法向和切向的最大位移值。

根据式(5)～式(8)，内聚力单元在法向和切向的应力-位移关系的双线性模型描述，如图5所示。

图5 内聚力单元本构模型Fig.5 Constitutive model of cohesive element

由图5可知，失效位移δs0和δsf可以通过断裂能来表示

δnf=2GIC/knδn0

(9)

δsf=2GIIC/ksδs0

(10)

式中，GIC和GIIC分别为法向和切向断裂能。

受荷过程中，内聚力单元除张拉破坏模式和剪切破坏模式外，还包括张拉-剪切复合模式，因此本文基于最大应力条件来考虑内聚力单元的复合损伤和失效。损伤起始条件为

(11)

失效条件可以表示为

(12)

式中，D为复合损伤变量。当D=Dn=1时，内聚力单元失效产生I型裂纹(张拉破坏)；
当D=Ds=1时，产生II型裂纹(剪切破坏)；
其他条件下的裂纹则属于I-II型复合裂纹。

2.2 内聚力参数的标定

内聚力模型中实体单元的力学参数可以直接通过试验获取，而内聚力单元的材料参数则需要借助单轴压缩和巴西劈裂试验进行标定。参数标定过程涉及数值试验和室内试验的对比，对网格质量有较高要求[23]，若大量的三角形边界在同一直线上，则对裂纹有很强的诱导作用，无法体现岩石材料的随机性。为真实地反映裂纹的扩展机制，本文采用非结构Delaunay算法生成网格，在保证网格质量的基础上避免了对裂纹扩展的过度人为干涉，如图6所示。

图6 非结构化网格划分Fig.6 Unstructured grid of coal sample

数值计算精度会随着网格尺寸l的减小而增加，经过试算，当网格尺寸l达到mm数量级后，计算结果基本稳定，随着网格尺寸的进一步减小，显式分析的计算精度虽然还会有小幅增加，但相应的稳定求解的时间步长也会减小，求解计算量显著增加，综合考虑模型精度和计算量，本文网格尺寸选取为1 mm。

内聚力单元的6个独立参数(法向、切向的刚度参数kn、ks，强度参数tn max、ts max和能量参数GIC、GIIC)，主要采用试错法进行标定，即当某组参数模拟得到的单轴压缩应力-应变曲线和抗拉强度均与试验结果一致，且破坏模式相符时，则认为该组参数是合理的，本文标定的模型参数如表2所示。

表2 内聚力单元参数Tab.2 Parameters of cohesive elements

模拟得到的单轴压缩应力-应变曲线如图7所示，设置的端面摩擦因数为0.1，试件为剪切破坏(双侧斜裂纹破坏)与室内试验基本一致。

图7 单轴压缩应力-应变曲线Fig.7 Stress-strain curves of uniaxial compressive test

本文模拟的巴西劈裂试验得到的煤岩抗拉强度为0.81 MPa，与室内试验相差小于5%，准确度较高，而破坏形态如图8所示。也与试验基本一致。

图8 巴西劈裂破坏模式Fig.8 Failure pattern of Brazilian disc test

可见表2中的参数在模拟煤岩准静态破坏时具有较高的可行性。

3.1 模型的建立

在LS-DYNA中建立的SHPB模型如图9所示。杆件系统的直径为50 mm，其中子弹长度为20 cm，入射杆长度为2.0 m，透射杆长度为1.2 m，均与试验一致，杆件间及杆件与试件的接触采用基于罚函数的面-面自动接触。杆件网格控制尺寸为1 cm，共有单元1 700个；
试件直径为48 mm，厚度为25 mm，网格控制尺寸为1 mm，共有实体单元1 694个，内聚力单元2 487个。

图9 SHPB试验的数值模型(LS-DYNA)Fig.9 Numerical model of SHPB test in LS-DYNA software

SHPB试验过程中，压杆系统始终处于弹性状态，采用线弹性材料进行模拟，弹性模量E0=210 GPa，泊松比υ0=0.3；
紫铜滤波片则采用附加Gruneisen状态方程的J-C(Johnson-Cook)模型进行模拟[24]，煤岩试件采用内聚力模型来模拟，参数在前文中已标定，按表1、表2选取。采用与试验相同的冲击速度，共进行了五组数值模拟。

3.2 破碎形态及应力-应变关系

数值模拟得到的试件最终破碎形态和应力应变曲线，如图10所示。

(a) v=1.98 m/s

(b) v=2.45 m/s

(c) v=3.02 m/s

(d) v=3.56 m/s

(e) v=4.16 m/s图10 数值模拟得到的应力-应变曲线Fig.10 Numerical model of SHPB test in LS-DYNA software

由图10可知，数值模拟得到的平均应变率和应力-应变曲线均与试验值十分接近。虽然内聚力模型本身没有与应变率相关的参数，但仍能通过惯性效应，即裂纹数量的增加来提高试件的宏观强度，说明本文数值模型在模拟煤岩动态破坏时也是可行的。

就最终的破坏形态而言，当冲击速度较低(v=2 m/s)时，试件在加载阶段于入射端角部和透射段中部出现少量剪切裂纹，加载过程结束后，试件仍具有较强的完整性，可以有效卸载，卸载过程无新裂纹生成，最终的破坏形态与室内SHPB试验结果一致。随着冲击速度增加，试件除角部的剪切裂纹数量有所增加，其破坏形态也发生了改变，当冲击速度达到2.5 m/s时，由于透射面附近反射的压力波与加载压力波叠加，产生张拉效应，试件于透射端出现了一条垂直于加载方向的张拉裂纹，进而导致试件出现层裂破坏，在煤岩的SHPB试验中也较常见。随着冲击速度的进一步增加，试件在入射端的剪切裂纹数量继续增加，并于试件中部产生沿加载方向的张拉裂纹。分析认为，轴向压力产生的泊松效应，导致试件在环向受拉，进而诱发张拉破坏。当冲击速度达到3.0 m/s后，试件的破坏形态基本稳定，破碎程度与冲击速度正相关。

3.3 破坏过程分析

数值模拟的最大优势在于可以对瞬时的破坏过程进行实时观测，记入射波前沿到达试件为T=0时刻，则根据数值模拟结果，冲击速度为3.0 m/s时，煤岩试件破坏过程，如图11所示。总历时1 200 μs。

图11 煤岩试件的破坏过程Fig.11 Failure process of coal specimen

相应的裂纹数-时间曲线如图12所示。

图12 裂纹数-时间曲线Fig.12 Crack number-time curve

在破坏类型上，煤岩的冲击破坏包括入射面附近的局部剪切破坏和整体张拉破坏两部分。在T=0时，入射波到达试件的入射端，T=300 ms时，入射波峰值进入试件，煤岩试件在入射端角部出现局部剪切破坏，岩石颗粒与母体脱离，向两侧飞出，其后随着应力波在试件内的继续传播，试件在入射杆侧的剪切裂纹数量继续增加，并出现少量张拉裂纹，裂纹处于稳定扩展阶段，此阶段内张拉裂纹尚未表现出扩展趋势；
在T=600 μs，透射波和反射波开始生成，试件出现一定程度的卸载，该阶段试件中部的多条张拉裂纹萌生，并迅速地向透射杆端延伸，裂纹进入快速扩展阶段，期间裂纹宽度也显著增加；
张拉裂纹在T=900 μs时刻贯穿整个试件，试件破坏为3个较大碎块及大量粉末，随后整体应力水平的降低，裂纹也进入缓慢扩展阶段，该阶段裂纹宽度继续增加，而裂纹数量基本稳定，仅在原有裂纹周边出现局部延伸；
最后在T=1 200 μs时裂纹不再增加，煤岩的冲击破坏过程结束。

煤岩的动态破坏过程与静载破坏过程有明显区别，在宏观上表现为强度的应变率相关性，和破坏形态的差异，在细观上则与煤岩结构相关。由图7和8可知，单轴压缩和巴西劈裂试验中煤岩的破坏主要表现为剪切裂纹和张拉裂纹的贯穿，产生的裂纹数量较少。而冲击破坏过程伴随着高密度的剪切破坏区和多条贯穿的张拉裂纹，单位体积内的裂纹总长度更多，因此煤岩冲击破坏时吸收的能量更多，于裂纹面释放的自由能更高，宏观上表现为应变率强化，同时试件的破碎也更严重。而弹性模量主要由基质材料决定，与裂纹面关系不大，因此煤岩弹性模量的应变率效应并不明显。

利用Φ50 mm SHPB装置进行了煤岩冲击压缩试验，研究了煤岩的动态破坏规律，并基于FDEM方法模拟了冲击破坏过程，主要得到以下四点结论：

(1) 煤岩的动态抗压强度与加载应变率间存在经验关系，应变率为98.05 s-1、119.22 s-1和135.85 s-1时其DIF分别为1.92、2.08和2.23；
煤岩动态弹性模量略大于静态弹性模量，但受应变率影响较小。

(2) FDEM可以模拟煤岩的裂纹扩展和破坏过程，通过煤岩的单轴压缩试验和巴西劈裂试验，可以对煤岩内聚力参数进行快速标定。由于惯性效应的存在，当网格尺寸合理时，FDEM模型也适用于煤岩动态破坏的模拟。

(3) 当冲击速度较低时，煤岩出现局部破坏后仍能保持一定的承载能力，卸载阶段有回弹现象；
冲击速度较高时，煤岩的冲击破坏包括加载段的局部剪切破坏和卸载阶段的整体张拉破坏，最终破碎为若干碎块和大量粉末。

(4) 煤岩的动态破坏过程与静载破坏过程有明显区别，在细观上，煤岩冲击破坏时裂纹密集，能量的吸收率和裂纹表面自由能释放率都较高，因此动态抗压强度在宏观上表现出较强的应变率相关性，同时动载作用下试件的破碎程度较准静态加载更严重。