近年来，含裂隙地质构造中的渗流和溶质运移过程分析在许多环境和工业工程应用中变得越来越重要，如二氧化碳封存、原油及天然气开采、核废料处理等^[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)

其中, v是单元内的局部流速矢量。忽略单元的储水作用，单元之间的公共边上，满足质量守恒方程:

$\sum\limits_{i = 1}^m {{q_i}} = 0. $

(2)

其中，m为与公共边相接的单元数量，q_i表示从公共边流向单元i的流量。

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)

其中: k_m是基质的渗透系数，ab为线段ab的长度(后文中所有上划线均表示线段长度)。注意到cd是Voronoi图的一条边，可以找到任意相互连接的过流断面(如图 1a中的加粗红色多段线)满足上述方程，使得EDFN模型表现出的宏观渗透性在网格最小尺度之上都是正确的。

接下来将模型扩展到三维空间。以图 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)

其中: v₁是三维有限元方法基质体单元内的流速矢量，v₂是DFEN方法中ijk面单元内的流速矢量。对于式(5)的左边，采用四节点四面体有限元格式，得到流经过流断面pml的流量，

$\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)

其中, h_i和h_j分别为i、j两节点的水头值。对于式(5)的右边，采用三节点三角形有限元格式，假设子单元ijo的导水系数为T_ijo，有

${Q_{op}} = \overline {op} \cdot {T_{ijo}} \cdot \frac{{{h_j} - {h_i}}}{{\overline {ij} }}. $

(7)

将式(6)、式(7)代入式(4)中，化简得

${T_{ijo}} = \frac{{\overline {ml} }}{2} \cdot {k_m}. $

(8)

注意到，式(8)中T_ijo只是线段ml的长度与基质渗透系数k_m的函数，而图 1e中另外2个过流断面rml和qml也共享了ml线段，因此整个三角形单元ijk上的导水系数是一致的，即

${T_{ijk}} = {T_{ijo}} = \frac{{\overline {ml} }}{2} \cdot {k_m}. $

(9)

1.3 模型的验证

对于复杂裂隙网络条件下的裂隙多孔介质渗流，Bogdanov等^[8]提出过一个经典算例。如图 2所示，在边长为L的立方体区域中布置正六边形裂隙，裂隙外接圆半径R=1/4 L，裂隙朝向和形心坐标随机均匀分布。边界条件如图 3所示，一对边界作为上下游，施加水头边界条件，其他边界为不透水边界。算例考察了模型在不同裂隙密度和不同裂隙-基质渗透性差异下的表现。对裂隙密度的描述采用无量纲密度ρ′ ^[11]，

$\rho ' = \frac{1}{2}\left\langle {AP} \right\rangle \rho . $

(10)

其中：〈AP〉表示裂隙面积与周长乘积的算术平均值，对本例中的等径正六边形裂隙有〈AP〉=33L/8，ρ是裂隙的数量密度。裂隙-基质渗透性差异以及最终的裂隙岩体渗透性也同样采用无量纲形式表示，将R作为特征长度，有

$\sigma ' = \frac{T}{{R{k_m}}}, K' = \frac{{{k_{{\rm{eff}}}}}}{{{k_m}}}. $

(11)

其中，k_eff为裂隙岩体的渗透系数，T为裂隙的导水系数。

算例使用11组不同密度的样本，ρ′分别为0.487, 0.974, 1.46, 1.95, 2.44, 2.92, 3.41, 3.90, 4.99, 5.97和7.92。σ′变化从10增长到10⁴。每组样本计算100次取渗透性平均值，无量纲化后记为〈K′〉, 结果如图 4所示，图中对比了3种模型的计算结果，数据分别来自Bogdanov等^[8]，Vu等^[9]和本文中的EDFN模型。结果表明，本文提出的EDFN模型与其他裂隙多孔介质模型计算结果吻合良好，该模型是有效的。

1.4 模型的网格依赖性

为考察该模型的网格依赖性，设计一组试验：固定裂隙几何不变，调整网格划分尺度及划分方式，计算裂隙岩体的有效渗透性。试验采用了3类代表性的裂隙岩体，分别为无裂隙、裂隙网络未贯通和裂隙网络贯通的岩体。在前处理中，设置网格节点在靠近裂隙边缘和相交处逐渐加密，因此无裂隙的样本为常密度网格，含裂隙的样本为变密度网格。计算结果如图 5所示，3类样本从节点数3 700左右的稀疏网格逐渐变化为节点数80 000左右的密集网格。可以看出，在3类裂隙网络中，本文的模型渗透性均不受网格密度的影响。

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×10²(σ′见式(10))。对每组参数(ρ′, R/L)计算1 000次有效渗透性。由于在低裂隙密度下可能同时存在贯通和非贯通两种裂隙网络，分别对两类网络下的渗透性取统计平均值，记做〈K′〉，结果如图 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′〉表示)会“夹住”大尺度样本的这两种平均渗透性。这意味着，小尺度样本的有效渗透性统计结果构成了大尺度样本的上下界，结论可以为推测无限大尺度下低裂隙密度岩体有效渗透性提供依据。