刘振等：循环加卸载下闪长玢岩蠕变特性及损伤本构模型 147 偏应力水平超过某一阈值时，试样产生的瞬时塑 介于固体和理想流体之间的元件.其模型的表达 性应变达到不容忽视的地步.为描述之一现象，将 式为： 一个弹性元件和一个塑性元件形成并联结构（如 d"s(t) r(0=刀drn 0≤n≤1 (9) 图4所示)来描述之一效果2四 式中：n为黏滞系数；n为Abel黏壶的分数阶导数 当=0时，该模型变成弹簧元件，当1时，该模 型又变成牛顿体.偏应力保持不变时，该元件三维 蠕变本构方程通过积分可得2： =-3. 0≤n≤1 (10) 3n I(n+1) 图4瞬时塑性模型 Fig.4 Transient plasticity model 将分数阶Abel黏壶与Kelvin体串联组成一个新 型黏弹性模型（如图5所示），其蠕变本构方程如下： 根据上文提到将静水压力造成的轴向应变分 (c1-c3r1-03 量视为不可恢复塑性应变的一部分，偏应力小于 3nceol(n+1)3Gce (11) 应力阈值S时岩石几乎不会产生瞬时塑性应变， 即塑性应变仅为静水压力引起的应变.偏应力大 于等于应力阈值S时，其三维本构方程如下： 5p-1m）-S0+o1=)-S0+ 9Kmp (6) 3Gmp 3Kmp 式中：S为应力阈值，即施加的应力超过该阈值时， 试样产生的塑性应变不容忽视；Kp为塑性体积模 图5新型黏弹性模型 量，可由式(3)确定；Gmp为塑性剪切模量，因为cmp Fig.5 New viscoelastic model 可由试验获得，所以剪切模量Gp可由下式确定： 式中：Ieeo为分数阶Abel黏壶黏滞系数；n介于 (σ1-03)-S0 0到1之间；r(n+1)为Gamma函数 Gmp= (7) (c1-03)+303-S0 3 Emp 为了验证该模型的优越性，利用1 stopt软件采 9Kme 用最小二乘法对其进行非线性拟合.表2给出了 2.3黏弹性模型 2模型的拟合参数，根据拟合结果可以看出改进黏 众所周知，在三轴压缩蠕变试验中，保持偏应 弹性体的拟合程度(R)均在0.995以上，拟合程度 力（σ一o?0)不变时，会产生随时间逐渐增大的蠕 更高.并且随着偏应力水平的增加Kelvin体的拟 变应变.而在卸载偏应力后，会产生黏弹性应变时 合程度越来越低，可见在高应力水平下Kelvin体 间恢复曲线.假设试样在黏弹性应变恢复期间，没 已经不能够承担描述试样黏弹性应变的任务，故 有损伤产生，那么在加载过程中黏弹性应变曲线 此本文选用新型黏弹性模型来描述黏弹性应变曲 可由黏弹性恢复曲线向上翻转得到3-2刈大多学 线.以偏应力90MPa为例，两模型拟合对比结果 者采用Kelvin模型来描述试样的黏弹性应变曲 如图6(a)所示. 线，Kelvin模型不仅能够很好地描述黏弹性应变曲 2.4黏塑性模型 线，而且还可以描述弹性后效.其三维本构方程如 蠕变应变为黏弹性应变和黏塑性应变之和， 式(8)所示： 所以用加载阶段的蠕变曲线减去黏弹性应变曲线 e= 1-3 (8) 即可得到黏塑性应变曲线，黏弹性应变曲线由上 3Gce 文提供的方法获得，又因为试验条件限制，本文试 式中：Gce表示Kelvin体中的剪切模量，ce表示Kelvin 验卸载阶段持续时间为24h,加载阶段持续时间 体中的黏滞系数. 为72h,故先通过式(11)与黏弹性应变曲线进行拟 但是，在高应力下岩石的黏弹性应变率呈现 合，求得黏弹性模型参数，然后再将所求参数带入 出明显的非线性现象，而Kelvin模型难以描述该 式(11)求出72h内的黏弹性应变，最后再根据上 现象.本文通过引进一个分数阶Abel黏壶来克服 述方法进行黏弹-塑性应变分离，分离后黏弹、塑 Kelvin模型的不足.众所周知分数阶Abel黏壶是 性应变曲线如图6(b)所示.偏应力水平超过某一阈值时，试样产生的瞬时塑 性应变达到不容忽视的地步. 为描述之一现象，将 一个弹性元件和一个塑性元件形成并联结构（如 图 4 所示）来描述之一效果[22] . Gmp S0 图 4    瞬时塑性模型 Fig.4    Transient plasticity model 根据上文提到将静水压力造成的轴向应变分 量视为不可恢复塑性应变的一部分，偏应力小于 应力阈值 S0 时岩石几乎不会产生瞬时塑性应变， 即塑性应变仅为静水压力引起的应变. 偏应力大 于等于应力阈值 S0 时，其三维本构方程如下： εmp = (σ1 −σ3)−S 0 3Gmp + (σ1 −σ3)−S 0 9Kmp + σ3 3Kmp （6） 式中：S0 为应力阈值，即施加的应力超过该阈值时， 试样产生的塑性应变不容忽视；Kmp 为塑性体积模 量，可由式（3）确定；Gmp 为塑性剪切模量，因为 εmp 可由试验获得，所以剪切模量 Gmp 可由下式确定： Gmp = (σ1 −σ3)−S 0 3 ( εmp − (σ1 −σ3)+3σ3 −S 0 9Kme ) （7） 2.3    黏弹性模型 众所周知，在三轴压缩蠕变试验中，保持偏应 力（σ1−σ3>0）不变时，会产生随时间逐渐增大的蠕 变应变. 而在卸载偏应力后，会产生黏弹性应变时 间恢复曲线. 假设试样在黏弹性应变恢复期间，没 有损伤产生，那么在加载过程中黏弹性应变曲线 可由黏弹性恢复曲线向上翻转得到[23−24] . 大多学 者采用 Kelvin 模型来描述试样的黏弹性应变曲 线，Kelvin 模型不仅能够很好地描述黏弹性应变曲 线，而且还可以描述弹性后效. 其三维本构方程如 式（8）所示： εce = σ1 −σ3 3Gce ( 1−exp( − Gce ηce t )) （8） 式中：Gce 表示Kelvin 体中的剪切模量，ηce 表示Kelvin 体中的黏滞系数. 但是，在高应力下岩石的黏弹性应变率呈现 出明显的非线性现象，而 Kelvin 模型难以描述该 现象. 本文通过引进一个分数阶 Abel 黏壶来克服 Kelvin 模型的不足. 众所周知分数阶 Abel 黏壶是 介于固体和理想流体之间的元件. 其模型的表达 式为： σ(t) = η d n ε(t) dt n 0 ⩽ n ⩽ 1 （9） 式中：η 为黏滞系数；n 为 Abel 黏壶的分数阶导数. 当 n=0 时，该模型变成弹簧元件，当 n=1 时，该模 型又变成牛顿体. 偏应力保持不变时，该元件三维 蠕变本构方程通过积分可得[25] ： ε = σ1 −σ3 3η · t n Γ(n+1) 0 ⩽ n ⩽ 1 （10） 将分数阶 Abel 黏壶与 Kelvin 体串联组成一个新 型黏弹性模型（如图 5 所示），其蠕变本构方程如下： εce = (σ1 −σ3)t n 3ηceoΓ(n+1) + σ1 −σ3 3Gce ( 1−exp( − Gce ηce t ))（11） Gce ηce ηceo n 图 5    新型黏弹性模型 Fig.5    New viscoelastic model 式中 ： ηceo 为分数 阶 Abel 黏壶黏滞系数 ； n 介 于 0 到 1 之间；Γ(n+1) 为 Gamma 函数. 为了验证该模型的优越性，利用 1stopt 软件采 用最小二乘法对其进行非线性拟合. 表 2 给出了 2 模型的拟合参数，根据拟合结果可以看出改进黏 弹性体的拟合程度（R 2 ）均在 0.995 以上，拟合程度 更高. 并且随着偏应力水平的增加 Kelvin 体的拟 合程度越来越低，可见在高应力水平下 Kelvin 体 已经不能够承担描述试样黏弹性应变的任务，故 此本文选用新型黏弹性模型来描述黏弹性应变曲 线. 以偏应力 90 MPa 为例，两模型拟合对比结果 如图 6（a）所示. 2.4    黏塑性模型 蠕变应变为黏弹性应变和黏塑性应变之和， 所以用加载阶段的蠕变曲线减去黏弹性应变曲线 即可得到黏塑性应变曲线，黏弹性应变曲线由上 文提供的方法获得. 又因为试验条件限制，本文试 验卸载阶段持续时间为 24 h，加载阶段持续时间 为 72 h，故先通过式（11）与黏弹性应变曲线进行拟 合，求得黏弹性模型参数，然后再将所求参数带入 式（11）求出 72 h 内的黏弹性应变，最后再根据上 述方法进行黏弹‒塑性应变分离，分离后黏弹、塑 性应变曲线如图 6（b）所示. 刘    振等： 循环加卸载下闪长玢岩蠕变特性及损伤本构模型 · 147 ·