关于孔隙-裂缝耦合模型中水力特性的研究，用于预测矿井底板水的涌入

《Frontiers in Earth Science》：Study on hydraulic properties in a pore-fracture coupled model for predicting water inrush from mine floors

　　摘要
采矿引发的破坏带是主要的地下水流动通道之一。由于裂缝分布的复杂性和裂缝中地下水流动的特性，目前计算采矿破坏带内裂缝的导水性是一个热门且困难的研究课题。为了同时模拟底板整个岩体的渗透性和裂缝的水力传导性，建立了一个孔隙-裂缝多孔介质的数值模拟模型，以研究不同裂缝结构参数对底板水力传导性的影响。在各种裂缝结构参数和渗透性条件下，研究了水压和水流速度的空间分布。模拟结果表明，底板压力分布与实际情况一致。裂缝位置的水压呈现下降趋势，裂缝内的压力梯度小于周围多孔岩体中的压力梯度。该模型不仅能够模拟和计算复杂孔隙-裂缝介质对流体的阻力，还能简化计算过程，为底板裂缝带的水力传导性模拟和计算提供了一种新方法。

1 引言
矿井水灾是中国煤矿安全生产面临的主要灾害之一（Gu等人，2026；Cao等人，2022；Sun等人，2022）。随着浅层煤炭资源的逐渐枯竭，近年来中国的煤炭开采深度不断加大。深埋的 lower 岩层位于地质环境复杂的区域，具有较高的原位应力和水压（Cao等人，2022；Wang等人，2022）。关于矿井底板突水机制的基础理论研究是一个跨学科的技术问题，涉及流体力学、岩石力学、水文地质学等多个领域（Liang等人，2020；Zhang等人，2017）。中国专家学者进行了大量深入研究，并提出了一系列矿井底板水灾预防和控制的理论体系，包括“下三带”理论（Yin等人，2019）、突水系数法（Yan等人，2021；Li等人，2017）以及关键层理论（Liang等人，2020；Li等人，2022）。这些理论成果为矿井底板突水的预测和预防奠定了坚实的基础，对含水层下煤炭资源的安全开采理论和技术发展做出了重要贡献。

矿井底板突水发生的两个必要条件是存在可靠的水源和有效的水力传导路径（Hao等人，2020）。承压含水层中的地下水是底板突水的主要水源（Shi等人，2019），而含水层水压决定了地下水是否能通过水力传导路径迁移到采场（Li等人，2020；Sun等人，2023）。地下岩体中的破坏带，包括断层和坍塌柱等构造裂缝以及采矿引发的破坏带，是主要的水力传导路径（Mu等人，2020；Liu等人，2023），其导水性与裂缝开口度、裂缝密度等裂缝结构参数密切相关。通过单轴压缩对含裂隙的双层复合岩样进行力学性质和微裂纹演化的数值研究（Ma等人，2023）。利用离散裂缝网络模型对节理岩体中的隧道突水灾害进行灌浆优化（Huang等人，2025）。裂缝成为主要的水力传导通道的根本原因在于，裂缝中的地下水流动阻力低于低渗透性隔水层多孔岩体中的阻力（Ogata等人，2018），这种阻力直接决定了裂缝的导水性。地下水在从含水层通过裂缝路径迁移到采场过程中需要克服自身的重力以及路径中的流动阻力，其迁移驱动力来自含水层水压（Zhang，2021）。因此，研究水压与裂缝结构参数之间的定量关系对于准确评估底板突水风险至关重要。

由于裂缝嵌入多孔岩基质中，必须同时计算裂缝及其周围多孔介质中的地下水流动参数（Xue等人，2019）。多孔介质中流体的复杂流动行为通常通过体积平均法进行分析，该方法将流体物理量转换为平均参数，而裂缝中的流体参数可通过求解动量守恒方程获得。由于两种计算方法存在差异，同时计算裂缝和多孔岩体中的流动阻力较为困难，导致关于煤层底板裂缝及周围岩体中水流速度和水压的空间分布特征的研究相对有限（Xue等人，2019；Sheng，2006）。

为了解决上述研究空白，本研究使用COMSOL数值模拟软件，将基于体积平均法的达西定律与流体动量守恒方程相结合，建立了一个孔隙-裂缝多孔介质数值模型，能够同时计算裂缝及其周围多孔介质中的地下水流动参数。通过与孔隙型多孔介质数值模型进行比较，系统研究了不同裂缝结构参数（裂缝数量、开口度、倾角）对煤层底板裂缝带导水能力的影响机制，并从流体流动阻力的角度揭示了底板突水的成灾机制，为准确预测和预防矿井底板突水提供了理论基础和数值方法。

