2025-11-23T03:55:16.576493

A finite element method using a bounded auxiliary variable for solving the Richards equation

Benfanich, Bourgault, Beljadid
The Richards equation, a nonlinear elliptic parabolic equation, is widely used to model infiltration in porous media. We develop a finite element method for solving the Richards equation by introducing a new bounded auxiliary variable to eliminate unbounded terms in the weak formulation of the method. This formulation is discretized using a semi-implicit scheme and the resulting nonlinear system is solved using Newton's method. Our approach eliminates the need of regularization techniques and offers advantages in handling both dry and fully saturated zones. In the proposed techniques, a non-overlapping Schwarz domain decomposition method is used for modeling infiltration in layered soils. We apply the proposed method to solve the Richards equation using the Havercamp and van Genuchten models for the capillary pressure. Numerical experiments are performed to validate the proposed approach, including tests such as modeling flows in fibrous sheets where the initial medium is totally dry, two cases with fully saturated and dry regions, and an infiltration problem in layered soils. The numerical results demonstrate the stability and accuracy of the proposed numerical method. The numerical solutions remain positive in the presence of totally dry zones. The numerical investigations clearly demonstrated the capability of the proposed method to effectively predict the dynamics of flows in unsaturated soils.
academic

A finite element method using a bounded auxiliary variable for solving the Richards equation

基本信息

  • 论文ID: 2510.13012
  • 标题: A finite element method using a bounded auxiliary variable for solving the Richards equation
  • 作者: Abderrahmane Benfanich (University of Ottawa), Yves Bourgault (University of Ottawa), Abdelaziz Beljadid (University of Ottawa & Mohammed VI Polytechnic University)
  • 分类: math.NA cs.NA physics.comp-ph
  • 发表时间: 2025年10月14日 (arXiv预印本)
  • 论文链接: https://arxiv.org/abs/2510.13012

摘要

Richards方程是一个非线性椭圆抛物型方程，广泛用于多孔介质中的渗透建模。本文通过引入一个新的有界辅助变量来消除弱形式中的无界项，开发了求解Richards方程的有限元方法。该方法采用半隐式格式离散化，使用Newton方法求解非线性系统。该方法消除了正则化技术的需要，在处理干燥和完全饱和区域方面具有优势。在所提出的技术中，使用非重叠Schwarz区域分解方法来模拟分层土壤中的渗透。数值实验验证了所提方法的有效性，包括完全干燥纤维片中的流动建模、完全饱和和干燥区域的两种情况，以及分层土壤中的渗透问题。

研究背景与动机

问题定义

Richards方程描述多孔介质中的非饱和流动，是水文地质学、环境工程和土壤科学中的基础模型。该方程存在多种形式：

  1. 压力头形式(Ψ)：适用于完全饱和和部分非饱和流动，但在S→0时数值不稳定
  2. 饱和度形式(S)：在干燥土壤中表现更好，但在饱和区域(S=1)存在数值困难
  3. 混合形式(S,Ψ)：改善质量守恒特性，但增加了复杂性

现有方法的局限性

  1. 主变量切换方法：在饱和/非饱和界面处可能产生非物理解
  2. 正则化技术：需要额外的建模参数选择
  3. 传统有限元方法：在处理完全干燥(S=0)和完全饱和(S=1)区域时存在数值困难

研究动机

开发一种统一的数值方法，能够：

  • 避免变量切换和正则化
  • 稳定处理干燥和饱和区域
  • 保持良好的质量守恒特性
  • 适用于异质多孔介质

核心贡献

  1. 提出新的u-形式：引入有界辅助变量u，消除弱形式中的无界项
  2. 开发稳定的数值格式：结合半隐式时间离散和Newton迭代
  3. 实现区域分解方法：处理分层土壤等异质介质问题
  4. 验证方法有效性：通过多个数值实验证明方法的稳定性和精度

方法详解

任务定义

求解Richards方程： ∂θ∂t+∇⋅q=0\frac{\partial\theta}{\partial t} + \nabla \cdot q = 0

其中θ是含水量，q是根据Darcy-Buckingham定律的水流通量： q=−Ks(x)Kr(x,S)∇(Ψ(x,S)+z)q = -K_s(x)K_r(x,S)\nabla(\Psi(x,S) + z)

新变量u的构造

对于所有模型，Leverett函数的导数可写为： J′(S)=C(x)S−a(1−Sc)−bJ'(S) = C(x)S^{-a}(1-S^c)^{-b}

