无网格法 · SPH 理论与实现

SPH 核函数、平滑长度与粒子间距

比较三次样条、五次样条与 Wendland 族核函数，说明平滑长度比、邻居数、离散一致性与配对不稳定之间的关系，并给出一个可复现的核函数选择算例。

源码知识库生成器更新于 2026/9/22

SPH核函数平滑长度一致性配对不稳定进阶

SPH 核函数、平滑长度与粒子间距

核函数 $W(r,h)$ 决定一个粒子如何"读取"邻域信息，它与平滑长度、粒子间距一起构成 SPH 的分辨率与稳定性支柱。本文比较常用紧支撑核，讲清平滑长度比、邻居数、离散一致性与配对不稳定之间的关系，并给出一个可复现的核函数选择算例。

1. 结论与适用场景

结论先行：紧支撑核的连续性质不保证离散一致性；配对不稳定与核的傅里叶变换是否非负相关； $h/\Delta p$ 与维度共同决定邻居数，参数不可跨维度照抄。

常见核族分成四类。三次样条核计算便宜、历史悠久，但在邻居数偏高时可能出现配对不稳定；五次样条核更平滑、傅里叶变换非负，代价是邻居数更多；Wendland 族（C2/C4/C6）支持半径小、在大邻域下稳定，是当前主流；高斯核没有紧支撑，必须截断，截断处导数不连续会让一致性与守恒性变差，通常不推荐用于有边界的工程模拟。

适用场景可以这样判断：当模拟出现压力振铃、粒子团聚或拉伸区的"带状"伪结构时，应优先改用傅里叶变换非负的核，而不是盲目减小 $\Delta p$ ；当界面清晰度优先、允许适度噪声时，可选支撑更窄的核并配合人工黏性。

需要提醒的是，核函数不是一个可以孤立调优的旋钮。它同时影响求和成本、压力噪声、界面厚度与稳定性，而且与边界处理、密度算法和人工黏性强烈耦合。合理的工作流是：先固定一套公认鲁棒的核（如 Wendland C2）与平滑长度比，把整体流程跑通并验证守恒与收敛，再回过头做核的敏感性对比。跨维度迁移参数时，必须重新推导归一化常数并重新统计邻居数。

2. 核近似与守恒离散

核近似把场量写成邻域加权和，核必须归一化、对称且紧支撑。连续层面的归一化与对称性容易满足，但离散到无序粒子后，求和精度取决于局部粒子分布的规则程度。衡量指标是离散归一化量 $S_i$ ，它偏离 1 直接对应一致性误差。

守恒性对核的要求是"距离对称"：成对求和时，只要 $W$ 仅依赖 $|\mathbf{x}_i-\mathbf{x}_j|$ ，就有 $\nabla_i W_{ij}=-\nabla_j W_{ij}$ ，压力项与黏性项的成对力天然反对称，离散动量守恒。若引入可变平滑长度 $h_i\neq h_j$ ，必须对称化（取算术或几何平均），否则反对称性被破坏，出现净动量漂移。

核与邻居表是一套系统：支持半径 $r_c=\kappa h$ 同时决定邻居规模、求和成本与配对稳定性。因此核函数选择必须与边界处理、密度计算和稳定化项共同验证，而不是单独调参。

边界与自由面是另一条主线。那里核被几何地截断，无论用哪种核， $S_i$ 都会低于 1。如果边界处理不当，常见表现是近壁粒子被"吸"向壁面或向外飞散；前者往往因为壁面压力估计偏低，后者常因壁粒子层数不足或修正缺失。因此，评估核函数的正确姿势不是只算内部规则粒子，而是把内部、壁面与自由面分别统计。

3. 核心公式与推导

归一化条件要求核的积分为 1：

\int_\Omega W(\mathbf{r},h)\,d\mathbf{r}=1

三维下即 $\int_0^{r_c} W(r,h)\,4\pi r^2\,dr=1$ ，由此定出各核的归一化常数 $\sigma_D/h^D$ ，其中 $D$ 为维数。用无量纲 $q=r/h$ 表示时，三次样条核在 $0\leq q<1$ 时 $W\propto 1-\frac{3}{2}q^2+\frac{3}{4}q^3$ ，在 $1\leq q<2$ 时 $W\propto \frac{1}{4}(2-q)^3$ ， $q\geq2$ 时为零。Wendland C2 核在三维下写作

