基于能耗理论的岩石三维蠕变本构模型及临界分段

2022-01-23卢功臣祝荃芃周林林

长江科学院院报 2022年1期

卢功臣，祝荃芃，周林林

(1.陕西省水利电力勘测设计研究院，陕西咸阳 712000； 2.江西师范高等专科学校，江西鹰潭 335000；3.辽宁工程技术大学土木工程学院，辽宁阜新 123000)

1 研究背景

岩石流变主要包括蠕变、松弛和弹性后效等现象[1]。蠕变是指材料在应力恒定的情况下，变形随着时间缓慢累积的现象，因其与工程岩体紧密关联而引起学术界广泛关注[2]。其中加速蠕变阶段在岩石蠕变中最为关键，如何有效预测加速蠕变行为是流变力学中的核心内容[3]。

现国内外关于岩石蠕变行为(尤其是加速蠕变)预测已有较多成果，宋飞[4]提出一个与时间和应力相关的可描述岩石加速蠕变阶段的非线性SP元件；Shibata等[5]通过研究蠕变应变率与流变参数之间的关系，建立可较好描述凝灰岩蠕变行为的本构模型；Sterpi等[6]构建同时考虑黏弹性和黏塑性蠕变的损伤模型；宋勇军等[7]基于硬化和损伤效应建立岩石蠕变损伤模型；周宏伟等[8]引入分数阶微积分，建立可较好描述岩盐加速蠕变阶段的流变本构模型；张泷[9]在Rice内变量热力学框架，通过给定比余能推导出岩石黏弹塑性蠕变本构模型；Ping等[10]改进传统Burgers模型并与连续损伤理论结合，引入非线性元件，建立适用于高应力区岩石的蠕变损伤模型。

目前对岩石加速蠕变阶段的判别和模拟还存在一定争议，多数模型都是在传统元件模型的基础上添加非线性元件，来描述岩石加速蠕变，虽然部分模型改进后能反映加速蠕变量值变化，但从蠕变机制而言无法合理解释。实际上，岩石蠕变过程中与外界能量不断交换，岩石材料能量耗散导致损伤劣化效应增强，以能量的角度能从某种意义上较好地反映岩石蠕变机制。谢和平等[11]系统分析了岩石变形破坏中的能量耗散；Sih[12]通过应变能理论解释了微裂纹扩展引起宏观应力应变过程；张萍等[13]探讨了页岩变形破坏过程中应变能的积聚和耗散，为岩石变形破坏过程中能量耗散分析提供有益参考。由此，本文基于能量耗散理论，结合蠕变速率特征，以耗散率值作为蠕变三阶段的控制阈值，对蠕变阶段临界分段，建立一个能够较为全面地反映岩石蠕变机制的非线性蠕变本构模型，并验证所建模型的合理性和适用性。

2 岩石三维蠕变模型

2.1 基础模型及黏塑形模型的选择

大量研究表明[1,4,7,10]，岩石蠕变全过程中包含明显的黏弹塑性行为。传统的元件模型如Maxwell模型、Kelvin模型、Bingham模型、Burgers模型等均由弹簧和牛顿体通过不同的串并联方式组合而成，无法描述岩石在加速蠕变阶段中的黏塑性变形。Cvisc模型(图1)在Burgers模型的基础上，串联了一个黏塑性体，其模型结构为H- N -H|N- N|S(H是弹性体，N是牛顿体，S是黏塑性体。“-”表示串联，“|”表示并联)。从模型组成上分析，包含弹性体、牛顿体和阻尼元件，具备一定的反映岩石弹性-黏性-黏弹性-黏塑性蠕变行为的能力。Cvisc模型的蠕变方程为

(1)

图1 Cvisc模型Fig.1 Cvisc model

式中：σ为蠕变应力；E1为瞬时弹性模量；E2为黏弹性模量；η1、η2和η3分别为模型第1、2、3部分的黏滞系数；ε为模型总应变；t为时间；σL为长期强度。

由于Cvisc模型中黏塑性阶段的应变率为常数，难以反映应力水平超过长期强度后岩石的加速蠕变，本文以Cvisc模型为基础模型进行改进。Perzyna黏弹塑性理论[14]被证明能较好地体现材料黏塑性行为的非线性特征，故本文用Perzyna黏弹塑模型[9]代替Cvisc模型第3部分中的黏塑性体，其黏塑性应变率方程为

(2)

式中：εvp为黏塑性应变，上标圆点为对时间t的导数；η为黏滞系数；{m}为黏塑性流动方向；“〈〉”是Macauley括号，括号内为分段函数；“{}”是表示向量的一种形式；φ(F)为屈服函数F的任意函数，可表示为[14-15]

(3)

式中：n为试验常数；σ0为初始屈服强度[7]。

2.2 基于能耗理论和蠕变速率的临界分段

