基于SHPB的复合岩样动态力学特性数值模拟

常玉林， 温 森，张建伟

(1.河南大学 岩土与轨道交通工程研究所，河南 开封 475004；2. 河南大学 土木建筑学院，河南 开封 475004)

1 研究背景

近年来， 随着经济建设的不断发展， 有关动态荷载下岩体的破坏问题越来越突出， 如工程爆破、 地震等动荷载作用下的地下工程问题也受到越来越多学者的关注。 因此， 研究动荷载下岩石动力学性质是非常有必要的[1-3]。 在岩石动破坏问题的研究中， 数值模拟方法受到许多学者的青睐。 与其他方法相比， 数值模拟方法不但能克服许多现有试验方法的缺点， 还可节省时间及试验经费等； 此外， 它还可充分考虑岩体的非均匀性及其赋存的复杂应力环境， 因此目前不少研究人员采用该方法对岩石的动力特性进行了研究。 杨岳峰[4]利用RFPA2D软件， 建立了岩石在动力作用下的细观损伤模型， 并系统阐述了该方法的基本思想、 理论基础和分析计算流程。 李超等[5]采用RFPA2D软件， 模拟了节理岩体的动态破坏过程， 研究了节理倾角、 节理条数、 节理贯通度及应力波峰值对岩体动态破坏过程的影响规律。 左宇军[6]采用RFPA2D对试样在不同应力波延续时间、 围压及应力峰值作用下的破坏进行了分析。 宫凤强[7]利用LS-DYNA软件， 分别从轴向、 径向传播2个方面分析了不同应力波在5种不同直径杆中的传播规律。 Bertholf等[8]和Zhang等[9]利用RFPA2D对SHPB试验系统进行了模拟， 对比分析了惯性及端面摩擦效应对试验结果的影响程度。 Liao等[10]采用ANSYS软件通过对不同冲击速度、 冲头形状与长度及不同试样尺寸计算结果的对比， 揭示了在SHPB试验中应力波形的影响因素。 Yan等[11]采用二维离散元软件UDEC分别分析了凝灰岩与玄武岩在巴西劈裂试验中的动态响应， 得出了岩石动态抗拉强度与应变率有着典型的相关性的结论， 同时还研究了不同性质的岩层之间结构面对其动态抗拉强度的影响。 Li等[12]和Zhu等[13]采用PFC2D模拟了异形冲头加载的SHPB试验，研究了材料动态屈服强度放大系数的影响因素。 Xu等[14]利用颗粒流程序PFC2D建立了一种半圆形弯曲试样的SHPB数值试验程序， 从能量变化的角度分析了动态拉伸破坏下岩石的能量耗散。 李晓峰等[15]在试验的基础上， 利用PFC2D软件模拟了应变率>103/s-1时岩石材料的动态力学响应及破坏特征。以上众多学者主要从不同角度、 考虑不同因素， 针对不同单一岩石的动态拉、 压力学及破坏特征进行研究， 而对工程中经常遇到的复合岩层的动态力学特性研究却鲜有报道。 基于此， 本文在前人研究的基础上， 采用PFC2D建立复合岩样SHPB数值试验模型， 研究不同冲击速度及岩层倾角对复合岩样动态力学特性的影响。

2 数值模型的建立

二维颗粒流软件PFC2D模型是由一系列离散的圆形颗粒聚集而成，颗粒运动遵循牛顿第二定律，试样的宏观力学特性通过颗粒间的接触变形体现。本文模拟采用能较好地反映岩石材料力学特性的平行粘结模型[16]。

本文在Tien等[17]的复合岩层试样单轴抗压试验的基础上，根据其试验结果进行材料细观参数的标定，研究复合岩层试样的动态力学特性。

2.1 模型的建立

根据国际岩石力学学会(ISRM)建议的试验方法，数值试验中岩样大小选取为50 mm×50 mm，复合岩样每层的厚度与Tien等[17]的静态试验保持一致，设定为5 mm。为提高计算效率，SHPB的杆件模型应在满足应力波测量的基础上尽可能小，故本模拟中入射杆、透射杆长度分别取为1.50 m和0.75 m(如图1(a)所示)。试样与杆件的接触面处，颗粒均匀排列(如图1(b))，以降低杆件的横向弥散效应[11]。采用异形冲头输入应力波，试验时在冲头右侧施加初速度，并通过改变不同的冲击速度来实现不同的应变率加载。

图1 数值模型及其局部细节Fig.1 Numerical model of SHPB and local detailsof the model

2.2 细观参数标定

杆件的颗粒半径范围选择为0.9～3.0 mm，根据质量守恒定律推算颗粒密度为7 894.7 kg/m3。颗粒的接触弹性模量和刚度比分别为343 GPa和2.8。冲击试验中杆件始终处于弹性状态，其粘结强度设定为10100MPa。

