考虑黏弹塑性应变分离的岩石复合蠕变模型研究

张亮亮，王晓健,2

(1.安徽理工大学土木建筑学院，安徽淮南，232001；2.安徽理工大学省部共建深部煤矿采动响应与灾害防控国家重点实验室，安徽淮南，232001)

岩石蠕变力学特性是岩石固有属性之一，尤其是深部赋存条件复杂的软岩，受高地应力、高温和高水压影响，其节理裂隙往往发育较好，宏观上表现为显著的蠕变变形，该变形特性对服务期限较长结构的安全性和长期稳定性至关重要[1-2]，因而，在进行深部地下结构(地下硐室群、核废料储存室及水电站地基等)设计时，不仅要考虑荷载作用下的瞬时变形，还要着重考虑结构随时间的蠕变变形。目前，关于岩石蠕变特性研究主要包括室内单、三轴分级加载压缩蠕变试验研究以及根据试验结果建立相应的元件(弹性元件、黏性元件和塑性元件)组合模型理论研究。如韩阳等[3]提出了一种非线性黏性元件，将其替换Burgers 模型中的2个线性黏性元件，建立了一种可描述岩石三阶段蠕变特性的非定常参数Burgers 模型；陈卫忠等[4]基于泥岩现场大型真三轴蠕变试验结果，提出了一种能够真实反映深部软岩不同应力下的流变特征的非线性经验幂函数型蠕变模型；HU 等[5]基于分数阶理论和损伤理论，提出了一种非线性损伤蠕变模型，该模型可描述硬岩加速蠕变阶段；刘开云等[6]在Maxwell 模型基础上串联1 个非线性黏塑性体，建立了可以描述岩石非稳态蠕变的参数非定常非线性黏弹塑性蠕变本构模型；刘东燕等[7]基于不同围压作用下重庆地区砂岩蠕变结果，将黏滞系数随时间按幂函数衰减的经验模型与Kelvin模型串联，建立了可描述不同围压作用下加速蠕变特征的改进黏塑性模型；ROBERTS 等[8]针对盐岩进行了三轴压缩、拉伸和循环拉压蠕变试验，研究了不同加载方式及应力路径对盐岩蠕变及扩容特性的影响；SINGH 等[9]基于Maxwell模型采用声发射技术获得了盐岩单轴压缩过程中的弹性模量和黏性系数，并根据该参数预测盐岩的单轴蠕变行为；BOUKHAROV 等[10]分析了脆性岩石流变曲线中的3个典型过程，并将应力触发的虎克体和应变触发的非线性黏壶元件串联用于描述非稳态蠕变过程。通过对上述研究成果及目前岩石蠕变力学发展过程的分析总结，本文作者认为仍存在以下2个方面不足：

1)缺乏岩石循环加卸载蠕变试验研究。在地下硐室开挖过程中，围岩原始应力平衡状态被打破，围岩由于卸载应力得到释放产生较大变形，随硐室进一步开挖，扰动围岩继续受开挖过程中爆破及机械振动等荷载，此过程相当于加载过程；再如大型水电站基础，在蓄水过程中，水压通过大坝将部分压力传递给坝基，而在放水过程中这部分压力又会随水位的降低逐渐释放。实际岩体工程多数受循环加卸载作用，而目前针对不同类型岩石的单、三轴分级加载试验研究成果较多，但关于循环加卸载作用下岩石变形响应的试验报道却相对较少。李娜等[11-13]进行了该方面的研究，但由于岩石种类繁多、地质条件复杂，所取得的成果应用受限，因此，在不同应力、不同围压及其他复杂应力路径下岩石循环加卸载蠕变试验尚需进一步研究。

2)岩石蠕变模型建立尚需完善。目前岩石蠕变模型理论最成熟、应用最广泛的是元件组合模型，该模型将弹性元件、黏性元件和塑性元件通过不同的串并联方式进行组合来反映岩石在不同应力下的蠕变特性。但研究发现经典元件组合模型(Bingham，Burgers和Nishihara模型等)只能描述岩石在低应力状态下发生的稳态蠕变，而不能描述岩石在高应力状态下发生的非稳态蠕变，尤其是非线性特征尤为明显的加速蠕变阶段。因此，目前大多数研究重点是基于非线性损伤理论，采用非线性黏性元件代替线性元件来对传统模型进行改进和完善，且多数研究都是建立加载条件下岩石的蠕变模型，针对岩石卸载条件下的蠕变特征研究及相应条件下的蠕变本构方程研究较少。同时，不少研究都采用弹性体来描述岩石在加载过程中产生的瞬时应变，认为该应变都为瞬时弹性应变，忽略了瞬时塑性应变，赵延林等[14]的研究表明瞬时塑性应变约占瞬时应变的一半且占比随应力的提高而提高，因此，仅考虑瞬时弹性应变忽视瞬时塑性应变所建立的蠕变模型不够严谨。

