0 引言

岩体工程的长期稳定性受岩石蠕变特性影响^[1-2]，常因在外荷载长期作用下导致岩石失稳，引发工程事故，然而岩体工程失稳破坏必然经历加速蠕变阶段^[3-4]，经典的流变模型只能描述岩石的稳态蠕变^[5]，建立能够描述岩石蠕变过程中非线性特征的蠕变模型是研究的热点和难点。目前，国内外的学者们对于岩石的非线性流变本构模型已取得了一些研究成果。

自1869年麦克斯韦提出材料是黏弹性的，便开启了流变学的研究。1939年Griggs^[6]首次对灰岩、页岩和砂岩等岩石进行了蠕变试验研究，指出当应力水平为破坏荷载的12.5%~80%时，岩石具有明显的流变特性。在国内，陈宗基、孙钧等学者首先开始对岩体流变学进行系统性的研究。徐卫亚等^[7]建立了能充分反映岩石加速流变非线性特性的河海模型。蒋树等^[8]将非线性流变模型进行二次开发应用，实现了边坡渐进性失稳破坏的数值模拟研究。刘文博等^[9]结合最小耗能原理提出了一种考虑岩石损伤的蠕变模型。曹平等^[10]在研究节理岩体蠕变曲线时，引入了一种应变触发元件描述加速蠕变破坏阶段的非线性特征。徐平等^[11]采用柔性板法进行了岩体蠕变试验研究，得到了相应的蠕变模型及模型参数，为工程设计提供了参考依据。韩阳等^[12]基于损伤力学理论对Burgers模型进行了改进，改进后的模型能很好地描述加速蠕变阶段的非线性特征。

HU Bo等^[13]进行了硬岩蠕变试验，研究了硬岩蠕变三阶段的变形随时间的变化关系，提出了一种含分数阶导数阻尼器的非线性黏弹塑性模型。Ningyu Zhao等^[14]利用TFD-2000岩石三轴流变试验系统进行蠕变试验，研究了不同应力水平条件下的岩体的时变变形规律。Marco Emanuele Discenza等^[15]通过岩体蠕变试验研究了岩体的黏塑性，表明岩体的黏性是随加载应力和时间变化的，提出了等效黏度的概念。

针对经验蠕变模型物理意义不甚明确, 以及一些流变模型过于复杂，仅适用于某一类岩石的不足，本研究将从岩石蠕变试验中黏滞系数随不同应力水平和不同蠕变阶段变化的角度，建立包含应力、时间以及加速蠕变启动条件的黏性元件，结合经典蠕变模型，得到一种新的蠕变模型, 包含应力阈值和时间阈值双门槛条件，能准确表征蠕变曲线在衰减蠕变阶段和加速蠕变阶段非线性特征，适用性广，模型参数易求解。

1 建立非线性黏弹塑性流变模型 1.1 岩石蠕变试验应变-时间曲线分析

如图 1所示，完整的岩石流变曲线包含瞬时弹性变形、衰减蠕变、等速蠕变及加速蠕变4个阶段，是弹性、黏性、塑性、黏弹性和黏塑性等多种变形共存的一个复杂过程。

在低应力状态下，岩样在加载的瞬间产生瞬时弹性应变，随后蠕变速率衰减较快，并进入稳态蠕变阶段，可以用经典的流变模型饱依丁-汤姆逊体描述，但由于该模型由线性元件组合而成，对于衰减蠕变阶段的非线性特征拟合度精度不高，需进一步改进。

在高应力水平作用下，岩石发生稳定蠕变后可能会发生非稳定蠕变变形。目前，主要采用应力阈值(长期强度)开关作为加速流变特征的触发途径。然而，通过对不同岩性的岩石在不同应力水平下的蠕变曲线进行分析，当加载应力接近或者大于长期强度时，岩石发生黏塑性破坏类型不同，从稳定蠕变阶段进入非稳定蠕变阶段的起始时间及持续时间不同，相应的应变-时间曲线在加速蠕变启动时刻表现为偏塑性突变型和偏黏塑性圆滑型，理论上当蠕变速率不为0时岩样最终会进入加速蠕变阶段发生黏塑性破坏，只是流变历时不同而己，因此以某个应力阈值作为加速蠕变启动条件难以准确描述进入加速阶段的时间和持续时间。

1.2 模型建立 1.2.1 初始瞬时弹性变形和衰减蠕变阶段

孙钧等^[16]通过对大量的岩石流变试验曲线进行分析，得到在低应力水平下，岩石发生初始瞬时弹性变形后进入减速蠕变阶段，黏性系数随着加载应力的增大而减小，随着时间的增加而增加，此阶段的黏性系数符合关于应力和时间的复合幂函数关系，建立新的黏性元件N₁，满足如下关系式：

