CN110069844B - 一种考虑岩石结构特征和矿物组成的细观模型生成方法 - Google Patents

Abstract

本发明公开了一种考虑岩石结构特征和矿物组成的细观数值模型生成方法。能够考虑岩石的细观结构特征和矿物组成，建立反映岩石结构特征和矿物组成的细观数值模型，能够更加真实进行岩石力学数值模拟。本发明包括以下步骤：对岩石试样进行岩石结构特征分析得到矿物颗粒的尺寸、形态、矿物类型等信息；将不规则形状的矿物颗粒简化为一尺寸缩小的圆球，采用粒径膨胀法使颗粒发生重排列，获得一个较为密实的结构；然后提取颗粒空间位置和尺寸等信息，对颗粒集合体进行基于Voronoi图的空间区域划分；最后将得到的Voronoi多晶体结构进行有限元网格划分，并在矿物颗粒边界和颗粒内部***无厚度界面单元，生成可反映岩石结构和岩相结构的有限元数值模型。

Description

一种考虑岩石结构特征和矿物组成的细观模型生成方法

技术领域

本发明属于岩石力学研究领域，具体涉及一种考虑岩石结构特征和矿物组成的细观数值模型生成方法。

背景技术

岩石是地壳的基本组成物质，大量出露于地表，构成山川峡谷，是人类工程活动的基本载体和环境。岩石是由一种或几种矿物组成的，具有稳定外形的固态集合体，而矿物又具一定的化学成分和结晶构造。岩石结构指的是组成岩石的物质的结晶程度、矿物颗粒的大小、矿物颗粒的相对大小、矿物颗粒的形状以及它们之间的相互关系所表现出来的特征。按岩石的结晶程度可将结构分为全晶质结构、半晶质结构和玻璃质结构三大类；按岩石中矿物颗粒的绝对大小可分出粗粒、中粒、细粒、微粒等级别的结构；按矿物颗粒的相对大小又可划分等粒结构、不等粒结构和斑状、似斑状结构；按岩石中矿物的自形程度，还可分出自形结构、半自形结构及他形结构。

在岩石的力学特性研究过程中，矿物颗粒的尺寸、形状、互锁程度和接触类型都具有不同的影响。从矿物颗粒大小看，在岩浆岩、变质岩和沉积岩中，等粒结构一般比非等粒结构强度高；在等粒结构中，细粒结构比粗粒结构强度高。在斑状结构中，细粒基质比玻璃基质强度高；粗粒具斑晶的酸性深成岩强度最低；细粒微晶而无玻璃质的基性喷出岩强度最高。从结构连结上看，大部分岩浆岩、变质岩及沉积岩中的化学岩，晶粒之间结合紧密，强度较高，但在化学岩中，以可溶性结晶连结的，强度较高，但抗水性差。固结粘土岩有一部分是再结晶连结，其强度比其它坚硬岩石差得多。

国内外的许多学者对岩石的宏观变形和强度特性与岩石结构的关系开展了大量的研究。

2013年，Coggan等人研究了英格兰西南部花岗岩的力学性能，认为高岭石化的增加是导致该地区花岗岩强度显著降低的主要原因；同时Sousa等人也对葡萄牙不同地区的花岗岩进行了对比分析，认为岩石中石英的纹理特征对岩石的力学行为有至关重要的影响。

2014年，Sajid和Arif等人调查了巴基斯坦西北部的花岗岩，他们发现岩石内部的岩相变化将会导致细粒花岗岩强度降低，同时还会导致岩石的吸水率增加。

近年来一些先进的观测技术，如声发射、扫描电子显微镜、X射线断层扫描技术、数字图像相关技术等被引入到岩石力学试验，研究岩石等准脆性材料在荷载作用下，微裂纹的成核、扩展和汇合的物理过程。由于岩石内部的矿物成分、晶体构成十分复杂，其力学性能尚不能通过室内试验很好的揭示。为了深入理解岩石微观结构对其力学性能的影响，已经有学者开始进行考虑岩石微观结构的数值模拟，建立相应的微观模型来反映岩石的结构特性。

2014年，Ghazvinian等人提出了一种三维随机Voronoi晶粒模型用于模拟脆性岩石的裂纹扩展，通过Voronoi镶嵌获得晶体模型，并用晶体边界模拟晶粒间的裂纹。

