充填节理岩石损伤度离散元数值模拟方法及系统

技术领域

本发明属于充填节理岩石损伤测定技术领域，更具体地，涉及一种充填节理岩石损伤度离散元数值模拟方法及系统。

背景技术

工程岩体中随机分布有大量的节理，而节理往往含有一定的充填物。充填节理岩石不同于完整岩石和非充填节理岩石，充填节理具有低强度和大变形特性，其力学强度和变形特性更差，充填节理对整个岩石的强度和稳定性有着很大影响。实际中充填节理岩石不可避免地会遭受到爆破等冲击作用，会产生一定的损伤，而损伤程度关乎到充填节理岩石的安全稳定性，所以衡量冲击作用下充填节理岩石的损伤程度至关重要。

现有技术中往往通过波速值、抗压强度以及体积变化率等特征的变化反映损伤度，公开号为CN111208198A的中国发明专利公开了一种岩体实时波速测定及质量评价的方法，通过在不同超声波发送位置、不同时间不断按动超声波发送装置的按钮，并且不断计算每次超声波传播到每个传感器的速度，根据波速、衰减系数、波形、频率、频谱及振幅等参数变化值可反演得到岩土体的物理力学指标及细观结构特征，如节理的变化情况，进而分析矿山岩体稳定情况。公开号为CN113109182B的中国发明专利公开了一种基于体积变化率确定岩石节理动态损伤的实验方法，分别基于单次冲击和循环冲击的一维分离式霍普金森压杆SHPB实验和三维节理形貌扫描技术，可以确定岩石节理在单次冲击或循环冲击下的动态损伤状态，其与传统的基于静态或准静态力学实验的超声检测相比，更适用于岩石节理在爆破或地震作用等动荷载作用下损伤程度的测量，直接采用岩石节理本身的形貌变化反映节理的损伤程度。

但考虑到充填节理岩石可能具有复杂的内部结构，包括填充物和岩石的层状分布，可能会干扰波速值和抗压强度的变化，难以准确反映实际损伤程度；此外，充填节理岩石可能在不同部分存在非均质性，不同区域的波速值和抗压强度可能不同，难以用统一的标准评估损伤。现有研究较少从冲击作用下充填节理岩石的裂纹发育情况角度评判损伤程度，而裂纹发育情况往往决定岩体是否失稳破坏，裂纹发育情况对评价岩体稳定性具有重要参考意义。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供一种充填节理岩石损伤度离散元数值模拟方法及系统，通过构建充填节理岩石数值模型与落锤数值模型，模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，以模型中模拟形成微裂纹的发育情况这一特征评估冲击作用下充填节理岩石的损伤度，更贴合实际情况，具有精度高的优点。

为了实现上述目的，按照本发明的第一方面，提供一种充填节理岩石损伤度离散元数值模拟方法，包括：

S100：建立充填节理岩石数值模型，通过与实际物理试验结果对比进行细观参数标定，赋予充填节理岩石数值模型通过标定后的细观参数，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子a

；

；

其中，emod

S200：建立落锤数值模型，并赋予落锤数值模型的位置坐标、密度、下落初始速度和冲击接触时间参数，并模拟落锤数值模型冲击充填节理岩石数值模型，记录冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据；

S300：根据微裂纹数目以及微裂纹出现的位置数据计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，当最大贯通度大于阈值时，冲击作用下充填节理岩石损伤度为最大贯通度数值；当充填节理岩石数值模型的最大贯通度小于阈值时充填节理岩石的损伤度为：

；

式中，D表示冲击作用下充填节理岩石损伤度，f

进一步地，步骤S300中，所述阈值的确定包括：

S301：对建立考虑矿物组分、含量及矿物排列位置的充填节理岩石数值模型进行单轴压缩模拟；

S302：输出单轴压缩过程中充填节理岩石数值模型微裂纹数目随时间的变化曲线；

S303：提取变化曲线中微裂纹数目增长速度变化量最大的数据点，该数据点时刻充填节理岩石数值模型所对应的贯通度即为阈值。

进一步地，所述阈值与充填节理岩石的矿物组分、含量以及矿物排列位置相关。

进一步地，步骤S300中，所述最大贯通度为任一平面微裂纹所形成的连线与充填节理岩石数值模型在连线方向的总长度的比值。

进一步地，步骤S100中，充填节理岩石数值模型包括充填节理层和岩石层的接触数目、细观参数及厚度之类的参数。