W(q)=\frac{21}{16\pi h^3}\left(1-\frac{q}{2}\right)^4\left(2q+1\right),\qquad 0\leq q\leq 2

其傅里叶变换在整个支持域内非负，因此对高邻居数稳定。作为对比，五次样条核在 $0\leq q\leq3$ 上定义并具有非负傅里叶谱，代价是 $\kappa=3$ ，邻居数增至三次样条的约 $(3/2)^D$ 倍，在三维中接近三倍的计算量。离散归一化量定义为

S_i=\sum_j V_j W_{ij},\qquad V_j=\frac{m_j}{\rho_j}

理想情况下 $S_i\approx1$ 。自由表面与壁面因核截断使 $S_i<1$ ，这是几何原因而非核函数缺陷。相应地，梯度一致性量应接近零：

\boldsymbol{\beta}_i=\sum_j V_j \nabla_i W_{ij}\approx\mathbf{0}

若 $\boldsymbol{\beta}_i$ 显著偏离，可用 Shepard 修正重规化核：

W_{ij}^{\text{corr}}=\frac{W_{ij}}{\sum_j V_j W_{ij}}

内部粒子邻居数的量级估计为

N_{2D}\approx \pi(\eta\kappa)^2,\qquad N_{3D}\approx \frac{4}{3}\pi(\eta\kappa)^3

其中 $\eta=h/\Delta p$ 、 $\kappa=r_c/h$ 。这解释了为何相同 $\eta$ 下三维邻居数远大于二维。配对不稳定的判据来自线性稳定性分析：离散扰动增长率与核的傅里叶变换 $\hat{W}(k)$ 相关，若 $\hat{W}$ 在某波数下变负，规则排布便失稳并配对成簇。三次样条在邻居数偏高时会越过该阈值；Wendland C2/C4/C6 与五次样条在整个支持域保持非负，因而更稳。

实践中可用一条经验法则：先按目标邻居数确定 $\eta$ 与 $\kappa$ ，再用扰动测试确认没有配对，最后才审阅压力噪声。若必须使用三次样条，应把邻居数控制在较低范围（二维约 20 以内、三维约 40 以内），并留意拉伸区是否出现成对结构；一旦出现，换核比继续降 $\Delta p$ 更有效。

4. 数值实现与参数（h/粒子间距、CFL、人工黏性）

平滑长度比 $\eta$ ：二维常取 $\eta\approx1.3$ 、三维 $\eta\approx1.2$ ，对应内部邻居数二维约 20~30、三维约 50~80。 $\eta$ 增大则更平滑、噪声更低，但界面更厚、局部峰值被抹平； $\eta$ 减小则分辨率提高，但压力噪声与配对风险上升。

支持半径倍数 $\kappa$ ：三次样条与 Wendland C2 通常取 $\kappa=2$ 。 $\kappa$ 与 $\eta$ 是独立参数： $\eta$ 定分辨率， $\kappa$ 定每颗粒子扫过的邻居层数。增大 $\kappa$ 提高稳定性与平滑度，却也直接抬高每步的邻居数，成本近似随 $\kappa^D$ 增长。

CFL：时间步与 $h$ 成正比、与最大声速成反比。换核后 $h$ 不变则时间步基本不变，但换核改变了有效邻域，需重新核查密度波动与稳定性。

人工黏性：窄支撑核更易产生压力噪声，必要时提高人工黏性系数，但应独立记录，避免把数值耗散算成物理黏度。

导数形式：优先使用差值形式的梯度近似以保住常数场一致性；核导数应解析给出，避免数值差分引入额外误差。

维度一致性：二维与三维的归一化常数不同，二维 SPH 通常按单位厚度处理，参数不能直接迁移。

邻居数的实证校验：解析估计只是量级参考，实际邻居数还受粒子分布、边界截断与支持半径与单元尺寸不匹配的影响。部署前应直接统计每个粒子的实际邻居数，取最小、平均与最大值，并单独记录零邻居粒子；内部粒子邻居数应落在目标区间，自由面与角点允许偏低但不能为零。

