0.   引言

大量地下工程由于围岩随时间推移产生较大蠕变变形而出现失稳破坏，富水环境下围岩更易发生蠕变变形^[1]，并导致渗透性变化，诱发渗漏水事故，因此研究蠕变过程中渗透率演化规律对防治地下工程灾害有重要理论与工程意义。

国内外学者从试验角度针对蠕变、渗透率演化开展了深入研究。Liu等^[2-3]、Zhang等^[4]、王如宾等^[5]、张玉等^[6]分别研究了泥岩、碎屑砂岩、变质角砾岩、破碎带砂岩蠕变过程中渗透率演化规律。蔡婷婷等^[7]研究了不同温度与应力水平条件下煤岩蠕变过程中渗流规律，指出温度升高引起煤岩蠕变起始阈值降低、并加剧煤岩蠕变变形。Xu等^[8]、Yang等^[9]研究指出蠕变过程中渗透率与体应变具有较好一致性，并将二者随时间变化过程划分了4个阶段。马丹等^[10]针对断层破碎带开展了蠕变-冲蚀耦合试验，研究结果表明突水可分为初始期、冲蚀主控期、蠕变主控期，蠕变主控期间岩体变形导致围岩失稳。王凯等^[11]开展了不同路径下含瓦斯复合煤岩体渗透率演化规律。岳少飞等^[12]开展了不同加载速率条件下无烟煤的分级加载蠕变试验，研究了硬化损伤特征。

在理论方面，于冰冰等^[13]根据时间与应力将蠕变曲线划分为瞬时蠕变、减速蠕变、等速蠕变及加速蠕变阶段，引入瞬时塑性体、村山体、CYJ体建立了一维及三维非线性蠕变模型以描述蠕变各阶段特征。程爱平等^[14-15]针对胶结充填体开展了分级蠕变试验，并监测了超声波变化规律，利用波速定义了损伤，构建了考虑硬化-损伤特征的蠕变本构模型。Liu等^[16]基于统计损伤理论，提出了真三轴各向异性蠕变模型。康红普等^[17]将改进的伯格斯体、黏塑性体和应变软化塑性体串联，以描述围岩稳定与非稳定蠕变、流变损伤和流变扩容特征。贾善坡等^[18]通过开展室内试验分析了蠕变过程中黏土岩渗透率变化规律，基于K-C方程提出了渗透率表达式，并将该公式应用至实际地下工程进行了数值计算。李祥春等^[19]研究了蠕变过程中煤岩渗透率变化规律，提出了以轴向应变为自变量的渗透率经验表达式，参数较少但缺少一定理论参考。Zhou等^[20]分析了煤岩卸围压蠕变过程中渗透率变化规律，基于K-C方程以体应变为自变量建立了渗透率表达式，具有较高精度，但由于以应变为自变量难以反映其对应的应力状态。Zhou等^[21]考虑了基质与裂隙间相互作用，并引入内膨胀系数，建立了黏弹性阶段渗透率模型。刘帅奇等^[22]考虑了裂缝蠕变变形对渗透率影响，得到非线性黏弹塑性模型下裂隙面刚度表达式，推导出渗透率随时间的演化规律。张雷等^[23]结合考虑体积蠕变的分数阶蠕变模型，考虑了黏塑性变形与裂隙开度关系，基于火柴棍模型，建立了考虑蠕变影响的渗透率模型。王路军等^[24]通过相似模拟试验与三轴循环加卸载渗流耦合试验，提出了突渗和拓扑的力学定义，建立了突渗的逾渗模型。亓宪寅^[25]提出了各向异性蠕变-渗透率模型，该模型将蠕变段视为初始值，且仅考虑了黏弹性阶段，适用范围有限。

现有蠕变-渗透率模型多基于各向同性假设，以孔隙度、轴向应变、体变为自变量。实际上三向应力下渗透率呈各向异性特征，忽略渗透率各向异性变化易造成渗流路径预测不准，诱发渗漏水事故。蠕变应力与时间是渗透率变化的间接原因，二者引起岩石内部孔/裂隙结构闭合或连接、贯通进而渗流通道演化是渗透率各向异性变化的直接原因。因此，本文首先建立各向异性蠕变损伤模型以描述岩石蠕变特征，再以应变为桥梁，从三向应力下裂隙开度变化出发，建立正交各向异性蠕变-渗透率模型，为预测渗流路径、渗流量提供理论基础，对地下工程设计、施工及长期稳定有重要意义。

1.   各向异性蠕变损伤模型

