基于背景网格增强的不可压SPH方法用于自由表面流动问题建模

《Engineering Analysis with Boundary Elements》：A background grid enriched incompressible SPH method for free surface flow problems modeling

　　

摘要我们提出了一种背景网格增强的不可压缩光滑粒子流体动力学方法（BG-ISPH），用于模拟自由表面流动问题。该方法通过背景网格方法计算压力泊松方程（PPE），同时采用SPH方法求解动量方程。所提出的BG-ISPH方法兼具两种数值方法的优势：网格方法能够实现压力和的高效精确计算；同时保留了SPH方法在处理大变形自由表面问题方面的优势。因此，与传统ISPH方法相比，该方法提高了稳定性和数值精度。我们将该模型应用于一系列基准问题的模拟，涉及不同类型的自由表面，并将数值结果与实验观测值及文献中的其他模拟结果进行了比较。此外，我们研究了背景网格尺寸对模型精度的影响，并提出了最优的网格与粒子尺寸比范围为1.0–1.75。结果表明，BG-ISPH模型在模拟自由表面流动问题时具有稳定性和精度。引言不可压缩自由表面流动建模是一个关键而具有挑战性的研究课题。一个典型的例子是溃坝问题，它涉及高度非线性的现象，如大界面变形、流体破碎、飞溅以及强流体-结构相互作用。传统的网格型计算流体力学（CFD）方法——如有限元方法（FEM）[1]、[2]、有限差分法（FDM）[3]、[4]和有限体积法（FVM）[5]、[6]——已被广泛应用于这些模拟中。这些研究显著推动了CFD的发展。然而，网格型数值方法在尝试准确高效地模拟具有大自由表面变形的流动时可能会遇到若干挑战。主要问题是欧拉框架中的网格畸变，这可能会严重损害计算的精度和稳定性，在某些情况下甚至导致模拟失败[7]。此外，为了准确追踪自由表面，还需要专门的界面捕捉技术。例如，Hirt [8]和Osher [9]分别提出了流体体积法（VOF）和水平集方法。这些开创性方法展现了良好的性能，但也引入了相对复杂的计算程序。为了解决这些挑战，拉格朗日型无网格方法作为鲁棒的替代方案应运而生。其固有的避免网格畸变的能力使其在模拟具有移动自由表面的不可压缩流动问题时具有优势[10]、[11]。其中，光滑粒子流体动力学（SPH）方法[12]、[13]已成为应用最广泛的无网格方法之一。作为一种基于粒子的数值技术，SPH方法使用携带物理属性（如质量、密度、压力、速度）的粒子来离散化计算域，并采用基于核函数的插值来逼近场变量及其导数。该框架自然地规避了网格方法中的界面追踪难题，能够稳健地模拟大变形和复杂的界面动力学。因此，SPH已成功应用于多种水动力学问题的建模，包括摇摆流动[14]、[15]、多相流动[16]、[17]、[18]、波浪冲击[19]、[20]、[21]等。这些应用凸显了SPH在模拟高度非线性、剧烈自由表面流动方面优于传统网格方法的能力。对于自由表面流动模拟，SPH提供了两种主要途径。一种常用方案是通过采用弱可压缩假设发展而来，因此被称为弱可压缩SPH（WCSPH）方法[22]。在WCSPH中，粒子压力通过状态方程（EoS）计算，该方程将压力描述为密度及其他状态变量（如声速和温度）的函数，从而实现显式时间积分。尽管其简单且具有计算效率，但WCSPH存在明显的缺点，例如无法严格满足不可压缩条件，导致不可压缩问题的体积守恒误差。此外，声速在CFL条件下限制了时间步长[23]，这限制了计算效率。WCSPH方法的另一个令人不快之处在于由高频压力波引起的非物理压力振荡[24]，特别是在模拟剧烈流动过程时。也有研究工作致力于缓解压力振荡，例如Marrone等人[25]发展的δ-SPH方法，该方法在动量方程中引入人工粘性，以及其各种改进版本[26]、[27]、[28]。这些方法取得了良好的效果，并在研究中被广泛应用。然而，正如Ferrari等人[29]所指出的，这些方法也会引入人为干扰和过度耗散。另一种方法，即由Cummins和Rudman[30]开创的不可压缩SPH（ISPH），在粒子离散化的域上求解压力泊松方程（PPE）以强制执行不可压缩条件。与WCSPH不同，ISPH保持恒定的粒子密度，严格保证体积守恒。由于压力是隐式获得的，与WCSPH相比可以使用更大的时间步长。此外，使用ISPH方案可以自然地产生平滑稳定的压力场[31]、[32]。然而，在每一步求解隐式PPE时实现高精度和稳定性仍然是一个重大挑战，尤其是当计算域被不规则分布的粒子离散化时。为解决这一问题，人们做出了许多尝试来提高ISPH方法的数值效率和精度。例如，开发了PPE源项和拉普拉斯算子的更高阶离散化方案[33]。然而，即使具有精确的差分算子模型，也很难保证精确的无散度速度场。Gotoh等人[34]提出了误差补偿源（ECS）方案，在PPE源项中引入两个误差缓解项来纠正瞬态和累积偏差，从而提高精度。其他尝试包括梯度校正（GC）[35]、二次半解析有限差分插值（QSFDI）[36]和优化泊松求解器[37]。这些研究在一定程度上推动了ISPH的发展；然而，所有方法都依赖于基于SPH离散化的PPE求解，这存在固有的局限性，包括粒子无序、核近似误差和可能的矩阵非对称性，这些都会降低压力求解的精度和稳定性。为解决上述问题，我们采用了一种先进的混合方案，在背景网格上求解PPE，这增强了传统ISPH方法的计算稳定性，同时保持压力计算的高精度。该方案的主要思想是通过执行插值将计算结果从一组离散化单元（如节点、粒子、体积）传输到另一组。实际上，将拉格朗日粒子与欧拉网格耦合的概念有着悠久的历史和坚实的基础。粒子-in-Cell（PIC）方法最初由Evans等人[38]提出，是最早将粒子属性映射到网格以进行高效场计算的技术之一。在此基础上，物质点方法（MPM）[39]将这种混合方法扩展到固体力学领域，其中背景网格用于求解运动方程，而粒子携带变形历史。这些成熟方法展示了结合网格效率和粒子灵活性的有效性。近年来，类似的想法获得了关注；例如，Bui等人[40]和Feng等人[41]分别开发了流体和固体粒子之间的两相显式ISPH和WCSPH插值模型，以模拟多孔介质流体-结构相互作用问题。后续工作将其扩展到网格-粒子映射，例如Wang等人[42]将运动粒子半隐式（MPS）方法[43]与背景网格集成，以增强基于投影的压力计算，在MPS框架中求解压力之前提高PPE源项的连续性和平滑性。类似地，You等人[44]将背景网格与SPH结合以细化基于EoS的源项，从而产生更精确的压力场。此外，Xu等人[45]采用ISPH-FVM耦合方法来提高液体中气泡流动模拟的精度和效率。所有这些背景网格辅助方案都显著改善了模拟性能，展示了将背景网格与粒子方法结合的巨大优势。上述方法主要集中在利用背景网格来平滑源项和辅助局部插值。此外，Hwang [46]开发了运动粒子压力网格（MPPM），该方法完全在背景网格上求解泊松方程。在MPS框架内，该方法采用固定的欧拉网格来求解PPE，显著改善了传统MPS方法的稳定性、收敛性和精度。在MPPM的基础上，提出了一些扩展方法，如运动粒子水平集（MPLS）[47]和非结构运动粒子压力网格（UMPPM）[48]，以进一步提高数值精度并实现对复杂几何边界的处理。在本文中，我们通过在ISPH方法中引入类似于MPPM方法的背景网格泊松方程求解器，提出了一种改进的ISPH方法，称为背景网格增强的不可压缩光滑粒子流体动力学（BG-ISPH）方法。通过采用这种背景网格求解器，与常规ISPH方法相比，压力和的计算效率都可以得到改善。与MPPM方法不同，动量方程完全通过采用SPH离散化来求解，并且所有边界条件都在粒子层面施加，以便处理移动边界，如造波器（见第3.6节）。在后续章节中，我们通过模拟各种基准自由表面问题来验证BG-ISPH模型。数值结果与实验数据和现有文献进行了比较，证明了该模型在精度和计算稳定性方面的稳健性能。我们提出的BG-ISPH模型取得了良好的效果，为复杂流体动力学挑战提供了一个有前景的解决方案。为解决上述问题，我们开发了一种背景网格增强的不可压缩光滑粒子流体动力学方法（BG-ISPH），其中PPE在固定的结构化背景网格上求解，而动量方程仍在SPH粒子框架中进行离散化和推进。核心动机是将压力计算与限制传统ISPH的不规则粒子拓扑结构解耦，同时保留SPH在追踪强烈变形的自由表面和移动界面方面的拉格朗日优势。这样，所提出的方法旨在将基于网格的压力求解的数值鲁棒性和效率与基于粒子的流动表示的几何灵活性结合起来。将拉格朗日粒子与欧拉网格耦合的想法并非新鲜事物。早期的混合方法，如粒子-in-Cell（PIC）[38]方法和物质点方法（MPM）[39]，证明了在粒子和网格之间传递信息可以有效利用两种描述方法的互补优势。遵循这一总体理念，一些研究在流体模拟的粒子方法中引入了背景网格。例如，一些背景网格辅助的SPH和MPS方法[43]被开发出来以平滑PPE源项、改善局部插值质量或增强压力评估。Wang等人[42]将背景网格纳入基于投影的MPS中，以在求解压力之前改善PPE源项的连续性。You等人[44]将背景网格与SPH耦合以细化基于状态方程的压力计算。Xu等人[45]进一步探索了ISPH-FVM耦合以用于气泡流动模拟。这些研究证实，引入背景网格可以改善数值性能；然而，在这些方法中，背景网格主要起辅助作用，例如在源项平滑、插值支持或局部校正方面，而不是作为主要的压力求解框架。一个更密切相关的研究方向由运动粒子压力网格（MPPM）[46]及其扩展方法[47]、[48]代表，其中泊松方程在MPS方法背景下直接在欧拉网格上求解。这些方法证明了将PPE求解从不规则的粒子转移到规则网格上可以显著改善稳定性和收敛性。然而，本工作与MPPM类方法在方法基础和目标应用方面都不同。首先，所提出的方法建立在SPH框架而非MPS框架内，因此动量方程完全基于粒子并与SPH离散化一致。其次，在本BG-ISPH方法中，自由表面和固体壁边界的处理完全在粒子层面进行，而不是在背景网格上重建和施加。这种设计避免了在欧拉节点上反复进行几何重建，因此对于涉及强烈变形自由表面、移动边界和复杂几何形状的问题特别方便。从这个角度来看，本工作的方法论创新不仅仅在于引入背景网格，而在于在网格和粒子表示之间建立了一个明确划分的数值责任分工。