2017年，Peng等人提出了一种晶粒模型模拟了武吉知马花岗岩的微观开裂行为，通过构造多边形晶粒模拟结晶岩石的微观结构。

但是，上述方法并不能完全反映岩石微观结构，仅能考虑岩石的矿物成分和矿物含量，不能考虑岩石矿物颗粒的不规则形状和排列结构，因此需要建立考虑岩石结构特征和矿物组成的细观数值模型。

发明内容

本发明是为了解决上述不足而设计的，目的在于提供一种能反映岩石结构特征和矿物组成的细观数值模型生成方法。基于岩石结构特征分析结果，如岩石矿物颗粒含量、矿物粒度分布、晶粒尺寸、颗粒形状和晶粒的空间位置，更加真实的反映岩石微观结构对其宏观力学性质的影响。

本发明为了实现上述目的，采用了以下方案。

一种考虑岩石结构特征和矿物组成的细观数值模型生成方法，其特征在于,包括以下步骤：

步骤1，通过岩石结构特征分析，得到矿物颗粒分布和颗粒几何信息：

步骤2，离散元模拟；

将由步骤1中将岩石结构特征分析得到的矿物颗粒分布和颗粒几何信息收集汇总，将不规则形状的矿物颗粒简化为一尺寸缩小的圆球，采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程，直至颗粒体系和真实岩石试样具有相同的结构特性即颗粒粒径是否达到矿物颗粒尺寸分布时，停止颗粒粒径膨胀及离散元数值模拟，完成岩石试样建模；

步骤3，生成数值模拟；

提取岩石试样模型中颗粒空间位置和尺寸信息，对颗粒集合体进行基于Voronoi图的空间区域划分，每个Voronoi元胞对应一个不规则形状的矿物晶粒；根据岩石中不同矿物类型的体积分数，给每个Voronoi元胞代表的矿物晶粒分配矿物类型；

步骤4，生成有限元模型；

将得到的Voronoi多晶体结构进行有限元网格划分，并在矿物颗粒边界和颗粒内部***无厚度界面单元，生成可反映岩石结构和岩相结构的有限元数值模型，完成细观数值模型生成。

作为优选，步骤1中，岩石结构特征分析具体方法为：

首先将岩石切片，通过显微观察得到岩石试样切片图，然后对岩石试样切片图进行数字图像处理，对二值化的图像采用分水岭算法分割出矿物颗粒，识别矿物颗粒的分布和颗粒几何信息。

作为优选，步骤1中，所述矿物颗粒分布信息包括矿物颗粒的矿物种类、不同矿物种类的含量、粒径大小、粒径分布、颗粒的形心位置、圆度以及面积，颗粒几何信息包括矿物颗粒的尺寸、实心度、长短轴之比以及圆度。

作为优选，所述圆度通过下面公式定义：

公式一中，S_grain是通过岩石结构特征分析得到的矿物颗粒面积，P_grain是矿物颗粒周长，π为圆周率。

作为优选，步骤2中，采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程中，颗粒在给定的边界内开始按照相同的速率体积膨胀，在膨胀过程中，受试样边界限制和颗粒间互不侵入条件的约束，颗粒在生长过程中发生重排列，趋于一个密实和各向同性结构，上述采用离散单元法模拟岩石试样的颗粒集合体的结构特性可以采用径向分布函数和键取向序参数度量，通过这两个参数判断颗粒粒径是否达到矿物颗粒尺寸分布，若不满足要求则继续进行粒径膨胀；当满足要求时停止离散元数值模拟。

作为优选，所述径向分布函数为：

公式二中，dn(r)表明在距离原点颗粒相距r处，且在宽度dr范围内所能找到的颗粒数， N为岩石试样中总矿物颗粒数，V为岩石试样建模的体积，g(r)的物理意义是与原点颗粒相距 r处单位体积的颗粒数密度，从而获得颗粒体系的结构信息，π为圆周率。

作为优选，所述键取向序参数函数为：

公式三中，

是球谐函数，其中l表示对称性的序参量，m取-l≤m≤l，

即为颗粒接触对的方向向量，而

和

为该方向向量的极角和方位角，π为圆周率，Q_l为键取向序参数值。

