基于平面黏结接触模型模拟平台巴西圆盘的动态响应和破坏机制

梁尔祝 耿 威 方定星 刘 鑫2

（1.鞍钢矿业爆破有限公司，辽宁 鞍山 114046；2.辽宁科技大学矿业工程学院，辽宁 鞍山 114051）

0 引言

采矿工程、隧道开挖、水利水电等工程涉及高应变率条件下岩石动态响应和破坏机制问题，了解岩石在冲击载荷作用下的裂纹萌生、演化、扩展和传播的细观特征，有利于掌握岩石动态破碎机理和控制岩石破碎效果。因此，深入探究高应变率条件下岩石的动态响应机制及破裂机理，对于岩石开挖工程的设计与施工具有一定研究意义。

霍普金森压杆试验技术已经是研究岩石在不同应变率条件下的力学性质和动态本构关系相对成熟的方法，其他相关的室内试验方面比如利用高速摄影记录岩石动态破坏全过程，借助数字图像相关方法分析岩石动态破坏的表面应变变化全过程[1-4]。在数值模拟方面，马江锋等[5]运用岩体裂纹扩展破坏二维分析程序DDARF模拟了在不同入射波加载情况下的大理石巴西圆盘试样破裂全过程，认为试样破裂是从加载端部起裂的；Mehrdad Imani等[6]利用二维颗粒流程序模拟了巴西圆盘试样在不同应变率条件下的破裂模式，阐述了巴西圆盘的4种破裂机制，提出剪切微破裂裂纹数量的增长可能是导致试件动态强度相应增加的原因。因此，岩石在动载荷作用下的破裂全过程和失效机理还需更详细的描述。

离散单元法是建立在微观力学原理上的、广义的颗粒流模型，被广泛应用于岩土材料、岩石材料、混凝土材料和沥青材料，颗粒间联结破坏可实现裂隙形成、贯通和发展，易于体现岩石破裂过程[7-8]，O'Sullivan[9]证实了离散单元法模拟动态岩石破坏方面的有效性。颗粒流模型是通过设置颗粒间接触的不同黏结模型实现相互黏结的，当作用力大于颗粒间的黏结力时，黏结断裂并以显示方式表示裂隙进而显示出来[8]。相比线性接触黏结模型和平行黏结接触模型的岩石的单轴压缩强度与拉伸强度的比值（即压拉强度比）范围，平面黏结接触模型的压拉强度比范围更广，更适合用于岩石材料[10]。

应变率对岩石破坏机理的影响是岩石领域最重要的方面之一，本研究利用平面黏结接触模型建立花岗岩平台巴西圆盘试样的数值模型，从细观角度深入探究高应变率载荷情况下花岗岩平台巴西圆盘试样的动态破裂模式和破坏机制。

1 SHPB试验简介

图1所示，SHPB试验系统主要由撞击杆、入射杆、透射杆、吸收杆、测速系统、动态应变仪、信号采集卡、计算机、氮气瓶等组成。本次SHPB试验采用50 mm直径的杆径，入射杆长度1 800 mm，透射杆长度1 800 mm，撞击杆长度300 mm，弹性模量210 GPa，密度7 800 kg/m3，入射杆端部贴厚度为1 mm和直径10 mm的紫铜作为波形整形器。SHPB试验的花岗岩试件来自鞍山齐大山铁矿，密度2 565 kg/m3，弹性模量41.86 GPa，泊松比0.19，抗压强度54.93 MPa，抗拉强度5.35 MPa。

2 数值模型

2.1 平面黏结接触模型

关于颗粒流模型的接触模型选择方面，通常使用接触黏结模型和平行黏结模型，平行黏结模型可实现颗粒的转动，相比于接触黏结模型模拟岩石的相关问题更为合理，但二者的压拉强度比范围较小、内摩擦角过低和强度包络线是线性的[11]，致使其又具有一定的局限性，更适合用于模拟低强度的脆性岩石[12]。

