一种用于求解非饱和土中Richards方程的h自适应局部离散Galerkin方法

《Journal of Computational and Applied Mathematics》：An h-adaptive local discontinuous Galerkin method for solving Richards’ equation in unsaturated soils

　　摘要：本文将h自适应局部离散Galerkin方法应用于基于Richards方程的非饱和渗透流问题。针对LDG离散化方法，给出了稳定性分析与误差估计。当水渗透到初始干燥的土壤中时，会产生较大的压力头梯度，这就需要较高的空间分辨率来准确描述湿润前沿的演变。为降低计算成本，本文采用了四叉树结构表示的自适应网格。通过评估单元格自身及其直接和间接相邻单元格内的含水率梯度来确定“问题单元格”指标。此外，对于远离湿润前沿的单元格则采用网格粗化处理。数值实验表明，将h自适应LDG方法与对角隐式Runge-Kutta时间积分方法结合使用，可在实际模拟中将计算时间缩短约三分之二。

引言：水在非饱和和饱和土层中的渗透是地下水水文学、环境工程以及页岩油开采等多个领域中的重要现象。在变饱和多孔介质中描述水运动的常用模型为Richards方程[1]。该方程结合了质量守恒方程与广义达西定律，其表达式为?θ?t=?·(K?H)??K?z，其中t为时间，θ为体积含水率，H为压力头，z为垂直坐标。在饱和土壤中，水力传导系数K被视为常数，而在非饱和区域，它则显著依赖于含水率θ或压力头H。为构成完整的数学模型，还需指定含水率与水力传导系数之间的本构关系。例如，Brooks-Corey模型[2]的表达式为Θ=θ?θrθs?θr={(HHb)?λ,H
章节节选：
空间离散化：本节介绍了针对richards方程（1.1）的ldg离散化方法。
时间离散化：假设在时间t时，元素kij上θh(x, z, t)的多项式展开式的系数为θij(t)，我们将这些系数集合表示为u(t)?(θ11(t),θ12(t),…,θn1n2(t))。通过ldg离散化（2.7）-（2.8）可得一组常微分方程组sdu(t)dt+l(t,u(t))=0，初始条件为u(t0)=u0。由于richards方程属于二阶偏微分方程，为保证数值稳定性，其显式时间积分的时间步长需满足δt=o(h2)的要求。在存在较大含水率梯度的情况下……
网格自适应：本节介绍了在存在较大梯度时，用于求解richards方程的自适应矩形网格的实现方法。首先描述了自适应网格的数据结构，接着阐述了选择需要细化或粗化的候选单元格的标准，最后介绍了网格自适应的具体流程。
数值实验：本节通过数值模拟展示了自适应ldg-dirk方法在求解richards方程（2.5）时的准确性及有效性。自适应网格技术在实际模拟中被广泛用于处理大梯度问题，同时还能节省计算资源，相关示例可见我们之前的研究[19]、[20]。表5展示了不同情况下的计算成本对比，结果明显体现了自适应数值方法的效率优势。
结论性论述：本文采用了h自适应ldg方法来求解richards方程，给出了ldg离散化的稳定性分析及最优收敛速率。数值实验验证了理论分析的正确性，实际模拟表明，与在均匀网格上计算相比，自适应ldg方案可将计算时间缩短约三分之二。
资金支持：第一作者得到了中国国家重点研发计划（2023yfa1009003）、山东省自然科学基金（zr2024ma036）以及中央高校基本科研业务费（24cx03005a）的支持。第三作者获得了国家自然科学基金（12131014、12371416）以及中央高校基本科研业务费（22cx03025a）的资助。最后一位作者则得到了simons基金会961585项目的支持。
利益冲突声明：作者们声明存在以下可能被视为潜在利益冲突的财务利益/个人关系：杨阳指出其差旅费用由simons基金会提供；田露露则表示其得到了中国国家重点研发计划的财务支持、行政支持、文章发表费用、设备、药品或物资支持以及差旅费资助。田露露还提及有其他形式的财务支持、行政支持、文章发表费用、设备、药品或物资支持。
田露露|林璐|郭辉|程娟|杨阳
中国石油大学（华东）理学院，中国青岛，266580 genuchten模型[4]，其表达式为θ=θ?θrθs?θr={(11+(α|H|)n)m,H<0,1,H≥0,K=KsΘ1/2(1?(1?Θ1/m)m)2，其中m=1?1/n。参数α和n需通过实验室实验确定。在这两种模型中，Ks为饱和水力传导系数，θr和θs分别表示剩余含水率和饱和含水率。由于这些本构函数具有高度非线性，且边界条件较为复杂，因此只有诸如Brooks和Corey模型[3]、[5]这类特定情况才能得到解析解。目前，学术界已采用多种方法对Richards方程进行数值求解，包括有限差分法[6]、[7]、有限元法[8]、有限体积法[9]、[10]、[11]、[12]，以及时空离散Galerkin方法[13]、[14]。近期，Congreve[15]对高阶局部离散Galerkin方法进行了误差分析。关于这些方法的详细综述可参见[16]。尽管大多数模拟都是在均匀网格上进行的，但使用传统数值方法且网格细化程度不足可能会导致结果不准确[17]。由于只有在压力头梯度较大的区域才需要高分辨率，因此采用自适应网格对于大幅降低计算成本至关重要[18]、[19]、[20]、[21]。网格自适应过程首先需要建立判断标准，以决定计算网格中的每个元素是应细化还是粗化。由于目前缺乏针对Richards方程的局部离散Galerkin离散化的可靠事后误差估计方法，本文采用了基于含水率梯度的启发式指标。作为采用元素级离散多项式空间的有限元方法家族的一员，DG（LDG）方法特别适合用于自适应网格构建。该方法最初由Reed和Hill[22]为中子输运方程所提出，随后Cockburn和Shu在其理论与数值发展方面做出了重要贡献[23]、[24]。受[25]中方法的启发，他们进一步提出了LDG方法[26]，用于处理含有二阶导数项的非线性对流-扩散方程。DG法和LDG法都具有诸多优点，如高精度、局部守恒性以及固有的并行性，因此被广泛用于各类高阶非线性偏微分方程的求解中[27]、[28]、[29]。本文重点研究了一种用于非饱和土壤中Richards方程的h自适应LDG方法，该方法能够自然实现质量守恒。我们的工作主要有三方面的贡献：首先，我们开发了一种稳健的自适应网格策略，能够根据含水率前沿的演变动态地调整单元格的细化与粗化程度。在含水率前沿存在剧烈梯度的实际模拟中，该方法会在关键区域细化网格，而在其他区域则进行粗化，从而提升计算效率。为实现网格的动态调整，需要一种高效的数据结构，为此我们为每个单元格采用了四叉树结构，该结构能够便捷地支持元素的插入与删除操作。自适应网格细化的另一个关键要素是“问题单元格”指标，我们通过整合当前单元格及其直接和间接相邻单元格的梯度信息来定义这一指标，从而提升了自适应模拟的稳定性与可靠性。与[12]、[30]的方法不同，我们的方法更加强调对远离湿润前沿的单元格进行动态粗化处理。其次，我们为应用于Richards方程的LDG方案建立了严格的稳定性分析与误差估计机制，证明了选择合适的数值通量对于保持算法稳定性至关重要。为处理误差估计中的对流项，我们利用了数值解的梯度与其对应的辅助梯度变量之间的重要关系。最后，我们将h自适应LDG方法与对角隐式Runge-Kutta（DIRK）方案相结合，与在均匀网格上的模拟相比，计算时间可减少约67%。本文的其余部分结构如下：第2节介绍相关符号，并给出Richards方程的稳定LDG离散化形式及相应的误差估计；第3节阐述DIRK时间积分方法；第4节介绍网格自适应框架，包括底层数据结构以及细化与粗化的规则；第5节通过数值实验展示h自适应LDG-DIRK方法的性能；第6节为结论性论述。 章节节选： 空间离散化：本节介绍了针对richards方程（1.1）的ldg离散化方法。 时间离散化：假设在时间t时，元素kij上θh(x, z, t)的多项式展开式的系数为θij(t)，我们将这些系数集合表示为u(t)?(θ11(t),θ12(t),…,θn1n2(t))。通过ldg离散化（2.7）-（2.8）可得一组常微分方程组sdu(t)dt+l(t,u(t))=0，初始条件为U(t0)=U0。由于Richards方程属于二阶偏微分方程，为保证数值稳定性，其显式时间积分的时间步长需满足Δt=O(h2)的要求。在存在较大含水率梯度的情况下…… 网格自适应：本节介绍了在存在较大梯度时，用于求解richards方程的自适应矩形网格的实现方法。首先描述了自适应网格的数据结构，接着阐述了选择需要细化或粗化的候选单元格的标准，最后介绍了网格自适应的具体流程。 数值实验：本节通过数值模拟展示了自适应ldg-dirk方法在求解richards方程（2.5）时的准确性及有效性。自适应网格技术在实际模拟中被广泛用于处理大梯度问题，同时还能节省计算资源，相关示例可见我们之前的研究[19]、[20]。表5展示了不同情况下的计算成本对比，结果明显体现了自适应数值方法的效率优势。 结论性论述：本文采用了h自适应ldg方法来求解richards方程，给出了ldg离散化的稳定性分析及最优收敛速率。数值实验验证了理论分析的正确性，实际模拟表明，与在均匀网格上计算相比，自适应ldg方案可将计算时间缩短约三分之二。 资金支持：第一作者得到了中国国家重点研发计划（2023yfa1009003）、山东省自然科学基金（zr2024ma036）以及中央高校基本科研业务费（24cx03005a）的支持。第三作者获得了国家自然科学基金（12131014、12371416）以及中央高校基本科研业务费（22cx03025a）的资助。最后一位作者则得到了simons基金会961585项目的支持。 利益冲突声明：作者们声明存在以下可能被视为潜在利益冲突的财务利益 个人关系：杨阳指出其差旅费用由simons基金会提供；田露露则表示其得到了中国国家重点研发计划的财务支持、行政支持、文章发表费用、设备、药品或物资支持以及差旅费资助。田露露还提及有其他形式的财务支持、行政支持、文章发表费用、设备、药品或物资支持。 田露露|林璐|郭辉|程娟|杨阳>
章节节选：
空间离散化：本节介绍了针对richards方程（1.1）的ldg离散化方法。
时间离散化：假设在时间t时，元素kij上θh(x, z, t)的多项式展开式的系数为θij(t)，我们将这些系数集合表示为u(t)?(θ11(t),θ12(t),…,θn1n2(t))。通过ldg离散化（2.7）-（2.8）可得一组常微分方程组sdu(t)dt+l(t,u(t))=0，初始条件为u(t0)=u0。由于richards方程属于二阶偏微分方程，为保证数值稳定性，其显式时间积分的时间步长需满足δt=o(h2)的要求。在存在较大含水率梯度的情况下……
网格自适应：本节介绍了在存在较大梯度时，用于求解richards方程的自适应矩形网格的实现方法。首先描述了自适应网格的数据结构，接着阐述了选择需要细化或粗化的候选单元格的标准，最后介绍了网格自适应的具体流程。
数值实验：本节通过数值模拟展示了自适应ldg-dirk方法在求解richards方程（2.5）时的准确性及有效性。自适应网格技术在实际模拟中被广泛用于处理大梯度问题，同时还能节省计算资源，相关示例可见我们之前的研究[19]、[20]。表5展示了不同情况下的计算成本对比，结果明显体现了自适应数值方法的效率优势。
结论性论述：本文采用了h自适应ldg方法来求解richards方程，给出了ldg离散化的稳定性分析及最优收敛速率。数值实验验证了理论分析的正确性，实际模拟表明，与在均匀网格上计算相比，自适应ldg方案可将计算时间缩短约三分之二。
资金支持：第一作者得到了中国国家重点研发计划（2023yfa1009003）、山东省自然科学基金（zr2024ma036）以及中央高校基本科研业务费（24cx03005a）的支持。第三作者获得了国家自然科学基金（12131014、12371416）以及中央高校基本科研业务费（22cx03025a）的资助。最后一位作者则得到了simons基金会961585项目的支持。
利益冲突声明：作者们声明存在以下可能被视为潜在利益冲突的财务利益/个人关系：杨阳指出其差旅费用由simons基金会提供；田露露则表示其得到了中国国家重点研发计划的财务支持、行政支持、文章发表费用、设备、药品或物资支持以及差旅费资助。田露露还提及有其他形式的财务支持、行政支持、文章发表费用、设备、药品或物资支持。
田露露|林璐|郭辉|程娟|杨阳
中国石油大学（华东）理学院，中国青岛，266580>