CN106813973A - 岩体幂函数型细观时效破裂三维模型的构建方法 - Google Patents

Abstract

本发明涉及一种岩体幂函数型细观时效破裂三维模型的构建方法，包括考虑弯扭贡献因子的岩体细观颗粒粘结应力三维模式、考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式、考虑弯扭贡献效应且带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则、以及考虑阻尼效应的细观颗粒线性接触三维模型的构建过程。本发明适应于三维应力空间条件下应力和裂纹扩展速度之间的关系符合幂函数型的这类岩体，对于这类深部岩体工程在三维应力条件下的围岩长期稳定性预测、评价以及优化设计提供技术支持。

Description

岩体幂函数型细观时效破裂三维模型的构建方法

技术领域

本发明涉及工程岩体三维细观时效破裂分析技术领域，具体涉及一种岩体幂函数型细观时效破裂三维模型的构建方法。

背景技术

深部岩体工程开挖后的失稳与破坏往往不是在开挖后立刻发生的，一般都存在着明显的变形破裂时效性和灾变(如岩爆、大变形等)的滞后性，严重危害工程的施工安全与长期运营。目前，在细观方面的时效力学研究成果相对较少。《深埋大理岩破裂扩展时间效应的颗粒流模拟》一文对锦屏大理岩破裂的时间效应进行了试验和二维数值分析(岩石力学与工程学报，2011，Vol.30 No.10:1989-1996)；《锦屏大理岩蠕变损伤演化细观力学特征的数值模拟研究》一文应用二维蠕变细观力学模型对锦屏大理岩短期和长期强度特征进行了数值研究(岩土力学，2013，Vol.34 No.12:3601-3608)。这类模型是以指数型构建的二维驱动应力和裂纹扩展速度间的关系，用来描述岩石细观层面上的二维时效破裂，适用于平面条件下应力和裂纹扩展速度之间符合指数表达方式的岩体。另外，这类模型还存在如下不足之处：(1)颗粒间的剪切破裂准则是一条与平行粘结正应力平行的水平直线，也即这种剪切破裂准则与平行粘结正应力状态无关，只要平行粘结剪切应力大于或等于固定平行粘结剪切破裂强度，颗粒间即可发生剪切破裂，无法体现岩体中不同平行粘结正应力具有不同平行粘结剪切破裂强度的客观事实；(2)没有考虑粘结力矩的差异作用对接触破坏的影响，将粘结力矩的贡献度对不同岩性的影响均视为一致；(3)对于应力和裂纹扩展速度之间不合符指数表达方式的岩体，这类模型缺乏适应性；(4)对于必须考虑三维应力条件下的深部岩体工程问题，该类二维模型因不能描述三维应力对岩体时效变形破坏的影响，也同样缺乏适应性。

发明内容

本发明的目的在于提供一种岩体幂函数型细观时效破裂三维模型的构建方法，本发明适应于三维应力空间条件下应力和裂纹扩展速度之间的关系符合幂函数型的这类岩体，对于这类深部岩体工程在三维应力条件下的围岩长期稳定性预测、评价以及优化设计提供技术支持。

为解决上述技术问题，本发明公开的一种岩体幂函数型细观时效破裂三维模型的构建方法，包括如下步骤：

步骤1：设定岩体细观颗粒粘结接触的三维几何参数量包括三维粘结面积、三维粘结惯性矩和三维粘结极惯性矩；其中，R^(a),R^(b)分别为三维粘结接触两端的颗粒半径，粘结单位厚度为1时的三维粘结面积、粘结单位厚度为1时的三维粘结惯性矩和粘结单位厚度为1时的三维粘结极惯性矩分别通过公式(2)、公式(3)、公式(4)来确定：

其中：为岩体细观颗粒三维粘结半径，为三维粘结直径乘数或半径乘数，A为三维粘结面积，I为三维粘结惯性矩，J为三维粘结极惯性矩；

步骤2：利用岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt，通过三维幂函数形式计算岩体细观颗粒粘结直径，公式(5)来确定；

其中：为判断三维岩体细观颗粒开始时效劣化衰减时的应力阀值，为岩体细观颗粒三维粘结拉伸强度，为考虑扭矩贡献因子的岩体细观颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力，β₁为控制幂函数整体变化的岩体内部细观颗粒三维粘结时效劣化系数，β₂为控制幂函数上标部分变化的岩体内部细观颗粒三维粘结时效劣化系数，为岩体细观颗粒三维粘结随时间劣化衰减的直径，为岩体细观颗粒三维粘结未衰减时的直径；

步骤3：根据步骤2中的公式(5)，设定岩体细观颗粒三维粘结直径的幂函数型时效衰减因子，见公式(6)：

其中：β为岩体细观颗粒三维粘结直径的时效衰减因子，A'、I'、J'、分别为岩体内部细观颗粒三维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数(粘结直径乘数指粘结直径(或粘结半径)与粘结两端最小颗粒直径(或半径)的比值)，Δt为岩体时效衰减劣化的时间增量， A、I、J、分别为岩体内部细观颗粒三维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数；

步骤4：将上述步骤1的公式(1)和步骤3中的公式(6)，代入步骤1中的公式(2)、公式(3)和公式(4)中得到岩体细观颗粒三维粘结几何参数时效劣化衰减模式，该岩体细观颗粒三维粘结几何参数时效劣化衰减模式，即是在三维情况下，岩体细观颗粒粘结直径随着时间增加而不断劣化衰减，三维粘结的面积、惯性矩和极惯性矩也随着时间增加而不断劣化衰减，分别见公式(7)、公式(8)和公式(9)；

其中：A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，A'、I'、J'分别表示为岩体细观颗粒三维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子；

步骤5：依次计算待构建三维模型中的第j个至第k个岩体细观颗粒粘结包含时间效应的三维粘结法向弯矩增量、切向扭矩增量，具体计算方法为，由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δt_c，通过如下公式(10)、公式(11)、公式(12)、公式(13)，确定三维岩体细观颗粒粘结法向增量位移三维岩体细观颗粒粘结切向st方向的增量位移三维岩体细观颗粒粘结切向ss方向的增量位移确定三维岩体细观颗粒粘结法向相对转角三维岩体细观颗粒粘结切向ss方向的相对转角三维岩体细观颗粒粘结切向st方向的相对转角再结合步骤4中的公式(8)和公式(9)以及步骤3中的公式(6)，确定三维岩体细观颗粒粘结切向st方向的扭矩增量、切向ss方向的扭矩增量以及三维岩体细观颗粒粘结法向弯矩增量，见如下公式(14)、公式(15)以及公式(16)；

其中：ff、j、k是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的最末标号值，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的绝对运动速度，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的角速度，n_n、n_ss、n_st分别为三维岩体细观颗粒粘结接触的法向单位向量、切向ss方向的单位向量、切向st方向的单位向量，ss和st为同一平面上相互垂直的两个方向的代号，分别为三维岩体细观颗粒粘结法向的位移增量、切向ss方向的位移增量、切向st方向的位移增量，I、J分别为岩体细观颗粒三维粘结未衰减时的惯性矩、极惯性矩，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，分别为三维岩体细观颗粒粘切向ss方向的扭矩增量值、切向st方向的扭矩增量值，为三维岩体细观颗粒粘法向弯矩增量值，三维岩体细观颗粒粘的弯矩和扭矩按右手螺旋法则，确定其矢量方向；

步骤6：根据步骤203中的公式(7)～公式(9)、步骤204中的公式(10)～公式(13)以及步骤202中的公式(6)，并通过公式(17)、公式(20)、公式(23)、公式(24)计算第i个岩体细观颗粒三维粘结接触的粘结法向力、切向力、法向弯矩、切向扭矩

第i个岩体细观颗粒三维粘结接触的粘结法向力：

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力：

第i个岩体细观颗粒三维粘结接触的粘结切向st方向力：

第i个岩体细观颗粒三维粘结接触的粘结切向合力：

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩：

第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩：

第i个岩体细观颗粒三维粘结接触的粘结法向弯矩：

第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩：

其中：为第i个岩体细观颗粒三维粘结接触的粘结法向力、为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向st方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向合力，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向弯矩，为第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向位移增量，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子，ff为包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号，+＝为加法自反运算符，-＝为减法自反运算符；

步骤7：考虑三维岩体细观颗粒粘结法向扭矩对岩体细观颗粒三维粘结正应力的贡献程度，在三维粘结正应力计算公式中设置扭矩贡献因子考虑三维岩体细观颗粒粘结切向弯矩对岩体细观颗粒三维粘结剪应力的贡献程度，在三维粘结剪应力计算公式中设置弯矩贡献因子根据岩体细观颗粒三维粘结正应力公式和岩体细观颗粒三维粘结剪应力公式同时将这两个公式中A、I、J以及用A'、I'、J'及替换，然后将步骤4中的公式(7)～公式(9)以及步骤3中的公式(6)代入，可获得包含幂函数型时间效应且考虑弯扭贡献效应的岩体细观颗粒三维粘结法向正应力和三维粘结剪应力计算公式，分别见公式(25)和公式(26)；

步骤8：将步骤7中包含幂函数型时间效应且考虑弯扭贡献效应的代入公式(27)，可确定带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则，该准则包含幂函数型时间效应和弯扭贡献效应，该准则用于判断岩体细观颗粒三维粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒三维粘结应力中包含了幂函数型时间效应和弯扭贡献效应；

其中：f^s为摩尔-库伦细观颗粒三维粘结剪切时效破裂准则，fⁿ为摩尔-库伦细观颗粒三维粘结拉伸时效破裂准则，为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒三维粘结剪应力，为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体细观颗粒三维粘结正应力，f^s表示岩体细观颗粒三维粘结剪切破裂准则，f^s大于等于0表示三维粘结剪切破裂，小于0表示三维粘结未发生剪切破裂；fⁿ表示岩体细观颗粒三维粘结拉伸破裂准则，fⁿ大于等于0表示三维粘结拉伸破裂，小于0表示三维粘结未发生拉伸破裂；

步骤9：如果步骤8中的公式(27)中的f^s或fⁿ大于等于0，表明三维粘结发生了破裂，此后岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达；如果步骤8中的公式(27)中的f^s和fⁿ都小于0，表明三维粘结未破裂，继续循环步骤2至8，计算、更新、判断岩体细观颗粒接触的三维粘结状态，直至岩体不产生新的三维粘结破裂或者三维粘结破裂加速发展而形成宏观破坏，循环终止。

