SPH 粒子法求解流体 — 从核函数到人工粘性

h 比几个粒子间距还短时，核求和会变稀；关掉人工粘性，冲击处粒子会叠在一起。

1992 年,Joe Monaghan 在一篇天体物理综述中抛出了一个挑衅性的设想:别再围绕恒星铺网格了。把气体撒成数千个相互重叠的粒子,让它们通过光滑的压力场相互作用。三十年后,Smoothed Particle Hydrodynamics(SPH)在天体物理、自由表面流和冲击模拟中占据了一席之地。撑住它的是四个方程,旁边是经常踩进去的坑。

在高斯、三次样条、poly6、spiky 之间切换核函数，钟形和中心斜率会一起变。

kernel

h = 1.00

Cubic spline has compact support on r < 2h. Poly6 / spiky vanish at r = h. Gaussian is non-compact (truncated at 2.5h here).

poly6 在中心附近斜率熄灭。spiky 在中心仍然尖。压力项用哪一个，就由这个差别决定。

没有网格的流体 — SPH 的起点#

欧拉法把坐标系钉死,看流动从上面流过。SPH 反过来,把坐标系绑在流场上,让每个粒子带着自己的速度和密度一起走。拉格朗日形式的连续方程和动量方程是:

\frac{D\rho}{Dt} = -\rho \nabla \cdot \mathbf{u}, \qquad \frac{D\mathbf{u}}{Dt} = -\frac{\nabla P}{\rho} + \mathbf{g}.

其中 $\rho$ 为密度, $\mathbf{u}$ 为速度, $P$ 为压力, $\mathbf{g}$ 为重力。沿粒子轨迹看,左边的时间导数就是普通的 $dt$ 。问题在右边的空间导数。粒子集合没有规则的邻居,该怎么计算 $\nabla P$ ?这就是 SPH 要回答的问题。

核函数 W(r, h) — 粒子构造的虚拟体积#

办法是卷积。在每个粒子周围放一个钟形函数 $W(r, h)$ ,把任意场 $A(\mathbf{x})$ 用邻居粒子的加权和近似:

\tilde{A}(\mathbf{x}) = \sum_{j} \frac{m_j}{\rho_j}\, A_j\, W(\mathbf{x} - \mathbf{x}_j, h).

$m_j$ 是粒子质量, $\rho_j$ 是密度, $h$ 是平滑长度。比值 $m_j/\rho_j$ 解释为粒子 $j$ 占据的体积。求导也照搬到核函数上 — $\nabla A \approx \sum_j (m_j/\rho_j) A_j \nabla W$ 。网格上的差分变成了核函数的导数。格子扭了以后还想保住梯度，另一套权重是移动最小二乘。

合格的核函数要满足两条:积分为 1, $h \to 0$ 时收敛到狄拉克函数。常用的形式都是紧支撑的 — 在 $r < h$ 或 $r < 2h$ 之外为零,每个粒子只与少数邻居打交道。高斯、三次样条、poly6、spiky 的支撑宽度和中心斜率各不相同。

三次样条核在 $r = 2h$ 处平滑趋零。poly6 适合密度估计,但它的梯度在 $r$ 接近 0 时变平,因此不适合压力计算。这正是图形 SPH 在压力项里改用 spiky 核的原因。高斯核没有紧支撑,代价是要把所有粒子对都算一遍。天体物理 SPH 几乎一律用三次样条核。

密度与压力 — 用求和定义的场

粒子 $i$ 的密度就是邻居核的和:

\rho_i = \sum_{j} m_j\, W(\mathbf{x}_i - \mathbf{x}_j, h).

有了密度,就用状态方程(EOS)算压力。弱压缩自由表面流通常用 Tait 状态方程:

P_i = K\left[\left(\frac{\rho_i}{\rho_0}\right)^{\gamma} - 1\right].

$\rho_0$ 是静止密度, $K$ 是刚度, $\gamma$ 通常取 7。这条公式让密度的微小变化产生很大的压力响应,以模拟不可压缩性。让 $K$ 取值使声速 $c_0 = \sqrt{\gamma K / \rho_0}$ 约为最大流速的十倍,马赫数就能保持很低。

