1.一种岩体幂函数型细观时效破裂三维模型的构建方法，包括如下步骤：
步骤1：设定岩体细观颗粒粘结接触的三维几何参数量包括三维粘结面积、三维粘结惯性矩和三维粘结极惯性矩；其中，R(a),R(b)分别为三维粘结接触两端的颗粒半径，粘结单位厚度为1时的三维粘结面积、粘结单位厚度为1时的三维粘结惯性矩和粘结单位厚度为1时的三维粘结极惯性矩分别通过公式(2)、公式(3)、公式(4)来确定：
其中：为岩体细观颗粒三维粘结半径，为三维粘结直径乘数或半径乘数，A为三维粘结面积，I为三维粘结惯性矩，J为三维粘结极惯性矩；
步骤2：利用岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt，通过三维幂函数形式计算岩体细观颗粒粘结直径，公式(5)来确定；
其中： 为判断三维岩体细观颗粒开始时效劣化衰减时的应力阀值， 为岩体细观颗粒三维粘结拉伸强度， 为考虑扭矩贡献因子的岩体细观颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力，β1为控制幂函数整体变化的岩体内部细观颗粒三维粘结时效劣化系数，β2为控制幂函数上标部分变化的岩体内部细观颗粒三维粘结时效劣化系数， 为岩体细观颗粒三维粘结随时间劣化衰减的直径， 为岩体细观颗粒三维粘结未衰减时的直径；
步骤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个岩体细观颗粒粘结包含时间效应的三维粘结法向弯矩增量、切向扭矩增量，具体计算方法为，由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δtc，通过如下公式(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端的角速度，nn、nss、nst分别为三维岩体细观颗粒粘结接触的法向单位向量、切向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)，可确定带拉伸截止限的摩尔库伦细观颗粒粘结时效破裂准则，该准则包含幂函数型时间效应和弯扭贡献效应，该准则用于判断岩体细观颗粒三维粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒三维粘结应力中包含了幂函数型时间效应和弯扭贡献效应；
其中：fs为摩尔-库伦细观颗粒三维粘结剪切时效破裂准则，fn为摩尔-库伦细观颗粒三维粘结拉伸时效破裂准则， 为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒三维粘结剪应力， 为第i个接触的含幂函数型时间效应且考虑扭矩贡献因子的岩体细观颗粒三维粘结正应力，fs表示岩体细观颗粒三维粘结剪切破裂准则，fs大于等于0表示三维粘结剪切破裂，小于0表示三维粘结未发生剪切破裂；fn表示岩体细观颗粒三维粘结拉伸破裂准则，fn大于等于0表示三维粘结拉伸破裂，小于0表示三维粘结未发生拉伸破裂；
步骤9：如果步骤8中的公式(27)中的fs或fn大于等于0，表明三维粘结发生了破裂，此后岩体细观颗粒的空间运动模式采用考虑阻尼效应的三维线性接触模型来表达；如果步骤8中的公式(27)中的fs和fn都小于0，表明三维粘结未破裂，继续循环步骤2至8，计算、更新、判断岩体细观颗粒接触的三维粘结状态，直至岩体不产生新的三维粘结破裂或者三维粘结破裂加速发展而形成宏观破坏，循环终止。
2.根据权利要求1所述的岩体幂函数型细观时效破裂三维模型的构建方法，其特征在于：所述步骤2中，岩体细观颗粒三维粘结时效衰减劣化的初始时间步长增量Δt的确定方法为：通过采用考虑弯扭贡献效应的三维粘结时效劣化衰减的幂函数型模式，由每个时间步内的三维粘结首次衰减破裂所损耗的时间来确定，也即通过第一个三维粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个三维粘结破裂所需要的计算循环次数来估 算 初 始 时 间 步 长 ，见 公 式 其 中 ，
为第i个接触的岩体内部细观颗粒粘结直径乘数，nc为第一个岩体
内部细观颗粒粘结破裂所需的循环计算的次数，βσ、βτ分别为岩体内部细观颗粒三维粘结拉伸强度和剪切强度状态下的时效劣化因子， 为判断岩体内部细观颗粒三维粘结开始时效劣化衰减时的应力阀值， 为岩体内部细观颗粒三维粘结拉伸强度， 为考虑弯扭贡献因子的颗粒三维粘结应力比，为岩体细观颗粒三维粘结应力。
3.根据权利要求2所述的岩体幂函数型细观时效破裂三维模型的构建方法，其特征在于：所述岩体内部细观颗粒三维粘结拉伸强度状态下的时效劣化因子βσ和岩体内部细观颗粒三维粘结剪切强度状态下的时效劣化因子βτ的确定包括如下步骤；其中，这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的三维粘结破裂标号；
步骤1000：由三维岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算步Δtc，通过公式 确定三维粘结接触的法向相对转角 通过公式
确定三维粘结切向ss方向的相对转角 通过公式
确定三维粘结切向st方向的相对转角 通过公式
确 定 三 维 粘 结 法 向 增 量 位 移 通 过 公 式
确定三维粘结切向ss方向的增量位移 通过公式
确定三维粘结切向st方向的增量位移 通过公式
确定三维粘结接触的弯矩增量，通过公式 确定
三维粘结切向st方向的扭矩增量，通过公式 确定三维粘结切向ss方
向的扭矩增量；
步骤1001：根据步骤1000中的公式 通过公式
确 定 三 维 粘 结 法 向 力 ；根 据 步 骤 1 0 0 中 的 公 式
和 公 式 通 过 公 式
和 确定三维粘结切向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为三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离， 为三维线性接触端a的坐标， 为三维线性接触端b的坐标；
步骤2001：在所构建的三维细观时效破裂模型中，岩体中颗粒间每个接触点的单位向量通过公式(29)计算，如果是颗粒与颗粒之间的接触，则利用步骤2000中得到的三维线性接触两端的中心点坐标，其中三维线性接触端a的坐标 三维线性接触端b的坐标 及中心距离d计算岩体中颗粒间每个接触点的单位向量；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触点的单位向量：
其中：ni为三维线性接触的单位矢量， 为三维线性接触端b的方向向量， 为三维线性接触端a的方向向量，nwall为约束墙的方向向量；
步骤2002：在所构建的三维细观时效破裂模型中，岩体破裂后，每一个接触点的接触重叠量U，通过步骤2000计算的三维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离d，以及三维线性接触两端，a端、b端的颗粒半径Ra、Rb，再利用公式(30)来确定；通过设定颗粒三维线性接触参考距离gr，并结合公式(31)，确定颗粒三维线性接触的距离gs：
gs＝|U|-gr   (31)
步骤2003：在所构建的三维细观时效破裂模型中，确定岩体中细观颗粒三维线性接触点法向、切向等效刚度，可利用接触两端颗粒实体或是墙体的刚度ka，kb等效代替为接触点的法向刚度和切向刚度，按公式(32)计算：
其中：Kn、Ks为等效的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度；
步骤2004：在所构建的三维细观时效破裂模型中，确定岩体中接触两端颗粒间的相对运动速度，可利用公式(33)、公式(34)来计算，其中epqz为Ricci指标置换符号，按照公式(35)来计算：
其中：Vp与Vq等价，Vp与Vq为岩体中细观颗粒三维线性接触两端颗粒间的相对运动速度，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(Ra,Rb)   (36)
ΔUp1＝Vp1Δt   (39)
Δδss＝Δδsnss   (42)
Δδst＝Δδsnst   (43)
其中：m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k平为岩体细观颗粒系统平移刚度，k转为岩体细观颗粒系统旋转刚度；ΔUp1为岩体细观颗粒三维线性接触的总位移增量，Δδn、 为岩体细观颗粒三维线性接触的法向位移增量，Δδs、 为岩体细观颗粒三维线性接触的切向位移增量，Vp1与Vq1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，Δδss、Δδst为切向位移Δδs在ss方向、st方向的分量，三者之间的关系为：
nss、nst为岩体细观颗粒三维线性接触面的切向ss方向、st
方向的单位向量，p1,q1为张量指标变换符号；
步骤2006：在所构建的三维细观时效破裂模型中，可由公式(31)判定岩体中颗粒表面接触允许存在的最大距离，通过公式(44)计算法向和切向位移更新因子，另外，岩体细观颗粒三维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量ss方向分量的更新是采用前一步的切向位移增量ss方向分量与更新因子α的乘积获得，岩体细观颗粒三维线性接触切向位移增量st方向分量的更新是采用前一步的切向位移增量st方向分量与更新因子α的乘积获得：
其中：(gs)0为模型计算初始时刻的表面接触距离，gs为岩体细观颗粒接触的距离，α为位移更新因子；
步骤2007：在所构建的三维细观时效破裂模型中，三维法向线性力的更新采取了相对矢量累加Ml＝1和绝对矢量累加Ml＝0模式，通过公式(45)计算，而切向线性力的更新采用了三维接触滑动来表示，通过公式(48)、公式(49)计算；
其中：kn、ks为三维线性接触法向线性刚度、切向线性刚度，gs为模型颗粒在一定荷载下的表面接触距离，Δδn和Δδs分别为三维线性接触法向位移增量和切向位移增量， 为三维线性接触的法向接触力， 为初始法向力增量值和切向力增量值， 为三维线性接触的切向接触力， 为三维线性接触切向线性力在st方向、ss方向的分量，三者之间的关系为: 为岩体细观颗粒未滑动时的静摩擦力, 为
颗粒滑动摩擦力，通过摩擦系数μ与 乘积得到，Δδst、Δδss分别为三维线性接触切向增量Δδs在st方向位移增量和ss方向位移增量，Δδs、Δδst、Δδss三者之间的关系为：
步骤2008：在所构建的三维细观时效破裂模型中，法向阻尼力采用全法向模式Md＝{0,
2}和无拉伸模式Md＝{1,3}两种，通过公式(50)、公式(51)计算；切向阻尼力采用全剪切模式Md＝{0,1}和滑剪模式Md＝{2,3}，按照公式(52)、公式(53)来计算；
其中： 分别为三维线性接触的法向和切向阻尼力，βn为三维线性接触的法向阻尼系数、βs为三维线性接触的切向阻尼系数，kn为三维线性接触的法向线性刚度、ks为三维线性接触的切向线性刚度， 分别为三维线性接触的法向速率和三维线性接触的切向速率，mc为等效颗粒质量，m(1)为颗粒与颗粒之间的第一接触端的颗粒质量，m(2)为颗粒与颗粒之间的第二接触端的颗粒质量，Fd为三维线性接触的总阻尼力， 分别为三维线性接触切向阻 尼在ss方向、st方向的分量，三者之间的 关系为:
为三维线性接触的法向接触力， 表示
三维线性接触切向ss方向的速率， 表示三维线性接触切向st方向的速率，三者之间的关系为：