基于局部泰勒展开的高精度无网格算法

High-Precision Meshless Algorithm Based on Local Taylor Expansion

摘要: 光滑粒子动力学方法(Smoothed Particle Hydrodynamics, SPH)所具备的无网格特性使其摆脱了网格划分过程，在处理复杂几何边界上具有潜在优势。然而，传统SPH受限于空间微分高阶一致性缺失，无法实现高精度求解。本文基于SPH粒子的局部泰勒展开，通过邻域粒子空间矩阵求解，实现了满足高阶一致性的SPH微分算子算法。通过一、二维函数证明了该算法在任意粒子分布中微分2、3、4阶一致性。随后，基于热传导方程、Burgers方程、纳维–斯托克斯方程等偏微分方程证明了该算法求解空间算子的精确性与稳定性。基于局部泰勒展开的SPH算法为无网格方法在复杂几何域的高精度数值求解提供了一种可行方法。

Abstract: The meshless nature of Smooth Particle Dynamics (SPH) method eliminates the need for mesh partitioning and has potential advantages in handling complex geometric boundaries. However, traditional SPH is limited by the lack of high-order consistency in spatial differentiation and cannot achieve high-precision solutions. This article is based on the local Taylor expansion of SPH particles, and solves the SPH differential operator algorithm that satisfies high-order consistency through the neighborhood particle space matrix. The algorithm has been proven to have 2nd, 3rd, and 4th order consistency in differentiation in any particle distribution through one-dimensional and two-dimensional functions. Subsequently, the accuracy and stability of the algorithm for solving spatial operators were demonstrated based on partial differential equations such as the heat conduction equation, Burgers equation, and Navier Stokes equation. The SPH algorithm based on local Taylor expansion provides a feasible method for high-precision numerical solutions in complex geometric domains.

1. 引言

高精度的偏微分方程数值解是解决诸多科学问题的核心需求。基于网格的数值方法长期占据主导地位，对于简单几何域求解而言，诸如十阶有限差分(FD) [1]或指数收敛谱方法[2]能够实现极高的数值精度。然而面对复杂几何域的求解，基于正交网格的高阶差分不再适用，全局的谱方法不再可行。同时，求解精度具有网格依赖性，高质量的网格划分需要昂贵的计算成本。

此时，能够实现计算域快速离散的无网格方法体现了巨大的应用潜力。无网格方法中自动生成计算节点的过程相对直接，同时不需要计算与存储粒子间额外的连通信息，能够大幅提升计算效率。光滑粒子动力学(Smoothed Particle Hydrodynamics)作为一种典型的无网格方法最初为天体物理模拟而开发的[3] [4]，并在不可压缩流动、自由表面流动、流体结构相互作用和固体力学的模拟中取得了巨大成功[5]。在光滑粒子动力学中，流体粒子的性质和空间微分由其周围的领域粒子的加权求和所决定。权重值一般由多项式函数、高斯函数等光滑核函数所确定[3] [6]。然而，无论选取何种形式的光滑函数，传统光滑粒子动力学在空间微分理论上仅具有2阶一致性[7]。在粒子非正交排布状态下空间微分甚至无法达到1阶一致性，并且收敛速率取决于核函数[8]，适用于广泛使用的Wendland C2核[9]，该方法在N-5/2的情况下接近2的理论收敛速度(以h为单位)。

为了改善由于粒子非正交排布所导致的空间微分一致性缺失，大量基于光滑粒子动力学的修正算法不断被提出[10]。重构核粒子方法(RKPM) [11]和修正光滑粒子动力学(CSPH) [12]通过在核函数中引入修正因子以确保微分一致性。Zhang与Liu独立开发了一种对SPH插值的校正，称为有限粒子法(FPM) [13] [14]。通过包含二阶和高阶核导数，并求解线性系统以获得函数及其梯度的近似来实现高阶。Liu提出了一种恢复粒子一致性的方法[15]。方法保留了传统的非负平滑函数，而不是重建新的平滑函数，提高了近似精度。Asprone进一步改进了FPM，利用多个单项式的线性组合构建了新核函数，并证明了其具有良好的二阶一致性[16] [17]。

