SPH 核近似、守恒离散与时间步

光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）把连续介质离散为一组携带质量、速度与密度的粒子，用核函数加权邻域求和近似任意场量及其导数。它天然适合自由表面、大变形、多相与固液耦合问题：射流破碎、溃坝、波浪爬高、齿轮搅油、血液流动、滑坡涌浪等用网格法处理时界面拓扑变化代价高昂的场景。本文从核近似出发，推导成对对称的动量离散，讲清一致性、离散守恒与显式时间步约束的工程含义，并给出一个可复现的一维静水柱算例。

1. 结论与适用场景

SPH 的结论可以先行归纳为三句话：核近似只在邻域采样足够时才接近连续积分；守恒性来自成对力的反对称性而非核函数本身；时间步由声速、加速度与黏性三个尺度共同约束。 三者缺一，结果往往"看起来像流体，物理量却不对"。

适用判据可归纳为三条：界面是否发生拓扑变化；是否存在大剪切与拉伸；是否接受以显式小时间步换取实现简洁。若问题以定常内流为主、几何规则、追求高精度壁面边界层，有限体积法通常更高效。SPH 的代价是计算量随粒子数近似线性增长，且自由面与壁面附近的一致性需要专门处理。

需要强调的是，弱可压缩 SPH（WCSPH）用状态方程把压力与密度联系起来，回避了压力泊松方程，代价是数值声速必须足够大，时间步因此变小。理解这一权衡，是判断该用 SPH 还是网格法的核心。

按压力的求解方式，SPH 可分为三类：弱可压缩（显式积分，实现最简单、并行扩展最好，但受声速时间步限制）、不可压缩 SPH（ISPH，用投影法每步求解压力泊松方程，时间步更大、压力噪声更低）、密度不变式 SPH（DFSPH，用迭代满足密度约束，常见于图形学实时流体）。另有 Reno 等基于 Renormalization 的弱可压缩变体。选择哪一类，取决于对不可压性、时间步与实现复杂度的取舍；本文其余部分以 WCSPH 为主线。

2. 核近似与守恒离散

核近似的思想是用满足归一化、紧支撑、非负与对称性的核 $W$ 去平滑 δ 函数。连续核积分满足一致性条件，但离散到无序粒子上时，只有当邻域求和恰当地近似体积元 $V_j=m_j/\rho_j$，才能保持零阶一致性。因此核函数的连续性质与粒子求和的实际精度是两回事。

守恒性必须在离散层面刻意构造。把压力加速度写成成对形式，使粒子 $i$ 与 $j$ 受到的内力大小相等、方向相反，这依赖核梯度的反对称关系

$$\nabla_i W_{ij} = -\nabla_j W_{ij}$$

只要核仅依赖距离，该关系自动成立，离散动量便严格守恒（浮点舍入除外）。能量守恒则要求压力做功与密度更新在时间上一致耦合。反过来，若使用核梯度的不同变体、非对称平滑长度，或人为截断邻居表，都会破坏对称性并引入虚假动量——这正是后面"常见坑"反复出现的根源。

邻域规模由支持半径 $r_c=\kappa h$ 决定。二维与三维的邻居数差异明显，相同比值 $\eta=h/\Delta p$ 下三维邻居数增长更快，参数不能跨维度照抄。核截断还使边界与自由面附近的求和缺少外侧贡献，形成"邻居不足"，需靠边界粒子层或修正规化补偿。

核函数须满足归一化 $\int W\,dV=1$、紧支撑（$r>r_c$ 时 $W=0$）、对称性，以及 $h\to0$ 时趋近 δ 函数。紧支撑是把邻居搜索限制在 $r_c$ 内的前提，也是计算量可控的关键。核的光滑阶数决定梯度近似的连续性，而稳定性还与邻居数相关，因此核选取不是孤立问题，需要与边界处理和密度计算一起考虑，详见姊妹篇《SPH 核函数、平滑长度与粒子间距》。

3. 核心公式与推导

标量场的核近似与粒子求和为

$$A(\mathbf{x}_i)\approx\int_\Omega A(\mathbf{x}')W(\mathbf{x}_i-\mathbf{x}',h)\,d\mathbf{x}'\approx\sum_j V_j A_j W_{ij}$$

其中 $W_{ij}=W(|\mathbf{x}_i-\mathbf{x}_j|,h)$，$V_j=m_j/\rho_j$。对梯度算子，把梯度作用到核上并分部积分，得到差值形式

