水岩耦合作用下碳质泥岩蠕变特性及本构模型

Study on Creep Behavior and Constitutive Model of Carbonaceous Mudstone Under Hydro-mechanical Coupling

摘要

为揭示复杂水岩耦合作用下碳质泥岩蠕变行为的时效机理，针对高渗透压与偏应力共同作用的复杂工况，开展系统的蠕变试验研究，旨在突破传统模型难以同时描述水岩耦合效应与非线性蠕变全过程的难点。通过开展碳质泥岩在不同渗透压力与偏应力共同作用下的蠕变试验，系统研究了其蠕变变形特性、应变时程规律和破坏机制，并构建了水岩耦合非线性黏弹塑性蠕变模型。结果表明：1）渗透压显著促进了碳质泥岩的蠕变变形，同时降低了其长期强度，致使其在较低应力水平下即发生蠕变破坏；2）渗透压对破坏形态具有重要影响，特别是在低围压条件下，会削弱围压对内部裂隙的径向约束作用，导致宏观破坏形态呈现随机性特征，出现复合型破坏模式；3）在渗透压作用下，碳质泥岩蠕变破坏的细观断口形态发生明显变化，拉裂破坏表现为鳞片状脆性断口、台阶花样和根状花样，剪裂破坏则呈现韧窝、穿晶裂纹和滑移分离等特征，细观与宏观破坏特征具有良好的对应关系；4）所建立的水岩耦合非线性黏弹塑性蠕变模型能准确表征不同围压和渗透压条件下的蠕变全过程特征。本文模型以围压和渗透压为变量，参数确定方法简便可靠，基于室内试验数据获得的计算结果与实测值吻合良好。

Abstract

Objective Carbonaceous mudstone exhibits low mechanical strength, water-induced disintegration, and reduced structural compactness. It is widely distributed in the mountainous strata of Southwestern China, where tunneling activities frequently traverse formations dominated by carbonaceous mudstone. The surrounding rocks are subjected to prolonged coupled effects of in-situ stress and groundwater seepage pressure, which complicates the creep behavior of carbonaceous mudstone. Constitutive models are essential for characterizing creep-related mechanical properties and deformation mechanisms. However, current models derived from elastoplastic theory inadequately address hydro-mechanical coupling effects during creep. Therefore, establishing a coupled hydro-mechanical creep constitutive model is imperative for accurate deformation prediction in geotechnical engineering. Methods This study conducted laboratory creep tests under hydro-mechanical coupling conditions on carbonaceous mudstone samples collected from secondary lining fracture zones in the Yanglin Tunnel. The experiments revealed the evolution of deformation characteristics and macro-micro fracture mechanisms throughout the creep process. A three-stage nonlinear viscoelastic-plastic creep model incorporating hydro-mec-hanical coupling was developed based on rheological and elastoplastic theories, and parameter identification methods were established. The theoretical curves showed strong agreement with the experimental data, accurately captured the complete creep behavior of carbonaceous mudstone, and demonstrated the model's validity for engineering applications. Results and Discussions 1) Creep curves exhibited stepwise progression, with deformation increasing significantly under higher osmotic pressures and reduced maximum bearing capacities. Failure deviatoric stress decreased exponentially with rising osmotic pressure. At constant confining pressure, elevated osmotic pressure shortened the total creep duration before failure. 2) Volumetric strain dilation occurred earlier under 23 MPa osmotic pressures compared to 1 MPa, which indicated accelerated crack initiation and unstable creep progression. 3) Accelerated creep rates manifested as nonlinear increases driven by coupled osmotic-deviatoric stress effects on crack damage evolution, which reflected macroscopic fracture propagation. 4) Under 7 MPa confining pressure, osmotic pressure reduced radial crack control and yielded stochastic failure modes such as composite fractures at 3 MPa. At 14 MPa confining pressure, macroscopic failure patterns became homogenized across osmotic pressures due to enhanced crack confinement. 5) SEM analysis revealed that tensile and shear fractures dominated microscale failure. Osmotic pressure raised intergranular reorganization, including fracture, refinement, and sliding, to form stress-adaptive microstructures. 6) Model-experiment consistency validated the applicability of the model across creep stages, including decay/steady-state creep at low stress and acceleration at high stress, which confirmed its capacity to characterize hydro-mechanical coupling effects. Conclusions 1) Osmotic pressure intensifies creep deformation and reduces long-term strength in carbonaceous mudstone, inducing failure under lower deviatoric stresses. 2) Under low confining pressures, osmotic pressure reduces radial crack confinement, increasing stochastic failure mo-des characterized by composite fracture patterns. 3) At the microscopic level, osmotic pressure alters fracture morphologies: tensile fractures exhibit scaly brittle surfaces and root-like patterns, whereas shear fractures display dimples and transgranular cracks. The microstructural fracture patterns correlate well with macroscopic failure modes. 4) The proposed nonlinear viscoelastic-plastic creep model effectively characterizes full-stage creep behavior under varying confining/osmotic pressures. The model parameters derived from laboratory tests produce theoretical curves that closely match the experimental data by considering confining pressure and osmotic pressure as variables.

Graphical abstract

关键词

碳质泥岩 / 水岩耦合 / 蠕变特性 / 破坏特征 / 蠕变模型

Key words

carbonaceous mudstone / hydro-mechanical coupling / creep characteristics / failure characteristics / creep model

本刊网刊

随着中国深部地下工程建设的快速发展，深部岩石力学问题日益受到学界关注^[1‒3]。其中，岩石蠕变特性直接影响隧道、边坡等工程的长期稳定性^[4‒6]。碳质泥岩作为西南山区常见岩层，具有强度低、遇水软化等特点，在工程中常面临地应力与渗透压的共同作用，导致其蠕变行为复杂化。因此，深入研究水岩耦合条件下碳质泥岩的蠕变特性，对保障地下工程安全具有重要意义。

室内试验是研究岩石材料蠕变特性的重要手段，现阶段学者们已经开展了广泛的研究。试验结果表明，影响岩石蠕变行为的重要因素包括岩石种类、地应力、围压，以及环境因素（水^[7‒8]、温度^[9‒11]和化学腐蚀^[12]）。张培森等^[7]通过温压控制试验揭示了深部砂岩渗透特性。李安润等^[8]基于Burgers模型建立了含水损伤蠕变本构模型。Bai等^[9]基于统计损伤理论建立了冻融循环损伤本构模型。张胜利等^[10]结合分数阶理论构建了考虑温度效应的盐岩蠕变损伤模型。陈国庆等^[11]通过三轴蠕变试验分析了不同砂岩在冻融作用下的孔隙演化与蠕变劣化规律。杨玉良等^[12]研究发现，钙芒硝盐岩蠕变应变随温度升高呈加速增长趋势。以上分析多针对温度影响，然而，水是影响隧道等地下工程施工最为关键的环境因素之一。现有研究普遍关注水对岩石蠕变特性的影响，但主要集中于静水压力条件下的蠕变行为，如不同含水率状态和干湿循环作用下的蠕变特性分析^[13‒15]。然而，工程实践表明，导致地下工程灾害频发的根本原因并非静水压力作用，而是工程地质环境中普遍存在的动态承压水系统（包括溶洞水、暗河及蓄水构造等）。这些承压水体的水压会随地下水位发生季节性波动，进而导致围岩产生显著变形，并引发支护结构长期开裂破坏等工程问题^[16‒17]。

