一种改进的二阶离散化方案，用于MPS方法中的自由表面及其附近粒子，该方案借助了表面网格来辅助离散化过程

提出基于表面网格的改进第二阶离散化方法，通过动态构建自由表面附近的局部网格补偿粒子分布不均问题，解决高阶离散化中粒子重叠导致的数值发散，验证表明该方法在压力计算精度、表面粒子分布均匀性和边界一致性方面显著优于传统虚拟粒子方法。

　　

李根|廖云龙|姚培涛

华南理工大学电力工程学院，中国广州510641

摘要

移动粒子半隐式（MPS）方法是一种用于处理涉及复杂变形、破碎和各种流体界面聚合的自由表面流动的新兴数值方法。然而，MPS方法的高阶离散化方案仍然不够完善。自由表面附近的不均匀粒子分布容易导致数值发散。传统的基于虚拟粒子的二阶离散化方案会将自由表面及其附近粒子的离散化降级为低阶格式，从而导致误差的累积和传播。为了解决这些问题，借助表面网格开发了一种改进的二阶离散化方案。在自由表面构建的表面网格为表面粒子的二阶离散化提供了所需的位置信息，并弥补了粒子数密度的不足。通过对表面网格分辨率进行了敏感性分析，并确定了粒子与网格尺寸的比例，以平衡计算效率和准确性。与之前的基于虚拟粒子的二阶方法相比，所提出的方法能够准确离散化自由表面粒子，防止了由不均匀粒子分布引起的误差累积。通过模拟四个基准案例（静止水池压力、溃坝流、椭圆滴体演化和方形滴体旋转）进行了验证。结果表明，所提出的基于表面网格的方法在压力计算准确性、自由表面粒子分布均匀性和表面一致性方面表现出更优的性能。

引言

自由表面流动在自然界和工程领域中非常普遍，其中的关键难点是准确解析各种流体界面的复杂变形、破碎和聚合。在解决这类问题时，传统的基于网格的方法需要额外的界面捕获和重构技术。然而，在欧拉框架中，对流项的数值耗散会阻碍基于网格的方法准确跟踪界面运动[1]。相比之下，移动粒子半隐式（MPS）方法[2]是一种新兴的数值方法。在这种方法中，不同的流体由不同类型的粒子表示，通过跟踪粒子相互作用来直接捕获相界面[[3]，[4]，[5]，[6]，[7]，[8]]。

MPS方法中的粒子相互作用模型（例如梯度、散度和拉普拉斯算子）是在假设粒子分布均匀的情况下推导出来的。然而，数值模拟中的粒子分布通常是不均匀的，这会导致这些相互作用模型中的误差累积。为了提高MPS方法的准确性，学者们进行了广泛的研究。值得注意的是，在守恒和稳定性方面，Khayyer和Gotoh[9]开发了修正MPS（CMPS）方法，以确保严格满足线性和角动量守恒。为了抑制非物理压力振荡，他们随后提出了MPS-HS（高阶源）[10]和MPS-HL（高阶拉普拉斯）[11]方法，分别显著提高了源项和拉普拉斯模型的准确性。此外，为了提高空间稳定性，Khayyer和Gotoh[12]开发了优化MPS方法来防止粒子聚集，而Khayyer等人[13]对粒子正则化方案的准确性和守恒特性进行了比较研究，并提出了一种优化的粒子移动方案。在处理非物理空洞和自由表面奇异性方面，Tsuruta等人[14]提出了空间势粒子（SPP）方案，该方案在空洞区域加入虚拟粒子以强制体积守恒的一致性。

与此同时，研究人员也致力于减少基本的离散化误差。Randles等人[15]在平滑粒子流体动力学（SPH）框架内引入了一个校正矩阵来减少离散化误差。这种方法随后被Suzuki等人[16]和Khayyer等人[12]应用于MPS方法。基于类似的数学基础，最小二乘MPS（LSMPS）方法[[17]，[18]，[19]]被开发出来，通过最小化泰勒级数近似的加权残差来实现高精度模拟。虽然LSMPS在理论上是严格的并且可以达到任意的一致性阶数，但它通常需要为每个时间步的每个粒子解决一个局部线性方程组来确定微分算子。这个过程会导致较高的计算成本和实现复杂性。相比之下，校正矩阵通过一个单一的、显式的修改向量来考虑不均匀的粒子分布，从而提供了更大的简化。利用这一优势，Duan等人[20]基于一阶校正矩阵将散度和拉普拉斯模型的准确性分别提高到了第一阶和零阶收敛。最近，Duan等人[21]基于二阶校正矩阵进一步提高了模型准确性，其中梯度模型和散度模型达到了二阶收敛，而拉普拉斯模型达到了第一阶收敛。这些进步有效地减少了由不均匀粒子分布引起的离散化误差。