(1)

式中，σ_v，t_v为标准参考应力和参考时间，取值为1;n为与应力水平相关的流变指数。

将N₁替换饱依丁-汤姆逊体中的牛顿黏壶，则有：

(2)

对于胡克体：

(3)

改进后的饱依丁-汤姆逊体本构方程为：

(4)

在恒定应力σ₀作用下，解上述微分方程，得到改进后的饱依丁-汤姆逊体蠕变方程：

(5)

当t=0时，黏性体相当于刚性体，并未产生瞬时应变：

(6)

将式(6)代入式(5)得：

(7)

1.2.2 等速蠕变及加速蠕变阶段

当σ₀≥σ_s时，不同岩性的岩石在不同应力水平条件下进入加速蠕变阶段的时间不同，但只要时间足够长，终将进入加速蠕变阶段，此阶段的黏性系数随着应力和时间的增加而减小^[16]，符合关于应力和时间的复合指数函数关系。同时引入加速蠕变时间阈值开关t^*，当t＜t^*时，处于稳态蠕变阶段，黏壶的黏性系数为常数；当t≥t^*时，进入加速蠕变阶段，黏性系数随应力、时间变化，应变随时间加速增加，由此建立非线性黏性元件N₂：

(8)

式中，m为与岩石岩性相关的材料参数。将N₂与塑性元件并联组成新的非线性黏塑性元件，其本构方程为：

(9)

当σ₀≥σ_s时，解上述微分方程，得改进后的非线性黏塑性元件蠕变方程为：

(10)

1.2.3 一维非线性黏弹塑性流变模型

如图 2所示，首先建立非线性部分的微分方程，然后通过积分得出该部分的蠕变方程，再通过叠加法得到新的非线性黏弹塑性流变模型的蠕变方程^[17]。

一维蠕变方程为：

(11)

2 理论验证 2.1 非线性黏弹塑性流变模型的稳定性分析

定性分析非线性黏弹塑性流变模型的稳定性，令：

(12)

式(5)简化为：

(13)

对式(14)求解关于时间变量的一阶导数和二阶导数，得：

(14)

(15)

由于A, B, C均大于0，则式(14)恒大于0，式(15)恒小于0，即式(5)为随时间递增的凸函数，可以描述衰减蠕变阶段变形值递增、蠕变速率递减的特征。

同理对式(10)关于时间的一阶导数和二阶导数进行分析，得：当t＜t^*时，其一阶导数为常数，属于等速蠕变；当t≥t^*时，式(10)为随时间递增的凹函数，能够描述加速蠕变阶段变形值及蠕变速率递增的特征。

2.2 模型参数的敏感性分析

如图 3、图 4所示，模型参数n主要影响衰减蠕变阶段的蠕变速率，取值范围为0＜n＜1，n越大蠕变速率衰减得越快；模型参数m＞1，主要影响加速蠕变阶段的变形值和蠕变速率，m越大蠕变变形值越大，蠕变速率增量越大。通过调整n，m的值可以很好地刻画蠕变全过程中衰减蠕变阶段和加速蠕变阶段的非线性特征。

3 模型参数辨识

(1) E₁，E₂，η₁，n的确定

当t=0时，对于初始瞬时弹性变形值和减速蠕变阶段：

(16)

当t→∞时：

(17)

联立式(16)、式(17)，得到E₁，E₂，再通过衰减蠕变阶段的试验值时即可求得n，η₁。

(2) η₂的确定

当σ₀≥σ_s且t＜t^*时，为等速蠕变阶段，应变-时间曲线的斜率k：

(18)

结合牛顿黏性体的蠕变方程得：

(19)

当σ₀≥σ_s且t≥t^*时，根据上述方法，采用分段拟合的形式，确定E₁，E₂，n，η₁，η₂，然后代入式(11)，根据加速蠕变阶段试验值即可求得模型参数m。

上述模型参数均可在Origin软件中通过建立自定义函数，采用最小二乘法，根据试验得到的x组(t_i，ε_i)试验值，进行非线性拟合确定。

4 加速蠕变启动条件

不同岩性的岩石具有不同的蠕变特性，例如软岩在应力水平接近长期强度值时，以几乎恒定的蠕变速率产生蠕变变形，而硬岩在稳态蠕变阶段应变基本不发生变化。为准确描述岩石蠕变破坏的非线性特征，确定加速蠕变阶段的起始时间尤为重要，且加速蠕变启动时间的选择会影响到各流变力学参数的计算结果。

(1) 根据蠕变速率突变点确定