由上述分析可知，水岩耦合蠕变试验研究对揭示地下工程围岩变形机理具有重要价值。朱珍德等^[18]发现渗透水压会显著促进脆性岩石的裂缝扩展和膨胀变形；Li等^[19]通过煤岩试验揭示了渗透率与损伤演化的相关性；Yan等^[20]在节理花岗岩试验中观察到渗透率与蠕变速率变化的同步性；Cerasi等^[21]探讨了页岩中CO₂封存过程的渗流愈合效应；Wang等^[22]研究发现循环水压会导致红砂岩蠕变曲线产生周期性波动。然而，当前岩石‒水力耦合蠕变研究仍存在明显不足：一方面，试验对象多局限于高渗透性软岩，对工程中广泛存在的碳质泥岩关注不足；另一方面，对于渗透压与应力耦合作用下的蠕变破坏机制，无论是宏观变形特征还是微观断裂过程，都缺乏系统深入的研究。

近年来，岩石蠕变本构模型研究取得了显著进展。作为表征岩石时变变形特性的重要工具，蠕变本构模型主要分为两类：一类是基于试验数据拟合的经验模型，多采用幂函数或指数函数形式^[23‒25]，这类模型虽能准确描述岩石的时效变形特征，但因参数过多而限制了其工程适用性；另一类是基于流变元件组合的经典模型（如Nishihara模型和Burgers模型）^[26‒28]，通过弹簧、阻尼器等元件的串联或并联来模拟岩石蠕变行为，但其仅能反映减速和稳态蠕变阶段，无法完整描述整个蠕变过程。齐亚静等^[26]通过引入应变触发的非线性黏壶对西原模型进行改进，建立了考虑三维蠕变的本构方程。丁靖洋等^[27]结合分数阶微积分和损伤理论，采用Abel黏壶替代传统Newton黏壶，构建了盐岩分数阶流变模型。苏腾等^[28]基于Scott-Blair元件和变系数分数阶元件，提出变阶分数阶非线性黏弹塑性蠕变模型，并实现了三维扩展。尽管上述学者对描述加速蠕变阶段的经典模型进行了改进，并取得了良好的效果，但这些模型大多描述了传统的蠕变过程，未将渗透压力作为一个重要指标引入到本构模型中，以描述水岩耦合作用下的整个蠕变过程。

基于此，本文以杨林隧道工程为背景，针对二衬开裂部位的碳质泥岩开展水岩耦合蠕变试验，系统研究其全过程的变形特征演化规律及宏微观断裂机制。基于流变力学与弹塑性力学理论，建立能完整描述3阶段蠕变特性的水岩耦合非线性黏弹塑性本构模型，并提出相应的参数确定方法。研究表明，试验数据与理论模型预测结果具有良好的一致性，验证了该模型在表征碳质泥岩蠕变行为方面的准确性与适用性，为类似工程条件下的围岩稳定性分析提供了理论依据。

1 试验及结果分析

1.1 工程背景

杨林隧道左幅ZK14+712.0～ZK14+720.5段二衬施工完成后一年多的时间内，在地应力及水压的共同作用下，围岩发生明显的流变变形。左侧拱肩至拱顶环向偏右约10 m范围内二衬混凝土发生不同程度的开裂和剥落，裂缝周边混凝土呈碎裂状，二衬钢筋弯曲外露，并有环向发展趋势。其中，ZK14+716.0里程处二衬开裂更加明显，掉块严重。专家决定对ZK14+712.0～ZK14+720.5段的初支及二衬全部凿除，并在进行支护优化后重新施作初支和二衬。考虑到二衬开裂段存在高水位和大水压，并且位于隧道拱顶偏左处，该区域主要由碳质泥岩组成，因此对水岩耦合作用下碳质泥岩的蠕变进行了试验研究。现场二衬开裂位置如图1所示。

1.2 岩样获取及试验设备

从杨林隧道二衬开裂处采集碳质泥岩样本，制备成直径为50 mm、高为100 mm的标准岩样。X射线衍射分析矿物组成（质量占比）：伊利石28.3%、绿泥石30.7%、石英24.4%、斜长石9.5%、白云石7.1%。试验设备采用微机岩石综合控制系统（图2（a）），渗透水压采用伺服加载装置（图2（b））。

1.3 试验荷载设置

根据杨林隧道地勘数据，试验围压设为7和14 MPa以匹配实测条件；试验渗透压设为1、2、3 MPa，以模拟工程渗透条件。首先，进行常规力学试验获得峰值应力（表1），为蠕变试验荷载分级提供依据。

依据《水利水电工程岩石试验规程》（SL/T 264—2020）分级加载原则，碳质泥岩常规及水岩耦合蠕变试验以峰值应力的90%为最大应力，荷载均分6级；每级加载240 min，未破坏则续加载至峰值应力。具体荷载分级见表2。

1.4 试验结果

1.4.1 蠕变特性

图3为碳质泥岩应变时程曲线。由图3可知，应变时程曲线呈现典型的阶梯式增长特征，整体光滑连续。在分级加载过程中，试样首先产生瞬时弹性变形，随后依次经历减速蠕变和稳态蠕变阶段。当荷载增至临界值时，碳质泥岩进入加速蠕变阶段并最终破坏。值得注意的是，在渗透压作用下，虽然应变曲线仍保持阶梯式上升趋势，但随着渗透压增大，各荷载水平下的蠕变量显著增加且破坏荷载明显降低。试验结果表明，渗透压不仅会加剧碳质泥岩的蠕变变形，还会显著降低其长期强度，使得碳质泥岩在相对较低的应力水平下即发生蠕变破坏。

图4（a）为碳质泥岩在不同围压下破坏偏应力随渗透压的变化规律。由图4（a）可知，破坏偏应力与渗透压呈显著负相关，其衰减趋势符合指数函数特征。以7 MPa围压工况为例，当渗透压由0递增至3 MPa时，对应的破坏偏应力从56依次降至50、44 MPa，整体降幅达21.4%。

图4（b）为围压和渗透压对蠕变时长t_c的耦合影响。由图4（b）可知，在相同围压条件下，渗透压的增大会显著缩短蠕变破坏历时。以7 MPa围压为例，渗透压从0增至3 MPa时，破坏时间由1 322.4 min锐减至729.6 min，时间缩短幅度达44.8%，充分体现了渗透压对碳质泥岩蠕变破坏的促进作用。

碳质泥岩的三轴蠕变行为与其内部裂隙的萌生、扩展及最终破坏特征密切相关，这一演化过程直接反映在体积应变ɛ_v的变化规律上，如图5（a）所示。由图5（a）可知：在7 MPa围压和1 MPa渗透压条件下，前3级荷载（至38 MPa）作用下，体积应变呈阶梯式增长，在第3级荷载结束时达到峰值7.99×10^-4（A点）。此后体积应变转为下降趋势，表明A点是体积变形由压缩转为膨胀的临界点，即扩容点。值得注意的是，当渗透压增至2 MPa和3 MPa时，扩容点出现时间提前，而且总蠕变时长显著缩短（由1 MPa时的1 322.4 min到3 MPa时的729.6 min），证实了渗透压会促进裂隙发育并加速蠕变破坏进程。图5（b）呈现相似的规律，不再赘述。

