基于Voronoi 模型的脆性岩石微观结构数值模拟1)

邱 鑫 林 缅, 郑思平 陈天宇

∗(中国科学院力学研究所，北京100190)

†(中国科学院大学，北京100049)

∗∗(东北大学深部金属矿山安全开采教育部重点实验室，沈阳110819)

岩石颗粒尺度的微观结构控制着岩石的微观力学行为，从而控制着岩石的宏观力学反应。在颗粒尺度上，岩石的晶粒形状、晶粒矿物刚度、矿物接触分布和其他缺陷等导致了岩石微观结构的非均质性[1-3]，这将引起内部应力分布的不均一性，影响断裂萌生和断裂扩展的力学行为[4-5]，从而控制岩石在载荷作用下的破坏、变形能力和裂纹模式。然而，对于非均质性如何影响脆性岩石破坏机制的理解仍不完整[6]。

随着计算机计算能力的不断提高，数值模拟技术不断发展并提供了模拟岩石微观结构和研究其岩石力学特性的可能性。目前在模拟岩石的颗粒结构方面，主要有以下四种方法：(1) 颗粒流 PFC 实现的圆盘状颗粒[7-8]；(2) RFPA2D 代码实现的正方形颗粒[9]；(3) 离散元 UDEC[10]以及混合有限离散代码 Y-Geo[11]和 ELFEN[12]实现的三角形颗粒；(4) UDEC 实现的多边形颗粒[13-14]。从岩石实际的微观结构观察来看，多边形颗粒结构更加接近于真实的情况。在UDEC 中，泰森多边形划分方法被应用于生成岩石的多边形颗粒结构，即UDEC–Voronoi模型。在 UDEC–Voronoi 模型中，多边形随机形状的岩石颗粒与颗粒间接触嵌合自锁，接触面的力学特性由法向刚度、切向刚度、内聚力、抗拉强度和摩擦角等参数决定。

UDEC–Voronoi 模型被开发后，许多学者对其进行了研究和应用。Fabjan 等[15]对该模型进行了广泛的微观参数敏感性分析，提供了有效岩石模型校准的指南。同时，通过与连续介质模型的比较，提出UDEC–Voronoi 模型可以更好地表征岩石的膨胀、裂缝形态和峰后行为。而孙博等[16]则对比了Trigon模型和Voronoi 模型在微观参数敏感性方面的差异，展现了多边形颗粒结构模型和三角形颗粒结构模型在表征岩石力学特征时存在的差异。目前，UDEC–Voronoi 模型已成功模拟了砂岩、火成岩、油页岩等岩石的微观结构，并在岩石的裂纹起裂与扩展规律、压缩强度和变形破坏特征等方面的研究中得到了应用[17-20]。

本研究介绍了UDEC–Voronoi 模型的构建方法，构建了致密岩石的微结构模型。该模型可以模拟致密岩石的单轴压缩实验，还可以监测不同脆性致密岩石在压缩过程中的裂纹数量、裂纹分布以及裂纹密度在压缩过程中的变化特征。这些特征有助于我们更好地理解致密岩石的微观结构破坏机制，从而更好地理解致密岩石的宏观破坏机制。

1 UDEC–Voronoi 模型

本研究采用的是 UDEC–Voronoi 建模方法，该方法利用 UDEC 内建的 Voronoi 镶嵌模式的自动生成器，将一个特定区域的模型细分为随机大小的多边形。通过这种方式，多边形可以代表颗粒、颗粒组合或完整岩石中随机扰动样本的其他缺陷[3,10](图 1)。

1.1 运动方程

根据式(1) 中的新速度确定新的颗粒位置为

式中，xi为颗粒坐标位置，θ为颗粒绕质心转动的角度。

在UDEC 的计算过程中，每个时间步长都会产生新的离散块体位置，从而产生新的块体之间的接触力，利用作用在离散块体上的合力和合力矩计算各块体的直线加速度和角加速度，离散块体的速度和位移则通过对时间增量的积分来确定。这个计算过程不断重复，直至达到离散块体系统的平衡状态或出现一个持续的失效结果。

1.2 本构定律

在低应力条件下岩石内部的破裂主要为沿晶破裂，且Voronoi 模型中颗粒不会发生破坏，破坏只沿着颗粒边界发生，因此在本研究中颗粒本构模型选用各向同性弹性本构模型。各向同性弹性本构模型具有卸载时可逆变形的特征，应力应变定律是线性的，与路径无关。颗粒模型的变形特征由体积模量和剪切模量表征。接触的本构模型选取不含残余强度的库仑滑移−面接触模型，如图2 所示。接触模型的变形特征由其法向刚度kn和切向刚度ks表征，强度特征由其摩擦角φ、黏聚力c和抗拉强度T来表征。图 2 中 ∆τs为接触面切向应力改变量，为切向位移增量的弹性部分，∆us为总切向位移增量。

2 致密岩石的微观力学行为模拟

图2 离散块体接触本构[21]