作为优选，步骤3中，对颗粒集合体进行基于Voronoi图的空间区域划分指的是一种在三维空间内分区的技术，对于D∈R³的区域内，R指的是长度量纲，R³表示在三维空间内给定的D中有若干种子点，D中种子点的坐标由粒径膨胀法获得，对于上述种子点i，满足{S_i(x_i)}(for i＝{1,…,N})，每一个种子点都会被分配一个Voronoi单元C_i,，其中C_i,定义如公式四所示：

公式四中，d(P,S_j)表示欧几里德距离。

作为优选，步骤4中，***无厚度的界面单元指的是：在将岩石试样划分为有限元网格后，在网格的实体单元之间新建采用共节点链接的界面单元，界面单元本身没有厚度，不影响岩石试样的几何特征，当界面单元的应力状态满足断裂准则后，界面单元失效并从模型中删除，可以模拟岩石的微观断裂。

发明的有益效果是：

(1)本发明提供了一种考虑岩石结构特征和矿物组成的细观数值模型生成方法，该方法不仅考虑了岩石内部的矿物成分和各向异性，还能够考虑岩石矿物颗粒含量、矿物粒度分布、晶粒尺寸、颗粒形状和晶粒的空间位置，能够更加精确的模拟岩石内部的微观结构。

(2)本发明所提供的方法可以生成有限元计算的数值模型，可以采用大型商业有限元计算软件进行力学行为的仿真模拟，例如单轴压缩试验、三轴压缩试验、单轴拉伸试验、三点弯折试验等。

(3)本发明所提供的方法计算效率较高，通用性好。

附图说明

图1是本发明实施例中涉及的数值模型生成方法的示意图；

图2是本发明实施例中涉及的岩石结构特征分析的部分示方法意图；

图中，11表示岩石试样切片，12表示切片显微观察，13表示矿物颗粒边界识别，14表示颗粒几何信息提取，15表示颗粒几何信息计算；

图3是本发明实施例中通过岩石结构特征分析得到了矿物颗粒尺寸分布规律示意图；

图4是本发明实施例中采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程示意图；其中，24表示圆球投放初期示意图，25表示圆球膨胀中期示意图，26表示圆球膨胀完成，达到密实分布示意图；

图5是本发明实施例中在膨胀过程中统计圆球堆叠***的径向分布函数示意图；

图6是本发明实施例中在膨胀过程中统计圆球堆叠***的键取向序参数示意图；

图7是本发明实施例中涉及的Voronoi图的空间区域划分方法示意图；

其中，31表示得到的圆球空间位置和尺寸分布示意图，32表示进行基于Voronoi图的空间区域划分，生成数值模型示意图；

图8是本发明实施例中涉及的生成有限元模型方法示意图；

其中，32表示本发明实施例中生成的数值模型示意图，41表示将数值模型进行网格划分得到网格模型示意图，42表示在矿物颗粒的边界***零厚度的界面单元示意图；

图9是本发明实施例中算法流程图。

具体实施方式

以下参照附图对本发明所涉及的一种考虑岩石结构特征和矿物组成的细观数值模型生成方法作详细阐述。以下实施例中未详细阐述的部分均属于现有技术。

<实施例>

如图1所示，本发明所涉及的考虑岩石结构特征和矿物组成的细观数值模型生成方法可分为4个流程：岩石结构特征分析1、粒径膨胀法2、基于Voronoi图的空间区域划分3以及生成有限元模型4。

如图1和2所示，岩石结构特征分析1具体可包括：岩石试样切片11、切片显微观察12、矿物颗粒边界识别13、颗粒几何信息提取14以及颗粒几何信息计算15。本实施例中，岩石试样切片11取直径5cm×高度10cm的花岗岩切片得到，在光学显微镜下得到切片的放大图片，采用开源软件ImageJ进行边界识别得到边界轮廓，颗粒几何信息计算15包括：矿物颗粒的尺寸(如公式五所示)、矿物颗粒的实心度(公式六所示)、矿物颗粒的长短轴之比(如公式七所示)以及矿物颗粒的圆度(如表达式八所示)。

公式五中，P_grain表示矿物颗粒的周长，π表示圆周率，R即表示矿物颗粒的尺寸。

公式六，S_grain表示矿物颗粒的面积，S_convex表示用面积最小的凸多边形包围矿物颗粒时，凸多边形的面积，Grain Solidity表示矿物颗粒的实心度。

