一种模拟巷道围岩蠕变失稳的离散元方法及系统

1.本发明涉及围岩支护技术领域，特别是涉及一种模拟巷道围岩蠕变失稳的离散元方法及系统。


背景技术：

2.隧道岩体破坏的外在表现形式常表现为开挖面岩体的开裂、剥落、屈曲、滑移及岩爆等，其内在原因则是岩体内部原生裂纹或次生裂纹在应力作用下的萌生、扩展、贯通过程，从而在岩体一定范围内形成开挖损伤区或开挖扰动区。可见，准确描述和预测在围岩开挖扰动作用下围岩的损伤情况，对于确保地下工程进行时围岩的安全与稳定有着重要意义。
3.目前，对岩体开挖引起的损伤和破坏的预测与描述常采用有限差分和离散元数值模拟方法，利用mohr-coulomb强度准则进行预测与计算围岩的损伤区范围。现有技术存在利用有限差分法建立分析模型，模拟洞室开挖诱导的能量释放过程和计算能量释放系数的大小，确定开挖过程中围岩损伤与围岩能量释放之间的关系，从而预测出围岩损伤范围。这种基于能量释放系数的围岩损伤预测方法适用于研究巷道围岩的宏观破坏，但不能揭示岩体裂隙拓展对巷道围岩破坏的本质。现有技术还存在利用udec离散元数值模拟软件从裂隙拓展的角度揭示了巷道围岩破坏的本质，但是，在udec离散元软件中有viscous模型、burgers模型和cvisc模型等8种蠕变模型用于岩石的蠕变特性分析，但以上蠕变模型均未考虑岩石的变形和应力随着时间因素的影响，无法实现岩石加速蠕变阶段和蠕变过程中的微裂隙扩展，使得模拟巷道围岩蠕变失稳的研究结果不准确。


技术实现要素：