本研究选取了 4 类典型的页岩岩样进行能量色散 X 射线光谱学分析，获取岩样的矿物组分，并采用 Jin 等[22]提出的基于矿物成分的脆性指数(式(3)) 计算了试样的脆性指数，整理成表1。接着按矿物组成随机分布4 组矿物颗粒，分别代表致密岩石的4 种矿物组分并赋予不同的弹性参数，不同矿物之间的接触被分别赋予不同的微观力学参数，受计算机能力的限制，构建的 Voronoi 模型的平均粒径为 1.5 mm (图 3)。

表1 不同脆性指数岩样地化参数

图 3 构建的 Voronoi 模型

式中，FQtz+Fsp+Cb为石英、长石和方解石的矿物含量总和，Ftotal为总体矿物含量，该值为100%。

通过单轴压缩模拟实验与岩样的单轴压缩实验比对(图4 中横轴正应变为轴向应变，负应变为环向应变；下同)，校准模型的细观力学参数，最终得到模型的细观力学参数如表2 和表3 所示。

图4 岩样和模型的单轴压缩应力−应变曲线

2.1 不同脆性岩石的压缩破坏特征

本研究接着构建了不同脆性岩石的单轴压缩模拟实验，4 个试样分别采取相同的加载条件，获取的应力−应变曲线如图5 所示。随着脆性的降低，岩石的单轴抗压强度从200 MPa 下降至150 MPa，弹性模量从34.5 GPa 减少至8.87 GPa，同时显现出更强的延性特征，峰前应力−应变曲线也从弹性类型转变为弹−塑性类型，应力−应变曲线的向下弯曲度不断提升，岩石的应变值也显著增加了。

表2 矿物弹性参数[23-24]

4 种模型的破坏形态分别如图6(a)∼图6(d) 所示。岩石出现的破坏主要为拉伸破坏，局部也可见一些剪切破坏。相较于A1 和 A2 高脆性岩石，A3 和A4 低脆性岩石破坏时部分区域出现滑移趋势，且试样的环向变形较大，侧帮和底部出现了较大程度的滑移和脱落。

表3 不同脆性指数模型矿物间接触力学参数

2.2 不同脆性岩石的裂纹数量变化

为了模拟不同脆性岩石在压缩过程中的裂纹扩展情况，本研究利用 FISH 编写了监测 Voronoi 模型内部裂纹变化信息的功能，并选取不同的参量进行分析。由于模型在压缩过程中主要发生拉伸破坏，故主要监测拉伸裂纹的数量变化。在轴向力不断施加的过程中，岩石的微裂纹数量变化情况如图 7∼图 10 所示。

压缩过程中微裂纹的数量呈现出分阶段的变化特点，根据微裂纹数量的变化特征，选取裂纹数量变化曲线上切线斜率水平达到2500 的特征点：裂纹起裂点的应力为σci，表征岩石内部开始出现微裂纹；选取切线斜率水平达到7500 的特征点：裂纹贯通点的应力为σcc，表征岩石内部微裂纹开始贯通，微裂纹急剧增长；σp为岩石的峰值强度。

图5 不同脆性岩石的单轴压缩应力−应变曲线

图6 岩石的破坏形态(红线为拉伸裂纹，蓝线为剪切裂纹)

图7 A1 试样的微拉伸裂纹变化

图8 A2 试样的微拉伸裂纹变化

图9 A3 试样的微拉伸裂纹变化

图10 A4 试样的微拉伸裂纹变化

不同脆性岩石的裂纹起裂点和裂纹贯通点显现出不同的特征。首先分析裂纹起裂点。以各试样的裂纹起裂点应力与峰值强度的百分比值作为衡量标准，脆性岩石的裂纹起裂点应力水平较低脆性岩石高约10%σp，裂纹数量则较低脆性岩石更少。其次分析裂纹贯通点。脆性岩石在峰前阶段存在明显的裂纹贯通点，而低脆性岩石微裂纹数量保持较为缓慢的增长趋势直到试样破坏，其裂纹贯通点与峰值应力点几乎重合，不存在明显的裂纹贯通点。此外，在峰前阶段，脆性岩石的微裂纹数量增加速度要明显低于低脆性岩石；到达应力峰值点时，低脆性岩石的微裂纹数量大约是高脆性岩石的1.5 倍，表明低脆性岩石需要大量的裂纹累积才能形成宏观贯通裂纹破坏。

2.3 不同脆性岩石的裂纹分布变化

为了直观地观察压缩破坏过程中微裂纹的扩展过程，绘制了岩石模型在不同加载阶段的微裂纹分布图(图11)，需要说明的是，图中的裂纹是该阶段的某一时刻的裂纹，所以存在加载初期的剪切裂纹随着加载的进行向拉伸裂纹转化的情况。可以发现，微拉伸裂纹倾向于与试样加载方向平行，而微剪切裂纹则与试样加载方向成不同角度。在裂纹起裂阶段，微拉伸裂纹和剪切裂纹的数量均较少，且裂纹的长度较短；到达裂纹贯通阶段后，裂纹的长度得到了显著的增长，裂纹开始相互贯通，短裂纹的数量明显减少，此时剪切裂纹为主导型裂纹；到达峰值应力阶段后，微剪切裂纹的比例开始降低，岩石开始出现破坏，此时拉伸裂纹成为主导型裂纹；到达峰后阶段，拉伸裂纹进一步增加，在拉伸、剪切裂纹共同的作用下，试样形成了更大尺度的破坏。

