CN115964901B - 基于离散单元法的水致岩石强度劣化的模拟方法及系统 - Google Patents

Abstract

本发明涉及一种基于离散单元法的水致岩石强度劣化的模拟方法及系统，涉及水致岩石强度劣化的模拟领域，包括：确定岩石试样矿物类型、含量信息和模拟块体参数；制备不同饱和度下的岩石试样；确定岩石试样的宏观力学参数；参数校准；对校准后的参数设置梯度；确定不同接触法向刚度、不同接触粘聚力、不同接触内摩擦角以及不同接触抗拉强度下数值模型的杨氏模量、单轴抗压强度、巴西抗拉强度、接触粘聚力以及接触内摩擦角；确定各类模拟参数与模拟所得宏观力学参数的关系；确定各模拟参数的预测值；确定数值模型宏观力学参数；对各模拟参数的预测值进行微调。本发明中的上述方案能够实现水对岩石强度劣化的数值模拟表征，提高模拟精度以及计算效率。

Description

基于离散单元法的水致岩石强度劣化的模拟方法及系统

技术领域

本发明涉及水致岩石强度劣化的模拟领域，特别是涉及一种基于离散单元法的水致岩石强度劣化的模拟方法及系统。

背景技术

岩石浸水后发生软化，强度降低，更易发生变形破坏，是岩体工程地质灾害的主要原因之一。数值模拟可基于岩体结构特性，再现内外因素作用下岩体工程地质体响应过程，为岩体工程地质灾害防治提供重要依据。因此，实现水致岩石强度劣化数值模拟对岩体工程实践具有重要意义。

在众多数值模拟手段中，离散单元法可较为真实地模拟节理岩体的非线性大变形特征，受到了岩体工程相关从业者的广泛青睐。在离散单元法中，岩体这类不连续介质由刚体或可变形的颗粒组合而成，颗粒之间通过接触粘结。用户通过定义颗粒和接触的属性及本构方程，可求解力与相对位移等物理量，计算得到颗粒的平移、转动等行为和接触的分离、滑移、压缩，进而实现岩体的非线性力学行为模拟，因此被广泛应用于滑坡、采矿、桥隧等的数值模拟。

得益于离散单元法在模拟岩体力学行为的显著优势，现已涌现出一些基于离散单元法的水致岩石强度劣化模拟方法，根据其原理可分为三类：

（1）假定模拟参数与岩石宏观力学参数等效（Wang等，2019）。诸多岩石吸水性试验结果表明，水对岩石强度劣化体现为岩石宏观力学参数的降低，具体包括内摩擦角、粘聚力、杨氏模量、单轴压缩强度随岩石饱和度增加而减小。通过假设模拟参数与上述岩石宏观力学参数等效，则可根据吸水性试验结果得到模拟参数与饱和度的关系，从而实现水致岩石强度劣化的模拟。然而，模拟参数与岩石宏观力学参数在物理意义和尺度上存在明显差异，因此二者间的等效性在原理上是错误的。

（2）模拟水致岩石强度劣化机理（Gu等，2020）。已有研究表明，水致岩石强度劣化是由于水岩接触发生力学、化学、物理作用下岩石内部胶结物溶解。借助离散单元法的基本原理，可通过假定颗粒间接触的参数遇水后降低，具体包括接触刚度、接触内摩擦角、接触抗拉强度、接触粘聚力随饱和度增加而减小，进而实现水致岩石强度劣化的模拟。然而，数值模拟参数受岩石宏观力学参数约束而非饱和度，因此直接构建数值模拟参数和饱和度的关系忽略了这一基本原理，所得结果真实性差。

（3）根据不同饱和度下岩石宏观力学参数对模拟参数进行标定（骆祚森等，2019）。根据离散单元法中模拟参数受岩石宏观力学参数约束的原理，可将试验所得不同饱和度下岩石宏观力学参数作为条件，通过传统的试错法不断调整模拟参数，包括块体弹性模量、块体泊松比、接触刚度、接触内摩擦角、接触抗拉强度、接触粘聚力等，直至模拟结果与各组试验结果相一致，最终确定各组模拟参数，从而表征不同饱和度下的岩石力学行为。该方法要求对模拟参数不断试错，这极大占用了计算资源。

