无论是静态还是动态，岩石的抗拉强度远低于岩石的抗压强度，二者之间通常相差了一个数量级^[1]，正是因为岩石材料的这种性质，导致其大概率会发生拉伸破坏。由于岩石试件的制作和试样实现单轴拉伸加载比较困难，致使直接测量岩石拉伸强度较为不易。因此，国内外普遍使用巴西圆盘劈裂法^[2]来间接获得岩石抗拉强度。

对于大理岩材料，众多学者采用试验或数值模拟的方法研究了其动态力学行为。在试验方面，如唐礼忠等^[3]、李伟等^[4]通过分离式霍普金森压杆（SHPB）进行了动态劈裂试验，研究结果表明在一定的应变率范围内，拉伸强度对应变率的敏感性要高于压缩强度，且材料的拉伸强度有着明显的应变率效应。王伟祥等^[5]利用扫描电镜对受动载作用的岩样断口形貌进行细观观察，分析了大理岩在动态劈裂下的断口形貌特征。宋小林等^[6]通过记录试样动态劈裂的应变路程，判断大理岩圆盘的起裂位置和时刻，以及裂纹的扩展顺序和速度。Wong等^[7]在Carrara大理石的动态巴西圆盘试验中借助高速摄像仪，观察到了在动态加载下，损伤的演化速度和范围普遍随径向压缩应变率的增大而增大。在数值模拟方面，韩东波等^[8]采用ANSYS/LS-DYNA有限元的方法对大理岩动态劈裂拉伸试验过程进行数值模拟，研究了不同形状的冲头对试验结果的影响。Lou等^[9]基于巴西圆盘劈裂试验，利用LS-DYNA软件确定了一套适用于高强脆性大理石的HJC模型参数。方新宇等^[10]利用LS-DYNA模拟SHPB得出系统破碎能与系统耗散能演化关系相似。

综上所述，国内外学者围绕大理岩动态劈裂行为开展了大量研究，但以往工作多是集中在试样的尺寸、应变率效应以及弹速影响等方面。对动态劈裂条件下岩石的裂纹扩展及各阶段的能量演化机制研究开展较少。本文在动态劈裂试验的基础上，先结合相应数值模拟，标定出锦屏大理岩试样合适的本构模型参数，随后从大理岩动态劈裂裂纹扩展与能量演化等方面进行多角度分析，力求得出有实际参考价值的结论。

1.   大理岩动态劈裂试验

1.1   岩样制备

本文的大理岩试样来自锦屏II级水电站深埋引水隧洞，所有试样均取自同一块完整的岩石，试样长径比为0.5，即Ф50 mm×25 mm，试样两端面不平整度误差均小于0.02 mm。本次试验安排3组不同的加载弹速（5.0，6.7，9.0 m/s），为提高数据的准确性，每组准备3～5个试样，试样的基本物理力学参数见表1。

表  1  大理岩基本物理力学参数

Table  1.  Basic physical and mechanical parameters of marble

参数	ρ0/（kg·m−3）	E/GPa	G/GPa	μ	fc/MPa
取值	2 763	33.32	12.43	0.34	120
注：ρ0为密度；E为弹性模量；G为剪切模量；μ为泊松比；f为抗压强度。

下载: 导出CSV 

| 显示表格

1.2   测试方案和原理

图1为SHPB试验装置，该系统主要由高压气腔、子弹、入射杆，试样、透射杆、缓冲装置以及相配套的计算机数据处理系统组成。其中，入射杆和透射杆以及子弹的材料密度均为7 900 kg/m³，弹性模量为210 GPa，泊松比为0.3。子弹、入射杆、透射杆和吸收杆的长度分别为400，2 400，1 400，1 000 mm，此外，在入射杆前端放置一个脉冲整形器，以保证试验过程中应力均匀与恒应变率加载，达到提高试验精度的目的^[11]。试样的基本物理参数及弹速见表2。

图  1  SHPB装置示意图

Figure  1.  Schematic diagram of SHPB device

下载: 全尺寸图片 幻灯片

表  2  试样几何尺寸和对应的子弹速度

