层状岩体在我国广泛分布，其力学特性受层理倾角的影响，各向异性显著^[1]. 在隧道建设的过程中，层状构造的围岩使隧道易发生拱顶坍塌、初期支护开裂、钢拱架扭曲等灾害^[2]，因此研究层状岩体的力学特性具有重要的工程意义.

根据现有的隧道工程研究可知，开挖后的岩体变形并非瞬时完成，而要持续一定的时间才趋于稳定，即岩体的流变性质和时效性特征. 岩体流变力学中的蠕变效应是国内外众多学者研究范围最广、内容最深入的一个方面^[3]. Xu等^[4]考虑结构效应，分析层状岩体的蠕变特性，结果表明层状岩体蠕变参数受结构效应的影响显著. 范庆忠等^[5]对红砂岩进行单轴分级加载蠕变试验，研究发现岩体的瞬时弹性模量受应力水平的影响较大. 李永盛^[6]对4种不同力学性质的岩体进行单轴压缩蠕变试验，发现岩体的蠕变速率呈现3个阶段，各阶段持续时长与岩体性质和应力水平有关. 徐卫亚等^[7]建立非线性黏弹塑性模型，实现了对岩体加速蠕变的描述. 胡涛涛等^[8]对含软弱夹层炭质板岩的单轴压缩力学行为进行数值模拟研究，总结了端部约束对抗压强度及塑性区发展的影响规律. 周晓飞等^[9]借助幂函数，对泥质夹层剪切蠕变曲线进行拟合，提出基于幂级次求解泥质夹层长期强度的方法. Li等^[10]建立岩体损伤蠕变模型，该模型能够反映裂隙岩石在衰减、稳态和加速阶段蠕变曲线的非线性特征. 宋洋等^[11]在传统Nishihara模型中引入GTN模型，建立可描述节理岩体剪切蠕变全过程的蠕变模型. 赵晨阳等^[12]建立非线性元件并将其串联到Burgers模型中，实现了对含软弱结构面岩体三阶段的模拟. Trzeciak等^[13]采用Burgers模型，对波兰北部波罗的海盆地试样进行黏弹性分析，得到Burgers模型参数的各向异性特征. Lv等^[14]将Kachanov损伤理论引入到经典蠕变本构模型中，提出新的岩体蠕变损伤本构模型建立的方法. Tian等^[15]推导了层状岩石的横观各向同性蠕变模型，总结了层理面刚度比对蠕变行为的影响规律. 张向东等^[16]对冻黏土进行三轴蠕变试验，建立符合冻黏土蠕变特性的对数形蠕变方程. 周传涛等^[17]对含不同层理倾角的黄砂岩开展常规三轴压缩试验，得到其峰值强度随倾角呈现先降低后升高的趋势. 刘开云等^[18]将三轴压缩过程中岩石弹性模量的衰减方程引入本构方程，推导得到反映岩石蠕变全过程的三维非线性黏弹塑蠕变本构模型. He等^[19]针对有结构面的大理岩，提出描述损伤时变特征的时变统计损伤本构模型，验证了模型的合理性. 目前，虽然学者们通过开展蠕变试验和建立蠕变本构模型，研究节理、软弱夹层及层理倾角等单一因素对岩体流变特性的影响，但对于炭质板岩这一特定岩体，考虑内部裂损和层理倾角这2个因素对蠕变特性的影响的研究较少，现有的蠕变损伤模型往往不能全面、准确地描述岩石蠕变的全过程.

为了研究考虑层理倾角的炭质板岩非线性蠕变损伤本构模型，本文采用炭质板岩开展分级加载蠕变试验. 基于试验结果，考虑层理结构对炭质板岩蠕变过程的影响以及Nishihara模型无法描述岩石加速蠕变的问题，建立考虑层理倾角和裂隙损伤的改进Nishihara非线性蠕变损伤本构模型，为解决隧道工程建设和运营过程中的变形稳定问题提供理论依据.

不同层理倾角炭质板岩的蠕变试验在长安大学RTX-1000伺服控制岩石三轴流变系统上开展. 试样均采集自渭武高速木寨岭隧道，通过岩石取芯机在采集的岩样中取得不同角度的圆柱体试样. 试样高 为100 mm，直径为50 mm，通过打磨机将试样两端面打磨平整，控制端面平整度小于0.02 mm ，断面垂直度小于0.25°. 分别对层理倾角θ = 0°、45°和90° 3组炭质板岩试样在给定围压等级5、10和15 MPa下进行三轴分级加载蠕变试验，试样破坏图片如图1所示. 试验时先将围压加载至指定值，再施加轴压. 轴压按分级加载的形式施加，根据常规三轴压缩强度σ_c 确定每一级加载的应力，并以100 N/s的加载速率加载至三轴压缩强度的 20%、40%、60%、80%及 100%后保持24 h，具体的应力分级如表1所示. 表中，p为围压，σ_t为轴向应力. 常规三轴试验的结果表明， 0°层理倾角试样的峰值应变随着围压的增大而减小，炭质板岩表现出贯穿层理面的剪切破坏. 45°层理倾角试样的峰值应变随着围压的增大而增大，表现出沿层理面产生破坏面的弱面剪切破坏. 90°层理倾角试样的峰值应变随着围压先减小后增大，破坏模式呈现张拉-剪切形态.

