, 邵玉龙, 黄帆, 周创兵南昌大学 建筑工程学院, 南昌 330031
收稿日期：2021-01-20
基金项目：国家自然科学基金资助项目（U1765207，52079062，41762020）；江西省自然科学基金资助项目（20192ACB2102）；江西省研究生创新专项资金项目（YC2019-B009）
作者简介：何忱(1994-), 男, 博士研究生
通讯作者：姚池, 教授, E-mail: chi.yao@ncu.edu.cn
摘要：裂隙岩体的渗透性是许多地下工程的重要参数之一。低裂隙密度下岩体渗透性特征与目前广泛研究的高裂隙密度显著不同，并且具有强烈的不确定性。该文首先提出了一种裂隙多孔介质渗流数值方法，该方法建立在等效离散裂隙网络模型之上，采用二维三角形单元网络等效描述三维岩体的孔隙基质和裂隙，推导了等效单元导水系数的解析表达式。该方法自然考虑了裂隙与基质之间的流量交换，且无需添加额外的参数。随后，对不同裂隙密度、不同研究尺度的裂隙岩体渗透性进行了Monte-Carlo模拟。结果表明：裂隙网络是否贯通研究区域对有效渗透性影响巨大；裂隙岩体的有效渗透性在低密度区间内存在独有的尺度效应；将岩体的渗透性依据网络是否贯通分类后，两类渗透性平均值依据样本尺度呈现规律变化，该规律能够用帮组建立不同尺度下岩体有效渗透性之间的关系。
关键词：岩石力学裂隙岩体等效离散裂隙网络渗流低裂隙密度
Effective permeability of three-dimensional fractured rock with low fracture densities
HE Chen, YAO Chi
, SHAO Yulong, HUANG Fan, ZHOU ChuangbingSchool of Civil Engineering and Architecture, Nanchang University, Nanchang 330031, China
Abstract: The permeability of fractured rock is an important parameter in many underground projects, but the permeability characteristics are normally not well known and differ significantly from the effective permeability of rocks with high fracture densities that have been widely studied in previous research. This paper presents a numerical method to model the flows in three-dimensional fractured porous media based on the equivalent discrete fracture network model. This method uses two-dimensional triangular elements to simulate the flows in the fractures and the permeable matrix in the rock which predicts the fluid exchange between the fractures and the matrix. The results are used to derive an analytical expression for the conductivity coefficient of the equivalent element. Then, the Monte-Carlo method is used to analyze the permeability of fractured rock with various fracture densities and scales. The results show that the presence of a fracture network in the study area greatly influences the effective permeability and that the effective permeability of fractured rock has a unique scale effect for low fracture densities. After classifying the effective permeabilities with and without fracture networks in the domain, the average effective permeabilities for the two types show regular variations depending on the scale of the domain. This rule can be used to establish a relationship between the effective permeability of various size rock masses.
Key words: rock mechanicsfractured rockequivalent discrete fracture networkseepagelow fracture density
近年来，含裂隙地质构造中的渗流和溶质运移过程分析在许多环境和工业工程应用中变得越来越重要，如二氧化碳封存、原油及天然气开采、核废料处理等[1-3]。实际工程中，岩体的破碎程度因地区而异[4]，研究者通常会依据裂隙的密度和连通性采用合适的计算模型[5]。在裂隙密度较大的岩体中，贯通整个研究域的裂隙网络形成了主要的渗流通道，岩石基质中的孔隙对岩体宏观渗透性贡献非常微小。基于离散裂隙网络(discrete fracture network, DFN)描述的方法假设岩石基质不透水，渗流仅发生在裂隙中[6-7]。这种简化方式很好地诠释了复杂裂隙岩体渗流受裂隙主导的特征。但是，对于裂隙密度较低的岩体，裂隙之间的连通性欠佳，孔隙基质对岩体渗透性的贡献通常不能忽略，岩体的有效渗透性也不是基质和裂隙单独渗透性贡献的简单叠加[8]。
目前针对低裂隙密度岩体的渗透性特征的研究较少。低裂隙密度的岩体渗流问题通常比高裂隙密度更加复杂，这一方面是由于裂隙网络局部的渗透性特征变异性带来了大量的不确定性。另一方面，裂隙网络极小的几何差异可能会带来显著的影响，如网络贯通时与未贯通时岩体渗透性差异巨大[9-10]。因此，低裂隙密度渗透性需要结合非确定性方法和裂隙网络逾渗理论共同研究[11]。
数值方法在解决岩体渗透性问题时具有较高的灵活性，如有限元法[12-13]、有限体积法[8, 14]等。但由于裂隙厚度相较研究尺度过于微小，并且裂隙网络的几何特征非常复杂，依据岩体的实际三维几何离散求解域会带来极大的计算负担，该问题在需要同时描述裂隙和基质的模型中尤为明显。目前比较普遍的解决方法是引入“降维描述”的思想。Bogdanov等[8]将裂隙降维成没有厚度的有界平面，以基质单元之间界面的形式剖分，采用有限体积法求解了裂隙岩体渗流问题。Ren等[15-16]将裂隙和基质同时降维成一维管单元构成网络，使用等效后的管单元网络模拟裂隙多孔介质内的渗流行为，该方法无需定义裂隙和基质间的流量交互项。何忱等[17]提出了等效离散裂隙网络方法，该方法将裂隙和基质降维成没有厚度的三角形面单元构成的网络，在自然考虑裂隙-基质流量交互的前提下还保留了裂隙内面流的特征。但该模型物理参数的标定依赖于常密度网格，限制了模型的使用。
本文首先提出了一种裂隙多孔介质渗流数值模型。该模型建立在文[17]提出的等效离散裂隙网络思想之上，推导了等效单元物理参数的解析表达式。相较于文[17]中的方法，该方法消除了模型的网格依赖性，无需额外标定单元的物理参数，并且适应于变密度的网格。这些特性增强了该类方法的可靠性，也提高了程序的计算效率，更加适合解决大型工程问题。随后，使用该模型研究了低裂隙密度下岩体有效渗透性特征，结果展现了低裂隙密度下岩体渗透性特有的尺度效应以及裂隙网络贯通状态带来的重要影响。
1 数值模型1.1 等效离散裂隙网络(EDFN)模型等效离散裂隙网络(equivalent discrete fracture network，EDFN)模型是一种基于Galerkin有限元的裂隙多孔介质渗流模型。如图 1a所示，模型采用约束Delaunay四面体网格离散整个研究区域，图中绿色三角形单元用于描述裂隙，蓝色三角形单元用于等效基质。需要强调，模型中并没有三维的体单元，所有的计算单元都是Delaunay四面体两两之间的三角形交界面。这样做的好处是，模型对基质和裂隙的描述仅需通过一套控制方程即可完成，裂隙-基质间的流量交互在两类单元的公共边上自然考虑。模型以节点水头h为未知量，在每个三角形单元内部满足如下连续性方程：
| $\nabla \cdot \mathit{\boldsymbol{v}} = 0.$ | (1) |
 |
