水利水电、隧道及采矿工程建设中，受工程扰动等因素的影响，实际工程岩体内产生大量非连续的微观或宏观裂隙。在自然界或实际工程中，如库岸边坡、地下水库坝体和防水煤柱等，具有初始损伤的煤（岩）体与水直接充分接触，承受荷载和静水压力的耦合作用。随着压力水–岩耦合作用时间增加，岩石强度力学特性损伤劣化现象更加明显，会导致出现与水相关的工程灾害，如Malpasset拱坝溃决^[1]、三峡库区蓄水后的库岸滑坡^[2]，以及高水压下隧洞围岩蠕变对锦屏二级水电站长期稳定性的影响^[3]等。因此，研究压力水下裂隙岩石的蠕变特性，对评价处于荷载与压力水耦合作用下岩土工程的长期稳定性具有重要的科学意义。

目前，国内外学者对裂隙岩石蠕变特性开展了大量研究工作。杨红伟等^[4]采用煤岩流变仪，进行细砂岩孔隙水压分级加载蠕变试验，发现细砂岩孔隙体积的演化规律与微孔洞的损伤破坏特征相吻合。Zhao等^[5]在石灰岩上预制裂隙，开展了三轴加卸载试验，探讨加卸载下裂纹岩石的蠕变特性变化规律。牛双建等^[6]采用岩石三轴流变仪，对损伤破裂砂岩进行单轴蠕变试验，分析了相应的蠕变特性。陈芳^[7]对裂隙辉绿岩进行了三轴蠕变试验，分别从轴向、侧向及体积变形等方面分析了蠕变变化规律，并探讨了岩石蠕变损伤阈值。Shao^[8]根据微观结构建立了流变模型，基于脆性岩石损伤特点，引入2阶损伤张量，构建了含损伤变量的蠕变模型。王其虎等^[9]基于应力与法向应变关系，提出裂隙岩石塑性蠕变元件，描述岩石的瞬时塑性变形。陆银龙等^[10]采用MATLAB程序与Comsol模拟软件相结合，模拟了含微裂纹岩石的蠕变损伤破坏全过程。

上述研究成果对认识裂隙岩石的蠕变力学特性具有重要意义。但现阶段水压影响下煤（岩）体的蠕变特性研究，主要以预制裂纹岩样为对象，采用试样端部施加水压的三轴试验或数值模拟的方法，研究方法和设备的模拟工况与实际工况存在一定偏差。基于此，人工制备含有初始损伤的非连续裂隙砂岩岩样，借助考虑荷载与静水压力耦合作用的多通道流固耦合岩石流变试验系统，进行轴压水压耦合作用下三轴蠕变试验，揭示轴压水压耦合作用下裂隙岩石的蠕变损伤劣化演变规律，构建考虑时效损伤的3维非线性蠕变本构模型，以期为评价处于轴向荷载与压力水耦合作用的岩土工程的长期稳定性提供理论指导。

1 裂隙砂岩制备 1.1 岩样成分

砂岩岩样取自川东石炭二叠纪的隆昌青砂岩，利用四川大学扫描电镜能谱仪对砂岩的主要成分进行分析，如图1所示。

由图1可知，砂岩试样中主要含有O、Si、Na、Mg和K等元素，其中，O元素的含量最高，其他依次为Si、Al、Fe等元素。颗粒成分主要为石英、方解石或长石等较坚硬的造岩矿物。

1.2 裂隙试样制备步骤

借助四川大学MTS试验机制备裂隙砂岩岩样。裂隙砂岩岩样在制备过程中需严格按照国际岩石力学学会室内试验岩石力学试验标准；加工时，尤其注意岩石的各向异性。制备步骤为：首先，采用轴力控制方式，以60 kN/min的加载速率加载；当轴向荷载接近岩样屈服强度时，转换为位移控制方式，以0.15 mm/min的加载速率加载；当应力–应变曲线达到峰值且具有下降趋势时，终止试验，进行卸载，完成裂隙砂岩岩样制备工作。制备好的部分裂隙砂岩试样如图2所示。其中，裂隙试样在单轴压缩条件下两端会局部产生端部效应，但裂隙岩样整体保持着较好的完整性，无宏观裂隙出现。

