利用改进的DEM-PNM耦合方法对渗流过程中堵塞机制进行多尺度研究

《Powder Technology》：Multi-scale investigation of clogging mechanisms in suffusion using an improved DEM-PNM coupling method

孔隙堵塞与流体迁移多尺度耦合机制研究。本文提出了一种改进的DEM-PNM耦合方法，通过区分粗细颗粒与流体的相互作用，减少网格更新频率，将计算效率提升6.1倍。研究揭示了孔隙堵塞的三态演化规律（严重堵塞、轻微堵塞、无堵塞），发现初始孔隙率和水力梯度对堵塞过程具有显著调控作用。数值模拟表明堵塞演化存在四阶段特征（初始、渗透率下降、回升、稳定），该方法在保持孔隙尺度精度前提下有效解决了传统DEM-CFD方法网格粗化导致的精度损失问题，为颗粒介质流变研究提供了新的计算范式。

　　

张东梅|陈思睿|谢晓创|陈慧豪|陈兆庚

同济大学土木工程学院岩土工程系，上海200092，中国

摘要

渗流过程伴随着孔隙堵塞现象，局部颗粒的积累最终会影响整体颗粒的运输和渗流模式。然而，由于计算要求和精度限制，在数值模拟中研究渗流过程中的孔隙尺度堵塞和流体阻塞仍然具有挑战性。本研究开发了一种改进的孔隙尺度流体-颗粒耦合方法，该方法对粗颗粒和细颗粒与流体的相互作用采用了不同的耦合方案，从而避免了频繁的网格更新。该方法将颗粒堵塞模型和孔隙密度流动模型集成到离散元方法（DEM）-孔隙网络方法（PNM）框架中，以保持孔隙尺度的准确性，同时提高计算效率。通过具有不同初始孔隙度和水力梯度的模拟，分析了孔隙和元素尺度上的堵塞效应。结果表明，孔隙堵塞可以分为三种状态：严重堵塞、轻微堵塞和解除堵塞，这些状态的变化受孔隙喉部拱形形成和颗粒积累的控制。初始孔隙度和水力梯度显著影响孔隙堵塞。孔隙堵塞主要从试样底部开始，并通过依赖性堵塞传播，改变局部流场并抑制颗粒迁移。在元素尺度上，堵塞降低了孔隙喉部的渗透性，渗透性和流体速度表现出四个阶段的变化：初始阶段、下降阶段、上升阶段和稳定阶段。与传统DEM-CFD方法相比，所提出的方法在捕捉孔隙尺度细节方面表现出更高的能力，并具有更优的计算效率。这些发现扩展了我们对多尺度渗流机制的理解。

引言

渗流是颗粒系统内部侵蚀的一种典型形式，对地质结构（包括大坝、基础、地下管道和盾构隧道）构成重大威胁[1]、[2]、[3]。特别是在级配不均的砂土中——这种土壤常见于大坝基础、沿海冲积层或某些冰川河流沉积物中——细颗粒可以从粗颗粒骨架上脱落，并在渗流作用下容易被侵蚀。细颗粒的侵蚀可能导致沉洞和土壤空洞[4]、[5]，这会显著改变地质结构的力学特性并降低其使用性能[6]、[7]。因此，研究渗流的微观机制对于地质工程中的灾害预防和缓解至关重要。

渗流过程始终伴随着孔隙堵塞现象，这对渗流过程有重要影响[8]、[9]。例如，孔隙堵塞会阻碍流体运动和颗粒迁移，而孔隙的解除堵塞则可能导致瞬时侵蚀和土壤渗透性的增加[10]、[11]。然而，由于固体和流体相之间的复杂相互作用，仍然难以用各种方法对渗流过程中的孔隙尺度堵塞和解除堵塞进行建模和分析。过去几十年中进行了大量实验室测试来研究这些机制。Nguyen等人[12]使用X射线计算机断层扫描（CT）识别了渗流过程中的局部堵塞和微观结构变化。Liu等人[9]利用高精度摄影技术跟踪了孔隙尺度的堵塞演变过程，并研究了颗粒的堵塞机制以及流速对堵塞行为的影响。Wang等人[13]使用扫描电子显微镜（SEM）观察了水注入过程中的颗粒迁移和颗粒重新分布机制。这些实验增强了我们对孔隙堵塞的理解。然而，实验中观察孔隙尺度现象需要精密仪器，因此成本高昂且耗时较长。此外，孔隙堵塞包括颗粒堵塞和流体阻塞，但无法通过实验监测孔隙尺度上的流体运动，导致对渗流过程中孔隙尺度行为的动态原理了解不足。

适当的数值方法可以有效地模拟固体和流体相中的渗流过程。使用离散元方法（DEM）进行数值模拟已被证明是模拟渗流过程中颗粒运动和堵塞的有效方法[14]、[15]、[16]。通常，由于DEM的计算成本极高，渗流的模拟和分析通常在元素尺度上进行，方法是通过对工程尺度模型进行降尺度处理（如图1所示）。已经提出了几种流体-颗粒耦合方案来解决流场建模的挑战。作为最广泛使用的流体-颗粒耦合方法，DEM-CFD（计算流体动力学）已被用于研究渗流过程中的颗粒运输、力链演变和力学响应[17]、[18]、[19]、[20]。然而，由于颗粒必须适应流体元素才能求解修正后的纳维-斯托克斯方程，传统的DEM-CFD方法需要相对粗糙的CFD网格。因此，这种方法难以研究由于细颗粒运动引起的孔隙尺度流体流动和渗透性变化等细节[21]。这一限制在考虑孔隙尺度堵塞或解除堵塞对元素尺度渗流行为的影响时带来了显著挑战，特别是在局部水动力变化与整体渗流模式演变之间的多尺度耦合效应方面。

