基于颗粒结构的岩石破裂演化细观模型、方法及电子设备

技术领域

本发明涉及离散元模型技术领域，具体涉及一种基于颗粒结构的岩石破裂演化细观模型、方法及电子设备。

背景技术

岩石是由多种矿物颗粒胶结而成的复杂地质材料，与金属类材料相比，岩石通常不具有连续性和各向同性等特征。另外，考虑到岩石中天然存在的孔隙和裂隙，以及构成岩石矿物颗粒类别的差异，岩石通常也不能被简单的视为均质性材料。因此，岩石在荷载作用下的力学响应十分复杂，具有明显的非线性特征，导致采用传统的弹性或弹塑性理论均不能很好的拟合其应力-应变关系。另一方面，岩石在宏观层面的力学特性究其根本是细微观因素的宏观反映。根据以往的研究，岩石矿物颗粒的构成、颗粒粒径大小、孔隙率、裂隙密度等细观因素是影响宏观力学特性的几个主要因素。

目前，考虑到理论模型构建以及获得解析解方面的困难，学术界通常采用数值模拟的方法进行岩石类材料的应力-应变关系研究，且主要分为有限单元法和离散单元法。由于离散单元法可以反映岩石在受载变形过程中的大变形、断裂等主要特征，成为了主流的模拟方法。同时考虑到岩石的细微观特征，建立能够同时反映上述细观特征的离散元模型是进行精细化模拟的前提。

在离散元理论中，岩石的宏观力学性质由离散单元和单元间接触面的力学参数直接控制。在计算迭代过程中，块体的位移、速度及加速度通过牛顿第二定律计算，接触面的应力-应变关系则通过库伦摩擦定律计算。

单个块体的运动由不平衡力矩的大小和方向以及作用在其上的力决定，牛顿第二定律可以表示为如下形式：

其中

由于力取决于位移，因此力-位移关系可以由上式计算得到。而对于由多个方向的力以及重力作用的二维情况，则可由式(3)计算：

其中

相邻块体之间的受力关系由接触面力学参数决定，并由库伦摩擦定律控制。在垂直于颗粒边界的方向，接触面上单位面积的力-位移关系是线性的，由法向刚度kn控制，如式(4)所示：

Δσ

式中Δσ

在剪切方向上，由剪切刚度k

如果颗粒间剪应力小于上限值：

则有：

Δτ

当剪应力大于上限值时：

采用残余粘聚力和残余摩擦角进行抗剪强度计算：

式中Δu

上述各类强度参数，如法向刚度、剪切刚度、粘聚力(初始、残余)、内摩擦角(初始、残余)及抗拉强度(初始、残余)可以通过对接触赋值的方式分配给各个离散块体。这些参数均称为细观或微观特性，且这些参数的组合控制着所模拟材料的宏观力学行为。

根据目前的研究成果，较为先进离散元模型应能够反映岩石的如下特征：1.考虑到岩石在损伤过程中矿物颗粒的破碎是一种典型特征，因此单个矿物颗粒应能够发生破碎和断裂；2.考虑到花岗岩是由多种矿物颗粒组成的集合体，由于不同类型的矿物颗粒具有不同的力学性能，因此模型应能够反映成分的异质性；3.考虑到颗粒层面的结构性，模型中的颗粒粒径分布应与真实岩石具有较好的一致性。

发明内容

针对现有技术的不足，本发明公开了一种基于颗粒结构的岩石破裂演化细观模型、方法及电子设备，为了探究岩石细观层面特征对宏观力学性质的影响规律。

本发明通过以下技术方案予以实现：

本发明公开了一种基于颗粒结构的岩石破裂演化细观模型，包括用于在指定区域内生成模拟晶体颗粒的密排多边形或多面体的T模块；用于进行单个颗粒的网格划分的M模块；用于进行模型的图像输出的V模块。还包括用于进行模型数据转换的Neper-UDEC接口程序；用于多边形网格精细化处理并构建反映岩石细观结构特征的python程序。所述模型为通过自编Python程序，结合Neper插件所生成的初步多边形网格建立可用于UDEC离散元计算且可体现单个矿物颗粒可破碎特性、矿物颗粒粒径分布特征以及岩石多组分特征的一种模型。

一种基于颗粒结构的岩石破裂演化细观方法，所述方法所述基于颗粒结构的岩石破裂演化细观模型实现，所述模型基于岩石细观组构特征建立，根据细观观测试验数据，采用Neper生成初步的多边形网格，该过程包含以下步骤：

S1利用模块T进入网格生成模式，通过-dim命令指定区域的维度，并通过-domain命令指定颗粒生成区域的具体范围；