为控制样品的离散性，消除因岩样性质差别导致的结果误差，首先，对试样进行超声波测试，剔除纵波波速异常的岩样；其次，测试样的密度，根据密度进行试样分组，保证同一组试样的密度差别不大。

2 轴压水压耦合作用下裂隙砂岩蠕变特性试验 2.1 试验设备

岩石流变试验系统（YSL–200）由静油压加载系统、静水压加载系统和位移变形测量系统组成（图3）。静油压加载系统通过压缩油缸对试样实施轴向加载；静水压加载系统通过加压承压筒内的水对试样进行围压加载；位移变形测量系统采用高精度光栅尺测量试样的轴向位移量。其中：轴压力范围为10～200 kN；水压范围为0～10 MPa；光栅尺量程为0～25 mm，测量精度为±0.2 μm。

该岩石流变系统对岩样的轴压水压耦合作用原理与MTS中的轴压孔隙水压联合作用原理不同，孔隙水压仅仅是在两端施加水压力，水从一端的孔中流入，从另一端的孔中流出，从而实现水–岩耦合作用，是一个动水压过程；岩样围压是通过控制油压完成。本文借助的岩石流变系统，真正实现了水与整块岩样充分接触，通过静水压加载系统施加水压力，实现荷载与压力水耦合作用，是一个静水压过程，相应的岩样受力示意图如图4所示。

2.2 试验方案及步骤 2.2.1 试验方案

从两个方面对轴压水压耦合作用下裂隙砂岩蠕变特性进行试验研究，即：1）不同初始轴荷载下单级加载的三轴蠕变特性分析（水压不变、轴压变）；2）不同水压力下单级加载的三轴蠕变特性分析（轴压不变、水压变）。根据研究内容，共设计了5组三轴蠕变试验方案，如表1所示。

2.2.2 试验步骤

试验中所用的试验岩样均符合ISRM试验标准。共选取15块制备好的裂隙砂岩岩样，分成5组，每组3个，按照试验方案分别开展相应的三轴蠕变试验。以试验1为例描述试验步骤。

首先，按照饱水试验规程对3块岩石进行浸泡48 h；其次，将浸泡好的岩样分别放到固定岩样装置内，封桶和注水饱和72 h（此时认为试样已达到完全饱和，水对岩样基质的软化作用默认完成）；然后，设置三轴蠕变试验参数，如岩样直径及高度、目标水压（3 MPa）等；最后，变形清零，运行试验，保存试验数据。试验结束后，取出岩样，整理蠕变试验数据，进行下一组试验。

2.3 试验结果及分析 2.3.1 蠕变试验结果

采用系统自带的TestPilot软件采集试验数据，该软件能够准确记录试验全过程的荷载、位移与时间等数据，并根据设置的试验参数，自动计算相应的应力、应变值。5组不同试验条件下的三轴蠕变试验结果见图5～8。

由图5～8可知，裂隙岩样蠕变曲线先后经历了减速蠕变、等速蠕变和加速蠕变阶段。对于蠕变试验1～5，裂隙砂岩岩样分别在经历了310、210、15、300、350 h后发生蠕变破坏现象，发生蠕变破坏时的最大轴向蠕变分别为410×10^–6、760×10^–6、1 000×10^–6、390×10^–6与280×10^–6。

2.3.2 水压对裂隙岩样蠕变特性的影响

为分析在相同初始荷载作用下，水压对裂隙砂岩轴向蠕变变化规律的影响，选取表1中的试验1、4、5，分别对这3组试验的最大轴向蠕变量与蠕变速率进行对比分析。

1）水压与最大轴向蠕变的关系

水压与裂隙砂岩岩样的最大轴向蠕变关系如图9所示。