公式七中，L_Major表示用椭圆形拟合矿物颗粒边界时，该椭圆形长轴长，L_Minor即为该椭圆形的短轴长，Aspect ratio表示矿物颗粒的长短轴之比。

公式八中，S_grain表示矿物颗粒的面积，L_Major表示用椭圆形拟合矿物颗粒边界时，该椭圆形长轴长，π表示圆周率，Grain Roundness表示矿物颗粒的圆度。

如图1、图3、图4、图5和图6所示，通过岩石结构特征分析得到了矿物颗粒尺寸分布规律21之后，采用开源软件LIGGGHTS按照该尺寸分布规律在圆柱形边界22中投放圆球 23，圆球在投放初期24呈现随机分布规律，随后圆球开始以相同的速率膨胀，膨胀中期25 的分布规律表现为圆球的堆叠逐渐开始密实，最终达到一种较为密实分布26，在膨胀过程中统计圆球堆叠***的径向分布函数27和键取向序参数28，这两个参数可以衡量整个***是否已经满足真实岩石试样中矿物颗粒的尺寸分布。本实例中，观察径向分布函数27中代表经过粒径膨胀之后的径向分布函数曲线(蓝色曲线)的峰值是否与代表真实岩石式样的径向分布函数曲线(黑色曲线)第一峰值吻合，径向分布函数中的第一个峰值表示整个颗粒体系的有序程度，第一个峰值吻合就表明数值方法生成的颗粒集合体与真实岩石中矿物颗粒集合体排列规则程度上是相吻合的；同时观察键取向序参数28中代表经过粒径膨胀之后的键取向序参数(黑色散点)是否接近HCP(体心立方排列结构)与FCC(面心立方排列结构)的交点。当上述参数满足要求时，粒径膨胀这一过程即可结束。

如图1和7所示，将得到的圆球空间位置和尺寸分布31提取之后，进行基于Voronoi(泰森多边形)图的空间区域划分，生成数值模型32。

如图1和8所示，为了方便后续进行力学性能仿真模拟，本方法将生成的数值模型32进行网格划分，划分为可用于有限元软件(比如ANSYS软件)计算的网格模型41，为了保证计算精度其中每个矿物颗粒45都会被离散成大约100个四面体单元，并在矿物颗粒的边界***零厚度的界面单元42，同时为了模拟矿物颗粒内部的开裂情况，在矿物颗粒内部也会***界面单元44。

如图9所示，本方法首先通过显微观察岩石试样切片，对岩石试样切片进行数字图像处理，对二值化的图像采用分水岭算法分割出矿物颗粒，进而得到矿物颗粒的尺寸、形态、矿物类型等信息；将不规则形状的矿物颗粒简化为一尺寸缩小的圆球，采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程，受试样边界限制和颗粒间互不侵入条件的约束，颗粒在生长过程中发生重排列，趋于一个较为密实和各向同性结构；判断颗粒粒径是否达到矿物颗粒尺寸分布，若不满足要求则继续进行粒径膨胀；当满足要求时停止离散元数值模拟，提取颗粒空间位置和尺寸等信息，对颗粒集合体进行基于Voronoi图的空间区域划分，每个Voronoi 元胞对应一个不规则形状的矿物晶粒；根据岩石中不同矿物类型的体积分数，给每个Voronoi 元胞代表的矿物晶粒分配矿物类型；将得到的Voronoi多晶体结构进行有限元网格划分，并在矿物颗粒边界和颗粒内部***无厚度界面单元，生成可反映岩石结构和岩相结构的有限元数值模型。

以上实施例仅是对本发明技术方案所做的举例说明。本发明所涉及的一种多晶岩石微观数值模型的生成方法并不仅仅限定于在以上实施例中所描述的过程，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明所要求保护的范围内。

Claims (8)

1.一种考虑岩石结构特征和矿物组成的细观数值模型生成方法，其特征在于,包括以下步骤：

步骤1，通过岩石结构特征分析，得到矿物颗粒分布和颗粒几何信息：

步骤2，离散元模拟；