仅实现微分的2阶一致性无法满足高精度数值求解的需求。因此，Lind和Stansby通过核函数的系数组合成功实现了微分的高阶一致性，建立了局部各向异性基函数法(LABFM)，为无网格方法在高精度数值计算的应用奠定了基础[18]-[20]。然而，这种基于粒子局部泰勒展开的高阶算法仍然存在一些不足。在求解核函数系数矩阵的过程中存在取决于基函数的矩阵病态问题，对最终系数的求解精度有着直接作用，这对于基函数的形式有着严格的先验要求。本文基于LABFM构建了一种不依赖于基函数形式的高精度粒子算法，通过对于领域粒子权重的直接矩阵求解获得了高阶一致性的空间微分算子。正交与非正交粒子分布下，该方法在一维、二维函数求导，偏微分方程求解均被证明能够实现2，3，4阶一致性，为无网格方法的高精度数值计算提供了一种全新的方案。

2. 数值模拟和误差分析

算法原理

以一维计算域为例，x_i是粒子i的位置信息，将粒子i周围存在N个邻域粒子j，作为粒子i的空间导数加权支撑粒子集合，h_ij表示粒子i与粒子j间距离。此时假定粒子i的n阶空间导数等于邻域粒子j的加权求和，即：

$f^{\left(n\right)}\left(x_i\right)={\displaystyle \sum _{j=1}^{N}f\left(x_j\right)W_{ij}};n=0,1,2\cdots ,k$ (1)

其中W_ij代表粒子i与j之间的领域粒子权重集合。随后将f(x_j)在i点进行泰勒展开，可得：

$\begin{array}{l}{f}^{\left(n\right)}\left(x_i\right)=f(x_i){\displaystyle \sum _{j=1}^N{W}_{ij}}+f^{\prime }(x_i){\displaystyle \sum _{j=1}^N(x_j-x_i)W_{ij}}+\cdots +\frac{1}{n!}{f}^{\left(n\right)}\left(x_i\right){\displaystyle \sum _{j=1}^N{\left(x_j-x_i\right)}^n{W}_{ij}}+R_n\left(x_j\right)\\ R_n\left(x_j\right)=\frac{1}{\left(n+1\right)!}{f}^{\left(n\right)}\left(\xi \right){\left(x_j-x_i\right)}^{n+1}\end{array}$ (2)

为使得f^(m)(x_i)在等式左右严格相等，必然存在W_ij集合满足以下线性方程：

$\left(\begin{array}{cccc}1& 1& \cdots & 1\\ x_1-x_i& x_2-x_i& \cdots & x_n-x_i\\ {\left(x_1-x_i\right)}^2& {\left(x_2-x_i\right)}^2& \cdots & {\left(x_n-x_i\right)}^2\\ \begin{array}{l}\cdots \\ {\left(x_1-x_i\right)}^m\\ \cdots \\ {\left(x_1-x_i\right)}^n\end{array}& \begin{array}{l}\cdots \\ {\left(x_2-x_i\right)}^m\\ \cdots \\ {\left(x_1-x_i\right)}^n\end{array}& \begin{array}{l}\cdots \\ \cdots \\ \cdots \\ \cdots \end{array}& \begin{array}{l}\cdots \\ {\left(x_n-x_i\right)}^m\\ \cdots \\ {\left(x_1-x_i\right)}^n\end{array}\end{array}\right)\left(\begin{array}{l}{w}_{i1}\\ w_{i2}\\ \cdots \\ \cdots \\ w_{im}\\ \cdots \\ w_{in}\end{array}\right)=\left(\begin{array}{l}0\\ 0\\ 0\\ \cdots \\ m!\\ \cdots \\ 0\end{array}\right)$ (3)

基于该线性方程组能够解得W_ij集合的唯一解，实现f^(m)(x_i)的精确求解。实际中f^(m)(x_i)的一致性取决于选取的邻域粒子j的个数N，对于一维问题而言，m阶一致性需要N个邻域粒子，二维问题与三维问题需要考虑偏导项，因此邻域粒子个数N与一致性阶数关系如下：

$\begin{array}{l}N=\frac{m^2+3m}{2};m=1,2\cdots ,k\\ N={\displaystyle \sum _{m=2}^{m=k+1}\left(\frac{m^2+m}{2}\right)};m=1,2\cdots ,k\end{array}$ (4)

基于对于领域粒子与计算粒子位置信息的构建能够实现光滑粒子动力学中的空间微分算子的高精度求解。显然，系数矩阵的状态对于求解精度有直接影响，因此，本文分别考虑正交与非正交两种粒子分布状态，其中非正交粒子由正交位置叠加随机位置扰动确定，噪声强度由δ决定，如图1所示，展示了δ = 0.5的粒子分布情况，其中横轴和纵轴坐标单位为同一度量，单位为1。在这种分布情况下，能够很好的模拟在实际问题中的粒子分布，并且噪声强度在0.5的大小是以往研究中未曾有过的，本文考虑这种随机的扰动，使算法面临的随机性增大，其结果更具有可靠性。本文的所有算例均以δ = 0.5的粒子分布为基础，而正交的粒子分布状态下，结果相应比非正交时更优异，不再展示。