由图9可知，在相同的初始荷载作用下，裂隙砂岩轴向蠕变应变随水压增大而减小，这是由于裂隙砂岩岩样的抗压强度随着水压升高而增大。随着水压的增大，裂隙岩样的侧向力会随之增大，进而增大裂隙岩样裂纹的横向扩展，束缚裂隙岩样裂纹的纵向扩展，从而减小裂隙岩样在轴向方向的蠕变变形量。

2）水压与轴向蠕变速率的关系

为分析水压对裂隙砂岩蠕变速率变化规律的影响，选取表1中的试验1、4、5的数据，计算得到3组裂隙砂岩蠕变试验的蠕变速率，如图10所示。

由图10可知：在相同初始荷载、不同水压下，裂隙岩样蠕变过程先后经历了减速蠕变、等速蠕变和加速蠕变阶段。当裂隙砂岩刚进入蠕变阶段时，岩石蠕变速率先达到最高点，但持续时间较短；随着蠕变时间增加，蠕变速率不断减小，此阶段为减速蠕变阶段；当蠕变速率趋于稳定（图10（c）中的曲线弯折是由于试样中局部裂隙扩展后又被压密所致），并在较长的时间内保持不变，说明裂隙砂岩在这段时间内进入了稳定蠕变阶段；最后，由于微裂隙的贯通与扩展，产生了主裂纹，蠕变速率开始从稳定值逐渐增大，并在裂隙岩样破坏之前达到最大，说明在这段蠕变变化周期内，裂隙岩样蠕变已进入加速蠕变阶段。

上述分析可知，裂隙岩样的最大轴向蠕变量随水压增大而减小。基于以上3组裂隙砂岩岩样蠕变速率的变化规律，定量分析水压对蠕变加速阶段蠕变速率的影响，如图11所示。

由图11可知：裂隙砂岩的蠕变加速阶段的最大蠕变速率随着水压力的升高而增大。这是由于水压束缚了裂隙岩样内部裂纹的横向贯通与扩展，间接促进裂纹在纵向与横向之间扩展，诱发较多的剪切微裂隙；待裂隙岩样达到加速蠕变阶段时，会形成主剪切裂纹；当张拉主裂纹与剪切主裂纹贯通后，岩样将快速发生蠕变破坏现象。

2.3.3 初始轴应力对裂隙岩样蠕变速率的影响

为分析在相同水压力下，初始轴应力对裂隙砂岩蠕变速率的影响变化规律，选取表1中的试验1、2、3数据进行研究，结果如图12所示。

大量试验研究表明，增大初始荷载会提高岩样蠕变量，这与现实相符^[11]，不再赘述相关研究内容，作者主要分析初始荷载对减速蠕变阶段（初始阶段）和加速蠕变阶段（破坏阶段）的裂隙岩样最大轴向蠕变速率的影响变化规律。

由图12可知，初始阶段和破坏阶段的最大轴向蠕变速率随初始荷载增大而增大，这和初始荷载与蠕变应变量的变化关系相吻合。这也说明在相同水压环境下，增大初始荷载的力学作用效应为先提高岩样的蠕变速率，最终由蠕变应变量呈现初始荷载的影响。

3 轴压水压耦合作用下裂隙砂岩蠕变模型

裂隙砂岩岩样中的裂隙呈现非连续性、非均一性，由于大量微裂纹的存在，不能采用常规的蠕变模型分析和研究轴压水压耦合作用下裂隙砂岩的蠕变特性，需要构建新的裂隙砂岩蠕变模型。新构建的裂隙砂岩蠕变模型不仅需要考虑裂隙岩样自身的损伤部分，也要考虑压力水环境对裂隙岩样的弱化作用，更重要的是要考虑裂隙岩样的时效损伤。因此，基于裂隙砂岩的蠕变特性，借助损伤力学理论，引进损伤因子，建立考虑时效损伤的非线性蠕变本构模型。

3.1 裂隙岩样损伤因子的确定

