基于图像数字技术的砂岩裂隙可视化渗流特性试验研究

刘杰; 唐洪宇; 杨渝南; 石谦; 李政; 黎照; 高进; 兰俊

doi:10.11779/CJGE202011007

基于图像数字技术的砂岩裂隙可视化渗流特性试验研究 English Version

刘杰^{1, 2,},
唐洪宇^{1, 2},
杨渝南^{1, 2, ,},
石谦^{1, 2},
李政^{1, 2},
黎照^{1, 2},
高进^{1, 2},
兰俊³

1.
三峡大学湖北省地质灾害防治工程技术中心,湖北　宜昌　443002
2.
三峡大学三峡库区地质灾害教育部重点试验室,湖北　宜昌　443002
3.
武汉天华华中建筑设计有限公司,湖北　武汉　430021

基金项目:

湖北省技术创新重点项目 2017ACA189

国家自然科学基金项目 51579138

湖北省自然科学基金杰出青年人才计划项目 2018CFA065

湖北长江三峡滑坡国家野外科学观测研究站开放基金项目 2018KTL08

湖北省自然科学基金一般面上项目 2020CFB584

详细信息

作者简介:
刘杰(1979—),男,教授,博士生导师,主要从事岩土工程和水工结构研究。E-mail：liujieea@126.com

通讯作者:
杨渝南, E-mail：yangyunan_11@163.com

中图分类号: TU411
计量
- 文章访问数: 322
- HTML全文浏览量: 44
- PDF下载量: 265
出版历程
- 收稿日期: 2020-03-08
- 网络出版日期: 2022-12-05
- 刊出日期: 2020-10-31

Experimental research on visible seepage of sandstone fissure using digital image-based method

LIU Jie^{1, 2,},
TANG Hong-yu^{1, 2},
YANG Yu-nan^{1, 2, ,},
SHI Qian^{1, 2},
LI Zheng^{1, 2},
LI Zhao^{1, 2},
GAO Jin^{1, 2},
LAN Jun³

1.
Engineering Technology Center for Geological Disaster Prevention of Hubei Province, Yichang 443002, China
2.
Key Laboratory of Geological Hazards in Three Gorges Reservoir Area, Ministry of Education, Yichang 443002, China
3.
Wuhan Tianhua Huazhong Architectural Design Co., Ltd., Wuhan 430021, China

摘要

摘要: 自主研发的可视化裂隙渗流试验装置,应用颜色示踪的图像数字化处理技术,在不同渗透压和法向应力分级加卸载耦合作用下,进行空间垂直角度的砂岩裂隙可视化渗流试验。通过毫秒级分帧技术捕捉裂隙渗流状态并对渗流面积色域分区,基于数字化图像技术进行二值化处理,通过图层二次叠加渲染,自识别渗流面积参数,建立渗透压与法向应力作用下的渗流面积扩散规律。研究裂隙渗流优势路径各断面的过流宽度变化特征,提出可视化渗流速度求取公式,指出渗流速度沿流径变化为骤增、骤降、均匀损耗3个阶段,建立最小断面宽度对应的峰值流速与渗透压和法向应力耦合作用下的函数模型。基于实测数据建立主渗流路径中法向应力、渗透压与雷诺数的幂函数关系,标定渗流状态的惯性作用区与黏性作用区的转变临界点,构建可视化裂隙渗流雷诺数预测模型,对裂隙渗流中真实渗流路径识别、流速矢量实时变化、流体状态判定等关键科学问题提出了新的研究理论及方法。
- 可视化渗流 /
- 图像数字化处理 /
- 渗流扩散面积 /
- 过流断面宽度 /
- 渗流速度 /
- 雷诺数
Abstract: The visual fracture seepage test devices independently developed apply the color tracer image digital processing technology to carry out the visual seepage tests on the sandstone fractures with vertical spatial angle under the coupling effects of loading and unloading with different osmotic pressure and normal stress gradations. By means of the millisecond frame division technology, the seepage state of the crack is captured and the color domain of the seepage area is partitioned. Based on the digital image technology, the binary processing is carried out. The parameters of the seepage area are identified by double overlay of layers, and the diffusion law of the seepage area under the action of osmotic pressure and normal stress is established. The characteristics of the change of the flow width of each section along the dominant seepage path of the fracture are studied, and the formula for calculating the visual seepage velocity is put forward. It is pointed out that the seepage velocity changes along the flow diameter into three stages: sudden increase, sudden decrease and uniform loss. A functional model for the peak velocity corresponding to the minimum section width is established under the coupling action of osmotic pressure and normal stress. Based on the measured data of the main seepage path, the power function for the stress normal, osmotic pressure and Reynolds number is established, the change point of seepage state in the inertia function area and viscous effect area is calibrated, and the prediction model for Reynolds number of visualized fissure flows is formulated. The new theories and research method for the key scientific issues such as identification of real seepage in fissure seepage path, real-time change of velocity vector and determination of fluid state are put forward.
- visualized seepage /
- digital image processing /
- seepage diffusion area /
- cross section width /
- seepage velocity /
- Reynolds number

HTML全文

0. 引言

在岩土工程建设与运营中，地基中土单元往往发生主应力轴循环旋转，比如离岸岩土结构受波浪荷载^[1]，路基受交通荷载作用等^[2]。大量试验表明，主应力轴循环旋转会使产生极其复杂的塑性效应，一方面，导致土体产生累积塑性变形，排水条件下产生累积塑性压缩体变^[3]，不排水条件下产生累积孔压^[4-5]；另一方面，主应力轴循环旋转将诱发显著的变形刚度弱化与非共轴性^[5-8]，即材料当前的主应力方向与塑性主应变率方向不一致。主应力轴循环旋转下的复杂塑性效应严重威胁着离岸岩土结构、路基等的安全建设与运营，因此，合理地考虑主应力轴循环旋转塑性效应具有重要的理论与实际意义。