岩石材料在外界应力场作用下发生变形的过程中，岩石不断从外界吸收能量，内部能量不断耗散。流变力学理论认为材料内部因素是决定蠕变力学行为的主要因素[1]，能量耗散理论[15]认为材料力学性能改变与内能变化有关，现目前大多数学者[1,4,10]仅通过研究岩石外部变量的变化规律来建立本构关系，虽然拟合效果较好但缺乏说服依据。实质上岩石蠕变变形过程中的能量耗散主要源于材料内部内能变化，耗散能不断地累积导致岩石损伤劣化效应增强，故而通过能耗理论结合蠕变速率，以“内外结合”的方式来反映岩石材料蠕变行为。为了进一步区分蠕变变形三阶段，采用耗散率值[16]来定义蠕变变形三阶段的控制阈值，示意图如图2所示。

图2 蠕变临界分段示意图Fig.2 Critical segmentsof creep

由图2可知，衰减、稳定、加速蠕变阶段耗散率控制阈值为γ1和γ2，其临界下限值分别为γ1*和γ2*。当耗散率超过γ1时，岩石从衰减蠕变阶段进入稳定蠕变阶段，当耗散率超过γ2时，岩石从稳定蠕变阶段进入加速蠕变阶段。

由能耗理论可知，耗散率是关于应力和不可逆应变率的函数，耗散率值γ表示为

(4)

假设在F的约束下有

(5)

式中λ为比例常数。

当γ≤γ1时，岩石处于衰减蠕变阶段，假设该阶段中蠕变变形包括可逆和不可逆的过程。当γ1<γ<γ2时，岩石处于稳定蠕变阶段，假设在该过程中(γ>γ1)蠕变变形不可逆，岩石从衰减蠕变阶段到稳定蠕变阶段在耗散率临界点经历初始屈服，结合式(3)、式(5)将φ(F)重新定义为

(6)

式中λ1为衰减和稳定蠕变阶段之间耗散率临界点的比例常数。

取衰减和稳定蠕变阶段之间耗散率临界点应变量值为εγ1，结合式(4)则有

(7)

由此可得临界下限值γ1*为

(8)

而{m}在一维应力状态下，其流动方向与塑形流动方向一致，故取{m}=1，于是将式(4)表示为

(9)

联立式(2)、式(6)—式(9)得到岩石在稳定蠕变阶段的黏塑性蠕变应变率模型，即

式中ηA为应变率模型中影响稳定蠕变阶段的黏滞系数。

ε(t)=a+blnt+ct+dtp。

(11)

式中a、b、c、d和p均为材料参数，在岩石未进入加速蠕变阶段时，p=0。

取p=0，对式(11)求导得到可用于稳定蠕变阶段的应变率，即

(12)

将式(12)代入式(10)可得

(13)

式(13)即为改进后的岩石稳定蠕变阶段黏塑性蠕变模型。

一般认为，岩石外界荷载超过长期强度时，岩石进入加速蠕变阶段[1,7,10]。在本文蠕变阶段划分中，当γ≥γ2时，岩石进入加速蠕变阶段，该阶段的耗散率确定方法与稳定蠕变阶段一致，参考式(8)可得稳定、加速蠕变两阶段之间的耗散率临界下限值为

(14)

对式(11)求导得到可用于加速蠕变阶段的应变率，即

(15)

参考式(10)的求解方法可得

(16)

式中：ηB为应变率模型中影响加速蠕变阶段的黏滞系数；λ2为稳定和加速蠕变阶段之间临界处的比例常数。

再将式(15)代入式(16)得

(17)

式(17)即为改进后的岩石加速蠕变阶段黏塑性蠕变应变率模型。

根据耗散能定义，耗散率值之间有如下关系，即

(18)

式(18)即为本文岩石蠕变三阶段临界分段的耗散率值关系。

2.3 基于能耗理论的岩石蠕变模型

在Cvisc模型的基础上，将改进后岩石稳定、加速蠕变阶段的黏塑性蠕变应变率模型应用到Cvisc模型中可得：

(19)

(20)

γ≥γ2*。

(21)

式(19)—式(21)即为基于能耗理论的岩石一维蠕变本构模型，通常取试验常数n=1[18]。当γ≤γ1时，该模型用于描述衰减蠕变阶段，所建蠕变模型退化为Cvisc模型中H-N-H|N结构；当γ1<γ<γ2时，该模型可用于描述稳定蠕变阶段，所建模型为该阶段黏塑性应变率模型与H-N-H|N结构的串联连接；当γ≥γ2时，该模型用可于描述加速蠕变阶段，所建模型为稳定、加速蠕变阶段黏塑性应变率模型与H-N-H|N结构的结合。

实际上在岩体工程中，岩石通常处于三向受力的状态，现室内压缩蠕变试验多为三轴试验，由此建立岩石三维应力状态下的蠕变模型十分必要。

假设岩石为各向同性体，根据广义Hooke定律，三维应力状态下的本构关系为

(22)

式中：σm为球应力张量；Sij为偏应力张量；K为体积模量；G为剪切模量;εm为体积应变张量;eij为偏应变张量。分解岩石内部的应力张量，于是有

σij=Sij+δijσm。

(23)

式中：σij为应力张量；δij为单位张量。

屈服函数可以取如下形式[19]，即

(24)

式中J2为应力偏量第二不变量。

