岩石材料的力学性能及变形特征对于岩土、采矿、隧道等工程施工设计至关重要。常见的岩石材料是由不同矿物成分构成的非均质、不连续介质^[1]。从微观角度看，岩石材料的力学性能及变形特征是矿物晶质内与晶质间断裂、滑移等力学行为累积效应的宏观表达^[2]。因此，从矿物晶质尺度描述外部载荷下岩石断裂损伤的细观特征和演化过程，对于揭示岩石的损伤破坏机制具有重要的现实意义。

数值模拟是从细观尺度揭示岩石断裂损伤破坏的有效手段。李晓锋等^[3]在考虑岩石晶质分布的基础上，利用颗粒流程序（particle flow code，PFC）对冲击载荷下多种岩石材料的力学性能及能量耗散规律进行了模拟。Zhou等^[4]通过PFC^2D对某花岗岩的动态冲击损伤过程进行了反演，深入分析了高应变率下岩石内部微裂纹的演化规律。赵翰卿等^[5]利用PFC对陶瓷巴西圆盘的分离式霍普金森压杆（split Hopkinson pressure bar，SHPB）实验过程进行了模拟分析，得出的模拟破坏模式与实验结果基本一致。基于PFC的模拟方法在岩石材料的动、静态力学性能模拟及岩石断裂破坏特征表征等方面展现出良好的应用前景。然而，上述研究采用均质的颗粒结构以及单一的晶质种类，在表征晶质岩石的非均质性、细观结构的差异性以及力学性能的各向异性等方面有所欠缺。

Potyondy^[6]为研究晶质结构岩石的细观结构特征对其力学性能的影响，提出了基于Voronoi结构的晶质岩石模型（bonded particle model，BPM）构建方法，并对该模型在多种实验力学中的断裂损伤细观发展过程及裂纹演化过程进行模拟分析，得到了与实验结果高度吻合的模拟结果。早期的PFC-等效晶质模型（grain based model，GBM）利用Voronoi空间结构及PFC颗粒集合体表征晶质材料的细观结构特征和力学特性，为从细观晶质尺度解决岩石力学问题提供了有效的模拟手段^[7–8]。随着实验手段的不断丰富，PFC-GBM逐渐成为与CT电子扫描技术^[9]、电镜探测技术^[10]等微观检测手段以及“多参数Weibull分布”^[11–12]等数理统计方法相结合的，能够良好地表征材料非均质性及细观结构特征的数值模拟方法。Saadat等^[13–15]基于PFC-GBM方法，对晶质岩石的单轴抗压、直接拉伸以及巴西劈裂等实验进行了大量模拟分析，得到了与室内实验相近的力学性能和破坏特征。Liu等^[16]利用GBM对非均质岩石的晶质尺寸效应和细观破坏特征进行了研究，发现晶间裂纹与晶内裂纹的贯通是导致岩石宏观破坏的主要原因。张涛等^[17]基于PFC-GBM方法构建了晶质岩石模型，研究了岩石非均质性及晶质尺寸对其动态力学性能的影响，验证了利用PFC-GBM从矿物晶质尺度研究岩石动力学性能的合理性。

综上所述，基于PFC-GBM的岩石细观结构能够有效地表达非均质岩石材料的晶质结构特征，但是存在对岩石非均质性以及矿物晶内及晶间断裂损伤模式和破坏机理表征不足、模型构建过程复杂、计算效率较低等问题，亟待优化和解决。目前，基于PFC^2D-FLAC^2D耦合方法和PFC^2D-GBM方法的非均质岩石断裂损伤细观特征模拟研究较为鲜见，结合二者的优势，本研究通过数值模拟对矿物晶间与晶内裂纹进行细化分组和量化统计，剖析非均质岩石动态断裂损伤的细观发展过程，为从矿物晶质尺度揭示非均质岩石断裂损伤机制提供合理、有效的分析方法。

1.   基于PFC-GBM的单轴破坏细观特征模拟

1.1   PFC-GBM构建