传统的弹塑性本构理论，绝大多数模型都隐含了共轴性(应变率与应力的主方向一致性)的假设，无法模拟主应力循环旋转诱发的非共轴性。为模拟非共轴性，一些模型引入了非共轴流动法则^[9-10]。但由于纯主应力轴循环旋转加载中主应力大小保持不变，在应力不变量表示的空间中不会产生屈服，因此，这类模型无法模拟纯主应力轴循环旋转诱发的塑性变形和非共轴性。为模拟主应力轴循环旋转塑性效应，需引入额外的力学机制和本构框架。Wang等^[11]在边界面亚塑性框架基础上构建了砂土本构模型，模型中塑性应变率方向与应力率方向相关，模拟了砂土在纯主应力轴循环旋转下的不排水有效应力路径。基于下加载面模型，Tsutsumi等^[12]引入了能产生非共轴性的各向异性和切向加载效应。Li等^[13]建议了一种各向异性临界状态为基础的本构理论框架，模型引入了与组构各向异性及非共轴流动相关的非比例加载机制。基于这种本构框架，Gao等^[14]通过引入组构及其演化规律，建立了三维本构模型模拟砂土主应力轴循环旋转下的非共轴性，童朝霞等^[15]则将主应力幅值变化以及应力主轴旋转产生的塑性变形单独加以考虑，建立了可考虑应力主轴循环旋转效应的砂土本构模型。而Lashkari等 ^[16]在Gutierrez等^[17]总结的应力空间内建立了相应的边界面，并引入基于主应力旋转试验规律的非共轴流动方向，通过考虑土体组构和非共轴流动方向，模拟了砂土的非共轴变形特征。Yang等^[18]在一般应力空间建立了一种包含切向加载的考虑主应力轴旋转的运动硬化模型。随后，Tian等^[19]从材料各向异性的角度将UH模型推广至模拟土的主应力轴循环旋转效应，通过引入组构及其演化规律考虑了固有和应力诱发各向异性，模拟了土的纯主应力轴循环旋转效应。另外，陈洲泉等^[20]在Qian等^[21]的基础上，重新定义了非共轴流动方法，并考虑了边界面塑性映射法则和组构各向异性，模拟了砂土主应力轴旋转的非共轴性。应该说，上述研究采用不同的力学机制实现了土的共轴变形特性的本构模拟，但研究对象仍局限于砂性土。相比较，最新试验观察发现^[5-8]主应力偏转路径下，天然黏性土的非共轴特性更为复杂，如纯主应力轴循环旋转下孔压循环波动的累积效应、剪切刚度的循环弱化效应等，并且这些非比例加载变形特性与非共轴性有一定关系^[6-7]，因此，为更准确地实现对这些复杂的力学响应的模拟，有必要对既有非共轴本构模型加以改进与完善。

基于以上认识，本文建立了能够描述主应力轴循环旋转塑性效应的边界面模型。模型通过引入可移动映射法则，描述主应力轴循环旋转卸载情况下的塑性变形。通过考虑固有各向异性弹性，描述天然软黏土主应力轴循环旋转下循环波动的塑性累积行为。同时，修正了非共轴流动法则考虑循环过程中非共轴性变化。最后，通过对比温州软黏土的模型预测和试验结果，验证了模型的有效性。

1. 边界面模型的建立

根据弹塑性力学的基本假定，总应变率 ${\dot \varepsilon _{ij}}$ 可分解为弹性部分和塑性部分，即

${\dot \varepsilon _{ij}}{\text{ = }}\dot \varepsilon _{ij}^{\text{e}} + \dot \varepsilon _{ij}^{\text{p}} \text{，}$

(1)

式中，弹性应变 $\dot \varepsilon _{ij}^{\text{e}}$ 可由弹性体积模量K和弹性剪切模量G来计算。根据剑桥模型的假定，K和G由当前的球应力p确定：

$K = \frac{p}{{{\kappa ^ * }}} \text{，} G = \frac{{3(1 - 2\nu )}}{{2(1 + \nu )}}K \text{，}$

(2)

式中， $\nu$ 为泊松比， ${\kappa ^*} = \kappa /(1 + {e_0})$ ， ${e_0}$ 为初始孔隙比， $\kappa$ 为e–lnp空间中回弹曲线的斜率。增量型弹性应力应变关系为

$\dot{\sigma}_{i j}=\boldsymbol{D}_{i j k l}^{\mathrm{e}} \dot{\varepsilon}_{k l}^{\mathrm{e}} \text{，}$

(3)

式中， $\boldsymbol{D}_{i j k l}^{\mathrm{e}}$ 为弹性刚度矩阵，可由K和G计算。

1.1 各向异性边界面方程

Ling等^[22]建议的单面边界面模型能够合理描述软黏土各向异性应力应变特性^[23-24]，如图 1所示，本文仍沿用这一方法来定义各向异性边界面方程。

图 1 各向异性边界面

Figure 1. Anisotropic bounding surface

下载: 全尺寸图片幻灯片

考虑各向异性的边界面方程如下：

$F = (\bar p - {p_{\text{c}}})\left( {\bar p + \frac{{R - 2}}{R}{p_{\text{c}}}} \right) + \frac{{{{(R - 1)}^2}}}{\chi }\bar q_\alpha ^2 = 0 \text{，}$

(4)

式中， $\bar p$ 和 ${\bar q_\alpha }$ 分别为平均应力和折减广义剪应力， ${\bar q_\alpha } = \sqrt {3\bar s_{ij}^\alpha \bar s_{ij}^\alpha /2}$ ，其中 $\bar s_{ij}^\alpha = {\bar s_{ij}} - {\bar \sigma _{kk}}{\alpha _{ij}}/3$ 为折减应力偏量； ${\alpha _{ij}}$ 为各向异性张量，表示土体各向异性的大小， $\alpha = \sqrt {3{\alpha _{ij}}{\alpha _{ij}}/2}$ 为p–q空间中屈服面的倾角（如图 1所示），M为临界状态应力比，R为形状参数，参数的定义及其物理意义参见黄茂松等^[23]和Huang等^[24]。

1.2 可移动映射法则

本文关注的纯主应力轴循环旋转加载，又称纯环剪，是非比例加载的一种，加载中主应力方向连续循环旋转，但广义剪应力q和主应力大小保持不变^[25]。由于主应力保持不变，因此，在π平面中表示为边界面上一个位置不动的点。为更好地表示这种纯主应力轴旋转的情况，本文将其表征在偏平面 ${\tau _{z\theta }}$ – ${\text{(}}{\sigma _z} -$ ${\sigma _\theta }{\text{)/2}}$ 上，其中 ${\sigma _z}$ ， ${\sigma _\theta }$ 和 ${\tau _{z\theta }}$ 分别为轴向应力、环向应力和剪切应力。在偏平面上，屈服面表示为绕坐标原点的圆，而纯主应力轴旋转应力路径如所示。加载过程中，应力增量方向与屈服面外法线方向夹角始终大于90°，即 ${\text{(}}\partial f/\partial {\sigma _{ij}}{\text{)d}}{\sigma _{ij}}$ < 0，无塑性变形产生，为卸载情况。值得一提的是，这种加载因子小于0是普遍存在的，如超固结土在循环荷载条件下的变形^[26-27]，为描述这种卸载情况下塑性变形的产生，本文借鉴Wang等^[11]和Li^[28]提出的可移动映射法则，其在π平面上的表示见Li^[28]，为更好地表示这种映射法则对纯主应力轴旋转的情况，本文将其表征在偏平面上。

