带有微分再生核的通量平衡平滑粒子流体动力学方法：多相流无网格相场建模

《Engineering Analysis with Boundary Elements》：Flux-balanced Smoothed Particle Hydrodynamics with differential reproducing kernels for meshfree phase-field modeling of multiphase flow

　　

摘要

本文提出了一种通量平衡的光滑粒子流体动力学（FBSPH）公式，用于多相流的无网格相场（PF）建模，其中成对反对称数值通量确保离开一个粒子的通量恰好等于进入其邻居粒子的通量。此外，使用微分再生核（DRK）来构建更高阶一致性精度核函数。FBSPH被应用于不混溶两相和三相流的保守Allen-Cahn相场方程，适用于固定和移动粒子。在该框架中，扩散通量和对流通量分别采用对称内部惩罚（SIP）公式和本地Lax-Friedrichs黎曼求解器进行处理。基准测试表明，通量平衡在固定粒子上可将相场序参量的守恒精度保持在接近机器精度。在移动粒子上，与使用相同核函数的标准SPH相比，通量平衡显著降低了守恒误差和相间质量传递。不同核函数阶次的对比显示，守恒性由通量处理方式决定，而精度则由核函数一致性决定。所提出的无网格PF公式具有伽利略不变性，能够处理高密度比，简化了三元体系中的表面张力施加，并抑制了游离粒子的形成。本文对全面的一系列基准测试给出了数学公式和数值结果，并进行了讨论。

引言

相场（PF）方法已成为一种强大的弥散界面技术，用于多相流的数值模拟，特别是在涉及复杂界面动力学的问题中。在该框架中，引入连续的相场序参量来区分不同流体相，界面被表示为序参量平滑变化的有限厚度区域。这种公式使得相属性在界面处能够平滑过渡，与尖锐界面模型[1]、[2]形成对比。

相场序参量的演化通常由四阶Cahn-Hilliard（CH）方程[3]、[4]或二阶Allen-Cahn（AC）方程[5]控制。两个方程均源自Ginzburg-Landau型自由能泛函，但在动力学上存在根本差异：CH方程对应于H?1梯度流，本质上对序参量具有全局保守性；而经典AC方程代表L2梯度流，除非经过修改，否则不具有保守性。

各种数值方法已被用于求解相场控制方程，包括有限差分方法[6]、[7]、[8]、[9]、[10]、[11]、[12]、[13]、[14]、[15]、[16]，有限元方法[17]、[18]、[19]、[20]、[21]，间断Galerkin（DG）方法[22]、[23]、[24]、[25]以及格子Boltzmann方法[26]、[27]、[28]。这些方法的共同特点是依赖于网格或晶格来显式定义离散节点之间的连接关系，并利用该关系构建求解控制方程所需的插值函数。

与基于网格的技术不同，无网格方法能够在散点集上求解偏微分方程，而无需显式定义连接关系。无网格方法在相场建模中的应用是相对较新的进展。已有文献[29]、[30]提出了一种采用局部最大熵基函数的拉格朗日无网格框架，用于生物膜动力学的相场模拟。其他研究利用径向基函数（RBFs）[31]、[32]、[33]、[34]、[35]、[36]和移动最小二乘（MLS）逼近[37]来构建无网格相场求解器。其他公式包括粒子差分法[38]和再生核粒子法（RKPM）[39]。

光滑粒子流体动力学（SPH）是一种替代的无网格方法，被广泛认为是最早的无网格方法，最初由Lucy[40]在20世纪70年代开发，Gingold和Monaghan[41]独立进行了相同工作，用于模拟天体现象。由于历史悠久，SPH是最成熟的无网格方法之一，因其局部化公式、易于并行化以及处理高度复杂几何形状的能力而广受欢迎。重要的是，其数学基础具有惊人的普适性，使其适用于有界或无界域，以及从介观到天体物理尺度。