本研究的岩石试样及相关物理力学参数均取自赵奎等^[18]的实验数据，首先使用PFC-GBM构建非均质岩石模型，结合实验数据，确定模型的细观参数。赵奎等^[18]选用的花岗岩材料中，长石、石英、云母的质量分数分别约为60%、35%、5%。设矿物晶质尺寸在1～2 mm内均匀分布。PFC-GBM模型构建步骤^[19]：(1) 首先在宽50 mm、长100 mm的范围内生成如图1(a)所示的基本球颗粒，共572个；(2) 利用球颗粒的空间位置进行Voronoi结构细分，生成如图1(b)所示的多边形集合区域；(3)删除基本球颗粒，生成粒径在0.400～0.664 mm范围内均匀分布的颗粒集合体，将模型伺服至平衡状态后，通过Voronoi结构将颗粒按所在区域进行分组，如图1(c)所示；(4) 按各矿物的质量分数在各组颗粒间随机选取物料，直至模型中各矿物的质量分数与实际一致，最终生成如图1(d)所示的PFC-GBM模型。

图  1  PFC-GBM的构建过程

Figure  1.  Construction process of PFC-GBM

下载: 全尺寸图片 幻灯片

为模拟材料内部的胶结状态，对模型中晶内与晶间颗粒添加黏结模型。如图2所示，有石英晶内黏结、云母晶内黏结、长石晶内黏结和晶质边界黏结4种黏结。为简化分析，对晶内及晶间黏结均添加Flatjoint黏结模型^[20]。Flatjoint黏结模型可简化为由一组具有恒定法向和切向刚度的“弹簧”元件以及具有一定抗拉和抗压能力的“黏壶”元件构成。当黏结键被激活时，颗粒间的张拉强度由法向刚度k_n及法向黏结强度S_n定义，压剪强度由切向刚度k_s、切向黏结强度S_s和内摩擦角Ф定义；当外部载荷超过黏结键张拉或压剪极限后，黏结键断裂，颗粒恢复至非黏结状态，此时颗粒间无抗拉强度，抗压能力由接触间隙g_s和法向刚度k_n定义，抗剪能力由切向刚度k_s和摩擦系数f定义。

图  2  黏结分布及黏结原理

Figure  2.  Bond distribution and bond principle

下载: 全尺寸图片 幻灯片

1.2   晶质岩石单轴压缩模拟

确定合理的细观参数是PFC模拟材料力学性能及变形特征的重要前提。利用无侧限压的单轴压缩模拟实验，通过试错法^[21]不断调整细观参数取值，获取与实验结果相近的宏观力学参数。在二维平面中，对于50 mm×100 mm的标准岩石试件，设置轴向加载速率为0.005 m/s，保证实验过程中准静态力学行为，直至试样破坏。图3显示了应力-应变曲线及试样破坏模式的模拟与实验结果^[16]对比。标定获取的细观参数见表1。表2列出了试样的宏观物理力学参数模拟结果和实验结果，其中：K₁为线性接触刚度比，E_f为线性接触模量，K₂为Flatjoint黏结刚度比，E为弹性模量，μ为泊松比，σ为单轴抗压强度。可以看出，模拟得到的宏观力学参数与实验结果极为接近，说明所选取的PFC细观参数合理。

图  3  试样应力-应变曲线及破坏模式的数值模拟与实验结果^[18]对比

Figure  3.  Comparison of numerical simulation and experimental results^[18] of stress-strain curves and failure modes

下载: 全尺寸图片 幻灯片

表  1  试样的细观参数

Table  1.  Microscopic parameters of the specimen

Material	f		Ф/(°)	K1	Ef/GPa	K2	Sn/MPa	Ss/MPa
Quartz	0.20		30	2.3	53	2.3	34.0	42.0
Feldspar	0.15		33	2.5	43	2.5	27.2	33.6
Mica	0.18		36	3.8	34	3.8	20.4	25.2
Grain boundary	0.50		42	5.0	22	5.0	5.1	6.3

下载: 导出CSV 

| 显示表格

表  2  试样的主要力学参数

Table  2.  Main mechanical parameters of the rock specimen

Method	E/GPa	μ	σ/MPa
Simulation	30.95	0.205	125.72
Experiment	30.58	0.210	126.57
Error/%	1.2	2.4	0.7

下载: 导出CSV 

| 显示表格