$$\nabla_i A_i \approx \sum_j V_j (A_j-A_i)\nabla_i W_{ij}$$

差值形式使常数场的梯度严格为零，因此比直接对 $A_j$ 求核梯度更稳健，是密度、速度散度等算子近似的首选。若直接近似 $\nabla A_i\approx\sum_j V_j A_j\nabla_i W_{ij}$，常数场会因核梯度和不为零而产生虚假梯度，这正是零阶一致性在离散层面失守的直接体现。

密度有两条路径。其一是求和式

$$\rho_i=\sum_j m_j W_{ij}$$

其二是连续性方程形式。求和式简单、并行友好，但在自由表面因核截断低估密度，产生所谓表面密度亏损，表现为假的表面张力与粒子飞散，工程上常用密度扩散（δ-SPH）或 Shepard 修正缓解。

动量方程取成对对称压力项：

$$\frac{d\mathbf{v}_i}{dt}=-\sum_j m_j\left(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2}\right)\nabla_i W_{ij}+\mathbf{g}+\mathbf{a}_{\nu,i}$$

层流黏性项常用对称形式

$$\mathbf{a}_{\nu,i}=\sum_j m_j\frac{2\mu_i\mu_j}{\mu_i+\mu_j}\frac{\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}}{|\mathbf{x}_{ij}|^2+\epsilon h^2}\nabla_i W_{ij}$$

其中 $\mathbf{v}_{ij}=\mathbf{v}_i-\mathbf{v}_j$，$\mathbf{x}_{ij}=\mathbf{x}_i-\mathbf{x}_j$，$\epsilon\approx0.01$ 用于防止粒子接近时除零。

压力由 Tait 状态方程给出：

$$p_i=\frac{\rho_0 c_s^2}{\gamma}\left[\left(\frac{\rho_i}{\rho_0}\right)^{\gamma}-1\right]$$

常用 $\gamma=7$。取 $c_s\geq 10\,v_{\max}$ 可把密度波动压到 1% 以内，声速越大则不可压性越好、时间步越小。显式时间步同时受三个尺度约束：

$$\Delta t \leq \min\left(C_c\frac{h}{c_s+v_{\max}},\ C_a\sqrt{\frac{h}{a_{\max}}},\ C_\nu\frac{h^2}{\nu}\right)$$

典型取 $C_c\approx0.25$、$C_a\approx0.25$、$C_\nu\approx0.125$，声学项通常主导。

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

平滑长度与粒子间距：常用 $h=\eta\Delta p$，$\eta$ 取 1.2~1.5（二维约 1.3，三维约 1.2）。$h$ 过大则过度平滑、界面模糊；过小则邻居不足、压力噪声大。应保证内部粒子邻居数在二维约 20~30、三维约 50~80。

CFL 与声速：先由时间步约束反推允许的最大 $c_s$，再确认密度波动是否落在目标范围（通常 1% 以内）。降低 $\Delta t$ 或提高 $c_s$ 都会增大计算量，二者是精度与成本的主要旋钮。

人工黏性：Monaghan 型人工黏性用于抑制冲击与粒子接近时的非物理飞散：

$$\Pi_{ij}=-\alpha_{\Pi}\frac{c_s h\,\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}}{\rho_{ij}(|\mathbf{x}_{ij}|^2+\epsilon h^2)}$$

仅当粒子相互接近（$\mathbf{v}_{ij}\cdot\mathbf{x}_{ij}<0$）时激活，$\alpha_{\Pi}$ 常取 0.01~0.1。它会耗散物理动能，必须单独记录系数并做敏感性分析，避免把它当成材料黏度。

边界处理：常用固定壁粒子（Wakui/Adami 型）或镜像粒子产生壁面压力与斥力。壁粒子层数要足够覆盖支持半径，否则近壁粒子一侧缺邻居，出现"贴壁飞散"或穿壁。固壁无滑移则需额外施加切向黏性约束。

稳定性辅助：密度扩散抑制压力噪声，粒子移位修复拉伸区的粒子分布。每一项都改变有效物理，应独立开关验证。

参数扫描顺序：先固定 $h/\Delta p$ 与核函数，扫描 $\alpha_{\Pi}$ 直到冲击稳定；再调 $c_s$ 把密度波动压到目标范围；最后按时间步分解项确认主导约束，避免盲目全局缩小 $\Delta t$。