本发明的有益效果：

(1)本发明中的模型结构包括考虑弯扭贡献因子的岩体细观颗粒粘结应力三维模式、考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式、考虑弯扭贡献效应且带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则、考虑阻尼效应的细观颗粒线性接触三维模型等，这四部分构成了完整的模型结构系统，本发明所构建的模型中设置了考虑弯扭贡献因子的岩体细观颗粒粘结应力三维模式，不仅在岩体细观颗粒粘结正应力三维计算公式中设置了弯矩贡献因子，而且在岩体颗粒粘结剪应力三维计算公式中设置了扭矩贡献因子。这种模型结构及构建方法，不仅考虑了弯矩对细观颗粒粘结正应力的贡献程度，考虑了扭矩对颗粒粘结剪应力的贡献程度，而且还考虑了弯矩和扭矩的贡献程度对岩体长期强度的影响，适合描述一类岩体的细观力学破裂三维空间问题。

(2)本发明中构建了考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式，包括在岩体细观颗粒粘结时效劣化衰减时，设置了幂函数型与考虑弯扭贡献因子的粘结应力相关的细观颗粒粘结三维劣化衰减模式，设置了细观颗粒粘结直径随时间逐步劣化衰减的三维与考虑弯扭贡献因子的粘结应力相关幂函数型模式，设置了细观颗粒粘结面积、惯性矩和极惯性矩等时效劣化衰减三维模式；按照这种三维幂函数型时效劣化衰减模式估算岩体细观颗粒粘结破裂的初始时间步长。这种三维幂函数型构建模式适合描述空间状态下一类深部岩体的三维细观力学时效破裂机制和响应规律。

(3)本发明中在所构建的幂函数型时效破裂三维模型中，嵌入考虑弯扭贡献效应且带拉伸截止限的的摩尔库伦细观颗粒粘结时效破裂准则。在岩体细观颗粒粘结时效破裂时，采用所嵌入的考虑弯扭贡献效应且带拉伸截止限的的摩尔库伦细观颗粒粘结时效破裂准则来判断；在该准则的细观颗粒粘结应力中包含了幂函数型时间效应，并在颗粒粘结正应力中设置了弯矩贡献因子，在颗粒粘结剪应力中设置了扭矩贡献因子。这种模型结构中粘结时效破裂准则的构建方法，不仅可以描述与细观颗粒粘结正应力相关时效剪切破裂强度的差异，还可以对细观时效拉伸破裂进行合理的表达，且考虑了弯矩和扭矩贡献程度对细观粘结时效破裂的影响，符合空间状态下一类岩体三维细观时效破裂模式。

(4)本发明中在所构建的幂函数型时效破裂三维模型中，嵌入考虑阻尼效应的细观颗粒线性接触三维模型结构，在岩体时效破裂后，通过指定三维接触参考距离设定岩体细观颗粒空间接触距离，设置考虑岩体细观颗粒空间受力变形的三维接触模式以及在岩体细观颗粒之间设置考虑三向滑动摩擦面力的耦合作用模式，同时设置三维接触的空间阻尼模式，可合理描述三维应力状态下一类深部工程岩体在时效破裂后的颗粒空间运动和受力特征。

本发明所提出的一种岩体幂函数型时效破裂三维模型及构建方法，适应于三维应力空间条件下应力和裂纹扩展速度之间的关系符合幂函数型的这类岩体，对于这类深部岩体工程在三维应力条件下的围岩长期稳定性预测、评价以及优化设计提供技术支持。

附图说明

图1为本发明模型中细观颗粒与颗粒接触示意图；

图2为本发明模型中细观颗粒与刚性墙接触示意图；

图3为本发明模型中细观颗粒空间重叠状态示意图；

图4为本发明模型中细观颗粒刚度计算三维示意图；

图5为本发明模型中细观颗粒粘结线性切向力和切向位移示意图；

图6为本发明模型中细观颗粒粘接触本构物理模型示意图；

图7为本发明模型中细观颗粒线性粘结三维结构示意图；

图8为本发明模型中考虑弯扭贡献效应且带拉伸截止限的的摩尔库伦细观颗粒粘结时效破裂准则示意图；

图9为本发明模型中细观颗粒粘结直径(或半径)时效劣化衰减示意图；

图10为本发明模型中细观颗粒三维接触面的力和力矩分布量示意图；

图11为本发明模型中细观颗粒三维接触面的法向和切向向量示意图；

图12为本发明模型构建流程示意图；

图13为本发明模型试样图；

图14为本发明模型蠕变位移与时间关系曲线图。

其中，其中：1—两颗粒的中心距离d，2—岩体细观颗粒间的半接触距离，3—岩体细观颗粒间的半参考距离g_r，4—岩体细观颗粒a的坐标，5—岩体细观颗粒b的坐标，6—岩体细观颗粒表面接触距离的中心坐标，7—岩体细观颗粒表面接触距离g_s，8—岩体细观颗粒间的单位接触法向量，9—岩体细观颗粒a的半径R^a，10—岩体细观颗粒b的半径R^b，11—岩体细观颗粒接触点的接触重叠量U，12—代表b(岩体细观颗粒或是边界墙体)的刚度(法向、切向刚度统称)k^b，13—代表a(岩体细观颗粒或是边界墙体)的刚度(法向、切向刚度统称)k^a，14—岩体细观颗粒接触点的等效刚度，15—总位移增量ΔU_i，16—初始法向力增量值，17—初始接触力矢量和，18—初始切向力增量值，19—所构建的幂函数型时效破裂三维模型法向位移增量Δδ_n，20—所构建的幂函数型时效破裂三维模型切向位移增量Δδ_s，21—岩体细观颗粒粘结拉伸强度值22—岩体细观颗粒粘结法向刚度23—岩体细观颗粒接触点的法向刚度K_n，24—岩体细观颗粒粘结切向刚度25—岩体细观颗粒粘结剪切强度，25.1—为岩体细观颗粒粘结粘聚力，25.2—岩体细观颗粒粘结内摩擦角26—岩体细观颗粒接触点的切向刚度K_s，27—岩体细观颗粒线性接触滑动摩擦系数，28—为岩体细观颗粒线性接触法向阻尼系数β_n，29—岩体细观颗粒线性接触切向阻尼系数β_s，30—为岩体细观颗粒粘结直径乘数31—岩体细观颗粒粘结直径32—考虑弯扭贡献效应且带拉伸截止限的摩尔库伦时效破裂准则，33—第i个接触的包含幂函数时间效应且考虑扭矩贡献因子的岩体细观颗粒粘结剪应力34—第i个接触的包含幂函数时间效应且考虑弯矩贡献因子的岩体细观颗粒粘结正应力35—岩体细观颗粒粘结时效衰减的半径36—岩体细观颗粒粘结时效衰减的直径37—岩体细观颗粒粘结未衰减时的直径38—岩体细观颗粒粘结未衰减时的半径39—岩体细观颗粒粘结弯矩方向矢量，40—岩体细观颗粒粘结扭矩方向矢量，41—岩体细观颗粒粘结切向方向的力矢量，42—岩体细观颗粒粘结法向方向的力矢量，43—岩体细观颗粒粘结直径，44—岩体细观颗粒粘结单位厚度(一般取值为1)，45—岩体细观颗粒粘结切向ss方向的分量，46—岩体细观颗粒粘结切向st方向的分量，47—岩体细观颗粒接触面的法向向量n_n，48—岩体细观颗粒线性粘结接触面。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

本发明所设计的岩体幂函数型细观时效破裂三维模型的构建方法，包括如下步骤：

步骤1：设定岩体细观颗粒粘结接触的三维几何参数量包括三维粘结面积、三维粘结惯性矩和三维粘结极惯性矩；其中，R^(a),R^(b)分别为三维粘结接触两端的颗粒半径，粘结单位厚度为1时的三维粘结面积、粘结单位厚度为1时的三维粘结惯性矩和粘结单位厚度为1时的三维粘结极惯性矩分别通过公式(2)、公式(3)、公式(4)来确定：

其中：为岩体细观颗粒三维粘结半径，为三维粘结直径乘数或半径乘数，A为三维粘结面积，I为三维粘结惯性矩，J为三维粘结极惯性矩；

步骤2：利用岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt，通过三维幂函数形式计算岩体细观颗粒粘结直径，公式(5)来确定；

其中：为判断三维岩体细观颗粒开始时效劣化衰减时的应力阀值，为岩体细观颗粒三维粘结拉伸强度，为考虑扭矩贡献因子的岩体细观颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力，β₁为控制幂函数整体变化的岩体内部细观颗粒三维粘结时效劣化系数，β₂为控制幂函数上标部分变化的岩体内部细观颗粒三维粘结时效劣化系数，为岩体细观颗粒三维粘结随时间劣化衰减的直径，为岩体细观颗粒三维粘结未衰减时的直径；

步骤3：根据步骤2中的公式(5)，设定岩体细观颗粒三维粘结直径的幂函数型时效衰减因子，见公式(6)：

其中：β为岩体细观颗粒三维粘结直径的时效衰减因子，A'、I'、J'、分别为岩体内部细观颗粒三维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数(粘结直径乘数指粘结直径(或粘结半径)与粘结两端最小颗粒直径(或半径)的比值)，Δt为岩体时效衰减劣化的时间增量， A、I、J、分别为岩体内部细观颗粒三维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数；

步骤4：将上述步骤1的公式(1)和步骤3中的公式(6)，代入步骤1中的公式(2)、公式(3)和公式(4)中得到岩体细观颗粒三维粘结几何参数时效劣化衰减模式，该岩体细观颗粒三维粘结几何参数时效劣化衰减模式，即是在三维情况下，岩体细观颗粒粘结直径随着时间增加而不断劣化衰减，三维粘结的面积、惯性矩和极惯性矩也随着时间增加而不断劣化衰减，分别见公式(7)、公式(8)和公式(9)；

其中：A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，A'、I'、J'分别表示为岩体细观颗粒三维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子；