本文采用一维、二维函数与偏微分方程为参考验证该高精度微分算子算法的2，3，4阶一致性，其中h代表目标粒子i与邻域粒子j之间的平均间距，H代表计算域特征长度。采用L2均方根误差(Root Mean Square Error, RMSE)表征计算值与解析值间的误差，即：

$L2-norm{(\cdot )}^d=\frac{{\left({{\displaystyle \sum \left[L_i^d(\cdot )-C^d\cdot D(\cdot )|{}_i\right]}}^2\right)}^{\frac{1}{2}}}{{\left({\displaystyle \sum C^d\cdot D(\cdot )|{}_i^2}\right)}^{\frac{1}{2}}}$ (5)

Figure 1. Particle distribution map in two-dimensional domain

图1. 二维域内粒子分布图

3. 算例数值分析

下面将从一维函数和偏微分方程，二维函数和偏微分方程的角度，通过测试在非正交分布时h/H变化过程中误差的收敛情况验证微分算子的2，3，4阶一致性。

3.1. 一维函数

考虑x∈[−H, H]的测试域，H = 1。定义测试函数为正弦函数：

$y=\sin 2\pi x$ (6)

其一阶导数形式如下：

$y=2\pi \cos 2\pi x$ (7)

·黑色实线为函数的真实曲线，红色散点为数值解，如图2所示：

Figure 2. Trigonometric function images

图2. 三角函数图像

其对应的二阶，三阶和四阶收敛阶如图3所示，其中虚线表示对应收敛阶的理论斜率下降趋势。测试结果显示其理论收敛率均达到预期效果。在某些特征宽度上，会出现波动情况，整体趋势是符合相应理论收敛要求。

Figure 3. Convergence order of trigonometric functions

图3. 三角函数收敛阶

3.2. 一维偏微分方程

Burgers方程是作为一个非线性偏微分方程，结合了对流和扩散效应，并且在流体动力学和输运理论中有着广泛的应用。Burgers方程可以看作是线性对流方程(Advection Equation)和线性扩散方程(Diffusion Equation)的组合，其中包含了非线性项，这使得它能够描述更加复杂的物理现象。

Burgers方程形式很多，这里考虑最简单的一维不可压缩Burgers方程，其偏微分方程形式为：

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{{\partial }^{2}u}{\partial x^2}$ (8)

其中，u (x, t)是流体速度，ν是粘性系数，x和t分别是空间和时间变量。左边的第一项表示速度的对流(或迁移)效应，第二项表示速度的非线性对流效应，右边的项表示速度的扩散效应。设置初始条件为一个反弦波：

$y=-\sin \left(\pi x\right)$ (9)

测试域x∈[−H, H]，H = 1，ν = 0.01/π，dt为单次迭代时间步长1/T，时间迭代次数为T = 20,000，采用欧拉时间推进。采用周期性边界条件，使边界上的粒子拥有足够的支撑邻域粒子。当粘度较小时，Burgers方程的解可能会形成陡峭的冲击波或不连续点。这是因为低粘度无法有效地扩散由非线性项引起的急剧变化，为了解决这些数值稳定性问题，引入人工粘度$\epsilon$ = 0.005，同时在边界处施加零边界条件，保证边界处的数值稳定。

图4展示了在正交粒子分布条件下，当粘度增大到ν = 0.05/π，扩散效应变得更强，冲击波变得更加平滑和宽广，这是因为更大的粘度能够更有效地扩散速度场中的梯度。同时，解的振荡减少，增大粘度后，这些振荡会被粘性扩散效应抑制，解趋于更加平滑和单调。同时速度场中的高频成分会被更快地衰减，这些成分在高粘度下会被迅速扩散和消除，导致解更加光滑。此时正交分布状态下，数值解与解析解在各种粘性系数下均吻合良好。

Figure 4. Introducing artificial viscosity and Burgers equation images

图4. 引入人工粘度下的Burgers方程图像

Burgers方程在粘度为ν = 0.01/π同样噪声δ = 0.5影响下有相应的2，3，4阶收敛阶，如图5所示：