开展分级加载蠕变试验，全程记录试样的轴向应变. 加载阶段的应变记录时间间隔为15 s，蠕变阶段的应变记录时间间隔为60 s，在试验结束后绘制全程时间-应变曲线.

在蠕变试验中，0°层理倾角试样内部的损伤会沿着某一倾斜角度发展，导致0°层理倾角炭质板岩发生剪切破坏. 随着切应力的堆积，45°倾角层理面逐渐开始破坏，最终形成贯通破坏面. 在损伤沿着某一倾斜角度发展最终形成贯通破坏面的同时，90°层理上产生的拉应力未能达到可以形成张拉破坏面的程度，只产生剪切破坏面，呈现出剪切破坏形态.

对0°、45°和90°层理倾角炭质板岩分别进行三轴蠕变试验，得到7种状态下的蠕变试验曲线，如图2~4所示. 图中，σ₁、σ₃分别为第一、三主应力，ε_t为轴向应变，t为蠕变时间. 可以看出，炭质板岩的蠕变特性明显受应力的影响，在较小应力下，炭质板岩只经历衰减蠕变阶段，随着应力的提高，炭质板岩开始发生稳态蠕变，在这2个过程之间存在明显的应力阈值σ_s. 通过总结规律发现，σ_s均出现在第2、3级应力之间. 在应力加载至破坏应力级后，炭质板岩发生加速蠕变.

炭质板岩的蠕变量随着应力的提高而逐渐增大，并在破坏应力级陡然增大，这与试样发生加速蠕变有关. 炭质板岩在蠕变过程中发生的总蠕变量随着围压的增大而减小，可见围压越大，炭质板岩越不容易发生蠕变，围压对炭质板岩的蠕变有抑制作用.

岩体介质在荷载作用下，内部产生不断扩展的微裂缝且部分变形不可恢复，称为弹塑性损伤. 为了考虑层理结构在炭质板岩蠕变过程中的影响，引入弹塑性层状损伤元件^[20-21]，如图5所示.

层状损伤元件的理论基础是连续损伤力学，其物理意义在于描述层状复合材料的损伤行为. 将上述弹塑性层理损伤元件看成由未损伤区域A₁及损伤区域A₂并联组成. 在某一应力下，该元件所产生的应变ε与A₁、A₂产生的应变ε₁、ε₂相等，其平衡关系为

(1)$ \sigma A = {\sigma _1}{A_1}. $

式中：σ为整个元件所受的应力，σ₁为未损伤区域所受的应力.

损伤前，损伤区域A₂未发生作用，该元件退化为弹性元件，其应力-应变关系满足胡克定律.

(2)$ {\sigma _1} = E{\varepsilon _1}. $

式中：E为弹性模量.

由此可知，层状损伤元件中损伤部位的损伤变量为D_c=A₂/A，则元件中应力之间的关系为

(3)$ \sigma = {\sigma _1}\left( {1 - {D_{\mathrm{c}}}} \right). $

根据未损伤元件的应力应变关系及元件各部分应变相等的原则，可得

(4)$ \sigma = E\varepsilon_1 \left( {1 - {D_{\mathrm{c}}}} \right). $

根据Kachanov蠕变损伤演化理论可知，定义的损伤有效应力如下：

(5)$ \tilde \sigma = \frac{\sigma }{{1 - D}}. $

假定在加速蠕变阶段，炭质板岩应变微元破坏服从Weibull分布，分布函数为

(6)$ F(t) = 1 - \varphi (t). $

$F(t)$是介于0~1.0的增函数，引入

(7)$ \varphi (t) = 1 - \frac{t}{{t_{\mathrm{c}}}}{\text{ }}{\text{.}} $

假设蠕变开始(t = t₀ = 0)时的损伤因子D_c = 0，蠕变破坏(t = t_c)时的损伤因子D_c = 1，则岩石经历加速蠕变后，定义层状损伤元件中的损伤演化方程为

