考虑弯矩贡献因子的二维时效破裂模型的构建方法.pdf

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

1.一种考虑弯矩贡献因子的二维时效破裂模型的构建方法，其特征在于：包括如下步
骤：
步骤1：设定岩体细观颗粒平行粘结接触的几何参数量包括平行粘结面积和平行粘结
惯性矩，R^a、R^b分别为二维平行粘结接触a端的颗粒半径、b端的颗粒半径，为岩体细观颗粒
平行粘结直径或半径乘数，在二维情况下，平行粘结单位厚度为1时的平行粘结面积A和平
行粘结惯性矩I分别通过公式(2)、公式(3)来确定：
$\stackrel{\‾}{R}=\stackrel{\‾}{\mathrm{\λ}}\min (R^a,R^b)---\left(1\right)$
$A=2\stackrel{\‾}{R}---\left(2\right)$
$I=\frac{2}{3}{\stackrel{\‾}{R}}^3---\left(3\right)$
其中：为岩体细观颗粒二维平行粘结半径，为岩体细观颗粒二维平行粘结直径乘
数，A为岩体细观颗粒二维平行粘结面积圆，I为岩体细观颗粒二维平行粘结惯性矩；
步骤201：利用岩体细观颗粒二维平行粘结时效衰减劣化的初始时间步长增量Δt，通
过指数型函数计算岩体细观颗粒二维平行粘结时效衰减劣化的直径，公式(4)来确定：
${\stackrel{\‾}{D}}^{\′}=\stackrel{\‾}{D}-\left\{\begin{array}{c}0,\stackrel{\‾}{\mathrm{\σ}}\<{\stackrel{\‾}{\mathrm{\σ}}}_{aa}\\ {\mathrm{\β}}_1{e}^{{\mathrm{\β}}_2(\stackrel{\‾}{\mathrm{\σ}}/{\stackrel{\‾}{\mathrm{\σ}}}_c)}\mathrm{\Δ}t,{\stackrel{\‾}{\mathrm{\σ}}}_{aa}\≤\stackrel{\‾}{\mathrm{\σ}}\<{\stackrel{\‾}{\mathrm{\σ}}}_c\end{array}\right.---\left(4\right)$
其中：为考虑弯矩贡献因子的二维平行粘结法向应力，为判断岩体细观颗粒二维
平行粘结开始时效劣化衰减时的应力阀值，为岩体细观颗粒二维平行粘结的拉伸强度，
为考虑弯矩贡献因子的二维平行粘结应力比，β₁、β₂分别为岩体细观颗粒平行粘结时
效劣化的第一控制参数、第二控制参数，为岩体细观颗粒二维平行粘结随时间劣化衰减
的直径，为岩体细观颗粒二维平行粘结未衰减时的直径，e为自然常数，Δt为岩体细观颗
粒二维平行粘结时效衰减劣化的时间增量；
步骤202：根据步骤201中的公式(4)，设定岩体细观颗粒二维平行粒粘结直径的指数型
时效衰减因子β，见公式(5)：
$\mathrm{\β}={\stackrel{\‾}{D}}^{\′}/\stackrel{\‾}{D}={\stackrel{\‾}{R}}^{\′}/\stackrel{\‾}{R}={\stackrel{\‾}{\mathrm{\λ}}}^{\′}/\stackrel{\‾}{\mathrm{\λ}}=1-\left\{\begin{array}{c}0,\stackrel{\‾}{\mathrm{\σ}}\<{\stackrel{\‾}{\mathrm{\σ}}}_{aa}\\ {\mathrm{\β}}_1{e}^{{\mathrm{\β}}_2(\stackrel{\‾}{\mathrm{\σ}}/{\stackrel{\‾}{\mathrm{\σ}}}_c)}\frac{\mathrm{\Δ}t}{\stackrel{\‾}{D}},{\stackrel{\‾}{\mathrm{\σ}}}_{aa}\≤\stackrel{\‾}{\mathrm{\σ}}\<{\stackrel{\‾}{\mathrm{\σ}}}_c\end{array}\right.---\left(5\right)$
其中：A'、I'、分别表示为岩体细观颗粒二维平行粘结随时间劣化衰减的平行
粘结直径、平行粘结半径、平行粘结面积、平行粘结惯性矩、平行粘结直径乘数，
A、I、为岩体细观颗粒平行粘结未衰减时的平行粘结直径、平行粘结半
径、平行粘结面积、平行粘结惯性矩、平行粘结直径乘数；
步骤203：根据步骤1中的公式(1)～公式(3)以及步骤202中的公式(5)，设定岩体细观
颗粒二维平行粘结几何参数时效劣化衰减模式；岩体细观颗粒二维平行粘结直径随着时间
增加而不断劣化衰减，二维平行粘结单位厚度为1时的平行粘结面积和平行粘结惯性矩也
随着时间增加而不断劣化衰减，分别见公式(6)、公式(7)：
$A^{\′}=2{\stackrel{\‾}{R}}^{\′}=\mathrm{\β}A---\left(6\right)$
$I^{\′}=\frac{2}{3}{\stackrel{\‾}{R}}^{\′3}={\mathrm{\β}}^{3}I---\left(7\right)$
其中：A'、I'分别表示为岩体细观颗粒二维平行粘结随时间劣化衰减的平行粘结半
径、平行粘结面积、平行粘结惯性矩，A、I为岩体细观颗粒平行粘结未衰减时的平行粘结面
积、平行粘结惯性矩；
步骤204：依次计算第j个至第k个的岩体细观颗粒包含时间效应的二维平行粘结法向
弯矩增量；在二维情况下，由平行粘结两端颗粒的速度、角速度和给定的循环计算时间步长
增量Δt_c，通过公式(8)、公式(9)、公式(10)确定第i个岩体细观颗粒二维平行粘结相对转
角二维平行粘结法向增量位移以及二维平行粘结切向增量位移再结合
步骤203中的公式(7)和步骤202中的公式(5)，确定第i个岩体细观颗粒包含时间效应的二
维平行粘结弯矩增量，具体见公式(11)：
${\left({\mathrm{\ΔU}}_i^n\right)}_{ff}=({\stackrel{\·}{X}}_i^b-{\stackrel{\·}{X}}_i^a)n_n{\mathrm{\Δt}}_c---\left(8\right)$
${\left({\mathrm{\ΔU}}_i^s\right)}_{ff}=({\stackrel{\·}{X}}_i^b-{\stackrel{\·}{X}}_i^a)n_s{\mathrm{\Δt}}_c---\left(9\right)$
${\left({\mathrm{\Δ\θ}}_i^n\right)}_{ff}=({\mathrm{\ω}}_i^b-{\mathrm{\ω}}_i^a){\mathrm{\Δt}}_c---\left(10\right)$
${\left(\mathrm{\Δ}{{\stackrel{\‾}{M}}_i}^s\right)}_{ff}={\stackrel{\‾}{k}}_n{\mathrm{\β}}^{3}I{\left({\mathrm{\Δ\θ}}_i^n\right)}_{ff}---\left(11\right)$
其中，ff、j、k是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩体细
观颗粒二维平行粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循环计
算中，包含时间效应的岩体细观颗粒二维平行粘结衰减后未破裂的最末标号值，i为第一个
至最后一个平行粘结颗粒标号值，分别为第i个岩体细观颗粒二维平
行粘结接触的a端和b端的绝对运动速度和角速度，n_n和n_s为岩体细观颗粒二维平行粘结接
触面的法向单位向量和切向单位向量，和分别为岩体细观颗粒二维平行粘结法
向位移增量和切向位移增量，为岩体细观颗粒平行粘结法向刚度，为岩体细观颗粒
二维平行粘结弯矩增；
步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和公式
(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒平行粘结未破裂
并包含时间效应的二维平行粘结法向力、切向力和切向弯矩；通过公式(12)、公式(13)、公
式(14)计算第i个岩体细观颗粒接触的二维平行粘结法向力、切向力和切向弯矩；在二维情
况下，通过公式(15)来确定岩体细观颗粒平行粘结法向弯矩：
法向力：
切向力：
切向弯矩：
法向弯矩：
其中：分别为第i个岩体细观颗粒包含时间效应的
平行粘结法向力、平行粘结切向力、包含时间效应的平行粘结法向弯矩、平行粘结切向弯
矩、平行粘结法向位移增量和平行粘结切向位移增量，为岩体细观颗粒二维平行粘结切
向刚度，+＝为加法自反运算符，-＝为减法自反运算符；
步骤206：设置弯矩贡献因子考虑弯矩对岩体细观颗粒平行粘结法向正应
力的贡献程度，根据二维平行粘结正应力计算公式和二维平行粘结剪应
力计算公式同时将这两个公式中A、I以及用A'、I'及替换，然后将步骤203中
的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型时间效应和弯矩贡献
因子的二维平行粘结正应力和二维平行粘结剪应力计算公式，分别见公式(16)和
公式(17)；
${{\left(\stackrel{\‾}{{\mathrm{\σ}}_i}\right)}^{\′}}_{ff}=\frac{1}{2\mathrm{\β}\stackrel{\‾}{R}}(-{\left({\stackrel{\‾}{F_i}}^n\right)}_{ff}+\frac{3\stackrel{\‾}{\mathrm{\β}}\left|{\left({\stackrel{\‾}{M_i}}^s\right)}_{ff}\right|}{\mathrm{\β}\stackrel{\‾}{R}})---\left(16\right)$
${{\left(\stackrel{\‾}{{\mathrm{\τ}}_i}\right)}^{\′}}_{ff}=\frac{\left|{\left({\stackrel{\‾}{F_i}}^s\right)}_{ff}\right|}{2\mathrm{\β}\stackrel{\‾}{R}}---\left(17\right)$
步骤207：将步骤206中包含时间效应的代入公式(18)，可确定考虑弯矩贡献
因子且带拉伸截止限的摩尔-库伦时效破裂准则，并且依次计算第j个至第k个的二维平行
粘结应力，用于判断岩体细观颗粒平行粘结是否破裂以及破裂模式；在该准则的岩体细观
颗粒二维平行粘结应力中包含了指数型时间效应和考虑弯矩贡献因子；
$\left\{\begin{array}{c}{f}^s={\left({\stackrel{\‾}{\mathrm{\τ}}}_i\right)}^{\′}-{\stackrel{\‾}{\mathrm{\τ}}}_c={\left({\stackrel{\‾}{\mathrm{\τ}}}_i\right)}^{\′}+{\left({\stackrel{\‾}{\mathrm{\σ}}}_i\right)}^{\′}\tan \stackrel{\‾}{\mathrm{\φ}}-\stackrel{\‾}{c}=\frac{\left|{\stackrel{\‾}{F_i}}^s\right|}{2\mathrm{\β}\stackrel{\‾}{R}}+\frac{1}{2\mathrm{\β}\stackrel{\‾}{R}}(-{\stackrel{\‾}{F_i}}^n+\frac{3\stackrel{\‾}{\mathrm{\β}}\left|{\stackrel{\‾}{M_i}}^s\right|}{\mathrm{\β}\stackrel{\‾}{R}})\tan \stackrel{\‾}{\mathrm{\φ}}-\stackrel{\‾}{c}\\ f^n={\left({\stackrel{\‾}{\mathrm{\σ}}}_i\right)}^{\′}-{\stackrel{\‾}{\mathrm{\σ}}}_c=\frac{1}{2\mathrm{\β}\stackrel{\‾}{R}}(-{\stackrel{\‾}{F_i}}^n+\frac{3\stackrel{\‾}{\mathrm{\β}}\left|{\stackrel{\‾}{M_i}}^s\right|}{\mathrm{\β}\stackrel{\‾}{R}})-{\stackrel{\‾}{\mathrm{\σ}}}_c\end{array}\right.---\left(18\right)$
其中：f^s、fⁿ分别为岩体细观颗粒二维平行粘结的时效剪切破裂准则、时效拉伸破裂准
则，为第i个接触的含指数型时间效应的二维平行粘结剪应力，为第i个接触的
含指数型时间效应且考虑弯矩贡献因子的二维平行粘结正应力，分别为岩体细观
颗粒二维平行粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维平行粘结的粘聚力，为
岩体细观颗粒二维平行粘结的内摩擦角；f^s大于等于0表示平行粘结剪切破裂，小于0表示
平行粘结未发生剪切破裂；fⁿ大于等于0表示平行粘结拉伸破裂，小于0表示平行粘结未发
生拉伸破裂；
步骤208：如果步骤207中的公式(18)中的f^s或fⁿ大于等于0，表明岩体细观颗粒的平行
粘结发生了破裂，此后岩体细观颗粒的运动模式采用考虑阻尼效应的二维线性接触模型来
表达；如果步骤207中的公式(18)中的f^s和fⁿ都小于0，表明平行粘结未破裂，继续循环步骤
201至207，计算、更新、判断岩体细观颗粒接触的平行粘结状态，直至岩体不产生新的平行
粘结破裂或者平行粘结破裂加速发展而形成宏观破坏，循环终止。
2.根据权利要求1所述的考虑弯矩贡献因子的二维时效破裂模型的构建方法，其特征在
于：所述岩体细观颗粒二维平行粘结时效衰减劣化的初始时间步长增量Δt的确定，是采用考
虑弯矩贡献因子的平行粘结时效劣化衰减的二维指数型模式，由每个时间步内的二维平行粘
结首次衰减破裂所损耗的时间来确定，也即通过第一个平行粘结按指数型模式进行衰减破裂
所历时的时间除以直至第一个平行粘结破裂所需要的计算循环次数来估算初始时间步长
增量Δt，见公式其中，
为第i个接触的岩体细观颗粒二维平行粘结直径乘数，n_c为第一个岩体细观颗粒二维平行
粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体细观颗粒二维平行粘结拉伸强度对应的
时效劣化因子、二维平行粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩
体细观颗粒平行粘结数，∞为无穷大。
3.根据权利要求2所述的考虑弯矩贡献因子的二维时效破裂模型的构建方法，其特征
在于：所述岩体细观颗粒二维平行粘结拉伸强度对应的时效劣化因子β_σ和岩体细观颗粒二
维平行粘结剪切强度对应的时效劣化因子β_τ的确定包括如下步骤，其中，这些步骤中包含
的公式下标1代表第一个按指数型模式进行时效衰减劣化的二维平行粘结破裂标号；
步骤211：在二维情况下，由岩体细观颗粒平行粘结两端颗粒的速度、角速度和给定的
循环计算时间步长增量Δt_c，通过公式确定平行粘结接触的相对转
角通过公式确定平行粘结法向增量位移通过公式
确定平行粘结切向增量位移通过公式
确定平行粘结接触的弯矩增量；
步骤212：根据步骤211中的公式通过公式
确定平行粘结法向力；根据步骤211中的公式通过公式
确定平行粘结切向力；根据步骤211中的公式和
公式通过公式确定平行粘结切向弯矩；通过公
式确定平行粘结法向弯矩，其中，+＝为加法自反运算符，-＝为减法自反运算符；
步骤213：在二维情况下，通过公式确定平行粘结正应力，通
过公式确定平行粘结剪应力，将这两个公式中A、I以及用A'、I'及替换，然
后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型时间
效应和弯矩贡献因子的二维平行粘结正应力计算公式和包含
指数型时间效应的二维平行粘结剪应力计算公式
步骤214：将代入公式并令β＝β_σ，将
代入公式并令β＝β_τ，据此，分别得到所述岩体细观颗粒二维平行
粘结拉伸强度对应的时效劣化因子以及岩体细
观颗粒二维平行粘结剪切强度对应的时效劣化因子
4.根据权利要求1所述的考虑弯矩贡献因子的二维时效破裂模型的构建方法，其特征
在于：所述岩体细观颗粒平行粘结发生破裂后，岩体细观颗粒的运动模式采用考虑阻尼效
应的二维线性接触模型来表达，用于描述岩体细观颗粒平行粘结时效破裂后细观颗粒的应
力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：
步骤301：通过Monte Carlo搜索算法，遍历寻找岩体细观颗粒每个二维线性接触端a、
二维线性接触端b，颗粒与颗粒或者颗粒与墙的中心坐标，在二维情况下，通过公式(19)计
算两者中心距离：
$d=\sqrt{{(x_i^b-x_i^a)}^2+{(y_i^b-y_i^a)}^2}---\left(19\right)$
其中：d为二维线性接触两端颗粒与颗粒或者颗粒与墙之间的中心距离，为二维
线性接触端a的坐标，为二维线性接触端b的坐标；
步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)计算，
如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标及距
离计算，如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触段的
单位向量：

其中：n_i为接触的单位矢量，为接触端b的方向向量，为接触端a的方向向量，n_wall
为约束墙的方向向量；
步骤303：岩体细观颗粒平行粘结破裂后，二维线性接触段的接触重叠量U，通过步骤
301计算颗粒间距d以及二维线性接触两端的颗粒半径R^a、R^b，再利用公式(21)来确定；通过
设定颗粒二维线性接触参考距离g_r，并结合公式(22)，确定颗粒二维线性接触的距离g_s：

g_s＝|U|-g_r (22)
步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体或是
墙体的刚度k^a，k^b等效代替为接触点的刚度，按公式(23)计算：
$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(23\right)$
其中：K_n、K_s为等效的法向刚度和切向刚度，为颗粒与颗粒或者颗粒与墙的接触
a端的法向和切向刚度，为颗粒与颗粒或者颗粒与墙的接触b端的法向和切向刚度；
步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)来计
算；其中e_pqz为Ricci指标置换符号，按照公式(26)来计算：
$V_p={\left({\stackrel{\·}{x}}_p^{\left(c\right)}\right)}_b-{\left({\stackrel{\·}{x}}_p^{\left(c\right)}\right)}_a=({\stackrel{\·}{x}}_p^{\left(b\right)}+e_{pqz}{\mathrm{\ω}}_q^{\left(b\right)}(x_z^{\left(c\right)}-x_z^{\left(b\right)}\left)\right)-({\stackrel{\·}{x}}_p^{\left(a\right)}+e_{pqz}{\mathrm{\ω}}_q^{\left(a\right)}(x_z^{\left(c\right)}-x_z^{\left(a\right)}\left)\right)---\left(24\right)$
$V_q={\left({\stackrel{\·}{x}}_q^{\left(c\right)}\right)}_b-{\left({\stackrel{\·}{x}}_q^{\left(c\right)}\right)}_a=({\stackrel{\·}{x}}_q^{\left(b\right)}+e_{qpz}{\mathrm{\ω}}_p^{\left(b\right)}(x_z^{\left(c\right)}-x_z^{\left(b\right)}\left)\right)-({\stackrel{\·}{x}}_q^{\left(a\right)}+e_{qpz}{\mathrm{\ω}}_p^{\left(a\right)}(x_z^{\left(c\right)}-x_z^{\left(a\right)}\left)\right)---\left(25\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端两个接触端；
步骤306：对于线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计最小
的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋于稳
定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增量和
切向位移增量：
R＝min(R_a,R_b) (27)
$J1=\frac{2}{5}{\mathrm{\πR}}^5---\left(28\right)$

ΔU_p1＝V_p1Δt (30)
${\mathrm{\Δ\δ}}_n={\mathrm{\ΔU}}_{p1}^{nl}=V_{q1}{n}_{q1}{n}_{p1}\mathrm{\Δ}t---\left(31\right)$
${\mathrm{\Δ\δ}}_s={\mathrm{\ΔU}}_{p1}^{sl}={\mathrm{\ΔU}}_{p1}-{\mathrm{\ΔU}}_{p1}^{nl}=V_{p1}\mathrm{\Δ}t-V_{q1}{n}_{q1}{n}_{p1}\mathrm{\Δ}t---\left(32\right)$
其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的转动
惯量；k^平为岩体细观颗粒系统平移刚度，k^转为岩体细观颗粒系统旋转刚度；ΔU_p1为岩体细观
颗粒二维线性接触的总位移增量，Δδ_n、物理意义相同，均表示岩体细观颗粒二维线性
接触的法向位移增量，Δδ_s、物理意义相同，均表示岩体细观颗粒二维线性接触的切向
位移增量，V_p1与V_q1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1,q1为张
量指标变换符号；
步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向和切
向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一步的
法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的更新
是采用前一步的切向位移增量与更新因子α的乘积获得；
$\mathrm{\α}=\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(33\right)$
其中：(g_s)₀为模型计算初始时刻的表面接触距离，g_s为岩体细观颗粒接触的距离，α为
位移更新因子；
步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加(M_l＝1)和绝对
矢量累加(M_l＝0)模式，通过公式(34)、(35)计算获得；岩体细观颗粒二维线性接触的切向
线性力采用岩体细观颗粒接触滑动来表示，通过公式(36)计算获得：
$F_n^l=\left\{\begin{array}{c}\left\{\begin{array}{c}{k}_n{g}_s,g_s\<0\\ 0,otherwlse\end{array}\right.,M_l=0\\ min({\left(F_n^l\right)}_o+k_n{\mathrm{\Δ\δ}}_n,0),M_l=1\end{array}\right.---\left(34\right)$
$F_s^{*}={\left(F_s^l\right)}_o-k_s{\mathrm{\Δ\δ}}_s---\left(35\right)$
$F_s^l=\left\{\begin{array}{c}{F}_s^{*},\left|\right|F_s^{*}\left|\right|\≤F_s^{\mathrm{\μ}}\\ F_s^{\mathrm{\μ}}(F_s^{*}/|\left|F_s^{*}\right|\left|\right),otherwise\end{array}\right.---\left(36\right)$
其中：分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向线性
力，k_n、k_s分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδ_n、Δδ_s
分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量，分别为岩
体细观颗粒二维线性接触的初始法向力增量值、切向力增量值，为岩体细观颗粒未滑动
时的静摩擦力，为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与乘积得到；
步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式M_d＝{0,2}和无拉伸模式
M_d＝{1,3}两种，通过公式(37)计算，其中m_c为等效颗粒质量，按公式(38)计算，岩体细观颗
粒线性接触的切向阻尼采用全剪切模式M_d＝{0,1}和滑-剪模式M_d＝{2,3}，按照公式(39)来
计算，
$F_n^d=\left\{\begin{array}{c}\left(2{\mathrm{\β}}_n\sqrt{m_c{k}_n}\right){\stackrel{\·}{\mathrm{\δ}}}_n,M_d=\{0,2\}\\ \min (F^{*},-F_n^l),M_d=\{1,3\}\end{array}\right.---\left(37\right)$

$F_s^d=\left\{\begin{array}{c}\left(2{\mathrm{\β}}_s\sqrt{m_c{k}_s}\right){\stackrel{\·}{\mathrm{\δ}}}_s,M_d=\{0,1\}\\ 0,M_d=\{2,3\}\end{array}\right.---\left(39\right)$
其中：分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼力，β_n为
岩体细观颗粒线性接触的法向阻尼系数，β_s为岩体细观颗粒线性接触的切向阻尼系数，k_n为
岩体细观颗粒线性接触的法向线性刚度，k_s为岩体细观颗粒线性接触的切向线性刚度，
分别为岩体细观颗粒线性接触的法向速率和切向速率，F^*为岩体细观颗粒线性接
触的全法向阻尼力，表达式为m_c为岩体细观等效颗粒质量，m⁽¹⁾为岩体
细观颗粒接触端1的细观颗粒质量，m⁽²⁾为岩体细观颗粒接触端2的细观颗粒质量，F^d为岩体
细观颗粒线性接触的总阻尼力。

考虑弯矩贡献因子的二维时效破裂模型的构建方法

技术领域

本发明涉及工程岩体细观时效破裂分析技术领域，具体地指一种考虑弯矩贡献因
子的二维时效破裂模型的构建方法。

背景技术

深部岩体工程开挖后的失稳与破坏往往不是在开挖后立刻发生的，一般都存在着
明显的变形破裂时效性和灾变(岩爆、大变形等)的滞后性，严重危害工程的施工安全与长
期运营。目前，在细观方面的时效力学研究成果相对较少。刘宁等对锦屏大理岩破裂的时间
效应进行了试验和数值分析(深埋大理岩破裂扩展时间效应的颗粒流模拟，岩石力学与工
程学报，2011,Vol.30No.10:1989-1996)；孙金山等应用蠕变细观力学模型对锦屏大理岩短
期和长期强度特征进行了数值研究(锦屏大理岩蠕变损伤演化细观力学特征的数值模拟研
究，岩土力学，2013,Vol.34No.12:3601-3608)。这些都是以离散元中平行粘结模型
(Parallel-Bonded Model，PBM)为基础，根据驱动应力和裂纹扩展速度间的关系来描述岩
石细观层面上的时效破裂。但是，这类平行粘结模型存在很多不足之处：首先，平行粘结断
裂之后，颗粒之间只考虑滑动摩擦力，没有考虑断裂颗粒间的接触方式，不符合岩体破裂后
的不同细观颗粒在外荷载下的运动方式；其次，颗粒间的剪切破裂准则是一条与平行粘结
正应力平行的水平直线，也即这种剪切破裂准则与平行粘结正应力状态无关，只要平行粘
结剪切应力大于或等于固定平行粘结剪切破裂强度，颗粒间即可发生剪切破裂，无法体现
岩体中不同平行粘结正应力具有不同平行粘结剪切破裂强度的客观事实；另外，对于岩体
这类摩擦粘结性材料，上述这种平行粘结模型并没有考虑粘结力矩的差异作用对接触破坏
的影响，将粘结力矩的贡献度对不同岩性的影响均视为一致，夸大了粘结力矩对岩体的破
坏作用。

发明内容

本发明的目的在于针对上述缺陷，提出了一种考虑弯矩贡献因子的二维时效破裂
模型的构建方法，本发明适应于应力和裂纹扩展速度之间的关系符合指数型的这类岩体，
对于平面状态下这类深部岩体工程的围岩长期稳定性预测、评价以及优化设计提供技术支
持。

本发明的目的是通过如下措施来达到的：一种考虑弯矩贡献因子的二维时效破裂
模型的构建方法，其特殊之处在于，包括如下步骤：

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

其中：为岩体细观颗粒二维平行粘结半径，为岩体细观颗粒二维平行粘结直径
乘数，A为岩体细观颗粒二维平行粘结面积圆，I为岩体细观颗粒二维平行粘结惯性矩；

步骤201：利用岩体细观颗粒二维平行粒粘结时效衰减劣化的初始时间步长增量
Δt，通过指数型函数计算岩体细观颗粒二维平行粘结时效衰减劣化的直径，公式(4)来确
定：

其中：为考虑弯矩贡献因子的二维平行粘结法向应力，为判断岩体细观颗粒
二维平行粘结开始时效劣化衰减时的应力阀值，为岩体细观颗粒二维平行粘结的拉伸强
度，为考虑弯矩贡献因子的二维平行粘结应力比，β₁、β₂分别为岩体细观颗粒平行粘结
时效劣化的第一控制参数、第二控制参数，为岩体细观颗粒二维平行粘结随时间劣化衰
减的直径，为岩体细观颗粒二维平行粘结未衰减时的直径，e为自然常数，Δ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是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩
体细观颗粒二维平行粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循
环计算中，包含时间效应的岩体细观颗粒二维平行粘结衰减后未破裂的最末标号值，i为第
一个至最后一个二维平行粘结颗粒标号值，分别为第i个岩体细观颗
粒二维平行粘结接触的a端和b端的绝对运动速度和角速度，n_n、n_s为岩体细观颗粒二维平行
粘结接触面的法向单位向量和切向单位向量，分别为岩体细观颗粒二维平行
粘结法向位移增量、切向位移增量，为岩体细观颗粒二维平行粘结法向刚度，为岩
体细观颗粒二维平行粘结弯矩增量；

步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和
公式(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒平行粘结未
破裂并包含时间效应的二维平行粘结法向力、切向力和切向弯矩；通过公式(12)、公式
(13)、公式(14)计算第i个岩体细观颗粒接触的二维平行粘结法向力、切向力和切向弯矩；
在二维情况下，通过公式(15)来确定岩体细观颗粒平行粘结法向弯矩：

法向力：

切向力：

切向弯矩：

法向弯矩：

其中：分别为第i个岩体细观颗粒包含时间效
应的平行粘结法向力、平行粘结切向力、包含时间效应的平行粘结法向弯矩、平行粘结切向
弯矩、平行粘结法向位移增量和平行粘结切向位移增量，为岩体细观颗粒二维平行粘结
切向刚度，+＝为加法自反运算符，-＝为减法自反运算符；

步骤206：设置弯矩贡献因子考虑弯矩对岩体细观颗粒平行粘结法向
正应力的贡献程度，根据二维平行粘结正应力计算公式和二维平行粘结
剪应力计算公式同时将这两个公式中A、I以及用A'、I'及替换，然后将步骤
203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型时间效应和弯
矩贡献因子的二维平行粘结正应力和二维平行粘结剪应力计算公式，分别见公式
(16)和公式(17)；

步骤207：将步骤206中包含时间效应的代入公式(18)，可确定考虑弯矩
贡献因子且带拉伸截止限的摩尔-库伦时效破裂准则，并且依次计算第j个至第k个的二维
平行粘结应力，用于判断岩体细观颗粒平行粘结是否破裂以及破裂模式；在该准则的岩体
细观颗粒二维平行粘结应力中包含了指数型时间效应和考虑弯矩贡献因子；

其中：f^s、fⁿ分别为岩体细观颗粒二维平行粘结的时效剪切破裂准则、时效拉伸破
裂准则，为第i个接触的含指数型时间效应的二维平行粘结剪应力，为第i个接
触的含指数型时间效应且考虑弯矩贡献因子的二维平行粘结正应力，分别为岩体
细观颗粒二维平行粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维平行粘结的粘聚力，
为岩体细观颗粒二维平行粘结的内摩擦角；f^s大于等于0表示平行粘结剪切破裂，小于0表
示平行粘结未发生剪切破裂；fⁿ大于等于0表示平行粘结拉伸破裂，小于0表示平行粘结未
发生拉伸破裂；

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

优选地，所述岩体细观颗粒二维平行粘结时效衰减劣化的初始时间步长增量Δt
的确定，是采用考虑弯矩贡献因子的平行粘结时效劣化衰减的二维指数型模式，由每个时
间步内的二维平行粘结首次衰减破裂所损耗的时间来确定，也即通过第一个平行粘结按指
数型模式进行衰减破裂所历时的时间除以直至第一个平行粘结破裂所需要的计算循环次
数来估算初始时间步长增量Δt，见公式其
中，为第i个接触的岩体细观颗粒二维平行粘结直径乘数，n_c为第
一个岩体细观颗粒二维平行粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体细观颗粒二
维平行粘结拉伸强度对应的时效劣化因子、二维平行粘结剪切强度对应的时效劣化因子，i
依次为第一个至最后一个岩体细观颗粒平行粘结数，∞为无穷大。

优选地，所述岩体细观颗粒二维平行粘结拉伸强度对应的时效劣化因子β_σ和岩体
细观颗粒二维平行粘结剪切强度对应的时效劣化因子β_τ的确定包括如下步骤，其中，这些
步骤中包含的公式下标1代表第一个按指数型模式进行时效衰减劣化的二维平行粘结破裂
标号；

步骤211：在二维情况下，由岩体细观颗粒平行粘结两端颗粒的速度、角速度和给
定的循环计算时间步长增量Δt_c，通过公式确定平行粘结接触的相
对转角通过公式确定平行粘结法向增量位移通
过公式确定平行粘结切向增量位移通过公式
确定平行粘结接触的弯矩增量；

步骤212：根据步骤211中的公式通过公式
确定平行粘结法向力；根据步骤211中的公式
通过公式确定平行粘结切向力；根据步骤211中的公式
和公式通过公式确定
平行粘结切向弯矩；通过公式确定平行粘结法向弯矩，其中，+＝为加法自反运算
符，-＝为减法自反运算符；

步骤213：在二维情况下，通过公式确定平行粘结正应
力，通过公式确定平行粘结剪应力，将这两个公式中A、I以及用A'、I'及替
换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型
时间效应和弯矩贡献因子的二维平行粘结正应力计算公式和
包含指数型时间效应的二维平行粘结剪应力计算公式

步骤214：将代入公式并令β＝β_σ；将
代入公式并令β＝β_τ，据此，分别得到所述岩体细观颗粒二维平行粘结拉伸强
度对应的时效劣化因子以及岩体细观颗粒二维
平行粘结剪切强度对应的时效劣化因子

优选地，所述岩体细观颗粒平行粘结发生破裂后，岩体细观颗粒的运动模式采用
考虑阻尼效应的二维线性接触模型来表达，用于描述岩体细观颗粒平行粘结时效破裂后细
观颗粒的应力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：

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

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

步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)
计算，如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标
及距离计算，如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触
段的单位向量：

其中：n_i为接触的单位矢量，为接触端b的方向向量，为接触端a的方向向量，
n_wall为约束墙的方向向量。

步骤303：岩体细观颗粒平行粘结破裂后，二维线性接触段的接触重叠量U，通过步
骤301计算颗粒间距d以及二维线性接触两端的颗粒半径R^a、R^b，再利用公式(21)来确定。通
过设定颗粒二维线性接触参考距离g_r，并结合公式(22)，确定颗粒二维线性接触的距离g_s：

g_s＝|U|-g_r (22)

步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体
或是墙体的刚度k^a，k^b等效代替为接触点的刚度(法向刚度和切向刚度的统称)，按公式(23)
计算：

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

步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)
来计算。其中e_pqz为Ricci指标置换符号，按照公式(26)来计算：

其中：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端两个接触端)。

步骤306：对于线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计
最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋
于稳定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增
量和切向位移增量：

R＝min(R_a,R_b) (27)

ΔU_p1＝V_p1Δt (30)

其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的
转动惯量；k^平为岩体细观颗粒系统平移刚度，k^转为岩体细观颗粒系统旋转刚度；ΔU_p1为岩体
细观颗粒二维线性接触的总位移增量，Δδ_n、物理意义相同，均表示岩体细观颗粒二维
线性接触的法向位移增量，Δδ_s、物理意义相同，均表示岩体细观颗粒二维线性接触的
切向位移增量，V_p1与V_q1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1,q1
为张量指标变换符号。

步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向
和切向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一
步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的
更新是采用前一步的切向位移增量与更新因子α的乘积获得。

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

步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加(M_l＝1)和
绝对矢量累加(M_l＝0)模式，通过公式(33)、(34)计算获得；岩体细观颗粒二维线性接触的
切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(35)计算获得：

其中：分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向
线性力，k_n、k_s分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδ_n、
Δδ_s分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量，分别
为岩体细观颗粒二维线性接触的初始法向力增量值、切向力增量值，为岩体细观颗粒未
滑动时的静摩擦力，为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与乘积得到。

步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式M_d＝{0,2}和无拉伸
模式M_d＝{1,3}两种，通过公式(39)计算，其中m_c为等效颗粒质量，按公式(40)计算，岩体细
观颗粒线性接触的切向阻尼采用全剪切模式M_d＝{0,1}和滑-剪模式M_d＝{2,3}，按照公式
(41)来计算，

其中：分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼
力，β_n为岩体细观颗粒线性接触的法向阻尼系数，β_s为岩体细观颗粒线性接触的切向阻尼系
数，k_n为岩体细观颗粒线性接触的法向线性刚度，k_s为岩体细观颗粒线性接触的切向线性刚
度，分别为岩体细观颗粒线性接触的法向速率和切向速率，F^*为岩体细观颗粒线性
接触的全法向阻尼力，表达式为m_c为岩体细观等效颗粒质量，m⁽¹⁾为岩
体细观颗粒接触端1的细观颗粒质量，m⁽²⁾为岩体细观颗粒接触端2的细观颗粒质量，F^d为岩
体细观颗粒线性接触的总阻尼力。

本发明提出的考虑弯矩贡献因子的二维时效破裂模型的构建方法，在平行粘结时
效破裂之后，通过在颗粒之间嵌入考虑阻尼效应的二维线性模型来表达颗粒间的接触方
式，不仅可以考虑滑动摩擦作用，而且还可以考虑颗粒间的变形特性，符合岩体颗粒在平面
状态下的运动规律；其次，在平行粘结应力计算中引入弯矩贡献因子，不仅考虑弯矩了对平
行粘结正应力贡献，而且还考虑了弯矩对岩体长期强度的影响；另外，在所构建的二维时效
破裂模型中嵌入考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦时效破裂准则，该准则平
行粘结应力中包含了指数型时间效应且增加了弯矩贡献因子，不仅可以描述与平行粘结正
应力相关时效剪切破裂强度的差异，还可以对时效拉伸破裂进行合理的表达，且考虑了弯
矩对平行粘结时效破裂的影响，符合岩体在平面状态下时效破裂的基本特征。本发明对于
平面状态下深部岩体工程围岩长期稳定性预测、评价以及优化设计提供直接的技术支持。

本发明所提出的一种考虑弯矩贡献因子的二维时效破裂模型的构建方法，其有益
效果和优势主要体现在：

(1)本发明中通过在二维平行粘结正应力计算公式中设置弯矩贡献因子，不仅考
虑了弯矩对平行粘结正应力的贡献程度，而且还考虑了弯矩对岩体长期强度的影响，同时
减小岩体破裂过程中产生过强的能量冲击波对邻近区域造成的二次损伤，适合描述平面应
力或平面应变条件下岩体的细观力学破裂行为。

(2)本发明中通过构建考虑弯矩贡献因子的二维平行粘结直径时效劣化衰减模
式，包括设置指数型与考虑弯矩贡献因子的平行粘结应力相关的二维平行粘结直径劣化衰
减模式，设置平行粘结直径随着时间增加逐步劣化衰减模式，同时设置平行粘结的面积和
截面惯性矩相应的时效劣化衰减模式。这种构建模式适合描述平面应力或平面应变条件下
深部岩体的细观力学时效破裂机制和响应规律。

(3)本发明中在所构建的二维时效破裂模型中，嵌入考虑弯矩贡献效应且带拉伸
截止限的摩尔-库伦时效破裂准则。在该准则平行粘结应力中包含指数型时间效应且增加
了弯矩贡献因子，不仅可以描述与平行粘结正应力相关时效剪切破裂强度的差异，还可以
对时效拉伸破裂进行合理的表达，且考虑了弯矩对平行粘结时效破裂的影响，符合平面条
件下岩体时效破裂模式。

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

附图说明

图1是本发明模型中岩体颗粒与颗粒接触示意图。

图2是本发明模型中岩体颗粒与刚性墙接触示意图。

图3是本发明模型中岩体颗粒重叠状态示意图。

图4是本发明模型中岩体颗粒刚度计算示意图。

图5是本发明模型中线性切向力和切向位移示意图

图6是本发明中所构建的二维时效破裂物理模型示意图。

图7是本发明模型中线性平行粘结结构示意图。

图8是本发明模型中的考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦破裂准则
示意图。

图9是本发明模型中直径时效衰减示意图。

图10是本发明模型直径时效衰减率对数与应力曲线示意图。

图11是本发明模型二维接触面的法向和切向单位向量示意图。

图12是本发明中模型结构构建流程示意图。

图13是本发明模型蠕变损伤分布图。

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

图15是本发明模型不同力矩贡献度的蠕变位移与时间关系曲线图。

图中：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代表为对数更新率ln(γ)，40代表判断岩体颗粒开始时效劣化衰减时的
应力阀值41代表拉伸强度42代表考虑弯矩贡献因子的平行粘结应力比43
代表控制岩体细观颗粒平行粘结时效劣化的第二控制参数β₂，44代表二维接触面的法向向
量n_n，45代表二维接触面切向单位向量n_s。

具体实施方式

下面结合附图和具体构建步骤及实施实例，对本发明模型进行详细的阐述。实例
说明仅是辅助用于本发明的理解，而不限定本发明的实际应用范围。在阅读了本发明以后，
本领域的技术人员对本发明的各种等价形式的修改均属于本发明所申请的权利要求所限
定的范围。

注：说明书中的所有的标号前面写明了公式，如公式(1)，均为公式标号。

如图1～图11所示，本发明考虑弯矩贡献因子的二维时效破裂模型适应于二维颗
粒离散元、二维颗粒不连续变形分析方法、二维颗粒流形元；二维时效破裂模型包括考虑弯
矩贡献因子的二维平行粘结应力模式、考虑弯矩贡献因子的二维平行粘结直径时效劣化衰
减模式、考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦时效破裂准则以及岩体时效破裂
后考虑阻尼效应的二维线性接触模型。

考虑弯矩贡献因子的二维平行粘结应力模式是指二维平行粘结正应力计算公式
中设置了弯矩贡献因子考虑弯矩对二维平行粘结正应力的贡献程
度；

第i个接触的岩体细观颗粒平行粘结法向力的计算方法为：第i
个接触的岩体细观颗粒平行粘结切向弯矩的计算方法为：

在上述公式中，为第i个岩体细观颗粒二维平行粘结的正应力，分别
为第i个接触的岩体细观颗粒平行粘结法向力、切向弯矩；为岩体细观颗粒二维平行粘结
半径，为弯矩贡献因子，I为岩体细观颗粒二维平行粘结的惯性矩，A为岩体细观
颗粒二维平行粘结面积，i依次为第一个至最后一个岩体细观颗粒平行粘结数，为岩体
细观颗粒二维平行粘结法向刚度，为岩体细观颗粒二维平行粘结法向位移增量，
为岩体细观颗粒二维平行粘结切向相对转角增量，+＝为加法自反运算符，-＝为减法自反
运算符，法向弯矩

考虑弯矩贡献因子的二维平行粘结直径时效劣化衰减模式包括考虑弯矩贡献因
子的二维平行粘结时效劣化衰减模式，在岩体细观颗粒二维平行粘结时效劣化衰减时，设
置了与考虑弯矩贡献因子的平行粘结应力相关的指数型模式，这种指数型模式中的二维平
行粘结直径随时间逐步劣化衰减，见平行粘结直径公式式
中，为考虑弯矩贡献因子的二维平行粘结法向应力，为判断岩体细观颗粒二维平行粘
结开始时效劣化衰减时的应力阀值，为岩体细观颗粒二维平行粘结的拉伸强度，为
考虑弯矩贡献因子的二维平行粘结应力比，β₁、β₂分别为岩体细观颗粒平行粘结时效劣化的
第一控制参数、第二控制参数，为岩体细观颗粒二维平行粘结随时间劣化衰减的直径，
为岩体细观颗粒二维平行粘结未衰减时的直径，e为自然常数，Δt为岩体细观颗粒二维平
行粘结时效衰减劣化的时间增量；设置了岩体细观颗粒二维平行粘结面积和惯性矩时效劣
化衰减模式，分别见平行粘结单位厚度为1时的平行粘结面积计算公式平行
粘结单位厚度为1时的惯性矩计算公式其中，β为岩体细观颗粒二维平行粘
结直径的指数型时效衰减因子，其计算公式见
式中，A'、I'、分别表示为岩体细观颗粒二维平行粘结随时间劣化衰减的平行粘
结直径、平行粘结半径、平行粘结面积、平行粘结惯性矩、平行粘结直径乘数，
A、I、为岩体细观颗粒平行粘结未衰减时的平行粘结直径、平行粘结半径、平行粘
结面积、平行粘结惯性矩、平行粘结直径乘数；

同时按照这种二维指数型时效劣化衰减模式估算岩体细观颗粒二维平行粘结破
裂的初始时间步长，见公式其中，
为第i个接触的岩体细观颗粒二维平行粘结直径乘数，n_c为第一个岩体细观颗粒二维平行
粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体细观颗粒二维平行粘结拉伸强度对应的
时效劣化因子、二维平行粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个岩
体细观颗粒平行粘结数，∞为无穷大。岩体细观颗粒二维平行粘结拉伸强度对应的时效劣
化因子β_σ和岩体细观颗粒二维平行粘结剪切强度对应的时效劣化因子β_τ的计算公式分别为

其中，分别为第i个颗粒接触的平行粘结法向力、平行粘结切向力、
平行粘结切向弯矩，为岩体细观颗粒二维平行粘结的粘聚力，为岩体细观颗粒二维平行
粘结的内摩擦角。

考虑弯矩贡献效应且带拉伸截止限的摩尔-库伦时效破裂准则是指在判断岩体细
观颗粒二维平行粘结时效破裂时，采用所嵌入的考虑弯矩贡献效应且带拉伸截止限的摩
尔-库伦时效破裂强度准则来判断，见公式

其中，f^s、fⁿ分别为岩体细观颗粒二维平行粘结的时效剪切破裂准则、时效拉伸破
裂准则，分别为岩体细观颗粒二维平行粘结拉伸强度、抗剪强度，分
别为第i个接触的含指数型时间效应且考虑弯矩贡献因子的岩体细观颗粒二维平行粘结正
应力、剪应力，

第i个接触的含指数型时间效应且考虑弯矩贡献因子的岩体细观颗粒二维平行粘
结正应力的计算公式为第i个接触的含指数型时间效应的岩
体细观颗粒二维平行粘结剪应力的计算公式为

在该准则的二维平行粘结应力中包含了指数型时间效应，见岩体细观颗粒二维平
行粘结直径的指数型时效衰减因子计算公式β₁、β₂分
别为控制岩体细观颗粒平行粘结时效劣化的第一控制参数、第二控制参数；

f^s大于等于0表示岩体细观颗粒二维平行粘结剪切破裂，小于0表示岩体细观颗粒
二维平行粘结未发生剪切破裂；fⁿ大于等于0表示岩体细观颗粒二维平行粘结拉伸破裂，小
于0表示岩体细观颗粒二维平行粘结未发生拉伸破裂。

岩体时效破裂后考虑阻尼效应的二维线性接触模型是指岩体细观颗粒平行粘结
时效破裂后，通过给定的二维线性接触参考距离g_r设置了细观颗粒二维线性接触距离g_s，见
岩体细观颗粒二维线性接触距计算公式其
中，为岩体内部颗粒与颗粒二维线性接触端a的坐标，为岩体内部颗粒与颗粒
二维线性接触端b的坐标，R^a、R^b分别为岩体细观二维线性接触端a的颗粒半径和二维线性接
触端b的颗粒半径；设置考虑岩体细观颗粒间受力变形的二维线性接触模式，在岩体颗粒之
间设置考虑二维滑动摩擦线力的作用模式，岩体细观颗粒间受力变形的二维线性接触法向
线性力计算公式取M_l＝1为相对矢量累加模式，取M_l＝0为
绝对矢量累加模式，岩体细观颗粒间受力变形的二维线性接触切向线性力计算公式为
其中，分别为岩体细观颗粒间受力变形的二维线性接触
法向线性力、切向线性力，k_n、k_s分别为岩体细观颗粒间受力变形的二维线性接触法向、切向
线性刚度，Δδ_n、Δδ_s分别为法向位移增量、切向位移增量，分别为初始法向力
增量值和切向力增量值，为颗粒未滑动时的静摩擦力，为岩体细观
颗粒滑动摩擦力，通过摩擦系数μ与乘积得到；同时设置二维接触的四种阻尼模式，其中
法向阻尼采用全法向模式M_d＝{0,2}和无拉伸模式M_d＝{1,3}两种，通过公式
计算，其中m_c为等效颗粒质量，按公式
计算，切向阻尼采用全剪切模式M_d＝{0,1}和滑-剪模式M_d＝{2,3}，按照公式
来计算，其中：分别为法向阻尼力、切向阻尼力，β_n为
法向阻尼系数、β_s为切向阻尼系数，为法向速率、切向速率，F^*为岩体细观颗粒线性
接触的全法向阻尼力，表达式为m_c为等效颗粒质量，m⁽¹⁾为二维线性接
触端1的岩体细观颗粒质量，m⁽²⁾为二维线性接触端2的岩体细观颗粒质量。

本发明考虑弯矩贡献因子的二维时效破裂模型的构建方法，包括如下步骤：

步骤1：设定岩体细观颗粒平行粘结接触的几何参数量包括平行粘结面积和平行
粘结惯性矩，R^a、R^b分别为二维平行粘结接触a端的颗粒半径、b端的颗粒半径，为岩体细观
颗粒平行粘结直径或半径乘数，在二维情况下，平行粘结单位厚度为1时的平行粘结面积A
和平行粘结惯性矩I分别通过公式(2)、公式(3)来确定：

其中：为岩体细观颗粒二维平行粘结半径，为岩体细观颗粒二维平行粘结直径
乘数，A为岩体细观颗粒二维平行粘结面积圆，I为岩体细观颗粒二维平行粘结惯性矩；

步骤201：利用岩体细观颗粒二维平行粘结时效衰减劣化的初始时间步长增量Δ
t，通过指数型函数计算岩体细观颗粒二维平行粘结时效衰减劣化的直径，公式(4)来确定：

其中：为考虑弯矩贡献因子的二维平行粘结法向应力，为判断岩体细观颗
粒二维平行粘结开始时效劣化衰减时的应力阀值，为岩体细观颗粒二维平行粘结的拉伸
强度，为考虑弯矩贡献因子的二维平行粘结应力比，β₁、β₂分别为岩体细观颗粒平行粘
结时效劣化的第一控制参数、第二控制参数，为岩体细观颗粒二维平行粘结随时间劣化
衰减的直径，为岩体细观颗粒二维平行粘结未衰减时的直径，e为自然常数，Δ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_c，通过公式(8)、公式(9)、公式(10)确定第i个岩体细观颗粒二维平行粘结相
对转角二维平行粘结法向增量位移以及二维平行粘结切向增量位移再
结合步骤203中的公式(7)和步骤202中的公式(5)，确定第i个岩体细观颗粒包含时间效应
的二维平行粘结弯矩增量，具体见公式(11)：

其中，ff、j、k是自然数，且2≤j≤ff≤k，j为每次循环计算中，包含时间效应的岩
体细观颗粒二维平行粘结衰减后未破裂的初始标号值，ff为中间某一个标号值，k为每次循
环计算中，包含时间效应的岩体细观颗粒二维平行粘结衰减后未破裂的最末标号值，i为第
一个至最后一个平行粘结颗粒标号值，分别为第i个岩体细观颗粒二
维平行粘结接触的a端和b端的绝对运动速度和角速度，n_n和n_s为岩体细观颗粒二维平行粘
结接触面的法向单位向量和切向单位向量，和分别为岩体细观颗粒二维平行粘
结法向位移增量和切向位移增量，为岩体细观颗粒平行粘结法向刚度，为岩体细观
颗粒二维平行粘结弯矩增量。

岩体细观颗粒二维平行粘结时效衰减劣化的初始时间步长增量Δt的确定，是采
用考虑弯矩贡献因子的平行粘结时效劣化衰减的二维指数型模式，由每个时间步内的二维
平行粘结首次衰减破裂所损耗的时间来确定，也即通过第一个平行粘结按指数型模式进行
衰减破裂所历时的时间除以直至第一个平行粘结破裂所需要的计算循环次数来估算初始
时间步长增量Δt，见公式其中，
为第i个接触的岩体细观颗粒二维平行粘结直径乘数，n_c为第一个岩体细观颗粒二维平
行粘结破裂所需的循环计算的次数，β_σ、β_τ分别为岩体细观颗粒二维平行粘结拉伸强度对应
的时效劣化因子、二维平行粘结剪切强度对应的时效劣化因子，i依次为第一个至最后一个
岩体细观颗粒平行粘结数，∞为无穷大。

岩体细观颗粒二维平行粘结拉伸强度对应的时效劣化因子β_σ和岩体细观颗粒二
维平行粘结剪切强度对应的时效劣化因子β_τ的确定包括如下步骤，其中，这些步骤中包含
的公式下标1代表第一个按指数型模式进行时效衰减劣化的二维平行粘结破裂标号：

步骤211：在二维情况下，由岩体细观颗粒平行粘结两端颗粒的速度、角速度和给
定的循环计算时间步长增量Δt_c，通过公式确定平行粘结接触的相
对转角通过公式确定平行粘结法向增量位移通
过公式确定平行粘结切向增量位移通过公式
确定平行粘结接触的弯矩增量；

步骤212：根据步骤211中的公式通过公式
确定平行粘结法向力；根据步骤211中的公式
通过公式确定平行粘结切向力；根据步骤211中的公式
和公式通过公式确定
平行粘结切向弯矩；通过公式确定平行粘结法向弯矩，其中，+＝为加法自反运算
符，-＝为减法自反运算符；

步骤213：在二维情况下，通过公式确定平行粘结正应
力，通过公式确定平行粘结剪应力，将这两个公式中A、I以及用A'、I'及替
换，然后将步骤203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型
时间效应和弯矩贡献因子的二维平行粘结正应力计算公式和
包含指数型时间效应的二维平行粘结剪应力计算公式

步骤214：将代入公式并令β＝β_σ，将
代入公式并令β＝β_τ，据此，分别得到所述岩体细观颗粒二维平行
粘结拉伸强度对应的时效劣化因子以及岩体细
观颗粒二维平行粘结剪切强度对应的时效劣化因子

步骤205：根据步骤203中的公式(6)和公式(7)、步骤204中的公式(8)、公式(9)和
公式(11)以及步骤202中的公式(5)，依次更新计算第j个至第k个岩体细观颗粒平行粘结未
破裂并包含时间效应的二维平行粘结法向力、切向力和切向弯矩；通过公式(12)、公式
(13)、公式(14)计算第i个岩体细观颗粒接触的二维平行粘结法向力、切向力和切向弯矩；
在二维情况下，通过公式(15)来确定岩体细观颗粒平行粘结法向弯矩：

法向力：

切向力：

切向弯矩：

法向弯矩：

其中：分别为第i个岩体细观颗粒包含时间效
应的平行粘结法向力、平行粘结切向力、包含时间效应的平行粘结法向弯矩、平行粘结切向
弯矩、平行粘结法向位移增量和平行粘结切向位移增量，为岩体细观颗粒二维平行粘结
切向刚度，+＝为加法自反运算符，-＝为减法自反运算符。

步骤206：设置弯矩贡献因子考虑弯矩对岩体细观颗粒平行粘结法向
正应力的贡献程度，根据二维平行粘结正应力计算公式和二维平行粘结
剪应力计算公式同时将这两个公式中A、I以及用A'、I'及替换，然后将步骤
203中的公式(6)和公式(7)以及步骤202中的公式(5)代入，获得包含指数型时间效应和弯
矩贡献因子的二维平行粘结正应力和二维平行粘结剪应力计算公式，分别见公式
(16)和公式(17)。

步骤207：将步骤206中包含时间效应的代入公式(18)，可确定考虑弯矩
贡献因子且带拉伸截止限的摩尔-库伦时效破裂准则，并且依次计算第j个至第k个的二维
平行粘结应力，用于判断岩体细观颗粒平行粘结是否破裂以及破裂模式；在该准则的岩体
细观颗粒平行粘结应力中包含了指数型时间效应和考虑弯矩贡献因子。

其中：f^s、fⁿ分别为岩体细观颗粒二维平行粘结的时效剪切破裂准则、时效拉伸破
裂准则，为第i个接触的含指数型时间效应的二维平行粘结剪应力，为第i个接
触的含指数型时间效应且考虑弯矩贡献因子的二维平行粘结正应力，分别为岩体
细观颗粒二维平行粘结的拉伸强度、抗剪强度，为岩体细观颗粒二维平行粘结的粘聚力，
为岩体细观颗粒二维平行粘结的内摩擦角；f^s大于等于0表示平行粘结剪切破裂，小于0表
示平行粘结未发生剪切破裂；fⁿ大于等于0表示平行粘结拉伸破裂，小于0表示平行粘结未
发生拉伸破裂；

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

岩体细观颗粒平行粘结发生破裂后，岩体细观颗粒的运动模式采用考虑阻尼效应
的二维线性接触模型来表达，用于描述岩体细观颗粒平行粘结时效破裂后细观颗粒的应
力、变形及运行规律，考虑阻尼效应的二维线性接触模型的构建包括如下步骤：

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

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

步骤302：二维平面状态下岩体细观颗粒间每个接触点的单位向量通过公式(20)
计算，如果是颗粒与颗粒之间的接触，利用步骤301中得到二维线性接触两端的中心点坐标
及距离计算，如果是颗粒与墙接触，直接利用墙体的法向量等效替换来计算，确定每个接触
段的单位向量：

其中：n_i为接触的单位矢量，为接触端b的方向向量，为接触端a的方向向量，
n_wall为约束墙的方向向量。

步骤303：岩体细观颗粒平行粘结破裂后，二维线性接触段的接触重叠量U，通过步
骤301计算颗粒间距d以及二维线性接触两端的颗粒半径R^a、R^b，再利用公式(21)来确定。通
过设定颗粒二维线性接触参考距离g_r，并结合公式(22)，确定颗粒二维线性接触的距离g_s：

g_s＝|U|-g_r (22)

步骤304：确定岩体细观颗粒接触点法向、切向等效刚度，利用接触两端颗粒实体
或是墙体的刚度k^a，k^b等效代替为接触点的刚度(法向刚度和切向刚度的统称)，按公式(23)
计算：

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

步骤305：确定岩体中接触两端颗粒间的相对运动速度，利用公式(24)、公式(25)
来计算。其中e_pqz为Ricci指标置换符号，按照公式(26)来计算：

其中：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端两个接触端)。

步骤306：对于线性接触模型的初始时间步长增量Δt的取值，通过公式(29)估计
最小的时间步长Δt，确保所构建模型的计算时间步长小于该值，即可保证系统积分计算趋
于稳定；通过公式(30)、公式(31)、公式(32)确定每个线性接触的总位移增量、法向位移增
量和切向位移增量：

R＝min(R_a,R_b) (27)

ΔU_p1＝V_p1Δt (30)

其中：R为岩体细观颗粒的等效半径，m为岩体细观颗粒质量，J1为岩体细观颗粒的
转动惯量；k^平为岩体细观颗粒系统平移刚度，k^转为岩体细观颗粒系统旋转刚度；ΔU_p1为岩体
细观颗粒二维线性接触的总位移增量，Δδ_n、物理意义相同，均表示岩体细观颗粒二维
线性接触的法向位移增量，Δδ_s、物理意义相同，均表示岩体细观颗粒二维线性接触的
切向位移增量，V_p1与V_q1为岩体细观颗粒接触两端的相对运动速度，n为单位法向量，p1,q1
为张量指标变换符号。

步骤307：由公式(22)判定岩体细观颗粒表面接触允许存在的最大距离，计算法向
和切向位移更新因子，另外，岩体细观颗粒二维线性接触法向位移增量的更新是采用前一
步的法向位移增量与更新因子α的乘积获得，岩体细观颗粒二维线性接触切向位移增量的
更新是采用前一步的切向位移增量与更新因子α的乘积获得。

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

步骤308：岩体细观颗粒二维线性接触的法向线性力采取相对矢量累加(M_l＝1)和
绝对矢量累加(M_l＝0)模式，通过公式(34)、(35)计算获得；岩体细观颗粒二维线性接触的
切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(36)计算获得：

其中：分别为岩体细观颗粒间受力变形的二维线性接触法向线性力、切向
线性力，k_n、k_s分别为岩体细观颗粒间受力变形的二维线性接触法向、切向线性刚度，Δδ_n、
Δδ_s分别为岩体细观颗粒二维线性接触的法向位移增量、切向位移增量，分别
为岩体细观颗粒二维线性接触的初始法向力增量值、切向力增量值，为岩体细观颗粒未
滑动时的静摩擦力，为岩体细观颗粒滑动摩擦力，其值可通过摩擦系数μ与乘积得到。

步骤309：岩体细观颗粒线性接触的法向阻尼采用全法向模式M_d＝{0,2}和无拉伸
模式M_d＝{1,3}两种，通过公式(37)计算，其中m_c为等效颗粒质量，按公式(38)计算，岩体细
观颗粒线性接触的切向阻尼采用全剪切模式M_d＝{0,1}和滑-剪模式M_d＝{2,3}，按照公式
(39)来计算，

其中：分别为岩体细观颗粒线性接触的法向线性阻尼力、切向线性阻尼
力，β_n为岩体细观颗粒线性接触的法向阻尼系数，β_s为岩体细观颗粒线性接触的切向阻尼系
数，k_n为岩体细观颗粒线性接触的法向线性刚度，k_s为岩体细观颗粒线性接触的切向线性刚
度，分别为岩体细观颗粒线性接触的法向速率和切向速率，F^*为岩体细观颗粒线性
接触的全法向阻尼力，表达式为m_c为岩体细观等效颗粒质量，m⁽¹⁾为岩
体细观颗粒接触端1的细观颗粒质量，m⁽²⁾为岩体细观颗粒接触端2的细观颗粒质量，F^d为岩
体细观颗粒线性接触的总阻尼力。

实验实例

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

步骤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、法向阻尼系数Apfan(图6)、切向阻尼系数Apfas(图6)以及内摩擦角pb_phi(图8)等19
个参数，参数具体值见表一。

步骤3：生成岩石样本模型

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

步骤4：精确标定岩石样本的细观参数

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

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

对岩体进行一系列不同应力强度比条件下的时效力学试验，得到不同应力强度比
条件下岩体变形随时间演化的曲线(图14)。通过参数拟合法，匹配实际岩体的时效变形过
程，确定岩体时效劣化指数β₁、β₂。

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

在荷载一定的条件下，分别设定不同的弯矩贡献因子，获得了弯矩贡献因子对岩
体时效变形和破坏的影响规律(图15)。

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

上述实施例中，公式的符号与图1～图11以及附图说明中的符号相互对应。

其它未详细说明的部分均为现有技术，以上所有参数均可通过查阅手册或计算得
到。本发明并不严格地局限于上述实施例。以上所述仅为本发明的特定实施例，并不用于限
制本发明。凡在本发明的精神和原则之内所做的任何修改、等同替换及改进等，都在本发明
的保护范围之内。