基于损伤统计的砂岩三维蠕变本构模型

Three-dimensional Creep Constitutive Model of Sandstone Based on Damage Statistics

摘要

岩体蠕变对工程长期稳定性至关重要，岩石是岩体的组成部分，本构模型最能反映岩石蠕变过程的力学特征和变形规律，利用MTS‒815岩石刚性试验机和STAC600‒600型岩石三轴流变试验机开展砂岩常规三轴压缩试验和分级蠕变试验，借助SMP准则和Weibull分布函数建立损伤统计模型，并验证模型准确性。采用损伤改进西原模型并引入裂隙塑性体，建立砂岩三轴蠕变条件下的本构模型，对模型参数进行优化并与其他本构模型进行对比。结果表明：1）利用SMP准则与Weibull分布函数推导的损伤变量能够有效描述砂岩三轴压缩与三轴蠕变条件下的损伤演化过程，且损伤演化参数u、m与围压有着明显函数关系；2）将围压对本构模型相关变量的影响纳入构建的损伤变量后，模型拟合结果与砂岩分级蠕变过程能够高度契合；3）本构模型中的相关参数（包括泊松比等）均推广至三维应用形式，且裂隙塑性体的引入使模型不单适用于完整岩石。该模型在工程中具有一定实用性，对其他类型岩石也具有参考价值。

Abstract

Objective The creep of rock mass is very important for long-term project stability, and the constitutive model best reflects the mechanical characteristics and deformation patterns of the rock creep process. In recent years, extensive studies on creep characteristics of rocks have been conducted, and many creep constitutive models have been derived from classical elastoplastic mechanics. These models reveal the creep behavior of specific rocks, but they often do not sufficiently take the influence of creep parameters into account, or they overlook the certain factorsduring researchs, so establishing a model that comprehensively incorporates parameters and accurately reveals the creep mechanical behavior of rocks is of significant value. Methods The Nishihara model in the classical elastoplastic framework was selected, and damage was introduced to compensate for its inability to reveal the nonlinear creep characteristics of rocks. The Weibull distribution function was applied as the damage statistical function, and the SMP criterion was adopted as the strength criterion of sandstone to derive the damage variables under three-dimensional conditions. A conventional triaxial compression test was conducted using the MTS-815 rock rigidity testing machine, and the applicability of the damage statistical model was verified by comparing the compression test results with the damage evolution curve. The yield function was applied to fit the yield surface of sandstone to determine the yield stress under each confining pressure, and a fracture plastic body was incorporated to represent the plastic deformation of the primary fracture. The components of each part of the model were then integrated to derive the constitutive equation in accordance with the mechanical boundary. The STAC600-600 rock rheology testing machine was utilized to perform the rock classification rheological test, and the fitting results of the creep model were analyzed. The model parameters were optimized to improve the fitting accuracy. Finally, the fitting results of this model were compared to those of other models to analyze the advantages and disadvantages of the approach. Results and Discussions The damage variables derived from the SMP criterion and Weibull distribution function showed strong consistency with the rock deformation behavior under conventional triaxial compression. Damage did not develop during the early loading stage, and when the deviator stress exceeded the yield stress, damage developed rapidly until destruction ofwhole rock was happened. The damage evolution trend indicated that the fracture unit of the rock accumulated rapidly once the deviator stress reached the yield limit, and the rate of damage accumulation varied under different confining pressures, with lower confining pressure leading to a faster accumulation rate. The improved Nishihara model was applied to fit the creep deformation curve and produced satisfactory results; however, the fitting results of the accelerated creep curve under 8, 12, and 20 MPa confining pressures all shifted forward. The Poisson ratio μ was refitted as a function related to confining pressure, and the accuracy improved significantly. An algebraic relationship existed between the key model parameters u and m and the confining pressure, and the constitutive parameters under each confining pressure were obtained after establishing this relationship. Compared to the fitting results of the other two models, the fitting accuracy of this model was higher, and the variable “time”, which often caused errors, did not participate in the parameter calculation. The damage parameters were directly related to stress and strain, which improved accuracy to a certain extent. Conclusions The results indicate that the damage statistical creep constitutive model accurately fits the creep deformation characteristics of rocks under three-dimensional conditions, expands the range of influencing factors for constitutive parameters, and provides practical guidance for related engineering applications.

Graphical abstract

关键词

岩石蠕变 / SMP准则 / 损伤变量 / Weibull分布 / 西原模型

Key words

rock creep / SMP criterion / damage variables / Weibull distribution / Nishihara model

本刊网刊

岩石蠕变是应力和温度不变的情况下，变形随时间增长的现象^[1]。深部工程中围岩的应力场和应变分布表现出显著的时间相关性。研究岩石蠕变特性具有重要意义，本构关系能直观反映岩石变形破坏和力学响应，建立适合的蠕变本构模型具有重要价值。

