多孔岩体蠕变诱发延迟压实局部化失稳特性

摘要

针对多孔岩体蠕变诱发的延迟压实局部化失稳问题，提出了一种基于弹黏塑性及可控性理论的压实局部化失稳判别准则。通过构建局部常微分方程系统，定义定常外部摄动条件下的弹黏塑性本构算子，用于识别自发传播压实局部化带内的不稳定加速变形。通过检查方程系统的可控性损失，得出了多孔岩体压实局部化的失稳指数。基于弹黏塑性本构模型，分析了平面应变条件下多孔岩体中的拟瞬时和延迟压实局部化，对所提出的失稳指数进行验证。进一步地，利用有限元数值模拟，验证了该失稳指数在分析边值问题中延迟压实局部化失稳的适用性。

关键词

压实局部化; 材料失稳; 蠕变; 弹黏塑性

自然界中多孔岩体（如砂岩）在受荷变形过程中可能发生局部化^［1］，其通常表现为塑性（或损伤）应变和位移场的局部化分布，并伴随着局部化带内孔隙塌陷，亦可称为“压实局部化失稳”^［2］。压实局部化变形模拟和预测对于地球物理和岩土工程问题具有重要意义，如斜坡、隧道的稳定性以及天然气回收、二氧化碳储存等都可能遇到与之相关的问题。

近年来有关加载诱发的局部化变形研究方面取得较大进展，室内试验中已成功观察到多孔岩体在压实过程中形成的压实带，其发生表现为孔隙塌陷和颗粒破碎为特征^［3-6］。Liu等^［7］的研究表明，多孔岩体中可出现具有非零倾角的剪切增强压实带和完全垂直于最大主应力方向的纯压实带。当前脆性岩土材料的延迟失稳也逐渐成为一个热门研究领域（如岩石^［8-9］及多孔岩体^［10］）。多孔岩体在蠕变过程中，由于其随时间变化的特性，往往会发生以加速变形为特征的自发性失稳^［10］。多孔岩体中微结构的存在导致了延迟压实局部化的形成，且伴随显著的带内应变加速演化现象。Heap等^［11］研究了蠕变试验中加速变形的时间演化和纯压实局部化带的起始和扩展。Brantut等^［10］的试验研究结果表明，具有非零倾斜角度的剪切增强型压实带也能被蠕变激活。

有关多孔岩体压实局部化的理论判别方面，通过弹塑性理论框架下的变形分叉理论或可控性分析，可预测加载导致的压实局部化起始点^［12-14］。然而，弹塑性框架不能反映多孔岩体随时间变化的力学特性，因此需要使用弹黏塑性理论。在弹黏塑性框架下，由于无法直接定义增量应力和增量应变关系的刚度矩阵，导致分叉理论和可控性理论不再适用。Pisano等^［15］通过多个常微分方程系统提供了一种识别弹黏塑性固体在稳态外部摄动（即蠕变条件）下应变加速（即应变不稳定演化）和减速（即应变稳定演化）的可能方法。结合特定的运动学约束，该方法已被成功用于分析饱和及准饱和黏土的延迟分散性型失稳^［16］以及多孔岩体中纯压实带的理论判别中^［14］。然而，现有的研究并没有考虑具有非零倾角的剪切增强压实局部化情形。

本文提出一个黏塑性压实局部化失稳指数，用于表征率相关多孔岩体中拟瞬时局部化和延迟局部化。基于黏塑性可控性理论，建立非线性常微分方程系统，将响应变量的加速度与响应变量率联系起来。推导该系统的失稳指数，用于识别自发传播的压实局部化带内不稳定加速变形响应，并通过本构模拟和有限元数值模拟对所提出失稳指数的合理性进行了验证。

1 蠕变条件下压实局部化触发条件

在定常外部摄动下（即蠕变），多孔岩体可发生以应变加速为特征的延迟压实局部化失稳。在这一过程中，压实带内的本构响应可表示为