针对以上传统的水致岩石强度劣化模拟方法存在原理不清、真实性差、计算繁琐的问题，本发明旨在提出一种新的基于离散单元法的水致岩石强度劣化的模拟方法及系统，简单可行地实现水对岩石强度劣化的数值模拟表征。

发明内容

本发明的目的是提供一种基于离散单元法的水致岩石强度劣化的模拟方法及系统，实现水对岩石强度劣化的数值模拟表征，提高模拟精度以及计算效率。

为实现上述目的，本发明提供了如下方案：

第一方面，本发明提供一种基于离散单元法的水致岩石强度劣化的模拟方法，所述模拟方法包括：

S1：利用X射线衍射仪确定岩石试样矿物类型和含量信息；

S2：根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block；

S3：制备不同饱和度下的岩石试样；

S4：利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

S5：选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont；

S6：分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度；

S7：分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

S8：通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系；

S9：将S4中所得不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值；

S10：将S9中所得不同饱和度下的各模拟参数的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，得到数值模型宏观力学参数；

S11：对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，得到对比结果；

S12：根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致；

S13：基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

可选的，所述选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数具体包括以下步骤：

S5.1：输入任意一组模拟参数，包括接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont，接触内摩擦角、接触抗拉强度σ_t ^cont；

S5.2：判断接触法向刚度和接触剪切刚度是否均小于或等于预设阈值；

S5.3：若否，则返回S5.1；

S5.4：若是，则执行下一步骤；

S5.5：判断接触剪切刚度和接触法向刚度的比值是否与试样剪切模量和试样杨氏模量的比值相等；

S5.6：若否，则返回S5.1；

S5.7：若是，则执行下一步骤；

S5.8：在离散单元法软件中进行单轴压缩数值模拟；数值模型尺寸和边界条件与单轴压缩试样和加载条件相同；

S5.9：判断单轴压缩数值模拟所得数值模型的杨氏模量是否与试验值相等；

S5.10：若否，则返回S5.1；

S5.11：若是，则确定接触剪切刚度和接触法向刚度的值，在后续步骤中即使返回S5.1仍保持接触剪切刚度和接触法向刚度值不变，并执行下一步骤；

S5.12：在离散单元法软件中进行巴西劈裂数值模拟；数值模型尺寸和边界条件与巴西劈裂试样和加载条件相同；

S5.13：判断巴西劈裂数值模拟所得的数值模型巴西抗拉强度是否与试验值相等；

S5.14：若否，则返回S5.1；

S5.15：若是，则确定接触抗拉强度的值，在后续步骤中即使返回S5.1仍保持接触抗拉强度值不变，并执行下一步骤；

S5.16：在离散单元法软件中进行单轴压缩数值模拟和三轴压缩数值模拟；数值模型尺寸和边界条件与用于单轴压缩试验和三轴压缩试验的试样和加载条件相同；

S5.17：判断单轴压缩数值模拟和三轴压缩数值模拟所得数值模型的粘聚力和内摩擦角是否与试验值相等；

S5.18：若否，则返回S5.1；

S5.19：若是，则确定接触粘聚力和接触内摩擦角，试错法结束，用于模拟选定饱和度下试样力学行为的模拟参数校准完成。

可选的，所述制备不同饱和度下的岩石试样具体包括：

依据国际规范制备饱和岩石试样，利用烘干箱制备干燥岩石试样，在真空容器中应用质量控制法使干燥岩石试样吸水，以得到不同饱和度下的岩石试样。

可选的，所述模拟块体表征矿物颗粒。

可选的，所述预设阈值的表达式如下：

其中，Δz_min表示相邻单元在垂直方向的最小宽度，max表示对于与接触相邻的所有单元的最大值，K^block表示模拟块体体积模量，G^block表示模拟块体剪切模量。

可选的，所述各类模拟参数与模拟所得宏观力学参数的关系具体包括：接触法向刚度-杨氏模量、接触抗拉强度-巴西抗拉强度、接触粘聚力-单轴抗压强度、接触粘聚力-粘聚力、接触内摩擦角-单轴抗压强度以及接触内摩擦角-内摩擦角的关系。

第二方面，基于本发明中的上述方法，本发明另外提供一种基于离散单元法的水致岩石强度劣化的模拟系统，所述系统包括：