已有的本构模型强调对岩石加速蠕变阶段的描述，主要分为3种。1）采用非线性元件代替线性元件进行描述，如康永刚等^[2]采用非定常黏壶替换模型中的定常黏壶，提出了一种非定常Burgers模型；王春萍等^[3]采用分数阶黏弹性元件替换Burgers模型中的黏性元件，并在黏弹性元件上叠加损伤，建立了裂隙花岗岩的非线性蠕变本构模型。2）组合非线性元件与线性元件来描述，如刘东燕等^[4]用串联的手段组合Kelvin模型与Burgers模型，构建一种新的黏弹塑性模型；齐亚静等^[5]通过在西原模型中串联一个具备应变触发特性的非线性黏壶，提出一种改进版的西原模型，并将其成功拓展应用于三维蠕变模型中；康永刚等^[6]提出一种改进的岩石蠕变模型，采用非牛顿黏壶替代传统Kelvin模型中的牛顿黏壶，更精确地刻画岩石的衰减蠕变和稳态蠕变特性。3）在元件中融入损伤变量使模型参数非线性化，如于冰冰等^[7]将以时间相关的损伤因子引入弹性和黏性元件中，并采用裂隙塑性体改进Bingham模型；宋战平等^[8]采用时效非线性损伤体改进Burgers模型，并通过反向传播（BP）神经网络对参数辨识过程进行优化；刘文博等^[9]对加速蠕变阶段的时间微分处理，推导了受时间、应力双重影响的蠕变参数时效性劣化型蠕变本构模型。

在上述模型中，损伤型模型最契合岩石力学特征，前人对损伤型本构模型的研究已取得了诸多成果，包括在水化学损伤^[10]、干湿循环损伤^[11]、高温损伤^[12‒13]等条件下开展的一系列研究，并建立了对应的本构模型。但较多模型考虑参数时不够全面，一些经验模型虽借助深度分析^[14]完成了结果的优化，但未消除模型本身由于忽略某些参数带来的影响。

为建立考虑岩石原生裂隙和围压影响的三维蠕变本构模型，本文依托砂岩常规三轴与分级蠕变试验结果分析了岩石损伤随压缩进程的关系，借助SMP准则推导了Weibull分布形式的损伤并叠加至西原模型元件，最后串联一个裂隙塑性体用于表征原生裂隙影响，构建了能描述岩石蠕变全过程的损伤本构模型，将泊松比影响纳入模型进行优化，并对比分析了其他模型。

1 概率密度函数与强度准则

1.1 Weibull分布函数

岩石强度失效是组成岩石的介质产生损伤并不断累积的结果，概率分布能表征岩石材料非均质强度特性。概率密度函数形式多样，各有优劣。正态分布函数^[15]能反映岩石变形全过程但模型参数需查询正态分布表；幂函数^[16]无法反映岩石的峰后变形特性；Weibull分布函数能反映后期变形特征却无法反映岩石加载初期的变形特征^[17‒19]，而岩石损伤也不存在于加载初期，因此它能够契合岩石的损伤演化规律。

岩石变形破坏过程中，内部的裂隙随机分布，根据损伤发育的随机性，采用Weibull分布函数^[20]描述岩石的强度特征：

P(F)=\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{F}{u}\right)}^m\right)（1）

p(F)=\frac{m}{u}{\left(\frac{F}{u}\right)}^{m-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{F}{u}\right)}^m\right)（2）

式中，F为岩石的强度，u为基元体力学属性，m为均匀性系数，p(F)为统计分布密度，P(F)为统计分布函数。

若定义损伤变量D为损伤微元数n占全岩石微元数N的比例，那么：

D=\frac{n}{N}（3）

岩石所受偏应力\sigma在区间(0，\sigma)上时，若破坏的微元数服从Weibull概率分布，根据积分，损伤微元数n可表示为：

n={\int }_{\mathrm{0}}^{\sigma }Np(\sigma )\mathrm{d}\sigma =NP(\sigma )（4）

联立式（1）、（3）、（4）可将损伤表示为：

D=\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{F}{u}\right)}^m\right)（5）

1.2 SMP准则

随着对岩石力学的深度探索，岩石研究中产生了许多新的强度准则，传统Mohr‒Coulomb准则^[21‒22]是以最大和最小主应力来构建的二元摩擦准则，Hork‒Brown准则^[23]可以反映出岩石与岩体的固有非线性破裂模式，并且适合解释各向异性的岩体，但这种方式没有考虑到中间主应力的作用。

SMP准则即Shao、Matsushima和Ping准则^[24‒25]，它考虑了中间主应力的影响，能更准确地描述岩石的三维力学特性。其表达式为：

\frac{I_{\mathrm{1}}{I}_{\mathrm{2}}}{I_{\mathrm{3}}}=\mathrm{8}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\phi +\mathrm{9}（6）

式中：I_{\mathrm{1}}、I_{\mathrm{2}}、I_{\mathrm{3}}为应力第1、2、3不变量，可通过最大、最小和中间主应力求取；\phi为内摩擦角。

Satake等^[26]拓展了SMP准则，提出新的黏结应力定义即{\sigma }_{\mathrm{0}}=Ccot φ，其中，{\sigma }_{\mathrm{0}}为黏结应力，C为黏聚力。

SMP准则强度公式为：

F=\frac{{\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{0}}}{{\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}}}（7）

式中，{\sigma }_{\mathrm{1}}、{\sigma }_{\mathrm{3}}分别为最大、最小主应力。

2 三轴条件下的岩石损伤统计

2.1 损伤参数求取

岩石内部未产生损伤时，可视为弹性体，根据广义HOOK定律弹性体本构关系三维形式，可将应变{\epsilon }_{ij}^{}表示为:

\begin{array}{c}{\epsilon }_{ij}^{}=\frac{\mathrm{1}}{\mathrm{2}G}{s}_{ij}+\frac{\mathrm{1}}{\mathrm{3}K}{\sigma }_{\mathrm{m}}{\delta }_{ij}=\frac{\mathrm{1}}{E}\left(\left(\mathrm{1}+\mu \right)s_{ij}+\left(\mathrm{1}-\mathrm{2}\mu \right){\sigma }_{\mathrm{m}}{\delta }_{ij}\right)\end{array}（8）