$$\left\{\begin{array}{c}{\dot{\mathit{\sigma }}}_{\alpha }\\ {\dot{\mathit{\sigma }}}_{\beta }\end{array}\right\}=\left[\begin{array}{cc}{\mathit{D}}_{\alpha \alpha }^{\mathrm{e}}& {\mathit{D}}_{\alpha \beta }^{\mathrm{e}}\\ {\mathit{D}}_{\beta \alpha }^{\mathrm{e}}& {\mathit{D}}_{\beta \beta }^{\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{\epsilon }}}_{\alpha }-{\dot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{v}\mathrm{p}}\\ {\dot{\mathit{\epsilon }}}_{\beta }-{\dot{\mathit{\epsilon }}}_{\beta }^{\mathrm{v}\mathrm{p}}\end{array}\right\}$$

（1）

式中：${\dot{\mathit{\sigma }}}_{\alpha }$及${\dot{\mathit{\epsilon }}}_{\beta }$为加载应力率和加载应变率；相应地，${\dot{\mathit{\epsilon }}}_{\alpha }$及${\dot{\mathit{\sigma }}}_{\beta }$为响应应变率和响应应力率；${\mathit{D}}^{\mathrm{e}}$为弹性矩阵；${\dot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{v}\mathrm{p}}$及${\dot{\mathit{\epsilon }}}_{\beta }^{\mathrm{v}\mathrm{p}}$为黏塑性应变率。

考虑到广义压实带以狭窄区域内密集的局部变形为特征，因此可利用简单剪切模式来描述带内的运动变形特性，如图1所示。

图1 压实局部化带内简单剪切运动学示意图

Fig. 1 Schematic diagram of simple shear kinematics inside compaction localization bands

从图1可看出，带内变形受到法向应力率${\dot{\mathit{\sigma }}}_{\mathrm{11}}$及切向应力率${\dot{\tau }}_{\mathrm{12}}$控制；与此同时，还需要满足平面应变条件${\dot{\mathit{\epsilon }}}_{\mathrm{33}}=\mathrm{0}$、${\dot{\mathit{\epsilon }}}_{\mathrm{13}}=\mathrm{0}$、${\dot{\mathit{\epsilon }}}_{\mathrm{23}}=\mathrm{0}$及简单剪切变形约束${\dot{\mathit{\epsilon }}}_{\mathrm{22}}=\mathrm{0}$。因此，压实带内的控制加载变量率和响应变量率可表述为

$$\begin{array}{l}{\left\{{\dot{\mathit{\sigma }}}_{\alpha },{\dot{\mathit{\epsilon }}}_{\beta }\right\}}^{\mathrm{T}}={\left\{{\dot{\sigma }}_{\mathrm{11}},{\dot{\tau }}_{\mathrm{12}},{\dot{\epsilon }}_{\mathrm{22}},{\dot{\epsilon }}_{\mathrm{33}},{\dot{\gamma }}_{\mathrm{13}},{\dot{\gamma }}_{\mathrm{23}}\right\}}^{\mathrm{T}}\\ {\left\{{\dot{\mathit{\epsilon }}}_{\alpha },{\dot{\mathit{\sigma }}}_{\beta }\right\}}^{\mathrm{T}}={\left\{{\dot{\epsilon }}_{\mathrm{11}},{\dot{\gamma }}_{\mathrm{12}},{\dot{\sigma }}_{\mathrm{22}},{\dot{\sigma }}_{\mathrm{33}},{\dot{\tau }}_{\mathrm{13}},{\dot{\tau }}_{\mathrm{23}}\right\}}^{\mathrm{T}}\end{array}$$

（2）

式（1）中，黏塑性应变率通过过应力方法计算得到

$${\dot{\mathit{\epsilon }}}^{\mathrm{v}\mathrm{p}}=\Phi \left(F\right)\frac{\partial Q}{\partial \mathit{\sigma }}$$

（3）

式中：$\Phi \left(F\right)$为黏性核函数，其大小取决于非弹性应力状态（即所谓的过应力）与屈服面之间的距离；$Q\left(\mathit{\sigma }\right)$为塑性势函数。为考虑具有任意倾角的压实带，需要进一步考虑旋转参考系统，旋转角度即压实带的倾角，可通过方向余弦矩阵T对参考系统施加旋转^［16］。相应地，旋转应力率、应变率、弹性矩阵可表示为