岩石试样矿物类型和含量信息确定模块，用于利用X射线衍射仪确定岩石试样矿物类型和含量信息；

模拟块体参数确定模块，用于根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block；

不同饱和度下的岩石试样制备模块，用于制备不同饱和度下的岩石试样；

不同饱和度岩石试样的宏观力学参数确定模块，用于利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

校准模块，用于选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括：接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont；

梯度设置模块，用于分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度；

第一数值模拟模块，用于分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

线性拟合模块，用于通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系；

预测模块，用于将不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值；

第二数值模拟模块，用于将所得不同饱和度下的各模拟参数的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，得到数值模型宏观力学参数；

比对模块，用于对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，得到对比结果；

调整模块，用于根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致；

不同饱和度下的岩石工况模拟模块，用于基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

第三方面，本发明提供一种电子设备，包括存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行上述的离散单元法的水致岩石强度劣化的模拟方法。

第四方面，本发明提供一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现上述的离散单元法的水致岩石强度劣化的模拟方法。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

不同于传统方法中试图建立模拟参数与饱和度的关系，本发明立足于“模拟参数受岩石宏观力学参数约束”这一离散单元法基本原理，通过构建模拟参数与岩石宏观力学参数的关系，进而将不同饱和度下的岩石宏观力学参数输入至该关系，输出相应的预测模拟参数并进行微调，最终确定不同饱和度下的模拟参数，从而实现水致岩石强度劣化的模拟。因此，该发明解决了传统的水致岩石强度劣化模拟方法原理不清的问题；得益于离散单元法在模拟岩体非线性大变形等方面的显著优势，本发明可真实再现不同饱和度下岩石的变形、破裂等力学行为；在上述离散单元法基本原理的基础上，本发明所得到的不同饱和度下模拟参数作为输入，输出模拟宏观力学参数与试验结果具有较高一致性，算例表明，本发明所得到的模拟宏观力学参数与试验结果误差在0.74%以下，应力-应变曲线和破裂模式与试验结果相似性高，印证了本发明可真实模拟水致岩石强度劣化过程；本发明借助模拟参数与岩石宏观力学参数的关系来确定不同饱和度下模拟参数，避免了以不同饱和度下岩石宏观力学参数为约束，大量重复性应用试错法对模拟参数进行标定，这极大节约了运算时间和内存占用。

附图说明

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

图1为本发明基于离散单元法的水致岩石强度劣化的模拟方法流程图；

图2为本发明单轴压缩模拟和巴西劈裂模拟数值模型示意图；

图3为本发明接触法向刚度-宏观杨氏模量参数关联示意图；

图4为本发明接触抗拉强度-巴西抗拉强度参数关联示意图；

图5为本发明接触粘聚力及接触内摩擦角-单轴抗压强度参数关联示意图；

图6为本发明试验与模拟所得宏观力学参数对比示意图；

图7为本发明试验与模拟结果对比示意图。

具体实施方式

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

本发明的目的是提供一种基于离散单元法的水致岩石强度劣化的模拟方法及系统，实现水对岩石强度劣化的数值模拟表征，提高模拟精度以及计算效率。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明基于离散单元法的水致岩石强度劣化的模拟方法流程图，如图1所述，本发明中的模拟方法包括：

S1：利用X射线衍射仪确定岩石试样矿物类型和含量信息。

S2：根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block。

其中，在离散单元法中矿物颗粒常用模拟块体表征。

S3：制备不同饱和度下的岩石试样。

具体的，是依据国际规范制备饱和岩石试样，利用烘干箱制备干燥岩石试样，在真空容器中应用质量控制法使干燥岩石试样吸水，以得到不同饱和度下的岩石试样。

S4：利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角。

其中，试验过程参考中华人民共和国行业标准《水利水电工程岩石试验规程》（SL264-2001），以得到不同饱和度岩石试样的宏观力学参数。

S5：选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括：接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont。

其中，在离散单元法软件中应用试错法校准本饱和度下试样的宏观力学参数所对应的模拟参数，具体步骤如下：

S5.1：输入任意一组模拟参数（也叫接触参数），包括接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont，接触内摩擦角、接触抗拉强度σ_t ^cont；

S5.2：判断接触法向刚度和接触剪切刚度是否均小于或等于预设阈值。

其中，预设阈值为：