PFC模型通过颗粒伺服^[22]的方法使颗粒集合体内部应力获得最大程度的释放，同时保证其具有较高的密实度，但这将导致PFC颗粒集合体在赋予黏结后具有较低的初始空隙率，难以展现受压载荷时的压密现象。从轴向应力-轴向应变曲线可以看出，从开始加载至起裂时刻，模拟得到的应力与应变表现出良好的线性关系，此阶段可认为是模型在单轴载荷下的弹性响应阶段。图4给出了晶质岩石在单轴压缩过程中的微裂纹演化过程。可以看出，晶质岩石的单轴压缩破坏过程主要分为初始状态、起裂、微裂纹稳定发育、裂纹贯通4个阶段。从初始状态至微裂纹稳定发育阶段，模型内部的损伤主要以晶间裂纹为主；在裂纹贯通阶段，晶间与晶内均出现较高程度的损伤，从而使岩石表现出宏观破坏。虽然不同时刻微裂纹在岩石内部的分布情况在一定程度上反映了晶质岩石单轴破坏的细观特征，但是仍然不足以定量地描述微裂纹种类、数量与加载破坏进度的关系。

图  4  非均质岩石单轴破坏的微裂纹演化过程

Figure  4.  Microcrack evolution of heterogeneous rock under uniaxial failure condition

下载: 全尺寸图片 幻灯片

为了更清晰地了解晶质岩石单轴破坏过程中的断裂损伤细观特征，通过自编Fish语言，统计得到各矿物晶内及晶间的微裂纹数随轴向应变的演化曲线，如图5所示。继续加载时，起裂后试样内部微裂纹稳定发育，微裂纹数逐渐增加，损伤破坏程度明显升高。Wong等^[23]基于GBM模型对晶质岩石断裂破坏的细观特征进行了深入研究，结果表明，晶质间的张拉破坏和晶间裂纹的逐步贯通是引起岩石最终发生屈服破坏的主要原因。王桂林等^[24]采用GBM模型对花岗岩细观破坏过程中的晶间裂纹与穿晶裂纹行为特点进行了模拟研究，结果显示，单轴条件下的模拟结果与实验结果基本吻合，且穿晶裂纹数明显低于晶界裂纹数。本研究的模拟结果与文献[19–20]的结论基本一致，说明模拟结果合理。

图  5  单轴破坏中的微裂纹演化过程

Figure  5.  Microcrack evolution in uniaxial failure

下载: 全尺寸图片 幻灯片

从图6中微裂纹细化分组情况可以看出，在晶质岩石的单轴破坏过程中，晶间裂纹数明显高于晶内裂纹数，说明晶间黏结破坏是晶质岩石断裂损伤破坏的主要原因。随着加载的进行，晶间裂纹数由多到少分别为长石晶内裂纹、石英晶内裂纹和云母晶内裂纹，与各矿物成分的质量分数对应。单轴载荷下，岩石内部各部分受力均匀，各类矿物晶质的含量与分布是岩石局部断裂破坏呈各向异性的重要原因。在破坏后期，石英晶内裂纹数超过长石晶内裂纹数。周喻等^[25]基于等效晶质模型对实验过程中晶质岩石断裂损伤细观特征进行了研究，也得出了类似的结果。因此，利用PFC-GBM从矿物晶质尺度研究非均质晶质岩石的断裂损伤过程及力学性能具有一定的合理性和优势。

图  6  微裂纹统计及细化分组

Figure  6.  Microcracks statistics and subdivision

下载: 全尺寸图片 幻灯片

2.   晶质岩石动态冲击破坏细观特征

2.1   耦合SHPB系统及等效载荷施加