S2根据岩石矿物颗粒的观测结果，通过-morpho命令进行颗粒几何特征的控制，通过diameq及sphericity参数并通过对数正态分布函数(lognormal)进行控制。

S3输出分别记录了每个颗粒的编号、中心坐标、边数及边的编号、顶点及顶点编号信息的.tess类型文件。

模型的建立以初步网格文件为基础，采用Python语言编写接口程序，建立Neper所生成.tess文件与UDEC程序的联系。

进一步的，单个所述多边形网格对应单个矿物颗粒，通过建立单个颗粒内部接触，实现颗粒可破碎特性。

进一步的，所述网格文件结合Python程序进行岩石不同矿物成分的赋值，具体为：

T1建立不同种类矿物块体属性、晶内接触属性；

T2在进行单个颗粒参数赋予前，建立(0，1)区间内的随机数，根据矿物含量将该区间分为多个子区间，随机数落在某区间内，则对该颗粒赋予相应的参数；

T3单个颗粒赋值完毕后，更新第T2步中的随机数，并对下一个颗粒重复上述操作，直到所有颗粒都被赋值后，结束循环。

一种电子设备，包括处理器以及存储有执行指令的存储器，当所述处理器执行所述存储器存储的所述执行指令时，所述处理器执行上述方法。

本发明的有益效果为：

本发明的模型可以用于模拟岩石类材料的受载破坏特征，可以联系材料的细观结构和宏观力学特征，并直观的展示了岩石内部的破裂损伤过程。根据模拟试验的结果，模型可以准确的再现岩石低拉压比特征，同时模型的破裂模式与实测结果吻合；验证了在0-50MPa围压作用下的单轴及三轴强度并采用摩尔-库伦准则进行了拟合，模拟结果和实测结果具有良好的一致性。能够准确的反映岩石的宏观力学特性，直观的再现岩石内部复杂的裂纹发展过程，适用于深部岩石的细观研究。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是颗粒可破碎特性设置方法示意图；

图2是颗粒参数赋值方法示意图；

图3是Neper-UDEC初步模型构建流程图；

图4是初步多边形模型粒径分布与实测数据对比图；

图5是基于颗粒结构的岩石破裂演化细观模型示意图；

图6是不同接触面抗拉强度模型的应力-应变曲线图；

图7是模型接触面抗拉强度与宏观抗拉强度关系图；

图8是模型k

图9是模型接触刚度值变化对杨氏模量和泊松比的影响示意图；

图10是接触面粘聚力与抗压强度关系图；

图11是接触面摩擦角与抗压强度关系图；

图12是单轴压缩模拟试验应力-应变曲线图；

图13是三轴压缩试验实测与模拟数据对比图；

图14是不同围压下裂纹形态差异示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

本实施例公开如图3所示的一种基于颗粒结构的岩石破裂演化细观初步模型的建立流程，所述方法包括以下步骤：

S1基于岩石细观试验数据，进入Neper程序，通过-dim命令指定区域的维度，并通过-domain命令指定颗粒生成区域的具体范围；

S2根据岩石矿物颗粒的观测结果，通过-morpho命令进行颗粒几何特征的控制，通过diameq及sphericity参数并通过对数正态分布函数进行控制；

S3输出分别记录了每个颗粒的编号、中心坐标、边数及边的编号、顶点及顶点编号信息的.tess网格文件。

由于Neper程序对本文使用的UDEC6.0没有相应的接口文件，因此无法直接生成UDEC可导入的数据文件。但上述.tess文件中记录了建立UDEC模型所需的各类数据，在了解.tess文件数据结构的情况下，可以采用编程的方式生成UDEC可读取的数据文件，将Neper生成的数据转换到UDEC中。

采用上述方法进行初步的模型建立，模型中颗粒粒径与统计值的对比如图4所示，可见两者具有较好的一致性。由于本方法在二维条件下研究，因此可以调节模型颗粒分布的自由度较低，三维情况下可以更精确的模拟真实的粒径分布。

实施例2

本实施例针根据文献的研究，穿晶断裂是岩石材料受载破坏过程中的一个重要特征，单个矿物颗粒的破坏特性能否再现，以及模型能否合理的反映沿晶破坏和穿晶破坏的差异对模拟结果有着显著的影响。公开一种颗粒离散元模型中实现单个颗粒可破碎特性的设置方法，如图1所示，根据Neper程序生成的多边形块体，由于已知多边形各顶点的坐标位置，则可以根据式(10)计算得到块体的几何中心：

其中p