在常规三轴压缩蠕变试验中，σ2=σ3，则有

(25)

一般认为，在三维应力下，岩石体积蠕变近乎为0，本文三维蠕变方程中不考虑体积蠕变，由此结合式(22)—式(23)和式(25)将式(19)—式(21)推广为三维应力状态，即：

(26)

(27)

γ≥γ2*。

(28)

将式(24)代入式(26)—式(28)即可得到本文所建基于能耗理论的岩石三维蠕变本构模型的轴向蠕变方程。

3 岩石蠕变试验

3.1 试验仪器及试样制备

蠕变试验采用RLW-2000型三轴流变试验系统(图3)。该系统由轴向加载系统、围压加载系统、孔压加载系统、伺服系统、控制系统、数据采集和自动绘图系统等部分组成。轴向加载系统和围压加载系统的控制部分采用全数字伺服控制器(EDC)，设备最大加载围压70 MPa，最大轴向荷载2 000 kN。岩样取自某水电工程高边坡滑带附近，制备成Φ50 mm×100 mm的圆柱样，砂岩岩样的典型破坏照片如图4所示。

图3 岩石三轴流变试验系统Fig.3 Triaxial rheological test system for rock

图4 典型破坏照片Fig.4 Typical damage photos

3.2 试验结果

蠕变试验采用逐级增量加载的方式，围压σ3设置为5 MPa，在三轴流变试验系统采集的数据中择取关键数据，通过玻尔兹曼[2,20]线性叠加处理得到分别加载蠕变曲线，将轴向应力荷载σ1标识于曲线上，如图5所示。

图5 分别加载蠕变曲线Fig.5 Creep curves under separate loading

由图5可看出，岩石在加载瞬间表现有瞬时弹性变形，接着进入衰减蠕变阶段，该阶段历时较短，蠕变速率逐渐降低，然后进入蠕变速率基本保持恒定的稳定蠕变阶段，当达到破坏应力水平时，岩石进入加速蠕变阶段，随即岩石屈服破坏。

引用张春阳等[21]的蠕变速率求解方法，其计算式为

(29)

式中：Δti为蠕变时间；n为蠕变试验数据个数。

通过式(29)的方法绘制蠕变速率曲线，如图6所示。

图6 蠕变速率曲线Fig.6 Curves of creep rate

由图6可知，曲线从左到右3个区域分别对应衰减、稳定和加速蠕变阶段，从曲线外侧到内侧，对应加载等级从第1级到第4级逐步提升，初始、稳态蠕变速率逐渐增大。

4 模型验证及参数求取

4.1 模型临界分段参数确定

通过最小二乘法用式(12)、式(15)对图6中蠕变速率点进行拟合，如图7所示。表1为拟合参数。

图7 蠕变速率拟合曲线Fig.7 Fitted curves of creep rate

表1 拟合参数Table 1 Fitting parameters

由图7和表1可知，式(12)和式(15)对蠕变速率点的拟合效果较好，平均R2达到0.967，通过数据拟合的方法确定了模型参数b、c、d和p。

表2 临界点参数Table 2 Parameters of thresholds

4.2 模型验证及参数求解

图8 等时应力-应变曲线Fig.8 Isochronousstress-strain curves

其余参数通过一般非线性最小二乘法拟合求解。为了对比预测效果，通过式(22)—式(25)的方法推广Cvisc模型到三维情形，根据室内砂岩蠕变试验数据验证本文所建模型和三维Cvisc模型，试验值与预测值对比曲线如图9所示，模型参数如表3所列。

图9 砂岩蠕变试验值与预测值对比曲线Fig.9 Comparison of sandstone creep betweentest and model prediction

表3 模型参数Table 3 Model parameters

由图9可看出，本文所建模型具有较强的预测能力，平均R2达到0.994 7，三维Cvisc模型平均R2为0.937 5。两模型在衰减、稳定蠕变阶段的辨识能力无显著差异，而在加速蠕变阶段中，本文所建模型具有明显的优越性。

为了验证本文基于能耗理论的岩石三维蠕变本构模型及临界分段方法描述岩石蠕变行为的适用性，引用相关文献中泥岩[20]和冻结软岩[22]的蠕变试验数据，通过本文模型和三维Cvisc模型进行辨识，预测结果如图10所示。模型临界分段参数、σL和σ0根据蠕变试验结果处理而得，限于篇幅，这里不再列出，其余模型参数值如表3所列。

图10 泥岩、冻结软岩蠕变试验值与预测值对比曲线Fig.10 Comparison of creep of mudstone and frozensoft rock between test and model prediction

结合图9、图10可看出，本文所建模型对砂岩、泥岩和冻结软岩蠕变数据拟合效果较好，平均R2分别达到0.987 9和0.990 2，能识别不同曲线形态的加速蠕变阶段，克服了传统元件模型难以反映加速蠕变阶段的困难，而传统三维Cvisc模型的平均R2分别仅有0.937 1和0.926 6。由此充分体现本文所建模型的优越性，证明本文基于能耗理论的岩石三维蠕变本构模型的合理性和适用性。