Potyondy[10]提出的平面黏结接触模型克服了压拉强度比范围较小的局限性，适用于模拟硬质岩石，图2为典型平面黏结接触模型示意图，平面黏结接触模型模拟2个抽象表面之间的交界面的行为，这些抽象表面与相应的颗粒紧密相连，因此抽象表面定义了每一个颗粒的有效表面，颗粒可以被认为是“有边缘的颗粒”，当接触面破坏以后“有边缘的颗粒”能提供连锁阻力和旋转阻力。破坏后的颗粒的旋转被抑制，使得压拉强度比显著增大。相比接触黏结模型和平行黏结模型压拉强度比在3～7，平面节理模型的压拉强度比可达到24[6]，所选花岗岩试样的压拉强度比是13.1。

2.2 接触处的相对运动方程

2.3 力—位移定律

2.4 细观参数校准

岩石材料的细观参数标定需通过多次设置颗粒流模型的细观参数进行单轴加载试验和误差检验的数值试验，根据试验结果获得与材料宏观性质（图1）大体一致的细观参数[13]。对于固定的岩石颗粒流模型的颗粒尺寸分布，本研究最大与最小颗粒粒径比为1.66，模型所含颗粒数为21 610个，调试细观力学参数的一般步骤[14]为：调试细观参数杨氏模量来确定宏观参数弹性模量；调试细观参数刚度比确定宏观参数泊松比；调试细观参数黏结接触力确定宏观参数抗压强度和抗拉强度；霍普金森压杆冲击试验检验和修正，此次校准的应变率是196.26 s-1，花岗岩宏观参数和数值模拟见表1，模型细观力学参数的选取见表2，单轴试验数值计算结果见图3。

?

?

3 数值计算

动态抗压强度（动态抗拉强度）与静态抗压强度（静态抗拉强度）的比值为动态增长因子，应变率在10-3～102s-1范围内，动态增长因子随着应变率的递增而缓慢增加，动态增长因子大约在1～1.5范围内；应变率在102～104s-1范围内变化，动态增长因子随着应变率的递增而快速增加，动态增长因子大约1.5～6范围内变化，说明了应变率与动态增长因子之间存在一定关系[15-17]。相关研究介绍了动载荷作用下巴西圆盘起裂模式和最终破坏形式，周子龙等[18]分析收集的动态耦合荷载加载巴西圆盘试样的破碎碎块，认为其是3种破裂模式，直径劈裂式、伴有楔形粉碎区的中心裂纹式、伴有压碎带的破坏式。周子龙等[19]结合高速摄影图像，认为其是2种破裂模式，直径劈裂式和伴有压碎带的破坏式。Mehrdad Imani等根据不同应变率范围内裂纹起裂位置，介绍了4种巴西圆盘破裂模式，如图4所示。

小于150 s-1属于拉伸劈裂破坏模式（图4（a）是应变率小于50 s-1，裂纹由试件中心起裂，随着载荷增加裂纹向试件两端扩展，最终试件破裂成两半；图4（b）是应变率小于150 s-1，更多的裂纹在试件中心产生，随着载荷增加裂纹向试件两端扩展，最终试件破裂成多半），150～600 s-1属于分支破坏模式，大于600 s-1属于粉碎破坏模式。因此，本文在10-1～104s-1分析平台巴西圆盘（为了避免加载点应力集中现象，岩石试件的平台角20°，直径50 mm）的起裂模式和最终破裂模式。图5为应变率分别是5 s-1、125 s-1、500 s-1、1 000 s-1的裂纹发展方式。

