用于求解偏微分方程（PDEs）的最小二乘有限体积粒子方法

《Computers & Mathematics with Applications》：Least square finite volume particle method for solving PDEs

　　张震|刘晓星
中法核工程与技术研究所，中山大学，中国广东省珠海市

**摘要**
尽管通过将泰勒级数展开与最小二乘技术结合，无网格粒子方法可以实现任意阶精度，但由于高阶矩阵的求逆，这种方法会带来较高的计算成本。在本研究中，我们开发了一种名为最小二乘有限体积粒子（LSFVP）的新无网格方法，用于高效求解偏微分方程（PDEs）。粒子在二维空间中被建模为正方形。使用LSFVP框架对每个有限体积粒子的PDEs的积分形式进行离散化。原始的二阶导数在概念上被转换为一阶导数。最小二乘方法用于估计粒子表面的通量，而高斯积分确保了评估表面通量积分时的整体二阶精度。通过几个数值示例验证了所提出的LSFVP方法。

**引言**
无网格方法在无需生成网格的情况下求解PDEs方面具有优势。典型的无网格方法包括平滑粒子流体动力学（SPH）方法[1]、广义有限差分方法（GFDM）[2]、径向基函数（RBF）方法[3]、再生核方法（RKM）[4]和移动粒子半隐式（MPS）方法[5]。在这些方法中，SPH和MPS已广泛应用于工程分析，特别是在机械、海洋和核工程领域[[3], [4], [5], [6], [7], [8]]。广义有限差分方法（GFDM）也被广泛用于求解各种类型的PDEs。该领域的最新发展包括将Krylov延迟校正技术与GFDM结合，用于稳定的长时间动态裂纹分析[9]，以及将这一高阶框架扩展到模拟动态载荷下双材料介质中的界面裂纹[10]。此外，Benito等人[11]利用GFDM求解抛物线和双曲线方程。Ure?a等人进一步将GFDM的应用扩展到对流-扩散方程[12]和二阶非线性抛物线PDEs[13]。

当粒子分布变得不规则时，MPS和SPH方法的精度和可靠性会显著降低，离散化误差也会大幅增加。在无网格移动粒子方法领域，高阶梯度和拉普拉斯模型的发展尤为重要。已经开发了几种技术来提高粒子方法的精度，包括校正矩阵方法[[14], [15], [16], [17], [18]]、最小二乘方法[19,20]和成对松弛方法[21,22]。无网格粒子方法的其他发展还包括高阶交错公式，这些公式将拓扑梯度算子与虚拟对偶网格上的广义移动最小二乘（MLS）散度结合，以实现椭圆问题的兼容离散化[23]，以及稳定的拉格朗日-欧拉配置方案，这些方案明确分析和强制局部或全局守恒，以提高不可压缩流动模拟的精度[24]。其中，最小二乘移动粒子半隐式（LSMPS）方法因其能够以任意精度离散化微分算子而脱颖而出[[25], [26], [27], [28]]。

LSMPS基于泰勒级数展开。通过对固定位置附近粒子的标量变量进行泰勒级数展开，可以获得一系列泰勒级数展开公式。然后可以通过最小化误差函数来获得误差函数，并在固定位置导出导数向量[20]。在这方面，LSMPS的数学公式与GFDM有显著相似之处。这两种方法都依赖于局部泰勒级数展开和加权最小二乘最小化来重构不规则点云上的微分算子，主要区别在于它们的具体实现细节和历史发展路径。LSMPS包括A型LSMPS和B型LSMPS。A型LSMPS假设固定位置的函数值已知，只有导数未知；而在B型中，函数值和导数都假设为未知。如果泰勒级数多项式展开到二阶项，则LSMPS中的二阶导数模型仅具有一阶精度。这涉及到在二维空间中求逆5阶矩阵（B型为6阶矩阵），在三维空间中求逆9阶矩阵（B型为10阶矩阵）。为了在LSMPS中获得二阶精确的二阶导数模型，泰勒级数多项式需要展开到三阶项。然而，高阶矩阵的求逆需要更多的计算成本和更多的邻居粒子。