合理的细观参数是模拟岩石材料力学特性的关键。复合岩层试样细观参数主要由3部分组成：材料A、材料B及层理面的细观参数。本次模拟细观参数标定步骤为：①室内单轴压缩试验标定材料A、材料B的细观参数；②室内巴西圆盘试验对细观参数进行校验与修正；③将A,B材料复合，通过复合材料的室内单轴压缩试验标定层理面的细观参数；④改变层理角度，与室内试验结果对比，修正层理细观参数；⑤动态试验校验并修正。本次数值试验所用材料A,B及对应的复合岩样的力学参数取自Tien等[17]的研究成果。表1列出了标定得到的细观参数，表2、表3分别列出了均质试样与复合岩层试样的室内试验、数值试验结果的对比。对比结果说明，数值模拟与室内试验结果具有较好的一致性，标定的材料细观参数合理。

表1 岩石材料细观参数Table 1 Meso-mechanical properties of rock

表2 均质试样室内试验、模拟参数对比Table 2 Comparison of parameters of homogeneoussample between laboratory test and numerical simulation

表3 复合岩层室内试验、模拟单轴抗压强度对比Table 3 Comparison of strength of mixed rock stratabetween lab test and numerical simulation

图2 入射杆与透射杆中应力波形Fig.2 Stress waveforms ofincident bar andtransmission bar

2.3 数值模拟程序验证

为了确保模拟程序的合理性，在研究前先对本文的程序进行验证。图2为数值模拟过程中入射杆、透射杆中的波形图。该波形是冲击速度V=10 m/s时产生的，模型的应力波形与实际SHPB试验中的半正弦波加载方式一致。利用三波法[17]进行验证，图3给出三波验证图。

图3 三波法验证应力平衡Fig.3 Verification of stress balance by three-wavemethod

由图3可以看出，透射波与反射波和入射波的叠加波基本重合；随着试样加载，两者的差别在一个较小范围内波动，说明试样在加载过程中满足应力平衡条件，数值模拟程序合理。

3 层理面倾角对复合岩样动态响应的影响

定义复合岩样层理面倾角为α，图4为试样在数值模拟中的放置图。应力波传播方向与层理面的夹角同α为互余关系。模拟中α依次取为0°，15°，30°，45°，60°，75°，90°，每层厚度为5 mm，冲击速度设定为10 m/s。

图4 岩样放置示意图Fig.4 Schematic diagram of rock sample placement

3.1 层理面倾角对动态力学参数的影响

图5给出了7种不同倾角的复合岩样的动态应力-应变曲线。数值模拟中没有考虑材料的初始微缺陷，因此峰前曲线接近线性变化，无材料的初始压密过程。α=0°，15°，30°所对应的3条曲线接近重合，表明α较小时，层理面倾角对复合岩层动态力学特性的影响较小，而当α≥45°时，倾角对动态力学性能产生的影响较大。

图5 不同α复合岩层试样的动态应力-应变曲线Fig.5 Dynamic stress-strain curves of mixed rockstrata with different values of dip angle α

为更加清晰地分析α的影响，绘制了试样动态峰值应力随α的变化图，如图6。从图6中可以看出，α=60°时动态峰值应力最低，为61.66 MPa；α=90°时动态峰值应力最高，达到81.82 MPa，两者相差32.69%。试样的峰值应力随着α的增大呈现一种先减小后增大的近“V”字型。

图6 动态峰值应力与α的关系Fig.6 Relationship between peak stress and α

定义峰值应力所对应的应变为εp，图7(a)为εp与α的关系曲线。复合岩样的εp在0.539%～0.615%之间，最大值与最小值相差12.3%。随着α的变化，εp呈现小幅波动变化。

复合岩层试样的动弹性模量与α的关系见图7(b)。图7(b)中规律与动态峰值应力随α变化规律基本一致，最大与最小动弹性模量之间相差27.42%。

图7 εp和动态弹性模量与α的关系Fig.7 Relationship of strain at peak stress and dynamicelastic modulus against α

3.2 层理面倾角对动态破坏模式的影响

为了分析α对动态裂纹的影响，数值试验中对裂纹数目与模式进行了统计。图8为根据统计结果绘制的柱状图。从图8可见，复合岩层的动态裂纹以拉伸裂纹为主，剪切裂纹为辅。随着岩层倾角的增大，裂纹的总数在2 466～2 670之间呈小幅度波动变化，表明在同样的冲击速度下，α对于复合岩层动态破坏的破碎程度的影响并不明显。α较小时(如α=0°，15°，30°)，试样剪切裂纹较少；随着α的增加，即从45°至75°时剪切裂纹数目逐渐增多；而在90°时，剪切裂纹数目又减少。

4 应变率对复合岩样动态响应的影响