标准SPH过程采用基于配点的方法，通过在核K_i中使用核加权求和直接近似粒子i处的场变量u的导数，来求解场变量u演化的控制方程。例如，场变量u在粒子i处的一阶梯度?u_i可近似为?u_i ≈ Σ_{j∈K_i}(u_j - u_i)?w_{ij}V_j，拉普拉斯算子通常近似为[42] ?2u_i ≈ Σ_{j∈K_i}(u_j - u_i)·2r_{ji}/|r_{ji}|·?w_{ij}V_j，其中V_j是邻近粒子j的体积，r_{ji} = x_i - x_j，?w_{ij}是核加权函数的梯度。

经典SPH使用反对称的核函数（?w_{ij} = -?w_{ji}）。具有该特性的核函数满足F_{ij} = -F_{ji}，其中F_{ij}表示粒子对(i, j)相关的通量。换言之，离开粒子i的通量被保守地传递给其邻近粒子j。然而，这通常以牺牲一致性为代价，经典核函数通常不能满足零阶一致性[43]、[44]。该性质对于提高精度和改善收敛性至关重要。尽管如此，一般SPH公式对核函数的选择施加的限制很少，除了一致性分解等基本性质外。这种灵活性激发了大量一致性修正变体的发展，包括修正SPH方法[43]、[45]，有限粒子方法[46]，核梯度无关方法[47]以及拉格朗日差分动力学[48]。通过引入径向基函数（RBF）[49]、[50]、移动最小二乘（MLS）[44]、[49]、[50]和再生核（RK）[51]、[52]进一步实现了改进。这些方法可被视为原始SPH方法的扩展，主要区别在于其核函数的构建或修正方式。

追求高阶一致性的一个主要缺陷是可能导致反对称性质的丧失。虽然强制执行该性质对局部守恒至关重要，但传统的"对称化"算子方法，例如显式平均核函数，通常会降低底层逼近的多项式再现性质。这导致一致性和精度的损失，特别是在核支撑被截断的边界附近。因此，大多数SPH公式通常以满足一个性质为代价牺牲另一个性质。Dilts[44]、[53]、Frontiere等人[54]、Vila[55]以及Lanson和Vila[56]、[57]的工作更详细地讨论了粒子云中一致性与守恒性之间的固有张力。

广泛用于流体物理模拟的有限体积（FV）方法通过在单元间界面进行通量平衡来强制执行局部守恒，确保离开一个单元的通量等于进入与其共享边或面的邻居单元的通量。在将基于有限体积的创新引入粒子方法方面已取得了若干重要成就。这包括Monaghan[58]的工作、Gudonov-SPH（GSPH）[59]、WENO-SPH[60]、MUSCL-SPH[61]、光滑体积积分（SVI）方法[62]、有限体积粒子法（FVPM）[63]、[64]、[65]以及成对松弛无网格方法[66]、[67]。此类方法的详细综述可参见文献[68]及其参考文献。它们的主要共同特征是采用了某种形式的粒子间通量平衡，类似于有限体积方法中那样，使得通量在粒子之间被保守地传递。这些方法主要开发用于激波捕获应用，而在固定和移动粒子分布上利用高阶核函数守恒多相流相场序参量方面的应用则非常有限（如果有的话）。

SPH已在多项工作中[69]、[70]、[71]、[72]、[73]用于多相流的数值模拟。由于其能够在无网格条件下追踪物质界面，因此天然适用于多相流。区分流体相的最常见方法是使用颜色函数，其中每个粒子被标记有相编号标识符。这实质上假设了尖锐的不连续界面，使得依赖颜色场梯度和曲率的表面张力力计算更容易受到噪声影响，并对粒子无序更加敏感。为计算表面张力，已采用RBF[74]和等值面算法[75]来近似相界面。另一种方法是在界面附近平滑颜色函数[72]、[73]。这种纯数值平滑缺乏热力学恢复力。因此，它无法防止"游离"粒子穿过界面迁移——即某一相的粒子游离到另一相中，导致非物理的相污染、数值不稳定以及界面区域噪声加剧，从而复杂化表面张力的施加。