通过利用有限积分和有限体积方法中固有的积分过程，原始的二阶导数可以在概念上转换为一阶导数。基于有限体积的类似想法已在多个粒子研究中采用[[29], [30], [31]]。尽管如此，这些有限体积粒子方法的精度与传统的SPH和MPS相当。本研究的目的是基于有限体积粒子概念开发一种精确的无网格方法。假设每个粒子在二维空间中占据一个正方形。通过采用高阶一阶导数模型来估计粒子表面的通量。显然，采用高斯积分确保了评估表面通量积分时的整体二阶精度。

与GFDM等直接离散化微分算子的强形式无网格方法不同，所提出的LSFVP方法基于控制方程的积分形式。LSFVP相对于GFDM的主要优势在于通过积分降低了导数的阶数。通过将二阶导数的体积积分转换为一阶导数，LSFVP避免了直接重构高阶导数，从而使用低阶多项式近似更容易实现高阶精度，从而提高了解决方案的精度和鲁棒性。然而，该方法也存在某些局限性。首先，虽然LSFVP解决了PDEs的守恒积分形式，但它不具备基于网格的有限体积方法中的严格离散守恒属性，因为缺乏网格拓扑结构导致共享边界上的通量无法严格抵消。其次，当前的公式和验证主要针对二维问题。

本文的其余部分组织如下。第2节概述了LSMPS方法，提供了A型LSMPS和B型LSMPS的详细信息。第3节介绍了本研究开发的最小二乘有限体积粒子方法。第4节展示了多个数值测试的数值结果。第5节总结了关键发现。

**LSMPS方法**
LSMPS由Tamai和Koshizuka[20]开发，理论上可以获得任意阶精度。LSMPS包括A型和B型方案。在A型方案中，围绕相邻粒子的标量变量的泰勒级数展开表示为：
$$
\phi_j = \phi_i + \sum_{m=1}^p \frac{1}{m!} \left( r_j \cdot \nabla \right)^m \phi_{|r_i + o(|r_j|p+1)}
$$
忽略高阶项，可以得到以下方程：
$$
\phi_j - \phi_i = \sum_{m=1}^p \frac{1}{r_m} \left( r_j \cdot r_s \cdot \nabla \right)^m \phi_{|r_i} = \pi_j \cdot \left( H_r s^{-1} \Lambda \phi_{|r_i} \right)
$$
其中 $\phi_i = \phi(r_i)$，$\phi_j = \phi(r_j)$，$r_j = r_j - r_i$，$r_s$

**最小二乘有限体积粒子（LSFVP）方法**
在LSFVP中，假设粒子在二维空间中为正方形。图1展示了这一概念。假设初始粒子间距为 $l_i$，控制体积是一个边长为 $l_i$ 的正方形。控制体积的表面包括正方形的四条边，分别用N、S、W、E表示。使用LSMPS估计每条边上的通量。需要注意的是，B型...

**数值测试**
为了定量评估精度，定义了离散L2误差范数：
$$
\|E\|^2 = \frac{1}{N} \sum_{i=1}^N \left( u_{num} - u_{exact} \right)^2
$$
其中 $N$ 是粒子总数。收敛阶数通过以下公式计算：
$$
\text{Order} = \frac{\ln\|E\|^2_{h_1}}{\ln\|E\|^2_{h_2}} = \frac{\ln(h_1)}{\ln(h_2)}
$$
其中 $h_1$ 和 $h_2$ 表示两个连续细化级别的粒子间距。

为了确保在不同流动状态下的数值稳定性，时间步长 $\Delta t$ 必须满足对流的Courant-Friedrichs-Lewy（CFL）条件和稳定性准则...

**结论**
在本文中，我们提出了LSFVP方法，这是一种新颖的无网格有限体积框架，它使用B型LSMPS算子在控制体积的边缘中点重构解和梯度值。通过单次最小二乘求解恢复边缘中点值，并将未知数限制在粒子体积中心，LSFVP方法实现了二阶空间精度，而无需引入额外的自由度。我们使用四个测试问题来验证该方法的空间精度...