图 2 偏平面τ_zθ–(σ_z-σ_θ)/2上纯主应力旋转应力路径

Figure 2. Stress path of pure principal stress rotation in τ_zθ–(σ_z-σ_θ)/2 space

下载: 全尺寸图片幻灯片

根据边界面塑性理论^{[24, 29]}可知，塑性加载因子定义为

$\dot{L}=\frac{1}{H_{\mathrm{p}}} \frac{\partial F}{\partial \bar{\sigma}_{i j}} \dot{\sigma}_{i j}=\frac{1}{\bar{H}_{\mathrm{p}}} \frac{\partial F}{\partial \bar{\sigma}_{i j}} \dot{\bar{\sigma}}_{i j} \text{，}$

(5)

式中 ${{\partial F} \mathord{\left/ {\vphantom {{\partial F} {\partial {{\bar \sigma }_{ij}}}}} \right. } {\partial {{\bar \sigma }_{ij}}}}$ 表示塑性加载方向，为边界面虚应力处外法线方向； ${H_{\text{p}}}$ 与 ${\bar H_{\text{p}}}$ 分别为当前应力点与虚应力点塑性模量，当前应力点与虚应力点重合时(即 ${\bar \sigma _{ij}} = {\sigma _{ij}}$ )， ${H_{\text{p}}} = {\bar H_{\text{p}}}$ 。

定义加载方向从加载转变为卸载（塑性加载因子 $\dot L$ 从正转变为负）时的上一步的应力点作为可移动的映射中心 ${\omega _{ij}}$ ，以连接映射中心 ${\omega _{ij}}$ 和当前应力点 ${\sigma _{ij}}$ 的直线与边界面的交点作为虚应力点 ${\bar \sigma _{ij}}$ ，如图 3所示，则构成了可移动映射法则。映射中心初始位置在偏平面的原点O。中所示为一般加载情况，塑性加载因子 $\dot L$ ＞0，则映射中心 ${\omega _{ij}}$ 保持不变。而对于一般卸载情况，如所示，塑性加载因子 $\dot L$ ＜0，此时，映射中心 ${\omega _{ij}}$ 移动至上一步的应力点位置，如所示，映射中心 ${\omega _{ij}}$ 移动后，以连接映射中心 ${\omega _{ij}}$ 和应力增量 ${\text{d}}{\sigma _{ij}}$ 的直线与边界面的交点作为虚应力点 ${\bar \sigma _{ij}}$ ，如中分析，塑性加载因子重新变为 $\dot L > 0$ 。即卸载时加载因子为负会导致映射中心 ${\omega _{ij}}$ 移动，而通过映射中心 ${\omega _{ij}}$ 移动会重新导致加载因子为正，这种可移动映射法则保证加卸载时有塑性变形产生。

图 3 偏平面τ_zθ–(σ_z-σ_θ)/2上加载映射法则

Figure 3. Mapping rule in τ_zθ–(σ_z-σ_θ)/2 space during general loading situation

下载: 全尺寸图片幻灯片

图 4 映射中心移动机制

Figure 4. Relocatable mechanism of projection center

下载: 全尺寸图片幻灯片

，中 $\rho$ 为映射中心 ${\omega _{ij}}$ 到当前应力点 ${\sigma _{ij}}$ 的距离， $\bar \rho$ 为映射中心 ${\omega _{ij}}$ 到虚应力点 ${\bar \sigma _{ij}}$ 的距离，映射法则表示为

${\bar \sigma _{ij}} = {\omega _{ij}} + \frac{{\bar \rho }}{\rho }({\sigma _{ij}} - {\omega _{ij}}) 。$

(6)

1.3 各向异性弹性的引入

由于天然软黏土的沉积使其呈现出横观各向同性弹性，又称各向异性弹性^[30-31]。主应力轴循环旋转下也会产生这种各向异性弹性效应，纯主应力轴循环旋转试验表明，不排水条件下的孔隙水压和排水条件下的体应变均呈波动形式累积^{[3-6, 8]}，如图 5所示。根据体应变可分为不可逆和可逆两种的定义^[32]，不排水试验中孔隙水压波动的下限对应于排水试验中永久累积塑性不可逆体应变，而上限和下限的差即为弹性波动可逆的孔隙水压，如图 5所示。早期的研究表明，孔隙水压波动这种塑性累积行为可能归因于天然软黏土弹性行为的各向异性^[30-31]。

图 5 纯主应力轴循环旋转孔隙水压累积中呈现的各向异性弹性

Figure 5. Anisotropic elasticity in accumulation of excess pore pressure with cycles of pure principal stress rotation

下载: 全尺寸图片幻灯片

为考虑这种波动塑性累积行为，模型中引入各向异性弹性，可表示为

${\left[ D \right]^{ - 1}} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{E_{\text{v}}}}}}&{ - \frac{{{\nu _{{\text{vv}}}}}}{{{E_{\text{v}}}}}}&{ - \frac{{{\nu _{{\text{vv}}}}}}{{{E_{\text{v}}}}}}&0&0&0 \\ { - \frac{{{\nu _{{\text{vv}}}}}}{{{E_{\text{v}}}}}}&{\frac{1}{{{E_{\text{h}}}}}}&{ - \frac{{{\nu _{{\text{vh}}}}}}{{{E_{\text{h}}}}}}&0&0&0 \\ { - \frac{{{\nu _{{\text{vv}}}}}}{{{E_{\text{v}}}}}}&{ - \frac{{{\nu _{{\text{vh}}}}}}{{{E_{\text{h}}}}}}&{\frac{1}{{{E_{\text{h}}}}}}&0&0&0 \\ 0&0&0&{\frac{1}{{2{G_{{\text{vh}}}}}}}&0&0 \\ 0&0&0&0&{\frac{1}{{2{G_{{\text{vh}}}}}}}&0 \\ 0&0&0&0&0&{\frac{{1 + {\nu _{{\text{vh}}}}}}{{{E_{\text{h}}}}}} \end{array}} \right] \text{，}$

(7)