SHPB是研究材料动力学特性最常用的实验方法之一。利用数值模拟方法模拟SHPB实验全过程，对于快速合理预测材料动力学性能、降低实验成本具有重要意义。为研究非均质岩石的动态破坏过程，构建如图7所示的二维耦合SHPB系统。利用PFC-GBM构建岩石试样，细观参数见表2。在FLAC^2D中构建杆件，杆件的细观参数如下：密度为7.89 g/cm³，入射杆和透射杆的长度均为1.8 m，弹性模量为220 GPa，泊松比为0.21。采用PFC分析动态力学问题时，需取合适的动态阻尼值。经模拟确定，颗粒法向阻尼比 ${\gamma _{\rm{n}}}$ 取0.32，颗粒切向阻尼比 ${\gamma _{\rm{s}}}$ 取0.44。基于FLAC网格节点构建的杆件模型优势在于：(1) 可以以任意“载荷-时间”函数关系，在杆件端部施加载荷；(2) 只用单一的、少量的建模单元即可实现实验过程的快速、合理模拟，极大地节约时间成本；(3) 良好的本构关系能够合理反映应力波在杆件中的均匀传递。另外，使用PFC-GBM方法构建的非均质岩石模型能够有效表征材料的损伤破坏特性。

图  7  耦合SHPB系统及应力波施加

Figure  7.  Coupling SHPB system and stress wave application

下载: 全尺寸图片 幻灯片

SHPB实验过程中，通过不同形状的子弹冲击入射杆产生入射波，关于子弹形状与入射波形之间的关系，人们已开展了大量研究^[26]。一般情况下，纺锤形子弹冲击杆件后产生正半弦波，圆柱形子弹冲击杆件后产生矩形波。在数值模拟中，常见的入射载荷施加方法有构造子弹冲击法^[27]和等效载荷施加法^[28]。其中，构造子弹冲击法可实现对实验过程的完整反演，但子弹的冲击速度v与冲击载荷p之间的关系需要标定，模拟过程较为复杂。等效载荷施加法是将设计的入射波形与加载时间的函数关系转换为入射载荷直接施加在入射端，可实现多种形式载荷的施加，模拟过程相对简便。

本研究选择等效载荷施加法实现正弦波入射载荷的施加。设入射载荷p与时间t满足函数关系$p = {p_{\text{m}}}\left[ {1 - \cos (2\text{π} f_0t)} \right]/2$，其中：p_m为入射载荷峰值，由于FLAC中拉应力为正，压应力为负，因此载荷p为负值；f₀为入射频率，当载荷施加时间为T时，有$T = 1/f_0$。

下面对等效载荷施加过程进行验证。T=200 μs、p_m=150 MPa时杆件的应力传播过程如图8所示。从图8中可以清晰地看到冲击应力在杆件中的入射、反射和透射过程。

图  8  耦合SHPB系统中杆件应力波的传播过程

Figure  8.  Stress wave propagation of the rod in coupling SHPB system

下载: 全尺寸图片 幻灯片

图9为杆件应力时程曲线，监测圆获取的杆件入射应力峰值同样为150 MPa，且杆件入射波形均与等效施加载荷一致。因此，该方法满足模拟要求。

图  9  杆件应力时程曲线

Figure  9.  Stress-time history curve of the rod

下载: 全尺寸图片 幻灯片

验证SHPB实验合理性的关键在于试件两端达到应力平衡过程是否符合一维应力波假设以及应力均匀性假设^[29]。图10给出了T=200 μs，p_m分别为150、200、250、300、350及400 MPa时的应力平衡关系曲线。可以看出，入射应力与反射应力之和与透射应力曲线的吻合程度较高，满足应力平衡条件。

下载: 全尺寸图片 幻灯片

图  10  不同冲击载荷下的应力平衡验证

Figure  10.  Uniformity of stress under different impact loading

下载: 全尺寸图片 幻灯片

2.2   动态冲击载荷下的微裂纹演化过程

为更清晰地展现冲击载荷下非均质晶质岩石的裂纹演化过程，选择破碎程度较高的p_m= 400 MPa的冲击过程进行分析。图11显示了在p_m=400 MPa的冲击载荷下不同节点时刻岩石内部微裂纹分布及宏观破坏模式。试样受到的整个冲击作用时间约为300 μs。应力波自318 μs左右经入射杆进入岩石，高强度的动态冲击载荷开始传入岩石内部；350 μs左右，岩石入射端首先产生微裂纹；382～443 μs，微裂纹数快速增长，在极短的时间内岩石内部的损伤程度显著升高，多处裂纹有贯通趋势；501 μs左右，岩石透射端出现损伤，岩石发生明显变形；532 μs左右，岩石内部多处微裂纹贯通，岩石变形程度及宏观破坏程度明显加大，但微裂纹数的增长速度明显放缓；614 μs左右，微裂纹数停止增长，试样不再发生进一步破坏，冲击破坏过程结束。