用代码里的数字把这条式子踩一遍。 $N=80$ ，区间从 0.5 到 1.0。间距是 $0.5/79\approx0.00633$ 。粒子质量 $m=0.5\times1000/80=6.25$ 。 $m/\rho_0=6.25/1000=0.00625$ ，平衡密度下的粒子体积几乎等于间距。 $h=0.018$ 时 $2h=0.036$ ，一侧邻居有 $0.036/0.00633\approx5.7$ 个，五个到六个。 $h$ 降到 0.03，支撑里还有粒子。再短，求和就变稀。

把 $\gamma=7$ 、 $K=200$ 代进 Tait，密度 1% 会在压力里放大多少一目了然。 $\rho/\rho_0=1.01$ 时 $(1.01)^7\approx1.072$ 。 $P=K\times0.072=14.4$ 。1% 变成压力 14。用同一组 $K$ 、 $\gamma$ 、 $\rho_0$ 写声速， $c_0=\sqrt{7\times200/1000}=\sqrt{1.4}\approx1.18$ 。代码里的人工粘性不用这个数，用的是 $c_0=20$ 。粘性项比状态方程的声速大约强一个量级。 $\alpha=0.5$ ，代码里没有 $\beta$ 项。弱坍塌靠 $-\alpha c_0\mu/\bar\rho$ 还能撑住。冲击一强， $\mu^2$ 项缺了，粒子又会叠到一起。

守住动量的对称化技巧

直接写 $\nabla P / \rho$ 的 SPH 形式不会守恒动量。关键的恒等式是

\nabla\!\left(\frac{P}{\rho}\right) = \frac{\nabla P}{\rho} + \frac{P}{\rho^2}\nabla\rho.

代入后得到对粒子对对称的压力项:

\frac{D\mathbf{u}_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}.

$j$ 对 $i$ 的力与 $i$ 对 $j$ 的力恰好反号,牛顿第三定律对每一对粒子都成立,总动量与角动量只随浮点误差漂移。这种对称化撑起了 SPH 稳定性的一半。

给冲击波加人工粘性 — Monaghan-Gingold 配方#

仅凭上面这些方程,粒子在冲击波处会互相穿透。Monaghan 和 Gingold(1983)提出只对靠近的粒子对生效的人工粘性:

\Pi_{ij} = \begin{cases} \dfrac{-\alpha\, \bar{c}_{ij}\, \mu_{ij} + \beta\, \mu_{ij}^2}{\bar{\rho}_{ij}} & \mathbf{u}_{ij} \cdot \mathbf{r}_{ij} < 0 \\ 0 & \text{otherwise} \end{cases}

\mu_{ij} = \frac{h\, \mathbf{u}_{ij} \cdot \mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2 + 0.01 h^2}.

通常取 $\alpha \approx 0.5$ – $1.0$ , $\beta = 2\alpha$ 。粒子远离时不作用,因此在拉伸区域不会粘住粒子。副作用是它在剪切层里也会触发,这促使后来出现了 Balsara 开关等改进。

用 Python 模拟一维列坍塌#

来看 80 个垂直堆叠的粒子在重力下坍塌的一维模型。

import numpy as np
 
def cubic_spline_1d(r, h):
    """一维三次样条核 W(r, h)."""
    q = abs(r) / h
    c = 2.0 / (3.0 * h)
    if q < 1:
        return c * (1 - 1.5 * q * q + 0.75 * q ** 3)
    if q < 2:
        return c * 0.25 * (2 - q) ** 3
    return 0.0
 
def cubic_spline_dW(r, h):
    """W(r, h) 的一阶导数,符号随 r."""
    q = abs(r) / h
    c = 2.0 / (3.0 * h * h)
    s = 1.0 if r > 0 else (-1.0 if r < 0 else 0.0)
    if q < 1:
        return c * (-3 * q + 2.25 * q * q) * s
    if q < 2:
        return c * (-0.75 * (2 - q) ** 2) * s
    return 0.0
 
def sph_density_pressure(x, m, h, rho0, K, gamma=7.0):
    n = len(x)
    rho = np.zeros(n)
    for i in range(n):
        for j in range(n):
            rho[i] += m * cubic_spline_1d(x[i] - x[j], h)
    p = K * ((rho / rho0) ** gamma - 1.0)
    return rho, p
 