图5（a）和图5（b）显示了拉伸型微破裂裂纹占主导，图5（c）和图5（d）显示了剪切型微破裂裂纹和拉伸型微破裂裂纹共同占主导作用，图 5（e）、图 5（f）、图5（g）、图5（h）显示了裂纹起裂、扩展位置随计算时步递增的发展历程，并结合拉伸型微破裂裂纹和剪切型微破裂裂纹最终产生形式，随着应变率的增加，剪切型微破裂裂纹数量不断增加。图5（e）和图5（f）显示了微破裂裂纹是由中心起裂的，图5（f）中心处产生的裂纹数量多于图5（e），即分别是中心规律微破裂裂纹分布的中心起裂的拉伸破裂和中心随机微破裂裂纹分布的起裂的破裂，与图4（a）和图4（b）描述的破裂模式一致，其破裂方式是试样中心起裂破裂成完整2块和试样中心起裂破裂成完整2块伴随碎屑。图5（c）是加载端压碎区共生拉伸型微破裂裂纹和剪切型微破裂裂纹，微破裂裂纹由压碎区前端呈放射性向试样中心发展，首先在试样中心区域产生较多随机分布拉伸型微破裂裂纹，随后剪切型微破裂伴生于拉伸型微破裂裂纹，连接拉伸型微破裂裂纹；图5（d）是加载端压碎区共生拉伸型微破裂裂纹和剪切型微破裂裂纹，共生微破裂裂纹由压碎区前端呈放射性向试样中心发展，最终导致试件破裂，结合图5（g）和图5（h），与图4（c）和图4（d）描述的破裂模式一致，其破裂方式是伴生压碎区的复合型微破裂裂纹共同作用的破裂。

图6为不同应变率条件下试样内诱发微破裂裂纹的数量情况，剪切型微破裂裂纹数目随着应变率的增加而增加，在应变率0.1～1 000 s-1范围内递增过程中，剪切型微破裂裂纹数占拉伸型微破裂裂纹数的百分比是3.42%～78.38%，其中在应变率在150 s-1和600 s-1相应的剪切型微破裂裂纹数占拉伸型微破裂裂纹数的百分比分别是21.68%和55.10%。由图对比还可以看出，随着应变率的增长，剪切型微破裂裂缝数目的增量趋势明显大于拉伸型微破裂裂纹数目的增加，最终趋于与拉伸型微破裂裂纹相同数量级的裂纹数，已有关于应变率与动态应力之间的变化关系研究，证明动态应力随应变率的增加而增加，由此可说明，剪切微破裂裂纹数目的增加可能是导致试件动态强度相应增加的原因[6]。

为分析剪切型微破裂裂纹数量变化与平台巴西圆盘试件应力之间的关系，对平台巴西圆盘动态试验模拟过程进行微破裂裂纹数量及应力数值实时监测，结果如图7所示。其中，图7（a）显示了第一次应力峰值（8.80 MPa）大于第二次应力峰值（7.76 MPa），拉伸型微破裂裂纹在第一次峰值前67.09%的应力峰值处，剪切型裂纹在达到第二次应力峰值前112.18%的第二次应力峰值处，说明此时的剪切型微破裂裂纹产生在第一应力峰值的下降阶段处，即98.87%的第一次应力峰值处。图7（b）显示了第二次应力峰值（8.57 MPa）略大于第一次应力峰值（8.49 MPa），拉伸型微破裂裂纹在第一次峰值前69.55%的应力峰值处，剪切型微破裂裂纹在达到第二次应力峰值前88.94%的第二次应力应力峰值处，即89.81%的第一次应力峰值处，计算时步显示了剪切型微破裂裂纹产生在第一次应力峰值之后的下降阶段。

图7（a）和图7（b）分别展现了2种应力与曲线类型，应力首先随着计算时步的增加呈线性增加到第一次峰值应力，应力达到第一次峰值后开始跌落，此时，拉伸型微破裂裂纹数量呈现突变现象，说明由试件中心起裂的微破裂裂纹向试件两端开始扩展。因为试件还未完全破裂致使其还具有一定的承载能力，所以应力跌落到一定程度后开始继续上升至第二次应力峰值，此时，在应力达到第二次应力峰值之前，试件加载端部产生剪切型微破裂裂纹，达到第二次应力峰值之后，剪切型微破裂裂纹数量趋于稳定。应力达到第二次应力后开始逐渐降低直至试件破裂，此阶段产生次生微破裂裂纹。