其中，Δz_min表示相邻单元在垂直方向的最小宽度，max表示对于与接触相邻的所有单元的最大值，K^block表示模拟块体体积模量，G^block表示模拟块体剪切模量（即相邻接触可能存在几种材料）。

S5.3：若否，则返回S5.1；

S5.4：若是，则执行下一步骤；

S5.5：判断接触剪切刚度和接触法向刚度的比值是否与试样剪切模量和试样杨氏模量的比值相等；

S5.6：若否，则返回S5.1；

S5.7：若是，则执行下一步骤；

S5.8：在离散单元法软件中进行单轴压缩数值模拟，数值模型尺寸和边界条件与单轴压缩试样和加载条件相同；

S5.9：判断单轴压缩数值模拟所得数值模型的杨氏模量是否与试验值相等；

S5.10：若否，则返回S5.1；

S5.11：若是，则确定接触剪切刚度和接触法向刚度的值，在后续步骤中即使返回S5.1仍保持接触剪切刚度和接触法向刚度值不变，并执行下一步骤；

S5.12：在离散单元法软件中进行巴西劈裂数值模拟；数值模型尺寸和边界条件与巴西劈裂试样和加载条件相同；

S5.13：判断巴西劈裂数值模拟所得的数值模型巴西抗拉强度是否与试验值相等；

S5.14：若否，则返回S5.1；

S5.15：若是，则确定接触抗拉强度的值，在后续步骤中即使返回S5.1仍保持接触抗拉强度值不变，并执行下一步骤；

S5.16：在离散单元法软件中进行单轴压缩数值模拟和三轴压缩数值模拟，数值模型尺寸和边界条件与用于单轴压缩试验和三轴压缩试验的试样和加载条件相同；

S5.17：判断单轴压缩数值模拟和三轴压缩数值模拟所得数值模型的粘聚力和内摩擦角是否与试验值相等；

S5.18：若否，则返回S5.1；

S5.19：若是，则确定接触粘聚力和接触内摩擦角，试错法结束，用于模拟选定饱和度下试样力学行为的模拟参数校准完成。

S6：分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度。

其中，梯度设置方向应视所选定的饱和度而定，理论上应从0至预期极大值均匀设定。

S7：梯度设置完毕后，分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，数值模型尺寸和边界条件与试验相同，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、接触粘聚力c以及接触内摩擦角。

S8：通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系。

具体包括接触法向刚度-杨氏模量、接触抗拉强度-巴西抗拉强度、接触粘聚力-单轴抗压强度、接触粘聚力-粘聚力、接触内摩擦角-单轴抗压强度、接触内摩擦角-内摩擦角的关系。

S9：将S4中所得不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值。

具体的，将S4中所得不同饱和度下实际岩石的宏观力学参数（杨氏模量、单轴抗压强度、巴西抗拉强度、粘聚力、内摩擦角）输入至上述各类模拟参数与数值模型宏观力学参数的关系（接触法向刚度-杨氏模量、接触抗拉强度-巴西抗拉强度、接触粘聚力-单轴抗压强度、接触粘聚力-粘聚力、接触内摩擦角-单轴抗压强度、接触内摩擦角-内摩擦角），得到不同饱和度下各模拟参数（接触法向刚度、接触抗拉强度、接触粘聚力、接触内摩擦角）的预测值。

S10：将S9中所得不同饱和度下的各模拟参数（接触法向刚度、接触抗拉强度、接触粘聚力、接触内摩擦角）的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，数值模型尺寸和边界条件与试验相同，得到数值模型宏观力学参数。

S11：对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，包括杨氏模量、单轴抗压强度、巴西抗拉强度、粘聚力、内摩擦角，得到对比结果。

S12：根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致。

S13：基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