为了区分岩石中的不同矿物成分，需要对由同一个多边形分割而成的多个三角形赋予相同的块体属性，并且需要区分颗粒间接触以及颗粒内部接触，对不同类型接触赋予不同的强度参数。UDEC中的changereg命令可以指定一个由四个坐标点组成的凸四边形，当块体中心点位于四边形内时，该块体的材料参数将会改变。同理，当接触的中心位于指定的四边形内部时，其材料参数也会随之改变。如图2所示，通过Python程序循环提取.tess文件中某颗粒的顶点坐标，并计算该颗粒的几何中心。以几何中心以及相邻的三个顶点坐标建立四边形区域，如图中OABC。根据A、B、C三点以及中心点的坐标，可以计算得到OA、OB、OC三个向量。以OA向量为例，设向量的长度为L，计算得到与向量OA方向相同，同样以O为原点的向量Oa，Oa的长度L

对于岩石的矿物成分组成，可以根据观测结果，结合Python程序容易的实现：

1.首先建立不同种类矿物块体属性、晶内接触属性；

2.在进行单个颗粒参数赋予前，建立(0，1)区间内的随机数，根据矿物含量将该区间分为多个子区间，随机数落在某区间内，则对该颗粒赋予相应的参数；

3.单个颗粒赋值完毕后，更新第2步中的随机数，并对下一个颗粒重复上述操作，直到所有颗粒都被赋值后，结束循环。

本实施例模型构建方面，首先结合Neper插件生成初步的多边形网格，通过建立单个颗粒内部接触的方法实现了颗粒的可破碎特性，并编写Python程序对颗粒块体及接触进行了精准赋值。所提出的颗粒参数赋值方法可以准确的对晶内接触及晶间接触进行区分，再现岩石破裂过程中的穿晶和沿晶破坏，并可以根据矿岩成分分析结果建立具有成分异质性的复杂模型，体现多类矿物颗粒强度差异对整体力学性质的影响，最终效果如图5所示。

实施例3

