用于三维饱和-非饱和渗流分析与自由表面追踪的自适应Picard-SBFEM水平集方法

《Journal of Hydrology》：An adaptive Picard-SBFEM-level set formulation for three-dimensional saturated-unsaturated seepage flow analysis and free-surface tracking

　　

•用于三维非饱和渗流的统一Picard-SBFEM级联框架
•有界Aitken松弛技术稳定了非线性Picard迭代
•通过水平集方法实现自由面的隐式表达
•计算效率高，自由度少且迭代次数少

**引言**

具有移动自由面的三维非饱和渗流是水文、岩土工程和环境工程中的基本问题，广泛控制着降雨入渗、边坡稳定性、大坝安全评估以及地下水-地表水相互作用等一系列自然和工程过程（Jiang et al., 2023, Biniyaz et al., 2022, Dai and Lei, 2025, Huang et al., 2021, Li et al., 2026）。非饱和多孔介质的水力行为本质上具有非线性特征，因为水力传导系数和水保持特性强烈依赖于压力水头，这极大地增加了控制流方程数值求解的复杂性。此外，自由面的位置通常无法事先确定，且随水力状态演变，由此引入了一个高度耦合的移动边界问题。当扩展到三维时，计算负担、收敛鲁棒性和精度要求变得更加严苛，尤其是对于非均质区域和复杂几何形状。因此，开发一个高效的、稳定的且精确的数值框架，能够同时解决强非线性问题并可靠地追踪三维非饱和渗流中的自由面，仍然是一项重要且具有挑战性的任务。

土壤中的非饱和流通常在多孔介质水力学的理论框架内进行描述，其中Darcy型流假设和Richards方程构成了模拟稳态和非稳态渗流过程的基础。这些公式中的一个关键本构要素是土水特征曲线（SWCC），它表征了基质吸力、体积含水量和水力传导系数之间高度非线性的耦合关系。在过去的几十年中，人们提出了大量本构模型来表示这一关系，例如Gardner（1958）提出的指数模型、Brooks-Corey模型（Brooks and Corey, 1964）、van Genuchten-Mualem公式（Van Genuchten et al., 1992）以及Fredlund-Xing模型（Fredlund and Xing, 1994）。尽管这些模型在数学结构、参数化和适用吸力范围方面有所不同，但它们都在水力传导系数与吸力或有效饱和度之间建立了函数联系，从而为非饱和渗流的分析发展和数值模拟形成了一个通用的理论基础。

在这一理论框架的基础上，早期研究致力于提高非饱和流预测的准确性。Papagianakis和Fredlund（1984）提出了一种稳态饱和-非饱和流模型，并证明考虑非饱和区中可变的水力传导系数比仅考虑饱和模型能提供更准确的预测。后续研究探讨了参数不确定性的影响，揭示出水力传导系数对渗流变异具有主导影响，而SWCC参数通常起次要作用（Calamak and Yanmaz, 2017, Nam et al., 2010, Pourbeirak et al., 2026, Pham et al., 2018, Bakhshi et al., 2023）。进一步的研究强调，SWCC模型的选择强烈影响非稳态渗流下的边坡稳定性预测，表明实验室拟合准则应与全面的非稳态分析互补。数值方法的进步使得非饱和流的模拟更加鲁棒和精确。Chen等人（2024）开发了一种数值流形法（NMM）结合后向Euler离散化方法，实现了稳定且质量守恒的解。Sun等人（2024）应用虚拟元法对Richards方程进行离散化，证明了该方法在具有复杂网格几何的均质和分层土壤中的准确解。Motaman等人（2018）和Krifa和Chapuis（2025）强调了时间离散化的重要性，表明小时间步长对于在非稳态非饱和流中实现物理上合理且数值上收敛的结果至关重要。这些建模进展已被广泛应用于工程领域。Tan等人（2017）和Jassam和Abdulrazzaq（2019）研究了土坝，证明水力参数和空间非均质性均显著影响流量和自由面位置。Yu等人（2023）研究了黄土填筑边坡内部界面的非饱和渗流，表明局部积水和孔隙动力学主导了浅层破坏。Hou等人（2022）开发了一个解析-数值框架来量化非饱和流和饱和流对隧道掌子面稳定性的贡献，发现当上覆非饱和区较厚时，非饱和流主要影响剪切强度，而饱和流主导孔隙压力分布。Li和Li（2019）表明裂隙和排水孔强烈影响岩体中的渗流模式和潜水面，而Ma等人（2022）采用耦合SPH框架模拟了非饱和/饱和多孔介质中的多相相互作用、侵蚀和大变形。Yu等人（2025）采用离散-连续混合方法研究降雨诱发滑坡，强调了软弱夹层和非线性入渗在破坏起始中的作用。

