一种用于三维稳态无限制渗流问题的自适应多面片等几何分析方法

《Engineering Analysis with Boundary Elements》：An adaptive multi-patch isogeometric analysis for three-dimensional steady-state unconfined seepage problems

　　习继华|吴晨|王琳|桑达拉拉詹·纳塔拉詹|方伟华|于天堂摘要传统用于自由面渗流分析的移动网格方法往往存在网格变形和数值不稳定的问题。而固定网格方法则在精度或计算效率方面存在妥协，尤其是在三维问题中。为克服这些缺陷，本文采用自适应等几何分析框架来模拟三维稳态无限制渗流问题，旨在获得高效且稳定的数值结果。文中引入了截断分层NURBS（TH-NURBS）用于局部细化，结合渗透率调整机制，该框架能够在无需持续重网格的情况下实现精确的自由面定位。通过多块贴片结合尼奇方法，可以对复杂的三维几何体进行建模与仿真。所提出的自适应策略依赖于事后误差估计器，在迭代过程中识别并细化关键区域，从而在无需牺牲数值精度的同时显著提升计算效率。该框架在多种渗流条件下的稳健性已得到验证，数值结果与参考解高度一致。研究表明，这一框架为未来多物理场耦合及渗流分析的实际工程应用奠定了基础。引言在实际的水利和岩土工程问题中，如土石坝、边坡以及地下岩土工程[1]、[2]、[3]，无限制渗流现象十分普遍。自由面的位置变化、渗流场的分布以及渗流面边界的确定，直接决定了结构的渗流量、渗透稳定性以及工程安全评估。因此，开展准确高效的数值模拟仍是无限制渗流分析领域的重要研究方向。有限元法（FEM）因其能够灵活处理复杂几何形状和边界条件，被广泛用于无限制渗流建模。然而，传统的移动网格方法[4]、[5]、[6]容易产生网格变形，且在描述水平分层和流固相互作用时存在困难。为解决这些问题，固定网格公式[7]逐渐被更多采用，即在求解渗流场时保持计算网格不变。现有的固定网格策略大致可分为两类：（i）启发式调整方法，包括单元渗透率矩阵调整法[8]以及流动调整法[9]、[10]；（ii）变分不等式方法[11]。尽管这些方法效果良好，但存在一些局限性，如计算成本相对较高、收敛性对初始猜测值敏感，以及在强非线性情况下的稳健性较差。此外，当地下水位与过渡单元相交时，渗透率的赋值可能产生歧义，从而导致迭代振荡和解决方案精度下降。近年来，等几何分析（IGA）[12]被引入到无限制渗流问题中[13]、[14]，它具有几何精确性、高阶连续性，并且无需手动生成网格。Nguyen等人[15]提出了一种基于IGA的饱和流问题数值方法。非均匀有理B样条（NURBS）基函数的应用使得场变量的近似精度更高，从而能更准确地捕捉湿润前沿。Shahrokhabadi和Vahedifard[16]通过结合IGA与随机场，提出了一种概率框架，用于数值模拟非饱和土壤中的渗流现象。Kazemzadeh-Parsi[17]则利用IGA解决了具有未知自由面、不规则域以及空间变化渗透率的稳态饱和无限制渗流问题。上述研究中均广泛使用了NURBS，它虽然具备全局高阶连续性，但缺乏局部细化的能力。这一局限促使人们开始使用可局部细化的样条，如T样条[18]、[19]、LR B样条[20]以及PHT样条[21]。尽管这些样条能够实现有效的局部网格细化，但其基函数在局部细化过程中可能存在线性无关性不足[22]、[23]的问题，或者出现性能退化现象[24]。层次化B样条（HB-spline）[25]、[26]以及截断层次化B样条（THB-spline）[27]、[28]提供了一种通过嵌套样条空间来解决上述问题的替代方法。此外，THB-spline还通过截断机制恢复了样条基函数所具有的单位分解特性。本文采用了其扩展形式TH-NURBS，它在保留所有前述优势的同时，还能实现对圆锥几何体的精确描述[29]。与此同时，IGA也逐渐与边界元法（BEM）相结合，形成了等几何边界元法（IGABEM），该方法已成为解决渗流和力学问题的重要手段[30]。由于IGABEM在局部网格处理方面的能力有限，自适应技术也成为其发展的关键方向。Neto和Leonel[31]提出了一种基于超奇异残差的自适应IGABEM策略，用于三维问题，该方法利用应变误差估计器驱动网格细化，对于均匀和各向异性弹性介质都具有良好的收敛性。Gantner和Praetorius[32]提出了一种基于层次化B样条的自适应IGABEM算法，分析了加权残差事后误差估计器的可靠性和最优收敛性，并证明了其在三维线性弹性问题中的有效性。Feischl等人[33]则为二维IGABEM提出了加权残差事后误差估计器以及一种可实现最优收敛速率的自适应算法。这些研究充分展示了IGA与自适应技术结合在各类工程领域的广泛应用，为本论文的研究内容提供了有力补充。尽管样条技术日益丰富，但复杂模型的参数化仍然是一个挑战。由于NURBS具有张量积性质，单个贴片无法表示任意复杂的几何形状，因此多贴片建模成为必然。关于非重叠贴片耦合的研究主要聚焦于四种主流方法：惩罚法[34]，虽然简单但易出现病态问题；拉格朗日乘数法[35]，需要满足LBB条件[36]，并且存在鞍点问题，不过其迭代版本（如IGATI[37]）可缓解这一问题；增广拉格朗日法[38]结合了两者的优点，但仍受LBB约束；以及尼奇方法[29]，通过在界面处施加弱积分形式来实现耦合。最后一种方法在计算精度和收敛性方面表现更优，且在耦合非协调网格时更具灵活性，因此本文采用了该方法。为在不牺牲精度的前提下减少计算时间，需要一种局部自适应网格。从计算域的初始粗略离散化开始，自适应方案通过事后误差估计逐步进行局部网格细化，直至达到预定精度。误差估计通常采用两种成熟的方法：基于恢复的误差估计器[39]和基于残差的误差估计器[40]。前者因其简单性、高计算效率以及易于实现而备受青睐。在这一框架下，Bekele等人[41]通过采用基于LR B样条的自适应IGA方法，解决了稳态地下水流动中的若干数值难题。最近，吴晨等人[42]提出了一种基于TH-NURBS的自适应IGA方法，用于二维稳态无限制渗流问题，该方法无需重构网格即可精确追踪自由面，同时动态细化其周边关键区域。本研究在参考文献[42]关于无限制渗流问题的基础上进行了拓展，将问题扩展到三维复杂域中。其主要贡献可归纳为以下几点：?将原本为二维稳态无限制渗流问题设计的基于TH-NURBS的自适应多贴片IGA方法扩展到了三维场景。?通过结合尼奇方法，多个贴片相互耦合，从而为复杂的三维几何体保持场变量的连续性。?通过三维基于恢复的误差估计来评估水力梯度场的误差，以此驱动自适应网格细化。?将固定网格方法与自适应细化相结合，从而在保证计算精度的同时降低计算成本。本文的其余结构如下：第2节详细阐述三维无限制渗流的基本理论；第3节介绍无限制渗流问题的数值实现方法，包括基于TH-NURBS的IGA离散化、自适应策略以及数值算法；第4节通过多个数值算例验证该框架的适用性与准确性；第5节对全文进行总结，并提出未来研究方向。问题描述图1示意性地展示了一个三维稳态无限制渗流问题，计算域被划分为两个相邻的子域Ω1和Ω2。这两个子域通过一个公共界面Γ?相连，界面上的单位外法向量分别记为n1和n2。假设水流仅发生在饱和带内，每个子域中的水流都会在干区与湿区之间形成自由面Γf。区域概述TH-NURBS是THB样条的合理延伸。它不仅继承了NURBS在描述圆锥曲线时的几何精确性，还具备THB样条的局部细化能力。给定三个参数方向上的初始结点向量Ξ0、H0和Z0，通过逐次插入结点，可以生成相应的嵌套结点向量序列{Ξ}a=0N、{H}a=0N、{Z}a=0N，其中N表示最大细化级别。相应地，也会形成嵌套的NURBS空间结构。数值算例为了展示所提方法的数值性能，给出了四个数值算例。第一个算例通过解析解进行验证，接下来的两个算例则在复杂条件下与已发表的二维数据对比，以验证方法的正确性，最后一个算例则通过与COMSOL的结果对比，来检验其在真实三维问题中的实用性。基于先前的研究[48]，所有数值算例均采用二阶样条基函数，以便在保持较低计算成本的同时获得较高的计算精度。结论本文在我们之前的工作[42]基础上，扩展了IGA框架，用于解决三维稳态无限制渗流问题。为对三维复杂结构进行建模，首先对域进行分解，然后采用尼奇方法来确保子域界面之间的耦合。为在保持解决方案精度的同时大幅降低三维问题的计算成本，本文引入了基于恢复的误差估计器来实现自适应网格细化。CRediT作者贡献说明习继华：撰写——原始草稿、可视化、软件、方法论、研究、正式分析、数据整理。吴晨：撰写——审阅与编辑、验证、软件、方法论、研究。王琳：撰写——审阅与编辑、验证、软件、方法论、研究。桑达拉拉詹·纳塔拉詹：撰写——审阅与编辑、验证、方法论、研究。方伟华：撰写——审阅与编辑、验证、研究。于天堂：撰写——审阅与编辑。人工智能使用声明在准备本文档的过程中，作者使用了DeepSeek来优化句子结构。使用DeepSeek之后，作者对内容进行了必要的审阅和修改，并对最终发表的内容负全责。利益冲突声明作者声明不存在任何可能影响本文研究的已知财务利益或个人关系。