然而，将离散化模型的准确性提高到二阶不可避免地会导致数值不稳定。Matsunaga等人[22]和Oger等人[23]指出，在高阶离散化方案中，粒子更容易聚集，导致由于粒子重叠而产生数值发散。Duan等人[24]分析了自由表面不稳定的原因，发现自由表面处的快速误差累积会导致单个粒子的物理量（例如速度和压力）与其邻近粒子的物理量显著偏离。这种累积的差异最终会引发粒子飞溅和计算终止。为了解决这些问题，Duan等人[24]引入了无量纲数< />_i、c_i和< />_i来增强自由表面粒子的识别。通过这种改进，可以消除自由表面的误差累积，从而确保计算的稳定性。然而，当自由表面发生复杂变形时，仍然会观察到非物理现象，如参考文献[24]中的图29所示。为了优化自由表面的粒子分布，Matsunaga等人[18]使用移动表面网格来表示自由表面边界，其中表面网格所围成的区域定义为液相，外部区域定义为气相。在计算中，通过求解应力平衡方程来确定表面网格节点的位置。虽然这种方法为自由表面提供了准确的边界条件，但在应用于三维模拟时变得极其复杂[25]。在Matsunaga的模型和张的方法[26]的基础上，Duan等人[25]建立了一种表面重构方法，其中重构曲线（或表面）被拟合到自由表面粒子的位置。然后使用表面法线粒子移动来调整自由表面粒子。总之，上述方法通过施加更准确的边界条件来确保稳定性。然而，自由表面及其附近的粒子仍然表现出偏分布，这可能导致较大的离散化误差，正如Duan等人[24]所分析的。

为了消除由自由表面粒子偏分布引起的模型误差，有必要补偿自由表面附近的粒子数密度（PND）不足。因此，虚拟粒子方案经常被应用于二阶离散化模型。与试图实现准确自由表面边界的稳定化方法相比，虚拟粒子方法可以从根本上解决误差累积问题，因为自由表面附近的粒子得到了完全补偿。正如Shibata等人[27]和Chen等人[28]指出的，虚拟粒子边界与幽灵粒子边界不同，因为虚拟粒子被包含在压力泊松方程（PPE）的计算中，但并不物理上位于模拟域内。此外，假设粒子分布均匀时，二阶校正矩阵成为一个单位向量。因此，即使没有定义位置信息，也可以将虚拟粒子应用于自由表面及其附近的粒子。然而，如图1所示，当粒子重叠或表现出偏分布时，均匀分布的假设就不再成立，即使修改后的PND等于初始PND也是如此。在这种情况下，二阶离散化模型需要虚拟粒子的位置信息。缺乏计算节点会阻止自由表面附近的粒子以二阶精度进行离散化，从而降低模型准确性。

为了解决由于重叠和偏分布的粒子导致二阶离散化无法有效实施的问题，本文提出了一种改进的稳定化方案，该方案借助表面网格辅助实现。使用辅助网格来增强压力稳定性是一种公认的策略，例如Wang等人[29]提出的背景网格（BM）方案。然而，BM方案旨在通过增强PPE中源项的连续性来减轻非物理压力噪声。本研究提出的表面网格（SM）方案则有所不同，它并不是为了修改PPE源项而设计的。相反，所提出的SM是一个局部且动态的网格，仅在自由表面附近每个时间步重新生成。SM的唯一目的是为二阶校正矩阵提供所需的位置信息，否则在粒子不足的区域，校正矩阵将无法完成计算。因此，自由表面及其附近的粒子可以借助这个完成的网格使用二阶校正矩阵进行离散化。通过对粒子与网格比例进行了敏感性分析，以最小化PND与标准PND的偏差。通过模拟基准案例并将结果与基于虚拟粒子的方法进行比较，证明了所提出方法的准确性和稳定性。

小节片段

控制方程和标准MPS离散化

不可压缩粘性流体流动的控制方程是连续性方程和纳维-斯托克斯方程。在拉格朗日框架中，这些方程表示如下：$\frac{D\rho }{Dt}=0$$DuDt=?\frac{1}{\rho }?P+\upsilon {?}^{2}u+g$其中ρ是流体密度，t是时间，u是速度向量，P是压力，υ是运动粘度，g是重力加速度向量。物质导数D/Dt表示跟随移动流体粒子的时间变化率。

改进的二阶离散化方案的验证

通过四个数值基准测试验证了所提出方法的准确性和稳定性：静止水池中的压力分布、溃坝流、椭圆滴体的演化和方形滴体的旋转。为了在比较原始结果和改进结果时保持清晰，以下图中故意省略了表面网格和SPP的可视化。

结论

本文提出了一种基于表面网格的改进稳定化方案，用于自由表面及其附近粒子的二阶离散化。表面网格补偿了这些粒子的PND不足，并为缺失的邻居粒子提供了位置信息。包括梯度模型、拉普拉斯模型和散度模型在内的粒子相互作用模型通过考虑表面网格的贡献并加以改进。

CRediT作者贡献声明

李根：撰写 – 审稿与编辑、监督、方法论、资金获取、概念化。廖云龙：撰写 – 原始草稿、软件、调查、数据管理。姚培涛：可视化、验证、形式分析。