显式不可压缩SPH方法的理论与数值评估

《Computers & Mathematics with Applications》：Theoretical and numerical assessment of explicit incompressible SPH method

　　

摘要

显式不可压缩SPH（EISPH）方法是在基于投影的不可压缩SPH（ISPH）框架基础上发展而来的一种完全显式方法，其设计初衷在于同时实现基于投影的ISPH方法所固有的数值稳定性和完全显式算法所特有的计算效率。尽管该方法的应用日益广泛，但其核心假设——将压力泊松方程（PPE）的迭代求解截断为单步——却从未得到严格的验证。本研究填补了这一空白：首先建立了EISPH与弱可压缩δ-SPH方法之间的理论类比关系。该分析表明，传统EISPH等价于一种具有非物理高可压缩性和过量数值扩散的δ-SPH方法，从而损害了其精度和稳定性。基于这些发现，我们针对CFL条件提出了适当的应用建议，并提出了一种稳定化方案以缓解其缺陷。随后，我们进行了基本的基准测试，以系统地评估并数值验证这些固有局限性以及所提出参数设置和方案的有效性。

引言

光滑粒子流体动力学（SPH）方法[1]是一种无网格拉格朗日方法，近年来已被应用于多种工程问题。其基于粒子的特性使其特别适合模拟复杂的自由表面流动、流固耦合以及涉及大变形等现象，而这些现象往往使传统的基于网格的方法面临巨大挑战（例如[2]、[3]、[4]、[5]、[6]）。

SPH方法在流体动力学领域的发展主要沿两条不同的路径展开。最初的框架称为弱可压缩SPH（WCSPH），由Monaghan[7]提出。该方法采用显式时间积分格式，并依靠状态方程将压力与密度关联起来，从而为方程组提供闭合条件。然而，初始WCSPH公式的一个显著缺点是其容易受到数值不稳定的影响，通常需要添加人工稳定化项，如人工粘性[8]，以维持模拟的稳定性。近年来，δ-SPH方法[9]在连续方程中引入数值密度扩散项，已成为显式SPH框架中的主流方法。相比之下，由Shao和Lo[10]发展的不可压缩SPH（ISPH）方法采用投影法[11]来严格执行不可压缩条件。ISPH方法本身具有稳定性，无需人工处理，并且由于其隐式算法的特性可以采用较大的时间步长。然而，ISPH的主要缺点在于求解压力泊松方程（PPE）隐式离散化产生的大型稀疏矩阵方程所带来的计算成本和并行效率降低。

第三种方法是显式不可压缩SPH（EISPH）方法，由Hosseini等人[12]提出。该方法旨在结合WCSPH和ISPH两者的优势。其公式基于ISPH方法，但通过引入压力场在单个时间步内不会显著变化的假设，使计算方案变为完全显式。这使得压力泊松方程矩阵的Jacobi迭代求解可以在仅一次迭代后截断。因此，EISPH据称既保留了隐式算法（ISPH）的鲁棒性和大步长能力，又受益于显式方案（δ-SPH）的实现简单性和高并行效率。表1总结了EISPH文献，表明该方法近年来获得了越来越多的关注和应用。

然而，尽管其使用日益增多，支撑EISPH核心假设的根本有效性在以往研究中从未得到严格检验。具体而言，在单步非迭代步骤后终止矩阵求解的不利影响在很大程度上被忽视了。此外，该方法基本数值特性的全面验证也尚未得到系统性考察。现有的验证研究数量有限，且通常局限于短周期基准问题，如溃坝问题，关于该方法长期精度和稳定性的研究尚未见报道。这种情况揭示了文献中的一个关键空白：该方法往往在缺乏充分数学和数值验证的情况下被应用于实际工程问题。