具体的，根据对比结果，对不同饱和度下的各预测模拟参数（接触法向刚度、接触抗拉强度、接触粘聚力、接触内摩擦角）进行微调。考虑到岩石具有非均质性，微调范围建议在线性拟合曲线的90-100%预测带与置信带内。进而重复单轴压缩数值模拟、巴西劈裂数值模拟、三轴压缩数值模拟，数值模型尺寸和边界条件与试验相同，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致，包括杨氏模量、单轴抗压强度、巴西抗拉强度、粘聚力、内摩擦角，最终可确定不同饱和度所对应的模拟参数（接触法向刚度、接触抗拉强度、接触粘聚力、接触内摩擦角），由此则可进行不同饱和度下的岩石工况模拟，例如：得到了不同饱和度对应的模拟参数，比如饱和度10%对应岩石宏观参数1和模拟参数a，饱和度50%对应岩石宏观参数2和模拟参数b，那么岩石宏观参数1到宏观参数2的变化就能用模拟参数a到模拟参数b的变化来表示，进而在后续应用时，如降雨入渗岩体导致滑坡，10%饱和度到50%的饱和度如何模拟，就可以应用所得的模拟参数a到模拟参数b的变化。

下面以一个具体的实施例对本发明作进一步解释说明：

选取湖北宜昌红砂岩为研究原型，借助X射线衍射分析方法确定砂岩矿物类型、含量、尺寸信息，据此确定模拟块体参数，包括块体密度、块体体积模量、块体剪切模量。

对岩样进行吸水试验，得到不同饱和度（9.23%，30%，60%，90%，100%）下的红砂岩试样。

利用WAW-300电液伺服刚性试验机对不同饱和度下的试样进行单轴压缩试验和巴西劈裂试验，得到各组试样的宏观力学参数，包括杨氏模量、泊松比、单轴抗压强度、巴西抗拉强度。

应用离散单元法软件UDEC（Universal Distinct Element Code）对试样进行建模，并选择天然饱和度（9.23%）下试样宏观力学参数作为约束，开展单轴压缩数值模拟和巴西劈裂数值模拟（见图2），应用试错法不断调整模拟参数，直至模拟所得宏观力学参数与试验所得结果一致。

基于上述天然饱和度对应的模拟参数，将接触法向刚度、接触粘聚力、接触内摩擦角、接触抗拉强度分别设为自变量，设置合理梯度并分别进行单轴压缩数值模拟和巴西劈裂数值模拟，得到接触法向刚度-宏观杨氏模量（如图3所示）、接触抗拉强度-巴西抗拉强度（如图4所示）、接触粘聚力-单轴抗压强度（如图5所示）、接触内摩擦角-单轴抗压强度的关系（如图5所示）。

将试验所得不同饱和度下的杨氏模量、巴西抗拉强度、单轴抗压强度输入至上述各类模拟参数与模拟宏观力学参数的关系，得到预测接触法向刚度、预测接触抗拉强度、预测接触粘聚力。由于接触内摩擦角-单轴抗压强度斜率绝对值较小，导致预测接触内摩擦角为负数；加之已有研究表明，接触内摩擦角主要影响内摩擦角，而内摩擦角受水致劣化程度极低，故接触内摩擦角设为恒定值。

将得到的各组预测模拟参数输入至UDEC，进行单轴压缩数值模拟和巴西劈裂数值模拟，得到各组预测模拟参数所对应的模型宏观力学参数。

对比试验和模拟所得宏观力学参数，包括杨氏模量、巴西抗拉强度、单轴抗压强度。其中，五组饱和度所对应的试验和模拟所得杨氏模量误差在0.69%以下，表明预测接触法向刚度无需调整，可确定为最终接触法向刚度（如图3中实心圆点）；模拟所得巴西抗拉强度相较试验结果偏低，经过对预测接触抗拉强度进行微调，五组饱和度所对应的试验和模拟所得巴西抗拉强度误差在0.23%以下，微调后的接触抗拉强度（如图4中实心圆点）在接触抗拉强度-巴西抗拉强度线性关系的95%预测带内，可确定为最终接触抗拉强度；模拟所得单轴抗压强度略低于试验值，对预测接触粘聚力进行微调，五组饱和度所对应的试验和模拟所得单轴抗压强度误差在0.74%以下，微调后的接触粘聚力（如图5中实心圆点）在接触粘聚力-单轴抗压强度线性关系的99%预测带内，可确定为最终接触粘聚力。

至此，五组饱和度所对应的接触法向刚度、接触粘聚力、接触内摩擦角、接触抗拉强度均已确定。将上述模拟参数输入至UDEC，进行单轴压缩数值模拟和巴西劈裂数值模拟，将模拟结果与试验结果进行对比。图6展示了模拟与试验所得单轴抗压强度、杨氏模量、巴西抗拉强度三类宏观力学参数的对比；图7中的a部分展示了模拟与试验所得轴向应力-应变曲线的对比；图7中的b部分展示了模拟与试验所得单轴压缩下破裂模式的对比；图7中的c部分展示了模拟与试验所得巴西劈裂下破裂模式的对比，可见本发明实现了水致岩石强度劣化的模拟。

