基于MLS-NMM方法的饱和土壤中浆液渗流的数值模拟

《Computers & Mathematics with Applications》：Numerical simulation of slurry infiltration in saturated soil based on MLS-NMM

　　摘要：本文在MLS-NMM框架内实现了一个现有的浆体渗透模型，该模型考虑了浆体颗粒传输与水力特性的耦合。首先，基于控制方程和边界条件构建问题的弱形式。然后，利用移动最小二乘法（MLS）中节点的影响域作为数学贴片，来构建数值流形方法（NMM）的数学覆盖。接着，采用两步向后有限差分方法进行时间离散化，从而建立了浆体渗透问题的MLS-NMM离散化模型。采用增量法和迭代法求解非线性控制方程。将数值结果与现有的一维实验室测试结果进行比较，验证了理论模型的合理性，并证明所提出的MLS-NMM结合了MLS和NMM的优势。最后，将所提出的方法应用于浆体沟槽的工程模型，成功地再现了流体压力和颗粒浓度的时空分布特性及演变趋势。

引言：自上个世纪以来，许多学者对浆体在多孔介质中的渗透行为进行了广泛研究[1,2]。浆体渗入地层的过程是一个复杂的过程，受多种因素的影响，主要包括浆体渗透速度、孔隙水压力、浆体质量浓度、浆体粘度和介质孔隙度等关键参数。在浆体通过多孔介质的渗透过程中，浆体所携带的颗粒会与介质孔隙内的滞留位置接触，颗粒的逐渐沉积会导致介质孔隙结构的物理堵塞[3]。受此影响，多孔介质的孔隙空间逐渐被占据，进而导致渗透区的渗透性显著降低[4]。在典型的土木工程场景中，浆体在饱和沙体中的渗透对施工稳定性有着极其重要的影响[5,6]。因此，与浆体颗粒沉积和渗透相关的问题已成为工程师关注的核心问题。准确量化这些物理量有助于提高工程设计施工的科学性和针对性。在形成致密且不透水的泥饼之前，浆体压力无法有效地作用于土壤骨架。当浆体压力大于地层的静水孔隙水压力时，由此产生的压力差会驱动浆体中的细小颗粒渗入土壤，导致土壤孔隙水压力持续增加。在浆体渗透过程中，携带沙子的颗粒逐渐沉积在土壤颗粒之间的孔隙通道中，最终在土壤表面形成由浆体颗粒和沙粒组成的致密防水混合层[7,8]。过滤饼形成后，渗入的浆体压力可以有效平衡作用在挖掘地层上的土压力和水压力[9,10]。当浆体的性质与地层参数不匹配时，浆体颗粒可能通过地层的孔隙通道，未能形成致密、低渗透性的过滤饼。这会导致浆体在地下形成大量通道，阻碍额外的推力传递到隧道掘进面，从而显著增加隧道施工事故的概率，包括掘进面不稳定、周围地面沉降不均、地面隆起过度甚至水涌入等危险[11][12][13][14]。基于此，对浆体渗透的复杂耦合作用机制及其核心影响因素进行系统研究（如图1所示）具有重要的工程实践价值和理论探索意义。研究发现，当浆体颗粒沿孔隙通道注入沙层时，其迁移模式可以分为对流（在孔隙水动力学的作用下，溶质整体随流体移动的过程）和扩散（迁移过程中的溶质空间扩散现象），最终表现为两种运动状态：在孔隙内沉积或继续通过沙介质迁移。

随着与浆体渗透相关的工程实践的不断进步和深化，理论分析、实验室测试和数值模拟已成为近年来科学领域对该过程进行系统研究的主流技术方法。传统的经验方法往往无法实现上述基本物理量的定量表征，因为它们不遵循严格的物理定律和数学推导逻辑。相比之下，理论研究方法在处理多尺度工程问题和实现变量定量分析方面具有独特优势。为了阐明浆体渗透过程中各种参数的演变规律，Moghadasi等人[15]假设渗透颗粒在模型构建过程中能够直接填充介质孔隙。然而，这项研究有两个局限性：首先，它没有考虑颗粒在多孔介质内的滞留现象；其次，它没有考虑颗粒堵塞行为对渗透场的干扰效应，导致模型与实际工程条件之间存在偏差。Reddi等人[16]将土壤介质简化为若干长度相等、直径不同的平行毛细管束，并基于Kozeny水力半径模型建立了土壤渗透系数与孔隙度之间的相关方程。然而，这样的物理堵塞模型无法与实验结果达到定量一致。关于浆体渗透的耦合效应问题，Alem等人[17]从地层过滤效应与颗粒堵塞行为之间的耦合关系入手，建立了过滤损失、孔隙度和渗透系数之间的定量表征方程，并分析了浆体流量对深层渗透过程的影响特性。然而，该模型在应用范围上存在局限性，无法系统解释渗透过程中多个变量的协同演变规律。