本研究通过从根本上调查EISPH方法的理论和数值特性来填补这一空白。首先，通过数学推导展示了EISPH与δ-SPH公式之间的理论相似性，从而揭示了EISPH方法固有的潜在问题（第3节）。其次，基于理论分析中发现的缺陷，我们为CFL条件及其系数值提出了适当的计算设置以保证稳定性；并开发了一种稳定化方案——动态压力拉普拉斯（DPL）模型——以改善自由表面稳定性（第4节和第5.2节）。第三，进行了一系列基本基准测试（第6节）。这些测试旨在系统评估方法的基本特性（包括精度、稳定性和守恒性），数值验证其理论局限性，并展示所提出计算设置和方案的有效性。

控制方程

EISPH和δ-SPH中使用的控制方程如下：

Dρ/Dt = -ρ?·u
Du/Dt = -1/ρ ?p + g + Π

其中ρ是密度，t是时间，u是速度向量，p是压力，g是重力加速度向量，Π是粘性耗散项。

SPH插值

光滑粒子流体动力学（SPH）方法建立在积分插值原理的基础上。任意场变量φ在位置r处的值通过对邻近位置ξ处的值在紧致支撑域上进行积分来近似。

压力更新公式

为阐明EISPH与δ-SPH之间的关系，对它们各自的压力更新方案进行比较。各方法中压力的时间演化可表示如下：

对于EISPH，由式(17)：

Δp?/Δt_EI = -h2/(A??·ρ·Δt_EI2)·Σ? u??·?w??·V? + 1/Δt_EI·(-p?? + A??/A??)

对于δ-SPH，由式(28)：

Δp?/Δt_WC = -c?2/ρ?·Σ? u??·?w??·V? + δc?A??/h·(-p?? + A??/A??)

其中时间步长Δt_EI和Δt_WC分别由其对应的CFL条件（式(18)和(29)）确定：

Δt_EI = C_EI·h/u_max
Δt_WC = C_WC·h/c?

需要注意的是

EISPH中适当时间步限制的调查

如第3节所示，传统EISPH方法等价于δ-SPH模型，允许约5%或更高的可压缩性。因此，本节探讨了改善方法稳定性的方法并确定EISPH的最优时间步限制。

EISPH的稳定化方案

在前一节中，通过数学分析与δ-SPH模型的比较，揭示了EISPH方法的关键和理论问题。另一方面，如表1所示，EISPH方法已在众多研究中被应用，这些工作报告的结果具有表面可接受的稳定性和精度。在本节中，我们回顾了那些已产生表面稳定结果的改进方案的文献。此外，我们提出了两种新方法来……

EISPH的数值验证

如第1节所述，EISPH方法的验证研究主要集中在溃坝等短期场景上。关键特性，包括长期稳定性和守恒特征，尽管对工程应用至关重要，但尚未得到充分调查。为填补这一空白，本节使用三个基本测试用例来验证数值方案：无粘静水柱、无粘驻波和振荡液滴。

结论

本研究对显式不可压缩SPH（EISPH）方法进行了严格的理论和数值重新评估，该方法旨在结合基于投影的不可压缩方法的稳定性与显式算法的效率。通过详细的数学推导，本工作确立了该方法的基本数值特性，并考察了其与弱可压缩δ-SPH公式之间的关系。此外，本文填补了一个关键空白……

生成式AI及写作过程中AI辅助技术的使用声明

在本文准备过程中，第一作者使用了Gemini进行写作校对。在使用该工具/服务后，作者对内容进行了审阅和编辑，并对出版物的内容承担全部责任。

CRediT作者贡献声明

Yuma Shimizu：概念构思、方法论、软件、验证、形式分析、调查研究、可视化、撰写原稿、撰写——审阅与编辑。

Dongfang Liang：方法论、监督、撰写——审阅与编辑。

Hitoshi Gotoh：方法论、监督、撰写——审阅与编辑。

利益冲突声明

作者声明不存在利益冲突。

Yuma Shimizu | Dongfang Liang | Hitoshi Gotoh

大阪大学苑大学（Momoyama Gakuin University）工学部，大阪府和泉市学之井1-1，邮编594-1198，日本