与二维或简化分析相比，理解三维流动能够更真实地评估水力响应和破坏机制。在此背景下，Boger（1998）开发了边坡三维非饱和流的解析解，强调了土壤各向异性和导水率方向对渗流模式的影响，为最小化入渗提供了见解。Li和Yang（2020）应用运动极限分析结合吸力-应力Mohr-Coulomb准则来评估非饱和条件下三维主动土压力，证明三维流动和非饱和效应均显著改变了压力分布。Chen和Zhang（2006）利用三维饱和-非饱和渗流分析研究了古后坝的水力入渗，表明与二维模型相比，三维效应、分层和边坡边界会加速渗流并增加破坏风险。在此基础上，Li和Yang（2020）将非稳态非饱和流纳入三维分析，揭示在三维和非稳态效应共同作用下主动土压力显著降低，其变化取决于墙型几何和土壤类型。Ng等人（2022）采用具有旋转各向异性随机场的三维随机有限元分析来评估降雨下的边坡可靠性，表明三维效应强烈影响稳定性，尤其在横层面边坡中尤为明显。Zhai等人（2021）将二维解析解推广到三维梯形堤坝，并通过数值和实验结果验证了该方法的适用性。Di等人（2022）开发了盾构隧道掌子面三维渗流理论模型，证明其能准确预测水头且计算时间减少。Wei等人（2021）引入具有有限元实现的等效裂隙网络模型来研究具有自由面的稳态渗流，捕捉了裂隙和非均质性对流动行为的影响。从更广泛的计算力学角度来看，近期发展的开裂单元法和非局部损伤公式为处理局部不连续性和演化内部特征提供了有用的参考。例如，Wang等人（2025）、Zhang和Mang（2020）、Zhang等人（2021）、Zhang等人（2020）开发了基于Galerkin的开裂单元方法，用于处理断裂和不连续性演化，无需重网格划分、增强或复杂的裂缝追踪，而非局部近场动力学公式进一步改善了演化局部损伤的处理，并减少了网格依赖性（Zhang et al., 2021）。尽管与渗流在物理机制上不同，但这些研究为稳定表示演化界面和局部过渡区域提供了方法论上的见解。

与Richards型饱和-非饱和公式不同，经典的饱和自由边界渗流通常将潜水面视为饱和域的一个未知边界。因此，人们已给予大量关注来准确确定无压渗流中的自由面，由于其位置未知且可能具有时间依赖性，这仍然是一个关键挑战。Jiang等人（2010）采用三维NMM模拟无压渗流，证明了无需重网格划分即可精确确定自由面位置和计算渗流力，提高了相对于传统FEM的计算效率。Jiang等人（2014）进一步基于抛物型变分不等式（PVI）的有限元程序研究了裂隙网络中的非稳态自由面渗流，强调了裂隙密度和方向对自由面演化的影响。Yuan和Zhong（2016）引入了一种弱形式数值积分单元方法，采用"自由点"插值来模拟土坝中的三维无压渗流，以较低的计算成本实现了高效且准确的自由面预测。Sharif等人（2001）将稳定化有限元公式与水平集方法相结合来模拟土石坝中的自由面渗流，在二维和三维情况下均能准确捕捉界面，展示了大规模应用潜力。Bardet和Tobita（2002）提出了一种实用的无压渗流有限差分方法，提供了一种易于实现的自由面确定方案，具有教育和工程意义。Rafiezadeh和Ataie-Ashtiani（2014）开发了一种边界元方法，结合广义运动边界条件来处理各向异性介质中的三维非稳态自由面渗流，强调了各向异性对潜水面形状的关键影响。Wang等人（2014）基于能量-功原理建立了无压渗流的NMM方法，能高效确定潜水面并适应复杂几何和边界条件，即使采用粗网格也能实现高精度和更快的收敛。自由面也可被视为由水力场和边界约束控制的演化内部界面。在相关的耦合输运问题中，Weber等人（2025）强调了加热混凝土中渗透率表征对孔隙水压力预测的重要性。此外，代理模型和数据辅助策略，如Zhang等人（2023、2022）和Wang等人（2024）开发的基于图像和神经网络方法，已显示出加速重复场预测的潜力。这些研究表明，基于可靠的高保真数值数据集，可能为快速参数化或重复三维渗流分析提供了未来方向。