| 图 1 等效离散裂隙网络模型与基质单元导水系数推导过程 |
| 图选项 |
其中, v是单元内的局部流速矢量。忽略单元的储水作用，单元之间的公共边上，满足质量守恒方程:
| $\sum\limits_{i = 1}^m {{q_i}} = 0.$ | (2) |
1.2 对基质等效的解析方法如文[15, 17]中提到的，使用“降维描述”最大困难在于如何把降维过后的单元对原介质进行等效。对于EDFN模型，文[17]中采用了一种基于节点密度的方式赋予每个等效基质单元相同的导水系数，最终达到模型宏观渗透性与原介质相近的结果。但这是一种经验方法，显然此类模型有较强的网格依赖性，当网格失去尺寸均质的性质时(如裂隙网络引起的网格畸形)，EDFN模型会一定程度偏离正确的结果。
本文提出了一种解析方法，基于网格几何，赋予每个等效基质单元独特的导水系数T。为了便于读者理解，首先从二维单元开始。如图 1b所示，使用二维Delaunay网格离散一个矩形区域，此时模型单元降维成图中黑色线段。假设水力梯度方向从左往右，图中蓝色箭头是单元内的流速方向。图中红色直线是与Delaunay图对偶的Voronoi图。如图 1c所示，现任意取出一个计算单元ab，线段cd是与其对偶的Voronoi线段，将其作为ab单元对应的过流断面。假设基质单元每个边(面)控制区域内其流动有所在边(面)的压力控制，与对向点的压力无关。现需要满足：在任意已知的节点水头场中，使用基于三节点三角形有限元方法计算得到的流过cd断面的流量，等同于使用EDFN模型中流过ab单元的流量。即
| $\int_{{l_{cd}}} {{k_m}\left( {\frac{{\partial h}}{{\partial x}} + \frac{{\partial h}}{{\partial y}}} \right){\rm{d}}l = {T_{ab}}\left( {\frac{{{h_a} - {h_b}}}{{\overline {ab} }}} \right)} .$ | (3) |
接下来将模型扩展到三维空间。以图 1d中的三角形单元ijk为例，此时有3个垂直于ijk的Voronoi平面，他们共享了垂直于ijk的线段lm。假设相邻的DFEN计算单元平分了它们之间的过流断面(如图 1d中的opsl被ijk和ijg沿对角线平分)，最终得到了基本等效单元图 1e。仿照二维情况中的做法，先取单元ijk的子区域ijo讨论，需要满足：在任意已知的节点水头场中，使用四节点四面体有限元方法算得流经过流断面lpm的流量，需要等同于使用DFEN方法算得流经线段op的流量，即
| $\begin{array}{*{20}{l}}{{Q_{pml}} = {\rm{ }}{Q_{opl}} + {Q_{opm}} = {\rm{ }}{Q_{op}}.}\end{array}$ | (4) |
| $\begin{array}{l}\int_{{A_{opl}}} {{\mathit{\boldsymbol{v}}_{1, opl}} \cdot \mathit{\boldsymbol{n}}\;{\rm{d}}A} + \int_{{A_{opm}}} {{\mathit{\boldsymbol{v}}_{1, opm}} \cdot \mathit{\boldsymbol{n}}\;{\rm{d}}A} = \\ \;\;\;\;\; \int_{{l_{op}}} {{\mathit{\boldsymbol{v}}_2} \cdot b \cdot \mathit{\boldsymbol{n}}\;{\rm{d}}l.} \end{array}$ | (5) |
| $\begin{array}{l}{Q_{pml}} = \left( {\frac{{\overline {op} \cdot \overline {ol} }}{2} + \frac{{\overline {op} \cdot \overline {om} }}{2}} \right) \cdot {k_m} \cdot \frac{{{h_j} - {h_i}}}{{\overline {ij} }} = \\\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\overline {op} \cdot \overline {ml} \cdot {k_m}({h_j} - {h_i})}}{{2 \times \overline {ij} }}.\end{array}$ | (6) |
| ${Q_{op}} = \overline {op} \cdot {T_{ijo}} \cdot \frac{{{h_j} - {h_i}}}{{\overline {ij} }}.$ | (7) |
| ${T_{ijo}} = \frac{{\overline {ml} }}{2} \cdot {k_m}.$ | (8) |
| ${T_{ijk}} = {T_{ijo}} = \frac{{\overline {ml} }}{2} \cdot {k_m}.$ | (9) |
| $\rho ' = \frac{1}{2}\left\langle {AP} \right\rangle \rho .$ | (10) |
 |