将由步骤1中将岩石结构特征分析得到的矿物颗粒分布和颗粒几何信息收集汇总，将不规则形状的矿物颗粒简化为一尺寸缩小的圆球，采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程，直至颗粒体系和真实岩石试样具有相同的结构特性即颗粒粒径是否达到矿物颗粒尺寸分布时，停止颗粒粒径膨胀及离散元数值模拟，完成岩石试样建模；

步骤3，生成数值模拟；

提取岩石试样模型中颗粒空间位置和尺寸信息，对颗粒集合体进行基于Voronoi图的空间区域划分，每个Voronoi元胞对应一个不规则形状的矿物晶粒；根据岩石中不同矿物类型的体积分数，给每个Voronoi元胞代表的矿物晶粒分配矿物类型；

步骤4，生成有限元模型；

将得到的Voronoi多晶体结构进行有限元网格划分，并在矿物颗粒边界和颗粒内部***无厚度界面单元，生成可反映岩石结构和岩相结构的有限元数值模型，完成细观数值模型生成；

步骤2中，采用离散单元法模拟岩石试样内颗粒粒径逐渐膨胀的过程中，颗粒在给定的边界内开始按照相同的速率体积膨胀，在膨胀过程中，受试样边界限制和颗粒间互不侵入条件的约束，颗粒在生长过程中发生重排列，趋于一个密实和各向同性结构，上述采用离散单元法模拟岩石试样的颗粒集合体的结构特性可以采用径向分布函数和键取向序参数度量，通过这两个参数判断颗粒粒径是否达到矿物颗粒尺寸分布，若不满足要求则继续进行粒径膨胀；当满足要求时停止离散元数值模拟。

2.根据权利要求1所述的细观数值模型生成方法，其特征在于：步骤1中，岩石结构特征分析具体方法为：

首先将岩石切片，通过显微观察得到岩石试样切片图，然后对岩石试样切片图进行数字图像处理，对二值化的图像采用分水岭算法分割出矿物颗粒，识别矿物颗粒的分布和颗粒几何信息。

3.根据权利要求2所述的细观数值模型生成方法，其特征在于：步骤1中，所述矿物颗粒分布信息包括矿物颗粒的矿物种类、不同矿物种类的含量、粒径大小、粒径分布、颗粒的形心位置、圆度以及面积，颗粒几何信息包括矿物颗粒的尺寸、实心度、长短轴之比以及圆度。

4.根据权利要求3所述的细观数值模型生成方法，其特征在于：所述圆度通过下面公式定义：

公式一中，S_grain是通过岩石结构特征分析得到的矿物颗粒面积，P_grain是矿物颗粒周长，π为圆周率。

5.根据权利要求1所述的细观数值模型生成方法，其特征在于：所述径向分布函数为：

公式二中，dn(r)表明在距离原点颗粒相距r处，且在宽度dr范围内所能找到的颗粒数，N为岩石试样中总矿物颗粒数，V为岩石试样建模的体积，g(r)的物理意义是与原点颗粒相距r处单位体积的颗粒数密度，从而获得颗粒体系的结构信息，π为圆周率。

6.根据权利要求1所述的细观数值模型生成方法，其特征在于：所述键取向序参数函数为：

公式三中，

是球谐函数，其中l表示对称性的序参量，m取-l≤m≤l，

即为颗粒接触对的方向向量，而

和

为该方向向量的极角和方位角，π为圆周率，Q_l为键取向序参数值。

7.根据权利要求1所述的细观数值模型生成方法，其特征在于：步骤3中，对颗粒集合体进行基于Voronoi图的空间区域划分指的是一种在三维空间内分区的技术，对于D∈R³的区域内，R指的是长度量纲，R³表示在三维空间内给定的D中有若干种子点，D中种子点的坐标由粒径膨胀法获得，对于上述种子点i，满足{S_i(x_i)}(for i＝{1,…,N})，每一个种子点都会被分配一个Voronoi单元C_i，其中C_i定义如公式四所示：

公式四中，d(P,S_j)表示欧几里德距离。

8.根据权利要求1所述的细观数值模型生成方法，其特征在于：步骤4中，***无厚度的界面单元指的是：在将岩石试样划分为有限元网格后，在网格的实体单元之间新建采用共节点链接的界面单元，界面单元本身没有厚度，不影响岩石试样的几何特征，当界面单元的应力状态满足断裂准则后，界面单元失效并从模型中删除，可以模拟岩石的微观断裂。

Images