在离散元模型中，模型的宏观抗拉强度受到接触面抗拉强度的影响。本实施例采用巴西劈裂实验模型进行模拟，研究接触面的抗拉强度与模型整体抗拉强度的相关性。模拟结果如图6所示，当接触面的抗拉强度增加时，模型的宏观抗拉强度相应增长，且两者呈现线性关系(图7，R

在本实施例模型中，接触面的法向刚度k

采用所提出的模型，通过调整两个刚度值探究其对材料泊松比和杨氏模量的影响。如图8所示，各个模型采用相同的法向刚度，晶间、长石、石英及黑云母法向接触刚度分别设为为2.3×10

实施例4

本实施例中对模型的接触面抗剪强度参数分析，接触面抗剪强度由库伦摩擦公式计算，其中接触面的粘聚力c和摩擦角

从图10及图11可以看出，接触面摩擦角及粘聚力均对三轴强度产生了显著的影响。与上述分析相同，当保持摩擦角恒定而增加接触面平均粘聚力时，材料的宏观强度所有增加。同时，随着围压的增加，三组模型的抗压强度增加幅度近乎一致，三组强度包络线呈现近乎平行的关系。另一方面，当保持接触面粘聚力相同但采用不同的摩擦角时，三组模型的包络线呈现差异性的斜率。分别将摩擦角的正切值增加和减小20％后，材料在50MPa围压下分别达到了365MPa、460MPa以及521MPa。随着内摩擦角的增加，颗粒间相对滑动时所需的剪切应力随之增加，因此更难以发生剪切破坏。而在高围压下岩石的断裂主要由剪切破坏控制，因此提高平均内摩擦角使得材料的三轴抗压强度快速增加。

通过调参过程，可以总结出细观参数对岩石宏观力学特性及破坏特征的影响规律，如表1所示。

表1细观参数对宏观力学性能的影响规律

1.当单独降低接触的抗拉强度时，材料的整体强度降低。由于接触的法向刚度k

2.当同比例降低接触的法向刚度k

3.当同比例降低接触面法向刚度k

4.同时增加接触面粘聚力c以及接触面内摩擦角

5.降低接触面粘聚力c的同时增加接触面内摩擦角

实施例6

本实施例阐述基于试验数据的验证算例，岩石模型采用1∶1比例建模，模型高100mm，直径50mm。模型上下设置加载板，下方加载板固定所有自由度，并在加载板内部设置轴向应力监测点。在模型两侧设置了5组10个横向位移监测点，用于监测模型的横向变形。为了消除端部效应引起加载板附近的复杂应力分布，岩石块体与加载板间的切向刚度、与抗剪强度相关的参数、与抗拉强度相关的参数均设置为极小值，用以模拟岩石与加载板间的光滑接触特性。通过对上方加载板施加恒定速度进行轴向加载，并在上方加载板底部设置位移监测点。本次实验中对试样轴向施加0.15m/s的恒定速度。

通过对模型在加载过程中应力及应变的监测，可以绘制出如图12所示的应力-应变曲线。图中黑色为轴向应变-应力曲线，红色为环向应变-应力曲线，蓝色为体应变-应力曲线。可以看出，模型可以较好的再现岩石单轴压缩实验中的各向特征：

1.加载初期呈现典型的线弹性特征。另外，离散元模型的一个主要优势是可以再现加载过程中的裂纹积累，以往研究表明岩石的启裂应力可以根据裂纹数量达到峰值处裂纹总数的1％进行判断。本文中计算得到的启裂应力为65.2MPa，约为峰值应力的39％，与实际较为相符。

2.随着加载的继续进行，岩石逐渐被压缩，体积应变持续增加。在达到峰值应力的85％处，曲线发生了明显的非线性增长，此时对应岩石的加速破裂阶段，随应变增加轴向应力增长缓慢。实测得到的损伤应力在峰值应力的80％左右，可见该模型能够较准确的反应岩石的破裂损伤情况。

3.达到峰值应力后，曲线发生迅速的跌落，呈现出典型的脆性破坏特征。此时，横向应变急剧增加，从0.2％快速增长至1％左右。此时岩石内部已经含有贯通的裂隙，岩石体积膨胀主要由裂隙开度的持续增加提供。因此对于体应力-应变曲线，其在峰值应力前为正值，代表岩石整体处于压缩状态。接近峰值应力时，体应力-应变曲线开始向负值方向发展，代表了岩石从压缩状态向扩容状态的转变。超过峰值应力后随环向应变的增加快速向负值方向发展，岩石发生体积的迅速增加。

实施例7

本实施例通过模拟巴西劈裂试验，验证模型的抗拉性能。模拟结果表明，当荷载达到峰值抗拉强度的85％时，模型开始发生接触面的断裂，最初出现的裂隙主要是晶粒间的拉伸破坏，分布在轴线附近且孤立存在。接近峰值强度时，轴线附近的张拉裂隙持续增加且相互连通。圆盘中心处首先出现与加载方向近似平行的贯通张拉裂纹。通过监测加载过程中拉伸裂纹和剪切裂纹，张拉裂隙在临近峰值阶段迅速积累，而剪切裂纹数量远低于张拉裂纹，且在峰值处才开始出现。同时轴线处形成了贯通的宏观张拉裂隙带，与实验结果高度吻合。计算得到的抗拉强度为8.35MPa，为单轴抗压强度的4.99％，可见模型很好的反映了花岗岩低拉压比的特征。

实施例8

本实施例进行三轴压缩试验模拟，采用一组三轴压缩试验验证模型模拟侧限条件下的性能。并采用摩尔-库伦进行强度包络线的拟合。从图13中可以看出，模型的包络线与测试曲线高度吻合。模拟得到的三轴强度与实测值相符，实测值与模拟值计算得到的粘聚力分别为39.28MPa以及36.44MPa，内摩擦角分别为44.86°及44.77°。如图14，在相对较小的20MPa围压作用下，主裂隙与加载方向夹角较小，同时局部出现了倾斜的剪切裂隙。整体上，破坏呈现出脆性特征，岩石裂纹数量较少且相对平直，没有出现明显的塑性破坏区，除两条主裂隙外，岩石其他部位颗粒基本无破碎。而在相对较大的50MPa围压作用下，显示出了明显的塑性特征。主裂隙与加载方向的角度增大，裂隙不再保持平直的形态，并且出现了明显的分叉。裂隙面两侧存在明显的晶体碎屑，且主裂面处发生了块体的剥落。从模型计算结果来看，两组围压下的模型较好的反映了上述特征。20MPa围压作用下模型的断裂由三条明显的剪裂隙和若干张拉裂隙组成。其中主剪裂隙从模型底部延伸至模型右上部，由张拉微裂纹和剪切微裂隙共同组成。50MPa围压作用的模型主要受一组共轭剪裂隙控制，剪裂隙相比20MPa模型更加倾斜，且剪切滑移裂隙的数量明显增多。在模型左下部出现了大量裂隙的集中，呈现出局部塑性破裂的特征。

本发明实施例结合三个验证算例说明了该模型在模拟花岗岩破裂演化过程的适用性。单轴压缩实验模拟中，该模型准确的反映了岩石的启裂应力及损伤应力，并直观的展示了岩石内部的破裂损伤过程；通过巴西劈裂实验模拟，再现了花岗岩低拉压比的特征，同时模型的破裂模式与实测结果吻合；验证了在10-50MPa围压作用下的三轴强度并采用摩尔-库伦准则进行了拟合，模拟结果和实测结果具有良好的一致性。综上，该模型能够准确的反映岩石的宏观力学特性，直观的再现岩石内部复杂的裂纹发展过程，适用于深部岩石的细观研究。

以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。