典型蠕变曲线一般包括三阶段，瞬时蠕变、稳定蠕变、加速蠕变，西原模型能描述蠕变各阶段特征而被广泛使用，本节基于西原模型，建立岩石各向异性蠕变损伤模型。对于应力引起的瞬时应变，根据广义胡克定律，瞬时应变为

式中：$\boldsymbol{\varepsilon} _{ij}^{\mathsf{ins}}$为瞬时应变张量；${\boldsymbol{C}_{ijkl}}$为弹性刚度矩阵；${\boldsymbol{\sigma} _{kl}}$为应力张量。

对于黏弹性阶段，提出黏弹性参数具有各向异性特征，黏弹性阶段蠕变应变张量为

式中：$\boldsymbol{\varepsilon} _{ij}^\mathsf{e}$为黏弹性应变张量；S_ij为应力偏张量；$G_{ij}^\mathsf{e}$，$\eta _{ij}^\mathsf{e}$分别为黏弹性阶段剪切模量和黏滞系数。

当蠕变应力水平较高甚至超过长期强度，蠕变曲线呈典型非线性特征，西原模型仅能描述黏塑性阶段线性变化，而Pernzyna过应力理论能较好描述蠕变非线性趋势^[26-27]，因此，本文采用Pernzyna过应力理论定义黏塑性应变以代替经典西原模型中黏塑性元件模型，其黏塑性应变率表达式为

式中，$\langle f\rangle= \begin{cases}\left(\frac{f}{f_0}\right)^a & f>0 \\ 0 & f \leqslant 0\end{cases}$，其中，f₀为屈服函数的初始值，可取为1；a为常数，且一般取1。f为过应力，与屈服函数相关，常用的屈服函数有主要有Mohr- Coulomb、Hoek-Brown(H-B)、Drucker-Prager(D-P)、Mogi-Coulomb等。其中，Mohr-Coulomb与H-B准则忽略了中间主应力的影响，D-P准则考虑了中间主应力但偏保守，Mogi-Coulomb的适用性较好^[28]，其具体表达式如下：

式中，${\tau _{{\text{oct}}}}$为八面体剪应力，

式（3）中$\frac{{\partial Q}}{{\partial {\sigma _{ij}}}}$体现了流动方向，Q为势函数，为便于计算采用关联流动法则，即Q=f。

在黏塑性阶段，岩石裂隙萌生、发育、汇聚、连通，损伤程度逐渐加剧，变形随之发展。因此引入蠕变损伤变量表示岩石内部损伤累积，考虑损伤的黏滞系数为：

式中：$\eta _{ij}^{\text{p}}$，$\eta _{ij}^{\text{p}}(D_{ij}^{\text{p}})$分别为初始及考虑损伤的黏滞系数。

统计损伤理论是通过引入代表性微元体概念，假设岩石内部孔裂隙呈非均匀分布，采用损伤微元体数量与微元体总量描述材料损伤程度。假设岩石微元体强度服从Weibull分布，概率密度函数为^[29-32]

式中：m，F₀为分布参数；F为微元体强度，可选取不同强度准则，考虑到采用应变表示的便利性，本文以应变作为分布变量定义损伤变量，概率密度函数表示为

微元体不断破坏引起了岩石材料损伤，将荷载作用下破坏的微元体数量记为${N_{\text{f}}}$，损伤变量可以表示为破坏的微元体数量与微元体总量的比：

式中：D为损伤变量；N为总微元体数量。

在任意区间[ε，ε+dε]，破坏的微元体数量为Nf(ε)dε，当加载至任意应力水平，破坏的微元体数量为

损伤变量表示为

三轴压缩条件下侧向应变为张拉应变，一般表示为负值，为保证分布变量为正值，以各向黏塑性应变的绝对值作为分布变量，具体表示为

式中：$D_{ij}^{\text{p}}$为蠕变损伤变量；${\text{ }}m_{ij}^{\text{p}}$，$\varepsilon _{ij0}^{\text{p}}$分别为蠕变分布参数。

式（12）不便于积分求解蠕变应变，文献[27]指出可采用多项式近似表示蠕变应变，各向蠕变应变表示为

式中：n为多项式幂指数；${a_{ijn}}$为常数。

对于稳定蠕变阶段，蠕变速率衰减至稳定，式（13）简化为

对于加速蠕变阶段，式（13）简化为

将式（6），（12），（15）代入式（3），则黏塑性应变率为

令式（16）中${a_{ij0}}/\varepsilon _{ij}^{\text{p}} = \alpha _{1,ij}^{}$，${a_{ij1}}/\varepsilon _{ij0}^{\text{p}} = \alpha _{2,ij}^{}$，式（16）可改写为