总体而言，这些研究表明三维非饱和渗流建模、演化界面处理和自由面预测取得了显著进展；然而，现有数值方法常存在收敛缓慢、自由面追踪精度有限和计算成本高企等问题，突显了对更高效和鲁棒求解策略的需求。

尺度边界有限元法（SBFEM）由Song和Wolf（1997）提出，是一种半解析的、以边界为重点的技术，通过单一缩放中心利用边界来建模整个域。通过在边界上集中离散化，该方法实现了增强的计算效率和精确解析局部特征（如角落、界面和陡峭梯度）的能力。与传统FEM相比，SBFEM将空间维度降低一维，能够高效建模无限或半无限域，并允许对域内部进行解析表示。与传统BEM不同，它避免了奇异积分，并能轻松处理各向异性材料。在过去几十年中，SBFEM已成功应用于结构力学（Chen et al., 2019, Ye et al., 2026, Huang et al., 2026, Huang et al., 2025）、弹动力学（Liu et al., 2025, Nie et al., 2025, Zou et al., 2019）、断裂分析（Yu et al., 2025, Natarajan and Ooi, 2024, Du et al., 2023）、岩土力学（Zhang et al., 2026, Gan et al., 2025, Tan et al., 2025, Han et al., 2024）和流体力学（Gao et al., 2025, Zang et al., 2024, Ye et al., xxxx）等问题，在处理复杂几何、捕捉奇异场和减少网格划分工作方面展现了持续的优势。在其基础公式的基础上，SBFEM已逐步扩展应用于日益复杂的渗流问题。Bazyar和Talebi（2015）将SBFEM应用于分区和各向异性土壤，证明仅边界离散化即可准确捕捉局部渗流特征。Liu等人（2018）引入了一种基于NURBS的等几何SBFEM用于多层域的渗流分析，改善了边界表示和解的精度。Huang等人（2025、2025）进一步推进了SBFEM在半空间和半无限域渗流中的应用，增强了数值稳定性和计算效率。在此基础上，Yi等人（2026）提出了一种自适应SBFEM用于渗流问题，将polytree型细化与嵌入式积分点和误差控制的精度相结合，克服了传统边界积分和重网格划分的局限性。Yan等人（2026）开发了一种多面体SBFEM，采用Wachspress形函数和多面体/八叉树网格，能够高效准确地模拟三维渗流（包括自由面流动），且无需重网格划分。

鉴于上述方法的广泛范围，表1总结了它们的方法论特征，并阐明了它们与本研究的交叉和差异。比较表明，Richards型公式主要贡献于非线性水力闭合的处理，FEM/NMM/水平集方法主要推进了自由面或界面更新，现有SBFEM研究主要改善了空间离散化和网格划分效率。本研究的方法论贡献定位在SBFEM、自适应Picard线性化和水平集追踪的耦合，统一于三维饱和-非饱和渗流公式框架中，其中非线性水力场和演化潜水面以相互耦合的方式求解。在该框架中，SBFEM和水平集方法发挥互补但相互作用的角色。