以α=30°为例，分别对冲头施加V=7.5，10.0，12.5，15.0，17.5，20.0 m/s 6种不同的冲击速度，研究复合岩层试样在不同应变率下的动态力学特性。

4.1 应变率对动态力学特性的影响

数值试验中，6种不同的冲击速度对应的试样应变率分别为66.72，99.41，115.48，140.23，168.90，188.98 s-1。应变率随着冲击速度的增大而增大，其对应的应力-应变曲线如图9所示。由图9可知，在该应变率范围内，曲线峰前基本保持弹性变化，其峰值应力与εp均随着应变率的增大而增大。

图9 不同应变率下复合岩层的动态应力-应变曲线Fig.9 Dynamic stress-strain curves of mixed rockstrata at different strain rates

为了更清楚地分析应变率对动态峰值应力、εp和动弹性模量的影响，将图9中不同应变率对应的峰值应力和εp分别提取出，绘制图10。图10(a)显示α=30°时，应变率从66.72 s-1增大到188.98 s-1的过程中，动态峰值应力从66.15 MPa增加到98.41 MPa，增幅达到48.76%；可见试样的峰值应力表现为明显的率相关性，并且其变化趋势接近于线性增长。图10(b)显示εp随着应变率的增大而增加，其数值从0.505%增大至0.729%，增幅为44.35%。从图10(c)可知，试样的动态弹性模量在12.81～13.52 GPa内有小幅波动，最大动弹性模量与最小动弹性模量之间相差5.5%，并且没有一个明显的变化趋势，说明在该范围内，应变率对复合岩层材料的动态弹性模量影响相对较小。

图10 动态峰值应力、εp、动态弹性模量与应变率的关系Fig.10 Relationship of dynamic peak stress,peak strain,and dynamic elastic modulus against strain rate

4.2 应变率对破坏模式的影响

为了分析应变率对动态裂纹的影响，根据统计结果绘制了图11。由图11可知，α=30°时，无论应变率如何变化，拉伸裂纹都占了绝对优势，表明试样的动态裂纹以拉伸裂纹为主；并且随着应变率的增大，裂纹总数也逐渐增多，表明随着应变率的增加，复合岩层试样动态破坏的破碎程度越来越严重。

图11 动态破坏裂纹数目与应变率的关系 Fig.11 Relationship between number of dynamiccrack and strain rate

5 关于破坏模式的讨论

前面讨论了岩层倾角、应变率对试样破坏模式的影响，分别统计了拉伸裂纹、剪切裂纹数目。然而实际工程中并不能完全以裂纹数目来判断宏观破坏模式。宏观破坏模式的判断除了看细观裂纹条数，还要观测裂纹的贯通情况及裂纹的位置。

PFC软件中的裂纹数目是指颗粒之间连接断开的数目，而颗粒间的连接数由颗粒的数量决定。以α=30°和α=60°、冲击速度为7.5 m/s为例，将其细观破坏模式与材料A均质试样进行对比，三者的细观裂纹扩展情况如图12所示。图中黑色为拉伸裂纹，红色为剪切裂纹，灰色颗粒代表材料A，白色颗粒代表材料B。可见，均质试样的细观破坏以沿应力波传播方向的拉伸裂纹为主；倾角为30°的试样以穿越层理面的拉伸裂纹为主，伴随少量层理面附近的剪切裂纹，虽然有少量裂纹沿层理面发展，但没有完全贯通，破坏模式与均质试样相似，说明此时层理面对试样破坏模式的影响有限。倾角为60°的试样既有穿越层理面的拉伸裂纹，又有沿层理面发展的剪切裂纹。剪切裂纹主要在层理面中出现，虽然其数目显著少于拉裂纹数，但却沿某些层理面实现了贯通。因此此时剪裂纹在破坏中起到了控制作用，而拉裂纹加剧了试样的破碎程度。上述分析说明试样动态破坏模式会随α的变化而有所不同，并且最终的破坏模式受细观裂纹种类、数量、位置以及贯通度的综合影响。

图12 不同试样的细观裂缝扩展Fig.12 Propagation of meso-cracks of different samples

6 结 论

(1)在动载荷作用下，岩层倾角α对复合岩层试样的力学特性有较大影响。随着α角度的增大，动态抗压强度、动态弹性模量整体先减小后增大、呈近“V”字型变化，在α=60°附近时达到最小值。

(2)复合岩层的动态力学特性有着明显的率相关性，其动态强度、动态峰值对应的应变均随着应变率的增加呈现近似线性的增长，而动态弹性模量受应变率变化的影响较小。

(3)复合岩层试样破坏裂纹以拉伸裂纹为主，α较大时剪切裂纹数目会随着α的增大逐渐增多；随着应变率的增大，试样的动态破坏破碎程度越来越严重。

(4)试样的细观破坏模式受岩层倾角的影响，并且最终宏观破坏模式受细观裂纹数目、位置及贯通情况的综合影响。