本发明公开了岩体二维细观时效破裂幂函数型模型的构建方法，该模型包括考虑弯矩贡献因子的岩体细观颗粒粘结应力二维模式、考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式、考虑弯矩贡献效应且带拉伸截止限的摩尔‑库伦细观颗粒粘结时效破裂准则、考虑阻尼效应的细观颗粒线性接触二维模型。本发明适应于应力和裂纹扩展速度之间的关系符合幂函数型的这类岩体，对于平面状态下这类深部岩体工程的围岩长期稳定性预测、评价以及优化设计提供技术支持。
1.岩体二维细观时效破裂幂函数型模型的构建方法，其特征在于，包括如下步骤：
步骤1：设定岩体内部细观颗粒粘结接触的几何参数量，包括颗粒粘结面积和颗粒粘结惯性矩，Ra、Rb分别为二维粘结接触a端的颗粒半径、b端的颗粒半径，为岩体细观颗粒粘结直径乘数或半径乘数，在二维情况下，粘结单位厚度为1时的颗粒粘结面积A和粘结惯性矩I分别通过公式(2)、公式(3)来确定：其中：为岩体细观颗粒二维粘结半径，为岩体细观颗粒二维粘结直径乘数或半径乘数，A为岩体细观颗粒二维粘结面积，I为岩体细观颗粒二维粘结惯性矩；
步骤201：利用岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt，通过幂函数形式计算岩体细观颗粒二维粘结时效衰减劣化的直径 由公式(4)来确定；
其中： 为考虑弯矩贡献因子的二维粘结法向应力， 为判断岩体细观颗粒二维粘结开始时效劣化衰减时的应力阈值， 为岩体细观颗粒二维粘结的拉伸强度， 为考虑弯矩贡献因子的二维粘结应力比，β1、β2分别为岩体细观颗粒粘结时效劣化的第一控制参数、第二控制参数， 为岩体细观颗粒二维粘结随时间劣化衰减的直径， 为岩体细观颗粒二维粘结未衰减时的直径，Δt为岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量；
步骤202：根据步骤201中的公式(4)，设定细观颗粒二维粘结直径的时效衰减因子β，见公式(5)：其中： A'、I'、 分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数， A、I、为岩体细观颗粒二维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数；
步骤203：根据步骤1中的公式(1)～公式(3)以及步骤202中的公式(5)，设定岩体细观颗粒二维粘结几何参数时效劣化衰减模式，在平面二维情况下，岩体细观颗粒二维粘结直径随着时间增加而不断劣化衰减，粘结单位厚度为1时的颗粒粘结的面积和惯性矩也随着时间增加而不断劣化衰减，分别见公式(6)、公式(7)；
其中： A'、I'分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩，A、I为岩体细观颗粒二维粘结未衰减时的粘结面积、粘结惯性矩；
步骤204：依次计算第j个至第k个的岩体细观颗粒包含时间效应的二维粘结法向弯矩增量，在二维情况下，由粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δt，通过公式(8)、公式(9)、公式(10)确定第i个岩体细观颗粒二维粘结相对转角 岩体细观颗粒二维粘结法向增量位移 以及岩体细观颗粒二维粘结切向增量位移 再结合步骤203中的公式(7)和步骤202中的公式(5)，确定第i个岩体细观颗粒包含时间效应的二维粘结弯矩增量，具体见公式(11)；
其中：ff、j、k、i是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的初始标号值，ff为中间标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的最末标号值，i为第一个至最后一个二维粘结颗粒标号值， 分别为第i个岩体细观颗粒二维粘结接触的a端和b端的绝对运动速度和角速度，nn、ns分别为岩体细观颗粒二维粘结接触面的法向单位向量和切向单位向量， 分别为岩体细观颗粒二维粘结法向位移增量和切向位移增量， 为岩体细观颗粒二维粘结法向刚度， 为岩体细观颗粒二维粘结弯矩增量；
步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和公式(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒粘结未破裂的粘结接触并包含时间效应的二维粘结法向力、切向力和切向弯矩，通过公式(12)、公式(13)、公式(14)计算第i个岩体细观颗粒二维粘结接触的粘结法向力、切向力和切向弯矩，在二维情况下，通过公式(15)来确定岩体细观颗粒粘结法向弯矩，法向力：
粘结法向力、粘结切向力、包含时间效应的粘结法向弯矩、粘结切向弯矩、粘结法向位移增量和粘结切向位移增量， 为岩体细观颗粒二维粘结切向刚度，+＝为加法自反运算符，-＝为减法自反运算符；
步骤206：设置弯矩贡献因子 考虑弯矩对岩体细观颗粒粘结法向应力的
贡献程度，根据岩体细观颗粒粘结二维正应力计算公式 和二维粘结剪应
力计算公式 同时将这两个公式中A、I以及 用A'、I'及 替换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含幂函数型时间效应和细观颗粒粘结考虑弯矩贡献因子的二维粘结正应力 计算公式和二维粘结剪应力 计算公式，分别见公式(16)和公式(17)，步骤207：将步骤206中包含时间效应的 代入公式(18)，确定考虑弯矩贡献因子且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则，并且依次计算第j个至第k个的岩体细观颗粒二维粘结应力，用于判断岩体细观颗粒粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒粘结应力中包含了幂函数型时间效应和考虑弯矩贡献因子，其中：fs、fn分别为岩体细观颗粒二维粘结的时效剪切破裂准则、时效拉伸破裂准则，为第i个接触的含幂函数型时间效应的二维粘结剪应力， 为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的二维粘结正应力， 分别为岩体细观颗粒二维粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维粘结的粘聚力，为岩体细观颗粒二维粘结的内摩擦角；fs大于等于0表示岩体细观颗粒粘结剪切破裂，小于0表示岩体细观颗粒粘结未发生剪切破裂；fn大于等于0表示岩体细观颗粒粘结拉伸破裂，小于0表示岩体细观颗粒粘结未发生拉伸破裂；
步骤208：当步骤207中的公式(18)中的fs或fn大于等于0时，表明岩体细观颗粒粘结发生破裂，此时岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达；当步骤207中的公式(18)中的fs和fn都小于0时，表明岩体细观颗粒粘结未破裂，继续循环步骤
201至207，计算、更新、判断岩体细观颗粒接触的粘结状态，直至岩体不产生新的细观颗粒粘结破裂或者岩体细观颗粒粘结破裂加速发展而形成宏观破坏，循环终止。
2.根据权利要求1所述的岩体二维细观时效破裂幂函数型模型的构建方法，其特征在于：所述岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt的确定，是采用考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式，由每个时间步内的细观颗粒二维粘结首次衰减破裂所损耗的时间来确定，即通过第一个细观颗粒粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个细观颗粒粘结破裂所需要的计算循环次数来估算初始时间步长增量Δt，见公式 其中， 为第i个接触的岩体细观颗粒二维粘结直径乘数，nc为第一个
岩体细观颗粒二维粘结破裂所需的循环计算的次数，βσ、βτ分别为岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子、二维粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩体细观颗粒粘结数，∞为无穷大。
3.根据权利要求2所述的岩体二维细观时效破裂幂函数型模型的构建方法，其特征在于：所述岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子βσ和岩体细观颗粒二维粘结剪切强度对应的时效劣化因子βτ的确定包括如下步骤，其中这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的细观颗粒二维粘结破裂标号：步骤211：在二维情况下，由岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δtc，通过公式 确定细观颗粒粘结接触的相对转角 通过公式 确定细观颗粒粘结法向增量位移 通
确定细观颗粒接触的粘结法向力；根据步骤211中的公式
确定细观颗粒接触的粘结切向弯矩；通过公式 确定细观颗粒
接触的粘结法向弯矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；
步骤213：在二维情况下，通过公式 确定细观颗粒接触的粘结
正应力，通过公式 确定细观颗粒接触的粘结剪应力，将这两个公式中A、I以及用A'、I'及 替换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得细观颗粒粘结包含幂函数型时间效应和考虑弯矩贡献因子的二维粘结正应力计算公式 和包含幂函数型时间效应的二维粘结剪应力计算公式步骤214：将 代入公式 并令β＝βσ；将
代入公式 并令β＝βτ，据此，分别得到所述岩体细观颗粒二维粘结拉
4.根据权利要求1所述的岩体二维细观时效破裂幂函数型模型的构建方法，其特征在于：所述岩体细观颗粒粘结发生破裂后，岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达，用于描述岩体细观颗粒粘结时效破裂后细观颗粒的应力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：步骤301：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个二维线性接触端a、二维线性接触端b的中心坐标，在二维情况下，通过公式(19)计算接触点a端，接触点b端的中心距离：其中：d为二维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离， 为二维线性接触端a的坐标， 为二维线性接触端b的坐标；
步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)计算，如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标及距离计算；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算；确定每个接触端的单位向量：其中：ni为接触的单位矢量， 为接触端b的方向向量， 为接触端a的方向向量，nwall为约束墙的方向向量；
步骤303：岩体细观颗粒粘结破裂后，每一个二维线性接触点的接触重叠量U，通过步骤
301计算颗粒间距d以及二维线性接触两端的颗粒半径Ra、Rb，再利用公式(21)来确定；通过设定颗粒二维线性接触参考距离gr，并结合公式(22)，确定颗粒二维线性接触的距离gs：gs＝|U|-gr (22)
步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体或是墙体的刚度ka，kb等效代替为接触点的刚度，按公式(23)计算：其中：Kn、Ks为岩体细观颗粒接触点等效的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度；
步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)来计算；其中epqz为Ricci指标置换符号，按照公式(26)来计算：其中：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端两个接触端；
步骤306：对于岩体细观颗粒线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增量和切向位移增量：R＝min(Ra,Rb) (27)
其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k平为岩体细观颗粒系统平移刚度，k转为岩体细观颗粒系统旋转刚度；ΔUp1为岩体细观颗粒二维线性接触的总位移增量，Δδn、 物理意义相同，均表示岩体细观颗粒二维线性接触的法向位移增量，Δδs、 物理意义相同，均表示岩体细观颗粒二维线性接触的切向位移增量，Vp1与Vq1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1,q1为张量指标变换符号；
步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向和切向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的更新是采用前一步的切向位移增量与更新因子α的乘积获得；
其中：(gs)0为模型计算初始时刻的表面接触距离，gs为岩体细观颗粒接触的距离，α为位移更新因子；
步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加Ml＝1和绝对矢量累加Ml＝0模式，通过公式(33)、(34)计算获得；岩体细观颗粒二维线性接触的切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(35)计算获得：其中： 分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向线性力，kn、ks分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδn、Δδs分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量， 分别为岩体细观颗粒二维线性接触的初始法向力增量值、切向力增量值， 为岩体细观颗粒未滑动时的静摩擦力， 为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与 乘积得到；
步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式Md＝{0,2}和无拉伸模式Md＝{1,3}两种，通过公式(36)计算，其中mc为等效颗粒质量，按公式(37)计算，岩体细观颗粒线性接触的切向阻尼采用全剪切模式Md＝{0,1}和滑-剪模式Md＝{2,3}，按照公式(38)来计算，其中： 分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼力，βn为岩体细观颗粒线性接触的法向阻尼系数，βs为岩体细观颗粒线性接触的切向阻尼系数，kn为岩体细观颗粒线性接触的法向线性刚度，ks为岩体细观颗粒线性接触的切向线性刚度，分别为岩体细观颗粒线性接触的法向速率和切向速率，F*为岩体细观颗粒线性接触的全法向阻尼力，表达式为 m c为岩体细观等效颗粒质量，m(1)为岩体细观颗粒接触端1的细观颗粒质量，m(2)为岩体细观颗粒接触端2的细观颗粒质量，Fd为岩体细观颗粒线性接触的总阻尼力。
岩体二维细观时效破裂幂函数型模型的构建方法
技术领域
[0001] 本发明涉及工程岩体细观时效破裂分析技术领域，具体地指岩体二维细观时效破裂幂函数型模型的构建方法。
背景技术
[0002] 深部岩体工程开挖后的失稳与破坏往往不是在开挖后立刻发生的，一般都存在着明显的变形破裂时效性和灾变(岩爆、大变形等)的滞后性，严重危害工程的施工安全与长期运营。目前，在细观方面的时效力学研究成果相对较少。《深埋大理岩破裂扩展时间效应的颗粒流模拟》一文对锦屏大理岩破裂的时间效应进行了试验和二维数值分析(岩石力学与工程学报，2011，Vol.30 No.10：1989-1996)；《锦屏大理岩蠕变损伤演化细观力学特征的数值模拟研究》一文应用二维蠕变细观力学模型对锦屏大理岩短期和长期强度特征进行了数值研究(岩土力学，2013，Vol.34 No.12：3601-3608)。这类模型是以指数型构建驱动应力和裂纹扩展速度间的关系，用来描述岩石细观层面上的时效破裂，适用于应力和裂纹扩展速度之间符合指数表达方式的这类岩体。另外，这类模型还存在如下不足之处：(1)颗粒间的剪切破裂准则是一条与平行粘结正应力平行的水平直线，也即这种剪切破裂准则与平行粘结正应力状态无关，只要平行粘结剪切应力大于或等于固定平行粘结剪切破裂强度，颗粒间即可发生剪切破裂，无法体现岩体中不同平行粘结正应力具有不同平行粘结剪切破裂强度的客观事实；(2)没有考虑粘结力矩的差异作用对接触破坏的影响，将粘结力矩的贡献度对不同岩性的影响均视为一致；(3)对于应力和裂纹扩展速度之间不合符指数表达方式的岩体，这类模型缺乏适应性。
发明内容
[0003] 本发明的目的在于针对上述缺陷，提出了一种岩体二维细观时效破裂幂函数型模型的构建方法，该模型结构包括考虑弯矩贡献因子的岩体细观颗粒粘结应力二维模式、考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式、考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则和考虑阻尼效应的细观颗粒线性接触二维模型。本发明适应于应力和裂纹扩展速度之间的关系符合幂函数型的这类岩体，对于平面状态下这类深部岩体工程的围岩长期稳定性预测、评价以及优化设计提供技术支持。
[0004] 本发明的目的是通过如下措施来达到的：岩体二维细观时效破裂幂函数型模型的构建方法，其特殊之处在于，包括如下步骤：
[0005] 步骤1：设定岩体内部细观颗粒粘结接触的几何参数量，包括颗粒粘结面积和颗粒粘结惯性矩，Ra、Rb分别为二维粘结接触a端的颗粒半径、b端的颗粒半径，为岩体细观颗粒粘结直径乘数或半径乘数，在二维情况下，粘结单位厚度为1时的颗粒粘结面积A和粘结惯性矩I分别通过公式(2)、公式(3)来确定：
[0006]
[0007]
[0008]
[0009] 其中：为岩体细观颗粒二维粘结半径，为岩体细观颗粒二维粘结直径乘数或半径乘数，A为岩体细观颗粒二维粘结面积，I为岩体细观颗粒二维粘结惯性矩；
[0010] 步骤201：利用岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt，通过幂函数形式计算岩体细观颗粒二维粘结时效衰减劣化的直径 由公式(4)来确定；
[0011]
[0012] 其中： 为考虑弯矩贡献因子的二维粘结法向应力， 为判断岩体细观颗粒二维粘结开始时效劣化衰减时的应力阀值， 为岩体细观颗粒二维粘结的拉伸强度， 为考虑弯矩贡献因子的二维粘结应力比，β1、β2分别为岩体细观颗粒粘结时效劣化的第一控制参数、第二控制参数， 为岩体细观颗粒二维粘结随时间劣化衰减的直径， 为岩体细观颗粒二维粘结未衰减时的直径，Δt为岩体细观颗粒二维粘结时效衰减劣化的时间增量；
[0013] 步骤202：根据步骤201中的公式(4)，设定细观颗粒二维粘结直径的时效衰减因子β，见公式(5)：
[0014]
[0015] 其中： A′、I′、分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数， A、I、为岩体细观颗粒二维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数；
[0016] 步骤203：根据步骤1中的公式(1)～公式(3)以及步骤202中的公式(5)，设定岩体细观颗粒二维粘结几何参数时效劣化衰减模式，在平面二维情况下，岩体细观颗粒二维粘结直径随着时间增加而不断劣化衰减，粘结单位厚度为1时的颗粒粘结的面积和惯性矩也随着时间增加而不断劣化衰减，分别见公式(6)、公式(7)；
[0017]
[0018]
[0019] 其中： A′、I′分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩，A、I为岩体细观颗粒二维粘结未衰减时的粘结面积、粘结惯性矩；
[0020] 步骤204：依次计算第j个至第k个的岩体细观颗粒包含时间效应的二维粘结法向弯矩增量，在二维情况下，由粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δt，通过公式(8)、公式(9)、公式(10)确定第i个岩体细观颗粒二维粘结相对转角 岩体细观颗粒二维粘结法向增量位移 以及岩体细观颗粒二维粘结切向增量位移再结合步骤203中的公式(7)和步骤202中的公式(5)，确定第i个岩体细观颗粒包含时间效应的二维粘结弯矩增量，具体见公式(11)；
[0021]
[0022]
[0023]
[0024]
[0025] 其中：ff、j、k、i是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的初始标号值，ff为中间标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的最末标号值，i为第一个至最后一个二维粘结颗粒标号值， 分别为第i个岩体细观颗粒二维粘结接触的a端和b端的绝对运动速度和角速度，nn、ns分别为岩体细观颗粒二维粘结接触面的法向单位向量和切向单位向量， 分别为岩体细观颗粒二维粘结法向位移增量和切向位移增量， 为岩体细观颗粒二维粘结法向刚度， 为岩体细观颗粒二维粘结弯矩增量。
[0026] 步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和公式(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒粘结未破裂的粘结接触并包含时间效应的二维粘结法向力、切向力和切向弯矩，通过公式(12)、公式(13)、公式(14)计算第i个岩体细观颗粒二维粘结接触的粘结法向力、切向力和切向弯矩，在二维情况下，通过公式(15)来确定岩体细观颗粒粘结法向弯矩，
[0027] 法向力：
[0028] 切向力：
[0029] 切向弯矩：
[0030] 法向弯矩：
[0031] 其中： 分别为第i个岩体细观颗粒包含时间效应的粘结法向力、粘结切向力、包含时间效应的粘结法向弯矩、粘结切向弯矩、粘结法向位移增量和粘结切向位移增量， 为岩体细观颗粒二维粘结切向刚度，+＝为加法自反运算符，-＝为减法自反运算符。
[0032] 步骤206：设置弯矩贡献因子 考虑弯矩对岩体细观颗粒粘结法向应力的贡献程度，根据岩体细观颗粒粘结二维正应力计算公式 和二维粘结剪应力计算公式 同时将这两个公式中A、I以及 用A′、I′及 替换，然后将步骤
203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含幂函数型时间效应和细观颗粒粘结考虑弯矩贡献因子的二维粘结正应力 计算公式和二维粘结剪应力 计算公式，分别见公式(16)和公式(17)，
[0033]
[0034]
[0035] 步骤207：将步骤206中包含时间效应的 代入公式(18)，确定考虑弯矩贡献因子且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则，并且依次计算第j个至第k个的岩体细观颗粒二维粘结应力，用于判断岩体细观颗粒粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒粘结应力中包含了幂函数型时间效应和考虑弯矩贡献因子，[0036]
[0037] 其中：fs、fn分别为岩体细观颗粒二维粘结的时效剪切破裂准则、时效拉伸破裂准则， 为第i个接触的含幂函数型时间效应的二维粘结剪应力， 为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的二维粘结正应力， 分别为岩体细观颗粒二维粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维粘结的粘聚力，为岩体细观颗粒二维粘结的内摩擦角；fs大于等于0表示岩体细观颗粒粘结剪切破裂，小于0表示岩体细观颗粒粘结未发生剪切破裂；fn大于等于0表示岩体细观颗粒粘结拉伸破裂，小于0表示岩体细观颗粒粘结未发生拉伸破裂；
[0038] 步骤208：当步骤207中的公式(18)中的fs或fn大于等于0时，表明岩体细观颗粒粘结发生破裂，此时岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达；当步骤207中的公式(18)中的fs和fn都小于0时，表明岩体细观颗粒粘结未破裂，继续循环步骤201至207，计算、更新、判断岩体细观颗粒接触的粘结状态，直至岩体不产生新的细观颗粒粘结破裂或者岩体细观颗粒粘结破裂加速发展而形成宏观破坏，循环终止。
[0039] 优选地，所述岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt的确定，是采用考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式，由每个时间步内的细观颗粒二维粘结首次衰减破裂所损耗的时间来确定，即通过第一个细观颗粒粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个细观颗粒粘结破裂所需要 的 计 算 循 环 次 数 来 估 算 初 始 时 间 步 长 增 量 Δ t ，见 公 式其中， 为第i个接触的岩体细观颗粒二维粘结直径乘数，nc为第一个岩体细观颗粒二维粘结破裂所需的循环计算的次数，βσ、βτ分别为岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子、二维粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩体细观颗粒粘结数，∞为无穷大。
[0040] 优选地，所述岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子βσ和岩体细观颗粒二维粘结剪切强度对应的时效劣化因子βτ的确定包括如下步骤，其中这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的细观颗粒二维粘结破裂标号：
[0041] 步骤211：在二维情况下，由岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δtc，通过公式 确定细观颗粒粘结接触的相对转角 通过公式 确定细观颗粒粘结法向增量位移
通过公式 确定颗粒粘结切向增量位移 通过公式
确定颗粒粘结接触的弯矩增量；
[0042] 步骤212：根据步骤211中的公式 通过公式确定细观颗粒接触的粘结法向力；根据步骤211中的公式
通过公式 确定细观颗粒接触的粘结切向
力；根据步骤211中的公式 和公式 通过公式
确定细观颗粒接触的粘结切向弯矩；通过公式 确定细观颗粒
接触的粘结法向弯矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；
[0043] 步骤213：在二维情况下，通过公式 确定细观颗粒接触的粘结正应力，通过公式 确定细观颗粒接触的粘结剪应力，将这两个公式中A、I以及 用A′、I′及 替换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得细观颗粒粘结包含幂函数型时间效应和考虑弯矩贡献因子的二维粘结正应力计算公式 和包含幂函数型时间效应的二维粘结剪应力计算公式
[0044] 步骤214：将 代入公式 并令β＝βσ；将代入公式 并令β＝βτ，据此，分别得到所述岩体细观颗粒二维粘结拉
伸强度对应的时效劣化因子 以及岩体细观颗粒
二维粘结剪切强度对应的时效劣化因子
[0045] 优选地，所述岩体细观颗粒粘结发生破裂后，岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达，用于描述岩体细观颗粒粘结时效破裂后细观颗粒的应力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：
[0046] 步骤301：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个二维线性接触端a、二维线性接触端b(颗粒与颗粒、颗粒与墙)的中心坐标，在二维情况下，通过公式(19)计算接触点a端，接触点b端的中心距离：
[0047]
[0048] 其中：d为二维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离， 为二维线性接触端a的坐标， 为二维线性接触端b的坐标。
[0049] 步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)计算，如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标及距离计算；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算；确定每个接触端的单位向量：
[0050]
[0051] 其中：ni为接触的单位矢量， 为接触端b的方向向量， 为接触端a的方向向量，nwall为约束墙的方向向量；
[0052] 步骤303：岩体细观颗粒粘结破裂后，每一个二维线性接触点的接触重叠量U，通过步骤301计算颗粒间距d以及二维线性接触两端的颗粒半径Ra、Rb，再利用公式(21)来确定；通过设定颗粒二维线性接触参考距离gr，并结合公式(22)，确定颗粒二维线性接触的距离gs：
[0053]
[0054] gs＝|U|-gr (22)
[0055] 步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体或是墙体的刚度ka，kb等效代替为接触点的刚度，按公式(23)计算：
[0056]
[0057] 其中：Kn、Ks为岩体细观颗粒接触点等效的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度。
[0058] 步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)来计算。其中epqz为Ricci指标置换符号，按照公式(26)来计算：
[0059]
[0060]
[0061]
[0062] 其中：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端两个接触端)。
[0063] 步骤306：对于岩体细观颗粒线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增量和切向位移增量：
[0064] R＝min(Ra，Rb) (27)
[0065]
[0066]
[0067] ΔUp1＝Vp1Δt (30)
[0068]
[0069]
[0070] 其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k平为岩体细观颗粒系统平移刚度，k转为岩体细观颗粒系统旋转刚度；ΔUp1为岩体细观颗粒二维线性接触的总位移增量，Δδn、 物理意义相同，均表示岩体细观颗粒二维线性接触的法向位移增量，Δδs、 物理意义相同，均表示岩体细观颗粒二维线性接触的切向位移增量，Vp1与Vq1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1，q1为张量指标变换符号。
[0071] 步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向和切向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的更新是采用前一步的切向位移增量与更新因子α的乘积获得。
[0072]
[0073] 其中：(gs)0为模型计算初始时刻的表面接触距离，gs为岩体细观颗粒接触的距离，α为位移更新因子。
[0074] 步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加(Ml＝1)和绝对矢量累加(Ml＝0)模式，通过公式(33)、(34)计算获得；岩体细观颗粒二维线性接触的切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(35)计算获得：
[0075]
[0076]
[0077]
[0078] 其中： 分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向线性力，kn、ks分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδn、Δδs分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量， 分别为岩体细观颗粒二维线性接触的初始法向力增量值、切向力增量值， 为岩体细观颗粒未滑动时的静摩擦力， 为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与 乘积得到。
[0079] 步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式Md＝{0，2}和无拉伸模式Md＝{1，3}两种，通过公式(36)计算，其中mc为等效颗粒质量，按公式(37)计算，岩体细观颗粒线性接触的切向阻尼采用全剪切模式Md＝{0，1}和滑-剪模式Md＝{2，3}，按照公式(38)来计算，
[0080]
[0081]
[0082]
[0083] 其中： 分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼力，βn为岩体细观颗粒线性接触的法向阻尼系数，βs为岩体细观颗粒线性接触的切向阻尼系数，kn为岩体细观颗粒线性接触的法向线性刚度，ks为岩体细观颗粒线性接触的切向线性刚度， 分别为岩体细观颗粒线性接触的法向速率和切向速率，F*为岩体细观颗粒线性接触的全法向阻尼力，表达式为 mc为岩体细观等效颗粒质量，m(1)为岩体细观颗粒接触端1的细观颗粒质量，m(2)为岩体细观颗粒接触端2的细观颗粒质量，Fd为岩体细观颗粒线性接触的总阻尼力。
[0084] 本发明所提出的一种岩体二维细观时效破裂幂函数型模型及构建方法，其有益效果和优势主要体现在：
[0085] (1)本发明中岩体细观颗粒粘结正应力二维计算公式中设置了弯矩贡献因子，不仅考虑了弯矩对细观颗粒粘结正应力的贡献程度，而且还考虑了弯矩对岩体长期强度的影响，适合描述平面应力或平面应变条件下岩体的细观力学破裂行为。
[0086] (2)本发明中构建了考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式，包括在岩体细观颗粒粘结时效劣化衰减时，设置了幂函数型与考虑弯矩贡献因子的粘结应力相关的细观颗粒粘结二维劣化衰减模式，设置了细观颗粒粘结随时间逐步劣化衰减的二维模式，设置了细观颗粒粘结的面积和截面惯性矩时效劣化衰减二维模式；同时按照这种时效劣化衰减模式估算岩体细观颗粒粘结破裂的初始时间步长。这种幂函数型构建模式适合描述平面状态下一类深部岩体的细观力学时效破裂机制和响应规律。
[0087] (3)本发明中在所构建的二维细观时效破裂幂函数型模型中，嵌入考虑弯矩贡献效应且带拉伸截止限的的摩尔-库伦细观颗粒粘结时效破裂准则。在岩体细观颗粒粘结时效破裂时，采用所嵌入的考虑弯矩贡献效应且带拉伸截止限的的摩尔-库伦时效破裂准则来判断；在该准则颗粒粘结应力中包含幂函数型时间效应且增加了弯矩贡献因子，不仅可以描述与颗粒粘结正应力相关时效剪切破裂强度的差异，还可以对时效拉伸破裂进行合理的表达，且考虑了弯矩对粘结时效破裂的影响，符合平面条件下一类岩体细观时效破裂模式。
[0088] (4)本发明中在所构建的二维细观时效破裂幂函数型模型中，嵌入考虑阻尼效应的细观颗粒线性接触二维模型结构，在岩体时效破裂后，通过指定二维接触参考距离设定岩体细观颗粒间接触距离，设置考虑岩体细观颗粒间受力变形的二维接触模式以及在岩体细观颗粒之间设置考虑二维滑动摩擦的作用模式，同时设置二维接触的阻尼模式，可合理描述平面条件下一类深部工程岩体在时效破裂后的颗粒运动和受力特征。
附图说明
[0089] 图1是本发明模型中细观颗粒与颗粒接触示意图。
[0090] 图2是本发明模型中细观颗粒与刚性墙接触示意图。
[0091] 图3是本发明模型中细观颗粒重叠状态示意图。
[0092] 图4是本发明模型中细观颗粒刚度计算示意图。
[0093] 图5是本发明模型中细观颗粒粘结线性切向力和切向位移示意图。
[0094] 图6是本发明模型中细观颗粒粘接触本构物理模型示意图。
[0095] 图7是本发明模型中细观颗粒线性粘结结构示意图。
[0096] 图8是本发明模型中考虑弯矩贡献效应且带拉伸截止限的的摩尔-库伦细观颗粒粘结时效破裂准则示意图。
[0097] 图9是本发明模型中细观颗粒粘结直径(或半径)时效劣化衰减示意图。
[0098] 图10是本发明模型中细观颗粒二维接触面的法向和切向单位向量示意图。
[0099] 图11是本发明模型的构建方法的流程示意图。
[0100] 图12是本发明模型基础模型试样
[0101] 图13是本发明模型蠕变位移与时间关系曲线图
[0102] 图中：1代表两颗粒的中心距离d，2代表岩体细观颗粒间的半接触距离，3代表岩体细观颗粒间的半参考距离gr，4代表岩体细观颗粒a的坐标，5代表岩体细观颗粒b的坐标，6岩体细观颗粒表面接触距离的中心坐标，7代表岩体细观颗粒表面接触距离gs，8代表岩体细观颗粒间的单位接触法向量，9代表岩体细观颗粒a的半径Ra，10代表岩体细观颗粒b的半径Rb，11代表岩体细观颗粒接触点的接触重叠量U，12代表b(岩体颗粒或是边界墙体)的刚度(法向、切向刚度统称)kb，13代表a(岩体颗粒或是边界墙体)的刚度(法向、切向刚度统称)ka，14代表岩体细观颗粒接触点的等效刚度，15代表总位移增量ΔUi，16代表初始法向力增量值，17代表初始接触力矢量和，18代表初始切向力 增量值，19代表所构建的二维时效破裂模型法向位移增量Δδn或 20代表所构建的二维时效破裂模型切向位移增量Δδs或 21代表岩体细观颗粒粘结的拉伸强度 22代表岩体细观颗粒粘结法向刚度 23代表岩体细观颗粒接触点的法向刚度Kn，24代表岩体细观颗粒粘结切向刚度 25代表岩体细观颗粒粘结剪切强度，25.1代表 为岩体细观颗粒粘结的粘聚力，25.2代表岩体细观颗粒粘结内摩擦角 26代表岩体细观颗粒接触点的切向刚度Ks，27代表岩体细观颗粒线性接触滑动摩擦系数，28代表岩体细观颗粒线性接触法向阻尼系数βn，29代表岩体细观颗粒线性接触切向阻尼系数βs，30代表为岩体细观颗粒粘结半径乘数 31代表岩体细观颗粒粘结直径 32代表考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦时效破裂准则，33代表第i个接触的包含时间效应的粘结剪应力 34代表第i个接触的包含时间效应且考虑弯矩贡献因子的粘结正应力 35代表岩体细观颗粒粘结时效衰减的半径 36代表岩体细观颗粒粘结时效衰减的直径 37代表岩体细观颗粒粘结未衰减时的直径 38代表岩体细观颗粒粘结未衰减时的半径 39代表岩体细观颗粒接触面的法向向量nn，40代表岩体细观颗粒接触面切向单位向量ns。
具体实施方式
[0103] 下面结合附图和具体构建步骤及实施实例，对本发明模型进行详细的阐述。实例说明仅是辅助用于本发明的理解，而不限定本发明的实际应用范围。在阅读了本发明以后，本领域的技术人员对本发明的各种等价形式的修改均属于本发明所申请的权利要求所限定的范围。
[0104] 注：说明书中的所有的标号前面写明了公式，如公式(1)，均为公式标号。
[0105] 如图1～图10所示，本发明所述的岩体二维细观时效破裂幂函数型模型，适应于二维颗粒离散元、二维颗粒不连续变形分析方法、二维颗粒流形元；幂函数型模型包括考虑弯矩贡献因子的岩体细观颗粒粘结应力二维模式、考虑弯矩贡献因子的岩体细观颗粒粘结时效劣化衰减的二维幂函数型模式、考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则以及考虑阻尼效应的细观颗粒线性接触二维模型。
[0106] 考虑弯矩贡献因子的岩体细观颗粒粘结应力二维模式是指岩体细观颗粒粘结正应力二维计算公式 中设置了弯矩贡献因子 考虑弯矩对岩体细观颗粒二维粘结正应力的贡献程度；在上述 公式中， 为第i个岩体细观颗粒二维粘结正应力，分别为第i个接触的岩体细观颗粒二维粘结法向力、切向弯矩， 为岩体细观颗粒二维粘结半径，为弯矩贡献因子， I为岩体细观颗粒二维粘结的惯性矩，A为岩体细观颗粒二维粘结面积，i依次为第一个至最后一个岩体细观颗粒粘结数。第i个接触的岩体细观颗粒二维粘结法向力 粘结切向弯矩 的计算方法为：法向力
式中， 为岩体细观颗粒二维粘结法向位移增量， 为岩体细观颗粒二维粘结法向刚度，切向弯矩 式中， 为岩体细观颗粒二维粘结切向相对转角增量，+＝为加
法自反运算符，-＝为减法自反运算符，法向弯矩
[0107] 考虑弯矩贡献因子的岩体细观颗粒粘结时效劣化衰减的二维幂函数型模式包括在岩体细观颗粒二维粘结时效劣化衰减时，设置了与考虑弯矩贡献因子的粘结应力相关的幂函数型模式，这种幂函数型模式中的细观颗粒二维粘结直径随时间逐步劣化衰减，见粘结直径公式，
[0108]
[0109] 式中， 为岩体细观颗粒二维粘结随时间劣化衰减的直径， 为岩体细观颗粒二维粘结未衰减时的直径， 为考虑弯矩贡献因子的二维粘结法向应力， 为判断岩体细观颗粒二维粘结开始时效劣化衰减时的应力阀值， 为岩体细观颗粒二维粘结的拉伸强度，为考虑弯矩贡献因子的二维粘结应力比，β1、β2为控制岩体细观颗粒二维粘结时效劣化的两个指定指数，Δt为岩体细观颗粒二维粘结时效衰减劣化的时间增量；设置了岩体细观颗粒粘结面积和面惯性矩时效劣化衰减二维模式，分别见粘结单位厚度为1时的岩体细观颗粒随时间劣化衰减的粘结面积计算公式 粘结单位厚度为1时的惯性矩I′计算公式 其中，β为岩体细观颗粒二维粘结直径的时效衰减因子，其计算公式见
[0110]
[0111] 其中， A′、I′、分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数， A、I、为岩体细观颗粒粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数；同时按照这种幂函数型时效劣化衰减模式估算岩体细观颗粒粘结破裂的初始时间步长Δt，见公式[0112]
[0113] 其中， 为第i个接触的岩体细观颗粒二维粘结直径乘数，nc为第一个岩体细观颗粒二维粘结破裂所需的循环计算的次数，βσ、βτ分别为岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子、二维粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩体细观颗粒粘结数，∞为无穷大。
[0114] 岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子βσ和二维粘结剪切强度对应的时效劣化因子βτ的计算公式分别为
[0115]
[0116]
[0117] 其中， 分别为第i个颗粒接触的粘结法向力、切向力和粘结切向弯矩， 为岩体细观颗粒二维粘结拉伸强度，为岩体细观颗粒二维粘结的粘聚力，为岩体细观颗粒二维粘结的内摩擦角。
[0118] 考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则是指在岩体细观颗粒二维粘结时效破裂时，采用所嵌入的考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦时效破裂准则来判断，见公式
[0119] 其中，fs、fn分别为岩体细观颗粒二维粘结的时效剪切破裂准则、时效拉伸破裂准则，分别为岩体细观颗粒二维粘结拉伸强度、抗剪强度， 分别为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒二维粘结正应力、剪应力。
[0120] 第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的岩体细观颗粒二维粘结正应力的计算公式为
[0121]
[0122] 第i个接触的含幂函数型时间效应的岩体细观颗粒二维粘结剪应力的计算公式为[0123]
[0124] 在该准则的二维粘结应力中包含了幂函数型时间效应，见岩体细观颗粒二维粘结直径的时效衰减因子β计算公式
[0125] β1、β2分别为控制岩体细观颗粒粘结时效劣化的第一控制参数、第二控制参数；fs大于等于0表示岩体细观颗粒二维粘结剪切破裂，小于0表示岩体细观颗粒二维粘结未发生剪切破裂；fn大于等于0表示岩体细观颗粒二维粘结拉伸破裂，小于0表示岩体细观颗粒二维粘结未发生拉伸破裂。
[0126] 考虑阻尼效应的细观颗粒线性接触二维模型是指在岩体细观颗粒粘结时效破裂后，通过给定的二维线性接触参考距离gr设置了细观颗粒二维线性接触距离gs，见岩体细观颗粒二维线性接触距计算公式
[0127]
[0128] 其中， 为岩体内部颗粒与颗粒二维线性接触端a的坐标， 为岩体内部颗粒与颗粒二维线性接触端b的坐标，Ra、Rb分别为岩体细观二维线性接触端a的颗粒半径和二维线性接触端b的颗粒半径；设置了考虑岩体细观颗粒间变形的二维线性接触模式，在岩体细观颗粒之间设置了考虑二维滑动摩擦线力的作用模式，岩体细观颗粒间受力变形的二维线性接触法向线性力计算公式 取Ml＝1为相对矢量累加模式，取Ml＝0为绝对矢量累加模式，岩体细观颗粒间受力变形的二维线性接触切向线性力计算公式为 其中， 分别为岩体细观颗粒间受力变形的二
维线性接触法向线性力、切向线性力，kn、ks分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδn、Δδs分别为法向位移增量、切向位移增量， 分别为初始法向力增量值和切向力增量值， 为颗粒未滑动时的静摩擦力，
为岩体细观颗粒滑动摩擦力，通过摩擦系数μ与 乘积得到；
[0129] 同时设置二维线性接触的阻尼模式，其中法向阻尼采用全法向模式Md＝{0，2}和无拉伸模式Md＝{1，3}两种，通过公式 计算，其中mc为等效颗粒质量，按公式 计算，切向阻尼采用全剪切模式Md＝{0，1}和滑-
剪模式Md＝{2，3}，按照公式 来计算，其中： 分别为法
向阻尼力、切向阻尼力，βn为法向阻尼系数、βs为切向阻尼系数，kn为法向线性刚度、ks为切向线性刚度， 为法向速率、切向速率，mc为等效颗粒质量。F*为岩体细观颗粒线性接触的全法向阻尼力，表达式为 m(1)为二维线性接触端1的颗粒质量，m(2)
为二维线性接触端2的颗粒质量。
[0130] 如图11所示，本发明所述的岩体二维细观时效破裂幂函数型模型的构建方法，包括如下步骤：
[0131] 步骤1：设定岩体内部细观颗粒粘结接触的几何参数量，包括颗粒粘结面积和颗粒粘结惯性矩，Ra、Rb分别为二维粘结接触a端的颗粒半径、b端的颗粒半径，为岩体细观颗粒粘结直径乘数或半径乘数，在二维情况下，粘结单位厚度为1时的颗粒粘结面积A和粘结惯性矩I分别通过公式(2)、公式(3)来确定：
[0132]
[0133]
[0134]
[0135] 其中：为岩体细观颗粒二维粘结半径，为岩体细观颗粒二维粘结直径乘数或半径乘数，A为岩体细观颗粒二维粘结面积，I为岩体细观颗粒二维粘结惯性矩；
[0136] 步骤201：利用岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt，通过幂函数形式计算岩体细观颗粒二维粘结时效衰减劣化的直径 由公式(4)来确定；
[0137]
[0138] 其中： 为考虑弯矩贡献因子的二维粘结法向应力， 为判断岩体细观颗粒二维粘结开始时效劣化衰减时的应力阀值， 为岩体细观颗粒二维粘结的拉伸强度， 为考虑弯矩贡献因子的二维粘结应力比，β1、β2分别为岩体细观颗粒粘结时效劣化的第一控制参数、第二控制参数， 为岩体细观颗粒二维粘结随时间劣化衰减的直径， 为岩体细观颗粒二维粘结未衰减时的直径，Δt为岩体细观颗粒二维粘结时效衰减劣化的时间增量；
[0139] 步骤202：根据步骤201中的公式(4)，设定细观颗粒二维粘结直径的时效衰减因子β，见公式(5)：
[0140]
[0141] 其中： A′、I′、分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数， A、I、为岩体细观颗粒二维粘结未衰减时的粘结直径、粘结半径、粘结面积、粘结惯性矩、粘结直径乘数；
[0142] 步骤203：根据步骤1中的公式(1)～公式(3)以及步骤202中的公式(5)，设定岩体细观颗粒二维粘结几何参数时效劣化衰减模式，在平面二维情况下，岩体细观颗粒二维粘结直径随着时间增加而不断劣化衰减，粘结单位厚度为1时的颗粒粘结的面积和惯性矩也随着时间增加而不断劣化衰减，分别见公式(6)、公式(7)；
[0143]
[0144]
[0145] 其中： A′、I′分别表示为岩体细观颗粒二维粘结随时间劣化衰减的粘结半径、粘结面积、粘结惯性矩，A、I为岩体细观颗粒二维粘结未衰减时的粘结面积、粘结惯性矩；
[0146] 步骤204：依次计算第j个至第k个的岩体细观颗粒包含时间效应的二维粘结法向弯矩增量，在二维情况下，由粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δt，通过公式(8)、公式(9)、公式(10)确定第i个岩体细观颗粒二维粘结相对转角 岩体细观颗粒二维粘结法向增量位移 以及岩体细观颗粒二维粘结切向增量位移再结合步骤203中的公式(7)和步骤202中的公式(5)，确定第i个岩体细观颗粒包含时间效应的二维粘结弯矩增量，具体见公式(11)；
[0147]
[0148]
[0149]
[0150]
[0151] 其中：ff、j、k、i是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的初始标号值，ff为中间标号值，k为每次循环计算中，包含时间效应的岩体细观颗粒二维粘结衰减后未破裂的最末标号值，i为第一个至最后一个二维粘结颗粒标号值， 分别为第i个岩体细观颗粒二维粘结接触的a端和b端的绝对运动速度和角速度，nn、ns分别为岩体细观颗粒二维粘结接触面的法向单位向量和切向单位向量， 分别为岩体细观颗粒二维粘结法向位移增量和切向位移增量， 为岩体细观颗粒二维粘结法向刚度， 为岩体细观颗粒二维粘结弯矩增量。
[0152] 其中，岩体细观颗粒二维粘结时效衰减劣化的初始时间步长增量Δt的确定，是采用考虑弯矩贡献因子的细观颗粒粘结时效劣化衰减的二维幂函数型模式，由每个时间步内的细观颗粒二维粘结首次衰减破裂所损耗的时间来确定，即通过第一个细观颗粒粘结按幂函数型模式进行衰减破裂所历时的时间除以直至第一个细观颗粒粘结破裂所需要的计算循环次数来估算初始时间步长增量Δt，见公式
[0153] 其中， 为第i个接触的岩体细观颗粒二维粘结直径乘数，nc为第一个岩体细观颗粒二维粘结破裂所需的循环计算的次数，βσ、βτ分别为岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子、二维粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩体细观颗粒粘结数，∞为无穷大。
[0154] 岩体细观颗粒二维粘结拉伸强度对应的时效劣化因子βσ和岩体细观颗粒二维粘结剪切强度对应的时效劣化因子βτ的确定包括如下步骤，其中这些步骤中包含的公式下标1代表第一个按幂函数型模式进行时效衰减劣化的细观颗粒二维粘结破裂标号：
[0155] 步骤211：在二维情况下，由岩体细观颗粒粘结两端颗粒的速度、角速度和给定的循环计算时间步长增量Δtc，通过公式 确定细观颗粒粘结接触的相对转角 通过公式 确定细观颗粒粘结法向增量位移
通过公式 确定颗粒粘结切向增量位移 通过公式
确定颗粒粘结接触的弯矩增量；
[0156] 步骤212：根据步骤211中的公式 通过公式确定细观颗粒接触的粘结法向力；根据步骤211中的公式
通过公式 确定细观颗粒接触的粘结切向
力；根据步骤211中的公式 和公式 通过公式
确定细观颗粒接触的粘结切向弯矩；通过公式 确定细观颗粒
接触的粘结法向弯矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；
[0157] 步骤213：在二维情况下，通过公式 确定细观颗粒接触的粘结正应力，通过公式 确定细观颗粒接触的粘结剪应力，将这两个公式中A、I以及 用A′、I′及 替换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得细观颗粒粘结包含幂函数型时间效应和考虑弯矩贡献因子的二维粘结正应力计算公式 和包含幂函数型时间效应的二维粘结剪应力计算公式
[0158] 步骤214：将 代入公式 并令β＝βσ；将代入公式 并令β＝βσ，据此，分别得到所述岩体细观颗粒二维粘结拉
伸强度对应的时效劣化因子 以及岩体细观颗粒
二维粘结剪切强度对应的时效劣化因子
[0159] 步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和公式(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒粘结未破裂的粘结接触并包含时间效应的二维粘结法向力、切向力和切向弯矩，通过公式(12)、公式(13)、公式(14)计算第i个岩体细观颗粒二维粘结接触的粘结法向力、切向力和切向弯矩，在二维情况下，通过公式(15)来确定岩体细观颗粒粘结法向弯矩，
[0160] 法向力：
[0161] 切向力：
[0162] 切向弯矩：
[0163] 法向弯矩：
[0164] 其中： 分别为第i个岩体细观颗粒包含时间效应的粘结法向力、粘结切向力、包含时间效应的粘结法向弯矩、粘结切向弯矩、粘结法向位移增量和粘结切向位移增量， 为岩体细观颗粒二维粘结切向刚度，+＝为加法自反运算符，-＝为减法自反运算符。
[0165] 步骤206：设置弯矩贡献因子 考虑弯矩对岩体细观颗粒粘结法向应力的贡献程度，根据岩体细观颗粒粘结二维正应力计算公式 和二维粘结剪应力计算公式 同时将这两个公式中A、I以及 用A′、I′及 替换，然后将步骤
203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含幂函数型时间效应和细观颗粒粘结考虑弯矩贡献因子的二维粘结正应力 计算公式和二维粘结剪应力 计算公式，分别见公式(16)和公式(17)，
[0166]
[0167]
[0168] 步骤207：将步骤206中包含时间效应的 代入公式(18)，确定考虑弯矩贡献因子且带拉伸截止限的摩尔-库伦细观颗粒粘结时效破裂准则，并且依次计算第j个至第k个的岩体细观颗粒二维粘结应力，用于判断岩体细观颗粒粘结是否破裂以及破裂模式，在该准则的岩体细观颗粒粘结应力中包含了幂函数型时间效应和考虑弯矩贡献因子，[0169]
[0170] 其中：fs、fn分别为岩体细观颗粒二维粘结的时效剪切破裂准则、时效拉伸破裂准则， 为第i个接触的含幂函数型时间效应的二维粘结剪应力， 为第i个接触的含幂函数型时间效应且考虑弯矩贡献因子的二维粘结正应力， 分别为岩体细观颗粒二维粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维粘结的粘聚力，为岩体细观颗粒二维粘结的内摩擦角；fs大于等于0表示岩体细观颗粒粘结剪切破裂，小于0表示岩体细观颗粒粘结未发生剪切破裂；fn大于等于0表示岩体细观颗粒粘结拉伸破裂，小于0表示岩体细观颗粒粘结未发生拉伸破裂；
[0171] 步骤208：当步骤207中的公式(18)中的fs或fn大于等于0时，表明岩体细观颗粒粘结发生破裂，此时岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达；当步骤207中的公式(18)中的fs和fn都小于0时，表明岩体细观颗粒粘结未破裂，继续循环步骤201至207，计算、更新、判断岩体细观颗粒接触的粘结状态，直至岩体不产生新的细观颗粒粘结破裂或者岩体细观颗粒粘结破裂加速发展而形成宏观破坏，循环终止。
[0172] 所述岩体细观颗粒粘结发生破裂后，岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来表达，用于描述岩体细观颗粒粘结时效破裂后细观颗粒的应力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：
[0173] 步骤301：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个二维线性接触端a、二维线性接触端b(颗粒与颗粒、颗粒与墙)的中心坐标，在二维情况下，通过公式(19)计算接触点a端，接触点b端的中心距离：
[0174]
[0175] 其中：d为二维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离， 为二维线性接触端a的坐标， 为二维线性接触端b的坐标。
[0176] 步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)计算，如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标及距离计算；如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算；确定每个接触端的单位向量：
[0177]
[0178] 其中：ni为接触的单位矢量， 为接触端b的方向向量， 为接触端a的方向向量，nwall为约束墙的方向向量；
[0179] 步骤303：岩体细观颗粒粘结破裂后，每一个二维线性接触点的接触重叠量U，通过步骤301计算颗粒间距d以及二维线性接触两端的颗粒半径Ra、Rb，再利用公式(21)来确定；通过设定颗粒二维线性接触参考距离gr，并结合公式(22)，确定颗粒二维线性接触的距离gs：
[0180]
[0181] gs＝|U|-gr (22)
[0182] 步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体或是墙体的刚度ka，kb等效代替为接触点的刚度，按公式(23)计算：
[0183]
[0184] 其中：Kn、Ks为岩体细观颗粒接触点等效的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触a端的法向刚度和切向刚度， 为颗粒与颗粒或者颗粒与墙的接触b端的法向刚度和切向刚度。
[0185] 步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)来计算。其中epqz为Ricci指标置换符号，按照公式(26)来计算：
[0186]
[0187]
[0188]
[0189] 其中：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端两个接触端)。
[0190] 步骤306：对于岩体细观颗粒线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增量和切向位移增量：
[0191] R＝min(Ra，Rb) (27)
[0192]
[0193]
[0194] ΔUp1＝Vp1Δt (30)
[0195]
[0196]
[0197] 其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的转动惯量；k平为岩体细观颗粒系统平移刚度，k转为岩体细观颗粒系统旋转刚度；ΔUp1为岩体细观颗粒二维线性接触的总位移增量，Δδn、 物理意义相同，均表示岩体细观颗粒二维线性接触的法向位移增量，Δδs、 物理意义相同，均表示岩体细观颗粒二维线性接触的切向位移增量，Vp1与Vq1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1，q1为张量指标变换符号。
[0198] 步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向和切向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的更新是采用前一步的切向位移增量与更新因子α的乘积获得。
[0199]
[0200] 其中：(gs)0为模型计算初始时刻的表面接触距离，gs为岩体细观颗粒接触的距离，α为位移更新因子。
[0201] 步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加(Ml＝1)和绝对矢量累加(Ml＝0)模式，通过公式(33)、(34)计算获得；岩体细观颗粒二维线性接触的切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(35)计算获得：
[0202]
[0203]
[0204]
[0205] 其中： 分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向线性力，kn、ks分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδn、Δδs分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量， 分别为岩体细观颗粒二维线性接触的初始法向力增量值、切向力增量值， 为岩体细观颗粒未滑动时的静摩擦力， 为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与 乘积得到。
[0206] 步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式Md＝{0，2}和无拉伸模式Md＝{1，3}两种，通过公式(36)计算，其中mc为等效颗粒质量，按公式(37)计算，岩体细观颗粒线性接触的切向阻尼采用全剪切模式Md＝{0，1}和滑-剪模式Md＝{2，3}，按照公式(38)来计算，
[0207]
[0208]
[0209]
[0210] 其中： 分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼力，βn为岩体细观颗粒线性接触的法向阻尼系数，βs为岩体细观颗粒线性接触的切向阻尼系数，kn为岩体细观颗粒线性接触的法向线性刚度，ks为岩体细观颗粒线性接触的切向线性刚度， 分别为岩体细观颗粒线性接触的法向速率和切向速率，F*为岩体细观颗粒线性接触的全法向阻尼力，表达式为 mc为岩体细观等效颗粒质量，m(1)为岩体细观颗粒接触端1的细观颗粒质量，m(2)为岩体细观颗粒接触端2的细观颗粒质量，Fd为岩体细观颗粒线性接触的总阻尼力。
[0211] 下面以深部岩体为实例，结合附图详述本发明模型的数值实现的具体过程，请参阅实例附图说明中的图12至图13以及模型附图说明中的图1至图10，来理解本发明模型的数值实现步骤及效果：
[0212] 步骤1：采用C++编程语言，并结合fish语言，根据本发明的模型结构构建流程图(图11)，在数值平台上实现了岩体二维细观时效破裂幂函数型模型。
[0213] 步骤2：初步确定岩体时效破裂模型的细观参数
[0214] 粒径比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、法向阻尼系数Apfan(图6)、切向阻尼系数Apfas(图6)以及内摩擦角pb_phi(图8)等19个参数，参数具体值见表一。
[0215] 步骤3：生成岩体模型
[0216] 根据高斯分布或weibull分布确定模型的粘结拉伸强度和粘聚力分布，通过均匀随机函数分布法确定颗粒的粒径分布；通过各向同性应力调整法及自适应动态膨胀法，调整颗粒的位置，减少颗粒重叠量；通过悬浮颗粒删除法，删除孤立颗粒，提高模型样本的整体性，减少缺陷模型的生成。最后赋予模型材料粘结强度参数，生成具有真实岩体结构的岩体模型。岩体模型的直径为50mm、高度为100mm(图12)。
[0217] 步骤4：精确标定本发明中模型的细观物理力学参数
[0218] 通过室内单轴和三轴压缩试验得到的应力-应变曲线，确定岩体的宏观弹性模量峰值强度σp，以及泊松比 通过优化方法，使岩体单、三轴压缩模型的应力-应变曲线与室内试验的应力-应变以及宏观变形参数和强度参数吻合，获得本发明所构建模型的细观物理力学参数。
[0219] 步骤5：岩体时效力学参数标定
[0220] 对岩体进行一系列不同应力强度比条件下的时效力学试验，得到不同应力强度比条件下岩体变形随时间演化的曲线。通过参数拟合法，匹配实际岩体的时效变形过程，确定岩体细观颗粒粘结时效劣化的第一控制参数β1、第二控制参数β2。
[0221] 步骤6：岩体时效力学数值试验
[0222] 在荷载一定的条件下，分别设定不同的弯矩贡献因子，获得了岩体时效变形破坏的演化规律(图13)。
[0223] 本发明模型的参数名称及值如表一所示。
[0224] 表一：本发明模型的参数名称及值
[0225]
[0226]
[0227] 上述实施例中，公式的符号与图1～图10以及附图说明中的符号相互对应。
[0228] 其它未详细说明的部分均为现有技术，以上所有参数均可通过查阅手册或计算得到。本发明并不严格地局限于上述实施例。以上所述仅为本发明的特定实施例，并不用于限制本发明。凡在本发明的精神和原则之内所做的任何修改、等同替换及改进等，都在本发明的保护范围之内。
