岩体二维细观时效破裂幂函数型模型的构建方法与流程

技术特征：

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{\left(\frac{\stackrel{\‾}{\mathrm{\σ}}}{{\stackrel{\‾}{\mathrm{\σ}}}_c}\right)}^{{\mathrm{\β}}_2}\mathrm{\Δ}t,{\stackrel{\‾}{\mathrm{\σ}}}_{aa}\≤\stackrel{\‾}{\mathrm{\σ}}\<{\stackrel{\‾}{\mathrm{\σ}}}_c\end{array}\right.---\left(4\right)$

其中：为考虑弯矩贡献因子的二维粘结法向应力，为判断岩体细观颗粒二维粘结开始时效劣化衰减时的应力阀值，为岩体细观颗粒二维粘结的拉伸强度，为考虑弯矩贡献因子的二维粘结应力比，β₁、β₂分别为岩体细观颗粒粘结时效劣化的第一控制参数、第二控制参数，为岩体细观颗粒二维粘结随时间劣化衰减的直径，为岩体细观颗粒二维粘结未衰减时的直径，Δ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{\left(\frac{\stackrel{\‾}{\mathrm{\σ}}}{{\stackrel{\‾}{\mathrm{\σ}}}_c}\right)}^{{\mathrm{\β}}_2}\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，通过公式(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(8\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、i是自然数，且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)计算接触点a端，接触点b端的中心距离：

$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_{qpz}{\mathrm{\ω}}_q^{\left(b\right)}(x_z^{\left(c\right)}-x_z^{\left(b\right)}\left)\right)-({\stackrel{\·}{x}}_p^{\left(a\right)}+e_{qpz}{\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)模式，通过公式(33)、(34)计算获得；岩体细观颗粒二维线性接触的切向线性力采用岩体细观颗粒接触滑动来表示，通过公式(35)计算获得：

$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(33\right)$

$F_s^{*}={\left(F_s^l\right)}_o-k_s{\mathrm{\Δ\δ}}_s---\left(34\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(35\right)$

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

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

$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(36\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(38\right)$

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