实验室实验是研究浆体渗透参数影响因素的可靠方法。许多学者使用圆柱形沙柱测试系统进行了实验室测试，旨在深入分析浆体在土壤中的传输特性和演变行为[18][19][20][21]。实验室测试的成果加深了人们对浆体渗透物理过程的认识。测试表明，渗透沙柱的水力特性、浆体质量浓度、浆体粘度参数和工作压力在决定浆体渗透和沉积模式方面起着决定性作用[12,13]。在采用浆体方法的浆体沟槽建设等土木工程项目中，所使用的浆体颗粒尺寸通常大于75微米[5,22]。Herzig等人[2]的研究表明，在渗透过程中，大颗粒（颗粒尺寸大于30微米）与胶体颗粒之间的核心差异体现在两个方面：堵塞机制和悬浮对流的影响。胶体颗粒的分子吸附起着关键作用，而浆体颗粒的沉积模式与分子吸附有本质不同。基于实验室测试结果，该研究团队总结了不同大小胶体颗粒的沉积模式，并建立了悬浮渗透过程的半解析半经验模型。在模型的控制方程中，经验系数被定义为阻力函数的函数。Bradford等人[23]发现了实验测量的沉积系数与各种物理参数之间的复杂相关性，并通过引入附加影响因素（如孔隙尺寸分布和颗粒尺寸效应）进一步修订和优化了控制方程。实验室模型测试可以高精度地模拟特定土壤的灌浆过程。然而，这种方法不仅需要大量的时间和经济成本，而且用于测试的土壤样本与实际原位土壤层存在一定的偏差。为了解析浆体渗透的微观机制以及一些变量在实验室测试中难以测量的问题，许多研究人员进行了浆体渗透的数值模拟工作。目前广泛使用的数值模拟方法主要包括有限元方法（FEM）、有限差分方法（FDM）、有限体积方法（FVM）和离散元方法（DEM）。Zhou等人[24]建立了一种计算饱和沙中浆体渗透的方法，该方法耦合了颗粒扩散、浆体渗透和土壤变形的多场相互作用，并采用FEM完成了数值解。然而，该模型的局限性在于没有考虑浆体颗粒的粘度效应，而粘度效应是关键影响因素之一。Huang等人[8]基于液相连续性条件修订和改进了现有模型，以解决质量守恒问题。Lu等人[25]和Liu等人[26]采用FDM对一维浆体渗透模型进行了数值求解。这种方法涉及的变量参数较为全面，但在处理不规则二维渗透问题时存在显著局限性。DEM在模拟浆体渗透方面表现出良好的能力[27][28][29][30]。其优势在于能够准确再现颗粒介质的微观结构特征和颗粒迁移的演变过程，但该方法对可模拟的颗粒类型、几何形状和数量有严格限制，从而在工程规模的应用中面临相当大的挑战。Liu等人[31]采用计算流体动力学-离散元方法（CFD-DEM）耦合方法，系统研究了不同类型浆体的渗透特性和主要控制条件，并基于浆体颗粒尺寸分布与地层颗粒之间的匹配关系，对不同类型浆体的渗透规律进行了深入分析。

NMM最初由 Shi Genhua 博士于1991年提出[32]。该方法基于Galerkin变分原理，创新性地结合了FEM和离散变形分析（DDA）的技术优势。通过建立两个相互独立的覆盖系统（即数学覆盖和物理覆盖），NMM具备了表征解域内不连续性的内在能力，从而能够统一解决连续和不连续的力学问题。移动最小二乘法（MLS）由Lancaster和Salkauskas于1981年提出[33]。作为无网格方法中的核心局部加权拟合算法，它被广泛用于构建场近似。与全局最小二乘法不同，MLS定义了一个以评估点为中心的紧凑支持域，并仅使用相邻的离散节点构建局部多项式基函数。MLS已被广泛应用于Lotka单性种群模型的数值模拟[34]、基于神经纤维位置的时滞时空神经场方程[35]、一类非矩形二维域上的非线性时变混合积分方程[36]以及由空气污染模型产生的延迟积分方程[37]。在NMM的传统应用中，数学覆盖的构建通常依赖于三角形或四边形有限元网格，这种NMM变体被称为基于有限元的数值流形方法（FE-NMM）。为了规避复杂的预处理步骤并提高计算效率，Zheng等人[38]提出了一种基于MLS的NMM，简称MLS-NMM。该方法利用MLS节点的形状函数作为数学覆盖贴片的权重函数，其中数学覆盖贴片由对应于MLS节点的矩形或圆形影响域构成。与FE-NMM相比，MLS-NMM无需生成流形元素，具有更高效的建模特性。MLS-NMM已被应用于各种工程和力学问题的数值模拟研究，包括裂纹传播和演变[38][39][40]、无限制渗流[41,42]和热传导[43,44]等。