另一方面，相场方法通过系统自由能的最小化来建模界面运动。因此，它可以立足于问题本身的物理基础。这使得表面张力效应和复杂界面动力学可以直接纳入相场公式中，无需显式的几何界面追踪或表面重构。此外，相场序参量自然平滑了界面处流体属性的过渡，有助于提高相界面处的数值稳定性。更重要的是，相场自由能泛函中固有的热力学恢复力主动驱使序参量返回其平衡体相值，防止游离粒子形成并显著减少界面区域噪声。将相场模型与SPH耦合的另一个优势在于，它既适用于欧拉固定粒子配置，也适用于拉格朗日移动粒子配置。因此，可以使用固定粒子进行多相模拟，这通常比拉格朗日方法在计算上高效得多。使用颜色函数时不可能实现这一点，因为颜色函数仅依赖于移动粒子框架。

相场方法已在多项SPH研究[49]、[52]、[76]、[77]、[78]、[79]中被采用，所有这些研究均使用标准核加权配点SPH方法求解控制相场方程。此外，关于使用无网格方法进行涉及超过两种共存流体的多相流相场建模的研究仍然非常有限。值得注意的是，文献[52]、[80]中进行了涉及三元流体的无网格相场模拟。在文献[52]中，采用增强再生核（RK）和微分再生核（DRK）的SPH求解了Cahn-Hilliard-Navier-Stokes（CHNS）方程的三元体系。该工作中注意到一个重要的观察结果：尽管CH方程在数学上本质上是守恒的，且使用了至少一阶可再现的核函数，但在核加权SPH-RK/DRK框架中仍观察到相场守恒误差。

本文的主要动机是改善在粒子云上使用高阶核函数求解相场方程时的守恒性能。这使得相场模型能够在拉格朗日和欧拉SPH框架中均被使用，而无需颜色函数或相界面重构。我们的主要贡献是开发和应用了FBSPH公式，用于不混溶两相和三相流CAC方程的无网格求解。为实现此目标，用反对称粒子间通量平衡取代了求解控制相场方程的标准SPH核加权配点方法，确保离开一个粒子的通量恰好等于进入其邻居的通量。同时，微分再生核被用来确保核函数至少具有第一阶一致性。

这种构建的本质特征在于守恒性和一致性被分别施加在两个不同的对象上。守恒所需的反对称性被施加在粒子间接口的数值通量上，而非核梯度上。因此，再生核保持完全不被修改，与历史上会降低多项式再现性质的核对称化策略形成对比。

在保守光滑粒子公式中使用再生核已在保守再生核SPH（CRKSPH）方法[54]中提出，其中再生核与反对称粒子间公式耦合以强制守恒质量和动量。然而，据我们所知，尚无先前工作建立了通量平衡SPH框架来专门求解涉及超过两种不混溶流体相的移动粒子上多相流的相场方程。为此需要克服多项挑战。首先，弥散界面区域对粒子无序高度敏感。其次，数值方案不应引入过度的数值耗散而使界面区域非物理地展宽，使其在整个模拟过程中保持紧凑且厚度相对恒定，即使在移动粒子条件下也是如此。与CRKSPH中提出的公式不同，本文通过在类似于有限体积和间断Galerkin公式中使用反对称数值通量近似守恒方程的面积积分来强制局部守恒。这允许引入黎曼求解器来处理迎风稳定化和高阶重构。

在这项工作中，采用了一种用于守恒定律的单调上游居中格式（MUSCL）并结合斜率限制器，来重构粒子界面两侧的场变量。此外，扩散通量使用对称内部惩罚（SIP）格式进行处理，而对于固定粒子的情况，则使用局部Lax–Friedrichs黎曼求解器进行迎风稳定化。