式中， ${E_{\text{h}}}$ 和 ${E_{\text{v}}}$ 分别为水平和竖向弹性模量， ${\nu _{{\text{vh}}}}$ 和 ${\nu _{{\text{vv}}}}$ 分别为水平和竖向泊松比， ${G_{{\text{vh}}}}$ 为剪切模量，满足 ${E_{\text{h}}} = {n^2}{E_{\text{v}}}$ ， ${\nu _{{\text{vh}}}} = n{\nu _{{\text{vv}}}}$ ， ${G_{{\text{vh}}}} = {{n{E_{\text{v}}}} \mathord{\left/ {\vphantom {{n{E_{\text{v}}}} {\left[ {2(1 + n{\nu _{{\text{vv}}}})} \right]}}} \right. } {\left[ {2(1 + n{\nu _{{\text{vv}}}})} \right]}}$ ，其中， $n$ 为各向异性因子。

对于应力控制的各向同性压缩，根据上述方程可知：

${\dot \varepsilon _{\text{v}}} = {\dot \varepsilon _{11}} + {\dot \varepsilon _{22}} + {\dot \varepsilon _{33}} = \left( {1 - 4{\nu _{{\text{vv}}}} + \frac{2}{{{n^2}}} - \frac{{2{\nu _{{\text{vv}}}}}}{n}} \right)\frac{p}{{{E_{\text{v}}}}}。$

(8)

考虑到体积模量K可根据式(2)计算得到，竖向弹性模量可表示为

${E_{\text{v}}} = \left( {1 - 4{\nu _{{\text{vv}}}} + \frac{2}{{{n^2}}} - \frac{{2{\nu _{{\text{vv}}}}}}{n}} \right)\frac{{1 + {e_0}}}{\kappa }p。$

(9)

不同于各向同性弹性模型回弹指数 $\kappa$ 、泊松比 $\nu$ ，在各向异性弹性中还需要一个各向异性因子 $n$ 。根据以上分析可知， ${n^2}$ 为水平与竖向刚度之比，当 $n$ =1时，材料为各向同性。

1.4 非共轴流动法则

Qian等^[21]将非共轴应力率在当前应力状态方向上进行正交投影，给出了广义应力空间内的非共轴应力率：

$\dot s_{ij}^n = {\dot s_{ij}} - \frac{{{{\dot s}_{kl}}{s_{kl}}}}{{{s_{mn}}{s_{mn}}}}{s_{ij}} - \frac{{{{\dot s}_{kl}}{S_{kl}}}}{{{S_{mn}}{S_{mn}}}}{S_{ij}} \text{，}$

(10)

式中， ${S_{ij}} = {s_{ik}}{s_{kj}} - {{2{J_2}{\delta _{ij}}} \mathord{\left/ {\vphantom {{2{J_2}{\delta _{ij}}} 3}} \right. } 3} - {{3{s_{ij}}{J_3}} \mathord{\left/ {\vphantom {{3{s_{ij}}{J_3}} {2{J_2}}}} \right. } {2{J_2}}}$ 。相应的非共轴应变率定义为

$\dot \varepsilon _{ij}^{pn} = \frac{1}{{{H_{\text{t}}}}}\dot s_{ij}^n \text{，}$

(11)

式中， ${H_{\text{t}}}$ 为塑性模量。

大量黏土试验结果表明^{[8, 33-34]}，随着当前剪切应力增大，非共轴应变率与应力方向不一致逐渐减小，至临界状态时趋于共轴，这种非共轴特性不能被方程（11）定义的非共轴特性反映^{[16, 20]}。另外，方程（11）中的非共轴参数为切向塑性模量，在循环剪切过程中，切向塑性模量呈现显著的弱化效应，若忽略其弱化效应则无法准确地预测非共轴塑性变形。因此，为更合理地模拟非共轴流动性，本文修正了Qian等^[21]非共轴流动理论，具体表述为

$\dot \varepsilon _{ij}^{pn} = \dot L{\left\langle {1 - \frac{\eta }{M}} \right\rangle ^\gamma }{k_{\text{n}}}n_{ij}^{nc}\text{，}$

(12)

式中， $n_{ij}^{nc}$ 为和 $\dot s_{ij}^n$ 具有相同方向的非共轴应力率，表述为

$n_{ij}^{nc}{\text{ = }}n_{ij}^{\dot s} - \frac{{n_{kl}^{\dot s}{s_{kl}}}}{{{s_{mn}}{s_{mn}}}}{s_{ij}} - \frac{{n_{kl}^{\dot s}{S_{kl}}}}{{{S_{mn}}{S_{mn}}}}{S_{ij}} \text{，}$

(13)

其中， $n_{ij}^{\dot s} = {{{{\dot s}_{ij}}} \mathord{\left/ {\vphantom {{{{\dot s}_{ij}}} {\left\| {{{\dot s}_{ij}}} \right\|}}} \right. } {\left\| {{{\dot s}_{ij}}} \right\|}}$ 为单位偏应力增量， $\eta$ 和 $M$ 分别为当前应力比和临界状态应力比，当 $\eta$ = $M$ 时， $\dot \varepsilon _{ij}^{np} = 0$ ， $\gamma$ ， ${k_{\text{n}}}$ 均为材料参数，控制非共轴塑性应变随当前应力比变化程度及量级， $\dot L$ 为塑性加载因子。

Macaulay括号内表示随着应力比增加，非共轴程度减小，至临界状态时，应变增量方向与应力方向共轴，这与试验结论是符合的。值得一提的是，与Qian等^[21]非共轴不同的是，方程（13）中非共轴应变率与共轴应变率用了相同的塑性加载因子，这是由于非共轴与共轴应变率均为不可恢复塑性变形，由相同的塑性机制驱动，与Li等^[35]讨论的一致，仅仅是方向不同。另外，虽然用相同的塑性加载因子，但非共轴与共轴应变率强度却不同， $\gamma$ 和 ${k_{\text{n}}}$ 共同控制非共轴塑性应变相对于共轴塑性应变的量级。通过以上分析可知，本文修正使非共轴应力率无量纲化，然后非共轴流动可以和共轴流动通过塑性加载因子耦合起来，并考虑应力比影响。

共轴塑性应变采用相关联流动法则，即塑性流动方向与塑性加载方向重合，且等于边界面在虚应力点处的外法线方向。考虑共轴塑性应变和非共轴塑性应变，可得总的流动法则：

$\dot \varepsilon _{ij}^{\text{p}} = \dot L\left[ {\frac{{\partial F}}{{\partial {{\bar \sigma }_{ij}}}} + {{\left\langle {1 - \frac{\eta }{M}} \right\rangle }^\gamma }{k_{\text{n}}}n_{ij}^{nc}} \right]。$

(14)

1.5 增量型应力应变关系