本文作者以文献[15]中石灰岩循环加卸载蠕变试验数据为基础，采用瞬时弹性体(IEB)、瞬时塑性体(IPB)、非线性黏弹性体(NVEB)和黏塑性体(VPB)分别描述瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变，并将四者串联建立复合蠕变模型。最后，根据非线性岩石流变理论推导岩石单、三轴蠕变及卸载方程，结合试验曲线与模型理论曲线对比结果验证模型的合理性。

1 石灰岩循环加卸载蠕变试验数据分析

赵延林等[14]采用MTS815岩石试验系统对石灰岩进行了不同围压下的循环加卸载蠕变试验，石灰岩基本力学参数及蠕变试验中偏应力分别见表1和表2。

表1 石灰岩基本力学参数Table 1 Basic mechanical parameters of limestone

围压分别为4 MPa 和8 MPa 的L2 和L7 试样的石灰试样的循环加卸载蠕变曲线如图1所示，图1中，σ3为第3 主应力。由图1可见：在前2 级应力下，蠕变曲线由瞬时应变通过衰减蠕变进入等速蠕变后，蠕变速率基本为0，卸载后瞬时恢复应变比加载瞬时应变小，说明加载过程中不仅有瞬时弹性应变，还有瞬时塑性应变；后3级应力下等速蠕变曲线斜率明显大于0，且卸载时的永久残余应变大于加载瞬时塑性应变，这表明在蠕变过程中，不仅有黏弹性应变还存在黏塑性应变。

图1 不同围压下石灰岩试样的循环加卸载蠕变曲线Fig.1 Cyclic loading unloading creep curves of limestone samples at different confining pressures

典型的岩石加卸载蠕变曲线示意图[16]如图2所示。

图2 岩石加卸载应变分离示意图Fig.2 Schematic diagram of strain separation of rock loading and unloading

根据图2及上述试验结果分析，可以将加载过程中总应变ε分为与时间无关的瞬时应变εm和与时间相关的蠕变应变εv，其中瞬时应变又可分为瞬时弹性应变εme和瞬时塑性应变εmp，蠕变应变又可分为可随时间逐渐恢复的黏弹性应变εve和不可恢复的黏塑性应变εvp，即

式中：εrp为残余塑性应变。

基于岩石蠕变与卸载时应变恢复特征，将蠕变试验曲线根据式(1)进行分解，得到不同应力下各应变，如表3所示。

表3 岩石瞬时弹、塑性应变和黏弹、塑性应变试验结果Table 3 Testing results of instantaneous elasto-plastic strain and viscoelasto-plastic strain

由表3可以看出：瞬时弹性应变、瞬时塑性应变、黏弹性应变和黏塑性应变都随应力的提高而提高。以试样L2 为例，5 种应力水平下瞬时塑性应变占瞬时应变的比例分别为33.11%，32.16%，31.28%，33.61%和34.50%，由于所占比例较大，因此，在建立描述瞬时应变的蠕变模型时，不能为了简化模型而忽略瞬时塑性应变，应将瞬时弹性应变和瞬时塑性应变分开考虑，分别建立相应的模型。此外，随偏应力逐级增大，5种应力下黏塑性应变占蠕变应变的比例分别为21.67%，28.70%，38.73%，42.74%和43.64%，这表明随偏应力增加，岩石不可恢复的黏塑性应变不断增加且逐渐累积储存于试样内部，当试验条件不能有效扼制住试样的黏塑性应变发展时就会发生破坏。在低应力下，岩石黏塑性应变所占比例较小，岩石的蠕变变形更大程度上表现为黏弹性变形；而在高应力下，岩石的蠕变变形表现为黏弹性和黏塑性共存的状态。

2 复合蠕变模型的建立

2.1 瞬时弹性体与瞬时塑性体

根据表1试验处理结果，绘制偏应力与瞬时弹性应变和瞬时塑性应变的关系曲线，结果如图3所示。

图3 瞬时弹、塑性应变与偏应力的关系曲线Fig.3 Relation curves among instantaneous elastic,plastic strain and deviator stress

由于应力加载阶段时间比后期应力维持加载时间短，因此，可认为瞬时弹性应变和瞬时塑性应变瞬间完成，两者量值只与应力有关而与时间无关。由图3可知：瞬时弹性应变、塑性应变与应力近似为线性增长关系，因此，可将建立瞬时弹性体与瞬时塑性体串联来描述，力学模型见图4。

图4 瞬时弹性体和塑性体模型Fig.4 Models of instantaneous elastic and plastic body

蠕变方程如下式所示：

式中：σ为应力；Ee为瞬时弹性体和塑性体的瞬时弹性模量；Ep为瞬时塑性体的瞬时弹性模量。

2.2 非线性黏弹性体