4.本发明的目的是提供一种模拟巷道围岩蠕变失稳的离散元方法及系统，可以增加模拟巷道围岩蠕变失稳研究结果的准确度。
5.为实现上述目的，本发明提供了如下方案：
6.一种模拟巷道围岩蠕变失稳的离散元方法，包括：
7.获取待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线；所述赋存状况包括：巷道尺寸、埋深、巷道周围采掘情况和设定范围内的岩体参数；所述岩体参数包括：岩体倾角、岩体厚度、岩体岩性和岩体的层理结构；所述物理力学参数包括：单轴压力试验应力
–
应变曲线、巴西劈裂实验应力
–
应变曲线和岩石蠕变曲线；
8.利用udec离散元模拟软件建立压缩劈裂模型并根据所述单轴压力试验应力
–
应变曲线和所述巴西劈裂实验应力
–
应变曲线对所述压缩劈裂模型的接触参数和块体参数进行反演得到合理接触面参数值和合理块体参数值；所述压缩劈裂模型包括单轴压缩数值模型和巴西劈裂数值模型；
9.根据所述合理块体参数值和所述合理接触面参数值利用udec离散元模拟软件建
立块体单元为cvisc蠕变模型的第一单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演得到合理块体蠕变参数值；
10.根据所述合理接触面参数值和所述合理块体蠕变参数值利用udec离散元模拟软件构建块体单元为改进cvisc蠕变模型的第二单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的块体蠕变参数退化方程进行反演得到合理块体蠕变参数退化方程；所述第二单轴蠕变数值模型的块体单元为改进cvisc蠕变模型；所述改进cvisc蠕变模型为退化方程为块体蠕变参数退化方程的cvisc蠕变模型；
11.根据所述待模拟巷道围岩的变形量随时间变化的曲线、所述赋存状况和所述合理块体蠕变参数退化方程构建巷道离散元数值仿真模型并采用离散法对所述巷道离散元数值仿真模型进行网格划分得到经离散网格化处理的巷道离散元数值仿真模型；
12.根据所述待模拟巷道围岩的赋存状况、所述待模拟巷道设定范围内岩层的物理力学参数和所述待模拟巷道围岩的变形量随时间变化的曲线对所述经离散网格化处理的巷道离散元数值仿真模型进行模拟开挖得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；所述累计张拉裂隙长度为所述网格内所有的张拉裂隙长度和；所述累计剪切裂隙长度为所述网格内所有的剪切裂隙长度和；
13.根据各所述网格的累计张拉裂隙长度、累计剪切裂隙长度和各所述网格的临界蠕变损伤值确定所述待模拟巷道围岩蠕变失稳的范围。
14.一种模拟巷道围岩蠕变失稳的离散元系统，包括：
15.获取模块，用于获取待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线；所述赋存状况包括：巷道尺寸、埋深、巷道周围采掘情况和设定范围内的岩体参数；所述岩体参数包括：岩体倾角、岩体厚度、岩体岩性和岩体的层理结构；所述物理力学参数包括：单轴压力试验应力
–
应变曲线、巴西劈裂实验应力
–
应变曲线和岩石蠕变曲线；
16.合理接触面参数值和合理块体参数值确定模块，用于利用udec离散元模拟软件建立压缩劈裂模型并根据所述单轴压力试验应力
–
应变曲线和所述巴西劈裂实验应力
–
应变曲线对所述压缩劈裂模型的接触参数和块体参数进行反演得到合理接触面参数值和合理块体参数值；所述压缩劈裂模型包括单轴压缩数值模型和巴西劈裂数值模型；
17.合理块体蠕变参数值确定模块，用于根据所述合理块体参数值和所述合理接触面参数值利用udec离散元模拟软件建立块体单元为cvisc蠕变模型的第一单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演得到合理块体蠕变参数值；
18.合理块体蠕变参数退化方程确定模块，用于根据所述合理接触面参数值和所述合理块体蠕变参数值利用udec离散元模拟软件构建块体单元为改进cvisc蠕变模型的第二单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的块体蠕变参数退化方程进行反演得到合理块体蠕变参数退化方程；所述第二单轴蠕变数值模型的块体单元为改进cvisc蠕变模型；所述改进cvisc蠕变模型为退化方程为块体蠕变参数退化方程的cvisc蠕变模型；
19.巷道离散元数值仿真模型确定模块，用于根据所述待模拟巷道围岩的变形量随时间变化的曲线、所述赋存状况和所述合理块体蠕变参数退化方程构建巷道离散元数值仿真
模型并采用离散法对所述巷道离散元数值仿真模型进行网格划分得到经离散网格化处理的巷道离散元数值仿真模型；
20.裂隙长度确定模块，用于根据所述待模拟巷道围岩的赋存状况、所述待模拟巷道设定范围内岩层的物理力学参数和所述待模拟巷道围岩的变形量随时间变化的曲线对所述经离散网格化处理的巷道离散元数值仿真模型进行模拟开挖得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；所述累计张拉裂隙长度为所述网格内所有张拉裂隙长度的和；所述累计剪切裂隙长度为所述网格内所有剪切裂隙长度的和；
21.失稳范围确定模块，用于根据各所述网格的累计张拉裂隙长度、累计剪切裂隙长度和各所述网格的临界蠕变损伤值确定所述待模拟巷道围岩蠕变失稳的范围。
22.根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明获取待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线；建立压缩劈裂模型并根据单轴压力试验应力
–
应变曲线和巴西劈裂实验应力
–
应变曲线对压缩劈裂模型的接触参数和块体参数进行反演得到合理接触面参数值和合理块体参数值；根据合理块体参数值和合理接触面参数值建立块体单元为cvisc蠕变模型的第一单轴蠕变数值模型并根据岩石蠕变曲线对第一单轴蠕变数值模型的块体蠕变参数进行反演得到合理块体蠕变参数值；根据合理接触面参数值和合理块体蠕变参数值构建块体单元为改进cvisc蠕变模型的第二单轴蠕变数值模型并根据岩石蠕变曲线对第二单轴蠕变数值模型的块体蠕变参数退化方程进行反演得到合理块体蠕变参数退化方程；根据待模拟巷道围岩的变形量随时间变化的曲线、赋存状况和合理块体蠕变参数退化方程构建巷道离散元数值仿真模型并采用离散法对巷道离散元数值仿真模型进行网格划分得到经离散网格化处理的巷道离散元数值仿真模型；根据待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线对经离散网格化处理的巷道离散元数值仿真模型进行模拟开挖得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；根据各网格的累计张拉裂隙长度、累计剪切裂隙长度和各网格的临界蠕变损伤值确定待模拟巷道围岩蠕变失稳的范围，本发明的蠕变参数退化方程可以实现加速蠕变阶段和蠕变过程中的微裂隙扩展，使得模拟巷道围岩蠕变失稳的研究结果更加准确。
附图说明
23.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
24.图1为本发明实施例提供的模拟巷道围岩蠕变失稳的离散元方法的流程图；
25.图2为本发明实施例步骤s1、s2中岩石单轴压缩试验中应力
–
应变曲线图；
26.图3为本发明实施例步骤s1、s2中岩石巴西劈裂试验中应力
–
应变曲线图；
27.图4为本发明实施例步骤s1、s2中岩石蠕变曲线图；
28.图5为本发明实施例步骤s3中岩石单轴蠕变退化试验时间
–
应变曲线图；
29.图6为本发明实施例步骤s4中工程离散元数值仿真模型图；
30.图7为本发明实施例步骤s4中巷道围岩离散化网格示意图；
31.图8为本发明实施例步骤s6中不同时间巷道围岩蠕变损伤演化曲线图；
32.图9为本发明实施例步骤s6中围岩损伤破坏范围图；
33.图10为本发明实施例提供的模拟巷道围岩蠕变失稳的离散元方法的大体流程图。
具体实施方式
34.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
35.为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
36.本发明公布了一种模拟巷道围岩蠕变失稳的离散元方法，如图1所示，具体步骤如下：
37.步骤101：获取待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线。所述赋存状况包括：巷道尺寸、埋深、巷道周围采掘情况和设定范围内的岩体参数；所述岩体参数包括：岩体倾角、岩体厚度、岩体岩性和岩体的层理结构；所述物理力学参数包括：单轴压力试验应力
–
应变曲线、巴西劈裂实验应力
–
应变曲线和岩石蠕变曲线，其中设定范围可以为巷道三倍宽度范围。
38.步骤102：利用udec离散元模拟软件建立压缩劈裂模型并根据所述单轴压力试验应力
–
应变曲线和所述巴西劈裂实验应力
–
应变曲线对所述压缩劈裂模型的接触参数和块体参数进行反演得到合理接触面参数值和合理块体参数值；所述压缩劈裂模型包括单轴压缩数值模型和巴西劈裂数值模型。
39.步骤103：根据所述合理块体参数值和所述合理接触面参数值利用udec离散元模拟软件建立块体单元为cvisc蠕变模型的第一单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演得到合理块体蠕变参数值。
40.步骤104：根据所述合理接触面参数值和所述合理块体蠕变参数值利用udec离散元模拟软件构建块体单元为改进cvisc蠕变模型的第二单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的块体蠕变参数退化方程进行反演得到合理块体蠕变参数退化方程；所述第二单轴蠕变数值模型的块体单元为改进cvisc蠕变模型；所述改进cvisc蠕变模型为退化方程为块体蠕变参数退化方程的cvisc蠕变模型。
41.步骤105：根据所述待模拟巷道围岩的变形量随时间变化的曲线、所述赋存状况和所述合理块体蠕变参数退化方程构建巷道离散元数值仿真模型并采用离散法对所述巷道离散元数值仿真模型进行网格划分得到经离散网格化处理的巷道离散元数值仿真模型。
42.步骤106：根据所述待模拟巷道围岩的赋存状况、所述待模拟巷道设定范围内岩层的物理力学参数和所述待模拟巷道围岩的变形量随时间变化的曲线对所述经离散网格化处理的巷道离散元数值仿真模型进行模拟开挖得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；所述累计张拉裂隙长度为所述网格内所有的张拉裂隙长度和；所述累计剪切裂隙长度为所述网格内所有的剪切裂隙长度和。
43.步骤107：根据各所述网格的累计张拉裂隙长度、累计剪切裂隙长度和各所述网格的临界蠕变损伤值确定所述待模拟巷道围岩蠕变失稳的范围。
44.在实际应用中，步骤101具体包括：
45.步骤1.1：收集待模拟巷道围岩赋存状况包括巷道尺寸、埋深、巷道周围采掘情况和巷道三倍宽度范围内的岩体倾角、厚度、岩性、层理结构。
46.步骤1.2：室内实验测量获得待模拟巷道三倍宽度范围内岩层的物理力学参数包括：岩石应力-应变曲线(单轴压力试验应力
–
应变曲线和巴西劈裂实验应力
–
应变曲线)、岩石蠕变曲线。
47.步骤1.3：观测待模拟巷道围岩变形量随时间的演化规律得到待模拟巷道围岩的变形量随时间变化的曲线。
48.在实际应用中，步骤102具体包括：
49.利用udec离散元模拟软件构建单轴压缩数值模型和巴西劈裂数值模型。使用迭代试错法根据所述单轴压力试验应力
–
应变曲线对所述单轴压缩数值模型的块体参数和接触面参数进行反演以使所述单轴压缩数值模型与所述单轴压力试验应力
–
应变曲线相匹配，根据所述巴西劈裂实验应力
–
应变曲线对所述巴西劈裂数值模型的块体参数和接触面参数进行反演以使所述巴西劈裂数值模型与所述巴西劈裂实验应力
–
应变曲线相匹配，得到合理接触面参数值和合理块体参数值(具体为：步骤2.1：利用udec离散元模拟软件建立直径为50mm高为100mm的单轴压缩数值模型和直径为50mm的巴西劈裂数值模型，模型块体采用三角形块体或者voronoi多边形块体，块体单元和接触面均采用mohr-coulomb模型，块体参数包括体积模量k、剪切模量g、内摩擦角θz、内聚力cz和抗拉强度tz，接触面参数包括法相刚度kn和剪切刚度ks、接触面内摩擦角θj、接触面内聚力cj和接触面抗拉强度tj，采用迭代试错的方法，反演块体参数和接触面参数，使得单轴压缩数值模型和巴西劈裂数值模型与步骤1.2获得的岩石应力-应变曲线(单轴压力试验应力
–
应变曲线和巴西劈裂实验应力
–
应变曲线)相匹配，从而获得合理模型块体参数和合理接触面参数)。
50.在实际应用中，步骤103具体包括：
51.利用udec离散元模拟软件建立块体单元为cvisc蠕变模型、块体参数为所述合理块体参数值，接触面参数为所述合理接触面参数值的第一单轴蠕变数值模型。使用迭代试错法根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演以使所述第一单轴蠕变数值模型与所述岩石蠕变曲线相匹配得到合理块体蠕变参数值(具体为：步骤2.2：建立直径为50mm高为100mm的第一单轴蠕变数值模型，模型块体形状和大小与步骤2.1相同，块体单元采用软件自带的cvisc蠕变模型，块体蠕变参数包括maxwell剪切模量gm、maxwell粘度系数ηm，kelvin剪切模量gk、kelvin粘度系数ηk、内摩擦角θz、内聚力cz和抗拉强度tz，接触面参数与步骤2.1相同，采用迭代试错的方法，反演块体蠕变参数，使得第一单轴蠕变数值模型的初始蠕变和等速蠕变阶段与步骤1.2获得的岩石蠕变曲线相匹配，获得合理块体蠕变参数值)。
52.采用软件自带的cvisc蠕变模型可以较好的模拟岩石初始蠕变阶段和等速蠕变阶段，但由于cvisc蠕变模型未考虑岩石在恒定载荷作用下其力学参数随时间不断劣化的过程，无法实现岩石蠕变过程中的加速蠕变阶段和加速蠕变阶段的微裂隙快速扩展。因此，需要在cvisc蠕变模型基础上自定义块体蠕变参数退化方程，所以在实际应用中，步骤104具
体包括：
53.利用udec离散元模拟软件建立块体单元为改进cvisc蠕变模型、接触面参数为所述合理接触面参数值，块体蠕变参数为所述合理块体蠕变参数值的第二单轴蠕变数值模型。使用迭代试错法根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的改进cvisc蠕变模型中的块体蠕变参数退化方程中的蠕变常数进行反演以使所述第二单轴蠕变数值模型与所述岩石蠕变曲线相匹配得到合理块体蠕变参数退化方程(具体为，步骤3：建立直径为50mm高为100mm的第二单轴蠕变数值模型，将步骤2.1标定的合理接触面参数值和步骤2.2标定的合理块体蠕变参数值进行赋参，将自定义的初始蠕变参数退化方程嵌入到主程序中，采用迭代试错的方法，反演块体蠕变参数退化方程，使得第二单轴蠕变数值模型的初始蠕变阶段、等速蠕变阶段和加速蠕变阶段与步骤1.2获得的岩石蠕变曲线相匹配，实现蠕变过程中的微裂隙快速扩展，从而获得合理蠕变参数退化方程)。
54.在实际应用中，初始蠕变参数退化方程为：
[0055][0056]
式中t为试件的寿命；t0为短期强度时间，由荷载为单轴抗强度σc的单轴蠕变破坏时间确定；cz(t)为蠕变时间t时刻的块体内聚力；σ
cd
为蠕变损伤阈值，约定于1/2倍的单轴抗强度；σ1为块体实际受到的最大主应力；k为蠕变常数，采用迭代试错的方法确定。
[0057]
在实际应用中，步骤105具体包括：
[0058]
利用udec离散元模拟软件根据所述赋存状况构建接触面参数为所述合理接触面参数值、块体参数为所述合理块体参数值、块体蠕变参数为所述合理块体蠕变参数值、退化方程为所述合理块体蠕变参数退化方程的初始仿真模型。对所述初始仿真模型进行模拟开挖并采用迭代试错法对所述初始仿真模型的接触面参数、块体参数和块体蠕变参数进行反演以使所述初始仿真模型与所述待模拟巷道围岩的变形量随时间变化的曲线相匹配得到巷道离散元数值仿真模型(具体为：步骤4.1：根据步骤1.1获取的巷道围岩赋存状况包括巷道埋深、巷道周围采掘情况和巷道三倍宽度范围内的岩体倾角、厚度、岩性、层理结构，利用udec软件建立符合巷道围岩赋存状况的巷道数值仿真几何模型，巷道三倍宽度范围内围岩采用三角形块体或者voronoi多边形块体，块体边长不超过1/10倍巷道宽度，施加边界条件，将步骤2.1确定的合理接触面参数值、合理块体参数值和步骤2.2确定的合理块体蠕变参数值对巷道数值仿真几何模型进行初次赋参，进行模型的初始平衡。步骤4.2：将合理块体蠕变参数退化方程嵌入到巷道数值仿真几何模型主程序中，对步骤4.1中初始平衡后的模型，进行模拟开挖，并采用试错迭代的方法，调整蠕变参数(块体蠕变参数、模型块体参数和接触面参数)，使得巷道数值仿真几何模型与步骤1.3的巷道围岩变形量随时间的演化规律相匹配，得到与工程实际相匹配的蠕变参数。步骤4.3：将步骤4.2调整后的蠕变参数对步骤4.1的模型进行重新赋参，得到重新赋值的巷道数值仿真模型。将巷道离散元数值仿真模型离散为规则网格(网格离散化处理是为了更好地分析围岩的破坏，将研究区域用人为抽象出来的规则网格进行划分，利用软件自带的fish语言进行编程可对网格范围内的岩体裂隙发育情况进行实时监测)，并对各网格进行唯一识别标记，将唯一识别数据(为了方便区分不同区域，将网格进行编号，每个网格的编号不重复，编号唯一)标记在各网格处备用)。
[0059]
在实际应用中，所述临界蠕变损伤值的确定过程具体包括：
[0060]
根据所述巷道离散元数值仿真模型的接触面参数值、块体参数值、块体蠕变参数值和块体蠕变参数退化方程构建第三单轴蠕变数值模型。采用分级加载方式，根据所述第三单轴蠕变数值模型进行单轴蠕变数值模拟，并以单轴蠕变数值模拟过程中试件发生加速蠕变时的临界蠕变损伤值作为所述临界蠕变损伤值(具体为：步骤5.1：建立直径为2m高为4m的第三单轴蠕变数值模型，模型块体形状和大小与步骤4.1相同，蠕变参数采用步骤4.2确定的参数，通过fish语言将合理蠕变参数退化方程嵌入到数值模型主程序中。步骤5.2：采用分级加载方式，进行单轴蠕变数值模拟，直至试件发生加速蠕变破坏，通过fish语言统计单轴蠕变模拟过程中的裂隙演化情况包括张拉裂隙数量、剪切裂隙数量、张拉裂隙长度和剪切裂隙长度随时间的变化情况；步骤5.3：以单轴蠕变数值模拟过程中试件发生加速蠕变时的临界蠕变损伤d
min
作为岩体蠕变失稳判据)。
[0061]
在实际应用中，步骤106具体包括：
[0062]
步骤6.1：根据步骤1.1收集的矿井实际情况，对步骤4.4经离散网格化处理的巷道离散元数值仿真模型，进行模拟开挖，通过fish语言统计巷道蠕变条件下各网格内张拉裂隙和剪切裂隙时空演化规律得到张拉裂隙长度和剪切裂隙长度。
[0063]
在实际应用中，步骤107具体包括：
[0064]
根据各所述网格的累计张拉裂隙长度和累计剪切裂隙长度确定各所述网格的岩石蠕变损伤。
[0065]
根据各所述网格的岩石蠕变损伤和各所述网格的临界蠕变损伤值d
min
确定所述待模拟巷道围岩蠕变失稳的范围。
[0066]
在实际应用中，根据各所述网格的累计张拉裂隙长度和累计剪切裂隙长度确定各所述网格的岩石蠕变损伤具体为：步骤6.2：将步骤6.1统计的蠕变张拉裂隙和剪切裂隙情况代入到公式计算所述网格的岩石蠕变损伤中，计算得到各网格内岩石蠕变损伤di，将岩石蠕变损伤di与d
min
进行比对，若di》d
min
，则该网格内岩体发生蠕变破坏，否则稳定，从而确定巷道围岩失稳时间和破坏范围，其中，di为岩石蠕变损伤，ls为累计张拉裂隙长度，ln为累计剪切裂隙长度，l为理论上的累计裂隙长度，由fish语言统计得到。
[0067]
本发明实施例还提供了一种模拟巷道围岩蠕变失稳的离散元系统，包括：
[0068]
获取模块，用于获取待模拟巷道围岩的赋存状况、待模拟巷道设定范围内岩层的物理力学参数和待模拟巷道围岩的变形量随时间变化的曲线。所述赋存状况包括：巷道尺寸、埋深、巷道周围采掘情况和设定范围内的岩体参数；所述岩体参数包括：岩体倾角、岩体厚度、岩体岩性和岩体的层理结构；所述物理力学参数包括：单轴压力试验应力
–
应变曲线、巴西劈裂实验应力
–
应变曲线和岩石蠕变曲线。
[0069]
合理接触面参数值和合理块体参数值确定模块，用于利用udec离散元模拟软件建立压缩劈裂模型并根据所述单轴压力试验应力
–
应变曲线和所述巴西劈裂实验应力
–
应变曲线对所述压缩劈裂模型的接触参数和块体参数进行反演得到合理接触面参数值和合理块体参数值；所述压缩劈裂模型包括单轴压缩数值模型和巴西劈裂数值模型。
[0070]
合理块体蠕变参数值确定模块，用于根据所述合理块体参数值和所述合理接触面参数值利用udec离散元模拟软件建立块体单元为cvisc蠕变模型的第一单轴蠕变数值模型
并根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演得到合理块体蠕变参数值。
[0071]
合理块体蠕变参数退化方程确定模块，用于根据所述合理接触面参数值和所述合理块体蠕变参数值利用udec离散元模拟软件构建块体单元为改进cvisc蠕变模型的第二单轴蠕变数值模型并根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的块体蠕变参数退化方程进行反演得到合理块体蠕变参数退化方程；所述第二单轴蠕变数值模型的块体单元为改进cvisc蠕变模型；所述改进cvisc蠕变模型为退化方程为块体蠕变参数退化方程的cvisc蠕变模型。
[0072]
巷道离散元数值仿真模型确定模块，用于根据所述待模拟巷道围岩的变形量随时间变化的曲线、所述赋存状况和所述合理块体蠕变参数退化方程构建巷道离散元数值仿真模型并采用离散法对所述巷道离散元数值仿真模型进行网格划分得到经离散网格化处理的巷道离散元数值仿真模型。
[0073]
裂隙长度确定模块，用于根据所述待模拟巷道围岩的赋存状况、所述待模拟巷道设定范围内岩层的物理力学参数和所述待模拟巷道围岩的变形量随时间变化的曲线对所述经离散网格化处理的巷道离散元数值仿真模型进行模拟开挖得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；所述累计张拉裂隙长度为所述网格内所有张拉裂隙长度的和；所述累计剪切裂隙长度为所述网格内所有剪切裂隙长度的和。
[0074]
失稳范围确定模块，用于根据各所述网格的累计张拉裂隙长度、累计剪切裂隙长度和各所述网格的临界蠕变损伤值确定所述待模拟巷道围岩蠕变失稳的范围。
[0075]
作为一种可选的实施方式，所述合理接触面参数值和合理块体参数值确定模块，具体包括：
[0076]
单轴压缩数值模型和巴西劈裂数值模型构建单元，用于利用udec离散元模拟软件构建单轴压缩数值模型和巴西劈裂数值模型。
[0077]
合理接触面参数值和合理块体参数值确定单元，用于使用迭代试错法根据所述单轴压力试验应力
–
应变曲线对所述单轴压缩数值模型的块体参数和接触面参数进行反演以使所述单轴压缩数值模型与所述单轴压力试验应力
–
应变曲线相匹配，根据所述巴西劈裂实验应力
–
应变曲线对所述巴西劈裂数值模型的块体参数和接触面参数进行反演以使所述巴西劈裂数值模型与所述巴西劈裂实验应力
–
应变曲线相匹配，得到合理接触面参数值和合理块体参数值。
[0078]
作为一种可选的实施方式，所述合理块体蠕变参数值确定模块，具体包括：
[0079]
第一单轴蠕变数值模型构建单元，用于利用udec离散元模拟软件建立块体单元为cvisc蠕变模型、块体参数为所述合理块体参数值，接触面参数为所述合理接触面参数值的第一单轴蠕变数值模型。
[0080]
合理块体蠕变参数值确定单元，用于使用迭代试错法根据所述岩石蠕变曲线对所述第一单轴蠕变数值模型的块体蠕变参数进行反演以使所述第一单轴蠕变数值模型与所述岩石蠕变曲线相匹配得到合理块体蠕变参数值。
[0081]
作为一种可选的实施方式，所述合理块体蠕变参数退化方程确定模块，具体包括：
[0082]
第二单轴蠕变数值模型构建单元，用于利用udec离散元模拟软件建立块体单元为改进cvisc蠕变模型、接触面参数为所述合理接触面参数值，块体蠕变参数为所述合理块体
蠕变参数值的第二单轴蠕变数值模型。
[0083]
合理块体蠕变参数退化方程确定单元，用于使用迭代试错法根据所述岩石蠕变曲线对所述第二单轴蠕变数值模型的改进cvisc蠕变模型中的块体蠕变参数退化方程中的蠕变常数进行反演以使所述第二单轴蠕变数值模型与所述岩石蠕变曲线相匹配得到合理块体蠕变参数退化方程。
[0084]
作为一种可选的实施方式，所述巷道离散元数值仿真模型确定模块，具体包括：
[0085]
初始仿真模型构建单元，用于利用udec离散元模拟软件根据所述赋存状况构建接触面参数为所述合理接触面参数值、块体参数为所述合理块体参数值、块体蠕变参数为所述合理块体蠕变参数值、退化方程为所述合理块体蠕变参数退化方程的初始仿真模型。
[0086]
巷道离散元数值仿真模型确定单元，用于对所述初始仿真模型进行模拟开挖并采用迭代试错法对所述初始仿真模型的接触面参数、块体参数和块体蠕变参数进行反演以使所述初始仿真模型与所述待模拟巷道围岩的变形量随时间变化的曲线相匹配得到巷道离散元数值仿真模型。
[0087]
作为一种可选的实施方式，所述失稳范围确定模块，具体包括：
[0088]
岩石蠕变损伤确定单元，用于根据各所述网格的累计张拉裂隙长度和累计剪切裂隙长度确定各所述网格的岩石蠕变损伤。
[0089]
失稳范围确定单元，用于根据各所述网格的岩石蠕变损伤和各所述网格的临界蠕变损伤值确定所述待模拟巷道围岩蠕变失稳的范围。
[0090]
作为一种可选的实施方式，模拟巷道围岩蠕变失稳的离散元系统还包括：
[0091]
第三单轴蠕变数值模型构建模块，用于根据所述巷道离散元数值仿真模型的接触面参数值、块体参数值、块体蠕变参数值和块体蠕变参数退化方程构建第三单轴蠕变数值模型。
[0092]
临界蠕变损伤值确定模块，用于采用分级加载方式，根据所述第三单轴蠕变数值模型进行单轴蠕变数值模拟，并以单轴蠕变数值模拟过程中试件发生加速蠕变时的临界蠕变损伤值作为所述临界蠕变损伤值。
[0093]
本发明实施例还提供了一种采用上述方法进行处理的更具体的过程，如图10所示，包括巷道围岩地质力学测试、标定岩石蠕变参数、反演块体蠕变退化方程、构建巷道离散元数值仿真模型、建立岩石蠕变失稳判据、巷道围岩蠕变失稳预测六个步骤，具体过程如下：
[0094]
工程背景：某矿主采3#煤层，33号工作面位于井下3#煤层。3#煤层平均埋深为350m，厚度4.50m～6.30m，平均5.90m；煤层倾角1
°
～9
°
，平均5
°
，为近水平煤层。直接顶厚度1.4～3.0m，平均2.0m，岩性多为灰白色粉砂岩。基本顶为层理较为发育的砂质泥岩，厚度2～10m，平均5.0m。直接底为平均厚度4.0m的砂质泥岩。33号进风顺槽服务于本工作面的进风任务，同时兼做下一工作面的辅助回风顺槽，属二次复用巷道。
[0095]
步骤s1：巷道围岩地质力学测试。
[0096]
步骤s1.1：以33进风顺槽为研究对象，巷道净高2.6m，净宽4.2m，沿煤层底板掘进，煤层平均埋深为350m、巷道位于东北部为采区大巷，东南部为3201工作面，西南部为井田边界，西北部为未采区，巷道三倍宽度范围内岩层赋存情况如表1所示。
[0097]
表1巷道三倍宽度范围内岩层赋存情况表
[0098][0099]
步骤s1.2：室内实验测量获得巷道三倍宽度范围内岩层的物理力学参数包括：岩石单轴压缩试验中应力
–
应变曲线图如图2所示、岩石巴西劈裂试验中应力
–
应变曲线图如图3所示、岩石蠕变曲线(岩石单轴蠕变试验时间
–
应变曲线)如图4。
[0100]
步骤s1.3：掘巷期间从巷道掘进开始每隔50m布置一个表面位移测站，用卷尺或测杆监测巷道相对变形量。
[0101]
步骤s2：标定岩石蠕变参数。
[0102]
步骤s2.1：利用udec离散元模拟软件建立直径为50mm高为100mm的单轴压缩数值模型和直径为50mm的巴西劈裂数值模型，模型块体采用软件自带trigon算法生成的三角形块体，三角形块体最大边长不超过5mm，块体单元和接触面均采用mohr-coulomb模型，块体参数：体积模量k＝1.2e9pa、剪切模量g＝0.78e9pa、内摩擦角θz＝30
°
、内聚力cz＝6.3e6pa和抗拉强度tz＝4.4e6pa，接触面参数：法相刚度kn＝2.56e11pa、剪切刚度ks＝1.02e11pa、接触面内摩擦角θj＝30
°
、接触面内聚力cj＝3.2e6pa和接触面抗拉强度tj＝1.2e6pa，通过单轴压缩和巴西劈裂的室内实验与数值模型对比，如图2和3所示，岩石弹性模量和抗压强度误均差小于10％，说明确定的块体参数和接触面参数比较合理。
[0103]
步骤s2.2：建立直径为50mm高为100mm的第一单轴蠕变数值模型，模型块体形状和大小与步骤s2.1相同，块体单元采用软件自带的cvisc蠕变模型，块体蠕变参数：maxwell剪切模量gm＝5.97e8pa、maxwell粘度系数ηm＝2.9e13，kelvin剪切模量gk＝0.78e9pa、kelvin粘度系数ηk＝1.0e3、内摩擦角θz＝30
°
、内聚力cz＝6.3e6pa和抗拉强度tz＝4.4e6pa，接触面模型及参数与步骤s2.1相同，通过单轴蠕变室内实验与数值模型对比，如图4所示，数值模拟的初始蠕变阶段和等速蠕变阶段与室内实验结果匹配情况较好，采用软件自带的cvisc蠕变模型可以较好的模拟岩石初始蠕变阶段和等速蠕变阶段，由于cvisc蠕变模型未考虑岩石的蠕变损伤过程即岩石在恒定载荷作用下其力学参数随时间不断劣化的过程，无法实现加速蠕变阶段。
[0104]
步骤s3：反演蠕变退化方程。建立直径为50mm高为100mm的第二单轴蠕变数值模型，将步骤s2标定的接触面参数和块体蠕变参数进行赋参，将自定义初始块体蠕变参数退化方程嵌入到主程序中，对比数值模型与步骤s1.2获得的岩石蠕变曲线，由图4可知，数值模拟再现了岩石的初始蠕变、等速蠕变和加速蠕变阶段，由图5可知，本方法也可实现蠕变过程中的微裂隙扩展，大致分为三个阶段：初始蠕变阶段的缓慢扩展、等速蠕变阶段的匀速扩展和加速蠕变阶段的快速扩展，当微裂纹不断扩展并贯通时，就形成了宏观的变形破坏行为即破裂失稳，大量的室内单轴蠕变声发射实验也验证了这一点。
[0105]
步骤s3中蠕变参数初始退化方程：
[0106][0107]
式中t为试件的寿命；t0为短期强度时间，反演后确定t0＝0.04s；cz(t)为蠕变时间t时刻的块体内聚力；σ
cd
为蠕变损伤阈值，反演后确定σ
cd
＝6.72mpa；σ1为块体实际受到的最大主应力；k为蠕变常数，反演后确定k＝0.05。
[0108]
步骤s4：构建巷道离散元数值仿真模型。
[0109]
步骤s4.1：根据步骤s1.1获取的巷道围岩赋存状况包括巷道埋深、巷道周围采掘情况和巷道三倍宽度范围内的岩体倾角、厚度、岩性、层理结构，建立巷道数值仿真几何模型如图6所示，数值模拟模型尺寸为：x方向250m，y方向50m，巷道围岩采用块体边长不超过0.2m的三角形块体细化，模型上部边界垂直应力按深度350m、容重25kn/m3考虑，为8.75mpa，将步骤s2确定的各岩层蠕变参数进行初次赋参，进行模型的初始平衡。
[0110]
步骤s4.2：将步骤s3获得的蠕变参数退化方程嵌入到数值模型主程序中，并采用试错迭代的方法，调整蠕变参数，使得数值模型与步骤s1.3的巷道围岩变形量随时间的演化规律相匹配，得到与工程实际相匹配的蠕变参数。
[0111]
步骤s4.3：步骤s4.2调整后的蠕变参数对步骤s4.1的模型进行重新赋参，将待研究巷道围岩离散为规则网格，并对各网格进行唯一识别标记，将唯一识别数据标记在各网格处备用如图7所示。
[0112]
步骤s5：建立岩体蠕变失稳判据。
[0113]
步骤5.1：建立直径为2m高为4m的第三单轴蠕变数值模型，模型块体形状和大小与步骤s4.1相同，蠕变参数采用步骤4.2确定的参数，通过fish语言将步骤3的蠕变参数退化方程嵌入到数值模型主程序中。
[0114]
步骤5.2：模型采用分级加载方式，试件单轴抗压强度σc＝12.5mpa，在同一试件上依次加载8.5、9.5、10.5、11.5、12.5mpa，每个应力水平加载为2h，直至试件发生加速蠕变破坏，通过fish语言统计单轴蠕变模拟过程中的裂隙演化情况包括张拉裂隙数量、剪切裂隙数量、张拉裂隙长度和剪切裂隙长度随时间的变化情况。
[0115]
步骤5.3：以蠕变过程中微裂隙发育密度表征岩石蠕变损伤(di)，以单轴蠕变数值模拟过程中试件发生加速蠕变时的临界损伤d
min
＝59.5％作为岩体蠕变失稳判据。
[0116]
步骤s6：巷道围岩蠕变失稳预测。
[0117]
步骤6.1：根据步骤s1.1收集的矿井实际情况，对步骤s4.3经离散网格化处理的巷道离散元数值仿真模型，进行模拟开挖，通过fish语言统计巷道蠕变条件下各网格内张拉裂隙和剪切裂隙时空演化规律；
[0118]
步骤s6.2：将步骤s6.1监测的蠕变张拉裂隙和剪切裂隙情况代入到公式计算得到各网格内岩石蠕变损伤di如图8所示，将岩石蠕变损伤di与蠕变失稳判据d
min
进行比对,若di》d
min
则该网格内岩体发生蠕变破坏，否者稳定，从而确定巷道围岩失稳时间和破坏范围如图9所示。
[0119]
本发明的有益效果：
[0120]
1、本发明基于udec离散元模拟软件，模拟了岩石的初始蠕变阶段、等速蠕变阶段和加速蠕变阶段，实现了岩石蠕变过程中的微裂隙扩展，使得从细观角度揭示巷道围岩蠕变失稳过程及机理成为现实，并通过建立岩体蠕变失稳判据，可对巷道围岩蠕变失稳时间和破坏范围进行定量预测，为高地应力软岩巷道蠕变控制提供理论基础。尤其适用于岩土工程、煤矿井下、隧道工程中确定围岩蠕变破坏程度。
[0121]
2、通过本发明可得到岩石的初始蠕变阶段、等速蠕变阶段和加速蠕变阶段，并从微裂隙扩展角度揭示巷道围岩蠕变失稳过程及机理，为巷道支护设计提供理论依据。
[0122]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0123]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。