图11 A2 试样各阶段裂纹分布(红线为拉伸裂纹；蓝线为剪切裂纹)

2.4 不同脆性岩石的裂纹密度变化

为了更好地理解岩石破裂与裂纹扩展之间的联系，根据 Bristow 提出的有关裂纹密度的定义[25](式4)，将裂纹密度作为衡量岩石内部裂纹扩展程度的指标。为了更好地探究不同脆性岩石的裂纹密度变化与起裂应力、峰值应力之间的规律，增加了脆性指数分别为 0.86，0.74 和 0.5 的数值试件，其地化参数如表4 所示，增加的数值试件的构建过程与A1∼A4 试样完全一致。

式中，C为拉裂纹密度，Li为统计区域内第i条裂纹的长度，S为统计区域的面积。

表4 增加的数值试件的地化参数

图12 展示了不同脆性岩石的起裂应力和起裂点裂纹密度，从图中可以看出，随着岩石脆性指数的增加，岩石的起裂应力呈现增加趋势，起裂点的裂纹密度呈现减少趋势，高脆性的A4 试样的起裂点裂纹密度约为低脆性的A1 试样的2 倍，而A1 试样的起裂应力约为A4 试样的1.7 倍，可见高脆性岩石是在较高的应力状态下和较低的裂纹损伤状态下达到裂纹起裂点。换句话说，相较于低脆性岩石，高脆性岩石会产生更多高能量的裂缝[26]。

图12 不同脆性岩石的起裂点应力和裂纹密度

图13 展示了不同脆性岩石的峰值应力和裂纹密度，从图中可以看出，致密岩石的峰值应力仍然符合随着脆性指数的增加而增加的规律，A1 试样的峰值强度约为A4 试样的1.3 倍。但裂纹密度随着脆性指数的增加呈现先增加后减小的趋势。经过数据拟合，峰值点裂纹密度和脆性指数符合

式中，C为裂纹密度，B为脆性指数。

图13 不同脆性岩石的峰值点应力和裂纹密度

在实际的岩样测试中，试样的脆性指数往往是容易测试的，然后根据经验公式 (5) 可以对试样在峰值应力处的裂纹损伤水平进行估计。另外，从图13 中可以看出，当脆性指数约为0.6 时，峰值点裂纹密度和脆性指数的变化曲线出现拐点，可知在脆性指数小于 0.6 时岩石的性质出现了改变，该脆性指数可以作为区分页岩脆性的经验值。

3 结论和讨论

本研究介绍了UDEC–Voronoi 模型的构建方法，建立了致密岩石的微结构模型，进行了不同脆性岩石的微观力学行为模拟。结果显示，通过微结构参数的校准，UDEC–Voronoi 模型能够很好地模拟岩石力学特性曲线，表征岩石的变形和强度特性；在模拟致密岩石的微观结构，监测不同脆性岩石微观裂纹的生成与扩展过程，探究致密岩石的微观破裂机制方面，该模型有良好的应用潜力。通过数值模拟获得的研究结论主要有：

(1)在压缩过程中，致密岩石内部的拉伸裂纹逐渐增加，最终成为主导性裂纹，导致岩石的破坏形态多为拉伸破坏。高脆性岩石的单轴抗压强度要更高，具有更高的弹性模量；而低脆性岩石的塑性特征更加明显，峰前应力−应变曲线从弹性类型转换为弹−塑性类型。

(2)随着岩石脆性指数的增加，岩石的起裂应力呈现增加趋势，而起裂点的裂纹密度呈现减少趋势，高脆性岩石是在较高的应力状态下和较低的裂纹损伤状态下达到裂纹起裂点。

(3) 岩石的峰值应力随着脆性指数的增加而增加，峰值点裂纹密度和脆性指数符合一拟合公式，根据该公式可以在已知岩石的脆性指数的情况下对峰值应力处岩石的裂纹损伤程度进行估计。可以把脆性指数0.6 作为区分页岩脆性的经验值。

必须指出，本研究介绍的 UDEC–Voronoi 模型还存在一些限制，主要有以下几点：

(1) 构建的岩石微结构模型中含有石英、长石、伊利石和方解石四种矿物颗粒，在计算模拟中只考虑了致密岩石的一般性特征，并未考虑岩石学、矿物学和微观结构等细微特征，且未考虑四种矿物颗粒大小的不同。

(2)由于计算机能力的限制，构建的岩石模型的粒径并不等于实际的岩石颗粒尺寸，而岩石颗粒尺寸对岩石的破裂机制具有一定的影响。随着计算机能力的发展，模拟实际岩石的颗粒尺寸将有可能实现。

(3) UDEC–Voronoi 模型为二维模型，未来的研究将致力于把Voronoi 模型扩展到三维的岩石模型中，这将使得构建的微结构模型更加接近于岩石材料的实际状况。