Table  2.  Sample geometric dimension and corresponding striking velocity

编号	M/g	l/mm	d/mm	ρ0/（kg·m−3）	v/（m·s−1）
C14	134.5	24.99	49.65	2 779.89	5.11
C15	137.6	25.38	49.72	2 793.13	4.93
C21	138.2	25.40	49.70	2 806.09	5.06
C1	136.8	25.10	49.70	2 810.12	6.87
C2	133.9	25.10	49.68	2 753.49	6.64
C18	138.6	25.54	49.77	2 789.81	6.66
C8	136.1	25.60	49.72	2 738.21	8.77
C9	130.0	25.13	49.66	2 671.19	9.30
C11	138.3	25.40	49.72	2 805.13	8.88
注：M为质量；l为高度；d为直径；ρ0为密度；v为弹速。

下载: 导出CSV 

| 显示表格

当高速子弹撞击入射杆瞬间，在入射杆前端会生成入射波，入射波在入射杆方向上进行一维传播，当入射波传至试样界面时，一部分透过试样以透射波的形式进入透射杆，而另外一部分被反射回来，在入射杆上形成反射波。因此，基于一维假设和应力均匀性假设，试样应变（ε_s）、应变率（$\dot{\varepsilon_{\mathrm{s}}}$）和应力（σ）之间的关系可以采用“三波法”计算得到^[12]：

式中：C₀——压杆波速/（m·s⁻¹）；

E——压杆的弹性模量/GPa；

A₀——压杆的截面积/m²；

A_s——试样的截面积/m²；

ε_i(t)、ε_r(t)、ε_t(t)——入射波、反射波与透射波的 应变时程。

2.   试验结果分析

2.1   应力平衡分析

图2是大理岩试样典型（编号C2）的动态波形平衡图。可以看出入射波与反射波波形升降平缓，未出现明显振荡。此外，可见在波峰前，入射波与反射波之和与透射波几乎相等，说明试样两端在加载前已经基本达到应力平衡，满足应力均匀性假设，即本文试验结果是真实有效的。

图  2  动态应力均衡图

Figure  2.  Diagram of dynamic stress equilibrium

下载: 全尺寸图片 幻灯片

2.2   动态劈裂强度与应变率

在SHPB测试中，使试样沿直径方向劈裂破坏，可间接测定岩石抗拉强度（也称劈裂强度）。图3显示了锦屏大理岩劈裂强度随着应变率的变化关系，作为对比，引用了唐礼忠等^[3]、李伟等^[4]、Padmanabha等^[13]、Zou等^[14]的数据。由于试样内部存在微孔洞、裂隙等缺陷^{[15 − 16]}，从而导致了拟合曲线出现一定离散，但在应变率5～35 s⁻¹范围内，两者整体上呈线性正相关。一般而言，应变率敏感度与岩石内节理发育程度相关——节理越发育，其应变率敏感性相对较低^[17]。本文取自锦屏深埋引水隧洞大理岩试样的节理发育可能相比于比其他地区的更为成熟，从而导致其应变率敏感度较低。因此，在该地区进行爆破等（高应变率条件）施工作业时，工程岩体更易产生劈裂拉伸破坏，进而引发工程地质灾害，应当注意适时支护。

图  3  劈裂强度与应变率关系比较

Figure  3.  Comparison between splitting strength and strain rate

下载: 全尺寸图片 幻灯片

2.3   动态劈裂宏观破坏

图4给出了试样代表性破坏形态图。可以看出，在试验的弹速范围内，沿着径向有一条贯穿整个试样的主裂纹。试样两端由于接触面积小、冲击力作用迅速，致使接触端面发生应力集中，导致在主裂纹贯通之前，试样两端的应力就已经超过材料的强度，从而使接触面迅速萌生裂纹并且快速发展，最终呈现一个类似于三角形粉碎区的破坏形态^[18]。随着弹速的增大，粉碎区面积也随之增加。一般情况下，入射杆端的粉碎区面积略大于透射杆端的粉碎区面积。

图  4  典型弹速下试样破坏比较

Figure  4.  Comparison of sample failure modes at typical velocities

下载: 全尺寸图片 幻灯片