图  11  冲击载荷下岩石内部微裂纹及宏观破坏模式演化过程

Figure  11.  Evolution process of microcracks and macroscopic failure modes in rock under impact loading

下载: 全尺寸图片 幻灯片

根据不同时刻晶质岩石冲击破坏细观发展过程，可直观了解微裂纹的产生及分布情况，但是仍不足以定量地获得微裂纹的时空演化模式。为了更直观地了解冲击载荷下晶质岩石内部的微裂纹演化过程，经统计得到如图12所示的各矿物晶内及晶间微裂纹数时程曲线。各矿物的晶内裂纹及晶间裂纹分布如图13所示。

图  12  微裂纹数时程曲线

Figure  12.  Time-history curves of microcrack number

下载: 全尺寸图片 幻灯片

图  13  微裂纹细化分组

Figure  13.  Micro-crack refinement grouping

下载: 全尺寸图片 幻灯片

从晶质岩石动态冲击的微裂纹演化时程曲线可以看出，各矿物晶质内及晶质间的微裂纹增长模式均存在萌生期、快速增长期、缓慢增长期、停止增长期4个阶段。裂纹萌生期及快速增长期前期，矿物晶间裂纹数明显高于晶内裂纹数，说明在冲击载荷作用下岩石内部的裂纹演化首先以晶间黏结的断裂损伤为主，与静态下晶质岩石的单轴破坏过程相似。在裂纹快速增长期的中后期，长石矿物晶质内裂纹数超过晶间裂纹数，结合不同节点时刻下微裂纹演化过程及试样的宏观破坏模式，发现晶间裂纹的产生对晶内岩石的损伤具有一定的导向作用。对比各矿物的晶内损伤情况，微裂纹数由高到低分别为长石裂纹、石英裂纹、云母裂纹，与矿物质量分数排序一致。张涛等^[30]基于PFC-GBM对不同晶质尺寸的花岗岩进行了动态力学性能模拟，所得的各矿物晶内和晶间微裂纹数及其演化规律与本研究的模拟结果基本一致。究其原因，矿物晶内裂纹数与初始各矿物晶内黏结数密切相关。在岩石充分破坏的情况下，黏结键的断裂即在对应的黏结位置产生微裂纹，这是导致矿物晶内微裂纹数与对应矿物含量呈正相关的主要原因之一。晶间和晶内裂纹经缓慢增长期后，增长态势逐渐平稳，直至裂纹停止产生。

2.3   不同冲击载荷下岩石损伤破坏特征

表4列出了不同冲击载荷下非均质岩石破坏模拟结果，其中：N_d为微裂纹数，σ_d为峰值强度，${\dot \varepsilon _{\rm d}}$为峰值应变率。可以看出，对于T=200 μs的固定加载周期，随着冲击载荷的增强，岩石内的微裂纹数逐渐增加，破碎程度逐步增大，动态峰值强度和峰值应变率也同步增大。

表  4  不同冲击载荷下试样破坏的模拟结果

Table  4.  Simulation results of specimen failure under different impact loading

T/μs	pm/MPa	Crack distribution	Fragmentation distribution	Nd	σd/MPa	${\dot \varepsilon _{\rm d}}$/s−1
200	200			1017	120.56	60.28
	250			2386	143.94	76.65
	300			3241	162.85	104.59
	350			3950	184.84	114.55
	400			4776	203.12	157.64

下载: 导出CSV 

| 显示表格

对于常见的岩石材料，其动载强度高于甚至远高于静载强度。模拟得到的试样单轴抗压强度p_m为125.72 MPa，由不同冲击载荷下试样的应力-应变曲线可知，p_m=200 MPa的试样的动态峰值强度σ_d为120.56 MPa，与轴向应力较为接近，但冲击结果显示试样内部只产生少量微裂纹，无明显破坏模式，说明动、静态力学实验中岩石的力学性能存在差异；当p_m=250 MPa时，σ_d为143.94 MPa，微裂纹数N_d明显增加，试样入射端与透射端的损伤程度升高；当p_m=300 MPa时，σ_d为162.85 MPa，N_d急剧增加，试样内部的损伤程度逐渐升高，但无明显破坏模式；当p_m=350 MPa时，σ_d为184.84 MPa，试样上、下两端沿冲击方向的微裂纹逐渐有贯通趋势，出现“剥落”式破坏，试样中间部位虽产生大量微裂纹，但未出现明显的破坏；当p_m=400 MPa时，σ_d为203.12 MPa，试样发生高度破碎，从块度分布情况来看，试样最终发生“留心”式破坏。