相比之下，完全解析的DEM-CFD方法和DEM-LBM（格子玻尔兹曼方法）能够在孔隙尺度上很好地捕捉流体-颗粒相互作用的物理过程。因此，这些方法已被用于对单个孔隙的详细分析[11]、[22]。虽然这些方法研究了流体流动和孔隙堵塞-解除堵塞机制，但渗流并不仅仅是一个局限于单个孔隙的孤立现象，而是源于颗粒在多孔介质中耦合传输过程中的复杂水动力相互作用（如图1所示的元素尺度）。然而，这些方法在计算效率上较低，因为需要复杂的流体动力学计算和精细的流体域离散化，因此在元素尺度上模拟渗流（尤其是在三维情况下）面临重大挑战。局部堵塞对整体渗流过程的影响尚不清楚。因此，需要一种新的方法，能够在三维模拟中平衡计算效率和处理孔隙流体，以分析具有堵塞的渗流过程。

此外，还出现了几种用于处理DEM中流体计算的数值算法，包括平滑粒子动力学（SPH）、浸没边界方法（IBM）和孔隙网络方法（PNM）。其中，PNM为孔隙尺度上的流体计算提供了优秀的框架，提高了计算效率并确保了结果的可靠性[23]。耦合的DEM-PNM方法比其他耦合方法更有效地表示了水力机械过程，现在可以处理多达350,000个颗粒的元素尺度问题[24]。Chareyre等人[23]开发了一种计算效率高的DEM-PNM模型，用于模拟应力作用下流体饱和颗粒材料的演变。他们的方法通过局部应用于颗粒堆积的规则Delaunay三角剖分来定义孔隙体及其相互连接。Wautier等人[15]随后使用DEM-PNM研究了应力状态和水力加载条件对颗粒运输起始的影响。Zhu等人[25]进一步将孔隙密度流动模型集成到DEM-PNM框架中，从而提高了渗流特性的准确性。此外，Gao等人[26]使用DEM-PNM研究了初始细颗粒含量不同的级配不均土壤的渗流后力学响应。然而，传统的DEM-PNM框架在模拟级配不均土壤中的动态渗流过程时存在显著限制[25]。粗颗粒和细颗粒在渗流过程中的运动行为差异会导致网络发生显著变形，需要频繁调整孔隙网络拓扑结构。这些固有的方法学限制导致Gao等人[26]仅模拟3.3秒的元素尺度渗流过程就需要15天的计算时间，凸显了该方法在动态场景下的高昂计算成本。

鉴于上述限制，本研究提出了一种基于DEM-PNM的改进孔隙尺度流体-颗粒耦合方案，用于研究级配不均土壤的渗流过程。该方法针对级配不均土壤的几何异质性和渗流过程中粗颗粒与细颗粒的差异运动行为，采用不同的耦合方案来处理流体与粗颗粒和细颗粒的相互作用，以提高计算效率并解决传统DEM-PNM的不足。这是通过引入孔隙密度流动模型和颗粒堵塞模型来实现的。流体方案使用MATLAB编程在DEM框架内实现[27]。利用这种方法，将进行更深入的研究，以阐明孔隙堵塞机制及其对元素尺度渗流过程的影响。首先描述并验证了所提出的方法，然后应用该方法模拟具有不同初始孔隙度和水力梯度的渗流过程。在结果部分，首先分析了孔隙尺度上的堵塞和解除堵塞，随后考虑不同情况进行了元素尺度分析。所提出的方法和发现可以加深我们对渗流过程中土壤物理机制和水力行为的理解。

流体-颗粒耦合方法方案

渗流本质上是细颗粒从粗颗粒骨架上脱落的过程。粗颗粒通常不会发生大规模迁移。本研究提出的方法旨在仅在粗颗粒之间构建PNM网络，并基于相变理论重新考虑细颗粒与流体之间的耦合，从而从根本上减少PNM网格的频繁更新需求。因此，需要解决以下挑战

介观参数校准

在本研究中，使用赫兹接触定律[35]模拟颗粒之间的法向接触，而弹簧简化了切向连接。根据现有关于砂土性质的研究[36]校准输入参数，以确保准确性。在50、100、200和300 kPa的法向应力下进行了数值直剪试验，以获得土壤的宏观参数，如图7所示。介观参数为

孔隙堵塞和解除堵塞的演变

在孔隙尺度上的监测可以提供对渗流机制的详细理解。在渗流过程中，细颗粒在不同的孔隙之间迁移。由于孔隙结构和细颗粒分布的差异，可能会发生不同程度的孔隙堵塞。为了量化堵塞程度，根据公式(18)开发了一个无量纲系数$C_p$$C_p$$k_{\mathrm{ei}}$$k_i$