图6为最后一级偏应力下的蠕变速率-时间曲线（对数坐标）。由图6可知，蠕变过程可分为3个典型阶段：初始减速阶段（速率递减）、中期稳态阶段（速率波动但总体稳定）和末期加速阶段（速率急剧上升）。其中：稳态阶段持续时间最长；加速阶段虽仅占全过程的8%～12%，却贡献了约60%的总变形量，其快速变形特征与工程中围岩失稳破坏现象高度相关。值得注意的是，加速蠕变本质上是岩体内部裂隙网络贯通演化的宏观表征，在渗透压-偏应力耦合作用下，裂隙损伤累积导致变形速率呈非线性增长。

对比不同围压条件下的蠕变总时长发现：7 MPa围压下，渗透压每增加1 MPa可使总蠕变时长缩短约30%；在14 MPa围压下该规律不明显，这是由于高围压的约束效应抑制了渗透压的弱化作用。但试验数据证实，渗透压始终会显著提升碳质泥岩的长期蠕变速率（增幅达40%～65%），并使其时效强度降低18%～25%。

1.4.2 裂纹宏观破裂特征

试验研究表明，碳质泥岩的宏观破坏模式受围压和渗透压的协同控制。在7 MPa围压下，试样呈现典型的陡倾角剪切破坏特征，表现为上下端面形成有明显倾斜裂缝。当围压提升至14 MPa时，破坏模式转变为拉剪复合型，主剪切面贯穿试样的同时伴随有垂直劈裂裂纹。图7（a）、（b）分别为在7、14 MPa围压下碳质泥岩的破坏特征演变。

由图7（a）可知：渗透压为1 MPa时，形成缓倾角剪切破坏并伴随端部劈裂；渗透压为2 MPa时，出现缓倾角与陡倾角剪切复合破坏；渗透压为3 MPa时，因围压约束不足导致试样整体破碎，呈现共轭剪切破坏特征。由图7（b）可知：渗透压为1 MPa时，为陡倾角剪切破坏；渗透压为2 MPa时，转变为缓倾角剪切破坏；渗透压为3 MPa时，则出现类似单轴条件下的张拉破坏，但破坏程度较单轴工况轻微，这源于高围压对渗透压弱化效应的抑制作用。

综上可知：低围压（7 MPa）条件下，渗透压会显著增加破坏模式的随机性，表现为多种复合破坏形态；高围压（14 MPa）条件下能有效约束渗透压的影响，使破坏模式保持单一性。

1.4.3 裂纹微观破裂特征

针对岩石蠕变破坏断口细观形貌观测的既有研究较少，水岩耦合蠕变破坏的细观研究更鲜见^[29‒30]。鉴于围压与渗透压对蠕变宏观破坏形态的显著影响，本文基于不同围压、渗透压下碳质泥岩蠕变的细观破坏特征，揭示其水岩耦合蠕变破坏机制。图8为不同围压、渗透压下碳质泥岩的断口细观形貌。由图8可知：拉裂形态呈现鳞片状脆性断口、台阶花样和根状花样；剪裂形态则表现为韧窝、穿晶裂纹及滑移分离。

图8（a）中：渗透压为1 MPa时，断口分布大量剪切韧窝，伴有穿晶裂纹和滑移分离，同时存在鳞片状脆性断口（拉裂特征），这表明宏观的缓倾角剪切破坏与拉剪复合破坏与细观形貌一致；渗透压为2 MPa时，断口全为剪切特征（韧窝簇、滑移分离、穿晶裂纹），与宏观破坏形态对应；渗透压为3 MPa时，断口同时存在剪切韧窝簇和张拉台阶花样，与宏观仅见的共轭剪切破坏不符，这是因为较低围压（7 MPa）对裂隙控制作用弱，加之渗透压增加了随机性，导致未能呈现完整剪切破坏，反而出现少量张拉破坏特征。

图8（b）中：渗透压为1 MPa时，断口呈现剪切韧窝、滑移分离及张拉台阶花样，与宏观的陡倾角剪切破坏不符；渗透压为2 MPa时，明显有滑移分离、穿晶裂纹和部分韧窝特征，与宏观缓倾角剪切破坏一致；渗透压为3 MPa时，断口兼具剪切韧窝簇和张拉根状花样，与宏观张拉破坏特征不符，这是因为14 MPa的围压限制了次生裂纹发展，导致宏观剪切特征不明显，但细观仍保留部分韧窝特征。

综上，拉伸与剪切断裂是碳质泥岩微尺度断裂的主要类型。SEM图像清晰显示：在长期应力作用下，其微观结构发生胶结破坏，颗粒分布持续调整；渗透压促进了矿物颗粒间的断裂、细化、滑移及充填效应，形成新的承重结构以承受外部载荷。

2 水岩耦合非线性黏弹塑性蠕变模型

2.1 流变本构关系

随着时间推移，作用在岩石上的荷载会使岩石内部产生损伤裂痕，当长期荷载持续不变，这些裂纹就会不断发展，直至岩石发生不可恢复的破坏。岩石类型与外界条件的不同，加速破坏发生的时间节点亦随之变化。加速蠕变开始时所对应的时间和应变对于模型的参数识别影响很大，如果全程采用非线性方程，模型的参数辨识会非常困难，所以既要不断推广非线性模型，也要使得模型对于各种不同情况下的加速蠕变有较好的适用性。

通过引入一个带应变触发的非线性黏壶来描述岩石加速蠕变阶段的变形^[31]，如图9所示。图9中，η_nl为非线性黏壶黏滞系数，n为指数，ɛ_a为屈服应变。

该元件本构关系如下：

（1）

式中，\sigma为应力，\epsilon为应变，ɛ_nl为非线性黏壶应变，t为时间。

将式（1）第1式进行Laplace变换，可得：

\tilde{\sigma }={\eta }_{\mathrm{n}\mathrm{l}}\frac{S^{n+\mathrm{1}}}{n!}{\tilde{\epsilon }}_{\mathrm{n}\mathrm{l}}(s)（2）

式中，\tilde{\sigma }为\sigma的拉氏变换，{\tilde{\epsilon }}_{\mathrm{n}\mathrm{l}}为的{\epsilon }_{\mathrm{n}\mathrm{l}}拉氏变换，S为拉式空间变换的复变量，s为复频率变量。

将式（2）进行Laplace逆变换，可得

{\epsilon }_{\mathrm{n}\mathrm{l}}=\frac{\sigma }{{\eta }_{\mathrm{n}\mathrm{l}}}{\tau }^n,\epsilon \ge {\epsilon }_{\mathrm{a}}（3）

式中，\tau =\frac{t-t_{\mathrm{F}}}{t_{\mathrm{u}\mathrm{n}}}，t_F为岩石进入加速蠕变阶段的时间，t_{\mathrm{F}}={\left.t\right|}_{\epsilon ={\epsilon }_{\mathrm{a}}}，t_un为用于消除量纲的恒定时间，本文取t_un=1 h。