基于本发明中的上述方法，本发明另外提供了一种基于离散单元法的水致岩石强度劣化的模拟系统，所述系统包括：

岩石试样矿物类型和含量信息确定模块，用于利用X射线衍射仪确定岩石试样矿物类型和含量信息。

模拟块体参数确定模块，用于根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block。

不同饱和度下的岩石试样制备模块，用于制备不同饱和度下的岩石试样；

不同饱和度岩石试样的宏观力学参数确定模块，用于利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角。

校准模块，用于选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括：接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont。

梯度设置模块，用于分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度。

第一数值模拟模块，用于分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角。

线性拟合模块，用于通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系。

预测模块，用于将不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值。

第二数值模拟模块，用于将所得不同饱和度下的各模拟参数的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，得到数值模型宏观力学参数。

比对模块，用于对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，得到对比结果。

调整模块，用于根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致。

不同饱和度下的岩石工况模拟模块，用于基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

本发明还提供一种电子设备，包括存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行上述的离散单元法的水致岩石强度劣化的模拟方法。

本发明还提供一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现上述的离散单元法的水致岩石强度劣化的模拟方法。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (9)

1.一种基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述模拟方法包括：

S1：利用X射线衍射仪确定岩石试样矿物类型和含量信息；

S2：根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block；

S3：制备不同饱和度下的岩石试样；

S4：利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

S5：选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont；

S6：分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度；

S7：分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

S8：通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系；

S9：将S4中所得不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值；

S10：将S9中所得不同饱和度下的各模拟参数的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，得到数值模型宏观力学参数；

S11：对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，得到对比结果；

S12：根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致；

S13：基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

2.根据权利要求1所述的基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数具体包括以下步骤：

S5.1：输入任意一组模拟参数，包括接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont，接触内摩擦角、接触抗拉强度σ_t ^cont；

S5.2：判断接触法向刚度和接触剪切刚度是否均小于或等于预设阈值；

S5.3：若否，则返回S5.1；

S5.4：若是，则执行下一步骤；

S5.5：判断接触剪切刚度和接触法向刚度的比值是否与试样剪切模量和试样杨氏模量的比值相等；

S5.6：若否，则返回S5.1；

S5.7：若是，则执行下一步骤；

S5.8：在离散单元法软件中进行单轴压缩数值模拟；数值模型尺寸和边界条件与单轴压缩试样和加载条件相同；

S5.9：判断单轴压缩数值模拟所得数值模型的杨氏模量是否与试验值相等；

S5.10：若否，则返回S5.1；

S5.11：若是，则确定接触剪切刚度和接触法向刚度的值，在后续步骤中即使返回S5.1仍保持接触剪切刚度和接触法向刚度值不变，并执行下一步骤；

S5.12：在离散单元法软件中进行巴西劈裂数值模拟；数值模型尺寸和边界条件与巴西劈裂试样和加载条件相同；

S5.13：判断巴西劈裂数值模拟所得的数值模型巴西抗拉强度是否与试验值相等；

S5.14：若否，则返回S5.1；

S5.15：若是，则确定接触抗拉强度的值，在后续步骤中即使返回S5.1仍保持接触抗拉强度值不变，并执行下一步骤；

S5.16：在离散单元法软件中进行单轴压缩数值模拟和三轴压缩数值模拟；数值模型尺寸和边界条件与用于单轴压缩试验和三轴压缩试验的试样和加载条件相同；

S5.17：判断单轴压缩数值模拟和三轴压缩数值模拟所得数值模型的粘聚力和内摩擦角是否与试验值相等；

S5.18：若否，则返回S5.1；

S5.19：若是，则确定接触粘聚力和接触内摩擦角，试错法结束，用于模拟选定饱和度下试样力学行为的模拟参数校准完成。

3.根据权利要求1所述的基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述制备不同饱和度下的岩石试样具体包括：