(8)$ {D}_{{\mathrm{c}}}=F(t)=1-\left(1-{t}/{{t}_{{\mathrm{c}}}}\right)^{\frac{1}{c+1}}. $

式中：c为与层理倾角相关的损伤参数.

在岩石蠕变过程中，超过某个时间阈值后，弹塑性损伤元件中的损伤区域A₂将开启，定义该阈值为t_s(即岩体裂缝开始产生)，则损伤演化方程可以改进为

(9)$ {D}_{{\mathrm{c}}}=\left\{\begin{array}{l}0\text{，}t\leqslant {t}_{{\mathrm{s}}};\\ 1-\left(1-{t}/{{t}_{{\mathrm{c}}}}\right)^{\frac{1}{c+1}}\text{，}t > {t}_{{\mathrm{s}}}.\end{array}\right.$

综上，可得层状损伤元件的蠕变本构方程为

(10)$ \varepsilon (t)=\left\{\begin{array}{l}\dfrac{\sigma }{E}\text{，}t\leqslant {t}_{{\mathrm{s}}};\\ \dfrac{\sigma }{E}\left(1-\dfrac{t}{{t}_{{\mathrm{c}}}}\right)^{\frac{1}{c+1}}\text{，}t > {t}_{{\mathrm{s}}}.\end{array}\right. $

黏性损伤是指岩体介质在外力作用下产生裂隙并随时间不断演化的现象，不可恢复是其主要特点. 为了解决Nishihara模型无法描述加速蠕变的问题，引入非线性黏性损伤元件. 该元件的理论基础包括非线性力学和损伤力学，物理意义是描述材料在黏性流动过程中的损伤行为. 该元件在黏性元件中引入可描述岩体内部损伤演化规律的损伤变量，损伤变量的演化受到非线性应力和黏性流动的影响，导致岩石内部产生微裂纹、空穴的损伤形态，最终导致岩石材料的失效. 元件示意图如图6所示.

Kachanov蠕变损伤演化理论可以用于加速蠕变阶段，能够有效地模拟岩石蠕变的全过程，描述材料的逐渐衰变. 根据该理论定义岩体内部损伤变量为D_n（为净截面积与无损截面积之比，即连续度），损伤演化方程为

(11)$ {D}_{{\mathrm{n}}}=\left\{\begin{array}{l}0\text{，}t\leqslant {t}_{{\mathrm{n}}};\\ 1-{\mathrm{exp}}\;{[-(t-{t}_{{\mathrm{n}}})^{d}]}\text{，}t>{t}_{{\mathrm{n}}}.\end{array}\right. $

式中： t_n为岩石开始产生加速蠕变的时间.

当t$ \leqslant $t_n时，损伤还未发生，岩石处于稳态蠕变阶段，D_n = 0；当t＞t_n时，岩石进入加速蠕变阶段，D_n逐渐趋向于极限1，岩石开始发生蠕变破坏. 将上述损伤因子引入黏性元件中，组成可描述岩石非线性蠕变阶段的黏性损伤元件，蠕变本构关系为

(12)$ \varepsilon (t)=\frac{\sigma }{\eta }{{\mathrm{exp}}}\;{[(t-{t}_{{\mathrm{n}}})^{d}]}t. $

将上述两元件引入Nishihara中，建立考虑层理损伤的改进Nishihara非线性蠕变损伤本构模型，示意图如图7所示. 图中，η₁、η₂、η₃分别为衰减、稳态和加速阶段的蠕变速率，E₁为蠕变阶段开始前的瞬时弹性模量，E₂为衰减蠕变阶段的弹性模量.

当t$ \leqslant $t_s且t$ \leqslant $t_n时，改进Nishihara模型中的损伤变量均未发生作用，此时改进Nishihara模型退化为传统Nishihara模型，如图8所示. 蠕变本构方程如下所示：

(13)$ \varepsilon \left(t\right)=\left\{\begin{array}{l}\dfrac{\sigma }{{E}_{1}}+\dfrac{\sigma }{{E}_{2}}\left(1-{{\mathrm{exp}}}\;{\left(-\dfrac{{E}_{2}}{{\eta }_{2}}t\right)}\right),\;\sigma <{\sigma }_{{\mathrm{s}}};\\ \dfrac{\sigma }{{E}_{1}}+\dfrac{\sigma }{{E}_{2}}\left(1-{{\mathrm{exp}}}\;{\left(-\dfrac{{E}_{2}}{{\eta }_{1}}t\right)}\right)+\dfrac{\sigma -{\sigma }_{{\mathrm{s}}}}{{\eta }_{2}}t,\;\sigma \geqslant {\sigma }_{{\mathrm{s}}}.\end{array}\right. $

式中：$ {\sigma _{\mathrm{s}}} $为发生稳态蠕变的应力临界阈值.

当t > t_s且t > t_n时，模型中的损伤变量开始发挥作用，此时模型中各元件的状态方程如下所示：

(14)$ \left.\begin{array}{l}{\sigma }_{1}={E}_{1}\left(1-{t}/{{t}_{{\mathrm{c}}}}\right){\varepsilon }_{1},\\ {\sigma }_{2}={E}_{1}{\varepsilon }_{2}+{\eta }_{1}{{\dot \varepsilon }_{2}},\\ {\sigma }_{3}={\sigma }_{{\mathrm{s}}}+{\eta }_{2}{{\dot \varepsilon }_{3}},\\ {\sigma }_{4}={\eta }_{3}{{\mathrm{exp}}}\;{[-{\left(t-{t}_{{\mathrm{n}}}\right)}^{d}]}{{\dot \varepsilon }_{4}},\\ \sigma ={\sigma }_{1}={\sigma }_{2}={\sigma }_{3}={\sigma }_{4},\\ \varepsilon ={\varepsilon }_{1}+{\varepsilon }_{2}+{\varepsilon }_{3}+{\varepsilon }_{4}.\end{array}\right\} $

综上，根据 Boltzmann叠加原理，建立改进Nishihara非线性损伤蠕变模型的本构方程：

(15)$ \varepsilon \left(t\right) = \left\{ \begin{array}{l}\dfrac{\sigma }{{{E}}_{1}}+\dfrac{\sigma }{{E}_{2}}\left(1-{\rm{exp}}\;{\left(-\dfrac{{E}_{2}}{{\eta }_{1}}t\right)}\right),\;\sigma < {\sigma }_{\rm{s}};\\\dfrac{\sigma }{{E}_{1}}+\dfrac{\sigma }{{E}_{2}}\left(1-{\rm{exp}}\;{\left(-\dfrac{{E}_{2}}{{\eta }_{1}}t\right)}\right)+\dfrac{\sigma -{\sigma }_{\rm{s}}}{{\eta }_{2}}t,\text{ }\sigma \geqslant {\sigma }_{\rm{s}}\text{，}\\\qquad t\leqslant {t}_{\rm{s}},\;t\leqslant {t}_{\rm{n}};\\\dfrac{\sigma }{{E}_{1}} \left(1 - \dfrac{t}{{t}_{\rm{c}}}\right)^{\frac{1}{c+1}} + \dfrac{\sigma }{{E}_{2}} \left(1 - {\rm{exp}}\;{\left(-\dfrac{{{E}}_{2}}{{\eta }_{1}}t\right)}\right) + \dfrac{\sigma -{\sigma }_{\rm{s}}}{{\eta }_{2}}t+\\\qquad \dfrac{\sigma }{{\eta }_{3}}{\rm{exp}}\;{\left[\left(t-{t}_{\rm{n}}\right)^{d}\right]}t,\text{ }\sigma \geqslant {\sigma }_{\rm{s}}\text{，}t > {t}_{\rm{s}},\;t > {t}_{\rm{n}}.\end{array} \right.$

天然状态下的岩体基本处于三维应力状态下. 只单纯地建立岩体的一维蠕变本构方程无法满足实际需要，故将一维本构方程扩展到三维具有重要意义^[22-25]. 将一维蠕变本构方程扩展到三维应满足如下假定. 1）材料的体应变为弹性变形，且与时间因素无关，即在加载瞬间完成. 2）球应力张量不会引起岩石蠕变，只有偏应力张量才能产生蠕变变形. 3）在蠕变过程中，泊松比ν 不随时间推移而变化. 在黏塑性变形中，仅通过屈服强度σ_s判断岩石的屈服是不正确的，岩石的屈服由屈服函数确定，变形特性与流动法则有关，因此岩石蠕变本构方程中黏塑性变形阶段应包括屈服函数f与塑性势函数g. 此外，简单地将一维应力σ换成三维应力状态的偏应力S_ij不太合理.

根据弹塑性理论，将改进Nishihara模型一维蠕变本构方程扩展到三维，建立改进Nishihara模型的三维蠕变本构方程.

将应力σ_ij分解为球应力${{{{\delta}} }_{{{{ij}}}}} $σ_m和偏应力S_ij，相应地将应变${{{{\varepsilon}} }_{{{{ij}}}}} $分解为球应变δ_ijε_m和偏应变e_ij，其平衡关系为

(16)$ \left.\begin{gathered} {{{{\sigma}} }_{{{{ij}}}}}{\; =\; }{{{{S}}}_{{{{ij}}}}}+{{{{\delta}} }_{{{{ij}}}}}{{{{\sigma}} }_{{\rm{m}}}} ,\\ {{{{\varepsilon}} }_{{{{ij}}}}}{ \;=\; }{{{{e}}}_{{{{ij}}}}}+{{{{\delta}} }_{{{{ij}}}}}{{{{\varepsilon}} }_{{\rm{m}}}}. \\ \end{gathered}\right\} $

式中：δ_ij为Kronecker函数，

(17)$ {{{{\delta}} }_{{{{ij}}}}} = \left\{ \begin{gathered} 0,\;i \ne j ; \\ 1,\;i = j . \\ \end{gathered} \right. $

根据弹性力学理论可知，三维状态下弹性体的应力-应变关系可以表示为

(18)$\left. \begin{gathered} {{{{\sigma}} }_{{\rm{m}}}} = 3K{{{{\varepsilon}} }_{{\rm{m}}}}, \\ {{{{S}}}_{{{{ij}}}}} = 2G{{{{e}}}_{{{{ij}}}}}. \\ \end{gathered} \right\}$

式中：K和G分别为弹性体的体积模量和切变模量，

(19)$ \left.\begin{gathered} K = \frac{E}{{3\left(1 - 2\nu \right)}} ,\\ G = \frac{E}{{2\left(1+\nu \right)}}. \\ \end{gathered} \right\} $

式中：$\nu $为泊松比.

综上所述，弹性体的三维本构关系可以表示为

(20)$ {{{\varepsilon }}}_{{{{ij}}}}^{{\rm{e}}} = \frac{1}{{2{G_1}}}{{{{S}}}_{{{{ij}}}}}+\frac{1}{{3K}}{{{{\sigma}} }_{{\rm{m}}}}{{{{\delta}} }_{{{{ij}}}}}. $

由于球应力张量不会引起蠕变变形，黏弹性体的三维本构关系可以表示为

(21)$ {{{\varepsilon}} }_{{{{ij}}}}^{{{\rm{ve}}}} = \frac{1}{{2{G_2}}}{{{{S}}}_{{{{ij}}}}}\left(1 - {\exp{\left( - \frac{{{G_2}}}{{{\eta _1}}}t\right)}}\right). $

对于黏塑性体，应变的开启与应力阈值有关. 应力阈值函数f_s与屈服函数f的确定方法由蠕变试验确定. 黏塑性体应变可以表示为

(22)$ {{{\varepsilon }}}_{{{{ij}}}}^{{{\rm{vp}}}} = \frac{1}{{{\eta _2}}}\left[ {\phi \left( {\frac{{{f_{\mathrm{s}}}}}{{{f_0}}}} \right)} \right]\frac{{\partial {g_{\mathrm{s}}}}}{{\partial {{{{\sigma}} }_{{{{ij}}}}}}}t. $

式中：$ \phi $函数可以取为幂函数的形式，f₀为岩石屈服函数的初始参考值. 采用相关流动法则，有f_s = g_s（g_s为应力达到阈值时的塑性势函数），则黏塑性体的三维本构关系可以表示为

(23)$ {{{\varepsilon}} }_{{{{ij}}}}^{{{{{\mathrm{vp}}}}}} = \frac{1}{{{\eta _2}}}\left[ {{{\left( {\frac{{{f_{\mathrm{s}}}}}{{{f_0}}}} \right)}^{m}}} \right]\frac{{\partial {f_{\mathrm{s}}}}}{{\partial {{{{\sigma}} }_{{{{ij}}}}}}}t. $

式中：m为试验常数，通常取1. [·]为越阶函数，则式(23)可以展开为

(24)$ {{{\varepsilon}} }_{{{{ij}}}}^{{{{{\mathrm{vp}}}}}} = \left\{ \begin{gathered} 0,\;{f_{\mathrm{s}}} < 0; \\ \frac{1}{{{\eta _2}}}\frac{{{f_{\mathrm{s}}}}}{{{f_0}}}\frac{{\partial {f_{\mathrm{s}}}}}{{\partial {{{{\sigma}} }_{{{{ij}}}}}}}t,\;{f_{\mathrm{s}}} \geqslant 0 . \\ \end{gathered} \right. $

非线性黏性体的三维本构关系可以表示为

(25)$ {{{\varepsilon}} }_{{{{ij}}}}^{{\rm{p}}} = \frac{1}{{{\eta _3}}}{{{{S}}}_{{{{ij}}}}}{{\mathrm{exp}}{\left[{{\left(t - {t_{\mathrm{n}}}\right)}^{{d}}}\right]}}t. $

与改进Nishihara模型的一维本构方程进行类比，采用Boltzmann叠加原理，将上述元件的三维本构关系进行叠加，得到改进Nishihara模型本构方程的三维形式.

(26)$ {{\textit{ε}}}_{{\textit{ij}}}=\left\{\begin{array}{l}\dfrac{1}{2{G}_{1}}{{\textit{S}}}_{{\textit{ij}}}+\dfrac{1}{3K}{{\textit{σ}}}_{ {\rm{m}}}{{\textit{δ}}}_{{\textit{ij}}}+\dfrac{1}{2{G}_{2}}{{\textit{S}}}_{{\textit{ij}}}\left(1-{\rm{exp}}\;{\left(-\dfrac{{G}_{2}}{{\eta }_{1}}t\right)}\right) ,\\\qquad {f}_{\rm{s}} < 0;\\ \dfrac{1}{2{G}_{1}}{{\textit{S}}}_{{\textit{ij}}}+\dfrac{1}{3K}{{\textit{σ}}}_{ {\rm{m}}}{{\textit{δ}}}_{{\textit{ij}}}+\dfrac{1}{2{G}_{2}}{{\textit{S}}}_{{\textit{ij}}}\left(1-{\rm{exp}}\; {\left(-\dfrac{{G}_{2}}{{\eta }_{1}}t\right)}\right)+\\\qquad\dfrac{1}{{\eta }_{2}}\dfrac{{f}_{\rm{s}}}{{f}_{0}}\dfrac{\partial {f}_{\rm{s}}}{\partial {{\textit{σ}}}_{{\textit{ij}}}} t,\; {f}_{\rm{s}}\geqslant 0,\;t\leqslant {t}_{\rm{s}},\;t\leqslant {t}_{{\mathrm{n}}};\\ \left(\dfrac{1}{2{G}_{1}}{{\textit{S}}}_{{\textit{ij}}}+\dfrac{1}{3K}{{\textit{σ}}}_{ {\rm{m}}}{{\textit{δ}}}_{{\textit{ij}}}\right)\left(1-\dfrac{t}{{t}_{{\mathrm{c}}}}\right)^{\frac{1}{c+1}}+\\\qquad \dfrac{1}{2{G}_{2}}{{\textit{S}}}_{{\textit{ij}}}\left(1-{\rm{exp}}\;{\left(-\frac{{G}_{2}}{{\eta }_{1}}t\right)}\right)+\\ \qquad \dfrac{1}{{\eta }_{2}}\dfrac{{f}_{\rm{s}}}{{f}_{0}}\dfrac{\partial {f}_{\rm{s}}}{\partial {{\textit{σ}}}_{{\textit{ij}}}}t+\dfrac{1}{{\eta }_{3}}{{\textit{S}}}_{{\textit{ij}}}{\rm{exp}}\;{{\left[\left(t-{t}_{{\mathrm{n}}}\right)^{d}\right]}}t， \\\qquad {f}_{\rm{s}}\geqslant 0,\;t > {t}_{\rm{s}},\;t > {t}_{{\mathrm{n}}}.\end{array}\right. $

常规三轴压缩试验中有第二主应力σ₂ = σ₃，则式（26）中的一些参数可以转换为如下形式：

(27)$\left.\begin{aligned} &{{\textit{σ}}}_{\rm{m}}=\frac{1}{3}\left({\sigma }_{1}+2{\sigma }_{3}\right),\\ &{{S}}_{11}={\sigma }_{1}-{{\textit{σ}}}_{{{\mathrm{m}}}}=\frac{2}{3}\left({\sigma }_{1}-{\sigma }_{3}\right).\end{aligned}\right\}$

对于屈服函数，有

(28)$ \left.\begin{gathered} {f_{\mathrm{s}}} = \sqrt {{{{{J}}}_{2}}} - \frac{{{\sigma _{\mathrm{s}}}}}{{\sqrt 3 }}, \\ \sqrt {{{{{J}}}_{2}}} = \frac{{{\sigma _1} - {\sigma _3}}}{{\sqrt 3 }}. \\ \end{gathered} \right\}$

式中：J₂为应力偏量的第二不变量.

将式(25)、(26)代入式(24)，取岩石屈服函数的初始参考值f₀ = 1，则得到改进Nishihara模型的轴向蠕变方程：

(29)$ \begin{split} &{{\textit{ε}}}_{\boldsymbol{ij}}=\\&\left\{\begin{array}{l}\dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{1}} + \dfrac{{\sigma }_{1}+2{\sigma }_{3}}{9K} + \dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{2}}\left(1-{{\mathrm{exp}}}\;{\left(-\frac{{G}_{2}}{{\eta }_{1}}t\right)}\right),\\\qquad {f}_{{\mathrm{s}}} < 0;\\\dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{1}} + \dfrac{{\sigma }_{1} + 2{\sigma }_{3}}{9K} + \dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{2}}\left(1-{{\mathrm{exp}}}\;{\left(-\frac{{G}_{2}}{{\eta }_{1}}t\right)}\right)+\\\qquad \dfrac{{\sigma }_{1}-{\sigma }_{3}-{\sigma }_{\text{s}}}{3{\eta }_{2}}t,\text{ }{f}_{{\mathrm{s}}}\geqslant 0,t\leqslant {t}_{{\mathrm{s}}},t\leqslant {t}_{{\mathrm{n}}};\\ \left(\dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{1}}+\dfrac{{\sigma }_{1}+2{\sigma }_{3}}{9K}\right)\left(1-\dfrac{t}{{t}_{{\mathrm{c}}}}\right)^{\frac{1}{c+1}}+\\\qquad \dfrac{{\sigma }_{1}-{\sigma }_{3}}{3{G}_{2}}\left(1-{{\mathrm{exp}}}\;{\left(-\frac{{G}_{2}}{{\eta }_{1}}t\right)}\right)+\\ \qquad \dfrac{{\sigma }_{1} - {\sigma }_{3} - {\sigma }_{\text{s}}}{3{\eta }_{2}}t\;{+}\;\dfrac{2\left({\sigma }_{1}-{\sigma }_{3}\right)}{3{\eta }_{3}}{{\mathrm{exp}}}\;{\left[\left(t-{t}_{{\mathrm{n}}}\right)^{d}\right]}t,\text{ }\\\qquad {f}_{{\mathrm{s}}}\geqslant 0,t > {t}_{{\mathrm{s}}}, t > {t}_{{\mathrm{n}}}.\end{array} \right.\end{split}$

