张 飞，张东宾，王宇伟

(中国矿业大学(北京)力学与建筑工程学院，北京 100083)

0 引 言

在地矿开采过程中，往往对岩体施加爆炸、冲击等强动荷载作用使其破碎，而岩体中富集的结构面影响其动态力学行为，给工程活动的顺利进行带来了诸多困扰。岩体中的结构面形式多样，如孔洞、节理、裂隙等，其几何特征对岩体的力学特性影响也各有不同。众多学者对含结构面岩体的动态力学特性开展了广泛的研究工作。

谭卓英等[1]对含圆形孔洞结构面岩体的冲击动态响应进行了实验模拟研究，结果表明孔洞结构对岩体冲击破坏形式、峰值应力水平等有一定影响。朱哲明等[2]对含缺陷岩体的爆炸动态响应进行了研究，发现不同缺陷种类的岩体动态强度和破坏范围有所差别。杨仁树等[3]通过实验分析了含贯通裂纹的类岩石介质动态破裂行为和预应力场中含层理类岩石介质爆破破裂的动态响应。王蒙等[4]基于SCSCC试件对岩石复合型裂纹的动态扩展特性进行了实验研究。YAVUZ等[5]对含多裂纹无限大平板进行分析，得出裂纹间距和长度对应力强度因子的影响规律，并给出应力强度因子计算式。GREGOIRE等[6]在进行冲击试验时发现动态裂纹扩展时有止裂又继续扩展的现象。ZHU等[7]对巷道围岩的动态响应问题进行了数值模拟和实验研究，并分析了应力波作用下的动态响应及加载方式对动态响应的影响。AVACHAT等[8]对岩体的动态断裂行为引入高速摄影法进行了影像学分析。FUNATSU等[9]在分析岩体断裂破坏时，考虑了压力和温度等因素对砂岩动态破坏的影响。魏晨慧等[10]和谢冰等[11]对节理角度和间距变化时的岩体爆炸动态响应进行了数值模拟研究。以上诸多研究从多个角度对岩体的动态力学行为进行了研究，分析了动态应力强度因子、裂纹扩展行为等动力学特性。本文基于ABAQUS数值分析平台，建立模型，分析岩体结构面几何特征变化时的冲击动态响应，以期找到结构面几何特征对岩体冲击动态响应的影响。

1 应力波衰减规律

应力波在岩体等非连续性介质中传播时，往往受地质条件、结构面几何特征等因素影响，其反射透射及衰减状况尤为复杂，众多学者提出了一些理论模型，如接触界面模型、完好黏结界面模型、位移不连续体模型等。基于这些理论模型，应力波经过非连续界面时的传播规律了较为完善的理论支撑，如李业学等[12]推导了应力波穿越节理时反射透射系数FR1、FT1的解析解，其受节理等效刚度、波阻抗、频率等因素影响，见式(1)和式(2)。

(1)

(2)

应力波在岩体中传播时，受到阻滞作用发生衰减，携带能量耗散；岩体某一质点位移表示为式(3)。

u(x,t)=Aei(ωt-Ux)

(3)

其中：

(4)

(5)

(6)

式中：A为质点振幅；λ为波长；G为剪切模量；η为黏滞系数；ρ为密度；ω为应力波频率。

联立式(4)～式(6)可得式(7)。

U=β-iα

(7)

其中：

α=

β=

将式(7)代入式(3)可得式(8)。

u(x,t)=Ae-αxei(ωt-βx)

(8)

由此可知，应力波在岩体中传播时以式(8)的函数关系衰减，并且与岩体的黏滞系数、剪切模量、密度、应力波频率等因素有关。

2 模型建立

岩体本构模型选择德鲁克-普拉格准则(以下简称“D-P准则”)作为失效准则。D-P准则同时量化了体积应力、剪应力和中间主应力对岩石弹塑性强度的影响，能理想地反映岩石的力学行为。