| 图 2 裂隙网络示意图 |
| 图选项 |
 |
| 图 3 边界条件 |
| 图选项 |
其中：〈AP〉表示裂隙面积与周长乘积的算术平均值，对本例中的等径正六边形裂隙有〈AP〉=33L/8，ρ是裂隙的数量密度。裂隙-基质渗透性差异以及最终的裂隙岩体渗透性也同样采用无量纲形式表示，将R作为特征长度，有
| $\sigma ' = \frac{T}{{R{k_m}}}, K' = \frac{{{k_{{\rm{eff}}}}}}{{{k_m}}}.$ | (11) |
算例使用11组不同密度的样本，ρ′分别为0.487, 0.974, 1.46, 1.95, 2.44, 2.92, 3.41, 3.90, 4.99, 5.97和7.92。σ′变化从10增长到104。每组样本计算100次取渗透性平均值，无量纲化后记为〈K′〉, 结果如图 4所示，图中对比了3种模型的计算结果，数据分别来自Bogdanov等[8]，Vu等[9]和本文中的EDFN模型。结果表明，本文提出的EDFN模型与其他裂隙多孔介质模型计算结果吻合良好，该模型是有效的。
 |
| 图 4 3种模型的对比 |
| 图选项 |
1.4 模型的网格依赖性为考察该模型的网格依赖性，设计一组试验：固定裂隙几何不变，调整网格划分尺度及划分方式，计算裂隙岩体的有效渗透性。试验采用了3类代表性的裂隙岩体，分别为无裂隙、裂隙网络未贯通和裂隙网络贯通的岩体。在前处理中，设置网格节点在靠近裂隙边缘和相交处逐渐加密，因此无裂隙的样本为常密度网格，含裂隙的样本为变密度网格。计算结果如图 5所示，3类样本从节点数3 700左右的稀疏网格逐渐变化为节点数80 000左右的密集网格。可以看出，在3类裂隙网络中，本文的模型渗透性均不受网格密度的影响。
 |