3.   损伤破裂模式数值模拟分析

3.1   模型与本构的确定

为了进一步探究动态劈裂过程中大理岩的裂纹扩展过程，本部分借助ANSYS/LS-DYNA数值软件进行分析。如图5所示，模型尺寸与实物1∶1建立，压杆和试样选用3D-Solid164单元进行网格划分。入射杆和透射杆轴向方向网格分别划分为240段和140段，径向方向划分为48段，试样轴向和径向方向上分别划分50段和120段，本次模拟选择映射网格划分。试样与压杆间采用面面侵蚀接触，由于两者刚度差距较大，经过反复试算，最终确定最佳接触刚度比例因子为0.04～0.06。

图  5  动态模拟中有限元网格划分

Figure  5.  Finite element meshing for dynamic simulation

下载: 全尺寸图片 幻灯片

压杆采用线弹性模型，大理岩试样则采用Cowper-Symonds本构模型^[19]，该模型可以提供材料的塑性变形模量，并考虑动态加载下的应变率效应，可以较好模拟出混凝土等材料的劈裂过程^[20]。在软件中其关键字为MAT_PLASTIC_KINEMATIC，本构方程见式（4）：

式中：${\sigma _{\text{d}}}$——动态屈服应力/MPa；

${\sigma _{0}}$——静态屈服应力/MPa；

E_p——塑性硬化模量/GPa；

$\stackrel{.}{\varepsilon }$——应变率/s⁻¹；

$\overline{\varepsilon }$——等效塑性应变；

C、p——应变率常数。

本构方程中应变率常数C和p的取值可以依据文献[21]中锦屏大理岩动态压缩试验数据进行拟合所得（图6），最终C取162.19，p取1.25。

图  6  Cowper-symonds本构中参数拟合曲线

Figure  6.  Parametric fitting curves of the Cowper-symonds model

下载: 全尺寸图片 幻灯片

考虑到试样的裂纹扩展拉伸破坏效果，建模中通过关键字MAT_ADD_EROSION来定义失效准则。最终得到大理岩Cowper-Symonds本构模型参数，见表3。

表  3  Cowper-Symonds本构参数

Table  3.  Constitutive parameters of the Cowper-Symonds model

参数	ρ0/（kg·m−3）	E/GPa	μ	C	p	β
取值	2 763	33.32	0.34	162.19	1.25	0.1
注：ρ0为密度；E为弹性模量；μ为泊松比；C、p为应变率常数；β为硬化参数。

下载: 导出CSV 

| 显示表格

图7是模拟的入射波波形与试验中的入射波波形对比。从图中可以看出被赋予一定速度的子弹冲击入射杆产生的波形与试验中入射波波形基本吻合，能实现与试验相同的冲击荷载，说明子弹以及杆件建模的合理性。

图  7  入射波形对比

Figure  7.  Comparison of incident waveforms

下载: 全尺寸图片 幻灯片

3.2   模拟结果破坏过程分析

图8和图9分别展示了数值模拟中试样裂纹扩展过程以及不同弹速下的试样最终破坏形态，本次模拟的3种弹速情况下裂纹扩展趋势差异不大，故取低弹速（v=5.0 m/s）下的试样裂纹扩展情况进行概括分析。模拟结果显示，在试样应力未达到平衡之前，靠近入射端（下端）且偏向中心位置的应力会率先达到材料的抗拉强度，所以试样会在偏向中心靠近入射端的方向起裂。试样裂纹扩展的顺序为：偏心起裂→双向扩展→入射端产生裂纹→透射端产生裂纹→完全破坏。整个过程仅仅持续了几十微秒，试样最终呈现出显著的脆性拉伸破坏。这与文[6]中的试验结果以及文[7]中的高速摄影结果相吻合，并且最终破坏形态与图4展示的试验结果相似，进一步验证了本构模型的合理性以及模型参数取值的准确性。

图  8  不同时刻裂纹的扩展过程 （v=5.0 m/s）

Figure  8.  Crack development process at different times

下载: 全尺寸图片 幻灯片

图  9  不同弹度下裂纹扩展最终形态