定义辅助变量： u=∫0S(1−sc)−bdsu = \int_0^S (1-s^c)^{-b} ds

使用不完全Beta函数表示： u=1cB(Sc,1c,1−b)u = \frac{1}{c}B(S^c, \frac{1}{c}, 1-b)

u-形式的Richards方程

在齐次介质中，u-形式为： ϕ∂S(u)∂t−∇⋅(KsKr(S(u))(hcapCS(u)−a∇u+ez))=0\phi \frac{\partial S(u)}{\partial t} - \nabla \cdot \left(K_sK_r(S(u))\left(h_{cap}CS(u)^{-a}\nabla u + e_z\right)\right) = 0

变分形式

寻找u ∈ L²(I;H¹(Ω))，使得： ∫Ωϕ∂S(u)∂tvdx+∫ΩKsKr(S(u))(hcapS(u)−a∇u+ez)⋅∇vdx+边界项=0\int_\Omega \phi \frac{\partial S(u)}{\partial t} v dx + \int_\Omega K_sK_r(S(u))\left(h_{cap}S(u)^{-a}\nabla u + e_z\right) \cdot \nabla v dx + \text{边界项} = 0

时间离散化

采用半隐式Euler格式： ∫ΩϕS(uhn+1)−S(uhn)Δtvhdx+∫ΩKsKr(S(uhn))(hcapS(uhn)−a∇uhn+1+ez)⋅∇vhdx=0\int_\Omega \phi \frac{S(u_h^{n+1}) - S(u_h^n)}{\Delta t} v_h dx + \int_\Omega K_sK_r(S(u_h^n))\left(h_{cap}S(u_h^n)^{-a}\nabla u_h^{n+1} + e_z\right) \cdot \nabla v_h dx = 0

Newton迭代

对于非线性项S(u^{n+1,m+1})，使用一阶Taylor展开： S(un+1,m+1)=S(un+1,m)+∂S∂u(un+1,m)(un+1,m+1−un+1,m)+O((un+1,m+1−un+1,m)2)S(u^{n+1,m+1}) = S(u^{n+1,m}) + \frac{\partial S}{\partial u}(u^{n+1,m})(u^{n+1,m+1} - u^{n+1,m}) + O((u^{n+1,m+1} - u^{n+1,m})^2)

区域分解方法

对于分层土壤，采用非重叠Schwarz方法：

  • 在每个子域独立求解Richards方程
  • 通过界面条件耦合：压力头连续性和流量守恒
  • 使用Robin型传输条件进行迭代求解

实验设置

数值实现

  • 软件平台：FreeFEM++和FENICSx
  • 线性求解器：UMFPACK, PETSc, MUMPS
  • 网格：一致三角网格，线性有限元(k=1)

评价指标

L²误差： Eh,Δt(T)2=∥u(⋅,T)−uh,Δt(⋅,T)∥L2(Ω)2E_{h,\Delta t}(T)^2 = \|u(\cdot,T) - u_{h,\Delta t}(\cdot,T)\|_{L^2(\Omega)}^2

收敛阶数：

  • 空间收敛阶：p=log⁡2(Eh,Δt(T)/Eh/2,Δt(T))p = \log_2(E_{h,\Delta t}(T)/E_{h/2,\Delta t}(T))
  • 时间收敛阶：q=log⁡2(Eh,Δt(T)/Eh,Δt/2(T))q = \log_2(E_{h,\Delta t}(T)/E_{h,\Delta t/2}(T))

测试案例

  1. 1D完全非饱和纤维片：验证干燥条件下的稳定性
  2. 1D饱和-非饱和混合区域：测试处理不连续初值的能力
  3. 制造解验证：计算收敛阶数
  4. 2D完全饱和流动：与Hydrus软件对比
  5. 分层土壤问题：验证区域分解方法
  6. 异质介质包含体：测试复杂几何的处理能力

实验结果

主要结果

1D纤维片测试

  • 网格：500个均匀网格点
  • 时间步长：Δt = 1s
  • 结果：u-形式与S-形式结果完全一致，与实验数据吻合良好

饱和-非饱和混合区域

  • 网格：5000个均匀网格点
  • 时间步长：Δt = 10⁻⁵ 小时
  • Newton收敛：平均1-5次迭代
  • 稳定性：解保持单调性，无数值振荡

制造解收敛性分析

空间收敛性（固定Δt = 5×10⁻⁵）：