$${\dot{\mathit{\sigma }}}^{\mathrm{R}}=\mathit{T}\cdot \dot{\mathit{\sigma }},{\dot{\mathit{\epsilon }}}^{\mathrm{R}}={\left({\mathit{T}}^{\mathrm{T}}\right)}^{-\mathrm{1}}\cdot \dot{\mathit{\epsilon }},{\mathit{D}}^{\mathrm{e},\mathrm{R}}=\mathit{T}\cdot {\mathbf{D}}^{\mathrm{e}}\cdot {\mathit{T}}^{\mathrm{T}}$$

（4）

式中

$$\mathit{T}=\left[\begin{array}{cccccc}{l}_{\mathrm{11}}^{\mathrm{2}}& l_{\mathrm{12}}^{\mathrm{2}}& l_{\mathrm{13}}^{\mathrm{2}}& \mathrm{2}{l}_{\mathrm{11}}{l}_{\mathrm{12}}& \mathrm{2}{l}_{\mathrm{12}}{l}_{\mathrm{13}}& \mathrm{2}{l}_{\mathrm{11}}{l}_{\mathrm{13}}\\ l_{\mathrm{21}}^{\mathrm{2}}& l_{\mathrm{22}}^{\mathrm{2}}& l_{\mathrm{23}}^{\mathrm{2}}& \mathrm{2}{l}_{\mathrm{21}}{l}_{\mathrm{22}}& \mathrm{2}{l}_{\mathrm{22}}{l}_{\mathrm{23}}& \mathrm{2}{l}_{\mathrm{21}}{l}_{\mathrm{23}}\\ l_{\mathrm{31}}^{\mathrm{2}}& l_{\mathrm{32}}^{\mathrm{2}}& l_{\mathrm{33}}^{\mathrm{2}}& \mathrm{2}{l}_{\mathrm{31}}{l}_{\mathrm{32}}& \mathrm{2}{l}_{\mathrm{32}}{l}_{\mathrm{33}}& \mathrm{2}{l}_{\mathrm{31}}{l}_{\mathrm{33}}\\ l_{\mathrm{11}}{l}_{\mathrm{21}}& l_{\mathrm{12}}{l}_{\mathrm{22}}& l_{\mathrm{13}}{l}_{\mathrm{23}}& \left(l_{\mathrm{11}}{l}_{\mathrm{22}}+l_{\mathrm{12}}{l}_{\mathrm{21}}\right)& \left(l_{\mathrm{12}}{l}_{\mathrm{23}}+l_{\mathrm{13}}{l}_{\mathrm{22}}\right)& \left(l_{\mathrm{11}}{l}_{\mathrm{23}}+l_{\mathrm{13}}{l}_{\mathrm{21}}\right)\\ l_{\mathrm{11}}{l}_{\mathrm{31}}& l_{\mathrm{12}}{l}_{\mathrm{32}}& l_{\mathrm{13}}{l}_{\mathrm{33}}& \left(l_{\mathrm{11}}{l}_{\mathrm{32}}+l_{\mathrm{12}}{l}_{\mathrm{31}}\right)& \left(l_{\mathrm{11}}{l}_{\mathrm{33}}+l_{\mathrm{13}}{l}_{\mathrm{31}}\right)& \left(l_{\mathrm{12}}{l}_{\mathrm{33}}+l_{\mathrm{13}}{l}_{\mathrm{32}}\right)\\ l_{\mathrm{21}}{l}_{\mathrm{31}}& l_{\mathrm{22}}{l}_{\mathrm{32}}& l_{\mathrm{23}}{l}_{\mathrm{33}}& \left(l_{\mathrm{21}}{l}_{\mathrm{32}}+l_{\mathrm{22}}{l}_{\mathrm{31}}\right)& \left(l_{\mathrm{21}}{l}_{\mathrm{33}}+l_{\mathrm{23}}{l}_{\mathrm{31}}\right)& \left(l_{\mathrm{22}}{l}_{\mathrm{33}}+l_{\mathrm{23}}{l}_{\mathrm{32}}\right)\end{array}\right]\mathit{l}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$$

（5）