本文采用了一种分裂求解器策略：相场（PF）方程使用FBSPH求解，而Navier–Stokes（NS）方程则使用核加权SPH-DRK方法求解[52]。两种求解器运行在相同的粒子分布上，并共享相同的再生核，确保了耦合系统中几何表示的一致性。由于NS方程是使用基于配点的核加权离散方法求解的，其中DRK核在不规则粒子分布上并非严格反对称，因此离散的压力和粘性力不是成对反对称的。结果是，总动量并不严格守恒。

在FBSPH和SPH-DRK两种配置中，对NS方程均使用相同的核加权SPH-DRK求解器，这是一种刻意的方法论选择。这确保了在PF守恒性能、收敛速度或界面保持方面观察到的任何差异，仅归因于PF求解器的FBSPH格式，这正是本文的核心论点和焦点，而不会受到流体求解器差异的干扰。虽然自20世纪90年代以来，NS方程的反对称通量格式已经得到了广泛开发，但相比之下，移动粒子上PF序参量的守恒几乎未受到关注，而这正是本文所解决的空白领域。

本文其余部分组织如下：首先介绍基于PF的N相系统的控制方程，随后详细阐述FBSPH格式的推导、实现细节以及在拉格朗日和欧拉设置下的一系列数值验证。

**章节片段**

**用于N相系统的守恒Allen–Cahn相场模型**

对于包含两种不混溶流体的域 ?? ? ???（??=2,3），令 ?? 为一个在 ??? 中取值为 ??=0、在 ??? 中取值为 ??=1 的平滑变化的PF序参量。弥散界面定义为 ??={??∈??∣???(??,??)=0.5}。相关的Ginzburg–Landau自由能泛函为 ??(??)=∫??(1??2???(??)+1/2?|???|2)????。

这些值之间过渡区的宽度由界面厚度的毛细宽度 ?? 来表征。项 ??(??)=1/2???2?(1???)2 是一个双阱函数，在体相流体中具有最小值 ??=0 和 ??=1。

**通用格式**

考虑一个计算域 ?? ? ???，其中 ??=2,3。定义场变量 ???(??,??) 的瞬态输运方程残差为 ??????/?????????，其中 ?? 为时间，?? 为通量矢量，由式(14)对CAC方程进行表达。将残差 ? 乘以试验函数 ??，在 ?? 上积分，并对通量项进行分部积分，得到 ∫?????????/????d??+∫?????????d????∫_???????????d??=0。

我们将 ?? 划分为一组粒子 ??，每个粒子定义一个控制体积 ????，使得 ??=∪_{??∈??}????，其中 ?? 为粒子总数。

**粒子偏移技术**

采用粒子偏移技术（PST）[115],[116],[117] 以避免SPH模拟中因粒子聚集或团聚而引发的众所周知的不稳定性[92],[118],[119],[120]。PST的主要思想是，粒子聚集区域的粒子浓度较高。因此，如果实施菲克定律，可以将粒子略微向低浓度区域偏移，以缓解聚集现象。这通常通过使用扩散方程来完成。

**纯对流下的数值耗散**

在我们继续求解PF方程之前，首先研究FBSPH方法求解不含CAC方程中扩散/反扩散项的纯对流问题。纯对流方程为 ????/????+??(?????)=0。

在尺寸为 ??=[?1,1]×[?1,1] 的计算域 ?? ? ?2 上进行了一系列数值实验，以评估FBSPH方法在界面平流问题上的性能。该域使用10,000个粒子以结构化方式离散化。

**结论**

在本工作中，我们证明了对于高阶再生核，用于求解PF方程的标准核加权SPH方法会导致PF序参量次优的守恒性能。为弥补这一缺陷，我们提出了一种用于求解移动和固定粒子上PF方程的通量平衡SPH格式。标准的核加权方法被替换为反对称的粒子间通量平衡，确保离开某粒子的通量恰好等于相邻粒子之间的通量。

**利益冲突声明**

作者声明不存在可能影响本文所报告工作的已知竞争性经济利益或个人关系。

**致谢**

作者衷心感谢匿名审稿人花费时间仔细审阅本手稿。他们富有洞察力的意见和建议帮助提高了论文的质量。

Adam Y. Ghoneim