为分析岩体结构面几何特征不同时冲击荷载作用下的动态响应，建立A组、B组、C组三组模型。平面应力状态如图1所示，外观尺寸150 mm×50 mm，上部中点受冲击荷载，速度30 m/s，下部及左右部设置为约束边界，x、y、z方向位移为0；网格划分时，单元形状为四边形，采用进阶算法(advancing front)生成均匀网格，单元属性为四结点双线性平面应力四边形单元；为便于分析，取分析点及分析路径，以A组模型为例，模型A1结构面左侧端部分析点记为A1L，右侧端部分析点记为A1R，下部中间分析点记为A1D，左右两侧以端部为起始点竖直向下方向上分别取分析路径记为path1和path2，其他模型同理。

图1 模型示意图Fig.1 Diagram of model

选用ABAQUS/Explicit，其本质上是基于动力学方程求解，对于爆炸、冲击的数值模拟具有很好的仿真效果。

3 结果分析

3.1 结构面厚度变化时的动态响应

A组模型的结构面厚度值分别为0.1 mm、0.2 mm、0.3 mm、0.4 mm，截取模型中间部分如图2所示。

图2 A组模型示意图Fig.2 Diagram of group A model

四个分析点A1D、A2D、A3D、A4D的Mises应力值随时间变化如图3所示。以冲击块体接触模型为零时刻，在70 μs左右时，Mises应力值出现峰值，分别为12.52 MPa、9.89 MPa、9.80 MPa、9.34 MPa。 模型A1结构面厚度为0.1 mm，分析点A1D的Mises应力峰值高于同组其他3个模型，但其他3个模型此时峰值差别不大，这可能是结构存在某一特定厚度使应力波没有穿透结构面所导致。随着时间的延后，分析点A1D的Mises应力值震荡下降，在200 μs之前，其值基本上大于其他3个模型分析点，200 μs左右时，其他3分析点应力峰值高于分析点A1D，这是因为应力波发生绕射，导致结构面下部分析点处的应力值又有所上升，200 μs之后，Mises应力值震荡下降。

图3 下部分析点应力值Fig.3 Stress of down monitoring point

以模型A1为例，在分析路径path1上，其Mises应力值随距离衰减曲线如图4所示。拟合函数为f(x)=ax-b，三个时刻的拟合函数关系式分别为：f1(x)=28.13x-0.43，f2(x)=23.72x-0.43，f3(x)=18.06x-0.40，符合应力波衰减规律。分析路径path2情况与此类似。

图4 应力值衰减曲线Fig.4 Curve fitting of stress attenuation

模型A1应力云图变化如图5所示。60 μs时结构面两侧以端部为中心，并主要沿下侧、外侧逐渐扩散，可以观察到结构面下部区域云图颜色偏冷色，这说明应力波绕射、透射到此区域较少。120 μs左右，应力云图持续保持以结构面两端为中心，主要向下侧、外侧扩展，向结构面下部中间集中程度并不大，这正是所取分析点Mises应力值震荡变化的直观现象，其根本原因是应力波的持续传播，之所以出现不只一次的峰值，是因为结构面尖端集聚能量释放的结果，180 μs左右也基本呈现此现象，云图颜色逐渐转冷色调。 实际岩体承受强动荷载时，弱面处出现应力集中，集聚较大应变能，并常从此起裂破坏。

图5 应力云图变化图Fig.5 Changing diagram of stress cloud

3.2 结构面长度变化时的动态响应

B组四个模型结构面长度分别为10 mm、20 mm、30 mm、40 mm，截取模型中间部分如图6所示。

图6 B组模型示意图Fig.6 Diagram of group B model

四个分析点B1D、B2D、B3D、B4D的Mises应力值随时间变化见图7。分析点B1D和分析点B2D变化规律相似，先震荡上升，200 μs左右达到峰值，分别为18.02 MPa和16.15 MPa，后者较小，随后震荡下降，这是因为随着结构面长度的增加，应力波较难以绕射到结构面下部区域。同样，分析点B3D的峰值和整体应力值水平更低。但随着结构面长度的进一步增加，情况有所不同，模型B4的结构面长度为40 mm，分析点B4D在100 μs左右时达到峰值19.23 MPa，随后有很短暂的下降后又快速上升，在160 μs左右时再次达到峰值18.86 MPa，随后震荡下降。

图7 下部分析点应力值Fig.7 Stress of down monitoring point