hL²误差收敛阶p计算时间(s)
1.43×10⁻²6.04×10⁻³-25.1
7.14×10⁻³8.99×10⁻⁴2.7547.2
3.57×10⁻³1.39×10⁻⁴2.6982.8
1.79×10⁻³3.48×10⁻⁵2.00151.0

时间收敛性（固定h = 1×10⁻³）：

ΔtL²误差收敛阶q计算时间(s)
1×10⁻²8.23×10⁻³-6.9
5×10⁻³3.90×10⁻³1.0811.1
2.5×10⁻³1.92×10⁻³1.0219.3
1.25×10⁻³9.60×10⁻⁴1.0033.3

2D饱和流动对比

  • 网格：95,142个三角单元，47,972个顶点
  • 时间步长：Δt = 0.001 小时
  • 结果：与Hydrus软件结果高度吻合

分层土壤测试

  • 网格：89,306个三角单元，45,064个顶点
  • 时间步长：Δt = 0.05 小时
  • 收敛：容差10⁻³，最多10次迭代，λ = 25
  • 结果：与文献20中的结果一致

方法对比

在异质介质包含体测试中：

  • u-形式：ε=0时仍稳定收敛，平均4次迭代
  • Ψ-形式：ε=10⁻⁵时需要>1000次迭代，残差容差仍为0.04

相关工作

主要研究方向

  1. 变量切换方法：Diersch & Perrochet (1999), Forsyth et al. (1995)
  2. 正则化技术：Schweizer (2007), Pop & Schweizer (2011)
  3. 混合形式：Celia et al. (1990), Maina & Ackerer (2017)
  4. 区域分解：Lions (1990), Manzini & Ferraris (2004)

本文优势

  1. 统一处理：无需变量切换或正则化
  2. 数值稳定：在S∈0,1全范围内稳定
  3. 物理一致：保持压力头连续性
  4. 计算效率：Newton快速收敛

结论与讨论

主要结论

  1. u-形式有效性：成功消除了传统形式中的无界项
  2. 数值稳定性：在完全干燥和饱和区域都保持稳定
  3. 收敛性能：空间二阶收敛，时间一阶收敛
  4. 实用性：适用于复杂的异质介质问题

局限性

  1. 模型依赖：需要满足条件lim⁡S→0K(S)S−a<∞\lim_{S\to 0} K(S)S^{-a} < \infty
  2. 非线性求解：饱和区域需要Newton迭代
  3. 界面处理：异质介质需要区域分解技术
  4. 参数敏感性：区域分解中的Robin参数λ需要调节

未来方向

  1. 理论分析：严格的数值分析和收敛性证明
  2. 多物理耦合：扩展到溶质和热传输
  3. 自适应技术：开发自适应网格和时间步长
  4. 并行计算：大规模并行实现

深度评价

优点

  1. 创新性强：u-变量的构造巧妙，理论基础扎实
  2. 数值稳定：避免了传统方法的数值困难
  3. 实验充分：涵盖了从1D到2D，从均质到异质的多种情况
  4. 工程实用：与商业软件Hydrus的对比增强了可信度

不足

  1. 理论分析不足：缺乏严格的收敛性和稳定性证明
  2. 计算复杂度：Newton迭代的计算成本分析不够详细
  3. 参数选择：Robin参数λ的选择缺乏理论指导
  4. 三维验证：缺乏三维问题的验证

影响力

  1. 学术价值：为Richards方程数值求解提供了新思路
  2. 实用价值：在水文地质和土壤科学中有广泛应用前景
  3. 可复现性：提供了详细的实现细节，便于复现

适用场景

  1. 多孔介质流动：地下水流动、土壤水分运移
  2. 工程应用：大坝渗流、边坡稳定性分析
  3. 环境问题：污染物迁移、地下水修复
  4. 材料科学：纤维材料中的流体渗透

参考文献

关键文献

  1. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1(5), 318-333.
  2. Celia, M.A., et al. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26(7), 1483-1496.
  3. van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892-898.
  4. Manzini, G., Ferraris, S. (2004). Mass-conservative finite volume methods on 2-d unstructured grids for the Richards' equation. Adv. Water Resour. 27, 1199-1215.

总体评价：这是一篇在Richards方程数值求解领域具有重要创新意义的论文。通过引入有界辅助变量u，巧妙地解决了传统方法在处理干燥和饱和区域时的数值困难。方法具有良好的理论基础和实用价值，数值实验充分验证了方法的有效性。虽然在理论分析方面还有待加强，但整体上是一项高质量的研究工作。