式中：G为剪切弹性模量；K为体积模量；s_{ij}为应力偏张量；{\sigma }_{\mathrm{m}}为平均应力；{\delta }_{ij}为Kronecker记号；E为弹性模量；\mu为泊松比；i、j为弹塑性力学指标符号，i、j=1,2,3。

常规三轴压缩条件下有以下关系：

{\sigma }_{\mathrm{2}}={\sigma }_{\mathrm{3}}（9）

{\sigma }_{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}({\sigma }_{\mathrm{1}}+\mathrm{2}{\sigma }_{\mathrm{3}})（10）

s_{ij}={\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{m}}=\frac{\mathrm{2}}{\mathrm{3}}({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}})（11）

式（9）中，{\sigma }_{\mathrm{2}}为中间主应力。

联立式（8）、（10）、（11）可得：

E{\epsilon }_{\mathrm{1}}^{}={\sigma }_{\mathrm{1}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}（12）

式中，{\epsilon }_{\mathrm{1}}^{}为轴向应变。

联立式（7）、（12），强度F可表示为：

F_{}=\frac{E{\epsilon }_{\mathrm{1}}^{}+\mathrm{2}\mu {\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}}}{{\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}}}（13）

根据Lematrie应变等价假说^[27]，假设岩石产生的损伤满足式（3）、（4）所作假设，则有效主应力{\sigma }_i^{\text{'}}与主应力{\sigma }_i满足以下关系：

{\sigma }_i^{\text{'}}=\frac{{\sigma }_i}{\mathrm{1}-D}（14）

将式（12）中主应力{\sigma }_{\mathrm{1}}、{\sigma }_{\mathrm{3}}替换为有效主应力并联立式（14）可获得损伤条件下的应力‒应变本构关系：

\frac{{\sigma }_{\mathrm{1}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}}{\mathrm{1}-D}=E{\epsilon }_{\mathrm{1}}^{}（15）

联立式（5）、（15）可得：

{\epsilon }_{\mathrm{1}}^{}E\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{F}{u}\right)}^m\right)={\sigma }_{\mathrm{1}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}（16）

将应力‒应变曲线峰值点({\sigma }_{\mathrm{1}\mathrm{c}},{\epsilon }_{\mathrm{1}\mathrm{c}})代入式（16），此时轴向应力{\sigma }_{\mathrm{1}}对轴向应变{\epsilon }_{\mathrm{1}}的偏导数为0，可表示为：

\underset{\begin{array}{c}{\epsilon }_{\mathrm{1}}^{}\to {\epsilon }_{\mathrm{1}\mathrm{c}}^{}\\ {\sigma }_{\mathrm{1}}\to {\sigma }_{\mathrm{1}\mathrm{c}}\end{array}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\partial ({\sigma }_{\mathrm{1}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}})}{\partial {\epsilon }_{\mathrm{1}}^{}}=\mathrm{0}（17）

将式（13）代入式（16），并在峰值点({\sigma }_{\mathrm{1}\mathrm{c}},{\epsilon }_{\mathrm{1}\mathrm{c}})处求偏导数，即：

\frac{\partial ({\sigma }_{\mathrm{1}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}})}{\partial {\epsilon }_{\mathrm{1}}^{}}=

E\left(\mathrm{1}-\frac{mE{\epsilon }_{\mathrm{1}\mathrm{c}}}{({\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}})u}{\left(\frac{F_{\mathrm{c}}^{}}{u}\right)}^{m-\mathrm{1}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{F_{\mathrm{c}}^{}}{u}\right)}^m\right)（18）

式中，F_{\mathrm{c}}为式（13）中应变取{\epsilon }_{\mathrm{1}\mathrm{c}}对应的强度值。

联立式（16）～（18）可解出u和m：

{\left(\frac{F_{\mathrm{c}}^{}}{u}\right)}^m=\frac{({\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}})F_{\mathrm{c}}^{}}{mE{\epsilon }_{\mathrm{1}\mathrm{c}}^{}}=\mathrm{l}\mathrm{n}\frac{E{\epsilon }_{\mathrm{1}\mathrm{c}}}{{\sigma }_{\mathrm{1}\mathrm{c}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}}（19）

m=\frac{({\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{0}})F_{\mathrm{c}}^{}}{E{\epsilon }_{\mathrm{1}\mathrm{c}}^{}\mathrm{l}\mathrm{n}\frac{E{\epsilon }_{\mathrm{1}\mathrm{c}}}{{\sigma }_{\mathrm{1}\mathrm{c}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}}}（20）

u={\frac{F_{\mathrm{c}}^{}}{\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{l}\mathrm{n}\left(\mathrm{l}\mathrm{n}\frac{E{\epsilon }_{\mathrm{1}\mathrm{c}}}{{\sigma }_{\mathrm{1}\mathrm{c}}-\mathrm{2}\mu {\sigma }_{\mathrm{3}}}\right)}{m}\right)}}^{}（21）

2.2 损伤演化规律验证

用MTS‒815岩石刚性试验机对砂岩试样进行常规三轴压缩试验，砂岩的基础物理力学参数见表1。