综上，岩石各向异性蠕变损伤模型为

2.   各向异性蠕变-渗透率模型

岩石由基质与内部大量孔/裂隙构成，有学者将其视为双重介质，认为其具有两种渗透率与孔隙率，也有学者认为基质渗透率太低，因此忽略基质渗透率，仅考虑裂隙渗透率。本节围绕渗透率各向异性特征，仅考虑裂隙渗透率，假设应力、应变、渗透率主轴重合，建立正交各向异性蠕变-渗透率模型。

如图 1所示，将岩石化简为立方体模型，基质被正交裂隙均匀分割为正方体形式的单元体，基质单元初始宽度为a，各向裂隙宽度分别为b_x，b_y，b_z。

图  1  岩石简化示意图

Figure  1.  Simplified schematic diagram of rock

下载: 全尺寸图片 幻灯片

岩石初始孔隙率可表示为

式中：${\phi _0}$为初始孔隙率，${b_{x0}}$，${b_{y0}}$，${b_{z0}}$为初始裂隙宽度。对于初始各向同性岩石，b_x0= b_y0=b_z0=b₀。

2.1   黏弹性阶段渗透率模型

以单裂隙为例推导三向应力下裂隙开度变化，σ_y为该裂隙的法向应力，σ_x与σ_z为侧向应力，此处定义裂隙压缩为正、拉伸为负，如图 2所示。

图  2  三维应力状态下裂隙渗流模型

Figure  2.  Seepage model for a fracture under 3D stress

下载: 全尺寸图片 幻灯片

在法向应力σ_y作用下裂隙开度为^[33]

式中：${b_{\text{y}}}$为在法向应力作用下的裂隙开度；${b_0}$为初始裂隙开度；$E_{}^{\text{f}}$为裂隙弹性模量；$\varepsilon _{\text{y}}^{\text{f}}$为法向应力作用下裂隙法向应变。将裂隙弹性模量与岩石弹性模量建立联系，表示为

式中：${\xi ^{\text{E}}}$为裂隙与岩石的弹性模量比。

式（20）中裂隙法向应变可以用岩石应变表示：

将该线性关系推广至蠕变阶段，黏弹性阶段裂隙法向应变为

式中：$\varepsilon _{\text{y}}^{{\text{fe}}}$为黏弹性阶段裂隙法向应变；$\xi _{\text{y}}^{\text{e}}$为黏弹性阶段裂隙与岩石应变比值；$\varepsilon _{\text{y}}^{\text{e}}$为黏弹性阶段岩石蠕变应变。

参考文献[34]的双重介质蠕变模型，可证明式（23）合理性。裂隙的黏弹性应变张量表示为

式中：$G_{ij}^{{\text{fe}}}$，$\eta _{ij}^{{\text{fe}}}$分别为裂隙黏弹性剪切模量及黏滞系数。

根据式（2），（24），裂隙与岩石蠕变应变比值表示为

对式（25）中指数项进行一阶泰勒展开，即

将其代入式（25），化简得

式（27）与式（23）一致，可知参数$\xi _{ij}^{\text{e}}$有明确物理意义，表示裂隙与岩石的黏滞系数比。

此处以y方向为例，式（27）可表示为

综上，黏弹性阶段法向应力作用下裂隙开度的变化量为

侧向应力作用下泊松效应使岩石产生法向膨胀变形，裂隙受到挤压、开度减小。学者对于侧向应力对裂隙开度影响尚未形成统一认识。有学者将侧向应力等效为法向应力，认为其直接作用于裂隙，可以根据应力引起基质、岩石应变求解得到裂隙应变^[35]，文献[36]指出侧向应力能够引起法向应变但不会产生等效法向应力。有学者根据弹性力学分别计算应力引起的基质与岩体应变，根据二者关系求解裂隙开度变化，这一方法目前比较常见，忽略基质泊松比与岩石泊松比的差异，得到渗透率随侧向应力增加而增加的趋势，引起较大误差并且难以解释试验现象。实际上，侧向应力增加引起渗透率下降。文献[37]认为侧向应力引起的岩石闭合类似于吸附膨胀引起的裂隙内闭合，引入影响系数以体现骨架变形对裂隙变形的影响。

综上，笔者认为侧向应力对裂隙开度影响与法向应力类似，引入表示侧向应力对裂隙开度的侧向影响参数，侧向应力下裂隙开度变化量为

式中：$\Delta b_{{\text{sy}}}^{\text{e}}$为黏弹性阶段侧向应力影响下裂隙开度变化量；$\beta _{\text{y}}^{\text{e}}$为侧向影响系数；$\varepsilon _x^{\text{e}}$和$\varepsilon _z^{\text{e}}$分别为x，y方向的黏弹性应变。