与传统的ISPH方法相比，本文提出的方法避免了在不规则的粒子分布上直接求解压力泊松方程（PPE），从而缓解了与粒子无序、核函数不一致以及非对称线性方程组相关的误差。与仅将背景网格用于平滑或插值支持的背景网格辅助SPH或MPS方法相比，本方法完全在背景网格上求解PPE。与MPPM类方法相比，本文提出的BG-ISPH方法保留了基于SPH的动量求解和粒子级别的边界执行，这既保留了SPH在自由表面流动模拟中的灵活性，又受益于稳健的基于网格的压力计算。因此，本研究被定位为针对不可压缩自由表面流的混合投影框架，而非对现有SPH求解器的简单辅助网格修正。除了提出该框架外，我们还考察了网格-粒子耦合中固有的一个重要实际问题，即网格与粒子尺寸比对求解精度的影响。这一方面在数值算例中进行了系统研究，并为该方法的实际实施提供了有用的指导。

在接下来的章节中，我们通过模拟各种基准自由表面问题来验证BG-ISPH模型。数值结果与实验数据和现有文献进行了比较，展示了该模型在精度和计算稳定性方面的稳健性能。我们提出的BG-ISPH模型已取得了良好的效果，为复杂流体动力学挑战提供了一种有前景的解决方案。