模型采用两个独立的硬化内变量 ${p_{\text{c}}}$ 以及 ${\alpha _{ij}}$ , 其中 ${p_{\text{c}}}$ 反映应力历史的影响， ${\alpha _{ij}}$ 定义屈服面倾角的大小，反映各向异性的程度。硬化规律的具体计算见Huang等^[24]。边界面虚应力点的塑性模量 ${\bar H_{\text{p}}}$ 可根据一致性条件求出：

${\bar H_{\text{P}}} = - \frac{{\partial F}}{{\partial {p_{\text{c}}}}}\frac{{\partial {p_{\text{c}}}}}{{\partial \varepsilon _{\text{v}}^{\text{p}}}}\frac{{\partial F}}{{\partial \bar p}} - \frac{{\partial F}}{{\partial {\alpha _{ij}}}}\left( {\frac{{\partial {\alpha _{ij}}}}{{\partial \varepsilon _{\text{v}}^{\text{p}}}}\frac{{\partial F}}{{\partial \bar p}} + \frac{{\partial {\alpha _{ij}}}}{{\partial \varepsilon _{\text{s}}^{\text{p}}}}\frac{{\partial F}}{{\partial \bar q}}} \right) 。$

(15)

当前应力点的塑性模量 ${H_{\text{p}}}$ 可由插值函数计算^[24]：

${H_{\text{P}}} = {\bar H_{\text{P}}} + \zeta {p_{\text{a}}}\left[ {{{\left( {\frac{{\partial F}}{{\partial \bar p}}} \right)}^2} + {{\left( {\frac{{\partial F}}{{\partial \bar q}}} \right)}^2}} \right]\left[ {{{\left( {\frac{{\bar \rho }}{\rho }} \right)}^\psi } - 1} \right] \text{，}$

(16)

$\psi = {\psi _0}\exp ( - \xi \varepsilon _{\text{s}}^{\text{p}}) 。$

(17)

式中， ${p_{\text{a}}}$ 为标准大气压（101.325 kPa）， $\varepsilon _{\text{s}}^{\text{p}} = \int {\dot \varepsilon _{\text{s}}^{\text{p}}}$ 为塑性偏应变的累积值， $\zeta$ ， ${\psi _0}$ ， $\xi$ 为模型参数。

增量型弹性应力应变关系见式（3）， $\boldsymbol{D}_{ijkl}^{\text{e}}$ 为弹性刚度矩阵，可由式（7）计算，将式（1），（3），（5），（14）代入一致性条件得

$\dot L = \frac{{\frac{{\partial F}}{{\partial {{\bar \sigma }_{st}}}}\boldsymbol{D}_{stkl}^{\text{e}}{{\dot \varepsilon }_{kl}}}}{{{H_{\text{P}}} + \frac{{\partial F}}{{\partial {{\bar \sigma }_{st}}}}\boldsymbol{D}_{stkl}^{\text{e}}\left[ {\frac{{\partial F}}{{\partial {{\bar \sigma }_{kl}}}} + {{\left\langle {1 - \frac{\eta }{M}} \right\rangle }^\gamma }{k_{\text{n}}}n_{kl}^{nc}} \right]}}。$

(18)

最终可得如下弹塑性应力应变增量关系：

${\dot \sigma _{ij}} = \boldsymbol{D}_{ijkl}^{{\text{ep}}}{\dot \varepsilon _{kl}} \text{，}$

(19)

其中， $\boldsymbol{D}_{ijkl}^{{\text{ep}}}$ 为弹塑性刚度矩阵，表达为

${\boldsymbol{D}}_{ijkl}^{\text{ep}}={D}_{ijkl}^{\text{e}}-\frac{{D}_{ijmn}^{\text{e}}\left[\frac{\partial F}{\partial {\overline{\sigma }}_{mn}}+{\langle 1-\frac{\eta }{M}\rangle }^{\gamma }{k}_{\text{n}}{n}_{mn}^{nc}\right]\frac{\partial F}{\partial {\overline{\sigma }}_{st}}{D}_{stkl}^{\text{e}}}{{H}_{\text{P}}+\frac{\partial F}{\partial {\overline{\sigma }}_{st}}{D}_{stkl}^{\text{e}}\left[\frac{\partial F}{\partial {\overline{\sigma }}_{kl}}+{\langle 1-\frac{\eta }{M}\rangle }^{\gamma }{k}_{\text{n}}{n}_{kl}^{nc}\right]}。$

(20)

本文所有计算都通过Bardet等^[36]提出的积分算法实现。

1.6 模型参数

本文模型共需要14个参数和初始状态参数 ${e_0}$ ，如表 1所示，模型相关参数根据其功能可分为6组。相对于黄茂松等^[23]和Huang等^[24]采用的边界面模型，模型中减少了结构性参数，引入了非共轴参数 $\gamma$ 、 ${k_{\text{n}}}$ 和各向异性弹性参数 $n$ ，除上述3个参数以外的其他参数，其标定方法可以参考黄茂松等^[23]和Huang等^[24]的论述，本文不再重复。非共轴参数 $\gamma$ 控制随着应力比增加非共轴程度减小的速率， ${k_{\text{n}}}$ 控制非共轴塑性应变相对于共轴塑性应变的大小。这两个非共轴特性参数 $\gamma$ 和 ${k_{\text{n}}}$ 可由两组 ${\alpha _\sigma }$ ≠0°和 ${\alpha _\sigma }$ ≠90°的主应力轴固定单调剪切试验来标定，根据式（12）可得

$\gamma = {{\ln \left( {\frac{{{\text{d}}\varepsilon _{z\theta }^{{\text{p1}}}}}{{{\text{d}}\varepsilon _{z\theta }^{{\text{p2}}}}}} \right)} \mathord{\left/ {\vphantom {{\ln \left( {\frac{{{\text{d}}\varepsilon _{z\theta }^{{\text{p1}}}}}{{{\text{d}}\varepsilon _{z\theta }^{{\text{p2}}}}}} \right)} {\ln \left( {\frac{{1 - {{{\eta _1}} \mathord{\left/ {\vphantom {{{\eta _1}} M}} \right. } M}}}{{1 - {{{\eta _2}} \mathord{\left/ {\vphantom {{{\eta _2}} M}} \right. } M}}}} \right)}}} \right. } {\ln \left( {\frac{{1 - {{{\eta _1}} \mathord{\left/ {\vphantom {{{\eta _1}} M}} \right. } M}}}{{1 - {{{\eta _2}} \mathord{\left/ {\vphantom {{{\eta _2}} M}} \right. } M}}}} \right)}} \text{，}$

(21)

表 1 模型参数

Table 1. Model parameters