最小二乘法具有参数辨识度精度高、方法简便的优点，因而成为蠕变本构模型参数辨识最常用的方法. Levenberg-Marquardt算法是基于最小二乘法提出的非线性优化算法，Origin分析软件将该算法集成于非线性拟合功能中. 利用Origin分析软件中的Levenberg-Marquardt算法，对改进Nishihara模型进行拟合及参数辨识. 当t＜t_s且t$ \leqslant $t_n时，改进Nishihara模型将退化为传统Nishihara模型. 以不同围压下不同层理倾角的炭质板岩为例，分别采用传统和改进Nishihara模型，对破坏应力级蠕变曲线进行拟合验证.

通过图9的拟合对比发现，改进Nishihara模型对蠕变曲线的拟合度较高，可以描述炭质板岩蠕变的三阶段，解决了传统Nishihara模型无法描述加速蠕变阶段的缺点，验证了改进Nishihara模型描述炭质板岩的蠕变规律的适用性.

采用Levenberg-Marquardt算法，对上述蠕变试验中不同层理倾角炭质板岩的蠕变参数进行辨识，辨识结果如表2~4所示.

从图10、11可以看出，内部损伤参数d随着围压逐渐减小，当围压为5 MPa时取得最大值，当围压为15 MPa时取得最小值. 在围压为10 MPa的情况下，d随层理倾角先增大后减小，当层理倾角为45°时取得最大值，当层理倾角为0°时取得最小值，最大值与最小值的比值为2.5. 在围压为15 MPa的情况下，d随层理倾角逐渐增大，当层理倾角为90°时取得最大值，当层理倾角为0°时取得最小值，最大值与最小值的比值为3.1，均表现出明显的各向异性.