断续裂隙岩样是完整岩样受高强度压缩所致，因此，裂隙岩样在进行蠕变试验前已有损伤。从损伤力学的角度，引入损伤因子D描述裂隙砂岩刚性流变参数随时间的变化情况，构建新的蠕变模型元件及含损伤的1维非线性蠕变模型。在损伤力学中，常通过弹性模量变化比例定义损伤因子^[12]，即：

$D(t) = \frac{{{E_0} - {E_t}}}{{{E_0}}}$

(1)

式中：D(t)为岩石试样在t时刻的损伤因子；E₀为初始时刻岩样的弹性模量，GPa；E_t为t时刻岩样的弹性模量，GPa。

张强勇等^[13]对在蠕变压缩试验中获得的岩石弹性模量进行拟合分析，发现岩石弹性模量在进行蠕变变形过程中具有时效弱化规律，与时间呈负指数关系，表达式为：

$D'(t) = \frac{{{E_0} - {E_\infty }}}{{{E_0}}}[1 - \exp( - \gamma t)]$

(2)

式中：γ为与损伤变量有关的系数；E_∞为岩样的长期弹性模量，GPa。

以式（2）为基础，在经典蠕变模型中引入损伤因子，建立非线性的裂隙蠕变模型。由于裂隙砂岩岩样存在初始损伤，需对式（2）进行完善，修正后的表达式为：

${\;\;\;\;\;\;\;D''(t)} = (1 - {D_0}) \cdot \frac{{{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)}}{{{E_0}}}$

(3)

式中，D₀为裂隙岩样在蠕变开始时的损伤因子。

对于裂隙砂岩非线性蠕变模型中初始损伤因子参数的辨识，需根据裂隙砂岩岩样进行蠕变试验前的应力应变试验数据确定。则求解损伤因子的公式为^[13]：

$D = \frac{{E{\varepsilon _1} - {\sigma _1} + \mu \left( {{\sigma _2} + {\sigma _3}} \right)}}{{E{\varepsilon _1}}}$

(4)

由于裂隙砂岩岩样在进入蠕变试验之前，先进行的是轴压水压耦合作用下的三轴压缩试验，所以式（4）需改写为：

${D_0} = \frac{{E{\varepsilon _1} - {\sigma _1} + 2\mu \cdot p}}{{E{\varepsilon _1}}}$

(5)

式中： ${\varepsilon _1}$ 为轴应变； $\;\mu $ 为泊松比； ${\sigma _1}$ 为轴应力，MPa； $p$ 为水压，MPa； $E$ 为常规压缩试验的弹性模量，GPa。

3.2 裂隙蠕变模型基本元件的参数非线性研究

裂隙蠕变模型基本元件中，同样包括弹性元件、黏性元件和塑性元件，从岩样自身角度出发，需要对基本元件中的线性元件进行非线性化处理，在基本元件中能够进行非线性处理的参数为黏滞系数和弹性模量。根据式（5），分别对黏滞系数和弹性模量进行非线性处理，处理后的表达式为：