进一步地，步骤S100中，所述充填节理岩石数值模型接触为平行黏结模型。

按照本发明的第二方面，提供一种充填节理岩石损伤度离散元数值模拟系统，包括：

损伤权重确定模块，用于建立充填节理岩石数值模型，通过与实际物理试验结果对比进行细观参数标定，赋予充填节理岩石数值模型通过标定后的细观参数，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子a

；

；

其中，emod

落锤数值模拟模块，用于建立落锤数值模型，并赋予落锤数值模型的位置坐标、密度、下落初始速度和冲击接触时间参数，并模拟落锤数值模型冲击充填节理岩石数值模型，记录冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据；

损伤度计算模块，用于根据微裂纹数目以及微裂纹出现的位置数据计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，当最大贯通度大于阈值时，冲击作用下充填节理岩石损伤度为最大贯通度数值；当充填节理岩石数值模型的最大贯通度小于阈值时充填节理岩石的损伤度为：

；

式中，D表示冲击作用下充填节理岩石损伤度，f

按照本发明的第三方面，提供一种终端装置，包括：

存储器，用于存储计算机程序；

处理器，用于执行所述计算机程序时实现充填节理岩石损伤度离散元数值模拟方法的步骤。

按照本发明的第四方面，提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现充填节理岩石损伤度离散元数值模拟方法的步骤。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1.本发明的方法，通过构建充填节理岩石数值模型与落锤数值模型，模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，以模型中模拟形成微裂纹的发育情况这一特征评估冲击作用下充填节理岩石的损伤度，更贴合实际情况，具有精度高的优点。

2.本发明通过引入损伤权重因子这一系数，基于充填节理层与岩石层细观参数弹性模量的比例关系以及厚度比，对充填节理的软弱特性进行定量分析，可以更准确地评估充填节理岩石整体的损伤程度。

3.本发明在确定冲击作用下充填节理岩石的损伤度时综合考虑微裂纹数目以及微裂纹分布特性的影响，采用对比最大贯通度与阈值的关系确定损伤度，相较于基于单一因素评估充填节理岩石损伤度的方法，所确定的损伤度具有更高的可靠性。

4.本发明通过使用颗粒流数值模拟软件，模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，可以直观呈现冲击作用下充填节理岩石的微裂纹数目和分布情况，具有良好的可视化效果。

附图说明

图1为本发明实施例一种充填节理岩石损伤度离散元数值模拟方法的流程图；

图2为本发明实施例一种充填节理岩石损伤度离散元数值模拟方法的原理图；

图3为本发明一个实施例中的落锤数值模型冲击充填节理岩石数值模型示意图；

图4为本发明一个实施例中的受冲击后充填节理岩石数值模型裂纹图；

图5为本发明一个实施例中的充填节理岩石数值模型微裂纹数目随时间的变化曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1-图5所示，本发明实施例提供一种充填节理岩石损伤度离散元数值模拟方法，该方法包括以下步骤：

S1：建立充填节理岩石数值模型，所述充填节理岩石数值模型的接触为平行黏结模型；记录充填节理层和岩石层的接触数目、细观参数及厚度，所述细观参数通过与实际物理试验结果比对进行标定；计算充填节理层和岩石层的损伤权重因子a

；

；

其中，emod

S2：建立落锤数值模型，建立所述落锤数值模型包括位置坐标、密度、下落初始速度和冲击接触时间参数的赋予，所述参数由CAD数据导入。

S3：模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度；所述最大贯通度为任一平面微裂纹所形成的连线与充填节理岩石数值模型范围内该平面总长度的比值。

S4：当充填节理岩石数值模型的最大贯通度小于阈值时，计算冲击作用下充填节理岩石的损伤度：

；

其中，D表示冲击作用下充填节理岩石损伤度，f

当最大贯通度大于阈值时，冲击作用下充填节理岩石损伤度为最大贯通度数值。

其中，所述阈值与充填节理岩石的矿物组分、含量以及矿物排列位置相关，确定所述阈值包括以下步骤：

步骤一，对建立考虑矿物组分、含量及矿物排列位置的充填节理岩石数值模型进行单轴压缩模拟；

步骤二，输出单轴压缩过程中充填节理岩石数值模型微裂纹数目随时间的变化曲线；