从图12、13可以看出，层状损伤参数c随围压先增大后减小，当围压为5 MPa时取得最小值，当围压约为10 MPa时取得最大值. 在不同围压下，c随层理倾角的变化规律均表现出先增大后减小的趋势，均在层理倾角约为45°时取得最大值，在层理倾角为0°时取得最小值，表现出明显的各向异性. c与层理倾角相关，c与层理倾角θ的变化规律满足二次函数，设c与θ满足关系式c(θ)=Aθ²+Bθ+C，使用该关系式对层状损伤参数c进行拟合，可得不同围压下c与θ的关系. 当围压为10 MPa时，c与θ的关系为

(30)$ c = - 6.42 \times {10^{ - 7}}{\theta ^2}+8.22 \times {10^{ - 5}}\theta +0.008\;6. $

当围压为15 MPa时，c与θ的关系为

(31)$ c = - 4.03 \times {10^{ - 6}}{\theta ^2}+3.10 \times {10^{ - 4}}\theta +0.003\;7. $

图14给出改进Nishihara模型的蠕变参数随围压的关系曲线. 图中，n为不同应力加载级次. 从图14(a)可以看出，不同应力水平下蠕变参数E₂呈现随围压增大而增大的趋势，在第4级应力水平下，E₂随围压增大呈现先增大后减小的趋势. E₂越大，则衰减蠕变阶段产生的应变越小. 从图14(b)可以看出，不同应力水平下η₁随围压呈逐渐增大的趋势，在第4级应力水平下η₁随围压先减小后增大，η₁越大，衰减蠕变阶段的蠕变速率越小. 从图14(c)可以看出，不同应力水平下η₂随围压总体上呈逐渐增大的趋势，η₂越大，稳态蠕变速率越小，由此可得不同应力水平下稳态蠕变速率随围压的增大而减小. 从图14(d)可知，η₃随围压逐渐减小，而η₃越大，加速蠕变阶段的蠕变速率越小，因此加速蠕变阶段的蠕变速率随围压逐渐增大.