将图9中的元件与塑性元件并联组成一个非线性加速流变启动元件，然后并入西原模型中，形成新非线性黏弹塑性模型，如图10所示。图10中，σ为总应力，σ₁、ε₁分别为弹性应力、应变，σ₂、ε₂分别为黏弹性应力、应变，σ₃、ε₃分别为黏塑性应力、应变，σ₄为非线性应力，E₀、E₁和η₁、η₂分别为该模型的弹性和黏性参数，σ_s1和σ_s2分别为黏塑性蠕变阈值和加速蠕变阈值。

整个模型状态方程如下：

\left\{\begin{array}{l}\sigma ={\sigma }_{\mathrm{1}}={\sigma }_{\mathrm{2}}={\sigma }_{\mathrm{3}}={\sigma }_{\mathrm{4}},\\ \epsilon ={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}+{\epsilon }_{\mathrm{3}}+{\epsilon }_{\mathrm{n}\mathrm{l}},\\ {\sigma }_{\mathrm{1}}=E_{\mathrm{0}}{\epsilon }_{\mathrm{1}},\\ {\sigma }_{\mathrm{2}}=E_{\mathrm{1}}{\epsilon }_{\mathrm{2}}+{\eta }_{\mathrm{1}}{\dot{\epsilon }}_{\mathrm{2}},\\ {\sigma }_{\mathrm{3}}-{\sigma }_{\mathrm{s}\mathrm{1}}={\eta }_{\mathrm{2}}{\dot{\epsilon }}_{\mathrm{3}},\\ {\sigma }_{\mathrm{4}}-{\sigma }_{\mathrm{s}\mathrm{2}}=\frac{{\eta }_{\mathrm{n}\mathrm{l}}{t}_{\mathrm{u}\mathrm{n}}^n}{(t-t_{\mathrm{F}}{)}^n}{\epsilon }_{\mathrm{n}\mathrm{l}}\end{array}\right.（4）

根据式（4）和图10可知流变模型本构方程表达如下：

1）当σ＜σ_s1时，模型为黏弹性三参量模型，流变本构为：

(E_{\mathrm{0}}+E_{\mathrm{1}})\sigma +{\eta }_{\mathrm{1}}\dot{\sigma }=E_{\mathrm{0}}{E}_{\mathrm{1}}\epsilon +E_{\mathrm{0}}{\eta }_{\mathrm{1}}\dot{\epsilon }（5）

2）当σ_s1\leσ＜σ_s2或ε\leε_a时，模型为黏弹塑性西原模型，流变本构为：

{\eta }_{\mathrm{1}}\ddot{\epsilon }+E_{\mathrm{1}}\dot{\epsilon }=\frac{{\eta }_{\mathrm{1}}}{E_{\mathrm{0}}}\ddot{\sigma }+\left(\mathrm{1}+\frac{E_{\mathrm{1}}}{E_{\mathrm{0}}}+\frac{{\eta }_{\mathrm{1}}}{{\eta }_{\mathrm{2}}}\right)\dot{\sigma }+\frac{E_{\mathrm{1}}(\sigma -{\sigma }_{\mathrm{s}\mathrm{1}})}{{\eta }_{\mathrm{2}}}（6）

3）当σ\geσ_s2且ε＞ε_a时，模型为非线性黏弹塑性，此时模型状态方程为式（4），其本构方程为：

\begin{array}{l}{\eta }_{\mathrm{1}}\ddot{\epsilon }+E_{\mathrm{1}}\dot{\epsilon }=\left[\frac{{\eta }_{\mathrm{1}}}{E_{\mathrm{0}}}+\frac{{\eta }_{\mathrm{1}}}{{\eta }_{\mathrm{n}\mathrm{l}}}{\left(\frac{t-t_{\mathrm{F}}}{t_{\mathrm{u}\mathrm{n}}}\right)}^n\right]\ddot{\sigma }+\\ \left[\mathrm{1}+\frac{E_{\mathrm{1}}}{E_{\mathrm{0}}}+\frac{{\eta }_{\mathrm{1}}}{{\eta }_{\mathrm{2}}}+\frac{E_{\mathrm{1}}}{{\eta }_{\mathrm{n}\mathrm{l}}}{\left(\frac{t-t_{\mathrm{F}}}{t_{\mathrm{u}\mathrm{n}}}\right)}^n\right]\dot{\sigma }+\frac{E_{\mathrm{1}}(\sigma -{\sigma }_{\mathrm{s}\mathrm{1}})}{{\eta }_{\mathrm{2}}}\end{array}（7）

2.2 一维蠕变方程

1）当σ＜σ_s1时，模型为黏弹性三参量模型，组成该模型满足的基本条件如式（8）所示：

\left\{\begin{array}{l}\sigma ={\sigma }_{\mathrm{1}}={\sigma }_{\mathrm{2}},\\ \epsilon ={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}},\\ {\sigma }_{\mathrm{1}}=E_{\mathrm{0}}{\epsilon }_{\mathrm{1}},\\ {\sigma }_{\mathrm{2}}=E_{\mathrm{1}}{\epsilon }_{\mathrm{2}}+{\eta }_{\mathrm{1}}{\dot{\epsilon }}_{\mathrm{2}}\end{array}\right.（8）

对式（8）进行Laplace变换得

\tilde{\epsilon }=\frac{\sigma }{E_{\mathrm{0}}S}+\frac{\sigma }{S(E_{\mathrm{1}}+{\eta }_{\mathrm{1}}S)}（9）

将式（9）进行Laplace逆变换得到蠕变公式为：

\epsilon =\frac{\sigma }{E_{\mathrm{0}}}+\frac{\sigma }{E_{\mathrm{1}}}\left(\mathrm{1}-{\mathrm{e}}^{-\frac{E_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)（10）

2）当σ_s1\leσ＜σ_s2或ε\leε_a时，模型为西原模型，组成该模型满足的基本条件如式（11）所示：

\left\{\begin{array}{l}\sigma ={\sigma }_{\mathrm{1}}={\sigma }_{\mathrm{2}}={\sigma }_{\mathrm{3}},\\ \epsilon ={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}+{\epsilon }_{\mathrm{3}},\\ {\sigma }_{\mathrm{1}}=E_{\mathrm{0}}{\epsilon }_{\mathrm{1}},\\ {\sigma }_{\mathrm{2}}=E_{\mathrm{1}}{\epsilon }_{\mathrm{2}}+{\eta }_{\mathrm{1}}{\dot{\epsilon }}_{\mathrm{2}},\\ {\sigma }_{\mathrm{3}}-{\sigma }_{\mathrm{s}\mathrm{1}}={\eta }_{\mathrm{2}}{\dot{\epsilon }}_{\mathrm{3}}\end{array}\right.（11）

对式（11）进行Laplace变换得

\tilde{\epsilon }=\frac{\sigma }{E_{\mathrm{0}}S}+\frac{\sigma }{S(E_{\mathrm{1}}+{\eta }_{\mathrm{1}}S)}+\frac{\sigma -{\sigma }_{\mathrm{s}\mathrm{1}}}{{\eta }_{\mathrm{2}}{S}^{\mathrm{2}}}（12）

将式（12）进行Laplace逆变换得到蠕变公式：

\epsilon =\frac{\sigma }{E_{\mathrm{0}}}+\frac{\sigma }{E_{\mathrm{1}}}\left(\mathrm{1}-{\mathrm{e}}^{-\frac{E_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)+\frac{\sigma -{\sigma }_{\mathrm{s}\mathrm{1}}}{{\eta }_{\mathrm{2}}}t（13）

3）当σ\geσ_s2且ε＞ε_a时，模型为非线性黏弹塑性，组成该模型满足的基本条件如式（4）所示，对式（4）进行Laplace变换得：

\tilde{\epsilon }=\frac{\sigma }{E_{\mathrm{0}}S}+\frac{\sigma }{S(E_{\mathrm{1}}+{\eta }_{\mathrm{1}}S)}+\frac{\sigma -{\sigma }_{\mathrm{s}\mathrm{1}}}{{\eta }_{\mathrm{2}}{S}^{\mathrm{2}}}+\frac{\sigma n!}{{\eta }_{\mathrm{n}\mathrm{l}}{S}^{n+\mathrm{1}}}（14）

将式（14）进行Laplace逆变换得到蠕变公式为：

\epsilon =\frac{\sigma }{E_{\mathrm{0}}}+\frac{\sigma }{E_{\mathrm{1}}}\left(\mathrm{1}-{\mathrm{e}}^{-\frac{E_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)+\frac{\sigma -{\sigma }_{\mathrm{s}\mathrm{1}}}{{\eta }_{\mathrm{2}}}t+\frac{\sigma }{{\eta }_{\mathrm{n}\mathrm{l}}}{\left(\frac{t-t_{\mathrm{F}}}{t_{\mathrm{u}\mathrm{n}}}\right)}^n（15）

2.3 三维蠕变方程

对于三维本构关系的推导，孙钧^[32]作出了如下假设：1）岩石的体积蠕变影响较小；2）将岩石视为各向同性材料；3）泊松比不随蠕变过程变化。

在三维应力状态下，岩石内部应力张量由球应力张量σ_m和偏应力张量s_ij 组成，表达式为：

\left\{\begin{array}{l}{\sigma }_{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}})=\frac{\mathrm{1}}{\mathrm{3}}{\sigma }_{kk},\\ s_{ij}={\sigma }_{ij}-{\delta }_{ij}{\sigma }_{\mathrm{m}}={\sigma }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\sigma }_{kk}\end{array}\right.（16）

式中，σ_kk 为正应力分量，δ_ij 为克朗内克符号，σ_ij 为应力张量。

所以

{\sigma }_{ij}=s_{ij}+{\delta }_{ij}{\sigma }_m（17）

式中，s_ij 为偏应力张量。

同理，岩石内部应变张量由球应变张量ε_m和偏应变张量e_ij 组成，表达式为：

\left\{\begin{array}{l}{\epsilon }_{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}({\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}+{\epsilon }_{\mathrm{3}})=\frac{\mathrm{1}}{\mathrm{3}}{\epsilon }_{kk},\\ e_{ij}={\epsilon }_{ij}-{\delta }_{ij}{\epsilon }_{\mathrm{m}}={\epsilon }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\epsilon }_{kk}\end{array}\right.（18）

所以

{\epsilon }_{ij}=e_{ij}+{\delta }_{ij}{\epsilon }_{\mathrm{m}}（19）

式中，{\epsilon }_{ij}^{}为应变张量。

对于三维应力条件下，该非线性模型的总应变表达式为：

\epsilon ={\epsilon }_{ij}^{\mathrm{e}}+{\epsilon }_{ij}^{\mathrm{v}\mathrm{e}}+{\epsilon }_{ij}^{\mathrm{v}\mathrm{p}}+{\epsilon }_{ij}^{\mathrm{n}\mathrm{l}}（20）

式中，{\epsilon }_{ij}^{\mathrm{e}}、{\epsilon }_{ij}^{\mathrm{v}\mathrm{e}}、{\epsilon }_{ij}^{\mathrm{v}\mathrm{p}}、{\epsilon }_{ij}^{\mathrm{n}\mathrm{l}}分别为三维弹性、黏弹性、黏塑性和非线性黏壶元件应变。

在三维本构关系中，用弹性模量为E₀和泊松比μ分别表示岩石剪切模量G₀和体积模量K₀，即：

G_{\mathrm{0}}=\frac{E_{\mathrm{0}}}{\mathrm{2}(\mathrm{1}+\mu )}（21）

K_{\mathrm{0}}=\frac{E_{\mathrm{0}}}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mu )}（22）

根据三维应力状态下的广义胡克定律，三维弹性体本构关系为：

\left\{\begin{array}{l}{e}_{ij}=\frac{s_{ij}}{\mathrm{2}{G}_{\mathrm{0}}},\\ {\epsilon }_{kk}=\frac{{\sigma }_{kk}}{\mathrm{3}{K}_{\mathrm{0}}}\end{array}\right.（23）

所以，三维状态下弹性体应变可以表示为：

{\epsilon }_{ij}^{\mathrm{e}}=\frac{s_{ij}}{\mathrm{2}{G}_{\mathrm{0}}}+\frac{{\sigma }_{\mathrm{m}}{\delta }_{ij}}{\mathrm{3}{K}_{\mathrm{0}}}（24）

根据第2.3节的假设，岩石流变过程的特性主要为剪切变形方面，也就是形状的变化，所以对于三维黏弹性本构关系有

{\epsilon }_{ij}^{\mathrm{v}\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}{G}_{\mathrm{1}}}\left(\mathrm{1}-{\mathrm{e}}^{-\frac{G_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)s_{ij}（25）

式中，G₁为黏弹性剪切模量，η₁为剪切黏滞系数。

三维黏塑性体本构关系为：

{\epsilon }_{ij}^{\mathrm{v}\mathrm{p}}=\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}〈\varphi \left(\frac{F}{F_{\mathrm{0}}}\right)〉\frac{\partial Q}{\partial {\sigma }_{ij}}t（26）

式中，ϕ为幂指数取1时的幂函数，F、F₀分别为岩石的屈服函数和初始屈服函数，Q为塑性势函数。

当采用相关联流动法则，式（26）可变为：

{\epsilon }_{ij}^{\mathrm{v}\mathrm{p}}=\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}\left(\frac{F}{F_{\mathrm{0}}}\right)\frac{\partial F}{\partial {\sigma }_{ij}}t（27）

室内三轴蠕变试验中，当岩石达到加速蠕变阶段，主要表现为轴向应变极速增加，因此，可以选择第一主应变ε₁₁作为岩石是否达到应变阈值ε_a的判断条件，即根据ε₁₁和ε_a的大小判断黏壶参与流变与否。故非线性黏壶三维本构关系为：

（28）

将式（24）、（25）、（27）、（28）代入式（20），可得非线性黏弹塑性模型三维本构方程为：

（29）

岩石屈服准则多采用Mohr-Coulomb和von Mises屈服准则，但Mohr-Coulomb准则不能考虑中间主应力σ₂；von Mises准则忽略了球形应力对岩石蠕变特性的影响，尤其是软岩蠕变特性。因此，本文采用D-P屈服准则，可以弥补前两个屈服准则的不足。其表达式为：

F=\lambda I_{\mathrm{1}}+\sqrt[]{J_{\mathrm{2}}}-k（30）

式中，I₁为应力张量第一不变量，J₂为偏应力张量的第二不变量，λ和k为岩石参数。λ和k的表达式分别为：

\lambda =\frac{\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )}（31）

k=\frac{\mathrm{6}c\cdot \mathrm{c}\mathrm{o}\mathrm{s}\phi }{\sqrt[]{\mathrm{3}}(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )}（32）

式中，c和φ分别为岩石的黏聚力和内摩擦角。

对于等围压三轴压缩试验，围压{\sigma }_{\mathrm{2}}={\sigma }_{\mathrm{3}}，所以

\left\{\begin{array}{l}{s}_{\mathrm{11}}={\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{m}}=\frac{\mathrm{2}}{\mathrm{3}}\left({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}\right),\\ \sqrt[]{J_{\mathrm{2}}}=\frac{\mathrm{1}}{\sqrt[]{\mathrm{3}}}({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}),\\ \frac{\partial F}{\partial {\sigma }_{\mathrm{11}}}=\frac{\sqrt[]{\mathrm{3}}+\mathrm{3}\lambda }{\mathrm{3}}\end{array}\right.（33）

将式（33）代入式（29）中并令F₀=1，得到非线性黏弹塑性蠕变方程为^[33]：

（34）

第2.3节建立了不考虑渗透压的蠕变方程。本文聚焦渗透压对碳质泥岩蠕变特性的影响，需在方程中引入渗透压参数。自然环境中，水对岩石的作用涵盖力学、物理、化学及生物效应，而室内试验仅需考虑力学效应。渗透压作用具体表现为孔隙水压力与渗透力。孔隙水压力作为静水压力，在岩体内均匀分布，仅导致体积变形而不改变形状；渗透力为动水压力，源于渗流通道内水力梯度对岩骨架的作用，会引起形状改变。在碳质泥岩蠕变过程中，蠕变应力由孔隙水压力（静水压力）和有效应力（岩骨架承载）共同承担，其中，有效应力是诱发蠕变变形的主导因素，故有：

{\sigma }^{\text{'}}=\sigma -ap（35）

式中：{\sigma }^{\text{'}}为有效应力；σ为轴向压力或围压；a为Biot系数，0<a<1；p为渗透水压。根据式（35）可得轴压和围压的有效应力分别为：

{\sigma }_{\mathrm{1}}^{\text{'}}=\left\{\begin{array}{l}{\sigma }_{\mathrm{1}}-ap,t\le t_{\mathrm{1}};\\ {\sigma }_{\mathrm{1}}-a_{\mathrm{1}}p,{\sigma }_{\mathrm{0}}\ge {\sigma }_{\mathrm{s}},t>t_{\mathrm{1}}\end{array}\right.（36）

{\sigma }_{\mathrm{3}}^{\text{'}}={\sigma }_{\mathrm{3}}-ap（37）

式（36）、（37）中，a和a₁为Biot系数，t₁为发生加速蠕变的时间点。

将式（36）和（37）代入式（34）中，可得水岩耦合非线性黏弹塑性蠕变方程为：

（38）

2.4 模型参数确定

优化分析法中的最小二乘法对于非线性问题的拟合效果并不理想，其对于初值的依赖性比较高^[34]。所以，本文采用对于初值依赖性不强且不易收敛到局部极小值的L‒M非线性最小二乘法进行模型的参数识别。具体参数识别过程如下：

1）式（38）中，当t=0时，瞬时弹性变形为：

{\epsilon }_{\mathrm{11}}=\frac{{\sigma }_{\mathrm{1}}+\mathrm{2}{\sigma }_{\mathrm{3}}-\mathrm{3}ap}{\mathrm{9}{K}_{\mathrm{0}}}+\frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{3}{G}_{\mathrm{0}}}（39）

根据试验结果取3种情况下的瞬时弹性变形，可求解式（39）中的参数a、K₀和G₀。

2）η₂的求解可根据试验确定岩石从衰减蠕变阶段进入等速蠕变阶段，该应变之后蠕变曲线斜率为k，即：

k=\frac{\mathrm{3}\lambda +\sqrt[]{\mathrm{3}}}{\mathrm{3}{\eta }_{\mathrm{2}}}（40）

根据式（40）可得

{\eta }_{\mathrm{2}}=\frac{\mathrm{3}\lambda +\sqrt[]{\mathrm{3}}}{\mathrm{3}k}（41）

3）除上述4个参数外，式（38）还包括参数η₁、G₁、a₁、η_nl和n，可采用L-M非线性最小二乘法进行非线性回归求解剩余参数。

2.5 长期强度的确定

蠕变模型分析表明：当应力水平低于σ_s时，岩石呈衰减蠕变和等速蠕变；当应力超过σ_s时，渗透压与应力的协同作用会促使内部裂隙发展，导致等速蠕变向加速蠕变转变。故精确获取蠕变参数需先确定σ_s值，再进行方程参数拟合。

岩体长期强度作为时间依赖性强度指标，表征无限时域蠕变引发破坏的最小应力阈值，其值显著低于瞬时强度。该参数对隧道工程长期稳定性评估至关重要。现行长期强度确定方法主要有3类：等时应力-应变曲线簇法、过渡蠕变法及稳态黏塑性蠕变速率-应力函数曲线法。其中：等时曲线簇法^[35]虽应用广泛，但其拐点辨识困难导致精度受限；过渡蠕变法^[35]仅能给出应力区间，难以精确量化。基于此，本文选用稳态黏塑性蠕变速率-应力函数曲线法进行测定。

限于篇幅，本文仅展示在7 MPa围压下常规三轴蠕变下的拟合曲线。根据蠕变试验结果，拟合得到稳态黏塑性蠕变速率-应力的函数关系，如图11所示，具体拟合公式如下：

（42）

根据式（42）计算得出，当围压σ₃=7 MPa时，碳质泥岩的长期强度σ_s为34 MPa。采用相同方法，获得了不同围压和渗透压条件下的长期强度值，见表3。从表3可以看出，长期强度与瞬时强度的比值（σ_s/σ_c）随渗透压增大而递减。具体而言，无渗透压时，σ_s/σ_c为0.520～0.580；施加渗透压后，该比值降至0.410～0.540，证实渗透压对长期强度具有显著弱化效应。碳质泥岩长期强度随渗透压升高呈指数型衰减趋势，如图12所示。该规律可通过散点数据的非线性拟合得以验证。

3 模型本构验证及分析

为了验证模型的正确性，结合上述不同水压下的蠕变试验，利用σ₃=7 MPa、p=2 MPa和σ₃=14 MPa、p=1 MPa的试验结果对参数进行确定。采用第2.5节给出的参数求取方案，分别通过计算得到参数a、K₀、G₀和η₂，然后通过L-M非线性最小二乘法进行非线性回归，得到η₁、G₁、a₁、η_nl和n。同时，将本文模型与文献[36]中所给出的水岩耦合西原模型进行对比，结果如图13所示，模型参数分别如表4和5所示。由图13可知，试验曲线和理论曲线的吻合程度非常高。这一结果得益于建立的水岩耦合非线性黏弹塑性蠕变模型能够全面反映岩石不同阶段的蠕变特性，并且光滑的试验曲线与模型理论曲线基本一致。同时，该模型也能准确描述低应力水平下碳质泥岩的衰减、等速蠕变特性，并且可以描述碳质泥岩在高应力状态下的加速蠕变特征，从而验证了该模型的适用性。同时，通过同传统的水岩耦合西原模型对比，可以看出，本文提出的模型计算值同试验值相关系数高达0.979，而西元模型相关系数仅为0.795，因此本文提出的模型可以很好地呈现加速蠕变过程，更具优势。

4 结 论

本文以杨林隧道在地应力和水压耦合作用下围岩流变变形引发的二衬开裂为工程背景，针对开裂区的碳质泥岩，开展室内水岩耦合蠕变试验，在模拟现场地应力和渗透压条件下，深入探讨了围压及渗透压对其蠕变行为与宏微观破裂规律的影响。基于渗透压和偏应力作用，构建了水岩耦合非线性黏弹塑性蠕变模型，并给出了参数确定方法。结论如下：

1）渗透压作用可显著增强碳质泥岩的蠕变变形，同时削弱其长期强度，导致在较低荷载水平下诱发蠕变破坏。

2）渗透压对碳质泥岩蠕变破坏形态的影响尤为显著。低围压时，其弱化对裂隙的径向约束作用，使宏观破坏形态更随机，呈现多种组合类型。高围压下，围压控制作用增强，渗透压对裂隙发育随机性的影响减弱，宏观破坏形态趋于单一。

3）渗透压改变了碳质泥岩蠕变的细观断口形态，包括拉裂与剪裂。拉裂主要表现为鳞片状脆性断口、台阶花样和须状花样。剪裂则包含韧窝、穿晶裂纹及滑移分离。细观断口特征与宏观破坏形态基本对应，但某些差异源于渗透压影响裂隙发育的随机性。

4）在长期应力作用下，碳质泥岩微观结构胶结破坏，颗粒分布持续调整。渗透压促进了矿物颗粒的断裂、细化、滑移与充填效应，形成新的承重结构以承受外部载荷。

5）所建水岩耦合非线性黏弹塑性蠕变模型能反映不同岩石在围压及渗透压下的蠕变全过程特征。模型以围压和渗透压为变量，参数确定简便且可通过室内试验获取。计算曲线与试验数据吻合良好，具备良好适用性。

参考文献

[1] Wang Bo, Tian Zhiyin, Ma Shiji,et al.Experimental study on micro-damage of rock mass in sensitive neighborhood under rheological disturbance[J].Chinese Journal of Rock Mechanics and Engineering,2024,43(S2):3820‒3831.

[2] 王波,田志银,马世纪,等.流变扰动条件下敏感邻域内岩体微观损伤试验研究[J].岩石力学与工程学报,2024,43(S2):3820‒3831.

[3] Miao Shengjun, Yang Pengjin, Wang Hui,et al.Study on fatigue rheological damage model of siltstone under cyclic loading[J].Engineering Mechanics,2022,39(7):70‒80. doi:10.6052/j.issn.1000-4750.2021.03.0235

[4] 苗胜军,杨鹏锦,王辉,等.循环荷载作用下粉砂岩疲劳流变损伤模型研究[J].工程力学,2022,39(7):70‒80. doi:10.6052/j.issn.1000-4750.2021.03.0235

[5] Xu Jiang, Yang Hongwei, Peng Shoujian,et al.Experimental study of mechanical property of sandstone under pore water pressure and confining pressure[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(8):1618‒1623.

[6] 许江,杨红伟,彭守建,等.孔隙水压力‒围压作用下砂岩力学特性的试验研究[J].岩石力学与工程学报,2010,29(8):1618‒1623.

[7] Chu Zhaofei, Wu Zhijun, Liu Quansheng,et al.Analytical solution for lined circular tunnels in deep viscoelastic Bu-rgers rock considering the longitudinal discontinuous excavation and sequential installation of liners[J].Journal of Engineering Mechanics,2021,147(4):04021009. doi:10.1061/(asce)em.1943-7889.0001912

[8] Chen Weizhong, Lu Chen, Yu Hongdan,et al.Progress in long-term mechanical properties of clay rocks under thermo-hydro-mechanical coupling conditions[J].Chinese Journal of Rock Mechanics and Engineering,2021,40(2):233‒247.

[9] 陈卫忠,卢琛,于洪丹,等.黏土岩温度‒渗流‒应力耦合作用的长期力学特性研究若干进展[J].岩石力学与工程学报,2021,40(2):233‒247.

[10] Zhang Liangliang, Wang Xiaojian.Viscoelastic‒plastic damage creep model for rock[J].Chinese Journal of Geotechnical Engineering,2020,42(6):1085‒1092.

[11] 张亮亮,王晓健.岩石黏弹塑性损伤蠕变模型研究[J].岩土工程学报,2020,42(6):1085‒1092.

[12] Zhang Peisen, Zhao Chengye, Hou Jiqun,et al.Experimental study on seepage characteristics of deep sandstone under high temperature and different hydraulic pressures[J].Chinese Journal of Rock Mechanics and Engineering,2020,39(6):1117‒1128. doi:10.13722/j.cnki.jrme.2019.1113

[13] 张培森,赵成业,侯季群,等.高温与不同水压下深部砂岩渗透特性试验研究[J].岩石力学与工程学报,2020,39(6):1117‒1128. doi:10.13722/j.cnki.jrme.2019.1113

[14] Li Anrun, Deng Hui, Wang Hongjuan,et al.Constitutive mo-del of water-damaged silty mudstone under water-rock interactions[J].Hydrogeology & Engineering Geology,2021,48(2):106‒113.

[15] 李安润,邓辉,王红娟,等.水‒岩作用下粉砂质泥岩含水损伤本构模型[J].水文地质工程地质,2021,48(2):106‒113.

[16] Bai Yao, Shan Renliang, Ju Yang,et al.Study on the mechanical properties and damage constitutive model of frozen weakly cemented red sandstone[J].Cold Regions Science and Technology,2020,171:102980. doi:10.1016/j.coldregions.2019.102980

[17] Zhang Shengli, Liang Weiguo, Xiao Ning,et al.Fractional viscoelastic‒plastic creep damage model of salt rock considering temperature[J].Chinese Journal of Rock Mechanics and Engineering,2022,41(S2):3198‒3209.

[18] 张胜利,梁卫国,肖宁,等.考虑温度的盐岩分数阶黏弹塑性蠕变损伤模型[J].岩石力学与工程学报,2022,41(S2):3198‒3209.

[19] Chen Guoqing, Wan Yi, Sun Xiang,et al.Research on creep behaviors and fractional order damage model of sandstone subjected to freeze‒thaw cycles in different temperature ranges[J].Chinese Journal of Rock Mechanics and Engineering,2021,40(10):3198‒3209.

[20] 陈国庆,万亿,孙祥,等.不同温差冻融后砂岩蠕变特性及分数阶损伤模型研究[J].岩石力学与工程学报,2021,40(10):1962‒1975.

[21] Yang Yuliang, Liang Weiguo, Yang Xiaoqin,et al.Temperature effect on creep of glauberite salt rock under multi-field coupling[J].Journal of China Coal Society,2020,45(3):1070‒1080. doi:10.13225/j.cnki.jccs.2019.0238

[22] 杨玉良,梁卫国,杨晓琴,等.钙芒硝盐岩多场耦合作用下蠕变的温度效应研究[J].煤炭学报,2020,45(3):1070‒1080. doi:10.13225/j.cnki.jccs.2019.0238

[23] Jiang T, Shao J F, Xu W Y,et al.Experimental investigation and micromechanical analysis of damage and permeability variation in brittle rocks[J].International Journal of Rock Mechanics and Mining Sciences,2010,47(5):703‒713. doi:10.1016/j.ijrmms.2010.05.003

[24] Wang H L, Xu W Y, Jia C J,et al.Experimental research on permeability evolution with microcrack development in sandstone under different fluid pressures[J].Journal of Geotechnical and Geoenvironmental Engineering,2016,142(6):04016014. doi:10.1061/(ASCE)GT.1943-5606.0001462

[25] Yang T H, Jia P, Shi W H,et al.Seepage‒stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway[J].Tunnelling and Underground Space Technology,2014,43:11‒19. doi:10.1016/j.tust.2014.03.005

[26] Yu Mingyuan, Liu Baoguo, Chu Zhaofei,et al.Permeability,deformation characteristics,and damage constitutive model of shale under triaxial hydromechanical coupling[J].Bulletin of Engineering Geology and the Environment,2022,81(3):85. doi:10.1007/s10064-022-02584-4

[27] Yu Mingyuan, Liu Baoguo, Liu Kaiyun,et al.Creep behavior of carbonaceous mudstone under triaxial hydraulic co-upling condition and constitutive modelling[J].International Journal of Rock Mechanics and Mining Sciences,2023,164:105357. doi:10.1016/j.ijrmms.2023.105357

[28] Zhu Zhende, Xu Weiya, Zhang Aijun.Mechanism analysis and testing study on damage and fracture of brittle rock[J].Chinese Journal of Rock Mechanics and Engineering,2003,22(9):1411‒1416. doi:10.3321/j.issn:1000-6915.2003.09.002

[29] 朱珍德,徐卫亚,张爱军.脆性岩石损伤断裂机理分析与试验研究[J].岩石力学与工程学报,2003,22(9):1411‒1416. doi:10.3321/j.issn:1000-6915.2003.09.002

[30] Li S C, Chen Z Q, Miao X X.Experimental study on the properties of time-dependent deformation-seepage in water-saturated broken sandstone[J].Journal of Mining and Safety Engineering,2011,28:542‒547.

[31] Yan Lei, Yi Wenhua, Liu Liansheng,et al.Blasting-induced permeability enhancement of ore deposits associated with low-permeability weakly weathered granites based on the split Hopkinson pressure bar[J].Geofluids,2018,2018:4267878. doi:10.1155/2018/4267878

[32] Cerasi P, Lund E, Kleiven M L,et al.Shale creep as leakage healing mechanism in CO₂ sequestration[J].Energy Procedia,2017,114:3096‒3112. doi:10.1016/j.egypro.2017.03.1439

[33] Wang Qiang, Hao Shengwang.Creep behavior and permeability evolution of red sandstone in three Gorges Reservoir area subjected to cyclic seepage pressure[J].Geomechanics for Energy and the Environment,2023,36:100498. doi:10.1016/j.gete.2023.100498

[34] Wang Guijun.A new constitutive creep-damage model for salt rock and its characteristics[J].International Journal of Rock Mechanics and Mining Sciences,2004,41:61‒67. doi:10.1016/j.ijrmms.2004.03.020

[35] Yang Chunhe, Daemen J J K, Yin Jianhua.Experimental investigation of creep behavior of salt rock[J].International Journal of Rock Mechanics and Mining Sciences,1999,36(2):233‒242. doi:10.1016/S0148-9062(98)00187-9

[36] Main I G.A damage mechanics model for power-law creep and earthquake aftershock and foreshock sequences[J].Geophysical Journal International,2000,142(1):151‒161. doi:10.1046/j.1365-246x.2000.00136.x

[37] Qi Yajing, Jiang Qinghui, Wang Zhijian,et al.3d creep constitutive equation of modified nishihara model and its parameters identification[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(2):347‒355. doi:10.3969/j.issn.1000-6915.2012.02.014

[38] 齐亚静,姜清辉,王志俭,等.改进西原模型的三维蠕变本构方程及其参数辨识[J].岩石力学与工程学报,2012,31(2):347‒355. doi:10.3969/j.issn.1000-6915.2012.02.014

[39] Ding Jingyang, Zhou Hongwei, Chen Qiong,et al.Characters of rheological damage and constitutive model of salt rock[J].Rock and Soil Mechanics,2015,36(3):769‒776. doi:10.16285/j.rsm.2015.03.022

[40] 丁靖洋,周宏伟,陈琼,等.盐岩流变损伤特性及本构模型研究[J].岩土力学,2015,36(3):769‒776. doi:10.16285/j.rsm.2015.03.022

[41] Su Teng, Zhou Hongwei, Zhao Jiawei,et al.A creep model of rock based on variable order fractional derivative[J].Chinese Journal of Rock Mechanics and Engineering,2019,38(7):1355‒1363. doi:10.13722/j.cnki.jrme.2018.1382

[42] 苏腾,周宏伟,赵家巍,等.基于变阶分数阶导数的岩石蠕变模型[J].岩石力学与工程学报,2019,38(7):1355‒1363. doi:10.13722/j.cnki.jrme.2018.1382

[43] Zhou F, Sun W, Shao J,et al.Experimental study on nano silica modified cement base grouting reinforcement materials[J].Geomechanics Engineering,2020,20(1):67‒73.

[44] Yang X X, Sun D K, Jing H W.Morphological features of shear-formed fractures developed in a rock bridge[J].Engineering Geology,2020,278:105833. doi:10.1016/j.enggeo.2020.105833

[45] Liu Kaiyun, Xue Yongtao, Zhou Hui.Study on 3D nonlinear visco‒elastic‒plastic creep constitutive model with parameter unsteady of soft rock based on improved Bingham model[J].Rock and Soil Mechanics,2018,39(11):4157‒4164.

[46] 刘开云,薛永涛,周辉.基于改进Bingham模型的软岩参数非定常三维非线性黏弹塑性蠕变本构研究[J].岩土力学,2018,39(11):4157‒4164.

[47] 孙钧.岩土材料流变及其工程应用[M].北京:中国建筑工业出版社,1999.

[48] Nara Y, Tanaka M, Harui T.Evaluating long-term strength of rock under changing environments from air to water[J].Engineering Fracture Mechanics,2017,178:201‒211. doi:10.1016/j.engfracmech.2017.04.015

[49] Zhou Hongwei, Liu Di, Lei Gang,et al.The creep-damage model of salt rock based on fractional derivative[J].Energies,2018,11(9):2349. doi:10.20944/preprints201807.0584.v1

[50] Wu Fei, Zhang Hao, Zou Quanle,et al.Viscoelastic‒plastic damage creep model for salt rock based on fractional derivative theory[J].Mechanics of Materials,2020,150:103600. doi:10.1016/j.mechmat.2020.103600

[51] Yu Mingyuan, Liu Baoguo, Liu Kaiyun,et al.Creep behavior of carbonaceous mudstone under triaxial hydraulic coupling condition and constitutive modelling[J].International Journal of Rock Mechanics and Mining Sciences,2023,164:105357. doi:10.1016/j.ijrmms.2023.105357

基金资助

中央高校基本科研业务费专项资金项目(BSKYZX202543)

国家自然科学基金项目(41961134032)

国家自然科学基金项目(42577234)