**SPH对控制方程的离散化**

用于模拟自由表面流的控制方程的连续和动量方程可以写成拉格朗日形式如下：

??u=0

dudt=?1ρ?p+g+ν?(??u)

其中u表示速度；t表示时间；ρ表示密度；p表示压力；g表示重力加速度；ν表示运动粘度。

光滑粒子流体动力学（SPH）的基本概念是通过连续插值或积分来逼近场函数。

**数值算例**

本节中，为了验证所提出的BG-ISPH耦合模型的稳健性、精度和通用性，我们进行了几个典型自由表面流动案例的数值模拟，包括静水试验、动水试验、流体-结构相互作用试验和移动边界试验。模拟结果与其他数值结果或实验数据进行了比较和分析。此外，我们还研究了网格尺寸对模拟结果的影响，并进行了分析。

**总结与结论**

本研究采用一种改进的网格-粒子混合方案，将背景网格上的压力场求解纳入ISPH框架中，以模拟不可压缩自由表面流，本文将其称为BG-ISPH方法。BG-ISPH充分利用了有限差分法（FDM）和SPH的互补优势，在每个时间积分步骤中计算流场。它保留了基于网格的方法在求解偏微分方程中的效率，同时利用了拉格朗日方法的优势。

**CRediT作者贡献声明**

傅一涵：撰写—审稿与编辑，撰写—初稿，可视化，调查研究，形式分析，数据整理。孙征：监督指导，资金获取。冯殿雷：撰写—审稿与编辑，监督指导，项目管理，方法论，资金获取。

**利益冲突声明**

作者声明，他们没有已知的可能影响本文所述工作的竞争性经济利益或个人关系。

**致谢**

冯殿雷感谢中国中央高校基本科研业务费项目提供的资金支持。冯殿雷和孙征感谢江西源头地区水生态保护重点实验室（编号：WEC-2024-01）提供的资金支持。特别感谢张雄教授进行的有益讨论。

傅一涵 | 孙征 | 冯殿雷