为了方便研究，选取最后一级应力水平下的蠕变参数，探讨蠕变参数与层理倾角的关系. 从图15可以看出，不同围压下改进Nishihara模型的蠕变参数均随着层理倾角先减小后增大. E₂在层理倾角为45°时取到最小值，在层理倾角为90°时取到最大值，各向异性比分别达到4.04与4.20，可以评价为高级各向异性. E₂随层理倾角先减小后增大表征衰减蠕变阶段的弹性模量先减小后增大，则衰减蠕变阶段产生的应变先增大后减小. η₁在层理倾角为45°时取到最小值，在层理倾角为90°时取到最大值，最大值与最小值的比值分别达到3.42与2.82，各向异性明显. η₁随层理倾角先减小后增大，证明衰减蠕变阶段的蠕变速率减小速率随层理倾角先减小后增大. η₂在层理倾角为45°时取到最小值，在层理倾角为0°时取到最大值，最大值与最小值的比值分别达到11.95与2.62，各向异性明显. η₂随层理倾角先减小后增大，表示在应力达到或超过屈服应力后，稳态蠕变阶段的蠕变速率随层理倾角先增大后减小. η₃在层理倾角为45°时取到最小值，在层理倾角为0°时取到最大值，最大值与最小值的比值分别达到10.54与9.0，各向异性明显. η₃随层理倾角先减小后增大，表明加速蠕变阶段的蠕变速率随层理倾角先增大后减小.