$\left\{ \begin{array}{l} \eta (t) = {\eta _0}(1 - D(t)) = \dfrac{{{\eta _0}}}{{{E_0}}} \cdot (1 - {D_0}) \cdot \dfrac{{ {{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)} }}{{\exp({\rm{ - }}\gamma t)}}, \\ E(t) = {\rm{E}}(1 - D(t)) = \dfrac{E}{{{E_0}}} \cdot (1 - {D_0}) \cdot \dfrac{{ {{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)} }}{{\exp({\rm{ - }}\gamma t)}} \end{array} \right.$

(6)

由式（6）可知，裂隙模型基本元件的非线性变形特征是通过引入损伤因子对黏滞系数和弹性模量进行弱化处理来体现的。根据处理后的非线性参数表达式，构建基本元件的非线性本构方程，并描述其变形特征。

1）弹性损伤元件

由文献[14]可知弹性元件的本构方程为：

$\varepsilon (t) = \frac{\sigma }{{E(t)}}$

(7)

结合弹性模量的非线性方程式（6），整理可得考虑时效损伤的弹性模量的非线性表达式为：

${\;\;\;\;\;\;\;\;\varepsilon (t)} = \frac{{{E_0}}}{{E \cdot (1 - {D_0})}} \cdot \frac{{{\sigma _0}\exp({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)}}$

(8)

由式（8）可知，弹性元件的应变量随时间的增加而增大。这是由于弹性模量是关于时间的函数，随时间的增加而减小，这说明弹性损伤元件能够呈现岩样的损伤时效劣化现象。

2）黏性损伤元件

由文献[14]可知，黏性元件的本构方程为：

$\sigma = \eta \frac{{{\rm{d}}\varepsilon }}{{{\rm{d}}t}}$

(9)

结合黏滞系数的非线性关系表达式（6），整理得到考虑损伤变量的黏滞系数的非线性本构方程为：

$ {\;\;\;\;\;\;\;\;\frac{{{\rm{d}}\varepsilon (t)}}{{{\rm{d}}t}}} = \frac{{{E_0}}}{{{\eta _0} \cdot (1 - {D_0})}} \cdot \frac{{{\sigma _0}\exp({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)}} $

(10)

经积分得：

$\varepsilon (t) = \frac{{{\sigma _0} \cdot {E_0}}}{{\gamma {\eta _0} \cdot (1 - {D_0}) \cdot {E_\infty }}}\ln \frac{{{E_0} - {E_\infty } + {E_\infty }\exp({\rm{ - }}\gamma t)}}{{{E_0}}}$

(11)

当式（11）中的γ=0时，黏滞系数非线性表达式演变为线性本构方程（式（9））。一方面，说明构建的非线性黏性元件本构方程是正确的；另一方面，说明了参数γ能够描述岩样的黏性损伤程度。

3.3 基本元件组合的非线性蠕变模型研究

基于提出的非线性参数表达式，在构建裂隙蠕变模型过程中，优先考虑时间效应，在构建蠕变本构方程时，须引入非线性参数表达式，建立相应的蠕变状态方程。

3.3.1 Kelvin损伤模型

Kelvin（K）力学蠕变模型由弹性元件与黏性元件并联组成，K模型的本构方程^[14]为：

$\sigma = E(t) \cdot \varepsilon + \eta (t) \cdot \dot \varepsilon $

(12)

将式（6）中的参数E(t)和η(t)表达式代入式（12）中，其中初始条件为：

$\varepsilon (t) = \left\{ \begin{array}{l} 0,\;t = 0; \\ \dfrac{{{\sigma _0}}}{{(1 - {D_0}){E_{\rm{K}}}}} \cdot \dfrac{{{E_0}}}{{{E_\infty }}},\;t \to \infty {\text{。}}\end{array} \right.$

整理得K损伤模型的蠕变方程为：

$\varepsilon (t) = \frac{{{\sigma _0}}}{{( (1 - {D_0}){E_{\rm{K}}}}}\left(\frac{{{E_0}\exp ({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp ({\rm{ - }}\gamma t)}} - \exp \left( - \frac{{{E_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t\right)\right)$

(13)

当系数γ=0时，式（13）演变为经典的K模型，这证明了改进的K损伤蠕变方程的准确性，也说明了与经典K模型的相容性，但由于损伤变量的存在，又有别于经典K蠕变模型。

3.3.2 非线性黏塑性模型

研究发现，岩石发生蠕变破坏的主要原因是岩石的黏滞系数急剧下降^[15-16]。因此，针对裂隙砂岩在轴压水压耦合作用的特殊试验环境，需构建新的非线性黏塑性模型。

在典型经验方程法中，幂函数方程能够较好地描述非线性蠕变问题^[17]，一般的幂函数蠕变方程为：

$\dot \varepsilon = R{\sigma ^a}{t^b}$

(14)

式中：t为蠕变时间，h；R、a和b为拟合蠕变参数。由试验结果分析可知，蠕变加速阶段黏滞系数η主要受初始荷载、水压力p和时间t影响。因此，可参照幂函数蠕变方程建立一个关于p、t和σ的η函数，该函数可表示为：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\eta (\sigma ,p,t)} = \frac{{\sigma \cdot {\eta _0}}}{{\left( {1 + n} \right) \cdot {{\left( {\sigma {{ + p}}} \right)}^m}{t^n}}}$

(15)

式中：η₀为初始黏滞系数，GPa·h；p为水压力，MPa；t为蠕变时间，h；m和n为拟合蠕变参数。

线性黏弹体本构方程^[14]为：

$ \dot \varepsilon = \frac{\sigma }{\eta } $

(16)

则非线性黏弹体的本构方程可写为：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\dot \varepsilon} = \frac{\sigma }{{\eta (\sigma ,p,t)}} = \frac{{{{\left( {\sigma {{ + p}}} \right)}^m}{t^n}}}{{{\eta _0}}}$

(17)

假设非线性黏塑性模型（MNY模型）由摩擦片和非线性黏壶并联而成。将式（15）积分，整理得到描述裂隙砂岩的非线性黏塑性蠕变方程为：

$ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} {\varepsilon} = 0,\;\sigma < {\sigma _{\rm{f}}};\\ \varepsilon = \dfrac{{{{({\sigma _0} - {\sigma _{\rm{f}}}{\rm{ + }}p)}^m}{t^n}}}{{{\eta _{{\rm{MNY}}}}}},\;\sigma \ge {\sigma _{\rm{f}}} \end{array} \right.} $

(18)

为了分析所得各参数对新非线性黏性元件蠕变力学特性的影响，对式（18）中的参数m、n逐一讨论，分析各参数对非线性黏性元件蠕变特性的影响变化规律，并对所提出的非线性黏性元件的适应性进行描述。

当分析参数m对蠕变特性的影响时，设定初始黏滞系数η₀=500 GPa·h，轴向应力σ₀=25 MPa，水围压p=2 MPa，参数n=1.5。当参数m的取值分别为–0.5、0.3、0.7、1.0、1.5时，加速阶段的蠕变随时间的变化规律如图13所示。

由图13可知：参数m不断增大时，模型的蠕变量及蠕变速率随之增大。当蠕变时间较短时，不同参数m下的蠕变量相差不大，且增幅较小；随着蠕变时间的增加，同一参数m下的蠕变速率呈递增趋势。

分析参数n对蠕变的影响时，假设黏滞系数η₀=500 GPa·h，轴向应力σ₀=25 MPa，水围压p=2 MPa，参数m=1。当参数n的取值分别为0.3、0.6、1.0、1.5时，加速阶段的蠕变随时间的变化规律见图14。

由图14可知：随着参数n取值不断增大，模型蠕变速率随之增大，物理材料在水压作用下，发生蠕变破坏的时间也越短，尤其当参数n取值大于1之后。当参数n小于1时，非线性模型元件可描述过渡蠕变；当参数n等于1时，非线性模型元件可描述等速蠕变；当参数n大于1时，非线性模型元件可描述加速蠕变。

3.4 轴压水压下考虑时效损伤的非线性裂隙砂岩蠕变模型 3.4.1 考虑时效损伤的1维非线性蠕变组合模型

基于裂隙砂岩在轴压水压耦合作用下的时效损伤蠕变变形特征，引入弹性元件（H模型）描述裂隙砂岩的瞬时变形特征；引入考虑时效损伤的非线性K模型描述裂隙砂岩的减速蠕变阶段和稳定蠕变阶段；引入MNY模型描述裂隙砂岩的加速蠕变阶段，如图15所示。构建考虑时效损伤的1维蠕变方程的具体过程如下。

根据蠕变模型的串联组合特征，可获得轴压水压下考虑时效损伤的非线性蠕变状态方程为：

1）当 $\sigma < {\sigma _{\rm{f}}}$ 时，考虑时效损伤的1维非线性蠕变模型由H模型与K模型组成，该蠕变模型能够描述轴压水压下裂隙砂岩的瞬时弹性变形和减速蠕变阶段。因为H模型描述的岩石弹性变形是瞬间完成的，不需要考虑时效损伤，这里只需要对K模型考虑时效损伤，因此，相应的蠕变方程可表示为：

$\varepsilon (t) = {\varepsilon _{\rm{H}}} + {\varepsilon _{\rm{K}}}$

(19)

式中： ${\varepsilon _{\rm{H}}}$ 为H模型在恒定荷载作用下发生的瞬时弹性应变， ${\varepsilon _{\rm{H}}} = \dfrac{{{\sigma _0}}}{{{{{E}}_{\rm{H}}}}}$ ； ${\varepsilon _{\rm{K}}}$ 为K模型在考虑时效损伤下的非线性蠕变， ${\varepsilon _{\rm{K}}} = \dfrac{{{\sigma _0}}}{{(1 - {D_0}){{{E}}_{\rm{K}}}}}\left(\dfrac{{{{{E}}_0}\exp ({\rm{ - }}\gamma t)}}{{{{{E}}_0} - {{{E}}_\infty } + {{{E}}_\infty }\exp ({\rm{ - }}\gamma t)}} -\right. \left.\exp \left( - \dfrac{{{{{E}}_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t\right)\right)$ 。

2）当 $\sigma \ge {\sigma _{\rm{f}}}$ 时，考虑时效损伤的1维非线性蠕变模型由H模型、K模型与MNY模型串联组成，能够描述轴压水压下裂隙砂岩蠕变曲线的瞬时变形、减速蠕变阶段和加速蠕变阶段，相应的非线性蠕变方程为：

$\varepsilon (t) = {\varepsilon _{\rm{H}}} + {\varepsilon _{\rm{K}}} + {\varepsilon _{{\rm{MNY}}}}$

(20)

式中， ${\varepsilon _{{\rm{MNY}}}}$ 为MNY模型在考虑时效损伤下的非线性蠕变， ${\varepsilon _{{\rm{MNY}}}} = \dfrac{{{{({\sigma _0} - {\sigma _{\rm{f}}}{\rm{ + }}p)}^{{m}}} \cdot {t^n}}}{{{\eta _{{\rm{MNY}}}}}}$ 。

3.4.2 考虑时效损伤的3维非线性蠕变组合模型

采用应力、应变张量的方法将考虑时效损伤的1维非线性蠕变模型转换为3维非线性蠕变模型。在三轴应力的作用下，相应的轴向蠕变关系表达式如下：

1）3维偏应变表达式

根据偏应变转换原理^[18-19]，将新构建的1维轴向蠕变方程转换为相应的3维偏应变方程。整理得到蠕变模型的3维偏应变表达式为：

当 $ {\sigma _0} < {\sigma _{\rm{f}}}$ 时，

$\begin{aligned}[b]{e_{ij}} =& \frac{{S_{ij}'}}{{2{G_{\rm{H}}}}} + \frac{{S_{ij}'}}{{2(1 - {D_0}) \cdot {G_{\rm{K}}}}}\cdot\\ &\left[ {\frac{{{E_0}\exp ({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp ({\rm{ - }}\gamma t)}} - \exp \left( { - \frac{{{G_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t} \right)} \right]\end{aligned}$

(21)

当 ${\sigma _0} \ge {\sigma _{\rm{f}}}$ 时，

$\begin{aligned}[b]{e_{ij}} =& \frac{{S_{ij}'}}{{2{G_{\rm{H}}}}} + \frac{{S_{ij}'}}{{2(1 - {D_0}) \cdot {G_{\rm{K}}}}}\cdot\\ &\left[ {\frac{{{E_0}\exp ({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp ({\rm{ - }}\gamma t)}} - \exp \left( { - \frac{{{G_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t} \right)} \right] +\frac{{{{({S_{ij}'} - {\sigma _{\rm{f}}}{\rm{ + }}p)}^m}{t^n}}}{{2{\eta _{{\rm{MNY}}}}}}\end{aligned}$

(22)

2）轴向蠕变方程

式（21）～（22）联合球应变张量 ${\delta _{ij}}{\varepsilon _{\rm{m}}}$ ，整理得到蠕变模型的3维轴向蠕变方程为：

当 ${\sigma _0} < {\sigma _{\rm{f}}}$ 时，

$\begin{aligned}[b]\varepsilon (t) =& \frac{{{\sigma _0} + 2p}}{{9K}} + \frac{{{\sigma _0} - p}}{{3{G_{\rm{H}}}}} + \frac{{{\sigma _0} - p}}{{3(1 - {D_0}) \cdot {G_{\rm{K}}}}}\cdot\\ &\left[ {\frac{{{E_0}\exp ({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp ({\rm{ - }}\gamma t)}} - \exp \left( { - \frac{{{G_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t} \right)} \right]\end{aligned}$

(23)

当 $ {\sigma _0} \ge {\sigma _{\rm{f}}}$ 时，

$\begin{aligned}[b]\varepsilon (t) = &\frac{{{\sigma _0} + 2p}}{{9K}} + \frac{{{\sigma _0} - p}}{{3{G_{\rm{H}}}}} + \frac{{{\sigma _0} - p}}{{3(1 - {D_0}) \cdot {G_{\rm{K}}}}}\cdot\\ &\left[ {\frac{{{E_0}\exp ({\rm{ - }}\gamma t)}}{{{E_0} - {E_\infty } + {E_\infty }\exp ({\rm{ - }}\gamma t)}} - \exp \left( { - \frac{{{G_{\rm{K}}}}}{{{\eta _{\rm{K}}}}}t} \right)} \right] + \frac{{{{({\sigma _0} - {\sigma _{\rm{f}}})}^m}{t^n}}}{{3{\eta _{{\rm{MNY}}}}}}\end{aligned}$

(24)

4 讨　论

为验证所提出的考虑时效损伤的非线性蠕变本构模型的准确性与合理性，以试验1、2的裂隙砂岩单级加载蠕变试验曲线为例，采用MATLAB软件中的牛顿插值迭代法，对提出的非线性蠕变模型参数进行辨识；然后，绘制相应的理论蠕变曲线，并与蠕变试验结果和经典西原蠕变模型预测的理论结果对比分析结果，如图16、17所示，相应的非线性蠕变模型辨识参数如表2所示。

由图16、17可知：单级加载下轴压水压耦合作用的裂隙砂岩蠕变试验结果与提出的蠕变模型预测的理论曲线基本相吻合；经典西原蠕变模型仅能很好地预测瞬时和稳定蠕变阶段，不能描述加速蠕变阶段。一方面说明MATLAB软件对非线性蠕变模型参数辨识比较可靠，另一方面说明针对裂隙砂岩蠕变特性提出的考虑时效损伤的非线性蠕变模型比较准确与合理，可为轴压水压耦合作用下裂隙岩石的长期蠕变特性研究提供理论参考。

5 结　论

1）借助MTS试验机，通过预加载的方式，制备了含有随机裂纹分布的裂隙砂岩岩样。采用特制的流固耦合岩石流变试验系统，开展了单级加载蠕变试验，真正实现了压力水–岩耦合作用。

2）裂隙砂岩岩样在轴压水压耦合作用下产生了蠕变破坏现象，先后经历了蠕变初始阶段、蠕变稳定阶段及蠕变加速阶段。轴向蠕变量随水压增大而减小；蠕变加速阶段的最大蠕变速率随水压升高而增大；初始阶段和破坏阶段的最大轴向蠕变速率随初始荷载增大而增大。

3）基于损伤理论，建立了考虑初始损伤的非线性基本元件；根据幂函数方程，提出本文的非线性黏塑性模型，用于描述蠕变加速破坏阶段；基于增量理论，构建了轴压水压下考虑时效损伤的3维非线性裂隙砂岩蠕变模型，能够较好地描述轴压水压耦合作用下裂隙砂岩蠕变特征，为评价处于压力水–岩耦合作用下岩土工程的长期稳定性提供了理论依据，有助于完善传统蠕变本构模型理论。