def sph_acceleration(x, v, rho, p, m, h, alpha, c0, g):
    n = len(x)
    a = np.full(n, g)
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            r = x[i] - x[j]
            if abs(r) >= 2 * h:
                continue
            grad = cubic_spline_dW(r, h)
            a[i] -= m * (p[i] / rho[i] ** 2 + p[j] / rho[j] ** 2) * grad
            v_ij = v[i] - v[j]
            if v_ij * r < 0:                       # 靠近的粒子对
                mu = h * v_ij * r / (r * r + 0.01 * h * h)
                rho_bar = 0.5 * (rho[i] + rho[j])
                Pi = -alpha * c0 * mu / rho_bar
                a[i] -= m * Pi * grad
    return a
 
# 一维列: 80 个粒子在 y in [0.5, 1.0] 内均匀分布
N = 80
x = np.linspace(0.5, 1.0, N)
v = np.zeros(N)
m = 0.5 * 1000.0 / N        # 列总质量 / 粒子数
h = 0.018
rho0, K = 1000.0, 200.0
alpha, c0 = 0.5, 20.0
g, dt = -9.81, 5e-5
 
for step in range(8000):
    rho, p = sph_density_pressure(x, m, h, rho0, K)
    a = sph_acceleration(x, v, rho, p, m, h, alpha, c0, g)
    v += a * dt
    x += v * dt
    below = x < 0.0                                 # 底壁反射边界
    x[below] = -x[below]
    v[below] = -0.5 * v[below]
 
print(f"final mean rho = {rho.mean():.1f}, top particle y = {x.max():.3f}")

跑完 8000 步后,平均密度在 $\rho_0$ 附近振荡,顶部粒子稳定在 $y \approx 0.55$ 处,形成一个更密的底层。这段小程序已经包含了 SPH 求解器的全部骨架 — 求和算密度,EOS 算压力,对称项算加速度,人工粘性维持稳定。

亲手拨动粒子

把平滑长度 $h$ 的滑条推到 0.07，每个粒子看到的邻居变多，颜色趋于均匀。降到 0.03，粒子之间拉开空隙，压力震荡跟着变大。

smoothing length h = 0.045stiffness K = 3.0viscosity mu = 3.5show kernel radius (first particle)

64 particles, poly6 density + spiky pressure kernel. Color encodes local density (blue = sparse, orange = compressed).

把刚度 $K$ 调高更接近不可压缩,但声速也会上升,时间步必须缩小。把粘性 $\mu$ 调到零会看到粒子互相穿插的条纹 — 直观地解释了人工粘性为什么必须存在。

SPH 不擅长的场景#

凡有亮面便有阴影。SPH 在强剪切流中数值粘性容易暴走,常常压过物理粘性。粒子聚团会破坏各向同性,使表面张力等微小效应难以捕捉。可变平滑长度 $h_i$ 带来自适应性,但需要在动量项里加上修正(Springel 和 Hernquist 2002 的 $f_i$ 项)。在自由表面附近密度计数骤减,压力可能为负,产生让表面颤动的 tensile instability — XSPH 和 $\delta$ -SPH 等方案就是为此设计的。

把 h 取得比粒子间距还短#

$h=0.018$ 大约是间距 0.00633 的三倍。支撑半径 $2h$ 的一侧邻居超过五个。这个宽度下密度求和还光滑。把 $h$ 收到一个间距附近， $W$ 几乎只看见自己， $\rho_i$ 就开始晃。 $\gamma=7$ 把这份晃动放大成压力。时间步 $5\times10^{-5}$ 走 8000 步，时间是 0.4 秒。这段里顶部粒子从 1.0 落到大约 0.55，是原文记下的结果。压力项把重力接住了。

拿掉对称项，粒子对的力不再抵消。跳过 $\nabla(P/\rho)$ ，总动量每一步都往一侧漏。人工粘性也有开关。只有粒子对互相靠近时 $\Pi_{ij}$ 才打开。代码只放了 $\alpha=0.5$ 的线性项。自由表面处密度掉下去，压力变成负的，粒子互相拉，表面就抖。那种成对靠拢写在固定核的拉伸不稳定里。

参考文献

Monaghan, J. J. (1992). Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics, 30, 543.
Monaghan, J. J., and Gingold, R. A. (1983). Shock simulation by the particle method SPH. Journal of Computational Physics, 52(2), 374.
Müller, M., Charypar, D., and Gross, M. (2003). Particle-based fluid simulation for interactive applications. SCA '03.
Springel, V. (2005). The cosmological simulation code GADGET-2. MNRAS, 364, 1105.