| 图 5 不同网格密度下样本的渗透性 |
| 图选项 |
2 低裂隙密度岩体有效渗透性工程中对裂隙密度的评价指标并不唯一。对于裂隙网络岩体渗流问题，无量纲密度ρ′(式(9))综合考虑了裂隙的数量密度、几何形状与尺寸分布，大量文献[14, 18-20]将其作为研究岩体渗透性时的裂隙密度度量。依据裂隙网络连通性的差异，可以将裂隙密度划分为4个区间：极低密度区间(ρ′ < 1)，低密度区间(1≤ρ′ < 4)，较高密度区间(4≤ρ′ < 8)，高密度区间(ρ′≥8)。当ρ′ < 1时，裂隙几乎相互独立，部分****认为岩体的有效渗透性是ρ′的一次函数[21]；另一方面，当ρ′≥4时，裂隙网络的连通性达到了一个较高的水平，目前有许多解析模型可以较好地预测岩体的有效渗透性[14, 22-23]。
但是，对于低密度区间(1≤ρ′ < 4)的裂隙岩体，其渗透性有大量不确定性。一方面，这是因为裂隙网络可能贯通了研究区域，也可能没有贯通；另一方面，即使是同一个裂隙岩体，研究区域的尺度不同，渗透性也可能差异很大。为了减少不确定性带来的影响，本文采用Monte-Carlo模拟的方法，对不同裂隙密度和不同研究区域尺度的影响展开试验。
试验设计如下：采用半径为R的圆盘形裂隙；裂隙朝向在空间内随机均匀分布；研究区域是边长为L的立方体，裂隙形心遵循泊松过程布置在研究区域内。变化无量纲密度ρ′从0逐渐增大到6，以及研究区域尺度R/L=1/4、1/6、1/8、1/12和1/16，裂隙-基质间的渗透性差异取σ′=4×102(σ′见式(10))。对每组参数(ρ′, R/L)计算1 000次有效渗透性。由于在低裂隙密度下可能同时存在贯通和非贯通两种裂隙网络，分别对两类网络下的渗透性取统计平均值，记做〈K′〉，结果如图 6所示。发现以下现象：
 |