蠕变曲线由瞬时应变进入等速蠕变之前会经历衰减蠕变阶段，该阶段蠕变曲线斜率会逐渐减小，表现出非常明显的非线性特征，且衰减过程、特征及岩石进入等速蠕变的时间都与岩石应力密切相关。随蠕变试验进行，在常应力条件下岩石内部弱胶结物会逐渐发生破坏，宏观表现为蠕变应变随加载时间逐渐增加，在该阶段，黏性参数会随时间逐渐变大并在岩石进入等速蠕变阶段后逐渐趋于稳定值[17]，说明蠕变已达稳定状态。传统Kelvin模型虽然能够较好地描述特定应力条件下蠕变衰减特征，但由于其黏性参数始终为不随时间变化的常数，因此，难以精确地反映不同应力下蠕变非线性特征。本文考虑Kelvin 模型中黏性系数的时间相关性，假定该系数在蠕变过程中与时间满足幂函数关系[18]，从而建立非线性黏弹性模型，如图5所示。

图5 非线性黏弹性体模型Fig.5 Model of nonlinear viscoelastic body

该模型微分型本构方程为

式中：εve为黏弹性应变；Eve和ηve分别为非线性黏弹性体的弹性模量和黏性系数；λ为常数；t为时间。

解式(3)微分方程得到非线性黏弹性体的蠕变方程为

图6所示为不同λ对应的非线性黏弹性体蠕变曲线。由图6可知：λ越大，岩石由衰减蠕变进入等速蠕变所需的时间越短，衰减蠕变阶段蠕变斜率衰减的速度越快，说明该模型能够反映不同蠕变时间和不同应力状态下岩石衰减蠕变曲线特征，适用性较传统Kelvin模型更广泛。

图6 不同λ岩石蠕变曲线Fig.6 Creep curves of rock at different λ

假设在t0时刻卸载，式(3)可表示为

解式(5)得

式中：A为积分常数。

当t=t0时，代入式(6)得到积分常数A为

把式(7)代入式(6)得到卸载阶段蠕变εul为

2.3 黏塑性体

当试验应力介于某一临界应力σcritical和长期强度σs之间时，岩石等速蠕变阶段蠕变曲线表现出斜率大于0且随时间呈近似线性增长关系。该阶段岩石试样受长期蠕变特性及试样非均质性影响，在应力长时间作用下试样内部一些软弱结构处发生了局部破坏，此破坏由卸载后永久残余应变大于加载瞬时塑性应变所致。根据岩石流变理论，采用黏塑性体可描述该蠕变阶段岩石黏塑性应变随时间的变化关系，其力学模型如图7所示。

图7 黏塑性体模型Fig.7 Model of viscoplastic body

黏塑性体的蠕变方程为

式中：ηvp为黏塑性体的黏性系数。

2.4 复合蠕变模型

至此，建立了瞬时弹性体、瞬时塑性体、非线性黏弹性体和黏塑性体，分别描述岩石瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变。根据叠加原理，将上述4 个模型串联，建立可以同时描述加载过程中的上述4种应变、完全卸载时瞬时恢复的瞬时弹性应变和随时间逐渐恢复的黏弹性应变的复合蠕变模型，该复合蠕变力学模型如图8所示。

图8 复合蠕变模型Fig.8 Model of composite creep

当σ＜σcritical时，岩石总应变由瞬时弹性应变、瞬时塑性应变和非线性黏弹性应变组成，其本构方程为

当σ≥σcritical时，岩石总应变由瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变组成，其本构方程为

根据拉普拉斯变换和逆变换，得到复合蠕变模型的蠕变方程为

由于卸载过程中瞬时塑性应变和黏塑性应变不可恢复，故复合蠕变模型的卸载方程为

3 三维蠕变本构方程的建立

由于实际岩体往往处于复杂的三维应力状态下，且本文引用数据为三轴蠕变数据，因此，为反映实际岩石受力情况以及更好地将模型理论结果与试验结果进行对比，需建立岩石在三维应力状态下复合蠕变模型的三维蠕变本构关系。假设岩石在三维应力状态下总应变为εij，瞬时弹、塑性应变分别为和，非线性黏弹性应变和黏塑性应变分别为和，根据叠加原理，有

根据岩石弹塑性理论得到岩石三维应力状态下瞬时弹性应变和瞬时塑性应变可表示为[19]

式中：Ge和Ke分别为瞬时弹性体的剪切模量和体积模量；Gp为瞬时塑性体的剪切模量；Sij为偏应力张量；σm为平均应力；δij为Kronecker张量。

主应力空间下偏应力张量可表示为

非线性黏弹性体的三维蠕变方程可以参考一维应力状态下的蠕变方程：

式中：Gve为非线性黏弹性体的剪切模量。