可变平滑长度的对称化：若某区域需要更小的 $h$ ，务必用对称化的核参数（如 $h_{ij}=(h_i+h_j)/2$ ）并乘上适当的体积因子，否则成对力不再反对称，长时间积分会出现动量漂移。这类漂移往往在短时间看不出，却在数千步后显现为整体平移。

5. 可复现示例（伪代码/小算例）

在二维规则格点上计算离散归一化 $S_i$ 与梯度误差，比较三次样条与 Wendland C2，并扫描 $\eta$ 。

输入: 核类型 kernel, eta, kappa, 格点间距 dp, 域 Lx*Ly
h = eta*dp;  rc = kappa*h;  V_j = dp*dp
for 核 in [cubic, wendlandC2]:
    S_min, S_max, S_mean, gmax = +inf, -inf, 0, 0
    for 内部粒子 i:
        s = 0; gx = 0; gy = 0
        for 邻居 j (r_ij < rc):
            W  = kernel(r_ij, h)            # 含维度归一化常数
            dW = kernel_grad(r_ij, h)       # 解析核导数
            s  += V_j * W
            gx += V_j * dW * dx_ij / r_ij   # 梯度一致性
            gy += V_j * dW * dy_ij / r_ij
        统计 S = s 与 |(gx, gy)|
    print 核, eta, 邻居数统计, S 范围, 最大梯度误差

验证步骤：其一，规则格点上 $S_i$ 应非常接近 1，两种核都不例外；其二，梯度误差之和应接近零；其三，扫描 $\eta=1.1,1.2,1.3$ 观察邻居数与噪声的变化；其四，对粒子施加随机扰动后观察是否出现配对团聚，比较两种核的稳健性——这正是配对不稳定在实践中的直接表现。

6. 常见坑与排查

跨维度照抄 $\eta$ ：二维与三维邻居数差异大，照抄会得到过大或过小的邻域。
忽略归一化常数：二维与三维的 $\sigma_D$ 不同，混用会让 $S_i$ 系统性偏离 1。
盲目增大 $\eta$ 降噪：会模糊界面、增大耗散，应优先换更稳的核或调人工黏性。
可变 $h$ 未对称化： $h_i\neq h_j$ 未做对称平均会破坏动量守恒。
使用截断高斯核：截断处或其导数不连续，一致性明显变差。
把表面 $S_i<1$ 当 bug：这是几何截断的结果，应通过边界粒子或规化处理，而非强行限制密度下限。
只看单个粒子：应统计整场的 $S$ 范围与分布，否则偶然性会掩盖问题。
支持半径与单元尺寸不匹配：近邻单元边长应不小于 $r_c$ ，否则一次查询覆盖的单元数与邻居数波动很大，影响性能与正确性。
换核后未重跑守恒检查：不同核的有效邻域与压力噪声不同，动量、能量与密度波动需重新验证。
把 $S_i$ 均值当唯一指标：均值接近 1 但最小值很低时，仍可能在局部产生伪力，应同时关注分位数。
忽略核的函数形式与差分算子匹配：梯度与拉普拉斯算子近似中常出现核的二次导数，需确认所用核在该阶导数连续（尤其样条分段处），否则高阶项会引入噪声。

7. 检查清单与参考

[ ] 归一化常数与维度匹配， $\sigma_D$ 取值正确。
[ ] 规则粒子上 $S_i\approx1$ ；边界与表面处记录偏离量。
[ ] 梯度一致性量 $\boldsymbol{\beta}_i\approx\mathbf{0}$ 。
[ ] 邻居数落在二维 20~30、三维 50~80 的目标区间。
[ ] 扫描 2~3 组 $\eta$ ，比较邻居数、压力噪声与界面厚度。
[ ] 扰动粒子、拉伸区与自由面附近无配对团聚。
[ ] 记录核名称、 $\kappa$ 、 $\eta$ 、导数形式、截断策略与对称化方式。

参考资料：

Dehnen W., Aly H., Improving convergence in smoothed particle hydrodynamics simulations without pairing instability, Monthly Notices of the Royal Astronomical Society, 2012.
Monaghan J.J., Smoothed Particle Hydrodynamics, Annual Review of Astronomy and Astrophysics, 1992.
Liu M.B., Liu G.R., Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments, Archives of Computational Methods in Engineering, 2010.
Wendland H., Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in Computational Mathematics, 1995.