临界状态参数	形状参数	硬化参数	非共轴参数	各向异性弹性参数	边界面插值参数
$\lambda$ =0.36 $\kappa$ =0.04 ${M_{\text{e}}}$ =0.81 ${M_{\text{c}}}$ =1.11 $\nu$ =0.2	R=3.2	$\mu$ =60 $\beta$ =0.7	$\gamma$ =0.5 ${k_{\text{n}}}$ =50	n=1.1	$\zeta$ =40.0 ${\psi _0}$ =12.0 $\xi$ =6.0

下载: 导出CSV

| 显示表格

式中，( ${\text{d}}\varepsilon _{z\theta }^{{\text{p1}}}$ , ${\eta _1}$ )和( ${\text{d}}\varepsilon _{z\theta }^{{\text{p2}}}$ , ${\eta _2}$ )为 ${\text{d}}\varepsilon _{z\theta }^{\text{p}}$ – $\eta$ 曲线上两个不同的点，参数 ${k_{\text{n}}}$ 可在 $\gamma$ 确定后通过拟合其中一条 $\varepsilon _{z\theta }^{\text{p}}$ – $\eta$ 曲线来标定。各向异性弹性参数 $n$ 为各向异性因子，是天然黏土的各向异性弹性参数。 ${n^2}$ 为水平与竖向刚度之比，可从固结后水平与竖向切片三轴试验或弯曲元试验获得的弹性模量或泊松比之比来标定。如果试验数据不足，也可通过试错法模拟可用数据，然后通过最佳拟合来确定，特别是通过模拟不排水试验孔压波动曲线或排水试验中体应变波动曲线来标定。

2. 模型验证

模型通过温州软黏土在纯主应力轴循环旋转下的不排水空心圆柱扭剪试验进行验证。Wang等^[8]对温州天然软黏土开展了一系列纯主应力轴循环旋转试验，首先对试样进行各向同性固结，初始固结压力p′=150 kPa，然后增加剪应力q至固定值，保持各向异性固结，随后进行纯主应力轴循环旋转不排水试验，如2.2节中描述，保持平均主应力p、中主应力系数b和剪应力q不变，连续旋转主应力方向，试验具体条件和方案见Wang等^[8]，试验初始条件为 ${p_0}$ =100 kPa和 ${e_0}$ =1.570。土体计算参数如所示，需要说明的是，该试验并未提供临界状态参数 $\lambda$ 和 $\kappa$ ，根据Nakase等^[37]的研究，临界状态参数 $\lambda$ 和 $\kappa$ 与土体物理性质指标液限w_L或塑性指数I_P存在线性关系，但其关系需根据土的特性来确定，参考王立忠等^[38]对温州黏土的研究，对比其土体物理性质指标液限w_L和塑性指数I_P，综合确定的温州软黏土临界状态参数见表 1。

纯主应力轴循环旋转下，温州饱和软黏土不排水试验的孔隙水压变化与模拟结果如图 6所示。由图 6可知，虽然每个循环广义剪应力q保持不变，但孔压随主应力轴旋转呈波动形式累积，本文本构模型能够合理模拟孔压波动这种塑性积累行为。由图 6还可看出，b=0的孔压累积小于b=0.5，1.0，在循环的初级阶段，孔压累积较快，随循环数增加，孔压累积变慢，而本文模型也能合理模拟中主应力系数和循环阶段对孔压累积的影响。图 6中虚线表示b=1.0条件下采用各向同性弹性n=1的模拟，即不考虑各向异性弹性，结果显示没有孔压波动。通过对比可知，孔压波动这种塑性累积行为的模拟要归结于本构模型中引入的各向异性弹性。

图 6 孔隙水压累积试验与模拟

Figure 6. Tests and simulation of accumulation of pore water pressure

下载: 全尺寸图片幻灯片

图 7为温州软黏土在纯主应力轴循环旋转下b=1.0系列的应变分量随主应力轴旋转角度的关系。由图 7可知，尽管广义剪应力保持不变，纯主应力轴旋转仍然产生了显著的塑性应变累积，模型合理的模拟这种塑性应变累积，这主要由于边界面模型中采用了可移动映射法则。

图 7 应变分量变化规律试验与模拟

Figure 7. Tests and simulation of variation of strain components

下载: 全尺寸图片幻灯片

纯主应力轴循环旋转b=1.0下，温州饱和软黏土不排水试验的剪切应力应变变化与模拟结果如图 8所示，图 8（a），8（b）分别为5个循环的试验和模拟结果。如图 8（a）所示，尽管广义剪应力保持不变，剪切应力应变曲线仍呈现滞回特性并伴随着累积塑性的发展，滞回圈刚度显著弱化且滞回圈形状呈开放型性发展，对比图 8（a），8（b）可知，尽管有一定差距，但模型合理地模拟了试验中这种剪切应力应变关系。

图 8 剪切应力应变关系试验与模拟

Figure 8. Tests and simulation of shear stress strain relationship

下载: 全尺寸图片幻灯片

图 9为温州软黏土在纯主应力轴循环旋转下b=0.5的第三个循环（N=3）应变增量矢量试验和模拟结果，相应的非共轴角随主应力轴旋转角变化的试验和模拟如所示。由可知，应变增量在循环加载过程中具有明显的分段特征，当 ${\alpha _\sigma }$ 在[0°，45°]和[90°，135°]变化时，非共轴角度减小，而当主应力轴旋转角 ${\alpha _\sigma }$ 在[45°，90°]和[135°，180°]变化时，非共轴角增大，由模拟结果可知，模型计算出的分段特征大致与试验数据吻合。由的对比结果可知，模型计算的非共轴角随循环变化规律也与试验结果一致，说明模型能合理模拟黏土随主应力轴循环旋转表现的非共轴性。还给出了不考虑非共轴流动法则 ${k_{\text{n}}}$ =0，即共轴模型的预测结果，由于采用的可移动映射法则的映射中心并不在坐标中心，因此塑性流动方向并不等于应力方向角 ${\alpha _\sigma }$ ，所以非共轴角并不等于0，但是几乎不随主应力轴旋转变化。因此，对比可知本文引入的修正非共轴流动法则在非共轴性预测中起显著作用，合理模拟了非共轴角随循环变化的规律。

图 9 非共轴特性试验与模拟(N=3)

Figure 9. Tests and simulation of noncoaxial behavior (N=3)

下载: 全尺寸图片幻灯片

图 10 非共轴角试验与模拟

Figure 10. Tests and simulation of noncoaxial angle

下载: 全尺寸图片幻灯片

3. 结论