对于黏塑性体，其三维本构方程不仅与岩石屈服函数F和塑性势函数Q相关，而且与流动法则相关，并不是简单地将偏应力张量替换一维蠕变方程中的应力[20]。根据流变力学理论得到黏塑性体三维蠕变方程为

式中：J2为应力偏量第2不变量；σ1为第1主应力。

在常规三轴压缩蠕变试验中，第2主应力σ2=σ3，则

结合式(14)～(20)，得到复合蠕变模型三维蠕变方程为

卸载过程中的三维蠕变方程为

4 模型验证

以试样L7 为例，根据式(21)得到复合蠕变模型各级偏应力下理论曲线与试验曲线的对比结果如图9所示。由图9可知：模型理论曲线与试验曲线吻合程度非常高，相关系数均在0.96 以上，尤其在卸载过程，两者基本一致，由此说明本文建立的复合蠕变模型不仅能够描述石灰岩加载过程中非线性蠕变特征，而且能描述卸载过程中应变随时间的变化规律，证明了该模型的合理性。

图9 试样L7在不同偏应力水平下模型理论曲线与试验曲线的对比结果Fig.9 Comparisons of theoretical curves and experimental curves of L7 specimen at different deviator stresses

为进一步验证复合蠕变模型的适用性及可行性，以试样L2 试验曲线为基础，采用复合蠕变模型理论曲线和广义Kelvin 蠕变模型理论曲线分别与试验曲线进行对比，结果如图10所示。广义Kelvin蠕变模型三维轴向蠕变方程为

式中：G′e和K′e分别为广义Kelvin 模型中弹性体的剪切模量和体积模量；G′ve和η′ve分别为广义Kelvin模型中黏弹性体的剪切模量和黏性系数。

相应的三维卸载蠕变方程为

由图10可知，广义Kelvin 蠕变模型和复合蠕变模型都能很好地描述加载过程中石灰岩的蠕变特征，但在衰减蠕变阶段，复合蠕变模型描述效果比广义Kelvin 蠕变模型的效果好。同时，复合蠕变模型卸载蠕变曲线基本上与试验曲线保持一致，而广义Kelvin 模型曲线明显比试验曲线低，这是由于广义Kelvin 蠕变模型不能将瞬时塑性应变从瞬时应变中分离，也不能将黏塑性应变从蠕变应变中分离，而是认为加载过程中产生的瞬时应变在卸载的一瞬间完全恢复，从而使得模型理论曲线与试验曲线产生较大误差。而复合蠕变模型综合考虑了瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变，认为卸载过程中瞬时弹性应变瞬间恢复，非线性黏弹性应变随蠕变时间逐渐恢复，而瞬时塑性应变和黏塑性应变为永久残余应变，因此，复合蠕变模型卸载过程中模型曲线与试验曲线高度吻合，相关系数均大于0.968，对比结果进一步验证本文建立的能够同时描述石灰岩分级加卸载蠕变特征的复合蠕变模型是合理、可行的。

图10 试样L2在不同偏应力水平下试验曲线、复合蠕变模型拟合曲线和广义Kelvin模型拟合曲线对比结果Fig.10 Comparison results of test curve,composite creep model fitting curve and generalized Kelvin model fitting curve of sample L2 at different deviator stresses

石灰岩试样L2和L7在不同偏应力下复合蠕变模型参数见表4。

表4 各级偏应力下复合蠕变模型参数Table 4 Parameters of composite creep model under different deviator stresses

5 结论

1)本文根据石灰岩循环加卸载蠕变曲线特征，将瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变从总应变中分离，采用瞬时弹性体、瞬时塑性体、非线性黏弹性体和黏塑性体对上述4 种应变分别进行描述，再将四者串联，建立了一种能够同时描述石灰岩加、卸载过程蠕变特征的复合蠕变模型。

2)根据叠加原理和非线性流变理论推导了岩石在一维、三维应力状态下的加卸载蠕变本构方程。石灰岩不同偏应力下加卸载蠕变试验曲线和复合蠕变模型理论曲线的对比结果表明两者高度吻合，该模型不仅能够描述石灰岩加载过程中非线性蠕变特征，而且能描述卸载过程中应变随时间的变化规律。

3)用复合蠕变模型理论曲线和广义Kelvin 蠕变模型理论曲线分别与试验曲线进行对比，结果表明广义Kelvin 蠕变模型和复合蠕变模型都能很好地描述加载过程中石灰岩的蠕变特征，但在衰减蠕变阶段及卸载蠕变阶段，复合蠕变模型描述效果要优于广义Kelvin 蠕变模型。复合蠕变模型综合考虑了瞬时弹性应变、瞬时塑性应变、非线性黏弹性应变和黏塑性应变，在卸载过程中模型曲线与试验曲线高度吻合，相关系数均大于0.968，对比结果验证了本文建立的复合蠕变模型是合理、可行的。