当应力水平大于加速蠕变应力阈值时，岩石由稳态蠕变阶段进入加速蠕变阶段，加速蠕变启动时间可以根据蠕变速率-时间曲线的突变点直接确定。例如文献[18]中岩石在8.78 MPa应力水平作用下，第15 h时，蠕变速率突然增大到1.68×10^-3/h，且呈加速持续增长的趋势。因此，判定在该级应力水平作用下加速蠕变启动时间t^*=15 h。

(2) 岩石蠕变破坏时间与应力水平的关系

一些学者^{[16, 19]}在研究岩石长期强度值时，通过对一定数量的岩石蠕变试验结果进行研究，发现快剪、快压试验得到的瞬时强度值最高，破坏荷载随流变时间(t)及加速流变历时(Δt)的增加呈减小的趋势，符合指数函数关系，如图 5所示：

(20)

(21)

联立式(20)、式(21)得加速蠕变破坏的起始时间为：

(22)

式(22)为岩石发生蠕变破坏的时间与应力水平的一般关系式，根据岩石蠕变试验拟合得到参数α，β值，通过式(24)可以推算出不同应力水平下岩石蠕变破坏的起始时间及流变历时，对于岩体工程稳定性预警、预测研究意义重大。

5 三维蠕变方程

实际岩体工程中岩石处于三向应力状态，对于三维应力状态下的弹性体，材料内部任一点的应力张量σ_ij可以分为球应力张量σ_m和偏应力张量S_ij；应变张量ε_ij分为球应变张量ε_m和偏应变张量e_ij，由：

(23)

可得：

(24)

式中δ_ij为克罗内克函数：

(25)

同理得到：

(26)

对于常规三轴压缩蠕变试验：

(27)

结合式(24)、(27)得：

(28)

对于三维应力状态下的胡克体：

(29)

式中，K为体积模量; G为剪切模量。结合一维蠕变方程，根据文献[20]采用类比法，得到弹性体三维蠕变方程：

(30)

黏性体的三维蠕变方程为：

(31)

最后，采用叠加法得到非线性黏弹塑性流变模型的三维本构模型：

(1) 当σ₀＜σ_s时：

(32)

(2) 当σ₀≥σ_s且t＜t^*时：

(33)

(3) 当σ₀≥σ_s且t≥t^*时：

(34)

(35)

式中F₀为屈服函数初始值，一般取1^[21]。

由于Mohr-Coulomb屈服准则不能考虑中间主应力对岩石屈服的影响，Von-Mises屈服准则忽视了球应力对岩石蠕变特性的影响，笔者采用Drucker-Prager屈服准则，其表达式为：

(36)

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

(37)

6 试验验证

笔者采用文献[18]及文献[22]中流变试验数据对上述建立的非线性黏弹塑性蠕变模型的适用性及正确性进行验证。

文献[18]为现场大尺寸压缩流变试验，试样为节理裂隙较为发育的微风化砂岩，试样尺寸为730 mm×730 mm×780 mm，采用逐级加载方式。根据试验结果，在应力水平为8.78 MPa时，取t^*=15 h，σ_∞=7.50 MPa。

文献[22]为室内小尺寸剪切流变试验，试样为块状构造斑状结构的灰白色花岗岩圆柱体，尺寸为50 mm×100 mm，采用逐级加载方式。根据试验结果，在应力水平为17.27 MPa时，t^*=70 h，σ_∞=15.57 MPa。

结果如图 6、图 7、表 1、表 2所示，拟合效果很好，在最高级应力水平作用下，文献[18]和文献[22]的蠕变曲线拟合度分别达到了0.997 8和0.978 4，这说明基于岩石蠕变试验中黏性系数在不同应力水平和不同蠕变阶段中的变化规律建立的岩石非线性黏弹塑性蠕变模型，能反映岩石蠕变全过程曲线特征，解决了传统流变模型在衰减蠕变阶段拟合精度不高、无法描述加速蠕变阶段非线性特征及适用性不高的缺陷，且流变破坏物理意义明确。

7 结论

(1) 基于岩石蠕变试验中黏性系数在不同应力水平和不同蠕变阶段中的变化规律建立了参数求解简单、流变破坏物理意义明确的五元件非线性黏弹塑性蠕变模型。

(2) 通过模型稳定性和参数敏感性分析，并结合现场大尺寸软岩压缩流变试验和室内小尺寸饱水花岗岩剪切流变试验，从理论和试验两方面验证了所建模型的正确性与适用性。

(3) 当应力水平接近或者大于岩体的长期强度时，岩石并不会立即进入加速蠕变阶段发生黏塑性破坏，采用应力阈值和时间阈值双开关能精确刻画岩石加速蠕变阶段的非线性特征，对于岩体工程中因岩石流变特性失稳预测方面的研究具有重大意义。