Figure 5. Burgers equation convergence order

图5. Burgers方程收敛阶

Burgers方程的L2误差相对于三角函数变大，但在空间域上仍然基本满足对应收敛速率。误差主要来源于时间维度上的一阶欧拉推进。从图5中可以看出，其收敛速率依然贴合对应的阶数。但会产生一定的数值波动，波动强度与噪声强度正相关。空间收敛速率整体趋势上仍然满足要求。然而，当计算域内粒子增多后，存在h/H的临界值。低于该临界值，公式3的系数矩阵愈发病态，无法精确实现权重集合W_ij的数值求解，导致相应空间微分计算误差增加。

随后，基于热对流方程进行验证。热对流方程是描述热量通过流体流动传递的数学方程。其偏微分方程形式如下：

$y=\tanh \frac{x-t}{5}$ (10)

同样设置测试域x∈[−H, H]，H = 1，固定边界条件，通过计算对应空间导数后迭代时间，可以得出在时间域下的图像，数值解与解析解对比如图6所示，对应收敛阶如图7所示。由此可以看出，一维热对流方程因其形式简单，算法收敛性相较Burger方程更好。

Figure 6. Thermal convection equation images

图6. 热对流方程图像

Figure 7. Convergence order of heat convection equation

图7. 热对流方程收敛阶

布鲁斯–格林纳克方程(Bruce Greenark equation)是流体动力学中描述湍流的一个模型方程。它考虑了湍流中的非线性相互作用和能量耗散。其偏微分方程形式如下：

$\frac{dx}{dt}=A-\left(B+1\right)x+x^{2}y$ (11)

$\frac{dy}{dt}=Bx-x^{2}y$ (12)

其中x和y分别代表反应物的浓度，A和B是方程中的常数参数。初始条件为x(0) = 1.0和y(0) = 2.0。布鲁斯–格林纳克方程通常用于描述化学反应动力学，而化学反应通常没有明确的空间边界。参数A = 1.0，B = 3.0，时间迭代次数T = 50。设置测试域为x ∈ [0, H]，H = 10，可以得到如图8所示的非正交粒子分布条件下的数值解与解析解对比图，对应收敛阶如图9所示：

Figure 8. Images of Bruce Greenark equation

图8. 布鲁斯–格林纳克方程图像

Figure 9. Convergence orders of Bruce Greenark equation

图9. 布鲁斯–格林纳克方程收敛阶

和Burgers方程类似，布鲁斯–格林纳克方程收敛性在加密粒子之后不会严格按照阶数方向收敛，甚至当δ较大的时候，出现误差相对增大的情况。对于此方程而言，临界值h/H比Burgers方程的更大，然而在临界值之前误差能够按对应收敛阶下降。

在测试的一维函数和偏微分方程中算法有着良好表现，下面将基于二维函数和偏微分方程测试算法的可行性和收敛性。

3.3. 二维函数

接下来，针对取二维抛物面函数进行测试，其原函数解析式为：

$f\left(x,y\right)=\frac{1}{3}{x}^3+\frac{1}{3}{y}^3$ (13)

抛物面形式为：

$\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}=x^2+y^2$ (14)

通过算法绘制出来的抛物面形式如图10所示，算法对二维微分同样适用。非正交状态下的微分收敛如图11所示。此时算法符合相应收敛速率，误差已达到计算机精度。

Figure 10. Images of paraboloid

图10. 抛物面图像

Figure 11. Convergence orders of paraboloid

图11. 抛物面收敛阶

3.4. 二维偏微分方程

接着基于典型抛物线方程——(齐次)非定常热方程上测试本算法。热传导方程是描述物体内部热量分布随时间变化的偏微分方程。它是一个基本的物理方程，在工程、物理和数学等领域都有广泛的应用。热传导方程的一般形式如下：

$\frac{\partial u}{\partial t}=\kappa {\nabla }^{2}u$ (15)

其中：u (x, t)表示温度分布函数，它是位置x和时间t的函数。$\kappa$ 是热扩散系数，通常表示材料传导热量的能力s。${\nabla }^{2}u$ 是拉普拉斯算子，表示温度的空间二阶偏导的和，描述了热量在空间中的分布情况。定义初始条件为以下表达式：

$u\left(x,y,t\right)=\sin \left(\frac{2\pi x}{H}\right)\sin \left(\frac{2\pi y}{H}\right)\exp \left(\frac{-8\kappa {\pi }^{2}t}{H^2}\right)$ (16)