依据国际规范制备饱和岩石试样，利用烘干箱制备干燥岩石试样，在真空容器中应用质量控制法使干燥岩石试样吸水，以得到不同饱和度下的岩石试样。

4.根据权利要求1所述的基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述模拟块体表征矿物颗粒。

5.根据权利要求2所述的基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述预设阈值的表达式如下：

其中，Δz_min表示相邻单元在垂直方向的最小宽度，max表示对于与接触相邻的所有单元的最大值，K^block表示模拟块体体积模量，G^block表示模拟块体剪切模量。

6.根据权利要求1所述的基于离散单元法的水致岩石强度劣化的模拟方法，其特征在于，所述各类模拟参数与模拟所得宏观力学参数的关系具体包括：接触法向刚度-杨氏模量、接触抗拉强度-巴西抗拉强度、接触粘聚力-单轴抗压强度、接触粘聚力-粘聚力、接触内摩擦角-单轴抗压强度以及接触内摩擦角-内摩擦角的关系。

7.一种基于离散单元法的水致岩石强度劣化的模拟系统，其特征在于，所述系统包括：

岩石试样矿物类型和含量信息确定模块，用于利用X射线衍射仪确定岩石试样矿物类型和含量信息；

模拟块体参数确定模块，用于根据岩石试样矿物类型和含量信息确定模拟块体参数；所述模拟块体参数包括：模拟块体密度ρ、模拟块体体积模量K^block以及模拟块体剪切模量G^block；

不同饱和度下的岩石试样制备模块，用于制备不同饱和度下的岩石试样；

不同饱和度岩石试样的宏观力学参数确定模块，用于利用伺服刚性试验机对所述不同饱和度下的岩石试样进行单轴压缩试验、巴西劈裂试验和三轴压缩试验，得到不同饱和度岩石试样的宏观力学参数；所述宏观力学参数包括杨氏模量E、泊松比μ、剪切模量G、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

校准模块，用于选择任一饱和度下岩石试样的一组宏观力学参数，并基于所述模拟块体参数校准该饱和度下岩石试样的宏观力学参数所对应的模拟参数；所述模拟参数包括：接触法向刚度k_n、接触剪切刚度k_s、接触粘聚力c^cont、接触内摩擦角以及接触抗拉强度σ_t ^cont；

梯度设置模块，用于分别对校准后的接触法向刚度k_n、接触粘聚力c^cont、接触内摩擦角、接触抗拉强度σ_t ^cont设置梯度；

第一数值模拟模块，用于分别应用不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角、不同接触抗拉强度σ_t ^cont进行单轴压缩数值模拟、巴西劈裂数值模拟和三轴压缩数值模拟，从而得到不同接触法向刚度k_n、不同接触粘聚力c^cont、不同接触内摩擦角以及不同接触抗拉强度σ_t ^cont下数值模型的杨氏模量E、单轴抗压强度UCS、巴西抗拉强度BTS、粘聚力c以及内摩擦角；

线性拟合模块，用于通过线性拟合得到各类模拟参数与模拟所得宏观力学参数的关系；

预测模块，用于将不同饱和度下实际岩石的宏观力学参数输入至所述各类模拟参数与数值模型宏观力学参数的关系中，得到不同饱和度下各模拟参数的预测值；

第二数值模拟模块，用于将所得不同饱和度下的各模拟参数的预测值作为输入，在离散单元法软件中进行单轴压缩数值模拟、巴西劈裂数值模拟以及三轴压缩数值模拟，得到数值模型宏观力学参数；

比对模块，用于对比所述数值模型宏观力学参数与岩石试样的宏观力学参数，得到对比结果；

调整模块，用于根据对比结果，对不同饱和度下的各模拟参数的预测值进行微调，直至所得数值模型宏观力学参数与实际岩石宏观力学参数一致；

不同饱和度下的岩石工况模拟模块，用于基于所述微调后的不同饱和度下的各模拟参数的预测值对不同饱和度下的岩石工况进行模拟；所述微调后的不同饱和度下的各模拟参数的预测值为最终值。

8.一种电子设备，其特征在于，包括存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行如权利要求1-6中任一项所述的离散单元法的水致岩石强度劣化的模拟方法。

9.一种计算机可读存储介质，其特征在于，其存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1-6中任一项所述的离散单元法的水致岩石强度劣化的模拟方法。

Images