本发明提供可避免畸形计算单元的岩石细观数值模型生成方法及系统，方法包括：步骤1.通过数值图像处理得到岩石试样的显微图像特征，在空间域内进行基于点饱和理论的随机撒点，提取矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，模拟矿物颗粒不规则形状的细观几何结构；步骤2.根据多边形的几何信息计算其质心坐标，采用质心迭代法，重新进行空间剖分，直到满足设定的迭代停止标准；将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型；步骤3.对模型进行有限元网格划分，并在颗粒边界和内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型。
1.可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于，包括以下步骤：步骤1.岩石细观特征提取和量化
通过数值图像处理技术得到岩石试样的显微图像特征，在空间域内进行基于点饱和理论的随机撒点，提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构；
将由步骤1所得Voronoi元胞的几何信息收集汇总计算其质心坐标，采用质心迭代法，重新进行基于Voronoi图的空间剖分；判断结果是否满足迭代停止标准，若不是则继续采用质心迭代法进行优化，直到满足设定的迭代停止标准；遍历所有的矿物颗粒，将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型；
质心迭代法指的是将下次Voronoi图的生成元点移至当前Voronoi元胞质心处，重新进行Voronoi剖分；矿物颗粒在质心迭代过程中不断重排列，趋于一个较为密实和各向同性结构，受岩石试样边界、颗粒间互不侵入和密排条件的约束，Voronoi元胞代表的矿物颗粒在迭代过程中发生重排列，形状畸形的矿物颗粒逐渐减少，矿物颗粒的形状质量相应提高；
将得到的岩石矿物颗粒模型进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型。
2.根据权利要求1所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤1中，在空间中进行基于点饱和理论的随机撒点指的是：在空间域内，控制2
Voronoi多边形的成长点间最短距离dmin，对于二维空间R ，该成长点间最短距离通过公式3
一计算得到，对于三维空间R，成长点间最短距离通过公式二计算得到：式中，a、b、c为模型的长、宽、高；N为点的总数。
3.根据权利要求1所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：n
其中，在步骤1中，Voronoi图的定义为，在给定的空间域R内，n为欧氏空间维度，存在集合S由n个互不相等的点构成，S＝{p1,p2,…,pn}，第i个Voronoi多边形集合体V(pi)表示为：
4.根据权利要求1所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤2中，Voronoi元胞质心的计算方法为：
对于空间域R内的多面体M，设多面体的边界由三角形的并集给出；定义Ai,i＝0,…,N‑
1为多面体的N个三角形面，其顶点坐标为(ai,bi,ci)，运用散度定理，多面体M的体积V如公式四所示：式中：对每个三角形面Ai，其面积为 ni为三角形面的单位法向量，x为外法向向量面元，x·ni为常数；
再次运用散度定理，在空间域R定义{e1,e2,e3}为基向量，得到多面体M的质心c如公式六所示：2
对于空间域R内的多边形Z，其质心计算式可简写为公式七所示：式中，x为多边形Z各顶点的横坐标，y为多边形Z各顶点的纵坐标，N为多边形Z的边数。
5.根据权利要求1所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤2中，迭代停止标准是指短边边长小于阈值的颗粒数量不再变化，趋于稳定。
6.根据权利要求5所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤2中，短边替换指的是，在质心迭代优化后的岩石细观模型中，将边长小于阈值的短边替换为短边中点，同时将与短边相邻的邻边上，靠近短边一端的端点替换成短边中点，重新生成邻边。
7.根据权利要求5或6所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤2中，阈值为相应Voronoi元胞中最长边边长的十分之一。
8.根据权利要求1所述的可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于：其中，在步骤3中，插入无厚度的界面单元指的是，在将岩石细观数值试样划分为有限元网格后，在网格的实体单元之间新建采用共节点链接的不同强度的界面单元，界面单元本身没有厚度，不影响岩石试样的几何特征，当界面单元的应力状态在满足断裂准则后，界面单元失效并从模型中删除，能够模拟岩石的断裂；
界面单元的不同强度指的是，界面单元的弹性模量、剪切模量、内摩擦角、摩擦因数、断裂能、单轴抗压强度、压拉比计算参数不同，从而模拟岩体矿物颗粒晶间破裂与穿晶破裂。
9.可避免畸形计算单元的岩石细观数值模型生成系统，其特征在于，包括：岩石细观特征提取和量化部，通过数值图像处理技术得到岩石试样的显微图像特征，在空间域内进行基于点饱和理论的随机撒点，提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构；
岩石细观结构优化部，与所述岩石细观特征提取和量化部通信相连，根据Voronoi元胞的几何信息计算其质心坐标，并采用质心迭代法重新进行基于Voronoi图的空间剖分；判断结果是否满足迭代停止标准，若不是则继续采用质心迭代法进行优化，直到满足设定的迭代停止标准；遍历所有的矿物颗粒，将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型；以及岩石细观数值模型生成部，与所述岩石细观结构优化部通信相连，将得到的岩石矿物颗粒模型进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型；
其中，质心迭代法指的是将下次Voronoi图的生成元点移至当前Voronoi元胞质心处，重新进行Voronoi剖分；矿物颗粒在质心迭代过程中不断重排列，趋于一个较为密实和各向同性结构，受岩石试样边界、颗粒间互不侵入和密排条件的约束，Voronoi元胞代表的矿物颗粒在迭代过程中发生重排列，形状畸形的矿物颗粒逐渐减少，矿物颗粒的形状质量相应提高。
可避免畸形计算单元的岩石细观数值模型生成方法及系统
技术领域
[0001] 本发明属于岩石力学研究领域，具体涉及可避免畸形计算单元的岩石细观数值模型生成方法及系统。技术背景
[0002] 岩石是自然界中经地质作用形成的具有一定结构构造的矿物集合体，由于受内部矿物晶粒、孔隙及胶结材料等细观结构的影响，岩石的力学特性及变形破坏规律呈现出其独特性和复杂性。大量试验和数值模拟表明，岩石的宏观力学特性与其微观结构密切相关，矿物颗粒的尺寸、形状、互锁程度和接触类型都能影响岩石的力学特性。从矿物晶粒大小看，在岩浆岩、变质岩和沉积岩中，等粒结构一般比非等粒结构强度高；在等粒结构中，细粒结构比粗粒结构强度高。在斑状结构中，细粒基质比玻璃基质强度高；粗粒具斑晶的酸性深成岩强度最低；细粒微晶而无玻璃质的基性喷出岩强度最高。从结构粘结上看，大部分岩浆岩、变质岩及沉积岩中的化学岩，晶粒之间结合紧密，强度较高，但在化学岩中，以可溶性结晶连结的，强度较高，但抗水性差。固结粘土岩有一部分是再结晶连结，其强度比其它坚硬岩石差得多。通过高倍显微镜观察岩石切片，发现通常岩石材料组成颗粒的形状接近于不规则多边形，其中细粒砂岩和花岗岩的组成颗粒排列紧凑，胶结面薄，其切片微观形状接近一种较为密实和各向同性结构的多边形网格拓扑结构。
[0003] 国内外的很多学者依次对岩石细观模型的生成开展了大量的研究。现有的岩石细观数值研究，或通过数字图像处理技术重构数值岩石试验，或采用细观数值模型随机生成算法。这两类生成岩石细观数值模型的方法，会产生大量的畸形单元体，影响了有限元等数值模拟方法的收敛性。2014年，Ghazvinian等人提出了一种三维随机Voronoi晶粒模型用于模拟脆性岩石的裂纹扩展，通过Voronoi镶嵌获得晶体模型，并用晶体边界模拟晶粒间的裂纹。2016年，Gui等人提出了一种二维Voronoi晶粒模型模拟了Gosford砂岩的微观开裂模型，并通过正则化迭代改善了晶体颗粒结构。2016年，Cantor等人借助Neper软件，提出一种可用于模拟颗粒材料开裂的三维中心Voronoi镶嵌颗粒模型生成方式。2017年，Li等人提出一种考虑颗粒大小分布和孔隙的二维随机Voronoi晶粒模型，模拟了Cottaer砂岩的微观开裂模型。2019年，专利CN110069844A公开了一种考虑岩石结构特征和矿物组成的三维Voronoi晶粒模型，采用粒径膨胀法能更真实的反映岩石试样的结构特性。
[0004] 但是，上述岩石细观模型生成方法，矿物颗粒数量一般需控制在200～2000之间，导致计算精度受限；且畸形的矿物颗粒较多，在网格划分过程易产生畸形单元和过密网格，进而导致计算效率和计算精度下降，甚至出现仿真模拟不收敛无法继续的现象。另一方面，上述大多方法仅考虑岩石颗粒晶体间开裂的情况，没有考虑岩石矿物晶粒穿晶开裂现象，因此需要建立一种避免畸形矿物颗粒，计算效率和计算精度较高，能够考虑岩石细观结构下真实破碎情况的细观数值模型。
发明内容
[0005] 本发明是为了解决上述问题而进行的，目的在于提供一种可避免畸形计算单元的岩石细观数值模型生成方法及系统。
[0006] 本发明为了实现上述目的，采用了以下方案：
[0007] <方法>
[0008] 本发明提供一种可避免畸形计算单元的岩石细观数值模型生成方法，其特征在于，包括：
[0009] 步骤1.岩石细观特征提取和量化
[0010] 通过数值图像处理技术得到岩石试样的显微图像特征，将不规则形状的矿物颗粒简化成一点，在空间域内进行基于点饱和理论的随机撒点，提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构；
[0011] 步骤2.岩石细观结构优化
[0012] 将由步骤1所得Voronoi元胞的几何信息收集汇总计算其质心坐标，采用质心迭代法，重新进行基于Voronoi图的空间剖分；判断结果是否满足迭代停止标准，若不是则继续采用质心迭代法进行优化，直到满足设定的迭代停止标准；遍历所有的矿物颗粒，将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型；
[0013] 步骤3.岩石细观数值模型生成
[0014] 将得到的岩石矿物颗粒模型进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型。
[0015] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤1中，基于Voronoi图的空间区域剖分指的是一种在空间区域内分区的技术，Voronoi图的定义为：在空间域内，控制Voronoi多边形的成长点间最短距离dmin，对2 3
于二维空间R ，该成长点间最短距离通过公式一计算得到，对于三维空间R ，成长点间最短距离通过公式二计算得到：
[0016]
[0017]
[0018] 式中，a、b、c为模型的长、宽、高；N为点的总数。
[0019] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以n具有以下特征：在步骤1中，Voronoi图的定义为，在给定的空间域R内，n为欧氏空间维度，存在集合S由n个互不相等的点构成，S＝{p1，p2，...，pn}，
[0020] 第i个Voronoi多边形集合体V(pi)表示为：
[0021] V(pi)＝∩i≠j{p|d(p，pi)＜dp，pj，i，j＝{1，2，…，n}} (公式三)[0022] 式中，d(p，pi)表示点p与点pi间的欧氏距离；V(pi)为一个凸多边形，被称为点pi的Voronoi结构，pi为第i个Voronoi多边形集合体V(pi)的成长点或者核；点集合S为Voronoi图的生成元。
[0023] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：
[0024] 在步骤2中，Voronoi元胞质心的计算方法为：
[0025] 对于空间域R3内的多面体M，设多面体的边界由三角形的并集给出；定义Ai，i＝0，...，N‑1为多面体的N个三角形面，其顶点坐标为(ai，bi，ci)，运用散度定理，多面体M的体积V如公式四所示：
[0026]
[0027]
[0028] 式中：对每个三角形面Ai，其面积为 ni为三角形面的单位法向量，x为外法向向量面元，x·ni为常数；
[0029] 对多面体M的质心c∈R3，由质心定义有：
[0030]
[0031] 再次运用散度定理，在空间域R3定义{e1，e2，e3}为基向量，得到多面体M的质心c如公式六所示：
[0032]
[0033]
[0034] 对于空间域R2内的多边形Z，其质心计算式可简写为公式七所示：
[0035]
[0036]
[0037] 式中，x为多边形Z各顶点的横坐标，y为多边形Z各顶点的纵坐标，N为多边形Z的边数。
[0038] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤2中，质心迭代法指的是将下次Voronoi图的生成元点移至当前Voronoi元胞质心处，重新进行Voronoi剖分；矿物颗粒在质心迭代过程中不断重排列，趋于一个较为密实和各向同性结构，受岩石试样边界、颗粒间互不侵入和密排等条件的约束，Voronoi元胞代表的矿物颗粒在迭代过程中发生重排列，形状畸形的矿物颗粒逐渐减少，矿物颗粒的形状质量相应提高。
[0039] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤2中，迭代停止标准是指短边边长小于阈值的颗粒数量不再变化，趋于稳定。
[0040] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤2中，短边替换指的是，在质心迭代优化后的岩石细观模型中，将边长小于阈值的短边替换为短边中点，同时将与短边相邻的邻边上，靠近短边一端的端点替换成短边中点，重新生成邻边。
[0041] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤2中，阈值为相应Voronoi元胞中最长边边长的十分之一。
[0042] 优选地，本发明提供的可避免畸形计算单元的岩石细观数值模型生成方法还可以具有以下特征：在步骤3中，插入无厚度的界面单元指的是，在将岩石细观数值试样划分为有限元网格后，在网格的实体单元之间新建采用共节点链接的不同强度的界面单元，界面单元本身没有厚度，不影响岩石试样的几何特征，当界面单元的应力状态在满足断裂准则后，界面单元失效并从模型中删除，能够模拟岩石的断裂；界面单元的不同强度指的是，界面单元的弹性模量、剪切模量、内摩擦角、摩擦因数、断裂能、单轴抗压强度、压拉比等计算参数不同，从而模拟岩体矿物颗粒晶间破裂与穿晶破裂。
[0043] <系统>
[0044] 进一步，本发明还提供了一种可避免畸形计算单元的岩石细观数值模型生成系统，其特征在于，包括：
[0045] 岩石细观特征提取和量化部，通过数值图像处理技术得到岩石试样的显微图像特征，在空间域内进行基于点饱和理论的随机撒点，提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构；
[0046] 岩石细观结构优化部，与岩石细观特征提取和量化部通信相连，根据岩石细观特征提取和量化部得到的Voronoi元胞的几何信息计算其质心坐标，采用质心迭代法，重新进行基于Voronoi图的空间剖分；判断结果是否满足迭代停止标准，若不是则继续采用质心迭代法进行优化，直到满足设定的迭代停止标准；遍历所有的矿物颗粒，将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型；以及
[0047] 岩石细观数值模型生成部，与岩石细观结构优化部通信相连，将得到的岩石矿物颗粒模型进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型。
[0048] 发明的作用与效果
[0049] 本发明所提供的可避免畸形计算单元的岩石细观数值模型生成方法及系统，通过基于点饱和理论的随机撒点、质心迭代法和短边替换法，模拟出自然界岩石在地壳运动不断挤压趋于一个较为密实和各向同性结构的过程，优化岩石细观数值模型，生成高质量的有限元网格进行开裂模拟，能够非常有效地避免出现畸形矿物单元，可高效、精确地生成相应的二维岩石细观模型和三维岩石细观模型，生成的岩石细观模型能够真实地反映岩石细观结构和开裂行为，并且可以采用连续离散耦合分析方法进行数值材料试验，例如单轴压缩试验、三轴压缩试验、单轴拉伸试验、三点弯折试验、巴西圆盘劈裂试验等。
附图说明
[0050] 图1为本发明实施例中涉及的可避免畸形计算单元的岩石细观数值模型生成方法的流程图；
[0051] 图2为本发明实施例中涉及的岩石的细观结构特征和矿物颗粒开裂示意图；其中，11，12分别为岩石试样切片显微观察下，矿物颗粒不规则形状的细观几何结构示意图，13为岩石破碎中的晶间破裂示意图，14为岩石破碎中的穿晶破裂示意图；
[0052] 图3为本发明实施例中涉及的基于点饱和理论的Voronoi图三维空间区域剖分过程示意图；其中，21为基于点饱和理论的随机撒点示意图，22为基于Voronoi图的空间区域剖分示意图；
[0053] 图4为本发明实施例中涉及的质心迭代方法二维示意图；其中，31为矿物颗粒质心计算示意图，32为生成元点移至Voronoi元胞质心示意图，33为质心迭代过程中短边数量变化规律示意图，34为质心迭代完成、达到迭代停止标准的岩石矿物颗粒模型示意图；
[0054] 图5为本发明实施例中涉及的短边替换方法二维示意图；其中，41为包含短边的Voronoi元胞短边替换重生成示意图，42为短边替换完成后生成的岩石矿物颗粒模型示意图；
[0055] 图6为本发明实施例中涉及的二维有限元模型生成过程示意图；其中，42为本发明实例中生成的岩石矿物颗粒模型示意图(数值模型示意图)，51为将数值模型拉伸为单位厚度、便于后期插入界面单元的假三维体结构示意图，52为将数值模型进行网格划分得到的网格模型示意图，53为在矿物颗粒的边界插入零厚度的界面单元示意图，54为在矿物颗粒的内部插入零厚度的界面单元示意图；
[0056] 图7为本发明实施例中涉及的算法流程图；
[0057] 图8为本发明对比例中涉及的本发明方法与专利CN110069844A方法的结果对比图；其中，61为采用本发明方法生成岩石细观模型，62为采用专利CN110069844A方法生成岩石细观模型，63为采用本发明与专利CN110069844A方法得到的岩石细观数值模型的矿物颗粒长宽比累计分布对比图。
具体实施方式
[0058] 以下结合附图对本发明涉及的可避免畸形计算单元的岩石细观数值模型生成方法及系统的具体实施方案进行详细地说明。
[0059] <实施例>
[0060] 如图1所示，本实施例中所提供的生可避免畸形计算单元的岩石细观数值模型生成方法包括以下步骤：
[0061] 步骤1.岩石细观特征提取和量化
[0062] 通过数值图像处理技术得到岩石试样的显微图像特征，在空间域内进行基于点饱和理论的随机撒点，提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构。
[0063] 本实施例中，以包括矿物颗粒的晶间破裂和穿晶破裂的Carrara大理石切片11和小浪底砂岩切片12为例，进行特征提取和数值模拟。考虑岩石的细观结构特征和矿物颗粒破裂特征具体包括，在光学显微镜下得到Carrara大理石切片11、小浪底砂岩切片12的放大图片中，矿物颗粒不规则形状的细观形状结构，以及岩石破碎过程中，矿物颗粒的晶间破裂13和穿晶破裂14。
[0064] 如图1和3所示，通过岩石结构特征分析得到了矿物颗粒细观结构特征，按照该形状特征，通过商业数学软件MATLAB，进行基于点饱和理论的随机撒点21，采用Voronoi图的空间区域剖分22来模拟岩石矿物颗粒不规则多边形形状。
[0065] 步骤2.岩石细观结构优化
[0066] 如图1和4所示，将得到的Voronoi矿物颗粒元胞质心位置31提取之后，将每个矿物颗粒的生成点位置移至多边形质心32，不断进行质心迭代重新生成Voronoi多边形，达到最优迭代次数33后停止迭代，得到一个较为密实和各向同性的结构34。
[0067] 如图1和5所示，扫描颗粒体系，将Voronoi元胞中边长小于某阈值的短边替换为短边中点41，重新生成易划分优质网格的Voronoi元胞42。
[0068] 步骤3.岩石细观数值模型生成
[0069] 如图1和6所示，为了方便后续进行力学性能仿真模拟，本方法将生成的Voronoi元胞43，拉伸为体状数值模型51，划分为可用于有限元软件(比如ANSYS软件)计算的网格模型52，为了保证计算精度，其中每个矿物颗粒都会被离散成大约6‑10个三棱柱单元，并在矿物颗粒的边界插入零厚度的界面单元53，同时为了模拟矿物颗粒内部的穿晶破裂，在矿物颗粒内部也会插入界面单元54。
[0070] 如图7所示，本方法首先通过数值图像处理技术，得到岩石细观结构特征及矿物颗粒内部破碎特征；通过基于点饱和理论的Voronoi空间区域剖分，模拟岩石结构特征；采用质心迭代法模拟自然界岩石在地壳运动不断挤压的过程，受试样边界限制和颗粒间互不侵入条件的约束，颗粒在生长过程中发生重排列，趋于一个较为密实和各向同性结构；采用短边替换法重新生成易划分优质网格的Voronoi元胞；将得到的Voronoi多晶体结构进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入不同强度的无厚度界面单元，生成可反映真实岩石细观结构特征和矿物颗粒破裂特征的有限元数值模型。
[0071] 进一步，本实施例还提供了一种能够自动实现以上可避免畸形计算单元的岩石细观数值模型生成方法的系统。具体地，该系统包括：岩石细观特征提取和量化部、岩石细观结构优化部、岩石细观数值模型生成部、输入显示部以及控制部。
[0072] 岩石细观特征提取和量化部按照步骤1所描述的方法采用数值图像处理技术对岩石试样进行处理，得到岩石试样的显微图像特征，并在空间域内进行基于点饱和理论的随机撒点，然后提取岩石试样模型中矿物颗粒空间位置，对颗粒集合体在岩石试样边界内进行基于Voronoi图的空间区域剖分，生成Voronoi多边形，每个Voronoi多边形为一个Voronoi元胞，每个Voronoi元胞对应一个形状不规则的矿物颗粒，模拟矿物颗粒不规则形状的细观几何结构。
[0073] 岩石细观结构优化部与岩石细观特征提取和量化部通信相连，按照步骤2所描述的方法，根据岩石细观特征提取和量化部得到的Voronoi元胞的几何信息计算其质心坐标，采用质心迭代法，重新进行基于Voronoi图的空间剖分；判断结果是否满足迭代停止标准，若不是则继续采用质心迭代法进行优化，直到满足设定的迭代停止标准；遍历所有的矿物颗粒，将矿物颗粒的所有短边替换为短边中点，并重新生成岩石矿物颗粒模型。
[0074] 岩石细观数值模型生成部与岩石细观结构优化部通信相连，按照步骤3所描述的方法，将得到的岩石矿物颗粒模型进行有限元网格划分，并在矿物颗粒边界和颗粒内部插入无厚度界面单元，生成能够模拟岩石开裂的有限元数值模型。
[0075] 输入显示部与岩石细观特征提取和量化部、岩石细观结构优化部、岩石细观数值模型生成部均通信相连，根据根据输入的指令，对岩石细观特征提取和量化部、岩石细观结构优化部、岩石细观数值模型生成部的数据进行相应显示。
[0076] 控制部与岩石细观特征提取和量化部、岩石细观结构优化部、岩石细观数值模型生成部、输入显示部均通信相连，控制它们的运行。
[0077] <对比例>
[0078] 如图8所示，采用本发明方法生成岩石细观模型61，采用专利CN110069844A方法生成岩石细观模型62，岩石试样长度为100mm，直径为50mm，包含2000个矿物颗粒。扫描整个颗粒体系，得到每个矿物颗粒最长边与最短边之比(定义为长宽比)，绘制岩石细观数值模型的矿物颗粒长宽比累计分布图63，和专利CN110069844A相比，本发明方法生成岩石试样的矿物颗粒长宽比均在65以下，更利于划分高质量的有限元网格，提高计算效率，保证计算的收敛性。
[0079] 如下表1所示，分别对本发明和专利CN110069844A生成的岩石试样进行控制单元边长度的网格划分，定义长宽比大于20、最大内角超过165度的单元为畸形单元，统计得到不同单元线长下岩石细观模型的畸形单元数目,并在网格的实体单元之间插入可模拟破碎的界面单元，进行单轴压缩试验，控制加载条件相同，计算均采用24个Intel Xeon2.4GHz处理器和16GB DDR31600MHZ RAM内存的硬件条件，统计得到各模型的计算时间，与专利CN110069844A相比，本发明生成岩石试样的细观数值模型畸形计算单元显著减少，计算收敛性更佳，计算效率更好。
[0080] 表1
[0081]
[0082] 以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的可避免畸形计算单元的岩石细观数值模型生成方法及系统并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。