以模型B4为例，在分析路径path1上，其Mises应力值随距离衰减曲线如图8所示。拟合函数为f(x)=ax-b，三个时刻的拟合函数关系式分别为：f1(x)=50.44x-0.59、f2(x)=42.47x-0.60、f3(x)=46.88x-0.60，符合应力波衰减规律。分析路径path2情况与此类似。

图8 应力值衰减曲线Fig.8 Curve fitting of stress attenuation

Mises应力云图变化如图9所示，90 μs左右，结构面下部开始出现“水滴”状颜色转暖色区域，随后快速消失，在150 μs左右时再次出现此区域；同时可以观察到结构面侧面和下侧面接触，这说明结构面的应变积累，结构面的竖向位移值超过了结构面厚度，相当于节理的等效刚度降低，透射系数增加，在模型建立阶段结构面的上下侧表面相互作用设置为硬接触，上表面将应力传递到下表面，应力波发生了透射现象。

图9 应力云图变化图Fig.9 Changing diagram of stress cloud

3.3 结构面曲率变化时的动态响应

C组模型结构面变化情况如图10所示。为直观表示，图中数值为结构面的曲率半径大小；由曲率半径r和曲率k之间的关系式k=1/r可知两者呈反比，即曲率越来越小，依次为500、333、250、200。

图10 C组模型示意图Fig.10 Diagram of group C model

取结构面左侧端部四个分析点C1L、C2L、C3L、C4L，其Mises应力值随时间变化情况如图11所示。当结构面曲率越来越小时，Mises应力值在60 μs之前的增加速度越来越大，结构面曲率减小，结构面趋于平直，应力波传播到端部的距离缩短，因此端部的应力值会在较短时间内开始快速增加。分析点C1L在120 μs左右时第一次达到较大值为35.78 MPa，随后震荡下降；分析点C2L在60 μs左右时首次达到峰值42.57 MPa，随后震荡下降，但随着结构面曲率的进一步减小；分析点C3L在50 μs左右时便首次出现峰值，大小为42.61 MPa；分析点C4L的峰值出现时间更早，基本上在50 μs左右，这和结构面几何特征是平直时的端部分析点出Mises应力值峰值出现时间基本一致，这说明随着曲率减小，结构面趋于平直，其端部Mises应力值在应力波传播初期有着相似的变化情况。分析点C3L在首次峰值出现后有短暂时间的缓慢下降，趋势较为平直，这可能是时间间隔不够短引起的误差，其下降趋势应该是震荡下降，分析点C4L出现不只一次峰值后下降。

图11 左侧分析点应力值Fig.11 Stress of left monitoring point

应力波在传播过程中携带的能量主要转化为应变能储存在介质内部，取分析点处的应变能随时间变化情况如图12所示。结构面左侧端部的应变能随时间变化趋势基本和Mises应力值随时间变化趋势保持协调。分析C1L在100 μs左右时应变能积累到最大值，随后趋于下降；分析点C2L的应变能在60 μs左右时达到峰值，随后震荡下降；分析点C3L和C4L的应变能在50 μs左右快速上升到峰值，随后趋于较平缓下降，在200 μs左右时快速降低；可见随着结构面曲率的减小，结构面端部积累应变能的情况有所不同，曲率大时，也即结构面越来越“圆”，此时应力波从结构面弧所对应的弦的中垂线方向入射，结构面反射较多的应力波，应力波携带的能量也就不容易在此聚集，随着结构面趋于平直，应力波反射情况趋同于直线型结构面。

图12 左侧分析点应变能Fig.12 Strain energy of left monitoring point

4 结 论

1) 在以结构面端部为起始点竖直向下的分析路径上，Mises应力值以函数关系式f(x)=ax-b衰减。

2) 结构面厚度增大，应力波难以发生透射，下部分析点应力值水平降低；结构面长度增加，上下表面容易发生闭合现象，应力波发生透射，下部分析点应力值会增大。

3) 结构面曲率减小，即结构面趋于“平直”时，结构面端部应力峰值出现时间早，说明更容易积累应变能；应变能与应力值变化情况保持协调。