本文在考虑各向异性边界面模型的基础上，通过引入可移动映射法则，处理主应力轴循环旋转卸载下的塑性效应。通过考虑固有各向异性弹性，模拟主应力轴循环旋转下循环波动的塑性累积行为。同时，修正了非共轴流动法则反映非共轴性的循环变化，建立了一种考虑主应力轴循环旋转塑性效应的本构模型。最后，对温州天然饱和软黏土在纯主应力轴循环旋转下的不排水行为进行了模拟，得到以下3点结论。

（1）在纯主应力轴循环旋转试验中，因为模型引入了可移动映射法则，使模型可描述这种卸载情况下的塑性变形。

（2）由于模型中引入了固有各向异性弹性，从而模型可以描述纯主应力轴循环旋转中循环波动这种塑性累积行为。

（3）修正的非共轴理论中考虑了应力比影响并将非共轴流动和共轴流动耦合起来，从而可以描述纯主应力轴循环旋转中非共轴性随循环变化规律。

图 1 预制V形状槽劈裂法

Figure 1. Prefabricated V-shaped slot splitting method

下载: 全尺寸图片幻灯片

图 2 硅胶翻模

Figure 2. Silica gel turning mould

下载: 全尺寸图片幻灯片

图 3 可视化试样制备

Figure 3. Preparation of visualized sample

下载: 全尺寸图片幻灯片

图 4 试样组装

Figure 4. Assembly of sample

下载: 全尺寸图片幻灯片

图 5 裂隙渗流可视化试验装置图

Figure 5. Installation diagram of visualized fracture seepage tests

下载: 全尺寸图片幻灯片

图 6 裂隙可视化渗流试验流程图

Figure 6. Flow chart of visualized tests of fracture

下载: 全尺寸图片幻灯片

图 7 分帧、分区、特征点提取图像示例

Figure 7. Example of image extraction by frame, partition and feature points

下载: 全尺寸图片幻灯片

图 8 渗流面积临界Ⅰ值求取示意图

Figure 8. Schematic diagram of critical Ⅰ value of contact area of fracture surface

下载: 全尺寸图片幻灯片

图 9 渗流面积Ⅰ值划分

Figure 9. Division of seepage area byⅠvalue

下载: 全尺寸图片幻灯片

图 10 识别盲区与误区

Figure 10. Identification blind and misunderstanding area of seepage area

下载: 全尺寸图片幻灯片

图 11 自识别正确

Figure 11. Correctness of self identification

下载: 全尺寸图片幻灯片

图 12 图层混合选项处理

Figure 12. Processing of layer blending options

下载: 全尺寸图片幻灯片

图 13 裂隙面图形总面积像素值

Figure 13. Pixel values of total area of fracture surface

下载: 全尺寸图片幻灯片

图 14 裂隙面图形非渗流面积像素值

Figure 14. Pixel values of non-seepage of fracture surface

下载: 全尺寸图片幻灯片

图 15 渗流面积与法向应力、渗透压的关系

Figure 15. Relationship among seepage area and normal stress and osmotic pressure

下载: 全尺寸图片幻灯片

图 16 k与法向应力的关系

Figure 16. Relationship betweenk and normal stress

下载: 全尺寸图片幻灯片

图 17 图像断面划分与宽度获取

Figure 17. Division of image section and width acquisition

下载: 全尺寸图片幻灯片

图 18 过流断面宽度变化

Figure 18. Variation of cross-section width

下载: 全尺寸图片幻灯片

图 19 法向应力0.1 MPa时各断面速度关系

Figure 19. Relationship between speed section at normal stress of 0.1 MPa

下载: 全尺寸图片幻灯片

图 20 法向应力0.6 MPa时各断面速度关系图

Figure 20. Relationship between speed sections at normal stress of 0.6 MPa

下载: 全尺寸图片幻灯片

图 21 各分区速度变化率

Figure 21. Change rates of various speed partitions

下载: 全尺寸图片幻灯片

图 22 渗透压与峰值流速v_max的影响关系

Figure 22. Relationship between osmotic pressure and maximum velocity

下载: 全尺寸图片幻灯片

图 23 增长率a随法向应力的变化

Figure 23. Change of growth rate with normal stress

下载: 全尺寸图片幻灯片

图 24 0.1 MPa不同渗透压下渗流路径长度与雷诺数的关系

Figure 24. Reynolds number of flow path under different osmotic pressures at normal stress of 0.1 MPa

下载: 全尺寸图片幻灯片

图 25 m与渗透压的关系

Figure 25. Relationship betweenm and osmotic pressure

下载: 全尺寸图片幻灯片

图 26 μ和法向应力的关系

Figure 26. Relationship betweenμ and normal stress

下载: 全尺寸图片幻灯片

图 27 n与渗透压的关系

Figure 27. Relationship betweenn and osmotic pressure

下载: 全尺寸图片幻灯片

图 28 κ与法向应力的关系

Figure 28. Relationship betweenκ and normal stress

下载: 全尺寸图片幻灯片

表 1 渗流面积分类与RGB值对应表

Table 1 Classification of seepage area and values of RGB

区域命名	渗流面积分类	颜色表观	RGB数值大小反映颜色的偏向			临界	归一化
区域命名	渗流面积分类	颜色表观	R	G	B	Ⅰ值	Ⅰ值
1区域	有效渗流面积	液体扩散区	70~105	40~85	35~85	0.32	0.11~0.32
2区域	临界渗流面积	液体微过渡区	101~130	51~90	51~90	—	0.32~0.342
3区域	非渗流面积	液体未扩散区	131~160	90~120	90~120	0.342	0.342~0.918

下载: 导出CSV

表 2 黑白高亮二值化图

Table 2 Black and white highlighted binarization diagram

法向应力/MPa	渗透压/MPa
法向应力/MPa	0.1	0.2	0.3	0.4	0.5	0.6
0.1
0.6

下载: 导出CSV

表 3 裂隙固有渗流参数δ值

Table 3 Values of intrinsic seepage parameter of crack

法向应力/MPa	0.1	0.2	0.3	0.4	0.5	0.6
δ值	49.73	47.17	45.93	41.90	39.17	23.43

下载: 导出CSV

表 4 主惯性与主黏性临界点

Table 4 Critical points of main inertia and main viscosity

临界点坐标	渗透压/MPa
临界点坐标	0.1	0.2	0.3	0.4	0.5	0.6
x（X）	8.33	7.52	7.65	7.93	8.08	8.52
y（Re）	0.71	1.15	1.18	1.74	1.86	2.27

下载: 导出CSV