2 方法
2.1 控制方程
地下煤岩体被归类为多孔介质，其中流体的流动非常复杂。实际上，孔隙尺度上的实际流动参数（如物理速度）并不直接分析；而是采用体积平均法将物理速度转换为渗流速度。从达西实验得出的渗流速度是一种典型的平均速度，表达式为（Sheng，2006）：
$$v = -\κ\mu(\nabla p - \rho g)$$
式中：$v$ 是流体的达西速度，$\κ$ 是多孔介质的渗透性，$\mu$ 是流体的动力粘度，$p$ 是孔隙水压，$\rho$ 是流体密度，$g$ 是重力加速度向量。

多孔介质中流体流动的连续性方程为（Kong等人，2005）：
$$\frac{\partial}{\partial t}(\phi\rho f) + \nabla \cdot (\rho f v) = 0$$
式中：$\phi$ 是多孔介质的孔隙度。

结合方程1和2，得出多孔介质渗流的控制方程：

地下煤岩体裂缝中的流体流动通常通过两种方法进行研究：
- 将裂缝视为孔隙并应用达西定律，与多孔介质流动分析一致；
- 直接描述裂缝中的实际流动状态，由流体动量守恒方程控制（见方程3，Xue等人，2019）：
$$\rho \frac{\partial v}{\partial t} + \rho (v \cdot \nabla p) = -\nabla p + \nabla \cdot \left(\mu [\nabla v + (\nabla v)^T \right] + \rho g$$
式中：$\mu$ 是裂缝内流体的物理速度。

对于多孔岩体，动量守恒方程简化为布林克曼方程（见方程4）：
$$\rho f \phi \left[ \frac{\partial v}{\partial t} + (v \cdot \nabla p) v \phi \right] = -\nabla p + \nabla \cdot \left(1 - \phi \mu (\nabla v + (\nabla v)^T \right) - (\mu k + \rho \beta |v|) v + \rho f g$$
在COMSOL中使用布林克曼物理场接口时，将裂缝区域定义为裂缝流，其他区域定义为多孔介质。为了计算煤层顶板和底板中的采矿诱发性破坏带，通常假设多孔介质具有弹性。基于广义胡克定律并考虑外部应力（如孔隙水压、热应力）， constitutive 方程可表示为（见方程5）：
$$\nabla \cdot (\sigma + \sigma_{\text{ext}}) + F = 0$$
式中：$F$ 是体力向量。

对于裂隙岩体，采用岩石损伤变量来量化岩石破坏的程度及其演化过程。岩石应力-应变曲线的不线性通常归因于加载下的渐进性损伤、微裂纹的形成和扩展——尤其是在拉应力作用下，脆性行为更为明显。因此，弹性损伤力学的 constitutive 关系适用于描述岩石细观单元的力学性质（Jin等人，2021）。

自20世纪80年代以来，全球学者提出了各种损伤模型（弹性损伤、弹塑性损伤等）。根据应变等效原理——即应力 $\sigma$ 在损伤材料上引起的应变等同于有效应力在完整材料上引起的应变——定义了以下损伤变量（见方程6–8，Zhu等人，2003）：
$$D_d = \begin{cases}
0 & F_1 < 0, F_2 < 0, 1 - \left|\varepsilon_{t0}\varepsilon_3\right| < n \\
F_1 = 0 & 1 - \left|\varepsilon_{c0}\varepsilon_1\right| < n \\
F_2 = 0 & F_1 < 0, 1 - \varepsilon_{c0}\varepsilon_1 < n
\end{cases}$$
其中：$D_d$ 是损伤变量，$F_1$ 和 $F_2$ 是基于最大拉应力准则和莫尔-库仑准则的损伤判断应力状态函数。当 $F_1 = 0$ 时发生拉应力损伤，当 $F_2 = 0$ 时发生剪切损伤。

将上述损伤变量纳入考虑后，受损完整岩层的渗透性可表示为（见方程9，Jin等人，2021）：
$$k_D = k \exp(\alpha)$$
式中：$\alpha$ 是由于损伤引起的渗透性影响系数（本研究中取 $\alpha = 5$）。

2.2 孔隙-裂缝型数值模拟模型
孔隙型多孔介质数值模型通过赋予不同的孔隙度和渗透性值来模拟煤层底板的孔隙特性和渗透性。当底板存在裂缝时，需要一个孔隙-裂缝模型来同时表示多孔介质的渗透性和裂缝的导水性。为了明确裂缝结构参数对底板水传导性的影响，并突出其与孔隙型模型的差异，同时建立了孔隙型和孔隙-裂缝多孔介质数值模型。

2.2.1 孔隙型多孔介质模型
2.2.1.1 几何模型
孔隙型模型通过孔隙度来表征孔隙结构，使用体积平均法将物理流体速度转换为渗流速度的乘积；因此无需显式的孔隙几何形状。如补充图S1所示，该模型由3米厚的煤层、20米厚的底板隔水层和5米厚的基岩含水层组成，模型总宽度为30米。

2.2.1.2 边界条件
为了量化不同渗流速度下隔水层中的流动阻力（所需水压），将含水层下边界设为速度入口（每次模拟方案中的速度可变），采场上边界设为压力出口（0帕），所有其他边界设为零速度边界。

2.2.1.3 材料参数
关键参数包括地下水密度（$1 \times 10^3 \, \text{kg/m}^3$）、动力粘度（$1 \times 10^{-3} \, \text{Pa}\cdot \text{s}$）、隔水层渗透性（$1 \times 10^{-15} \, \text{m}^2$）和含水层渗透性（$1 \times 10^{-12} \, \text{m}^2$）。

2.2.2 孔隙-裂缝型多孔介质模型
孔隙-裂缝模型包含多孔岩体和离散裂缝。多孔介质由孔隙度定义（无显式孔隙几何形状），裂缝穿透隔水层连接含水层和采场。如补充图S2所示，模型尺寸为30米（宽）× 28米（高），包含3米厚的煤层、20米厚的隔水层和5米厚的含水层。裂缝结构参数根据模拟方案进行调整。

2.2.3 孔隙-裂缝型多孔介质模型
图1a和b展示了完整（无裂缝）和裂隙底板的空间压力分布（裂缝开口度 = 20毫米，入口速度 = $3 \times 10^{-5} \, \text{m/s}$）。裂隙案例模拟了单个裂缝（开口度20毫米，倾角18°）嵌入渗透性与孔隙型模型相匹配的多孔介质中。
**图1**：压力和速度。(a) 无裂缝条件下的压力分布；(b) 有裂缝条件下的压力分布；(c) 无裂缝条件下的速度分布；(d) 有裂缝条件下的速度分布。
- 水平方向：完整底板上的压力均匀分布，但裂隙底板上的压力呈凹陷向下分布（裂缝处压力低，周围岩石处压力高）。这是由于在均匀隔水层中流动阻力均匀，而在裂缝中的阻力较低。
- 垂直方向：完整底板上的压力梯度均匀，但在裂缝中的压力梯度较小。均匀的渗透性使得完整底板上的垂直流动阻力一致，而裂缝中的阻力降低导致局部压力梯度减小。

图1c和d展示了完整和裂隙底板的速度分布（单个裂缝，开口度200毫米，倾角18°，入口速度 = $3 \times 10^{-5} \, \text{m/s}$）。在完整底板中，由于渗透性均匀、压力梯度一致和质量守恒，渗流速度在空间上均匀分布。在断裂的岩石中，流体在裂缝内的流速远高于周围岩石中的流速，因为地下水更倾向于通过阻力较小的裂缝路径流动。在x=2米的测量线（图2a）上，含水层中的渗流速度大于裂缝附近多孔介质中的渗流速度。根据质量守恒原理，来自含水层的总流入量等于从隔水层流出的总量；大部分水流集中在裂缝中，导致含水层中的流速高于周围多孔岩石中的流速。

图2 渗流速度。（a）沿测量线的渗流速度；（b）不同裂缝开口下的压力-流速关系。为了单独研究裂缝的阻力，在分析中排除了静水压力。如图2b所示（入口速度从3×10^-8 m/s到4×10^-5 m/s，裂缝开口从5 mm到200 mm），裂缝入口处的平均水压在5 mm和200 mm开口之间存在大约六个数量级的差异。较大的裂缝开口可以降低沿路径的流动阻力，从而降低所需的含水层压力。入口速度与裂缝入口压力呈正相关。裂缝区域的水传导性与裂缝密度密切相关。3条裂缝（开口40 mm，入口速度1×10^-5 m/s）的压力分布分别绘制在补充图S3和S4中。与单裂缝行为类似，裂缝处的压力呈现出向下凹陷的分布，中间区域由于相邻裂缝之间的水力相互作用也显示出压力降低。随着裂缝数量从1条增加到5条（图3），含水层入口处的平均压力呈单调下降趋势。裂缝密度的增加扩大了总的横截面积，降低了裂缝区域的整体流动阻力。这证实了更多的裂缝发育能够提高水传导性，减少给定流量所需的压力，并增加涌水风险。

图3 含水层入口边界处的平均压力与裂缝数量之间的关系。

2.3 孔隙-随机裂缝水力-力学-损伤（HMD）耦合模型
采矿引起的矿顶和矿底岩层的变形和破坏形成了导水路径（例如，矿顶的导水裂缝区和矿底由采矿引起的破坏区）。地下水在裂缝中的流动类似于管道流动，而在周围岩石中的流动则遵循多孔介质的渗流规律。为了准确模拟这两种流动机制，必须同时求解动量守恒方程和多孔介质渗流方程。上述建立的孔隙-裂缝模型能够同时捕捉到裂缝管道流动和完整岩石的渗流，真实地再现了采矿受损区域中的地下水流动情况。为了进一步模拟采矿引起的损伤，将应力方程和损伤方程纳入模型中。引入随机裂缝网络以更好地表示煤层顶板和底板中的天然裂缝分布。因此，开发了一种孔隙-随机裂缝水力-力学-损伤（HMD）耦合模型，并以煤炭矿山工作面作为工程案例。

2.3.1 工作面概述
工作面的平均垂直深度为617米。第2号煤层的平均厚度为9.3米，倾斜角度为4°–6°。

2.3.2 建立孔隙-随机裂缝HMD耦合模型
如补充图S5所示，根据工程地质条件构建了一个几何模型：长度745米，高度662米，煤层厚度9.3米，工作面倾斜长度241米，并为矿顶和底板分配了随机裂缝网络。矿顶和矿底中的采矿引起损伤区在补充图S6中展示。在水力-力学-损伤耦合作用下，损伤沿着裂缝传播并连接离散的裂缝，形成连续的导水路径。在含水层压力分别为0、1、2和3 MPa（出口压力=0 MPa）的情况下，渗流速度分布如图4所示。当压力为0 MPa时，渗流速度为零。随着含水层压力的增加，地下水更倾向于通过相互连接的受损裂缝流动，最大渗流速度从5.51×10^-4 m/s上升到1.67×10^-3 m/s。模型两侧采用滚轴支护，在下边界施加固定约束，在上边界施加均匀分布的荷载。该模型采用非结构化网格。

图4 孔隙-裂缝模型的渗流速度分布。
对于孔隙控制模型（入口压力3 MPa，图4），最大渗流速度为3.02×10^-4 m/s，显著低于孔隙-裂缝模型的结果。这是因为孔隙-裂缝模型将裂缝视为离散的管道流动路径（低阻力），而不是均匀的多孔介质（高阻力），从而在受损区域获得了更高的整体渗透性和流速。因此，孔隙-裂缝框架更准确地反映了现场地质和水力条件。

3 结论与讨论
基于多孔介质渗流和流体力学，开发了一个孔隙-裂缝数值模型，用于同时模拟裂缝管道流动和周围多孔岩石的渗流。该模型量化了在不同裂缝结构参数下的地下水压力和速度的空间分布。模拟结果与现场观察一致：周围岩石中的流动阻力大于裂缝中的流动阻力，形成了向下凹陷的水平压力剖面以及裂缝内部较小的垂直梯度。裂缝中的流速远高于完整岩石中的流速。水传导性与裂缝开口大小和数量相关，随着裂缝发育程度的增加，给定流量所需的压力降低，涌水风险增加。所提出的模型将达西定律与流体动量守恒方程结合起来，提高了对完整岩石和裂缝区域的物理真实性。它准确地再现了复杂孔隙-裂缝介质中的流动阻力，同时简化了双流动机制的计算，为评估煤层底板裂缝区的导水能力提供了一种新的数值方法。该模型假设煤层顶板和底板的破坏区域等同于随机裂缝，未来的研究应根据实际情况进行建模。