用带噪声的笛卡尔分布离散域，进行时间积分，时间步长为dt = 0.05h × h/κ，设κ = 1。计算域H = 2。在周期边界条件下，所有计算粒子都有完整的邻域。当然在测试中发现，如果用齐次狄利克雷条件(边界上u = 0)代替周期边界条件，使得边界附近和边界上的粒子邻域不完整。狄利克雷条件(齐次或非齐次)通过对边界上的所有粒子规定u，并仅对内部粒子求解(15)来施加。当k ≤ 4时，方案是稳定的，能再次看到接近k的收敛速率。因此，对于非高阶收敛问题，即使是非周期性边界条件仍然适用。可以绘制出温度分布情况，如图12所示：

Figure 12. Temperature distribution

图12. 温度分布情况

对应的2，3，4阶收敛图如图13所示：

Figure 13. Convergence orders of heat conduction equation

图13. 热传导方程收敛阶

在二维的热传导方程中可以看到，其收敛的稳定性很好。即使随着时间域迭代，其收敛速率只有在低于临界值后会出现逐渐缓平，其误差还是相对很小。

最后，测试此方法求解不可压缩Navier-Stokes方程的效果。动量方程是非线性双曲–抛物方程，不可压缩性约束表现为椭圆PDE。Navier-Stokes方程可以用涡度来表示。泰勒–格林涡(Taylor-Green vortex)是一种在流体力学中研究的涡流现象，它描述了一个周期性的、二维的、不可压缩的流体在低雷诺数下的流动。其偏微分方程形式为：

$\frac{\partial \omega }{\partial t}+u\nabla \omega =\frac{1}{\mathrm{Re}}{\nabla }^2\omega$ (17)

需要初始化速度场：

$u\left(x,y,t=0\right)=\sin \left(x\right)\cos \left(y\right)\exp \left(-2\nu t\right)$ (18)

初始化温度场：

$T\left(x,y,t=0\right)=-\sin \left(x\right)\sin \left(y\right)$ (19)

初始化涡度场：

$w\left(x,y,t=0\right)=-\sin \left(x\right)\sin \left(y\right)\exp \left(-2\nu t\right)$ (20)

初始条件均取t = 0时刻，设置周期性边界条件，时间域不断迭代下一时刻的值，时间步长为dt = 0.05h × h/κ，κ = 1。动力粘度v = 0.1，雷诺数Re = 1.0。

Figure 14. Vorticity and velocity field at time t = 0

图14. t = 0时刻涡度及速度场

Figure 15. Vorticity and velocity field at time t = 5

图15. t = 5时刻涡度及速度场

Figure 16. Convergence orders of heat Navier-Stokes equation

图16. Navier-Stokes方程收敛阶

速度场和涡度场都随时间指数衰减，衰减率由运动学粘度控制。这反映了粘性耗散的效果。速度场和涡度场在空间上是周期性的。在图14和图15中，方块颜色代表涡度大小，箭头代表速度场方向。如图16可以看出，泰勒格林涡方程的收敛在δ = 0.5的情况下，正如预期的那样，误差近似地收敛于k阶。虽然只考虑了无界域中的Navier-Stokes方程，但这些结果清楚地表明，此方法作为一种可行的无网格计算流体动力学模拟方法，在复杂的几何形状中，这里选择的涡流函数公式需要仔细(且计算代价高昂)的边界处理，而对于三维流动则更加昂贵，因为必须为流函数向量场的每个分量求解泊松方程。

4. 结论

本文提出了一种不依赖于核函数的高阶一致性光滑粒子动力学方法。对领域粒子向目标粒子的泰勒展开建立线性方程组，通过求解所得的权重集合加权获得相应高阶一致性的微分算子。采用一维、二维函数与偏微分方程针对算法进行了2，3，4阶收敛一致性验证。通过这些算例的结果分析，得出以下结论：

通过算例结果，证明了不依赖于基函数形式，仅通过对于领域粒子权重的直接矩阵求解的方法仍然有效，并且结果表明无论在粒子正交抑或是非正交分布下，算法均能符合相应的收敛速率，证明了算法的稳定性与可行性，并且降低了算法的复杂度。同时，由于本文在算例的随机性上有更高的要求，与以往的研究相比，本文的算例粒子分布的随机性更大，这对算法的精度提出了更高的要求，而结果也证明了算法符合相应的精度要求。该高精度无网格算法在处理复杂几何域数值求解问题中具有潜在优势。

基金项目

国家自然科学基金(12002212)。