在侧向应力下裂隙的法向应变为

根据立方定律，得z方向黏弹性阶段渗透率表达式为

将式（28），（31）代入式（32），可得黏弹性阶段各向渗透率表达式为

2.2   黏塑性阶段渗透率模型

当蠕变进入黏塑性阶段，岩石内部裂隙萌生、发育、连通。部分裂隙与主裂隙间连通性较差，如图 3中的死裂隙，对形成有效渗流通道贡献较少；而一些裂隙与其他裂隙连通，如图 3中融合裂隙，形成有效连通裂隙，是主要渗流通道。因此，引入修正系数表示裂隙开度对渗透通道的影响，基于岩石黏塑性应变与裂隙开度的密切关系^[23]，建立黏塑性阶段渗透率模型。

图  3  不同类型裂隙

Figure  3.  Different types of cracks

下载: 全尺寸图片 幻灯片

黏塑性阶段裂隙开度的变化量表示为

式中：$\omega _y^{\text{p}}$为修正系数；$\varepsilon _y^{\text{p}}$为岩石黏塑性应变；$\Delta b_y^{\text{p}}$为黏塑性阶段裂隙开度变化量。

黏塑性阶段裂隙应变表示为

式中：$\omega _y^{\text{p}}$为修正系数；$\varepsilon _y^{{\text{fp}}}$为黏塑性阶段裂隙应变。

根据立方定律，黏塑性阶段的渗透率表示为

3.   渗透率模型应用分析

3.1   渗透率模型验证

为了验证各向异性蠕变-渗透率模型有效性，笔者以砂岩为研究对象，开展了真三轴条件下蠕变-渗流试验，采用逐级加载方式，渗透压力为3 MPa，具体应力水平设置、蠕变损伤模型参数详见文献[16]。

利用通用全局优化算法，根据蠕变曲线确定蠕变损伤模型参数，再根据渗透率曲线确定其他参数，如表 1所示。为了对比本文提出的各向异性渗透率模型与现有渗透率模型优劣，将本文模型与经典的K-C模型^[18]及文献[25]中考虑蠕变的各向异性渗透率模型对比，以$\sigma _{\text{3}}^{}$=5 MPa，$\sigma _{\text{2}}^{}$=20 MPa条件下最大主应力方向渗透率k₁、中间主应力方向渗透率k₂为例，模型对比结果如图 4。在黏弹性、黏塑性阶段，根据K-C模型与文献[25]考虑蠕变的各向异性渗透率模型计算得到的理论值与试验值偏差较大，而本文提出的各向异性蠕变-渗透率模型能较好描述黏弹性、黏塑性阶段渗透率响应变化，体现了本文提出的各向异性蠕变-渗透率模型优势显著。

表  1  σ₃=5 MPa，σ₂=20 MPa条件下渗透率模型参数

Table  1.  Parameters of permeability model under σ₃=5 MPa and σ₂=20 MPa

蠕变应力比	ξe	$\beta _1^{\text{e}}$	$\beta _2^{\text{e}}$	$\beta _3^{\text{e}}$	$\omega _1^{\text{p}}$	$\omega _2^{\text{p}}$	$\omega _3^{\text{p}}$
0.4	18	9.81	9.81	9.81
0.5		4.18	4.95	4.95
0.6		2.06	4.30	4.30
0.7		0.94	3.37	3.37
0.8					11.01	13.85	13.85
0.9					8.72	8.72	8.72
0.95					0.81	2.11	2.11

下载: 导出CSV 

| 显示表格

图  4  不同渗透率模型对比图

Figure  4.  Comparison among different permeability models

下载: 全尺寸图片 幻灯片

