关于粗糙单一裂隙中非达西效应与水力开度关系研究的学术报告

本文旨在向国内学术界介绍一篇于2021年发表于*Water Resources Research*期刊上的重要研究论文，题为“Characterizing the relationship between non-Darcy effect and hydraulic aperture in rough single fractures”。该研究由合肥工业大学资源与环境工程学院的邢坤、钱家忠（通讯作者）、马海春、罗乾坤和马雷（通讯作者）共同完成。

一、 研究的学术背景

本研究隶属于地质工程与水文地质学交叉领域，核心关注点是岩体裂隙中的流体流动问题。裂隙作为地壳中广泛存在的地质介质，是地下水、污染物、气体及二氧化碳等物质在低渗透性岩体中运移的主要通道。准确描述裂隙中的流体流动规律，对于二氧化碳封存、核废料处置、地热开采及地下水污染防治等工程至关重要。

传统的裂隙流理论基于立方定律（Cubic Law），该定律将裂隙简化为理想光滑的平行板模型。然而，自然界中的真实裂隙具有强烈的非均质性和复杂的几何结构，导致在高流速、高雷诺数（Reynolds number）或复杂几何条件下，达西定律（Darcy’s law）和立方定律往往失效，出现流速（或流量）与水动力梯度（或压力梯度）之间的非线性关系，这种现象被称为非达西效应（non-Darcy effect）。长期以来，研究者们通常采用临界雷诺数（critical Reynolds number, Rec）或非达西效应因子（non-Darcy effect factor, e）来量化流动从达西流向非达西流的转变起点。然而，这些参数的确定存在争议，例如特征长度的选择不统一，且传统的判断标准（如e=0.1）受到几何形状的影响而波动，其背后的物理机制也未完全统一。

水力开度（hydraulic aperture, eh）是表征裂隙几何和水力特性的关键参数，通常由立方定律反算得出。以往研究多默认水力开度仅取决于裂隙的几何属性（如机械开度、粗糙度），而与流体的惯性效应无关。但近期一些研究暗示水力开度可能随流动状态（惯性效应）变化。基于此，本研究提出了一个创新性设想：水力开度作为一个同时耦合了介质几何特性和流体惯性效应的参数，可能为表征非达西效应提供一种更简便、更直接的新方法。因此，本研究的主要目标是系统地探索利用水力开度来表征粗糙单一裂隙中非达西效应的可行性。

二、 详细研究流程

本研究主要采用基于有限元方法的计算流体动力学（CFD）软件COMSOL Multiphysics®进行二维数值模拟。整个工作流程严谨而系统，包含以下几个核心步骤：

1. 物理模型构建：研究选择构建二维规则粗糙单一裂隙模型而非三维模型，主要基于计算效率和对特定物理过程进行理想化研究的考虑。研究指出，大量前期工作已证明二维Navier-Stokes方程用于研究裂隙流是可靠的，并且可以作为理解更复杂三维模型的理论基础。为了探究粗糙度的影响，研究构造了三种具有规则粗糙元（凸起）的裂隙模型：矩形、梯形和三角形。通过设定不同的凸起高度（δ，0.4-1.0 mm）和最大开度（b_max，1.2-2.0 mm），生成了具有不同相对粗糙度（r/rsd）的一系列裂隙样本。所有粗糙元通过周期函数生成，确保了粗糙元密度的一致性。平均机械开度（em）和相对粗糙度通过附录中描述的统计方法计算。

2. 控制方程与边界条件设定：流体流动的控制方程为不可压缩牛顿流体的Navier-Stokes方程和连续性方程。模型的左端设为压力入口，右端设为压力出口（压力为0 Pa）。为获得精确的数值解，研究进行了细致的网格独立性分析。通过对比常规、细化、更细化、极度细化和超细化五种网格划分方案下的流速-水力梯度关系曲线，确定了当网格达到“极度细化”级别时，数值解趋于稳定。因此，最终采用“极度细化”网格进行所有模拟，在保证结果准确性的同时兼顾了计算效率。具体的边界条件、流体物性参数（水）和求解器设置（稳态求解）均在文中表格详细列出。

3. 模拟方案与数据分析流程：通过逐渐增大入口压力（即增大水力梯度J），模拟从低流速（达西流区）到高流速（非达西流区）的完整流动过程。对于每个流动状态，软件计算出流场分布（速度、压力）和总流量（q）。随后，基于立方定律（公式1）反算出该状态下的水力开度（eh）。同时，利用Forchheimer方程（公式4）对模拟得到的压力梯度（∇p）与流量（q）数据进行拟合，获得粘性项系数（a1）和惯性项系数（b1），进而计算每个状态对应的非达西效应因子（e，公式5）。e的物理意义是惯性项造成的压力损失占总压力损失的比例。此外，研究还利用流场数据，通过内置探针和自行编写的MATLAB代码，量化了低流速区（速度低于平均速度一个数量级）的体积以及回流区（recirculation zones, RZs）的体积，以探究非达西效应的微观机制。用于识别回流区的MATLAB代码基于速度分量方向与主流方向相反的特性开发，并在附录E中给出。

三、 主要研究结果

研究通过系统的数值模拟，获得了一系列揭示水力开度与非达西效应之间内在联系的深刻结果。

1. 水力开度表征非达西效应起始点：结果清晰地表明，在达西流动阶段（低e值），水力开度保持恒定，该值被定义为初始水力开度（eh0）。eh0仅与裂隙的几何属性（平均机械开度em、相对粗糙度r）相关，与惯性效应无关。这与传统观点部分吻合。然而，当惯性效应增强到某一临界点时，水力开度开始从eh0处下降。这个转折点被定义为非达西效应的起始点，对应的临界非达西效应因子记为e_d。研究发现，利用水力开度的下降来判定非达西效应起始，比传统标准e=0.1更为敏感和直接，因为e_d通常远小于0.1（例如在示例中e_d=0.0317）。这意味着非达西效应在惯性压力损失占比还很小时就已经开始。