从图14可以看出，与准静态加载下的单轴应力-应变曲线类似，对于固定的入射周期T=200 μs，不同冲击载荷下试样的应力-应变曲线在初始阶段均具有良好的线性关系。随着冲击峰值载荷p_m的不断增大，峰值应变率${\dot \varepsilon _{\rm d}}$逐渐增大，动态峰值强度σ_d和微裂纹数N_d同步增加，试样损伤破坏程度逐渐提高，均表现出良好的应变率相关性。此外，各曲线的弹性模量几乎一致，没有随应变率的升高发生明显变化。当p_m为200、250、300 MPa时，应力-应变曲线在下降阶段表现出明显的回弹现象，可能是由于试样的破坏程度较低，弹性形变较高；但当p_m增加至350和400 MPa后，回弹现象消失，峰后曲线快速下降，试样表现出良好的脆性破坏趋势。

图  14  不同冲击载荷下的应力-应变曲线

Figure  14.  Stress-strain curves under different impact loading

下载: 全尺寸图片 幻灯片

需要说明的是，冲击峰值载荷p_m的变化直接影响冲击时试样的峰值应变率${\dot \varepsilon _{\rm d}}$及动态峰值强度σ_d。图15给出了${\dot \varepsilon _{\rm d}}$和σ_d随p_m的变化关系，其中：${\dot \varepsilon _{\rm d}}$与p_m的拟合方程为${\dot \varepsilon _{\rm d}}$= 0.465 p_m–36.83，拟合精度R²为0.959，σ_d与p_m的拟合方程为σ_d= 0.412p_m+39.45，R²为0.998，二者均表现出良好的线性关系。这有助于快速确认冲击峰值载荷p_m与峰值应变率${\dot \varepsilon _{\rm d}}$及动态峰值强度σ_d的关系，对于降低模拟实验的时间成本等具有重要的现实意义。

图  15  σ_d和${\dot \varepsilon _{\rm d}}$关于p_m的线性拟合

Figure  15.  Linear fitting of σ_d and ${\dot \varepsilon _{\rm d}}$ with respect to p_m

下载: 全尺寸图片 幻灯片

3.   结　论

采用PFC-GBM方法构建了能够反映微观结构特征的非均质晶质岩石模型，并以此对晶质岩石在静态单轴载荷下的微裂纹演化过程进行了模拟分析；通过PFC-FLAC耦合SHPB系统，采用等效载荷施加方法，对晶质岩石在冲击载荷下的晶内及晶间裂纹演化过程和断裂损伤细观特征进行了研究，得到如下主要结论。

(1) 在单轴压缩条件下，晶间裂纹数高于晶内裂纹数，晶质岩石的断裂损伤主要以晶间裂纹为主导，进而引起晶内损伤破坏。

(2) 不同冲击载荷下，峰值应变率${\dot \varepsilon _{\rm d}}$与冲击峰值载荷p_m的拟合方程为${\dot \varepsilon _{\rm d}}$=0.465p_m–36.83，R²=0.959；动态峰值强度σ_d与p_m的拟合方程为σ_d=0.412p_m+39.45，R²=0.998，二者均表现出良好的线性关系。

(3) 从晶质岩石冲击破坏的微裂纹演化时程曲线及宏观破坏模式演化来看，各矿物晶质内及晶间的微裂纹增长模式均存在萌生期、快速增长期、缓慢增长期、停止增长期4个阶段。与静态单轴压缩下岩石破坏细观特征相似，冲击加载初期，晶间裂纹数明显高于各矿物晶内裂纹数，岩石主要发生沿晶损伤破坏；在岩石充分破坏的情况下，最终各矿物晶内裂纹数与对应的矿物含量呈正相关。