与现有研究相比，本文在MLS-NMM框架内实现了饱和沙体的两相渗透理论，考虑了浆体渗流和颗粒传输的耦合行为。与传统单场MLS-NMM公式不同，所提出的方法完全耦合了浆体流动和颗粒沉积过程，引入了由颗粒沉积引起的孔隙度和渗透性的动态演变，并实现了不可压缩浆体渗透过程中压力场和浓度场的迭代求解。该方法基于颗粒质量守恒定律，明确考虑了沉积和悬浮颗粒迁移的动态特性，推导出一组非线性控制方程，将颗粒扩散方程与渗流速度场耦合起来，从而系统地揭示了浆体渗透和沉积的耦合作用机制。通过构建控制方程的惩罚弱形式，建立了一个用于浆体渗透问题的离散MLS-NMM模型。该模型能够有效描述浆体渗透过程中孔隙水压力、介质孔隙度、渗流速度和颗粒沉积速率的动态演变规律。通过与现有的一维实验室测试结果比较，验证了所提理论模型的可靠性。最后，将这种方法应用于一个浆体沟渠工程案例的数值模拟中，再现了颗粒沉积量、流体压力和颗粒浓度的时空分布特征及演变趋势，进一步证明了本文提出方法的正确性和工程适用性。

**浆体渗透模型的建立**
浆体渗透的物理机制主要包含两个方面：一是浆体在多孔介质中的渗流过程，二是浆体颗粒的迁移与传输过程。相关物理机制的示意图如图2所示。

**MLS-NMM简介**
由于NMM能够统一处理连续和非连续问题，因此具有广泛的应用前景。该方法中结合了两种类型的“覆盖”：数学覆盖和物理覆盖。在大多数应用场景中，NMM的数学覆盖通常由规则几何形状的有限元构建。对于二维平面问题，常用的单元形式主要是三角形和四边形。这种构建方法的核心优势在于变分原理。

**浆体渗透模拟的控制方程**
根据第2节所述，浆体渗透模拟的相关控制方程如下：
$$
\left\{
\begin{array}{ll}
\mathcal{D}\left(\vec{u}+\mathbf{c}_p\right), & 0 \\
\vec{q} + (\mathbf{v}_p - \mathbf{v}_s), & 0 \\
\mathbf{P}_s - \mathbf{P}_e, & 0 \\
\end{array}
\right.
$$
其中：
$$
\mathcal{D}(\vec{u}+\mathbf{c}_p) = \frac{\partial}{\partial t}(\vec{u}+\mathbf{c}_p) - \nabla \cdot (\mathbf{F} \cdot \vec{u}+\mathbf{g})
$$
$$
\mathbf{q} = \nabla \cdot \vec{q} + \mathbf{S}
$$
$$
\mathbf{P}_s = \int_{\Sigma} p d\mathbf{v} \, dS
$$
$$
\mathbf{P}_e = \int_{\Sigma} p \, d\mathbf{u} \, dS
$$
$$
\mathbf{F} = \rho \nabla \mathbf{v}_p + \mathbf{g}
$$
$$
\mathbf{g} = \rho \mathbf{a} \cdot \vec{u}
$$

**一维浆体渗透过程的数值建模**
一维测试是研究浆体渗透机制的典型方法，可用于验证所提方法的合理性。本部分旨在评估所建立控制方程在捕捉和再现Lu等人[19,25]发表的实验结果方面的有效性。

**讨论**
本研究重点关注饱和多孔介质中浆体渗流和颗粒扩散现象，主要需要解决的物理过程包括浆体渗流和颗粒扩散。所提模型的一个关键假设是，浆体颗粒的渗透和沉积完全由物理行为主导，未考虑颗粒与天然有机物或聚合物相互作用对沉积速率的影响。

**总结**
本文基于MLS-NMM框架，提出了一种用于饱和砂中浆体渗透的双相场耦合数学模型的数值分析方法。首先，根据浆体颗粒交换的质量守恒原理，推导出了与渗流速度耦合的颗粒扩散控制方程。该方程能够动态描述浆体流动与颗粒迁移的耦合过程，准确反映了对流和扩散现象。

**作者信息**
潘志伟（Zhiwei Pan）| 郑宏（Hong Zheng）
教育部城市安全与灾害工程重点实验室，北京工业大学，北京100124，中国