2. 非达西效应机制与临界粗糙度的发现：通过分析流场，研究发现非达西效应的产生与粗糙元导致的低流速区密切相关。随着惯性效应增强，低流速区的体积逐渐增大。当低流速区发展到一定规模时，流场中高速主流区与低流速区之间的水力特性差异加剧，导致速度场非均质性增强，流动阻力增大，从而引发非达西效应（表现为eh下降）。研究还首次发现了一个有趣的“临界相对粗糙度（rc）”现象。当相对粗糙度r小于rc时，裂隙中始终存在一条最窄的优势流道（narrowest preferential flow path, NPFPS）。随着r增大，NPFPS变窄，低流速区与优势流道间的差异减小，反而抑制了非达西效应的早期发生，因此e_d随r增大而增大。然而，当r超过rc后，优势流道完全被低流速区阻塞，粗糙度的进一步增加直接导致低流速区体积显著扩大，使得非达西效应更容易发生，因此e_d随r增大而减小。这一发现修正了以往“粗糙度越大，越易发生非达西流”的简单认知，揭示了粗糙度影响的双阶段性。

3. 高惯性效应下水力开度的演变阶段：在非达西效应开始后，随着惯性效应（e值）持续增强，水力开度总体呈下降趋势，但下降过程并非单调平滑，而是呈现出多阶段性，可分为：稳定下降区、波动区、再稳定区和再波动区。通过分析回流区的形成与发展，研究将这种波动归因于回流区的动态演化。在稳定下降区，回流区尚未形成。当e达到一定值，回流区开始产生并进入活跃发展阶段（对应波动区），其体积增长率不断变化，扰动了水力开度的下降轨迹。当回流区发展充分、体积增长率达到峰值后开始下降时，水力开度进入再稳定下降区。后续可能因新的回流区在其他位置形成，再次引发波动（再波动区）。这一结果表明，水力开度的波动本身反映了流场结构的复杂变化，但其总体下降趋势仍稳健地表征了非达西效应的增强。

4. 基于水力开度的非达西效应量化模型：为了直接利用水力开度量化非达西效应的强度，研究对水力开度的衰减进行了归一化处理，定义了归一化水力开度衰减量（(eh0 - eh)/eh0）。研究发现，该衰减量与e之间存在良好的线性关系。通过对模拟数据进行拟合，得到了经验公式：(eh0 - eh)/eh0 = 0.684e （R² = 0.984）。为了验证该关系的普适性，研究引用了其他学者（Rong et al., 2020; Zhang et al., 2019）的室内实验数据进行检验，拟合结果（R² = 0.812）支持了这一线性关系。此外，研究还基于自身数据，建立了初始水力开度eh0与几何参数（em, r）的经验关系：eh0 = 1.756em / (1 + 2.34r)。将两者结合，可以得到一个直接基于几何参数和eh测量值来预测e的模型。

四、 研究结论与价值

本研究得出核心结论：水力开度（eh）是一个能够有效表征粗糙单一裂隙中非达西效应的关键参数。它耦合了裂隙的几何特性和流体的惯性效应。利用eh从恒定值（eh0）开始下降的点，可以更灵敏地判定非达西效应的起始。非达西效应的微观机制与粗糙裂隙表面诱导的低流速区（及其特殊形式——回流区）的发育密切相关。研究首次发现了临界相对粗糙度（rc）的存在，它通过影响优势流道和低流速区的发育格局，以双阶段模式影响非达西效应的起始。此外，研究还建立了基于水力开度衰减来量化非达西效应强度（e）的线性预测模型，并经实验数据初步验证。

本研究的科学价值在于：1）为理解非达西效应的物理机制提供了新的视角（低流速区发育），并揭示了粗糙度影响的复杂非线性特征（临界粗糙度）；2）提出并验证了一种基于水力开度的、更简便直观的非达西效应表征与量化新方法，避免了传统方法中特征长度选择争议和Forchheimer系数拟合困难等问题；3）支持了“水力开度具有惯性效应依赖性”这一前沿观点，深化了对这一经典参数的认识。

其应用价值在于：所提出的方法和模型有望用于现场裂隙介质渗透性和导水性的更准确评估，避免因忽略非达西效应而高估这些参数，从而为地热开发、二氧化碳地质封存、核废料处置库安全评估及地下水污染修复等工程中的水文地质计算提供更可靠的理论依据。

五、 研究亮点

1. 观点新颖：首次系统提出并证实了“水力开度可作为表征非达西效应的核心参数”这一创新性设想。
2. 发现重要：首次在数值模拟中观察到并阐述了“临界相对粗糙度”现象及其对非达西效应起始的双阶段影响机制，是对传统认识的深化和修正。
3. 方法整合：巧妙地将计算流体力学模拟、流场精细后处理（低流速区、回流区量化）和参数化建模相结合，研究流程完整、逻辑严谨。
4. 模型实用：建立的基于水力开度的非达西效应量化模型形式简洁，并经过了初步的实验数据交叉验证，显示了良好的应用潜力。

六、 其他有价值内容

研究在附录中详细给出了平均机械开度、相对粗糙度的计算方法，三种规则粗糙裂隙的周期函数表达式，最窄优势流道的识别方法，以及用于计算回流区体积的MATLAB核心代码。这些内容增强了研究的可重复性和透明度，为其他研究者提供了宝贵的工具和参考。此外，研究也坦承了其基于二维简化模型的局限性，并指出未来需要在更复杂的三维模型中进行验证，体现了科学的严谨性。