将试验所得应力‒应变关系曲线与损伤演化曲线（图1）进行对比，发现损伤发育与否与施加轴向偏应力和屈服应力的关系密切，施加轴向偏应力大于屈服应力时，损伤才发展，并在接近峰值应力时加剧，最终导致岩石破坏。损伤统计模型能有效描述三轴压缩条件下砂岩的累进性破坏规律，其模型参数见表2。

3 考虑原生裂隙的岩石塑性变形

大多数岩石在形成、演化过程中及改造作用下会形成一定的天然裂隙。图2反映了裂隙的受压变形示意图。图2中，对岩石施加荷载时，其裂隙表面并不平整且伴随剪切位移产生^[28]，即便应力释放后，裂隙也无法完全闭合。

文献[29]揭示了非充填粗糙结构面的法向应力{\sigma }_{\mathrm{n}}与法向变形{\delta }_{\mathrm{n}}的关联，描述为：

{\sigma }_{\mathrm{n}}=\frac{{\delta }_{\mathrm{n}}}{a-b{\delta }_{\mathrm{n}}}（22）

式中，a和b分别代表结构面几何特征和岩石力学属性。岩石未遭受完全破坏时，其内部所有裂隙具有非贯通性。如想保持一定的张开状态，需要克服张开回弹力{\sigma }_{\mathrm{e}}，当法向应力{\sigma }_{\mathrm{n}}达到闭合门槛值{\sigma }_{\mathrm{D}}时，有以下关系：

{\sigma }_{\mathrm{n}}+{\sigma }_{\mathrm{e}}=\frac{{\delta }_{\mathrm{n}}}{a-b{\delta }_{\mathrm{n}}},{\sigma }_{\mathrm{n}}⩾{\sigma }_{\mathrm{D}}（23）

对式（23）进行变形处理，裂隙闭合产生的法向塑性变形表示为：

{\delta }_{\mathrm{n}}=\frac{a{\sigma }_{\mathrm{e}}+a{\sigma }_{\mathrm{n}}}{\mathrm{1}+b{\sigma }_{\mathrm{e}}+b{\sigma }_{\mathrm{n}}}（24）