积分与近邻：显式积分常用预测—校正或辛格式（如 kick-drift-kick）。每步从 Verlet 表读取邻居，位移超过阈值时重建，且所有求和应保持确定性顺序，或定义统计容差。

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

一维静水柱：$N=200$ 粒子等间距排在 $[0,L]$，密度 $\rho_0=1000$，采用 Tait 状态方程，底壁由固定壁粒子表示，重力 $g=9.81$。运行到 $t=0.5$ s，检查静水压分布是否线性、自由面是否稳定。

输入: N, dx, rho0, cs, gamma, g
初始化 x[i]=i*dx, v=0, rho=rho0, m=rho0*dx
while t < 0.5:
    # 1. 近邻
    if 位移 > skin/2: 重建 Verlet 表(半径 rc+skin)
    # 2. 密度与压力
    rho[i] = sum_j m_j * W(r_ij, h)
    p[i]   = rho0*cs^2/gamma * ((rho[i]/rho0)^gamma - 1)
    # 3. 加速度（成对求和）
    for pair (i,j) with r_ij < rc:
        pterm = (p_i/rho_i^2 + p_j/rho_j^2) * gradW
        Pi    = 人工黏性（仅接近时）
        a[i] += -m_j*(pterm - Pi)*gradW
    加上重力与壁面斥力
    # 4. 显式积分 (kick-drift-kick)
    dt = min(dt_cfl, 0.25*h/(cs+vmax), 0.25*sqrt(h/amax))
    v += a*dt/2;  x += v*dt;  v += a*dt/2
    t += dt

验证步骤：其一，常数密度场下 $\rho_i/\rho_0$ 应接近 1，自由面附近允许略低；其二，静水柱底部压力应等于 $\rho_0 g H$，误差随 $h$ 减小而收敛；其三，无外力时净动量应保持，成对受力逐对抵消；其四，粒子数加倍后压力剖面按一阶收敛。

6. 常见坑与排查

• 表面密度亏损：自由面附近求和式密度偏低，产生假的表面张力与飞散。排查时绘制 $S_i=\sum_j V_j W_{ij}$ 在表面处的取值，再用密度扩散或修正公式。
• 把人工黏性当物理黏度：若改变 $\alpha_{\Pi}$ 显著改变结果，说明黏性主导，需要分离材料黏性与数值黏性。
• 时间步只满足声学 CFL：漏掉加速度或黏性约束，会在强加速度区突然发散。
• 邻居不足：拉伸区粒子稀疏、邻居数骤降会导致核和错误，应监控最小与平均邻居数。
• 非对称平滑长度：各粒子用不同 $h$ 时若不对称处理，会破坏动量守恒。
• 浮点归约顺序：并行下逐位不可复现，应比较物理统计量而非逐位相等。
• 声速过小：密度波动超过 1% 甚至出现负压，表现为体积不可压性丧失。
• 时间步自适应只取全局最大量：应在计算前扫描局部最大速度与加速度，按粒子而非全局统一取步长，否则局部爆点会污染全局。
• 初始分布过疏或过密：初始构型未做松弛（如按格点排布后让壁面压力平衡）会在起步阶段喷射，应先跑一段稳定化再开始统计。
• 把边界当自由面：边界粒子未参与邻居统计会让近壁密度偏低，需统一把边界纳入求和，或用修正规化补偿缺失体积。

7. 检查清单与参考

• [ ] 常数场：$\rho_i/\rho_0\approx1$，梯度近似为零。
• [ ] 线性场：一阶导数误差随 $h$ 收敛。
• [ ] 成对受力：大小相等、方向相反，净内力为零。
• [ ] 分区域统计邻居数（内部/表面/壁面）与零邻居粒子数。
• [ ] 监测密度范围、最小粒子间距、最大加速度与时间步分解项。
• [ ] 加密时同步检查 $h/\Delta p$、$\Delta t$ 与边界粒子层数。
• [ ] 守恒量（动量、能量）漂移在可接受范围内。

参考资料：

1. Monaghan J.J., Smoothed Particle Hydrodynamics, Annual Review of Astronomy and Astrophysics, 1992.
2. Monaghan J.J., Smoothed Particle Hydrodynamics and Its Diverse Applications, Annual Review of Fluid Mechanics, 2012.
3. Liu G.R., Liu M.B., Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, 2003.
4. Koschier D. et al., SPH Techniques for the Physics Based Simulation of Fluids and Solids, Eurographics Tutorial, 2019.