步骤5：依次计算待构建三维模型中的第j个至第k个岩体细观颗粒粘结包含时间效应的三维粘结法向弯矩增量、切向扭矩增量，具体计算方法为，由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δt_c，通过如下公式(10)、公式(11)、公式(12)、公式(13)，确定三维岩体细观颗粒粘结法向增量位移三维岩体细观颗粒粘结切向st方向的增量位移三维岩体细观颗粒粘结切向ss方向的增量位移确定三维岩体细观颗粒粘结法向相对转角三维岩体细观颗粒粘结切向ss方向的相对转角三维岩体细观颗粒粘结切向st方向的相对转角再结合步骤4中的公式(8)和公式(9)以及步骤3中的公式(6)，确定三维岩体细观颗粒粘结切向st方向的扭矩增量、切向ss方向的扭矩增量以及三维岩体细观颗粒粘结法向弯矩增量，见如下公式(14)、公式(15)以及公式(16)；

其中：ff、j、k是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的最末标号值，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的绝对运动速度，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的角速度，n_n、n_ss、n_st分别为三维岩体细观颗粒粘结接触的法向单位向量、切向ss方向的单位向量、切向st方向的单位向量，ss和st为同一平面上相互垂直的两个方向的代号，分别为三维岩体细观颗粒粘结法向的位移增量、切向ss方向的位移增量、切向st方向的位移增量，I、J分别为岩体细观颗粒三维粘结未衰减时的惯性矩、极惯性矩，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，分别为三维岩体细观颗粒粘切向ss方向的扭矩增量值、切向st方向的扭矩增量值，为三维岩体细观颗粒粘法向弯矩增量值，三维岩体细观颗粒粘的弯矩和扭矩按右手螺旋法则，确定其矢量方向；

步骤6：根据步骤203中的公式(7)～公式(9)、步骤204中的公式(10)～公式(13)以及步骤202中的公式(6)，并通过公式(17)、公式(20)、公式(23)、公式(24)计算第i个岩体细观颗粒三维粘结接触的粘结法向力、切向力、法向弯矩、切向扭矩

第i个岩体细观颗粒三维粘结接触的粘结法向力：

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力：

第i个岩体细观颗粒三维粘结接触的粘结切向st方向力：

第i个岩体细观颗粒三维粘结接触的粘结切向合力：

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩：

第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩：

第i个岩体细观颗粒三维粘结接触的粘结法向弯矩：

第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩：

其中：为第i个岩体细观颗粒三维粘结接触的粘结法向力、为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向st方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向合力，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向弯矩，为第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向位移增量，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子，ff为包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号，+＝为加法自反运算符，-＝为减法自反运算符；

步骤7：考虑三维岩体细观颗粒粘结法向扭矩对岩体细观颗粒三维粘结正应力的贡献程度，在三维粘结正应力计算公式中设置扭矩贡献因子考虑三维岩体细观颗粒粘结切向弯矩对岩体细观颗粒三维粘结剪应力的贡献程度，在三维粘结剪应力计算公式中设置弯矩贡献因子根据岩体细观颗粒三维粘结正应力公式和岩体细观颗粒三维粘结剪应力公式同时将这两个公式中A、I、J以及用A'、I'、J'及替换，然后将步骤4中的公式(7)～公式(9)以及步骤3中的公式(6)代入，可获得包含幂函数型时间效应且考虑弯扭贡献效应的岩体细观颗粒三维粘结法向正应力和三维粘结剪应力计算公式，分别见公式(25)和公式(26)；

步骤8：将步骤7中包含幂函数型时间效应且考虑弯扭贡献效应的代入公式(27)，可确定带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则，该准则包含幂函数型时间效应和弯扭贡献效应，该准则用于判断岩体细观颗粒三维粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒三维粘结应力中包含了幂函数型时间效应和弯扭贡献效应；

其中：f^s为摩尔-库伦细观颗粒三维粘结剪切时效破裂准则，fⁿ为摩尔-库伦细观颗粒三维粘结拉伸时效破裂准则，为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒三维粘结剪应力，为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体细观颗粒三维粘结正应力，f^s表示岩体细观颗粒三维粘结剪切破裂准则，f^s大于等于0表示三维粘结剪切破裂，小于0表示三维粘结未发生剪切破裂；fⁿ表示岩体细观颗粒三维粘结拉伸破裂准则，fⁿ大于等于0表示三维粘结拉伸破裂，小于0表示三维粘结未发生拉伸破裂；

步骤9：如果步骤8中的公式(27)中的f^s或fⁿ大于等于0，表明三维粘结发生了破裂，此后岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达；如果步骤8中的公式(27)中的f^s和fⁿ都小于0，表明三维粘结未破裂，继续循环步骤2至8，计算、更新、判断岩体细观颗粒接触的三维粘结状态，直至岩体不产生新的三维粘结破裂或者三维粘结破裂加速发展而形成宏观破坏，循环终止。

上述技术方案的步骤2中，岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt的确定方法为：通过采用考虑弯扭贡献效应的三维粘结时效劣化衰减的幂函数型模式，由每个时间步内的三维粘结首次衰减破裂所损耗的时间来确定，也即通过第一个三维粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个三维粘结破裂所需要的计算循环次数来估算初始时间步长，见公式其中， 为第i个接触的岩体内部细观颗粒粘结直径乘数，n_c为第一个岩体内部细观颗粒粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体内部细观颗粒三维粘结拉伸强度和剪切强度状态下的时效劣化因子，为判断岩体内部细观颗粒三维粘结开始时效劣化衰减时的应力阀值，为岩体内部细观颗粒三维粘结拉伸强度，为考虑弯扭贡献因子的颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力。

上述技术方案中，所述岩体内部细观颗粒三维粘结拉伸强度状态下的时效劣化因子β_σ和岩体内部细观颗粒三维粘结剪切强度状态下的时效劣化因子β_τ的确定包括如下步骤；其中，这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的三维粘结破裂标号；

步骤1000：由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δt_c，通过公式确定三维粘结接触的法向相对转角通过公式确定三维粘结切向ss方向的相对转角通过公式确定三维粘结切向st方向的相对转角通过公式确定三维粘结法向增量位移通过公式确定三维粘结切向ss方向的增量位移通过公式确定三维粘结切向st方向的增量位移通过公式确定三维粘结接触的弯矩增量，通过公式确定三维粘结切向st方向的扭矩增量，通过公式确定三维粘结切向ss方向的扭矩增量。

步骤1001：根据步骤1000中的公式通过公式确定三维粘结法向力；根据步骤100中的公式和公式通过公式和确定三维粘结切向st方向力、切向ss方向力，并通过确定三维粘结切向合力；根据步骤1000中的公式和公式通过公式确定三维粘结法向弯矩；根据步骤1000中的公式和公式以及和公式通过公式和公式确定三维粘结切向st方向扭矩、切向ss方向扭矩，并通过确定三维粘结切向合扭矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；

步骤1002：通过公式确定三维粘结法向正应力，通过公式确定三维粘结剪应力，将这两个公式中A、I、J以及用A'、I'、J'及替换，然后将步骤,4中的公式(7)～公式(9)以及步骤3中的公式(6)代入，可获得包含幂函数型时间效应和扭矩贡献因子的三维粘结法向正应力计算公式和包含幂函数型时间效应和弯矩贡献因子的三维粘结剪应力计算公式步骤1003：将代入公式并令β＝β_σ；将代入公式并令β＝β_τ，据此，可分别由公式和根据牛顿迭代法或斯蒂芬森加速迭代方法或二等法求解这两个方程，可分别得到对应拉伸强度状态下的三维粘结时效劣化因子β_σ以及对应剪切强度状态下的三维粘结时效劣化因子β_τ。

上述技术方案中，岩体细观颗粒三维粘结发生破裂后，岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达，用于描述岩体时效破裂后细观颗粒的三维应力和三维变形及空间运动规律，考虑阻尼效应的三维线性接触模型的构建包括如下步骤：

步骤2000：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个线性接触端a、二维线性接触端b(颗粒与颗粒、颗粒与墙)的中心坐标，在三维情况下，通过公式(28)计算两者中心距离：

其中：d为三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离，为三维线性接触端a的坐标，为三维线性接触端b的坐标；

步骤2001：在所构建的三维细观时效破裂模型中，岩体中颗粒间每个接触点的单位向量通过公式(29)计算，如果是颗粒与颗粒之间的接触，则利用步骤2000中得到的三维线性接触两端的中心点坐标(其中三维线性接触端a的坐标三维线性接触端b的坐标)及中心距离d计算岩体中颗粒间每个接触点的单位向量；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触点的单位向量：

其中：n_i为三维线性接触的单位矢量，为三维线性接触端b的方向向量，为三维线性接触端a的方向向量，n_wall为约束墙的方向向量；

步骤2002：在所构建的三维细观时效破裂模型中，岩体破裂后，每一个接触点的接触重叠量U，通过步骤2000计算的三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离d，以及三维线性接触两端(a端、b端)的颗粒半径R^a、R^b，再利用公式(30)来确定；通过设定颗粒三维线性接触参考距离g_r，并结合公式(31)，确定颗粒三维线性接触的距离g_s：

g_s＝|U|-g_r (31)

步骤2003：在所构建的三维细观时效破裂模型中，确定岩体中细观颗粒三维线性接触点法向、切向等效刚度，可利用接触两端颗粒实体或是墙体的刚度k^a，k^b等效代替为接触点的法向刚度和切向刚度，按公式(32)计算：

其中：K_n、K_s为等效的法向刚度和切向刚度，为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度，为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度；

步骤2004：在所构建的三维细观时效破裂模型中，确定岩体中接触两端颗粒间的相对运动速度，可利用公式(33)、公式(34)来计算，其中e_pqz为Ricci指标置换符号，按照公式(35)来计算：