表 5 雷诺数预测值的平均相对误差

Table 5 Average relative errors of predicted values of Reynolds number

法向应力/MPa	0.1	0.2	0.3	0.4	0.5	0.6
平均相对误差ε	0.2659	0.1442	0.1811	0.1872	0.2059	0.1667

下载: 导出CSV

参考文献(18)

[1]	LOUIS C. Rock Hydraulics in Rock Mechanics[M]. New York: Springer-Verlag, 1974.
[2]	速宝玉, 詹美礼, 赵坚. 仿天然岩体裂隙渗流的试验研究[J]. 岩土工程学报, 1995, 17(5): 19-24. doi: 10.3321/j.issn:1000-4548.1995.05.004 SU Bao-yu, ZHAN Mei-li, ZHAO Jian. Study on fracture seepage in the imitative nature rock[J]. Chinese Journal of Geotechnical Engineering, 1995, 17(5): 19-24. (in Chinese) doi: 10.3321/j.issn:1000-4548.1995.05.004
[3]	许光祥, 张永兴, 哈秋舲. 粗糙裂隙渗流的超立方和次立方定律及其试验研究[J]. 水利学报, 2003, 34(3): 74-79. doi: 10.3321/j.issn:0559-9350.2003.03.014 XU Guang-xiang, ZHANG Yong-xing, HA Qiu-ling. Super- cubic and sub-cubic law of rough fracture seepage and its experimental study[J]. Journal of Hydraulic Engineering, 2003, 34(3): 74-79. (in Chinese) doi: 10.3321/j.issn:0559-9350.2003.03.014
[4]	SINGH K K, SINGH D N, RANJITH P G. Laboratory simulation of flow through single fractured granite[J]. Rock Mechanics and Rock Engineering, 2015, 48(3): 987-1000. doi: 10.1007/s00603-014-0630-9
[5]	BABADAGLI T, REN X, DEVELI K. Effects of fractal surface roughness and lithology on single and multiphase flow in a single fracture:an experimental investigation[J]. International Journal of Multiphase Flow, 2015, 68: 40-58. doi: 10.1016/j.ijmultiphaseflow.2014.10.004
[6]	BRUSH D J, THOMSON N R. Fluid flow in synthetic roughwalled fractures:Navier-Stokes, Stokes, and local cubic law simulations[J]. Water Resources Research, 2003, 39(4): 1085-1100.
[7]	OLSSON R, BARTON N. An improved model for hydro mechanical coupling during shearing of rock joints[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(3): 317-329. doi: 10.1016/S1365-1609(00)00079-4
[8]	胡昱, 源新, 刘光廷, 等. 多轴应力作用下砂砾岩单裂隙渗流规律试验研究[J]. 地下空间与工程学报, 2007, 3(6): 1009-1013. doi: 10.3969/j.issn.1673-0836.2007.06.006 HU Yu, YUAN Xin, LIU Guang-ting, et al. Experiment research on the laws of seepage in calcirudite rock within single fracture under multiaxial stresses[J]. Chinese Journal of Underground Space and Engineering, 2007, 3(6): 1009-1013. (in Chinese) doi: 10.3969/j.issn.1673-0836.2007.06.006
[9]	刘才华, 陈从新. 三轴应力作用下岩石单裂隙的渗流特性[J]. 自然科学进展, 2007, 17(7): 989-994. doi: 10.3321/j.issn:1002-008X.2007.07.022 LIU Cai-hua, CHEN Cong-xin. Seepage characteristics of single fracture under triaxial stress[J]. Advances in Natural Science, 2007, 17(7): 989-994. (in Chinese) doi: 10.3321/j.issn:1002-008X.2007.07.022
[10]	王志良, 申林方, 徐则民, 等. 岩体裂隙面粗糙度对其渗流特性的影响研究[J]. 岩土工程学报, 2016, 38(7): 1262-1268. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201607013.htm WANG Zhi-liang, SHEN Lin-fang, XU Ze-min, et al. Influence of roughness of rock fracture on seepage characteristics[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(7): 1262-1268. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201607013.htm
[11]	QIAN J, ZHAN H, LUO S, et al. Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture[J]. Journal of Hydrology, 2007, 339(3/4): 206-215.
[12]	RANJITH P G, DARLINGTON W. Nonlinear single-phase flow in real rock joints[J]. Water Resources Research, 2007, 43(9): 146-156.
[13]	ZHANG Z, NEMCIK J. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures[J]. Journal of Hydrology, 2013, 477(16): 139-151.
[14]	胡少华, 周佳庆, 陈益峰, 等. 岩石粗糙裂隙非线性渗流特性试验研究[J]. 地下空间与工程学报, 2017, 13(1): 48-56. https://www.cnki.com.cn/Article/CJFDTOTAL-BASE201701008.htm HU Shao-hua, ZHOU Jia-qing, CHEN Yi-feng, et al. Laboratory research on nonlinear flow behavior of rough fractures[J]. Chinese Journal of Underground Space and Engineering, 2017, 13(1): 48-56. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BASE201701008.htm
[15]	鞠杨, 谢和平, 郑泽民, 等. 基于3D打印技术的岩体复杂结构与应力场的可视化方法[J]. 科学通报, 2014, 59(32): 3109-3119. https://www.cnki.com.cn/Article/CJFDTOTAL-KXTB201432002.htm JU Yang, XIE He-ping, ZHENG Ze-ming, et al. Visualization of the complex structure and stress field inside rock by means of 3D printing technology[J]. Chinese Science Bull, 2014, 59(32): 3109-3119. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-KXTB201432002.htm
[16]	刘建军, 汪尧, 宋睿, 等. 基于透明岩土材料的可视化渗流实验及其应用前景[J]. 地球科学, 2017, 42(8): 1287-1295. https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201708004.htm LIU Jian-jun, WANG Yao, SONG Rui, et al. Visual seepage experiment based on transparent rock-soil material and its application prospect[J]. Earth Science, 2017, 42(8): 1287-1295. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201708004.htm
[17]	盛金昌, 刘继山, 赵坚. 基于图像数字化技术的裂隙岩体非稳态渗流分析[J]. 岩石力学与工程学报, 2006, 25(7): 1402-1407. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200607015.htm SHENG Jin-chang, LIU Ji-shan, ZHAO Jian. Analysis of transient fluid flow in fractured rock masses with digital image-based method[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(7): 1402-1407. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200607015.htm
[18]	FOXRW , MCDONALD A T, PRITCHARD P J. Introduction to Fluid Mechanics[M]. 6th ed. Hoboken: John Wiley and Sons, 2004: 348.