对于裂隙岩石，假设原生裂隙对应的损伤为D₀，根据Lematrie应变等价假说，作用于岩石的应力可用有效应力{\sigma }^{\text{'}}表示：

{\sigma }^{\text{'}}=\frac{\sigma }{\mathrm{1}-D_{\mathrm{0}}}（25）

根据文献[30]，裂隙面的法向应力和裂隙岩石沿载荷方向的塑性应变分别表示为：

{\sigma }_{\mathrm{n}}={\sigma }^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha（26）

{\epsilon }_{\mathrm{p}}=\frac{{\delta }_{\mathrm{n}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{V_{\mathrm{R}}b}（27）

式中，\alpha为裂隙面法向与受载方向的夹角，V_R为裂隙体积，{\epsilon }_{\mathrm{p}}^{}为原生裂隙闭合产生的塑性变形。

联立式（25）～（27），{\epsilon }_{\mathrm{p}}^{}可表示为：

{\epsilon }_{\mathrm{p}}^{}=\frac{a\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{V_{\mathrm{R}}b}\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}+b{\sigma }_{\mathrm{e}}+{\sigma }^{\text{'}}b\mathrm{c}\mathrm{o}\mathrm{s}\alpha }\right)（28）

令Q=\frac{a\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{V_{\mathrm{R}}b},W=b{\sigma }_{\mathrm{e}},R=\frac{b\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{\mathrm{1}-D_{\mathrm{0}}}，则式（28）简化为：

（29）

4 岩石蠕变损伤本构模型

4.1 蠕变模型结构与蠕变方程

西原模型^[31]综合考虑了岩石的弹性、黏弹性和黏塑性，能有效地刻画岩石蠕变过程的前两个阶段，但在描述岩石的非线性变形特性时存在局限。可利用损伤对西原模型进行改进并串联一个裂隙塑性体，使模型能够描述岩石的整个蠕变过程。

西原模型由虎克体、开尔文体和理想黏塑性体组成，如图3所示。图3中，E_{\mathrm{0}}、E_{\mathrm{1}}分别为弹性、黏弹性模量，{\eta }_{\mathrm{1}}、{\eta }_{\mathrm{2}}为黏滞系数，{\sigma }_{\mathrm{Y}}为屈服应力。西原模型一维条件下的本构方程为式（30），式中t为加载时间。

（30）

三维条件下，采用由E_{\mathrm{0}}和\mu确定的弹性剪切模量G_{\mathrm{0}}和体积模量K替代弹性模量E_{\mathrm{0}}，采用黏弹性剪切模量G_{\mathrm{1}}替代黏弹性模量E_{\mathrm{1}}，替换后的西原模型本构方程为：

（31）

采用损伤对西原模型进行改进并增加一个裂隙塑性体，本文最终的本构方程可表示为：

（32）

4.2 模型特性

本构模型结构如图4所示。该结构由3部分组成，第1、2部分为损伤改进的西原模型，第3部分为带应力开关的裂隙塑性体。

当，本构模型简化至广义开尔文模型，该模型会展现出瞬时弹性变形与稳定蠕变的特性。其中，瞬时应变与偏应力之间呈现出线性相关关系。

当{\sigma }_{\mathrm{D}}\le \sigma \le {\sigma }_{\mathrm{Y}}，损伤变量D失效，该情形下，裂隙塑性体会参与到蠕变过程中。此时的模型会出现瞬时塑性变形，并呈现稳定的蠕变状态。同时，瞬时应变和偏应力之间存在非线性关系。

当\sigma >{\sigma }_{\mathrm{Y}}，损伤变量D和黏塑性体生效，由于损伤变量为强度F的函数，F由应力\sigma和应变\epsilon决定，黏塑性变形随时间增加，损伤增加并作用于各元件，模拟加速蠕变。

4.3 模型参数确定

模型参数包括K、G_{\mathrm{0}}、G_{\mathrm{1}}、{\eta }_{\mathrm{1}}、{\eta }_{\mathrm{2}}、W、R、Q，根据蠕变曲线特点判断施加应力所处区间后，通过以下方式分别求取：

1）当瞬时应变与偏应力呈线性关系时，，联立两组蠕变曲线可求解K、G_{\mathrm{0}}。随后岩石发生黏弹性变形，通过函数拟合可获得G_{\mathrm{1}}、{\eta }_{\mathrm{1}}。

2）当瞬时应变与偏应力呈非线性关系且不发生加速蠕变时，{\sigma }_{\mathrm{D}}\le \sigma \le {\sigma }_{\mathrm{Y}}。蠕变曲线的起始点，有瞬时应变{\epsilon }_{\mathrm{i}\mathrm{1}}^{}，{\epsilon }_{\mathrm{i}\mathrm{2}}^{}，{\epsilon }_{\mathrm{i}\mathrm{3}}^{}，…，{\epsilon }_{\mathrm{i}l}^{}，l为恒载等级，这些瞬时应变可被视为独立的瞬时弹性应变和塑性应变，减去弹性应变后，存在塑性应变{\epsilon }_{\mathrm{p}\mathrm{2}}^{}，{\epsilon }_{\mathrm{p}\mathrm{2}}^{}，{\epsilon }_{\mathrm{p}\mathrm{3}}^{}，…，{\epsilon }_{\mathrm{p}l}^{}，参数W、R、Q可以通过式（29）采用最小二乘法估计。

3）当岩石发生蠕变破坏时，\sigma >{\sigma }_{\mathrm{Y}}。取蠕变曲线直线段，其斜率倒数即为{\eta }_{\mathrm{2}}。蠕变破坏仅能在试验后获得，本文选取能表征弹塑性模型塑性屈服面的屈服准则（式（33））判定加速蠕变发生的应力区间。

f({\sigma }_{ij},{\alpha }_{\mathrm{p}})=q-{\alpha }_{\mathrm{p}}({\gamma }_{\mathrm{p}})R_{\mathrm{c}}\sqrt[]{A\left(C_{\mathrm{s}}+\frac{p}{R_{\mathrm{c}}}\right)}=\mathrm{0}（33）

式中：p为平均应力；q为偏应力，即σ₁‒σ₃；α_p(γ_p)为塑性硬化函数，值域为[1，2)；γ_p为塑性应变；R_{\mathrm{c}}为单元的强度；C_{\mathrm{s}}为内聚力系数；A为内摩擦有关系数；{\sigma }_{ij}为主应力。

5 模型验证与本构参数优化

5.1 蠕变试验设计

为了确认模型的有效性，采用STAC600‒600型号的岩石三轴流变试验机（图5）进行分级蠕变试验。该设备通过设置5阶段的轴向压力值来模拟相当于不同围压下峰值强度50%~90%的偏应力条件，即设置5级恒定轴向荷载，每级荷载持续时间48 h以上，试验前后的试件状态如图6所示。

5.2 蠕变试验结果与特征分析

图7（a）和（b）分别为0 和8 MPa围压下的分级加载结果。围压和偏应力提升时，砂岩的累积蠕变变形显著增加；相较于0围压，8 MPa围压下岩石最大破坏上升58.94%；围压升高时，岩石破坏特征从脆性向延性转变，且破裂面与荷载方向夹角增大，围压为0时夹角为13°，围压为8 MPa时夹角为36°。

5.3 模型应用与验证

根据三轴压缩试验数据，借助屈服准则表达式（34）拟合获得屈服面（图8），拟合参数见表3，通过屈服函数可获取不同围压下的屈服应力，结合第4.3节判断各级蠕变曲线所在应力区间。

分别采用加入裂隙塑性体前后的改进西原模型A、B对试验曲线进行拟合，结果见表4、5。表4中相同围压下体积模量K和剪切弹性模量G₀存在一定波动，该现象可能与岩石内部的天然裂隙有一定联系，试验所用砂岩为微风化砂岩，原生裂隙影响较小，对比不强烈，对风化程度为强风化甚至全风化的岩石，裂隙塑性体作用会更显著。

采用最终确定的蠕变损伤本构模型（对应图5的模型结构和式（32）的模型方程）对不同围压下发生蠕变破坏的蠕变变形进行拟合，结果见图9。围压为0时拟合效果较好，但在围压为8、12、20 MPa条件下破坏阶段拟合曲线有明显前移迹象，可能是损伤建立过程忽略了不同围压下泊松比的改变，进而影响损伤的计算结果并导致拟合偏差。因此，建立损伤统计模型时需要将泊松比也视作变量，进一步优化模型参数。

5.4 模型参数优化

为精确模拟各种围压条件下的试验结果，需要获取砂岩泊松比与围压之间的关系曲线，各围压下泊松比的研究确定方法参照文献[32]。不同围压下砂岩泊松比如图10所示。图10中，x为当前围压，随着围压的提升，砂岩泊松比呈现阶梯式下降。

确定单轴压缩条件下的泊松比函数后，围压与泊松比的关系可表示为：

{\mu }_x^{\mathrm{*}}=\frac{{\mu }_{\mathrm{0}}^{*}}{h+\frac{g{\sigma }_{\mathrm{3}}}{\mathrm{4}}}（34）

式中：{\mu }_x^{*}为某围压下三轴压缩条件的泊松比函数；{\mu }_{\mathrm{0}}^{*}为单轴压缩条件下的泊松比函数；g、h为拟合参数，本文取值分别为1和0.06。考虑泊松比影响之后的修正模型拟合结果如图11所示，该模型对于三轴条件下的砂岩加速蠕变曲线拟合效果有明显改善，8、12和20 MPa围压下拟合度R²分别较修正前提高了0.06、0.02和0.04。

优化的实质是采用函数{\mu }_x^{*}替代式（13）F表达式中的常量\mu，凡是与围压相关的参数（如屈服应力、泊松比、弹性模量等）均被作为损伤变量D的内变量或影响因素，通过影响损伤变量进而影响模型元件达到拟合非线性变形的目的，该过程是通过SMP准则和Weibull分布函数共同实现的。

根据图12损伤参数拟合曲线，三轴条件下计算损伤的重要参数u、m与围压σ₃分别呈二次函数和直线关系，u=7.44-0.09σ₃+σ₃²，m=24.18-0.70σ₃，u和m均能通过拟合关系式获得。因此，获得了岩石的基本物理参数后，该模型对蠕变过程的拟合较为简便。

5.5 模型应用对比

本文所建立的蠕变损伤本构模型适用性获得了验证。为进一步证明该本构模型相较于其他本构模型的优缺点，选取另外两种结构的蠕变本构模型在已有试验的基础上进行对比分析。

1）带非线性黏壶的改进Burgers模型，即在Burgers模型上串联一个非线性黏塑性体，该本构模型结构如图13所示，本构方程为式（35）。

(35)

式中，v为拟合参数。

在参数求取与拟合过程中，需要将时间分段并获取加速流变的起始时间t₂，这需要首先分析流变速率s与偏应力q的代数关系，拟合并找到流变速率曲线的转折点，该转折点的确定有很大的主观性。因此，该模型对加速蠕变阶段拟合精度较低，但处理过程较为简便。

2）时效损伤改进的经典弹塑性本构模型，该类模型将损伤直接与时间联系，属经验型模型，同样对时间依赖性极强。以时效损伤的西原模型为例，其示意图如图14所示，本构方程为式（36）、（37）。

(36)

D=\mathrm{1}-w{\sigma }^{zt}（37）

式中，w和z为拟合参数。

该模型的缺点与带非线性黏壶的改进Burgers模型相同，但运用随机森林、神经网络等算法能将拟合精度提升至一定水平。

两种模型的拟合结果如图15所示，在不采用其他手段修正的情况下，拟合精度相较于本文所建立的本构模型偏低。总的来说，本文建立的本构模型优点在于将时间变量剔出了模型参数的求取过程，减少误差的产生来源，并且将一维参数都拓展为三维形式，获得基本参数后具备完整的拟合系统，除特定的应用场景外，该模型考虑的影响因素较全面。缺点在于损伤变量求取及本构模型的建立过程比较复杂，前期参数计算比较烦琐。

6 结 论

本研究综合运用室内实验与理论分析的方法，得出以下主要结论：

1）将SMP准则作为强度准则，结合Weibull分布函数构建损伤变量能够有效描述砂岩常规三轴压缩过程与蠕变过程的损伤演化过程，表明了该准则下该形式损伤变量的适用性。

2）损伤变量中的重要参数u、m与围压有着明显的函数关系，由于本文损伤变量是基于强度F建立的，F又为主应力与黏结应力的表达式，因此获取任意岩石的u、m与围压关系曲线后，能直接采用应力表达损伤。

3）采用加入裂隙塑性体的改进西原模型对砂岩蠕变过程进行拟合时，效果较好，但围压不为0时，加速变形曲线前移，考虑泊松比的影响能有效解决这一问题。

4）本构模型中的相关参数（如泊松比、屈服应力等）在三维应力状态下得到了扩展，同时裂隙塑性体的引入使模型不仅适用于完整岩石，还能够描述带有裂隙或结构面的岩石。因此，除特殊环境，该模型适用于一些常见的工程条件，且对其他岩石类型也具有重要的参考价值，但某些物理量，如泊松比随着压缩进程也会发生改变，应当继续本研究并建立更加科学严谨的模型。

5）对比3种形式的本构模型拟合结果发现，本文所建立的泊松比函数修正的蠕变损伤本构模型拟合精度优于另外两种本构模型，且损伤与时间独立，规避了人为处理时间因素的误差，但前期处理较烦琐。

参考文献

[1] Sun Jun.Rock rheological mechanics and its advance in engineering applications[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(6):1081‒1106.

[2] 孙钧.岩石流变力学及其工程应用研究的若干进展[J].岩石力学与工程学报,2007,26(6):1081‒1106.

[3] Kang Yonggang, Zhang Xiue.A non-stationary model for rock creep based on Burgers model[J].Rock and Soil Mechanics,2011,32(S1):424‒427.

[4] 康永刚,张秀娥.基于Burgers模型的岩石非定常蠕变模型[J].岩土力学,2011,32(S1):424‒427.

[5] Wang Chunping, Liu Jianfeng, Chen Liang,et al.Nonlinear creep constitutive model of fractured granite and its parameters analysis[J].Advanced Engineering Sciences,2024,56(4):141‒149.

[6] 王春萍,刘建锋,陈亮,等.裂隙花岗岩非线性蠕变本构模型及参数分析[J].工程科学与技术,2024,56(4):141‒149.

[7] Liu Dongyan, Xie Linjie, Xiaofeng Tuo,et al.Study on creep characteristics and nonlinear viscoelastic-plastic model of sandstone under different confining pressures[J].Chinese Journal of Rock Mechanics and Engineering,2017,36(S2):3705‒3712.

[8] 刘东燕,谢林杰,庹晓峰,等.不同围压作用下砂岩蠕变特性及非线性黏弹塑性模型研究[J].岩石力学与工程学报,2017,36(S2):3705‒3712.

[9] 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.

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

[11] Kang Yonggang, Zhang Xiue.Nonstationary parameter fractional Burgers model of rock creep[J].Rock and Soil Mechanics,2011,32(11):3237‒3241.

[12] 康永刚,张秀娥.岩石蠕变的非定常分数伯格斯模型[J].岩土力学,2011,32(11):3237‒3241.

[13] Yu Bingbing, Li Qing, Zhao Tongde,et al.Full-time nonlinear creep damage model of fractured rock mass based on stress-time double threshold[J].Chinese Journal of Rock Mechanics and Engineering,2023,42(8):1928‒1944.

[14] 于冰冰,李清,赵桐德,等.基于应力与时间双阈值的岩石全时态非线性蠕变损伤模型[J].岩石力学与工程学报,2023,42(8):1928‒1944.

[15] Song Zhanping, Wang Bowen, Fan Shengyuan,et al.Nonlinear creep intrinsic modeling of hard rock considering aging damage[J].Chinese Journal of Rock Mechanics and Engineering,2024,43(10):2368‒2380.

[16] 宋战平,王博文,范胜元,等.考虑时效损伤的硬岩非线性蠕变本构模型研究[J].岩石力学与工程学报,2024,43(10):2368‒2380.

[17] Liu Wenbo, Zhang Shuguang, Sun Boyi.Time-dependentness of creep parameter degradation and variable parameter creep model[J].Advanced Engineering Sciences,2021,53(1):104‒112.

[18] 刘文博,张树光,孙博一.蠕变参数劣化的时效性及变参数蠕变模型[J].工程科学与技术,2021,53(1):104‒112.

[19] Chen Youliang, Chen Qijian, Xiao Peng,et al.A true triaxial creep constitutive model for rock considering hydrochemical damage[J].Chinese Journal of Theoretical and Applied Mechanics,2023,55(1):159‒168.

[20] 陈有亮,陈奇键,肖鹏,等.考虑水化学损伤的岩石真三轴蠕变本构模型[J].力学学报,2023,55(1):159‒168.

[21] Hou Di, Xu Zhen, Mo Zhikun,et al.Research and application of rock creep constitutive model under different dry-wet cycles[J].Yangtze River,2023,54(12):218‒224.

[22] 侯迪,许珍,莫至坤,等.干湿循环条件下岩石蠕变本构模型研究及应用[J].人民长江,2023,54(12):218‒224.

[23] Meng Lubo, Zhang Shan, Li Tianbin,et al.Study on the creep constitutive model of layered rockconsidering anisotropic and damage factors after hightemperature exposure[J].International Journal of Damage Mechanics,2025,34(3):520‒541. doi:10.1177/10567895241292748

[24] Shi Sheng, Zhu Fengjin, Zhu Jiancai.Investigation on damage creep constitutive model of rock under the coupled effect of freeze‒thaw cycles and loading[J].Fatigue & Fracture of Engineering Materials & Structures,2024,47(4):1281‒1299. doi:10.1111/ffe.14244

[25] Wang Junbao, Song Zhanping, Zhao Baoyun,et al.A study on the mechanical behavior and statistical damage constitutive model of sandstone[J].Arabian Journal for Science and Engineering,2018,43(10):5179‒5192. doi:10.1007/s13369-017-3016-y

[26] Ma Qiufeng, Qin Yueping, Zhou Tianbai,et al.Constitutive model of rock compaction stage based on contact theory[J].Journal of Central South University (Science and Technology),2019,50(8):1941‒1948.

[27] 马秋峰,秦跃平,周天白,等.基于颗粒接触理论的岩石压密阶段本构模型[J].中南大学学报(自然科学版),2019,50(8):1941‒1948.

[28] Ji Ming, Sun Zhongguang, Liu Wenpeng,et al.Study on sandstone damage constitutive model based on power function distribution[J].Journal of Xi’an University of Architecture & Technology (Natural Science Edition),2023,55(3):324‒331.

[29] 季明,孙中光,刘文朋,等.基于幂函数分布的砂岩损伤本构模型研究[J].西安建筑科技大学学报(自然科学版),2023,55(3):324‒331.

[30] Zheng Dong.Research on statistical damage constitutive model of rocks under high-temperature and high confining pressure conditions based on Weibull distribution of element strength[D].Beijing:Beijing University of Civil Engineering,2024.

[31] 郑栋.基于单元强度Weibull分布的高温高围压岩石统计损伤本构模型研究[D].北京:北京建筑大学,2024.

[32] Bai Yu, Zhao Kangpu, Liu Yi,et al.Shale dynamic damage constitutive model based on Weibull distribution under variable confining pressure[J].Engineering Blasting,2023,29(6):61‒67.

[33] 白宇,赵康朴,刘依,等.基于Weibull分布页岩围压下的动态损伤本构模型研究[J].工程爆破,2023,29(6):61‒67.

[34] Jiang Haopeng, Jiang Annan, Yang Xiurong.Statistical damage constitutive model of high temperature rock based on Weibull distribution and its verification[J].Rock and Soil Mechanics,2021,42(7):1894‒1902.

[35] 蒋浩鹏,姜谙男,杨秀荣.基于Weibull分布的高温岩石统计损伤本构模型及其验证[J].岩土力学,2021,42(7):1894‒1902.

[36] Zhang Zongtang, Gao Wenhua, Zhang Zhimin,et al.Evolution of particle disintegration of red sandstone using Weibull distribution[J].Rock and Soil Mechanics,2020,41(3):877‒885.

[37] 张宗堂,高文华,张志敏,等.基于Weibull分布的红砂岩颗粒崩解破碎演化规律[J].岩土力学,2020,41(3):877‒885.

[38] Gu Qiuxin, Zhang Qiang, Dai Wanli,et al.A novel statistical damage constitutive model of rock joints considering normal stress and joint roughness[J].International Journal of Damage Mechanics,2025,34(2):326‒351. doi:10.1177/10567895241277681

[39] Zhang Chao, Yu Jin, Bai Yun,et al.Statistical damage constitutive model of rock brittle-ductile transition based on strength theory[J].Chinese Journal of Rock Mechanics and Engineering,2023,42(2):307‒316.

[40] 张超,俞缙,白允,等.基于强度理论的岩石脆延转化统计损伤本构模型[J].岩石力学与工程学报,2023,42(2):307‒316.

[41] Wei Chao, Zhao Cheng, Zhao Chunfeng,et al.Statistical damage constitutive model for rocks considering temperature effects and its validation[J].Journal of Central South University (Science and Technology),2024,55(3):1056‒1067.

[42] 魏超,赵程,赵春风,等.考虑温度作用的岩石统计损伤本构模型及验证[J].中南大学学报(自然科学版),2024,55(3):1056‒1067.

[43] Shi Lei, Hou Kepeng, Sun Huafen.Damage softening model of freeze-thaw loaded rock based on SMP criterion[J].Science Technology and Engineering,2023,23(4):1658‒1664.

[44] 史磊,侯克鹏,孙华芬.基于SMP准则的冻融受荷岩石损伤软化模型[J].科学技术与工程,2023,23(4):1658‒1664.

[45] Zhu Jianming, Cheng Haifeng, Yao Yangping.Statistical damage softening model of fractured rock based on SMP failure criterion and its application[J]. Chinese Journal of Rock Mechanics and Engineering, 2013,32(S2):3160‒3168.

[46] 朱建明,程海峰,姚仰平.基于SMP准则的破裂岩统计损伤软化模型及其应用研究[J].岩石力学与工程学报,2013,32(S2):3160‒3168.

[47] Satake M.Stress-deformation and strength characteristics of soil under three difference principle stresses(Discussion)[J]. Proc. of Japan Society of Civil Engineers,1976,246:137‒138. doi:10.2208/jscej1969.1976.246_137

[48] Chen Song, Cao Xiuling, Yang Zhao.Three dimensional statistical damage constitutive model of rock based on Griffith strength criterion[J].Geotechnical and Geological Engineering,2021,39(8):5549‒5557. doi:10.1007/s10706-021-01844-6

[49] Wang Lei, Liang Weiguo.Experimental study on cracking initiation and fracture propagation of supercritical CO₂/clean water fracturing coal[J]. Chinese Journal of Rock Mechanics and Engineering, 2019,38(S1):2680‒2689

[50] 王磊,梁卫国.超临界CO₂/清水压裂煤体起裂和裂缝扩展试验研究[J].岩石力学与工程学报,2019,38(S1):2680‒2689.

[51] Barton N, Bandis S, Bakhtar K.Strength,deformation and conductivity coupling of rock joints[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,1985,22(3):121‒140. doi:10.1016/0148-9062(85)93227-9

[52] Wang Qihu, Ye Yicheng, Liu Yanzhang,et al.A creep constitutive model of rock considering initial damage and creep damage[J].Rock and Soil Mechanics,2016,37(S1):57‒62.

[53] 王其虎,叶义成,刘艳章,等.考虑初始损伤和蠕变损伤的岩石蠕变全过程本构模型[J].岩土力学,2016,37(S1):57‒62.

[54] Jia Baoxin, Chen Guodong, Liu Fengpu.Damage constitutive model of rock under high temperature and its verification[J].Rock and Soil Mechanics,2022,43(S2):63‒73.

[55] 贾宝新,陈国栋,刘丰溥.高温下岩石损伤本构模型及其验证[J].岩土力学,2022,43(S2):63‒73.

[56] Zhao Pan.Research and application of basic performance parameters of sandstone under confining pressure[D].Chengdu:Southwest Petroleum University, 2016.

[57] 赵攀.围压下砂岩基本性能参数扩展研究及应用[D].成都:西南石油大学,2016.

基金资助

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

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

国家重点实验室自主课题(SKLGP2019Z017)

四川省科技计划项目(2019YJ0403)