| 图 6 低裂隙密度岩体平均有效渗透性 |
| 图选项 |
1) 在低密度区间内(1≤ρ′ < 4)，将岩体有效渗透性依据裂隙网络是否贯通研究区域分开后，两者渗透性相差较大，且表现出明显的尺寸效应。
2) 随着研究域尺寸L的增大，贯通网络的〈K′〉逐渐降低，非贯通网络的〈K′〉逐渐升高。这意味着小尺寸样本的两类〈K′〉会“夹住”大尺寸样本的〈K′〉，这在推导无限大尺寸岩体的〈K′〉时能够作为一个有用的参考，如本例中无限大尺寸下的〈K′〉应在粉色虚线附近。
3) 上述尺度效应是在低裂隙密度下特有的。随着ρ′逐渐增大，裂隙网络连通性变强，尺度效应将逐渐消失。当ρ′>5时，所有尺度下的裂隙网络都贯通了研究域，〈K′〉都逐渐收敛在一条曲线上。
3 结论针对低裂隙密度下岩体渗流问题，本文首先提出了一种基于等效离散裂隙网络(EDFN)的三维裂隙岩体渗流模型，该模型能够自然考虑岩体渗流问题中裂隙、基质、裂隙-基质流量交换三者对有效渗透性的贡献。网格使用约束Delaunay四面体剖分方法得到，此类网格剖分方法具有快速、稳定和质量高等优势。针对EDFN模型中最困难的等效单元物理参数定义的问题，本文提出了一种解析模型，利用Delaunay-Voronoi图之间对偶的几何性质，消除了EDFN模型的网格依赖性。该模型能够处理任意复杂的裂隙网络，并且适应于变密度的网格。
随后，基于Monte-Carlo方法研究了裂隙密度区间(1≤ρ′ < 6)内裂隙岩体的有效渗透性。结果表明：在该密度区间内，岩石基质以及裂隙-基质间流量交换对渗透性的贡献都不能忽略；裂隙网络可能会贯通研究区域，也可能不会，2种情况下岩体渗透性表现出截然不同的特征，应当分开考虑。研究还发现了低裂隙密度区间内特有的尺度效应，其特征为：小尺度样本中将贯通网络和非贯通网络分开后得到的2种平均渗透性(文中通过〈K′〉表示)会“夹住”大尺度样本的这两种平均渗透性。这意味着，小尺度样本的有效渗透性统计结果构成了大尺度样本的上下界，结论可以为推测无限大尺度下低裂隙密度岩体有效渗透性提供依据。
参考文献
| [1] | XU T F, APPS J A, PRUESS K. Numerical simulation of CO2 disposal by mineral trapping in deep aquifers[J]. Applied Geochemistry, 2004, 19(6): 917-936. DOI:10.1016/j.apgeochem.2003.11.003 |
| [2] | TSANG C F, BERNIER F, DAVIES C. Geohydromechanical processes in the excavation damaged zone in crystalline rock, rock salt, and indurated and plastic clays[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(1): 109-125. DOI:10.1016/j.ijrmms.2004.08.003 |
| [3] | ZHANG H, FALCONE G, TEODORIU C. Modeling fully transient two-phase flow in the near-wellbore region during liquid loading in gas wells[J]. Journal of Natural Gas Science and Engineering, 2010, 2(2-3): 122-131. DOI:10.1016/j.jngse.2010.04.005 |
| [4] | 汪小刚, 贾志欣, 张发明, 等. 岩体结构面网络模拟原理及其工程应用[M]. 北京: 中国水利水电出版社, 2010. WANG X G, JIA Z X, ZHANG F M, et al. Rock mass structural plane network simulation principle and its engineering application[M]. Beijing: China Water Conservancy and Hydropower Press, 2010. (in Chinese) |
| [5] | LIU R C, LI B, JIANG Y J, et al. Review: Mathematical expressions for estimating equivalent permeability of rock fracture networks[J]. Hydrogeology Journal, 2016, 24(7): 1623-1649. DOI:10.1007/s10040-016-1441-8 |
| [6] | 叶祖洋, 姜清辉, 姚池, 等. 岩体裂隙网络非稳定渗流分析与数值模拟[J]. 岩土力学, 2013, 34(4): 1171-1178. YE Z Y, JIANG Q H, YAO C, et al. Formulation and simulation of non-steady seepage flow through fracture network in rock masses[J]. Rock and Soil Mechanics, 2013, 34(4): 1171-1178. (in Chinese) |
| [7] | 姚池, 姜清辉, 叶祖洋, 等. 裂隙网络无压渗流分析的初流量法[J]. 岩土力学, 2012, 33(6): 1896-1903. YAO C, JIANG Q H, YE Z Y, et al. Initial flow method for unconfined seepage problems of fracture networks[J]. Rock and Soil Mechanics, 2012, 33(6): 1896-1903. DOI:10.3969/j.issn.1000-7598.2012.06.045 (in Chinese) |
| [8] | BOGDANOV I I, MOURZENKO V V, THOVERT J F, et al. Effective permeability of fractured porous media in steady state flow[J]. Water Resources Research, 2003, 39(1): 257-263. |
| [9] | VU M N, POUYA A, SEYEDI D M. Effective permeability of three-dimensional porous media containing anisotropic distributions of oriented elliptical disc-shaped fractures with uniform aperture[J]. Advances in Water Resources, 2018, 118: 1-11. DOI:10.1016/j.advwatres.2018.05.014 |
| [10] | HYMAN J D, KARRA S, CAREY J W, et al. Discontinuities in effective permeability due to fracture percolation[J]. Mechanics of Materials, 2018, 119: 25-33. DOI:10.1016/j.mechmat.2018.01.005 |
| [11] | ADLER P M, THOVERT J F, MOURZENKO V V. Fractured porous media[M]. Oxford: Oxford University Press, 2013. |
| [12] | YAO C, HE C, YANG J H, et al. A novel numerical model for fluid flow in 3D fractured porous media based on an equivalent matrix-fracture network[J]. Geofluids, 2019, 2019: 9736729. |
| [13] | LANG P S, PALUSZNY A, ZIMMERMAN R W. Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map predictions[J]. Journal of Geophysical Research: Solid Earth, 2014, 119(8): 6288-6307. DOI:10.1002/2014JB011027 |
| [14] | MOURZENKO V V, THOVERT J F, ADLER P M. Permeability of isotropic and anisotropic fracture networks, from the percolation threshold to very large densities[J]. Physical Review E, 2011, 84(3): 036307. DOI:10.1103/PhysRevE.84.036307 |
| [15] | REN F, MA G W, WANG Y, et al. Unified pipe network method for simulation of water flow in fractured porous rock[J]. Journal of Hydrology, 2017, 547: 80-96. DOI:10.1016/j.jhydrol.2017.01.044 |
| [16] | REN F, MA G W, WANG Y, et al. Pipe network model for unconfined seepage analysis in fractured rock masses[J]. International Journal of Rock Mechanics and Mining Sciences, 2016, 88: 183-196. DOI:10.1016/j.ijrmms.2016.07.023 |
| [17] | 何忱, 姚池, 杨建华, 等. 基于等效离散裂隙网络的三维裂隙岩体渗流模型[J]. 岩石力学与工程学报, 2019, 38(S1): 2748-2759. HE C, YAO C, YANG J H, et al. A 3D model for flow in fractured rock mass based on the equivalent discrete fracture network[J]. Chinese Journal of Rock Mechanics and Engineering, 2019, 38(S1): 2748-2759. (in Chinese) |
| [18] | KOUDINA N, GARCIA R G, THOVERT J F, et al. Permeability of three-dimensional fracture networks[J]. Physical Review E, 1998, 57(4): 4466-4479. DOI:10.1103/PhysRevE.57.4466 |
| [19] | KHAMFOROUSH M, SHAMS K, THOVERT J F, et al. Permeability and percolation of anisotropic three-dimensional fracture networks[J]. Physical Review E, 2008, 77(5): 056307. DOI:10.1103/PhysRevE.77.056307 |
| [20] | BOGDANOV I I, MOURZENKO V V, THOVERT J F, et al. Effective permeability of fractured porous media with power-law distribution of fracture sizes[J]. Physical Review E, 2007, 76(3): 036309. DOI:10.1103/PhysRevE.76.036309 |
| [21] | SHAFIRO B, KACHANOV M. Anisotropic effective conductivity of materials with nonrandomly oriented inclusions of diverse ellipsoidal shapes[J]. Journal of Applied Physics, 2000, 87(12): 8561-8569. DOI:10.1063/1.373579 |
| [22] | 周创兵, 熊文林. 双场耦合条件下裂隙岩体的渗透张量[J]. 岩石力学与工程学报, 1996, 15(4): 338-344. ZHOU C B, XIONG W L. Permeability tensor for jointed rock masses in coupled seepage and stress field[J]. Chinese Journal of Rock Mechanics and Engineering, 1996, 15(4): 338-344. DOI:10.3321/j.issn:1000-6915.1996.04.006 (in Chinese) |
| [23] | 陈益峰, 周创兵, 胡冉, 等. 大型水电工程渗流分析的若干关键问题研究[J]. 岩土工程学报, 2010, 32(9): 1448-1454. CHEN Y F, ZHOU C B, HU R, et al. Key issues on seepage flow analysis in large scale hydropower projects[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(9): 1448-1454. (in Chinese) |