其中：V_p与V_q等价，V_p与V_q为岩体中细观颗粒三维线性接触两端颗粒间的相对运动速度，p、q为指标等价符号，p＝1，q＝1表示颗粒与颗粒接触，p＝2，q＝2时表示颗粒与墙接触，为颗粒与颗粒或是颗粒与墙的接触b端单元的速度，为颗粒与颗粒或是颗粒与墙的接触a端单元的速度，为颗粒与颗粒或是颗粒与墙的接触a端单元的角速度，为颗粒与颗粒或是颗粒与墙的接触b端单元的角速度，为颗粒与颗粒或是颗粒与墙的接触a端的位移，为颗粒与颗粒或是颗粒与墙的接触b端的位移，为位移指标变换的中间过渡符号，表示指标符号为p时颗粒-颗粒或者颗粒-墙的接触a端的速度，表示指标符号为q时颗粒-颗粒或者颗粒-墙的接触a端的速度，表示指标符号为p时颗粒-颗粒或者颗粒-墙的接触b端的速度，表示指标符号为q时颗粒-颗粒或者颗粒-墙的接触b端的速度(只有a端和b端两个接触端)；

步骤2005：在所构建的三维细观时效破裂模型中，对于时间步长Δt的取值，可以通过公式(38)得到最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳定，通过公式(39)、公式(40)、公式(42)、公式(43)确定每个线性接触的总位移增量、法向位移增量和切向位移增量：

R＝min(R_a,R_b) (36)

ΔU_p1＝V_p1Δt (39)

Δδ_ss＝Δδ_sn_ss (42)

Δδ_st＝Δδ_sn_st (43)

其中：m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k^平为岩体细观颗粒系统平移刚度，k^转为岩体细观颗粒系统旋转刚度；ΔU_p1为岩体细观颗粒三维线性接触的总位移增量，Δδ_n、为岩体细观颗粒三维线性接触的法向位移增量，Δδ_s、为岩体细观颗粒三维线性接触的切向位移增量，V_p1与V_q1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，Δδ_ss、Δδ_st为切向位移Δδ_s在ss方向、st方向的分量，三者之间的关系为：n_ss、n_st为岩体细观颗粒三维线性接触面的切向ss方向、st方向的单位向量，p1,q1为张量指标变换符号。

步骤2006：在所构建的三维细观时效破裂模型中，可由公式(31)判定岩体中颗粒表面接触允许存在的最大距离，通过公式(44)计算法向和切向位移更新因子，另外，岩体细观颗粒三维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量ss方向分量的更新是采用前一步的切向位移增量ss方向分量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量st方向分量的更新是采用前一步的切向位移增量st方向分量与更新因子α的乘积获得：

其中：(g_s)₀为模型计算初始时刻的表面接触距离，g_s为岩体细观颗粒接触的距离，α为位移更新因子；

步骤2007：在所构建的三维细观时效破裂模型中，三维法向线性力的更新采取了相对矢量累加(M_l＝1)和绝对矢量累加(M_l＝0)模式，通过公式(45)计算，而切向线性力的更新采用了三维接触滑动来表示，通过公式(48)、公式(49)计算；

其中：k_n、k_s为三维线性接触法向线性刚度、切向线性刚度，g_s为模型颗粒在一定荷载下的表面接触距离，Δδ_n和Δδ_s分别为三维线性接触法向位移增量和切向位移增量，为三维线性接触的法向接触力，为初始法向力增量值和切向力增量值，为三维线性接触的切向接触力，为三维线性接触切向线性力在st方向、ss方向的分量，三者之间的关系为: 为岩体细观颗粒未滑动时的静摩擦力,为颗粒滑动摩擦力，通过摩擦系数μ与乘积得到，Δδ_st、Δδ_ss分别为三维线性接触切向增量Δδ_s在st方向位移增量和ss方向位移增量，Δδ_s、Δδ_st、Δδ_ss三者之间的关系为：

步骤2008：在所构建的三维细观时效破裂模型中，法向阻尼力采用全法向模式M_d＝{0,2}和无拉伸模式M_d＝{1,3}两种，通过公式(50)、公式(51)计算；切向阻尼力采用全剪切模式M_d＝{0,1}和滑剪模式M_d＝{2,3}，按照公式(52)、公式(53)来计算；

其中：分别为三维线性接触的法向和切向阻尼力，β_n为三维线性接触的法向阻尼系数、β_s为三维线性接触的切向阻尼系数，k_n为三维线性接触的法向线性刚度、k_s为三维线性接触的切向线性刚度，分别为三维线性接触的法向速率和三维线性接触的切向速率，m_c为等效颗粒质量，m⁽¹⁾为颗粒与颗粒之间的第一接触端的颗粒质量，m⁽²⁾为颗粒与颗粒之间的第二接触端的颗粒质量，F^d为三维线性接触的总阻尼力，分别为三维线性接触切向阻尼在ss方向、st方向的分量，三者之间的关系为: 为三维线性接触的法向接触力， 表示三维线性接触切向ss方向的速率，表示三维线性接触切向st方向的速率，三者之间的关系为：

本发明所构建的岩体幂函数型细观时效破裂三维模型，所述三维模型包括考虑弯扭贡献因子的岩体细观颗粒粘结应力三维模式、考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式、考虑弯扭贡献效应且带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则、以及考虑阻尼效应的细观颗粒线性接触三维模型。

上述技术方案中，所述岩体幂函数型细观时效破裂三维模型适用于三维颗粒离散元分析方法、三维颗粒不连续变形分析方法、三维颗粒流形元分析方法。

上述技术方案中，考虑扭矩对岩体内部细观颗粒三维粘结法向正应力的贡献程度，所述弯扭贡献因子的岩体细观颗粒粘结应力三维模式为岩体内部细观颗粒粘结正应力三维计算公式中设置了扭矩贡献因子

考虑弯矩对岩体内部细观颗粒三维粘结剪应力的贡献程度在扭贡献因子的岩体细观颗粒粘结应力三维模式中的岩体内部细观颗粒粘结剪应力三维计算公式中设置了弯矩贡献因子在上述公式中，为岩体内部细观颗粒三维粘结半径，为用于确定扭矩在应力中的贡献程度的扭矩贡献因子， 为用于确定弯矩在应力中的贡献程度弯矩贡献因子，I为岩体内部细观颗粒三维粘结的惯性矩，J为岩体内部细观颗粒三维粘结的极惯性矩，A为岩体内部细观颗粒三维粘结面积，为第i个接触的岩体内部细观颗粒三维粘结正应力，为第i个接触的岩体内部细观颗粒三维粘结剪应力，分别为第i个接触的岩体内部细观颗粒三维粘结法向力、切向合力、切向合扭矩和法向弯矩，其中，岩体内部细观颗粒三维粘结法向力式中，为岩体内部细观颗粒三维粘结法向的位移增量，为岩体内部细观颗粒三维粘结法向刚度，+＝符号为加法自反运算符；

切向合力式中，分别为岩体内部细观颗粒粘结切向ss方向力和岩体内部细观颗粒粘结切向st方向力，其中，ss和st为同一平面上相互垂直的两个方向的代号；

切向合扭矩式中，分别为岩体内部细观颗粒粘结切向ss方向扭矩和岩体内部细观颗粒粘结切向st方向扭矩，ss和st为同一平面上相互垂直的两个方向的代号；

法向弯矩式中，为岩体内部细观颗粒粘结法向的相对转角增量，为岩体内部细观颗粒粘结切向刚度，-＝符号为减法自反运算符；

上述技术方案中，所述考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式包括在岩体细观颗粒粘结时效劣化衰减时，设置了幂函数型与考虑弯扭贡献因子的粘结应力相关的细观颗粒粘结三维劣化衰减模式，见幂函数型更新率：

式中，为判断岩体内部细观颗粒三维粘结开始时效劣化衰减时的应力阀值，为岩体内部细观颗粒三维粘结拉伸强度，为考虑弯扭贡献因子的颗粒三维粘结应力比，β₁为控制幂函数整体变化的岩体内部细观颗粒三维粘结时效劣化系数，β₂为控制幂函数上标部分变化的岩体内部细观颗粒三维粘结时效劣化系数，为岩体细观颗粒三维粘结应力；

在考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式中设置了三维与考虑弯矩贡献因子的粘结应力相关幂函数型模式，该幂函数型模式中的岩体内部细观颗粒粘结直径随时间逐步劣化衰减，见粘结直径公式式中，为岩体内部细观颗粒粘结随时间劣化衰减的直径，为岩体内部细观颗粒三维粘结未衰减时的直径，Δt为岩体时效衰减劣化的时间增量；

在考虑弯扭贡献因子的细观颗粒粘结时效劣化衰减的三维幂函数型模式中设置了岩体内部细观颗粒粘结面积、惯性矩和极惯性矩的时效劣化衰减三维模式，分别见粘结单位厚度为1时的三维粘结面积计算公式粘结单位厚度为1时的三维粘结惯性矩计算公式和粘结单位厚度为1时的三维粘结极惯性矩计算公式其中，β为岩体内部细观颗粒三维粘结直径的时效衰减因子，β的计算公式为其中，A'、I'、J'、分别为岩体内部细观颗粒三维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数，Δt为岩体时效衰减劣化的时间增量， A、I、J、分别为岩体内部细观颗粒三维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数；

同时按照这种三维幂函数型时效劣化衰减模式估算岩体内部细观颗粒三维粘结破裂的初始时间步长增量，见公式其中， 为第i个接触的岩体内部细观颗粒粘结直径乘数，n_c为第一个岩体内部细观颗粒粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体内部细观颗粒三维粘结拉伸强度和剪切强度状态下的时效劣化因子，∞为无穷大，岩体细观颗粒三维粘结拉伸强度状态下的时效劣化因子β_σ和岩体细观颗粒三维粘结剪切强度状态下的时效劣化因子β_τ可分别由公式：

根据迭代法(牛顿迭代法或斯蒂芬森加速迭代方法)或二等法求解这两个方程获得，其中为岩体内部细观颗粒三维粘结拉伸强度，为岩体内部细观颗粒三维粘结的粘聚力，为岩体内部细观颗粒三维粘结的内摩擦角,F_σ为β_σ的函数，F_τ为β_τ的函数，π为圆周率。

上述技术方案中，所述考虑弯扭贡献效应且带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则包括公式：