$${\mathit{D}}^{\mathrm{e}}=\left[\begin{array}{cccccc}\left(K+\frac{\mathrm{4}}{\mathrm{3}}G\right)& \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \left(K+\frac{\mathrm{4}}{\mathrm{3}}G\right)& \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \left(K-\frac{\mathrm{2}}{\mathrm{3}}G\right)& \left(K+\frac{\mathrm{4}}{\mathrm{3}}G\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& G& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& G& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& G\end{array}\right]$$

（6）

式中：$\theta$为压实带倾角；$K$和$G$分别为弹性体积模量及剪切模量。

显然，式（1）-（4）描述了多孔岩体压实带内的本构特性。利用上述公式推导出响应变量率${\dot{\mathit{\epsilon }}}_{\alpha }$及${\dot{\mathit{\sigma }}}_{\beta }$为

$$\left\{\begin{array}{c}{\dot{\mathit{\epsilon }}}_{\alpha }\\ {\dot{\mathit{\sigma }}}_{\beta }\end{array}\right\}=\mathit{A}\left\{\begin{array}{c}{\dot{\mathit{\sigma }}}_{\alpha }\\ {\dot{\mathit{\epsilon }}}_{\beta }\end{array}\right\}+\Phi \mathit{B}\left\{\begin{array}{c}\frac{\partial Q}{\partial {\mathit{\sigma }}_{\alpha }}\\ \frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }}\end{array}\right\}$$

（7）

$$\begin{array}{l}\mathit{A}=\left[\begin{array}{cc}{\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}& -{\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}{\mathit{D}}_{\alpha \beta }^{\mathrm{e},\mathrm{R}}\\ {\mathit{D}}_{\beta \alpha }^{\mathrm{e},\mathrm{R}}{\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}& {\mathit{D}}_{\beta \beta }^{\mathrm{e},\mathrm{R}}-{\mathit{D}}_{\beta \alpha }^{\mathrm{e},\mathrm{R}}{\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}{\mathit{D}}_{\alpha \beta }^{\mathrm{e},\mathrm{R}}\end{array}\right]\\ \mathit{B}=\left[\begin{array}{cc}{\mathit{I}}_{\alpha \alpha }& {\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}{\mathit{D}}_{\alpha \beta }^{\mathrm{e},\mathrm{R}}\\ {\mathit{O}}_{\beta \alpha }& -\left({\mathit{D}}_{\beta \beta }^{\mathrm{e},\mathrm{R}}-{\mathit{D}}_{\beta \alpha }^{\mathrm{e},\mathrm{R}}{\left({\mathit{D}}_{\alpha \alpha }^{\mathrm{e},\mathrm{R}}\right)}^{-\mathrm{1}}{\mathit{D}}_{\alpha \beta }^{\mathrm{e},\mathrm{R}}\right)\end{array}\right]\end{array}$$

（8）

式中：${\mathit{I}}_{\alpha \alpha }$和${\mathit{O}}_{\beta \alpha }$分别为单位矩阵及零矩阵。进一步地，可推导得到响应变量的演化加速度

$$\left\{\begin{array}{l}{\ddot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{R}}={\mathit{F}}_{\epsilon }^{\mathrm{a}}+\frac{\mathrm{d}\Phi }{\mathrm{d}t}\left(\frac{\partial Q}{\partial {\mathit{\sigma }}_{\alpha }^{\mathrm{R}}}+{\mathit{B}}_{\alpha \beta }\frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\right)+\\ \Phi \left(\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial Q}{\partial {\mathit{\sigma }}_{\alpha }^{\mathrm{R}}}+\frac{\mathrm{d}}{\mathrm{d}t}{\mathit{B}}_{\alpha \beta }\frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\right)\\ {\ddot{\mathit{\sigma }}}_{\beta }^{\mathrm{R}}={\mathit{F}}_{\sigma }^{\mathrm{a}}+\frac{\mathrm{d}\Phi }{\mathrm{d}t}{\mathit{B}}_{\beta \beta }\frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}+\\ \Phi \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{B}}_{\beta \beta }\frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\end{array}\right.$$

（9）

式中：${\mathit{F}}_{\epsilon }^{\mathrm{a}}$和${\mathit{F}}_{\sigma }^{\mathrm{a}}$为体力项，在蠕变条件下为零。考虑简单剪切条件（即${\dot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{R}}=\mathrm{0}$）及定常外部摄动条件（即蠕变，${\dot{\mathit{\sigma }}}_{\beta }^{\mathrm{R}}=\mathrm{0}$），式（9）可转化为

