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

doi:10.6052/1000-0879-20-433

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

邱鑫^*^,, 林缅^,^*^,^,²⁾, 郑思平^*^,, 陈天宇^**

^*中国科学院力学研究所,北京 100190

中国科学院大学,北京 100049

^**东北大学深部金属矿山安全开采教育部重点实验室,沈阳 110819

NUMERICAL SIMULATION OF BRITTLE ROCK MICROSTRUCTURE USING A VORONOI MODEL ¹⁾

QIU Xin^*^,, LIN Mian^,^*^,^,²⁾, ZHENG Siping^*^,, CHEN Tianyu^**

^*Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China

University of Chinese Academy of Sciences, Beijing 100049, China

^**Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China

通讯作者: 2)林缅,研究员,主要研究方向为非常规油气勘探开发。E-mail:linmian@imech.ac.cn

责任编辑: 王永会

收稿日期: 2020-10-12 修回日期: 2021-01-19 网络出版日期: 2021-04-20

基金资助:

1)国家自然科学基金重点项目.  42030808
国家自然科学基金重大项目.  41690132
中国科学院战略性先导科技专项.  XDA14010304
国家自然科学基金面上项目.  41872163

Received: 2020-10-12 Revised: 2021-01-19 Online: 2021-04-20

作者简介 About authors

摘要

基于离散元软件UDEC中的Voronoi模型,构建了致密岩石的标准压缩数值试件;应用FISH语言编译相关程序,监测不同脆性致密岩石的裂纹数量和密度在压缩过程中的变化特征。研究发现,随着岩石脆性指数的增加,岩石的起裂应力逐渐增加,而起裂点裂纹密度逐渐减小,高脆性岩石则是在较高的应力状态和较低的裂纹损伤状态下达到裂纹起裂点。岩石的峰值应力随着脆性指数的增加而增加,峰值点裂纹密度和脆性指数符合一拟合公式,根据该公式可以对峰值应力处岩石的裂纹损伤程度进行估计。可以把脆性指数0.6作为区分页岩脆性的经验值。

关键词： 离散元方法 ; Voronoi模型 ; 脆性岩石 ; 岩石微结构 ; 裂纹密度

Abstract

Numerical compression specimens with different brittleness are constructed based on the UDEC-Voronoi model. The number of cracks and the crack density are monitored by the FISH programs. It is shown that the crack initiation stress increases with the increase of the brittleness index. However, the crack density at the crack initiation stress decreases. That is to say, the crack initiation of high brittle rocks is at a higher stress and with a fewer cracks as compared with the low brittle rocks. Furthermore, the peak stress increases as the brittleness index increases. A fitting formula is established to estimate the crack density at the peak stress against the brittleness index, and the brittleness index 0.6 can be used as an empirical value for judging the shale brittleness.

Keywords： distinct element method ; Voronoi model ; brittle rocks ; rock microstructure ; crack density

PDF (14428KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

邱鑫, 林缅, 郑思平, 陈天宇. 基于Voronoi模型的脆性岩石微观结构数值模拟 ¹⁾. 力学与实践[J], 2021, 43(2): 244-251 DOI:10.6052/1000-0879-20-433

QIU Xin, LIN Mian, ZHENG Siping, CHEN Tianyu. NUMERICAL SIMULATION OF BRITTLE ROCK MICROSTRUCTURE USING A VORONOI MODEL ¹⁾. MECHANICS IN ENGINEERING[J], 2021, 43(2): 244-251 DOI:10.6052/1000-0879-20-433

岩石颗粒尺度的微观结构控制着岩石的微观力学行为,从而控制着岩石的宏观力学反应。在颗粒尺度上,岩石的晶粒形状、晶粒矿物刚度、矿物接触分布和其他缺陷等导致了岩石微观结构的非均质性^[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

新窗口打开| 下载原图ZIP| 生成PPT

图1 Voronoi随机节理生成算法

1.1 运动方程

在UDEC中,颗粒在二维空间中受多种作用力$F_i^{(t)}$和力矩$M^{(t)}$,其速度方程为

(1) $\begin{eqnarray} \left.\begin{array}{l} \dot{u}_{i}^{(t+\Delta t / 2)}=\dot{u}_{i}^{(t-\Delta t / 2)}+\left(\dfrac{\Sigma F_{i}^{(t)}}{m}+g_{i}\right) \Delta t\\ \dot{\theta}^{(t+\Delta t / 2)}=\dot{\theta}^{(t-\Delta t / 2)}+\dfrac{\Sigma M^{(t)}}{I} \Delta t \end{array}\right\} \end{eqnarray}$

式中,$\dot{u}_{i}$为速度,$g_i$为重力加速度,$\dot{\theta}$为颗粒绕质心转动的角速度,$I$为转动惯量。

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

(2) $\begin{eqnarray} \left. \begin{array}{l} x_{i}^{(t+\Delta t)}=x_{i}^{(t)}+\dot{u}_{i}^{(t+\Delta t / 2)} \Delta t \\ \theta^{(t+\Delta t)}=\theta^{(t)}+\dot{\theta}^{(t+\Delta t / 2)} \Delta t \end{array}\right\} \end{eqnarray}$

式中,$x_{i}$为颗粒坐标位置,$\theta$为颗粒绕质心转动的角度。

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

1.2 本构定律

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

图2

新窗口打开| 下载原图ZIP| 生成PPT

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

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

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

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

新窗口打开| 下载CSV

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 构建的Voronoi模型

(3) $\begin{eqnarray} B_{\text {compo }}=\frac{F_{\rm Q t z+F s p+C b}}{F_{\text {total }}} \end{eqnarray}$

式中,$F_{\rm Q t z+F s p+C b}$为石英、长石和方解石的矿物含量总和,$F_{\text {total }}$为总体矿物含量,该值为100%。

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

图4

新窗口打开| 下载原图ZIP| 生成PPT

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

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

新窗口打开| 下载CSV

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

新窗口打开| 下载CSV

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

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

图5

新窗口打开| 下载原图ZIP| 生成PPT

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

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

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

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

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

图7

新窗口打开| 下载原图ZIP| 生成PPT

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

图8

新窗口打开| 下载原图ZIP| 生成PPT

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

图9

新窗口打开| 下载原图ZIP| 生成PPT

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

图10

新窗口打开| 下载原图ZIP| 生成PPT

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

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

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

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

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

图11

新窗口打开| 下载原图ZIP| 生成PPT

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

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

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

(4) $\begin{eqnarray} C=\Sigma L_{i}^{2} / S \end{eqnarray}$

式中,$C$为拉裂纹密度,$L_{i}$为统计区域内第$i$条裂纹的长度,$S$为统计区域的面积。

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

新窗口打开| 下载CSV

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

图12

新窗口打开| 下载原图ZIP| 生成PPT

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

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

(5) $\begin{eqnarray} C=-0.418 3 B^{2}+0.527 1 B-0.035 4 \end{eqnarray}$

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

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

图13

新窗口打开| 下载原图ZIP| 生成PPT

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

3 结论和讨论

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

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

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

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

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

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

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

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

责任编辑: 王永会

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Lan

, Martin

, Hu

Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading

Journal of Geophysical Reseach, 2010,115:B01202