图16给出不同损伤参数d对模型加速蠕变阶段的影响规律. 可以看出，不同的d对加速蠕变阶段产生的蠕变量及加速蠕变时间有着明显的控制作用，d越大，加速蠕变阶段产生的蠕变量越多，加速蠕变持续时间越短. 不同的d对加速蠕变阶段的蠕变速率和蠕变速率增量有着明显的影响，d越大，则加速蠕变速率和加速蠕变速率增量越大. 从试验中试样的破坏方式角度分析可知，随着d的增大，试样从延性破坏发展为脆性破坏. 本文中d随围压先增大后减小，但总体上呈现减小的趋势，即随着围压的增大，炭质板岩更趋向于延性破坏.

图17给出不同蠕变速率η₃对模型加速蠕变阶段的影响规律. 可以看出，η₃越大，加速蠕变阶段产生的蠕变量越少，η₃对加速蠕变开启时间无明显影响. 不同的η₃对加速蠕变阶段的蠕变速率和蠕变速率增量有着明显的影响，η₃越大，则加速蠕变速率和加速蠕变速率增量越小.

（1）炭质板岩的蠕变过程存在明显的阈值σ_s. 通过总结规律发现，σ_s均出现在第2、3级应力水平之间. 当应力小于σ_s时炭质板岩只发生衰减蠕变；当应力达到或大于σ_s时，开始发生稳态蠕变；当应力达到或超过破坏应力时，炭质板岩发生加速蠕变并发生蠕变破坏. 在蠕变过程中，层理倾角为0°的炭质板岩发生剪切破坏，45°层理倾角面上形成贯通破坏面，而90°层理倾角面上产生了剪切破坏面，呈现出剪切破坏形态.

（2）模型参数反演辨识结果表明，改进Nishihara模型能够较好地拟合各阶段的蠕变曲线，解决了传统Nishihara模型无法描述加速蠕变阶段的缺点，验证了改进Nishihara本构模型的正确性.

（3）基于参数辨识结果，探讨模型参数与围压和层理倾角的关系，得到损伤参数c与层理倾角的关系式. 归纳了损伤参数d与黏滞系数η₃对加速蠕变阶段的敏感程度. d越大，则加速蠕变阶段产生的蠕变量越多，加速蠕变速率和增量越大；η₃越大，则加速蠕变阶段产生的蠕变量越少，加速蠕变速率和增量越小，但η₃对加速蠕变开启时间无明显影响.