其中，f^s为摩尔-库伦细观颗粒三维粘结剪切时效破裂准则，fⁿ为摩尔-库伦细观颗粒三维粘结拉伸时效破裂准则，为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒三维粘结剪应力，为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体细观颗粒三维粘结正应力，分别为岩体内部细观颗粒粘结拉伸强度、抗剪强度，为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体内部细观颗粒粘结正应力，的计算公式为 的计算公式为在所述考虑弯扭贡献效应且带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则的细观颗粒三维粘结应力中包含了幂函数型时间效应，见岩体内部细观颗粒粘结直径的时效衰减因子计算公式在岩体内部颗粒三维粘结剪应力中设置了弯矩贡献因子f^s大于等于0为岩体内部细观颗粒三维粘结剪切破裂，小于0为岩体内部细观颗粒三维粘结未发生剪切破裂；fⁿ大于等于0为岩体内部细观颗粒三维粘结拉伸破裂，fⁿ小于0为岩体内部细观颗粒三维粘结未发生拉伸破裂。

上述技术方案中，所述考虑阻尼效应的细观颗粒线性接触三维模型是指在岩体细观颗粒时效破裂后，通过三维接触参考距离g_r设定了岩体内部细观颗粒空间接触距离，见岩体内部细观颗粒空间接触距计算公式其中，为接触端a的坐标，为接触端b的坐标，R^a、R^b分别为岩体内部细观接触端a的颗粒半径和接触端b的颗粒半径；

在考虑阻尼效应的细观颗粒线性接触三维模型中设置了考虑岩体内部细观颗粒空间变形的线性三维接触模式，在岩体内部细观颗粒之间设置了考虑三向滑动摩擦面力的耦合作用模式，岩体内部细观颗粒空间变形的线性三维接触法向线性力计算公式取M_l＝1为相对矢量累加模式，取M_l＝0为绝对矢量累加模式，岩体内部细观颗粒空间变形的线性三维接触切向线性力计算公式为和其中，k_n、k_s为岩体内部细观颗粒空间变形的线性三维接触法向、切向线性刚度，Δδ_n为岩体内部颗粒线性接触的法向位移增量(Δδ_s为岩体内部颗粒线性接触的切向位移增量)，为岩体内部颗粒线性接触的初始法向力增量值和切向力增量值，为切向线性力在岩体内部细观颗粒粘结切向ss方向和岩体内部细观颗粒粘结切向st方向的分量，为颗粒滑动摩擦力，通过摩擦系数μ与乘积得到，为颗粒未滑动时的静摩擦力，Δδ_st、Δδ_ss分别为岩体内部细观颗粒粘结切向ss方向的位移增量和岩体内部细观颗粒粘结切向st方向的位移增量；

同时设置三维接触的空间阻尼模式，其中法向阻尼采用全法向模式M_d＝{0,2}和无拉伸模式M_d＝{1,3}两种，通过公式计算，其中，F^*为岩体内部颗粒线性接触的全法向阻尼力，表达式为 为岩体内部细观颗粒空间变形的线性三维接触法向线性力，m_c为等效颗粒质量，按公式计算；

切向阻尼采用全剪切模式M_d＝{0,1}以及滑动和剪切模式M_d＝{2,3}，根据

公式

公式

进行计算，其中，为ss方向的速率，ss为岩体内部颗粒三维线性接触面内部的某一个方向(粒与颗粒之间三维粘结破裂后，颗粒与颗粒之间的接触变为三维线性接触)，为st方向的速率，st为岩体内部颗粒三维线性接触面内部的另一个方向，该方向与ss方向相互垂直，为岩体内部颗粒线性接触的法向阻尼力，β_n为岩体内部颗粒线性接触的法向阻尼系数、β_s为岩体内部颗粒线性接触的切向阻尼系数，k_n为岩体内部颗粒线性接触的法向线性刚度、k_s为岩体内部颗粒线性接触的切向线性刚度，为岩体内部颗粒线性接触的法向速率(为的合速率，称为线性接触的切向速率。三者之间的关系为：)，m_c为等效颗粒质量，m⁽¹⁾为岩体内部颗粒与颗粒接触的第一接触端的颗粒质量，m⁽²⁾为岩体内部颗粒与颗粒接触的第二接触端的颗粒质量，F^d为总阻尼力，为线性接触的法向阻尼力，为线性接触的切向阻尼力，F^d为的合力，称为线性接触总阻尼力。三者关系为： 为切向阻尼在岩体内部细观颗粒粘结切向ss方向和岩体内部细观颗粒粘结切向st方向的分量。

下面以深部岩体为实例，结合附图详述本发明模型的数值实现的具体过程，请参阅实例附图说明中的图13至图14以及模型附图说明中的图1至图12，来理解本发明模型的数值实现步骤及效果：

步骤1：采用C++编程语言，并结合fish语言，根据本发明的模型结构构建流程图(图12)，在数值平台上实现了岩体幂函数型时效破裂三维模型。

步骤2：初步确定岩体时效破裂模型的细观参数

粒径比Rratio、线性接触法向刚度kn(图6)、线性接触切向刚度ks(图6)、颗粒密度ba_rho、颗粒接触模量b_Ec、粘结法向刚度pb_kn(图6)、粘结切向刚度pb_ks(图6)、粘结模型pb_Ec、颗粒的摩擦系数ba_fric、粘结拉伸强度pb_sn_mean、粘结拉伸强度的标准差pb_sn_sdev、粘聚力平均值pb_coh_mean、粘聚力标准差pb_coh_sdev、粘结半径乘数gamma(图7)、粘结弯矩贡献因子beta_sigma、粘结扭矩贡献因子beta_shear、法向阻尼系数Apfan(图6)、切向阻尼系数Apfas(图6)以及内摩擦角pb_phi(图8)等19个参数，参数具体值见表一。

步骤3：生成岩体模型

根据高斯分布或weibull分布确定模型的粘结拉伸强度和粘聚力分布，通过均匀随机函数分布法确定颗粒的粒径分布；通过各向同性应力调整法及自适应动态膨胀法，调整颗粒的位置，减少颗粒重叠量；通过悬浮颗粒删除法，删除孤立颗粒，提高模型样本的整体性，减少缺陷模型的生成。最后赋予模型材料粘结强度参数，生成具有真实岩体结构的岩体结构图模型。岩体模型的直径为50mm、高度为100mm(图13)。

步骤4：精确标定本发明中模型的细观物理力学参数

通过室内单轴和三轴压缩试验得到的应力-应变曲线，确定岩体的宏观弹性模量峰值强度σ_p，以及泊松比通过优化方法，使岩体单、三轴压缩模型的应力-应变曲线与室内试验的应力-应变以及宏观变形参数和强度参数吻合，获得本发明所构建模型的细观物理力学参数。

步骤5：岩体时效力学参数标定

对岩体进行一系列不同应力强度比条件下的时效力学试验，得到不同应力强度比条件下岩体变形随时间演化的曲线。通过参数拟合法，匹配实际岩体的时效变形过程，确定控制岩体细观颗粒粘结时效劣化的两个参β₁、β₂。

步骤6：岩体时效力学数值试验

在荷载一定的条件下，进行岩体三维时效力学数值模拟试验，获得了岩体时效变形破坏的演化规律(图14)。

表一：本发明模型的参数名称及值

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (4)

1.一种岩体幂函数型细观时效破裂三维模型的构建方法，包括如下步骤：

步骤1：设定岩体细观颗粒粘结接触的三维几何参数量包括三维粘结面积、三维粘结惯性矩和三维粘结极惯性矩；其中，R^(a),R^(b)分别为三维粘结接触两端的颗粒半径，粘结单位厚度为1时的三维粘结面积、粘结单位厚度为1时的三维粘结惯性矩和粘结单位厚度为1时的三维粘结极惯性矩分别通过公式(2)、公式(3)、公式(4)来确定：

$\stackrel{\‾}{R}=\stackrel{\‾}{\mathrm{\λ}}min(R^{\left(a\right)},R^{\left(b\right)})---\left(1\right)$

$A=\mathrm{\π}{\stackrel{\‾}{R}}^2---\left(2\right)$

$I=\frac{1}{4}\mathrm{\π}{\stackrel{\‾}{R}}^4---\left(3\right)$

$J=\frac{1}{2}\mathrm{\π}{\stackrel{\‾}{R}}^4---\left(4\right)$

其中：为岩体细观颗粒三维粘结半径，为三维粘结直径乘数或半径乘数，A为三维粘结面积，I为三维粘结惯性矩，J为三维粘结极惯性矩；

步骤2：利用岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt，通过三维幂函数形式计算岩体细观颗粒粘结直径，公式(5)来确定；