步骤三，提取变化曲线中微裂纹数目增长速度变化量最大的数据点，该数据点时刻充填节理岩石数值模型所对应的贯通度即为阈值。本发明在确定冲击作用下充填节理岩石的损伤度时综合考虑微裂纹数目以及微裂纹分布特性的影响，采用对比最大贯通度与阈值的关系确定损伤度，相较于基于单一因素评估充填节理岩石损伤度的方法，所确定的损伤度具有更高的可靠性。

本发明的方法，通过构建充填节理岩石数值模型与落锤数值模型，模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，以模型中模拟形成微裂纹的发育情况这一特征评估冲击作用下充填节理岩石的损伤度，更贴合实际情况，具有精度高的优点。

此外，本发明实施例还提供一种充填节理岩石损伤度离散元数值模拟系统，包括：

损伤权重确定模块，用于建立充填节理岩石数值模型，通过与实际物理试验结果对比进行细观参数标定，赋予充填节理岩石数值模型通过标定后的细观参数，根据充填节理层与岩石层细观参数的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子a

；

；

其中，emod

落锤数值模拟模块，用于建立落锤数值模型，并赋予落锤数值模型的位置坐标、密度、下落初始速度和冲击接触时间参数，并模拟落锤数值模型冲击充填节理岩石数值模型，记录冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据；

损伤度计算模块，用于根据微裂纹数目以及微裂纹出现的位置数据计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，当最大贯通度大于阈值时，冲击作用下充填节理岩石损伤度为最大贯通度数值；当充填节理岩石数值模型的最大贯通度小于阈值时充填节理岩石的损伤度为：

；

式中，D表示冲击作用下充填节理岩石损伤度，f

实施例1

在颗粒流数值模拟软件PFC中建立充填节理岩石数值模型：充填节理岩石数值模型的尺寸为50mm*35mm；赋予充填节理岩石数值模型接触为平行黏结模型并记录充填节理层和岩石层的接触数目c

表1 细观参数表

表中：pb_emod为平行粘结弹性模量、pb_krat为平行粘结刚度比、pb_ten为平行粘结抗拉强度、pb_coh为平行粘结强度、pb_fa平行粘结摩擦角、emod为弹性模量、krat为刚度比。

岩土是一种地质材料，它经受长期的地质构造作用，在一定的地质环境中形成一定的结构，显现出宽广和多变的材料响应范围。与一般工程材料相比，表现出：（1）结构上的不连续性、不均匀性和各向异性；（2）物理力学性质上的非线性；（3）随机性、模糊性和不确定性；（4）地质条件和工程条件的复杂性；（5）整个岩土工程系统的非线性。

基于岩石结构与力学性质上的不连续性和非线性，采用颗粒流数值模拟（particle follow code，PFC）进行建模，模拟岩石这种非连续介质材料的力学行为。颗粒流数值模拟是一种离散单元法，通过圆形颗粒介质的运动及其相互作用来模拟颗粒材料的力学特性。在这种颗粒单元研究的基础上，通过非连续的数值方法来解决包含复杂变形模式的实际问题。在具有颗粒结构特性的充填节理岩石中的应用，就是从其细观力学特征出发，将材料的力学响应问题从物理域映射到数学域内进行数值求解。与此相应，物理域内实物颗粒被抽象为数学域内的颗粒单元，并通过颗粒单元来构建和设计任意几何性状的试样，颗粒间的相互作用通过接触本构关系来实现，数值边界条件的确定和试样的若干应力平衡状态通过迭代分析进行，直到使数值试样的宏观力学特性逼近材料的真实力学行为或者工程特性。

PFC程序中内置接触模型有多种：线性模型、接触黏结模型、平行黏结模型、平直节理模型和光滑节理模型等，其中平行黏结模型可以很好地模拟岩石材料的力学行为。本发明通过使用颗粒流数值模拟软件，模拟落锤数值模型冲击充填节理岩石数值模型，记录并输出冲击作用下充填节理层和岩石层的微裂纹数目以及微裂纹出现的位置数据，可以直观呈现冲击作用下充填节理岩石的微裂纹数目和分布情况，具有良好的可视化效果。

本实施例中，根据充填节理层与岩石层细观参数弹性模量（emod）的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子a

；

其中：emod

本实施例中，构建落锤数值模型，落锤数值模型特征由CAD数据导入，落锤数值模型如图3所示，赋予落锤数值模型的位置坐标为（0，0.2m）、密度为2.7×10

