SPH基础(一): 核函数与导数近似

欢迎来到“SPH 系列教程”系列！

光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）是一种无网格的拉格朗日粒子法。与传统的基于网格的方法不同，SPH 通过一系列离散的粒子来代表连续的流体，极大地简化了对大变形、自由表面等问题的处理。

本系列的第一篇文章，我们将从 SPH 的两个最核心的概念——核函数近似与导数近似——开始。

什么是核函数近似？

在 SPH 中，任何一个连续场 $A$ 在空间任意一点 $\mathbf{r}$ 的值，都可以通过一个对邻近粒子的加权求和来近似：

$$ A(\mathbf{r}_i) \approx \sum_{j} A(\mathbf{r}_j) W(\mathbf{r}_i - \mathbf{r}_j, h) V_j $$

这里的：

$A(\mathbf{r}_i)$ 是我们想要求的粒子 $i$ 的场量值。
求和遍历粒子 $i$ 的所有邻居粒子 $j$。
$W$ 是核函数，一个根据距离分配权重的函数。
$h$ 是光滑长度，定义了核函数的影响范围。
$V_j$ 是粒子 $j$ 的体积，通常等于其质量除以密度 ($V_j = m_j / \rho_j$)。

这个公式的核心思想是：一个点的属性，可以由其周围点的属性加权平均得到。

如何近似一个场的导数？

SPH 的真正威力在于它也能方便地近似一个场的导数（如梯度、散度），这是构建物理控制方程（如流体力学的纳维-斯托克斯方程）的关键。场 $A$ 在粒子 $i$ 处的梯度 $\nabla A(\mathbf{r}_i)$ 可以通过以下对称形式来近似：

$$ \nabla A(\mathbf{r}_i) \approx \sum_{j} [A(\mathbf{r}_j) - A(\mathbf{r}_i)] \nabla_i W(\mathbf{r}_i - \mathbf{r}_j, h) V_j $$

其中 $\nabla_i W$ 是核函数对粒子 $i$ 坐标的梯度。这个形式具有良好的数值稳定性，并且保证了一个常数场的梯度为零。

与有限元法（FEM）的类比: 如果读者对有限元法有了解，会发现这两种方法在哲学上是相通的。无论是SPH还是FEM，它们都巧妙地将对未知场函数的微分操作，转移到了已知的、解析的基函数（SPH中的核函数，FEM中的形函数）上。这样做最大的好处是避免了对离散数据点进行直接的数值差分，因为后一种方法对粒子/节点的无序性和噪声非常敏感，容易导致数值不稳定和精度损失。

常用核函数及其梯度

一个好的核函数需要满足归一性、紧支撑性等性质。下面介绍几种在 SPH 中常用的核函数。

注意：以下公式中的归一化常数 $\alpha$ 均是三维空间下的值。在不同维度下，这些常数会发生变化。

1. 三次样条核 (Cubic Spline Kernel)

这是SPH中最经典和广泛使用的核函数之一，因其良好的稳定性和近似二阶高斯函数的特性而备受青睐。

其数学表达式为：

$$ W(R,h) = \alpha_d \times \begin{cases} \frac{2}{3} - R^2 + \frac{1}{2}R^3, & 0 \le R < 1; \\ \frac{1}{6}(2-R)^3, & 1 \le R < 2; \\ 0, & R \ge 2. \end{cases} $$

其中，归一化常数 $\alpha_d$ 在不同维度下分别为：

一维: $\alpha_d = \frac{1}{h}$
二维: $\alpha_d = \frac{15}{7\pi h^2}$
三维: $\alpha_d = \frac{3}{2\pi h^3}$

2. 二次光滑核 (Quadratic Kernel for Impact Problems)

根据Johnson等人 (1996b) 的研究，在模拟高速冲击问题时，可以采用以下的二次光滑函数。

其数学表达式为：

$$ W(R,h) = \alpha_d \left(\frac{3}{16}R^2 - \frac{3}{4}R + \frac{3}{4}\right), \quad 0 \le R \le 2. $$

其中，归一化常数 $\alpha_d$ 在不同维度下分别为：

一维: $\alpha_d = \frac{1}{h}$
二维: $\alpha_d = \frac{2}{\pi h^2}$
三维: $\alpha_d = \frac{5}{4\pi h^3}$

Matlab 实验：一维函数与导数近似

让我们通过一个更全面的 Matlab 实验来感受这些概念。下面的代码将：

在一个一维域上近似函数 $y = \sin(x)$ 及其导数 $y’ = \cos(x)$。
提供多种核函数（cubic, poly6, spiky）供选择。
可视化整个域上的近似结果和误差。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
% SPH_APPROXIMATION_1D_ADVANCED.m
%
% 描述:
% 这是一个更高级的1D SPH实验，用于近似一个函数及其一阶导数。
% 1. 使用多种核函数（Cubic Spline, Poly6, Spiky）。
% 2. 近似整个函数域，而不仅仅是一个点。
% 3. 可视化函数近似和导数近似的误差。

clear; clc; close all;

%% 1. 参数设置
L = 2 * pi;             % 定义域长度 [0, 2*pi]
N = 100;                % 粒子数量
dx = L / N;             % 粒子间距 (每个粒子的“体积”)

% --- 可调参数 ---
h_factor = 2.0;         % 光滑长度因子 h = h_factor * dx
kernel_choice = 'cubic'; % 可选: 'cubic', 'poly6', 'spiky'
% ... (代码其余部分和上面提供的一样)