${\stackrel{\&OverBar;}{D}}^{\&prime;}=\stackrel{\&OverBar;}{D}-\left\{\begin{array}{c}0,\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{aa}\\ {\mathrm{\&beta;}}_1{\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c}\right)}^{{\mathrm{\&beta;}}_2}\mathrm{\&Delta;}t,{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{aa}\&le;\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c\end{array}\right.---\left(5\right)$

其中：为判断三维岩体细观颗粒开始时效劣化衰减时的应力阀值，为岩体细观颗粒三维粘结拉伸强度，为考虑扭矩贡献因子的岩体细观颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力，β₁为控制幂函数整体变化的岩体内部细观颗粒三维粘结时效劣化系数，β₂为控制幂函数上标部分变化的岩体内部细观颗粒三维粘结时效劣化系数，为岩体细观颗粒三维粘结随时间劣化衰减的直径，为岩体细观颗粒三维粘结未衰减时的直径；

步骤3：根据步骤2中的公式(5)，设定岩体细观颗粒三维粘结直径的幂函数型时效衰减因子，见公式(6)：

$\mathrm{\&beta;}={\stackrel{\&OverBar;}{D}}^{\&prime;}/\stackrel{\&OverBar;}{D}={\stackrel{\&OverBar;}{R}}^{\&prime;}/\stackrel{\&OverBar;}{R}={\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}^{\&prime;}/\stackrel{\&OverBar;}{\mathrm{\&lambda;}}=1-\left\{\begin{array}{c}0,\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{aa}\\ {\mathrm{\&beta;}}_1{\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c}\right)}^{{\mathrm{\&beta;}}_2}\frac{\mathrm{\&Delta;}t}{\stackrel{\&OverBar;}{D}},{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{aa}\&le;\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c\end{array}\right.---\left(6\right)$

其中：β为岩体细观颗粒三维粘结直径的时效衰减因子，A'、I'、J'、分别为岩体内部细观颗粒三维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数(粘结直径乘数指粘结直径(或粘结半径)与粘结两端最小颗粒直径(或半径)的比值)，Δt为岩体时效衰减劣化的时间增量， A、I、J、分别为岩体内部细观颗粒三维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩、粘结直径乘数；

步骤4：将上述步骤1的公式(1)和步骤3中的公式(6)，代入步骤1中的公式(2)、公式(3)和公式(4)中得到岩体细观颗粒三维粘结几何参数时效劣化衰减模式，该岩体细观颗粒三维粘结几何参数时效劣化衰减模式，即是在三维情况下，岩体细观颗粒粘结直径随着时间增加而不断劣化衰减，三维粘结的面积、惯性矩和极惯性矩也随着时间增加而不断劣化衰减，分别见公式(7)、公式(8)和公式(9)；

$A^{\&prime;}=\mathrm{\&pi;}{\stackrel{\&OverBar;}{R}}^{\&prime;2}={\mathrm{\&beta;}}^{2}A---\left(7\right)$

$I^{\&prime;}=\frac{1}{4}\mathrm{\&pi;}{\stackrel{\&OverBar;}{R}}^{\&prime;4}={\mathrm{\&beta;}}^{4}I---\left(8\right)$

$J^{\&prime;}=\frac{1}{2}\mathrm{\&pi;}{\stackrel{\&OverBar;}{R}}^{\&prime;4}={\mathrm{\&beta;}}^{4}J---\left(9\right)$

其中：A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，A'、I'、J'分别表示为岩体细观颗粒三维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子；

步骤5：依次计算待构建三维模型中的第j个至第k个岩体细观颗粒粘结包含时间效应的三维粘结法向弯矩增量、切向扭矩增量，具体计算方法为，由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δt_c，通过如下公式(10)、公式(11)、公式(12)、公式(13)，确定三维岩体细观颗粒粘结法向增量位移三维岩体细观颗粒粘结切向st方向的增量位移三维岩体细观颗粒粘结切向ss方向的增量位移确定三维岩体细观颗粒粘结法向相对转角三维岩体细观颗粒粘结切向ss方向的相对转角三维岩体细观颗粒粘结切向st方向的相对转角再结合步骤4中的公式(8)和公式(9)以及步骤3中的公式(6)，确定三维岩体细观颗粒粘结切向st方向的扭矩增量、切向ss方向的扭矩增量以及三维岩体细观颗粒粘结法向弯矩增量，见如下公式(14)、公式(15)以及公式(16)；

${\left({\mathrm{\&Delta;U}}_i^n\right)}_{ff}=({\stackrel{\&CenterDot;}{X}}_i^b-{\stackrel{\&CenterDot;}{X}}_i^a)n_n{\mathrm{\&Delta;t}}_c---\left(10\right)$

${\left({\mathrm{\&Delta;U}}_i^{ss}\right)}_{ff}=({\stackrel{\&CenterDot;}{X}}_i^b-{\stackrel{\&CenterDot;}{X}}_i^a)n_{ss}{\mathrm{\&Delta;t}}_c---\left(11\right)$

${\left({\mathrm{\&Delta;U}}_i^{st}\right)}_{ff}=({\stackrel{\&CenterDot;}{X}}_i^b-{\stackrel{\&CenterDot;}{X}}_i^a)n_{st}{\mathrm{\&Delta;t}}_c---\left(12\right)$

$\left\{\begin{array}{c}{\left({\mathrm{\&Delta;\&theta;}}_i^n\right)}_{ff}=({\mathrm{\&omega;}}_i^b-{\mathrm{\&omega;}}_i^a)n_n{\mathrm{\&Delta;t}}_c\\ {\left({\mathrm{\&Delta;\&theta;}}_i^{ss}\right)}_{ff}=({\mathrm{\&omega;}}_i^b-{\mathrm{\&omega;}}_i^a)n_{ss}{\mathrm{\&Delta;t}}_c\\ {\left({\mathrm{\&Delta;\&theta;}}_i^{st}\right)}_{ff}=({\mathrm{\&omega;}}_i^b-{\mathrm{\&omega;}}_i^a)n_{st}{\mathrm{\&Delta;t}}_c\end{array}\right.---\left(13\right)$

${\left(\mathrm{\&Delta;}{{\stackrel{\&OverBar;}{M}}_i}^{st}\right)}_{ff}={\stackrel{\&OverBar;}{k}}_n{\mathrm{\&beta;}}^{4}I{\left({{\mathrm{\&Delta;\&theta;}}_i}^{st}\right)}_{ff}---\left(14\right)$

${\left(\mathrm{\&Delta;}{{\stackrel{\&OverBar;}{M}}_i}^{ss}\right)}_{ff}={\stackrel{\&OverBar;}{k}}_n{\mathrm{\&beta;}}^{4}I{\left({{\mathrm{\&Delta;\&theta;}}_i}^{ss}\right)}_{ff}---\left(15\right)$

${\left(\mathrm{\&Delta;}{\stackrel{\&OverBar;}{M}}_i^n\right)}_{ff}={\stackrel{\&OverBar;}{k}}_s{\mathrm{\&beta;}}^{4}J{\left({{\mathrm{\&Delta;\&theta;}}_i}^n\right)}_{ff}---\left(16\right)$

其中：ff、j、k是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒粘结衰减后未破裂的最末标号值，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的绝对运动速度，分别为第i个三维岩体细观颗粒粘结接触的a端和b端的角速度，n_n、n_ss、n_st分别为三维岩体细观颗粒粘结接触的法向单位向量、切向ss方向的单位向量、切向st方向的单位向量，ss和st为同一平面上相互垂直的两个方向的代号，分别为三维岩体细观颗粒粘结法向的位移增量、切向ss方向的位移增量、切向st方向的位移增量，I、J分别为岩体细观颗粒三维粘结未衰减时的惯性矩、极惯性矩，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，分别为三维岩体细观颗粒粘切向ss方向的扭矩增量值、切向st方向的扭矩增量值，为三维岩体细观颗粒粘法向弯矩增量值，三维岩体细观颗粒粘的弯矩和扭矩按右手螺旋法则，确定其矢量方向；

步骤6：根据步骤203中的公式(7)～公式(9)、步骤204中的公式(10)～公式(13)以及步骤202中的公式(6)，并通过公式(17)、公式(20)、公式(23)、公式(24)计算第i个岩体细观颗粒三维粘结接触的粘结法向力、切向力、法向弯矩、切向扭矩

第i个岩体细观颗粒三维粘结接触的粘结法向力：

${\left({\stackrel{\&OverBar;}{F}}_i^n\right)}_{ff}+={\stackrel{\&OverBar;}{k}}_n{\mathrm{\&beta;}}^{2}A{\left({\mathrm{\&Delta;U}}_i^n\right)}_{ff}---\left(17\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力：

${\left({{\stackrel{\&OverBar;}{F}}_i}^{ss}\right)}_{ff}+={\stackrel{\&OverBar;}{k}}_s{\mathrm{\&beta;}}^{2}A{\left({\mathrm{\&Delta;U}}_i^{ss}\right)}_{ff}---\left(18\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向st方向力：

${\left({{\stackrel{\&OverBar;}{F}}_i}^{st}\right)}_{ff}+={\stackrel{\&OverBar;}{k}}_s{\mathrm{\&beta;}}^{2}A{\left({\mathrm{\&Delta;U}}_i^{st}\right)}_{ff}---\left(19\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向合力：

${\left({{\stackrel{\&OverBar;}{F}}_i}^s\right)}_{ff}=\sqrt{{\left({\left({{\stackrel{\&OverBar;}{F}}_i}^{st}\right)}_{ff}\right)}^2+{\left({\left({{\stackrel{\&OverBar;}{F}}_i}^{ss}\right)}_{ff}\right)}^2}---\left(20\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩：

${\left({{\stackrel{\&OverBar;}{M}}_i}^{ss}\right)}_{ff}-={\stackrel{\&OverBar;}{k}}_n{\mathrm{\&beta;}}^{4}I{\left({{\mathrm{\&Delta;\&theta;}}_i}^{ss}\right)}_{ff}---\left(21\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩：

${\left({{\stackrel{\&OverBar;}{M}}_i}^{st}\right)}_{ff}-={\stackrel{\&OverBar;}{k}}_n{\mathrm{\&beta;}}^{4}I{\left({{\mathrm{\&Delta;\&theta;}}_i}^{st}\right)}_{ff}---\left(22\right)$

第i个岩体细观颗粒三维粘结接触的粘结法向弯矩：

${\left({{\stackrel{\&OverBar;}{M}}_i}^n\right)}_{ff}-={\stackrel{\&OverBar;}{k}}_s{\mathrm{\&beta;}}^{4}J{\left({{\mathrm{\&Delta;\&theta;}}_i}^n\right)}_{ff}---\left(23\right)$

第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩：

${\left({{\stackrel{\&OverBar;}{M}}_i}^s\right)}_{ff}=\sqrt{{\left({\left({{\stackrel{\&OverBar;}{M}}_i}^{st}\right)}_{ff}\right)}^2+{\left({\left({{\stackrel{\&OverBar;}{M}}_i}^{ss}\right)}_{ff}\right)}^2}---\left(24\right)$

其中：为第i个岩体细观颗粒三维粘结接触的粘结法向力、为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向st方向力、为第i个岩体细观颗粒三维粘结接触的粘结切向合力，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向弯矩，为第i个岩体细观颗粒三维粘结接触的粘结切向合扭矩，为第i个岩体细观颗粒三维粘结接触的粘结法向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向ss方向位移增量，为第i个岩体细观颗粒三维粘结接触的粘结切向st方向位移增量，为三维岩体细观颗粒粘结法向刚度，为三维岩体细观颗粒粘切向刚度，A、I、J分别为岩体细观颗粒三维粘结未衰减时的粘结面积、粘结惯性矩、粘结极惯性矩，β为岩体细观颗粒三维粘结直径的时效衰减因子，ff为包含时间效应的岩体细观颗粒粘结衰减后未破裂的初始标号，+＝为加法自反运算符，-＝为减法自反运算符；

步骤7：考虑三维岩体细观颗粒粘结法向扭矩对岩体细观颗粒三维粘结正应力的贡献程度，在三维粘结正应力计算公式中设置扭矩贡献因子考虑三维岩体细观颗粒粘结切向弯矩对岩体细观颗粒三维粘结剪应力的贡献程度，在三维粘结剪应力计算公式中设置弯矩贡献因子根据岩体细观颗粒三维粘结正应力公式和岩体细观颗粒三维粘结剪应力公式同时将这两个公式中A、I、J以及用A'、I'、J'及替换，然后将步骤4中的公式(7)～公式(9)以及步骤3中的公式(6)代入，可获得包含幂函数型时间效应且考虑弯扭贡献效应的岩体细观颗粒三维粘结法向正应力和三维粘结剪应力计算公式，分别见公式(25)和公式(26)；

${{\left(\stackrel{\&OverBar;}{{\mathrm{\&sigma;}}_i}\right)}^{\&prime;}}_{ff}=\frac{1}{{\mathrm{\&pi;\&beta;}}^2{\stackrel{\&OverBar;}{R}}^2}(-{\left({\stackrel{\&OverBar;}{F_i}}^n\right)}_{ff}+\frac{4{\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{\mathrm{\&sigma;}}\left|{\left({\stackrel{\&OverBar;}{M_i}}^s\right)}_{ff}\right|}{\mathrm{\&beta;}\stackrel{\&OverBar;}{R}})---\left(25\right)$

${{\left(\stackrel{\&OverBar;}{{\mathrm{\&tau;}}_i}\right)}^{\&prime;}}_{ff}=\frac{1}{{\mathrm{\&pi;\&beta;}}^2{\stackrel{\&OverBar;}{R}}^2}\left(\right|{\left({\stackrel{\&OverBar;}{F_i}}^s\right)}_{ff}|+\frac{2{\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{\mathrm{\&tau;}}\left|{\left({\stackrel{\&OverBar;}{M_i}}^n\right)}_{ff}\right|}{\mathrm{\&beta;}\stackrel{\&OverBar;}{R}})---\left(26\right)$

步骤8：将步骤7中包含幂函数型时间效应且考虑弯扭贡献效应的代入公式(27)，可确定带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则，该准则包含幂函数型时间效应和弯扭贡献效应，该准则用于判断岩体细观颗粒三维粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒三维粘结应力中包含了幂函数型时间效应和弯扭贡献效应；

$\left\{\begin{array}{c}{f}^s={\left({\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_i\right)}^{\&prime;}-{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_c={\left({\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_i\right)}^{\&prime;}+{\left({\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_i\right)}^{\&prime;}\tan \stackrel{\&OverBar;}{\mathrm{\&phi;}}-\stackrel{\&OverBar;}{c}=\\ \frac{1}{{\mathrm{\&pi;\&beta;}}^2{\stackrel{\&OverBar;}{R}}^2}\left(\right|{\stackrel{\&OverBar;}{F_i}}^s|+\frac{2{\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{\mathrm{\&tau;}}\left|{\stackrel{\&OverBar;}{M_i}}^n\right|}{\mathrm{\&beta;}\stackrel{\&OverBar;}{R}})+\frac{1}{{\mathrm{\&pi;\&beta;}}^2{\stackrel{\&OverBar;}{R}}^2}(-{\stackrel{\&OverBar;}{F_i}}^n+\frac{4{\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{\mathrm{\&sigma;}}\left|{\stackrel{\&OverBar;}{M_i}}^s\right|}{\mathrm{\&beta;}\stackrel{\&OverBar;}{R}})\tan \stackrel{\&OverBar;}{\mathrm{\&phi;}}-\stackrel{\&OverBar;}{c}\\ f^n={\left({\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_i\right)}^{\&prime;}-{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c=\frac{1}{{\mathrm{\&pi;\&beta;}}^2{\stackrel{\&OverBar;}{R}}^2}(-{\stackrel{\&OverBar;}{F_i}}^n+\frac{4{\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{\mathrm{\&sigma;}}\left|{\stackrel{\&OverBar;}{M_i}}^s\right|}{\mathrm{\&beta;}\stackrel{\&OverBar;}{R}})-{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_c\end{array}\right.---\left(27\right)$

其中：f^s为摩尔-库伦细观颗粒三维粘结剪切时效破裂准则，fⁿ为摩尔-库伦细观颗粒三维粘结拉伸时效破裂准则，为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒三维粘结剪应力，为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体细观颗粒三维粘结正应力，f^s表示岩体细观颗粒三维粘结剪切破裂准则，f^s大于等于0表示三维粘结剪切破裂，小于0表示三维粘结未发生剪切破裂；fⁿ表示岩体细观颗粒三维粘结拉伸破裂准则，fⁿ大于等于0表示三维粘结拉伸破裂，小于0表示三维粘结未发生拉伸破裂；

步骤9：如果步骤8中的公式(27)中的f^s或fⁿ大于等于0，表明三维粘结发生了破裂，此后岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达；如果步骤8中的公式(27)中的f^s和fⁿ都小于0，表明三维粘结未破裂，继续循环步骤2至8，计算、更新、判断岩体细观颗粒接触的三维粘结状态，直至岩体不产生新的三维粘结破裂或者三维粘结破裂加速发展而形成宏观破坏，循环终止。

2.根据权利要求1所述的岩体幂函数型细观时效破裂三维模型的构建方法，其特征在于：所述步骤2中，岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt的确定方法为：通过采用考虑弯扭贡献效应的三维粘结时效劣化衰减的幂函数型模式，由每个时间步内的三维粘结首次衰减破裂所损耗的时间来确定，也即通过第一个三维粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个三维粘结破裂所需要的计算循环次数来估算初始时间步长，见公式其中， 为第i个接触的岩体内部细观颗粒粘结直径乘数，n_c为第一个岩体内部细观颗粒粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体内部细观颗粒三维粘结拉伸强度和剪切强度状态下的时效劣化因子，为判断岩体内部细观颗粒三维粘结开始时效劣化衰减时的应力阀值，为岩体内部细观颗粒三维粘结拉伸强度，为考虑弯扭贡献因子的颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力。

3.根据权利要求1所述的岩体幂函数型细观时效破裂三维模型的构建方法，其特征在于：所述岩体内部细观颗粒三维粘结拉伸强度状态下的时效劣化因子β_σ和岩体内部细观颗粒三维粘结剪切强度状态下的时效劣化因子β_τ的确定包括如下步骤；其中，这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的三维粘结破裂标号；

步骤1000：由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δt_c，通过公式确定三维粘结接触的法向相对转角通过公式确定三维粘结切向ss方向的相对转角通过公式确定三维粘结切向st方向的相对转角通过公式确定三维粘结法向增量位移通过公式确定三维粘结切向ss方向的增量位移通过公式确定三维粘结切向st方向的增量位移通过公式确定三维粘结接触的弯矩增量，通过公式确定三维粘结切向st方向的扭矩增量，通过公式确定三维粘结切向ss方向的扭矩增量；

步骤1001：根据步骤1000中的公式通过公式确定三维粘结法向力；根据步骤100中的公式和公式通过公式和确定三维粘结切向st方向力、切向ss方向力，并通过确定三维粘结切向合力；根据步骤1000中的公式和公式通过公式确定三维粘结法向弯矩；根据步骤1000中的公式和公式以及和公式通过公式和公式确定三维粘结切向st方向扭矩、切向ss方向扭矩，并通过确定三维粘结切向合扭矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；

步骤1002：通过公式确定三维粘结法向正应力，通过公式确定三维粘结剪应力，将这两个公式中A、I、J以及用A'、I'、J'及替换，然后将步骤,4中的公式(7)～公式(9)以及步骤3中的公式(6)代入，可获得包含幂函数型时间效应和扭矩贡献因子的三维粘结法向正应力计算公式和包含幂函数型时间效应和弯矩贡献因子的三维粘结剪应力计算公式

步骤1003：将代入公式并令β＝β_σ；将代入公式并令β＝β_τ，据此，可分别由公式和

根据牛顿迭代法或斯蒂芬森加速迭代方法或二等法求解这两个方程，可分别得到对应拉伸强度状态下的三维粘结时效劣化因子β_σ以及对应剪切强度状态下的三维粘结时效劣化因子β_τ。

4.根据权利要求1所述的岩体幂函数型细观时效破裂三维模型的构建方法，其特征在于：岩体细观颗粒三维粘结发生破裂后，岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达，用于描述岩体时效破裂后细观颗粒的三维应力和三维变形及空间运动规律，考虑阻尼效应的三维线性接触模型的构建包括如下步骤：

步骤2000：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个线性接触端a、二维线性接触端b(颗粒与颗粒、颗粒与墙)的中心坐标，在三维情况下，通过公式(28)计算两者中心距离：

$d=\sqrt{{(x_i^b-x_i^a)}^2+{(y_i^b-y_i^a)}^2+{(z_i^b-z_i^a)}^2}---\left(28\right)$

其中：d为三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离，为三维线性接触端a的坐标，为三维线性接触端b的坐标；

步骤2001：在所构建的三维细观时效破裂模型中，岩体中颗粒间每个接触点的单位向量通过公式(29)计算，如果是颗粒与颗粒之间的接触，则利用步骤2000中得到的三维线性接触两端的中心点坐标(其中三维线性接触端a的坐标三维线性接触端b的坐标)及中心距离d计算岩体中颗粒间每个接触点的单位向量；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触点的单位向量：

其中：n_i为三维线性接触的单位矢量，为三维线性接触端b的方向向量，为三维线性接触端a的方向向量，n_wall为约束墙的方向向量；

步骤2002：在所构建的三维细观时效破裂模型中，岩体破裂后，每一个接触点的接触重叠量U，通过步骤2000计算的三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离d，以及三维线性接触两端(a端、b端)的颗粒半径R^a、R^b，再利用公式(30)来确定；通过设定颗粒三维线性接触参考距离g_r，并结合公式(31)，确定颗粒三维线性接触的距离g_s：

g_s＝|U|-g_r (31)

步骤2003：在所构建的三维细观时效破裂模型中，确定岩体中细观颗粒三维线性接触点法向、切向等效刚度，可利用接触两端颗粒实体或是墙体的刚度k^a，k^b等效代替为接触点的法向刚度和切向刚度，按公式(32)计算：

$K_n=\frac{k_n^a{k}_n^b}{k_n^a+k_n^b},K_s=\frac{k_s^a{k}_s^b}{k_s^a+k_s^b}---\left(32\right)$

其中：K_n、K_s为等效的法向刚度和切向刚度，为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度，为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度；

步骤2004：在所构建的三维细观时效破裂模型中，确定岩体中接触两端颗粒间的相对运动速度，可利用公式(33)、公式(34)来计算，其中e_pqz为Ricci指标置换符号，按照公式(35)来计算：

$V_p={\left({\stackrel{\&CenterDot;}{x}}_p^{\left(c\right)}\right)}_b-{\left({\stackrel{\&CenterDot;}{x}}_p^{\left(c\right)}\right)}_a=({\stackrel{\&CenterDot;}{x}}_p^{\left(b\right)}+e_{pqz}{\mathrm{\&omega;}}_q^{\left(b\right)}(x_z^{\left(c\right)}-x_z^{\left(b\right)}\left)\right)-({\stackrel{\&CenterDot;}{x}}_p^{\left(a\right)}+e_{pqz}{\mathrm{\&omega;}}_q^{\left(a\right)}(x_z^{\left(c\right)}-x_z^{\left(a\right)}\left)\right)---\left(33\right)$

$V_q={\left({\stackrel{\&CenterDot;}{x}}_q^{\left(c\right)}\right)}_b-{\left({\stackrel{\&CenterDot;}{x}}_q^{\left(c\right)}\right)}_a=({\stackrel{\&CenterDot;}{x}}_q^{\left(b\right)}+e_{pqz}{\mathrm{\&omega;}}_p^{\left(b\right)}(x_z^{\left(c\right)}-x_z^{\left(b\right)}\left)\right)-({\stackrel{\&CenterDot;}{x}}_q^{\left(a\right)}+e_{pqz}{\mathrm{\&omega;}}_p^{\left(a\right)}(x_z^{\left(c\right)}-x_z^{\left(a\right)}\left)\right)---\left(34\right)$

其中：V_p与V_q等价，V_p与V_q为岩体中细观颗粒三维线性接触两端颗粒间的相对运动速度，p、q为指标等价符号，p＝1，q＝1表示颗粒与颗粒接触，p＝2，q＝2时表示颗粒与墙接触，为颗粒与颗粒或是颗粒与墙的接触b端单元的速度，为颗粒与颗粒或是颗粒与墙的接触a端单元的速度，为颗粒与颗粒或是颗粒与墙的接触a端单元的角速度，为颗粒与颗粒或是颗粒与墙的接触b端单元的角速度，为颗粒与颗粒或是颗粒与墙的接触a端的位移，为颗粒与颗粒或是颗粒与墙的接触b端的位移，为位移指标变换的中间过渡符号，表示指标符号为p时颗粒-颗粒或者颗粒-墙的接触a端的速度，表示指标符号为q时颗粒-颗粒或者颗粒-墙的接触a端的速度，表示指标符号为p时颗粒-颗粒或者颗粒-墙的接触b端的速度，表示指标符号为q时颗粒-颗粒或者颗粒-墙的接触b端的速度(只有a端和b端两个接触端)；

步骤2005：在所构建的三维细观时效破裂模型中，对于时间步长Δt的取值，可以通过公式(38)得到最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳定，通过公式(39)、公式(40)、公式(42)、公式(43)确定每个线性接触的总位移增量、法向位移增量和切向位移增量：

R＝min(R_a,R_b) (36)

$J1=\frac{2}{5}{\mathrm{\&pi;R}}^5---\left(37\right)$

ΔU_p1＝V_p1Δt (39)

${\mathrm{\&Delta;\&delta;}}_n={\mathrm{\&Delta;U}}_{p1}^{nl}=V_{q1}{n}_{q1}{n}_{p1}\mathrm{\&Delta;}t---\left(40\right)$

${\mathrm{\&Delta;\&delta;}}_s={\mathrm{\&Delta;U}}_{p1}^{sl}={\mathrm{\&Delta;U}}_{p1}-{\mathrm{\&Delta;U}}_{p1}^{nl}=V_{p1}\mathrm{\&Delta;}t-V_{q1}{n}_{q1}{n}_{p1}\mathrm{\&Delta;}t---\left(41\right)$

Δδ_ss＝Δδ_sn_ss (42)

Δδ_st＝Δδ_sn_st (43)

其中：m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k^平为岩体细观颗粒系统平移刚度，k^转为岩体细观颗粒系统旋转刚度；ΔU_p1为岩体细观颗粒三维线性接触的总位移增量，Δδ_n、为岩体细观颗粒三维线性接触的法向位移增量，Δδ_s、为岩体细观颗粒三维线性接触的切向位移增量，V_p1与V_q1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，Δδ_ss、Δδ_st为切向位移Δδ_s在ss方向、st方向的分量，三者之间的关系为：n_ss、n_st为岩体细观颗粒三维线性接触面的切向ss方向、st方向的单位向量，p1,q1为张量指标变换符号；

步骤2006：在所构建的三维细观时效破裂模型中，可由公式(31)判定岩体中颗粒表面接触允许存在的最大距离，通过公式(44)计算法向和切向位移更新因子，另外，岩体细观颗粒三维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量ss方向分量的更新是采用前一步的切向位移增量ss方向分量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量st方向分量的更新是采用前一步的切向位移增量st方向分量与更新因子α的乘积获得：

$\mathrm{\&alpha;}=\left\{\begin{array}{c}\frac{g_s}{g_s-{\left(g_s\right)}_0},{\left(g_s\right)}_0>0,g_s<0\\ 1,otherwise\end{array}\right.---\left(44\right)$

其中：(g_s)₀为模型计算初始时刻的表面接触距离，g_s为岩体细观颗粒接触的距离，α为位移更新因子；

步骤2007：在所构建的三维细观时效破裂模型中，三维法向线性力的更新采取了相对矢量累加(M_l＝1)和绝对矢量累加(M_l＝0)模式，通过公式(45)计算，而切向线性力的更新采用了三维接触滑动来表示，通过公式(48)、公式(49)计算；

$F_n^l=\left\{\begin{array}{c}\{\begin{array}{c}{k}_n{g}_s,g_s<0\\ 0,otherwise\end{array},M_l=0\\ min({\left(F_n^l\right)}_o+k_n{\mathrm{\&Delta;\&delta;}}_n,0),M_l=1\end{array}\right.---\left(45\right)$

$F_{st}^{*}={\left(F_s^l\right)}_o-k_s{\mathrm{\&Delta;\&delta;}}_{st}---\left(46\right)$

$F_{ss}^{*}={\left(F_s^l\right)}_o-k_s{\mathrm{\&Delta;\&delta;}}_{ss}---\left(47\right)$

$F_{st}^l=\left\{\begin{array}{c}{F}_{st}^{*},\left|\right|F_{st}^{*}\left|\right|\&le;F_s^{\mathrm{\&mu;}}\\ F_s^{\mathrm{\&mu;}}(F_{st}^{*}/|\left|F_{st}^{*}\right|\left|\right),otherwise\end{array}\right.---\left(48\right)$

$F_{ss}^l=\left\{\begin{array}{c}{F}_{ss}^{*},\left|\right|F_{ss}^{*}\left|\right|\&le;F_s^{\mathrm{\&mu;}}\\ F_s^{\mathrm{\&mu;}}(F_{ss}^{*}/|\left|F_{ss}^{*}\right|\left|\right),otherwise\end{array}\right.---\left(49\right)$

其中：k_n、k_s为三维线性接触法向线性刚度、切向线性刚度，g_s为模型颗粒在一定荷载下的表面接触距离，Δδ_n和Δδ_s分别为三维线性接触法向位移增量和切向位移增量，为三维线性接触的法向接触力，为初始法向力增量值和切向力增量值，为三维线性接触的切向接触力，为三维线性接触切向线性力在st方向、ss方向的分量，三者之间的关系为: 为岩体细观颗粒未滑动时的静摩擦力,为颗粒滑动摩擦力，通过摩擦系数μ与乘积得到，Δδ_st、Δδ_ss分别为三维线性接触切向增量Δδ_s在st方向位移增量和ss方向位移增量，Δδ_s、Δδ_st、Δδ_ss三者之间的关系为：

步骤2008：在所构建的三维细观时效破裂模型中，法向阻尼力采用全法向模式M_d＝{0,2}和无拉伸模式M_d＝{1,3}两种，通过公式(50)、公式(51)计算；切向阻尼力采用全剪切模式M_d＝{0,1}和滑剪模式M_d＝{2,3}，按照公式(52)、公式(53)来计算；

$F_n^d=\left\{\begin{array}{c}\left(2{\mathrm{\&beta;}}_n\sqrt{m_c{k}_n}\right){\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_n,M_d=\{0,2\}\\ \min (F^{*},-F_n^l),M_d=\{1,3\}\end{array}\right.---\left(50\right)$

$F_{ss}^d=\left\{\begin{array}{c}\left(2{\mathrm{\&beta;}}_s\sqrt{m_c{k}_s}\right){\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_{ss},M_d=\{0,1\}\\ 0,M_d=\{2,3\}\end{array}\right.---\left(52\right)$

$F_{st}^d=\left\{\begin{array}{c}\left(2{\mathrm{\&beta;}}_s\sqrt{m_c{k}_s}\right){\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_{st},M_d=\{0,1\}\\ 0,M_d=\{2,3\}\end{array}\right.---\left(53\right)$

其中：分别为三维线性接触的法向和切向阻尼力，β_n为三维线性接触的法向阻尼系数、β_s为三维线性接触的切向阻尼系数，k_n为三维线性接触的法向线性刚度、k_s为三维线性接触的切向线性刚度，分别为三维线性接触的法向速率和三维线性接触的切向速率，m_c为等效颗粒质量，m⁽¹⁾为颗粒与颗粒之间的第一接触端的颗粒质量，m⁽²⁾为颗粒与颗粒之间的第二接触端的颗粒质量，F^d为三维线性接触的总阻尼力，分别为三维线性接触切向阻尼在ss方向、st方向的分量，三者之间的关系为: 为三维线性接触的法向接触力， 表示三维线性接触切向ss方向的速率，表示三维线性接触切向st方向的速率，三者之间的关系为：