从图 5中可以看出，黏弹性阶段蠕变应力比为0.4时，渗透率理论值k₁，k₂分别从2.77×10^-18，3.54×10^-18 m²降低至2.28×10^-18，3.07×10^-18 m²；当蠕变应力比为0.7时，渗透率理论值k₁，k₂分别从1.8×10^-18，2.6×10^-18 m²降低至1.6×10^-18，2.5×10^-18 m²，侧向影响系数$\beta _{\text{1}}^{\text{e}}$从9.81降低至0.94，$\beta _{\text{1}}^{\text{e}}$从9.81降低至3.37。当f > 0，进入黏塑性蠕变阶段，中主应力方向上黏塑性应变相对较小，对渗透率贡献较小，假设该方向黏塑性应变对渗透率影响参数与最小主应力方向影响参数相同，即$\omega _{\text{2}}^{\text{p}} = \omega _{\text{3}}^{\text{p}}$。当蠕变应力比为0.8时，渗透率理论值k₁，k₂分别由1.75×10^-18，2.62×10^-18 m²增加至1.93×10^-18，2.74×10^-18 m²，加速蠕变阶段时，裂隙迅速发展、形成宏观破裂面，渗透率理论值k₁、k₂急剧增加，由2.11×10^-18，2.85×10^-18 m²增加至1.7×10^-17，2.3×10^-17 m²，增大了一个数量级，修正系数$\omega _{\text{1}}^{\text{p}}$，$\omega _{\text{2}}^{\text{p}}$，$\omega _{\text{3}}^{\text{p}}$随着蠕变应力比增加而降低。本模型可以描述黏弹性阶段由于裂隙闭合致使渗透率减小的趋势以及黏塑性阶段也由于裂隙连通引起的渗透率激增趋势，并且模型也能够描述渗透率各向异性特征。

图  5  σ₃=5 MPa，σ₂=20 MPa蠕变条件下渗透率曲线

Figure  5.  Permeability curves under of σ₃=5 MPa and σ₂=20 MPa

下载: 全尺寸图片 幻灯片

3.2   渗透率模型参数敏感性分析

以σ₃=5 MPa，σ₂=20 MPa条件为例，分析渗透率模型黏弹性阶段参数ξ^e、$\beta _{\text{1}}^{\text{e}}$黏塑性阶段修正参数$\omega _{\text{1}}^{\text{p}}$对渗透率演化曲线影响，如图 6，7所示。

图  6  黏弹性阶段参数ξ^e，$\beta _{\text{1}}^{\text{e}}$对渗透率演化曲线影响

Figure  6.  Effects of parameters ξ^e and $\beta _{\text{1}}^{\text{e}}$ on permeability curves at viscoelastic stages

下载: 全尺寸图片 幻灯片

图  7  黏塑性阶段修正系数对渗透率演化曲线影响

Figure  7.  Effects of correction parameter on permeability curves at viscoplastic stages

下载: 全尺寸图片 幻灯片

从图 6（a）中可以看出，随着ξ^e增加，裂隙黏滞系数越小，在三向应力下裂隙越容易发生变形，促使裂隙闭合程度增加，抑制了岩石中水的流动, 渗透率衰减越快，蠕变末时刻稳定渗透率从2.32×10^-18 m²减小至1.85×10^-18m²。本节以侧向影响系数$\beta _{\text{1}}^{\text{e}}$为例讨论侧向应力对渗透率演化曲线的影响，从图 6（b）可以看出$\beta _{\text{1}}^{\text{e}}$越高渗透率初始衰减速率越高，稳定渗透率越低，表明侧向压力对渗透率的影响越大，随着$\beta _{\text{1}}^{\text{e}}$从7增加至12，稳定渗透率从2.41×10^-18 m²降低至2.15×10^-18 m²。在黏塑性阶段，以修正系数$\omega _{\text{1}}^{\text{p}}$为例分析$\omega _{\text{1}}^{\text{p}}$对渗透率影响，结果如图 7所示。$\omega _{\text{1}}^{\text{p}}$越高，稳定蠕变阶段渗透率增加幅度不大，随着岩石进入加速蠕变阶段渗透率开始大幅增加，随着$\omega _{\text{1}}^{\text{p}}$从7.81增加至11.81破坏时渗透率从7.2×10^-18 m²增加至17.5×10^-18 m²。

4.   结论

通过结合各向异性蠕变损伤模型、立方定律，建立了正交各向异性蠕变-渗透率模型，并验证模型的准确性，分析模型关键参数影响，得出以下3点结论。

（1）黏弹性阶段分别考虑法向应力与侧向应力对裂隙开度的影响，引入黏滞系数比、侧向应力影响系数，模型能够描述渗透率因裂隙闭合呈衰减至稳定的趋势。

（2）黏塑性阶段考虑裂隙连通程度影响有效渗流通道形成，引入修正系数，模型能够描述蠕变进入黏塑性阶段缓慢上升的趋势及加速蠕变阶段渗透率突增趋势。

（3）对各向异性蠕变-渗透率模型参数进行了敏感性分析。黏弹性阶段黏滞系数比增加，渗透率衰减越快，渗透率稳定值越低；侧向影响系数越高，渗透率初始衰减速率越高，渗透率稳定值越低。黏塑性阶段修正系数越大，加速蠕变阶段渗透率突增趋势越剧烈。