本实施例中，模拟落锤数值模型冲击充填节理岩石数值模型，编制命令流分别记录并输出落锤数值模型冲击充填节理岩石数值模型下充填节理层和岩石层的微裂纹数目f

本实施例中，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，最大贯通度为任一平面微裂纹所形成的连线与充填节理岩石数值模型范围内该平面总长度的比值，其值为0.387；

本实施例中，计算冲击作用下充填节理岩石数值模型的阈值；

所述阈值与充填节理岩石的矿物组分、含量及矿物排列位置相关，所述阈值确定的步骤为：

步骤一、对建立考虑矿物组分、含量及矿物排列位置的充填节理岩石数值模型进行单轴压缩模拟；

步骤二、输出单轴压缩过程中充填节理岩石数值模型微裂纹数目随时间的变化曲线，如图5所示；

步骤三、提取变化曲线中微裂纹数目增长速度变化量最大的数据点，所选数据点如图5所示，该数据点时刻充填节理岩石数值模型所对应的贯通度即为阈值，该阈值为0.7。

本实施例中，所述最大贯通度小于阈值，冲击作用下充填节理岩石损伤度计算公式为：

故本实施中例冲击作用下充填节理岩石的损伤度为0.425。

实施例2

在颗粒流数值模拟软件PFC中建立充填节理岩石数值模型：充填节理岩石数值模型的尺寸为60mm*40mm；赋予充填节理岩石数值模型接触为平行黏结模型并记录充填节理层和岩石层的接触数目c

本实施例中，根据充填节理层与岩石层细观参数弹性模量（emod）的比例关系和充填节理层与岩石层的厚度比确定充填节理层和岩石层的损伤权重因子a

；

其中：emod

本实施例中，构建落锤数值模型，落锤数值模型特征由CAD数据导入，赋予落锤数值模型的位置坐标为（0，0.25m）、密度为3.2×10

本实施例中，模拟落锤数值模型冲击充填节理岩石数值模型，编制命令流分别记录并输出落锤数值模型冲击充填节理岩石数值模型下充填节理层和岩石层的微裂纹数目f

本实施例中，计算冲击作用下产生的微裂纹在充填节理岩石数值模型任一平面的最大贯通度，最大贯通度为任一平面微裂纹所形成的连线与充填节理岩石数值模型范围内该平面总长度的比值，其值为0.633；

本实施例中，计算冲击作用下充填节理岩石数值模型的阈值；

所述阈值与充填节理岩石的矿物组分、含量及矿物排列位置相关，所述阈值确定的步骤为：

步骤一、对建立考虑矿物组分、含量及矿物排列位置的充填节理岩石数值模型进行单轴压缩模拟；

步骤二、输出单轴压缩过程中充填节理岩石数值模型微裂纹数目随时间的变化曲线；

步骤三、提取变化曲线中微裂纹数目增长速度变化量最大的数据点，该数据点时刻充填节理岩石数值模型所对应的贯通度即为阈值，该阈值为0.6。

本实施例中，所述最大贯通度大于阈值，冲击作用下充填节理岩石损伤度为最大贯通度，其数值为0.633。

故本实施中例冲击作用下充填节理岩石的损伤度为0.633。

实施例3

本申请还提供了一种终端装置，该终端装置可包括：

存储器，用于存储计算机程序；

处理器，用于执行计算机程序时可实现如上述任意一种确定冲击作用下充填节理岩石损伤度方法的步骤。

终端装置可以包括：处理器、存储器、通信接口和通信总线。处理器、存储器、通信接口均通过通信总线完成相互间的通信。

在本申请实施例中，处理器可以为中央处理器（Central Processing Unit，CPU）、特定应用集成电路、数字信号处理器、现场可编程门阵列或者其他可编程逻辑器件等。

处理器可以调用存储器中存储的程序，具体的，处理器可以执行异常IP识别方法的实施例中的操作。

存储器中用于存放一个或者一个以上程序，程序可以包括程序代码，程序代码包括计算机操作指令，在本申请实施例中，存储器中至少存储有用于实现确定冲击作用下充填节理岩石损伤度方法功能的程序。

实施例4

本发明的一个实施例还公开了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如本发明第一方面所述的充填节理岩石损伤度离散元数值模拟方法的步骤。

本发明实现上述实施例方法中的全部或部分步骤，可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储的介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器（ROM，Read-OnlyMemory）、随机存取存储器（RAM，Random Access Memory）、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括是电载波信号和电信信号。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。