Figure  9.  Final cracking shapes under different velocities

下载: 全尺寸图片 幻灯片

3.3   模拟结果能量分析

试样破坏过程往往伴随着能量的吸收与释放，当子弹以不同速度撞击入射杆时，子弹的动能会以其他能量的形式输入到整个系统中。LS-DYNA后处理软件LS-PREPOST中可以利用软件内部的能量处理手段自动跟踪能量转化过程，能够更直观地理解整个实验。本次模拟以入射杆-试样-透射杆为一个系统，所涉及到的能量如下^{[22 − 23]}：

①系统外部的总能量E_t，其包括系统中所有微元物理能量的积累，公式如下：

式中：Q——最后一个时间步的外部总能量/J；

N_p——微元数/个；

m——微元质量/kg；

g_i——重力加速度/（m·s⁻²）；

F_i——外部载荷/N；

ΔU_i——位移增量；

M_j——外部力矩/（N·m）；

Δθ_j——角度增量/（°）。

②加载期间系统的动能E_k，其定义为系统中所有微元平移和旋转所产生的能量，计算式如下：

式中：m_i——单元粒子的广义质量/kg；

ν_i——单元粒子的广义速度/（m·s⁻¹）。

③加载期间系统所储备的内能E_i，主要以弹性应变能、塑性应变能的形式存在。

为了减少杆件之间的摩擦，在试样与杆件接触面处涂抹凡士林，故系统中因摩擦而消耗的能量可以忽略不计。假设系统内部与外部没有热量的交换，试样耗散能E_a主要用于试样裂纹的扩展，其计算式为：

图10是SHPB试验过程中的能量变化。在子弹撞击之初，冲击应力波迅速施加到入射杆端，系统能量持续增加，一直增加到A时刻，系统总能量累积至最大，随后保持不变，此时系统总能量等于系统内能与系统动能之和。冲击波在B时刻抵达试样，一部分以反射拉伸波的形式返回入射杆，另一部分则以透射压缩波的形式继续传输，在AB时间间隔内，系统的内能几乎与动能相等。系统的内能与动能在B时刻分离，此时岩样中微元的应力超过其本身的强度，内部的微裂纹被激活。接着，系统储存的内能开始释放，系统耗散能增加，裂纹开始扩展。直至C时刻，试样完全破坏，系统动能和耗散能达到最大，内能处于最小。随后，动能减小，内能增加，耗散能则保持不变，整个系统的能量突变时刻对应着图8中的裂纹扩展时刻。可见，试样的破坏过程伴随着系统能量间的相互转换，也印证了数值模拟中用能量变化来解释试样的破坏规律是可行的。图11是不同弹速下系统内能和动能的关系。弹速的增大导致子弹冲击波的强度随之增加，系统内能和动能的数值也相应增大。与此同时，系统内部能量的转化速率加快，从而加速试样的破坏。

图  10  系统能量时程曲线

Figure  10.  System energy with time

下载: 全尺寸图片 幻灯片

图  11  系统内能与动能时程曲线

Figure  11.  System energy and kinetic energy with time

下载: 全尺寸图片 幻灯片

4.   结论

（1） 冲击荷载下大理岩动态劈裂破坏表现为明显的拉伸形式的脆性破坏，其劈裂强度具有明显的应变率效应。在应变率5～35 s⁻¹范围内，劈裂强度与应变率呈线性正相关。相比较其他地区的大理岩，锦屏大理岩劈裂强度的应变率敏感性略低。

（2） 考虑应变率效应的Cowper-Symonds本构能够较好地表征锦屏大理岩动态劈裂拉伸破坏形态，其模拟结果与试验观测现象比较吻合，表明了文中本构参数选用的合理性，可用来模拟分析锦屏水电站深部洞室围岩损伤破坏乃至失稳的全过程。

（3） 系统总能量维持稳定后，如出现系统内能急剧下降现象，则表明试样开始发生破坏；当试样完全破坏时，系统内能降至最低。系统能量的变化能较好解释试样从裂纹开始萌生到完全破坏的过程；随着弹速的增加，系统内能和动能随之增大，从而加速试样的破坏。