图 7（c）、图 7（d）、图 7（e）、图 7（f）、图 7（g）显示了拉伸型微破裂裂纹分别产生在应力达到应力峰值前60.82%、50.43%、21.74%、28.42%、42.18%的应力峰值处，剪切型微破裂裂纹分别产生在应力达到应力峰值前93.80%、67.75%、57.29%、32.75%、43.61%的应力峰值处，图7（h）、图7（i）、图7（j）显示拉伸型微破裂裂纹和剪切型微破裂裂纹同时产生，分别在达到应力峰值前45.83%、47.68%、47.67%的应力峰值处，图7（h）起始阶段拉伸型微破裂裂纹数量与剪切型微破裂裂纹数量基本一致，随后拉伸型微破裂裂纹数量多于剪切型微破裂裂纹数量，剪切型微破裂裂纹数量随计算时步呈现线性增长阶段，非线性增长阶段之后趋于平缓。图7（i）、图7（j）显示剪切型微破裂裂纹数量大于拉伸型微破裂裂纹数量，图7（j）更加明显呈现这一特性，随后拉伸型微破裂裂纹数量迅速增长超越剪切型微破裂裂纹数量，分析这一原因是应变率较大导致加载端部产生压破裂，远离加载端部位置主要是拉伸型微破裂裂纹占主要地位，当应变率增长到一定程度之后，试件主要受压破坏产生主裂纹，次生裂纹由拉伸型微破裂裂纹产生。

图7（c）～图7（j）展示了同一种应力曲线类型，应力随着计算时步增加至应力峰值后逐渐降低直至试件破裂，随着应变率的增加，剪切型微破裂裂纹的产生在由应力达到应力峰值之后转变为产生在应力达到应力峰值之前，即当在低应变率条件下，剪切微破裂裂纹数量在应力峰值下降阶段呈现突变现象，在中应变率条件下，剪切微破裂裂纹数量接近于应力峰值点处呈现突增现象，在高应变率条件下，剪切微破裂裂纹数量在应力—时步曲线的上升阶段呈现突增现象，且这一时步点随着应变率的递增而越早产生。

4 结论

利用平面黏结接触模型模拟花岗岩平台巴西圆盘试样在高应变率加载情况下的动态响应和破坏机制，从细观角度深入探究花岗岩平台巴西圆盘试样的动态破裂模式，其研究结论如下：

（1）随着应变率的增长，花岗岩平台巴西圆盘的破裂模式主要分为2种破裂模式，微破裂裂纹由中心起裂的破裂模式和伴生压碎区的复合微破裂裂纹共同作用的破裂模式，其中，中心起裂破裂模式还可以细分为中心规律微破裂裂纹分布的中心起裂的拉伸破裂和中心随机微破裂裂纹分布的起裂的破裂。当应变率较高时，已不满足拉伸试验基本原理，但是对于高应变率的工程问题具有一定指导意义。

（2）剪切型微破裂裂纹数目随着应变率的增加而增加，剪切型微破裂裂纹数目的增量趋势明显大于拉伸型微破裂裂纹数目的增加，剪切型微破裂裂纹数占拉伸型微破裂裂纹数的百分比是随应变率的增长呈显著递增趋势。

（3）在低应变率情况时，应力与计算时步曲线存在2个应力峰值，即峰后强化和峰后弱化，剪切型微破裂裂纹产生在第一应力峰值的下降阶段处；在中高应变率情况时，应力与计算时步曲线存在一个应力峰值，即单峰值，随着应变率的增加，剪切型微破裂裂纹的产生在由应力达到应力峰值之后转变为产生在应力达到应力峰值之前。