$$\left\{\begin{array}{c}{\ddot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{R}}\\ {\ddot{\mathit{\sigma }}}_{\beta }^{\mathrm{R}}\end{array}\right\}=\mathit{Z}\left\{\begin{array}{c}{\dot{\mathit{\epsilon }}}_{\alpha }^{\mathrm{R}}\\ {\dot{\mathit{\sigma }}}_{\beta }^{\mathrm{R}}\end{array}\right\}+\mathit{F}$$

（10）

$$\mathit{Z}=\left[\begin{array}{cc}-\Phi \frac{\partial \Phi }{\partial f}{H}_{\mathrm{i}\mathrm{n}}{\mathit{I}}_{\alpha \alpha }& \frac{{\partial }^{\mathrm{2}}Q}{\partial {\mathit{\sigma }}_{\alpha }^{\mathrm{R}}\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}+{\mathit{B}}_{\alpha \beta }\frac{{\partial }^{\mathrm{2}}Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\\ {\mathit{O}}_{\beta \alpha }& -\Phi \frac{\partial \Phi }{\partial f}{H}_{\mathrm{i}\mathrm{n}}{\mathit{I}}_{\beta \beta }-{\mathit{A}}_{\beta \beta }\frac{{\partial }^{\mathrm{2}}Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\end{array}\right]$$

（11）

在蠕变条件下有

$$H_{\mathrm{i}\mathrm{n}}=H-H_{\mathrm{\chi }},H_{\mathrm{\chi }}=-{\left\{\frac{\partial F}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}\right\}}^{\mathrm{T}}{\mathit{A}}_{\beta \beta }\frac{\partial Q}{\partial {\mathit{\sigma }}_{\beta }^{\mathrm{R}}}$$

（12）

式中：$H$为硬化模量；$H_{\mathrm{\chi }}$与式（2）、式（3）中的控制加载变量及响应变量的具体形式相关，亦可称为可控性模量。

最终，控制响应变量变速演化规律的常微分方程系统得以建立。式（10）通过依赖时间变化的矩阵$\mathit{Z}$（即弹黏塑性本构算子）控制响应变量的变速演化规律。若矩阵$\mathit{Z}$的特征值为负，则响应变量将呈现减速演变，可视为稳定响应；然而若矩阵$\mathit{Z}$的特征值为正，将导致响应变量速率的加速演变，可视为不稳定响应。式（10）表明，矩阵$\mathit{Z}$的特征值由${\mathit{Z}}_{\alpha \alpha }$和${\mathit{Z}}_{\beta \beta }$的特征值组成。${\mathit{Z}}_{\alpha \alpha }$为对角矩阵且特征值的正负号由$H_{in}$控制；而${\mathit{Z}}_{\beta \beta }$为一个对角矩阵和一个海森矩阵之和。由于塑性势函数为凸函数，海森矩阵的特征值为恒正，因此${\mathit{Z}}_{\beta \beta }$的特征值永远小于${\mathit{Z}}_{\alpha \alpha }$的特征值。这就意味着当$H_{\mathrm{i}\mathrm{n}}>\mathrm{0}$时，响应变量呈现减速演化趋势，即该系统是恒稳定的；相反，当$H_{\mathrm{i}\mathrm{n}}<\mathrm{0}$时，响应变量速率有加速发展的潜在可能，意味着系统是潜在不稳定的。因此，$H_{\mathrm{i}\mathrm{n}}>\mathrm{0}$是弹黏塑性材料稳定的充分条件，而$H_{\mathrm{i}\mathrm{n}}<\mathrm{0}$则对应潜在压实局部化失稳。

2 多孔岩体压实局部化理论预测

2.1　多孔岩体弹黏塑性模型

将上述理论与一种应变硬化弹黏塑性本构模型相结合，该模型由Nova等^［17］提出，包含了硬化和软化的相互竞争机制。屈服函数（当$h=F$时，代表屈服函数）及塑性势函数（当$h=Q$时，代表塑性势函数）表示为

$$h=\left[{\left(\mathrm{1}+\frac{{\eta }^{\mathrm{*}}}{K_{\mathrm{1}h}{M}_h}\right)}^{K_{\mathrm{1}h}/\left(\mathrm{1}-{\mu }_h\right)\left(K_{\mathrm{1}h}-K_{\mathrm{2}h}\right)}\cdot {\left(\mathrm{1}+\frac{{\eta }^{\mathrm{*}}}{K_{\mathrm{2}h}{M}_h}\right)}^{K_{\mathrm{2}h}/\left(\mathrm{1}-{\mu }_h\right)\left(K_{\mathrm{1}h}-K_{\mathrm{2}h}\right)}\right]p-p_c$$

（13）

$$K_{\mathrm{1}h,\mathrm{2}h}=\frac{{\mu }_h\left(\mathrm{1}-{\alpha }_h\right)}{\mathrm{2}\left(\mathrm{1}-{\mu }_h\right)}\left(\mathrm{1}\pm \sqrt[]{\mathrm{1}-\frac{\mathrm{4}{\alpha }_h\left(\mathrm{1}-{\mu }_h\right)}{{\mu }_h{\left(\mathrm{1}-{\alpha }_h\right)}^{\mathrm{2}}}}\right)$$

（14）

式中：$p={\sigma }_{kk}/\mathrm{3}$为球应力；$q=\sqrt[]{\mathrm{3}{s}_{ij}{s}_{ij}/\mathrm{2}}$为广义等效剪应力；$s_{ij}={\sigma }_{ij}-p{\delta }_{ij}$为偏应力；${\eta }^{\mathrm{*}}=q/\left(p+p_{\mathrm{t}}\right)$为应力比；$p_{\mathrm{c}}=p_{\mathrm{s}}+p_{\mathrm{m}}+p_{\mathrm{t}}$为各向同性压缩条件下的屈服应力，$p_{\mathrm{t}}=\kappa p_{\mathrm{m}}$，其中$p_{\mathrm{s}}$和$p_{\mathrm{m}}$为硬化变量，分别反映了骨架堆积和晶间胶结的贡献，$\kappa$控制了屈服面或塑性势面在拉力域内的扩展。

通过标定参数$M_h$、${\mu }_h$及${\alpha }_h$，可灵活控制式（13）中给出的屈服面和塑性势面的形状。模型中的硬化变量率${\dot{p}}_{\mathrm{s}}$和${\dot{p}}_{\mathrm{m}}$可表达为

$${\dot{p}}_{\mathrm{s}}=\frac{p_{\mathrm{s}}}{B_{\mathrm{p}}}{\dot{\epsilon }}_{\mathrm{v}}^{\mathrm{p}},{\dot{p}}_{\mathrm{m}}=-{\rho }_{\mathrm{m}}{p}_{\mathrm{m}}\left(\left|{\dot{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}\right|+{\xi }_{\mathrm{m}}{\dot{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}\right)$$

（15）

式中：$B_{\mathrm{p}}$、${\rho }_{\mathrm{m}}$及${\xi }_{\mathrm{m}}$为材料参数。$p_{\mathrm{s}}$可用于捕捉由于孔隙变化导致的硬化或软化机制，其中塑性体积应变反映孔隙的变化。多孔岩体压缩变形时，颗粒孔隙会明显减小，这一物理过程可由$p_{\mathrm{s}}$反映，而$p_{\mathrm{m}}$则同时通过体积塑性应变和等效塑性剪应变引入软化机制。

利用式（3），该模型可转变为弹黏塑性模型，黏塑性应变可通过式（16）求得：

$${\dot{\epsilon }}_{ij}^{\mathrm{v}\mathrm{p}}=\Phi \left(F\right)\frac{\partial Q}{\partial {\sigma }_{ij}}=\frac{\mathrm{1}}{\omega }\left[\frac{F}{p_{\mathrm{c}\mathrm{0}}}\right]\frac{\partial Q}{\partial {\sigma }_{ij}}$$

（16）

式中：$p_{\mathrm{c}\mathrm{0}}$为p_c的初始值；$\omega$为黏性参数。模型的弹性部分采用线弹性模型。

根据上述模型，Marinelli和Buscarnera^［18］基于三轴试验结果及变形分叉理论，标定了多种多孔岩体材料参数，标定的Berea砂岩材料参数如表1所示。基于所建立的本构模型开展有限元数值模拟，在不同围压（p₀=75 MPa、150 MPa、200 MPa以及 300 MPa）下模拟得到的力学响应与三轴试验结果对比如图2所示，两者较为符合。

表1 Berea砂岩材料参数^［18］

Tab. 1 Material parameters of Berea sandstone^[18]

$K$/MPa	G/MPa	${\rho }_{\mathrm{m}}$	${\xi }_{\mathrm{m}}$	$B_{\mathrm{p}}$	${\mu }_{\mathrm{f}}$	${\alpha }_{\mathrm{f}}$	$M_{\mathrm{f}}$	$\kappa$	${\mu }_g$	${\alpha }_g$	$M_g$	$p_{s0}$/MPa	$p_{m0}$/MPa	$\omega$/（MPa·s）
9 550	7 000	3.7	0.5	0.04	1.3	0.2	1.02	0.025	2	0.2	1.88	80	295	100

图2 Berea砂岩力学响应

Fig. 2 Mechanical response of Berea sandstone

2.2　平面应变单元压实局部化分析

进一步地，对平面应变状态下延迟压实局部化进行分析。将平面应变单元视为一个积分点，验证所提出的黏塑性压实局部化失稳指数的适用性。分析包括：①各向同性压缩，使模型具有初始围压；②平面应变剪切模拟，使模型进入非弹性状态；③改变约束条件，使模型进入蠕变变形状态。对于分析步②，当应力状态达到初始屈服面时，对应于某个角度的失稳指数可能变为零乃至负值（即$H_{\mathrm{i}\mathrm{n}}\left(\theta \right)\le \mathrm{0}$），此角度即为可能发生的压实局部化带倾角。当存在$H_{\mathrm{i}\mathrm{n}}\left(\theta \right)\le \mathrm{0}$时，$H_{\mathrm{i}\mathrm{n}}\left(\theta \right)$最小值对应的角度为最可能出现的局部化带倾角，即$A\left({\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)=\mathrm{m}\mathrm{i}\mathrm{n}A\left(\theta \right)$。在各种可能的局部化变形中，$A\left(\theta =\mathrm{0}\right)<\mathrm{0}$时发生的变形称为压实局部化变形。这种情况下，倾斜角度为0°的局部化带成为一种潜在可能。这种压实局部化变形又分为2种情况，其一为$A\left({\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)<A\left(\theta =\mathrm{0}\right)$，具有非0°倾斜角的剪切增强压实带；其二为 $A\left({\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)<A\left(\theta =\mathrm{0}\right)$，具有0°倾斜角的纯压实带。在随后的分析步③中，将约束条件更改为与该倾角对应方向施加简单剪切蠕变变形。应注意到，分析步③所施加的控制条件体现了压实变形带内的变形运动特征。

为研究当应力路径达到不同屈服面区域后的模型响应及失稳指数变化规律，选择p₁=105 MPa、p₂=190 MPa这2种围压开展Berea砂岩模拟，计算结果如图3所示。应力路径与屈服面的交点分别为屈服点Ⅰ和屈服点Ⅱ，由于黏塑性理论允许应力状态超越屈服面，采用围压p₁和p₂计算所得应力状态与屈服面变化的不匹配造成了变形特性的差异。当应力状态达到屈服面时，应力状态增长速度大于屈服面扩展速度，此阶段为不稳定蠕变变形阶段，可产生由加速变形造成的蠕变压实局部化失稳；而当应力状态增长速度逐渐小于屈服面扩展速度时，不稳定蠕变响应转化为稳定蠕变响应，此时出现的压实局部化为稳定变形。针对屈服点Ⅰ和Ⅱ的所有倾角失稳指数计算结果如图4所示，从图中可看出，基于所提出失稳指数的压实局部化理论预测与基于传统分叉理论的预测结果相匹配。这2种理论所得结果表明：屈服点Ⅰ点的应力状态对应倾角为±21°的剪切增强压实带；屈服点Ⅱ点预测得到倾角为0°的纯压实带。蠕变阶段轴向应变的演化如图5a所示，轴向蠕变应变最开始呈现为加速不稳定发展，随后应变率逐渐降低至趋于零，轴向应变达到稳定值。根据图3结果可知，应力状态和屈服面演化速率不匹配造成了应变加速（或减速）演化，因而应力状态增长速度等于屈服面扩展速度时对应的时间点即为从不稳定蠕变到稳定蠕变的转折点（实心点标出）。图5b表明，围压p₁和p₂条件下的失稳指数最初为负值，亦即式（8）中弹黏塑性算子的特征值为正，表明轴向应变呈加速变化趋势，即发生不稳定蠕变。随着蠕变变形发展，特征值由负值转化为正值，即逐渐进入减速变形阶段。图5b中失稳指数正负号的转折点与图5a不稳定蠕变转折点一致，说明所提出的黏塑性压实局部化失稳指数能有效识别蠕变的加速失稳和减速稳定阶段。

图3 材料点应力路径分析

Fig. 3 Stress paths of material point analysis

图4 分叉理论与失稳指数计算结果对比

Fig. 4 Comparison of the bifurcation theory and the instability index

图5 Berea砂岩剪切蠕变分析

Fig. 5 Analysis of simple shear creep of Berea sandstone

3 多孔岩体压实局部化有限元模拟

进一步开展Berea砂岩局部化形成过程的平面应变有限元数值模拟。分析模型及网格如图6a所示，模型长宽比为2：1，由12 800个C3D8单元组成。底部边界施加全约束，顶部按应变速率10^-5s^-1施加荷载，以使区域内每个高斯点的应力状态达到屈服面，然后将顶部的牵引力保持恒定，以确保试样可以在固定的外部扰动下随时间蠕变。值得注意的是，这里有限元计算模拟所采用的参数、轴压、控制条件以及加载量与第3节点积分理论分析完全相同。为了触发应变局部化，在模型中间设置了一个薄弱区域（将p_c取为原值的95%）。

图6 Berea砂岩平面应变压缩试验模拟

Fig. 6 Numerical simulations of plane strain compression test for Berea sandstone

Berea砂岩有限元数值模拟结果如图6b。由图可见，一旦应力状态进入塑性区（T / T_c = 0），在围压p₁和p₂条件下计算得到的塑性体积应变场呈压实局部化变形模式，且得到的局部化带倾角与材料点预测结果完全一致。图6b还体现了恒定轴向应力作用下体积塑性应变的空间传播模式。在围压p₁和p₂条件下所得的体积塑性应变首先从薄弱区域出现，然后在模型内部形成多条压实带，这种压实带的形成意味着可能发生局部化失稳。从位于压实局部化带内的高斯积分点获得的局部响应如图7所示，围压p₁和p₂计算结果中提取的高斯点（图6中灰色星号示出）显示出最初的加速趋势（即不稳定）以及随后的减速趋势（即稳定）。塑性体应变区域内的稳定和不稳定应变响应分别与对应倾角的失稳指数的正负号有关（即倾角为21°和0°分别对应$H_{\mathrm{i}\mathrm{n}}\left(\theta ={\mathrm{21}}^{\mathrm{o}}\right)$ 和$H_{\mathrm{i}\mathrm{n}}\left(\theta ={\mathrm{0}}^{\mathrm{o}}\right)$）。分析结果表明，所提出的黏塑性压实局部化失稳指数能很好地用于率相关多孔岩体中压实局部化的分析。

图7 Berea砂岩有限元高斯积分点分析

Fig. 7 Finite element Gaussian point analysis of Berea sandstone

4 结语

基于黏塑性假设及压实带内变形运动学特性，建立了描述多孔岩体蠕变条件下应变局部化响应变量随时间演化的常微分方程组，并提出了一个判别多孔岩体蠕变诱发延迟局部化的黏塑性压实局部化失稳指数。进一步地，采用能捕捉不同应变局部化模式的弹黏塑性应变硬化及服从非关联流动法则的本构模型，开展了材料点失稳分析及有限元模拟，结果表明：①通过建立的理论，在加载阶段不同压实应变局部化模式的预测结果与率无关弹塑性变形分叉理论预测结果一致；②蠕变过程模拟中响应变量的加速和减速是由超过弹性域的应力状态与屈服面之间的相对运动引起的，而失稳指数的符号变化（由负到正）对应了响应变量由加速向减速演化的转折点。这些结果表明，文中提出的失稳指数能够合理判别多孔岩体蠕变压实局部化现象。

作者贡献声明

薛大为：理论推导、程序编写及稿件撰写。

吕玺琳：提出研究思路和方法、稿件审核。

任中俊：结果校核、稿件撰写。

刘泳钢：方案设计、稿件修订。