SBFEM提供了以边界为核心的空间离散化以及缩放方向上的半解析表示，而水平集函数则在固定离散化网格中隐式地表征未知的自由浸润面。对于给定的自由面构型，压力水头相关的渗透系数被纳入SBFEM系数矩阵中，非线性渗流场通过自适应Picard迭代进行求解。由此产生的渗流残差随后用于修正水平集函数，更新后的界面改变局部饱和–非饱和状态和进入下一轮SBFEM求解的渗透系数。这种嵌套的渗流-几何交互使得渗流场与自由面位置能够一致收敛，而无需显式界面重构或局部重新划分网格。

因此，所提出的耦合方法的主要优势在于将SBFEM的建模能力扩展至具有未知且演化的浸润面的三维强非线性饱和–非饱和渗流问题。这一扩展涉及SBFEM框架内可处理的非线性自由边界问题的范围，而非对基础物理控制模型的修改。同时，以边界为核心的SBFEM离散化方法和固定网格水平集处理方法通过减少空间离散化需求并避免自由面演化过程中的重复网格重构，保留了额外的计算优势。因此，所提出的方法在单一耦合求解框架内结合了增强的问题求解能力与计算效率。

本文其余部分组织如下。第2节介绍了用于处理控制方程强非线性的自适应线性化Picard迭代理论。第3节发展了Picard-SBFEM方法，包括问题陈述、SBFEM离散化、非线性控制方程及相应的求解过程。第4节引入了基于水平集的自由面追踪方法及其与Picard-SBFEM求解器的耦合。第5节提供了若干三维数值算例以验证所提方法的精度、鲁棒性和计算效率，并进行了详细讨论。最后，第6节总结本文的主要发现及未来工作展望。

## 章节摘要

### 自适应线性化Picard迭代理论

Picard迭代是一种求解离散化控制方程所产生的非线性代数方程的经典方法。它特别适用于非线性来源于状态相关材料特性的问题，因为它通过使用前一次迭代的解来评估非线性项从而线性化系统（Huang等人，2026；Lott等人，2026；Saeed，2017；Han等人，2026）。非线性算子方程在渗流分析中经常出现。

### 数学问题陈述

多孔介质域$V$中的非饱和渗流由Richards方程控制，该方程由非饱和条件下质量守恒与Darcy定律的组合推导得出。在本方法中，假设孔隙水不可压缩且多孔骨架为刚性；因此不考虑固体变形和孔隙弹性耦合，并且由于分析仅限于稳态渗流，与流体或固体压缩性相关的储水效应不予考虑。

### 基于耦合Picard-SBFEM-水平集方法的自由面追踪

水平集方法提供了一种隐式且几何一致的框架，用于在固定空间描述中表示和更新浸润面，如图3所示。该隐式表示避免了显式界面重构，自然地适应复杂几何形状和拓扑变化，并在固定网格上保持数值鲁棒性。在本饱和–非饱和Richards方程中，浸润面被表示为界面。

### 矩形断面坝中的渗流

为了检验所提方法对具有空间变化渗透性的无压渗流问题的适用性，考虑了矩形坝中的无压渗流问题，如图4所示。坝高10米，坝顶长度9米。上游边界承受对应水库水位的10米恒定水头，下游水位维持在2米。除…外，其余边界假定为不透水边界。

### 结论

本研究发展了一种用于三维稳态饱和–非饱和渗流分析及自由面追踪的半解析数值框架。该方法结合了SBFEM、自适应Picard迭代和水平集修正，用于求解具有压力水头相关渗透系数和未知自由面的非线性渗流问题。基于理论公式推导和数值验证，可得出以下结论。

1. 不同情况下的精度

### CRediT作者贡献声明

胡超豪：撰写——初稿，软件，方法论。刘军：撰写——审阅与编辑，概念化。王福明：撰写——审阅与编辑，方法论。郭成超：调查，正式分析。甘磊：调查，概念化。叶文彬：软件。薛冰涵：正式分析。林高：监督指导。

### 利益冲突声明

作者声明不存在任何已知的可能影响本文所述工作的竞争性经济利益或个人关系。

### 致谢

本研究得到国家自然科学基金委长江流域水科学研究联合基金重点项目（资助号U2340231）和国家自然科学基金委（资助号52309149）的资助，作者对此表示诚挚感谢。

胡超豪|刘军|王福明|郭成超|甘磊|叶文彬|薛冰涵|林高