SPH平滑粒子流体模拟原理与实现

可视化演示

动态演示SPH平滑粒子流体模拟原理与的核心机制，通过总览/传输/转换/输出四种模式展示关键概念与交互流程。

SPH概述：不是网格，而是一锅食材

假设你是大厨，想精确掌握锅里每一滴水的去向。传统网格流体模拟的做法是给整口锅盖上一层密密麻麻的渔网，然后记录每个网眼里的平均速度和水位——这就是“欧拉视角”的网格法。现在换个思路：你不再盯着网格，而是盯住汤里的每一粒食材。胡萝卜丁、葱花、肉末各自以“微小一团流体代言人”的身份在锅里游走，它们带着自己的质量、速度、密度去碰撞、聚散，最终你观察到的仍然是整锅汤，但描述它的语言变成了“粒子”——这便是 SPH 所采用的“拉格朗日视角”。

SPH 全称 Smoothed Particle Hydrodynamics，中文叫“光滑粒子流体动力学”，1977 年由 Gingold 与 Monaghan 在研究天体物理（如恒星形成）时提出，后来被图形学工程师引入流体模拟。它不直接求解抽象的偏微分方程，而是用离散粒子的位置与属性重建速度场、密度场。每个粒子的行为遵循两条规则：一是“受周围粒子影响”（通过核函数加权），二是“贡献自己的质量与压力”影响周围。

为什么可以用一堆离散点代替连续流体？因为 SPH 有一个“插值魔法”：任意一点的物理量 f，都可以由附近粒子携带的值加权拼出来。这个拼贴公式就是 SPH 的基石：

$$\displaystyle f(\mathbf{r}) = \sum_{j} \frac{m_j}{\rho_j} f_j \, W(\mathbf{r} - \mathbf{r}_j, h)$$

1. 目标物理量 $f(\mathbf{r})$：这是要重建的连续场在位置 r 处的值，密度、压力、速度都适用。

$$\displaystyle f(\mathbf{r})$$

2. 粒子求和 $\sum_j$：遍历 r 附近光滑长度 h 内的所有粒子，逐个累加它们的贡献。

$$\displaystyle \sum_{j}$$

3. 体积权重 $\dfrac{m_j}{\rho_j}$：粒子 j 的质量除以密度等于它占据的“小体积”，体现它在空间里占了多少份额。

$$\displaystyle \frac{m_j}{\rho_j}$$

4. 粒子属性 $f_j$：粒子 j 自身携带的那个物理量的数值。

$$\displaystyle f_j$$

5. 核函数 $W(\mathbf{r}-\mathbf{r}_j, h)$：衡量粒子 j 对位置 r 的影响强度，距离越近权重越大，超过 h 就归零。

$$\displaystyle W(\mathbf{r} - \mathbf{r}_j, h)$$

回到汤里——这个公式就像“按勺子大小决定每人盛多少汤”的分配规则。每粒食材都在向周围宣告自己的存在，而任意一点能“尝到”的味道，正是附近食材贡献的加权平均。把这一思想铺开，就是 SPH 的全部起点。

核函数：香料影响一圈圈变淡

往汤里丢一粒花椒，它的麻味并不会整锅均匀扩散，而是离得越近越冲、越远越淡，超过某一圈距离就彻底尝不到了。这粒花椒的“影响曲线”就是核函数 W 的直观形象。光滑长度 h 就是那圈“尝不到”的边界——它决定了一粒食材能影响的半径，也决定了汤里的“社交圈”大小。

在 SPH 里，核函数是粒子之间传递信息的“规则”。它必须满足三个基本要求：归一化（整条影响曲线围出的面积等于 1，保证总量不凭空增加）、紧致支持（超过 h 就取 0，让计算只发生在局部）、非负且平滑（避免出现负的“反味道”和数值突变）。不同的核函数就像不同香料的扩散方式：有的香气锐利，有的绵长。

最常用的核函数之一是 Poly6，它由 Müller 等人于 1996 年引入图形学，形式优美、计算简单：

$$\displaystyle W_{poly6}(r, h) = \frac{315}{64\pi h^{9}} \left( h^{2} - r^{2} \right)^{3}, \quad 0 \le r \le h$$

1. 核函数 $W_{poly6}$：Poly6 是求密度时最常用的权重函数，名字源于表达式里 $(h^2-r^2)^3$ 的六次结构。

$$\displaystyle W_{poly6}(r, h)$$

2. 归一化系数 $\dfrac{315}{64\pi h^{9}}$：这个常数让核函数在全空间积分为 1，保证质量在加权求和时不增不减。

$$\displaystyle \frac{315}{64\pi h^{9}}$$

3. 光滑长度 $h$：影响半径，r 超过 h 时核函数直接取 0，邻域搜索只关心半径 h 内的粒子。

$$\displaystyle h$$

4. 中心距 $r$：目标点与粒子之间的距离，距离越小核函数值越大。

$$\displaystyle r$$

5. 幂次项 $(h^2-r^2)^3$：随距离增大，影响按三次方平滑坠落到 0，边缘无尖角，数值稳定性好。

$$\displaystyle \left( h^{2} - r^{2} \right)^{3}$$

花椒放得太多，整锅都是麻味；放得太少，又尝不出来。h 选大了，流体细节糊成一团；选小了，粒子各自为政、互相找不到邻域。核函数与光滑长度，就是这锅汤的“调料尺度”，决定了 SPH 模拟的细腻程度和计算成本。

密度计算：一勺汤的浓稠度

你要判断汤的哪个位置浓稠，最直接的办法就是舀一勺，看这一勺里食材挤了多少。SPH 也是这样——想知道粒子 i 所在处的密度，就把它邻域内的“食材”按质量数一遍，离得越近贡献越大，这就是密度求和公式的来源。

在连续介质力学里，密度通常由连续性方程（质量守恒的偏微分方程）求解。而 SPH 的聪明之处在于：密度是“数”出来的，压根不需要解方程。每个粒子的质量为 m_j，它像一颗小葡萄干，把质量“摊”在核函数的范围内，于是粒子 i 处的密度就等于邻域内所有葡萄干摊到 i 处的质量之和：

$$\displaystyle \rho_i = \sum_{j} m_j \, W(\mathbf{r}_i - \mathbf{r}_j, h)$$

1. 目标密度 $\rho_i$：粒子 i 所在位置的流体密度，SPH 里每个粒子每一帧都要重新计算。

$$\displaystyle \rho_i$$

2. 邻域求和 $\sum_j$：对光滑长度 h 范围内的所有邻居粒子遍历并累加。

$$\displaystyle \sum_{j}$$

3. 粒子质量 $m_j$：每个邻居粒子贡献多少“质量”，SPH 中质量通常恒定不变。

$$\displaystyle m_j$$

4. 邻距核权重 $W(\mathbf{r}_i-\mathbf{r}_j, h)$：距离越近质量分摊越多，距离超过 h 的粒子不参与求和。

$$\displaystyle W(\mathbf{r}_i - \mathbf{r}_j, h)$$

细节上需要注意两件事。首先，这个求和天然满足质量守恒：所有粒子的质量固定，加在一起总量不变，流体体积变化自然反映为密度升降，无需额外修正。其次，密度必须每帧重算——粒子位置时刻在变，“哪一勺浓、哪一勺稀”也跟着变。实际操作中，初始密度可能偏离目标密度 ρ0，这会在下一节产生压力，推动粒子重新分布。

就像大厨每搅拌一下都要重新判断哪里浓稠一样，SPH 每一帧都要重新“舀一勺”。密度不是一次性设定好的常数，而是粒子分布状态的实时写照——它是压力计算的前置条件，也是整个 SPH 循环的“温度计”。

压力与力：拥挤人群会推人

想象商场里人山人海：你周围挤满了人就会感到被往四面八方推搡，而稀疏的角落却没有人挤你。流体中的压力就是这种“微观推挤”的宏观描述。SPH 里压力不是凭空设定的，它由密度偏差驱动——哪个区域密度超过了理想值 ρ0，那里就产生了“拥挤压力”。

SPH 常用一个简单的状态方程（Tait 方程的近似）：p_i = k(ρ_i - ρ0)。k 是刚度系数，ρ0 是静水密度。密度低于 ρ0 时压力为负（轻微拉扯），高于 ρ0 时压力为正（向外推挤）。光有压力还不够，粒子感受到的是压力的“变化率”——压力梯度，就像人群的压力差让你被迫移动。离散化后，粒子 i 受到的压力梯度力为：

$$\displaystyle \mathbf{F}_i^{\text{pressure}} = -\sum_{j} m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W(\mathbf{r}_i - \mathbf{r}_j, h)$$

1. 压力梯度力 $\mathbf{F}_i^{\text{pressure}}$：作用在粒子 i 上、由邻域压力不均产生的合力，驱动流体从高压流向低压。

$$\displaystyle \mathbf{F}_i^{\text{pressure}}$$

2. 负号 $-$：力指向压力减小的方向——密度高的地方把粒子推往稀疏处，对应人群散开的直觉。

$$\displaystyle -$$

3. 对称压力项 $\left( \dfrac{p_i}{\rho_i^2} + \dfrac{p_j}{\rho_j^2} \right)$：同时考虑 i、j 双方的压力，交换下标后数值不变，天然满足牛顿第三定律。

$$\displaystyle \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}$$

4. 核梯度 $\nabla W$：核函数的梯度给出推力的方向和强弱，指向邻居方向且随距离衰减。

$$\displaystyle \nabla W(\mathbf{r}_i - \mathbf{r}_j, h)$$

5. 质量加权 $\sum_j m_j$：每个邻域粒子对力的贡献按质量累加，质量越大影响越大。

$$\displaystyle \sum_{j} m_j$$

除了压力梯度力，SPH 还叠加三种力。一是粘性力，相当于汤的“粘稠感”，它取邻域粒子速度差的加权平均，让相邻粒子速度趋向一致，抑制抖动；二是重力，最简单的外力，像地心引力让汤往锅底沉；三是边界力，粒子逼近容器壁时由惩罚函数产生推力，防止穿帮。把所有力求和除以质量，就得到粒子 i 的加速度 a_i = F_i / m_i。

汤太浓的地方，食材互相顶撞，把同伴推向清汤寡水的区域——这是压力的作用；汤太稀的地方，粒子又被“拉”回来填补空隙。压力梯度力、粘性力、重力与边界力四者配合，粒子们就在锅里自动演绎出翻涌、沉降、回旋的流体姿态。

邻域搜索：只数勺子里的食材

前面所有计算都有同一个前缀——“邻域内”。那么问题来了：怎么快速知道每个粒子的邻域是谁？最笨的办法是让每粒食材和其他所有食材两两打个招呼：一大锅汤有 N 粒食材，就要比 N² 次关系，这就是 O(N²) 的朴素遍历。当 N 从几百涨到几万，计算量会爆炸式增长。

聪明的大厨有自己的做法：先把锅底划成整齐的小方格，每个格子大约一个“光滑长度”见方。然后让每粒食材报出自己所在的格子编号，查找邻居时，只需要检查自己所在格子以及紧挨着的相邻格子就够了。这套思路在 SPH 中叫“空间哈希网格”（Spatial Hashing）。粒子的格子坐标可以表示为：

$$\displaystyle \text{cell}(\mathbf{r}) = \left( \left\lfloor \frac{x}{h} \right\rfloor ,\; \left\lfloor \frac{y}{h} \right\rfloor ,\; \left\lfloor \frac{z}{h} \right\rfloor \right)$$

1. 网格单元 $\text{cell}(\mathbf{r})$：空间位置 r 所归属的格子坐标，是整个邻域搜索的地图索引。

$$\displaystyle \text{cell}(\mathbf{r})$$

2. 向下取整 $\left\lfloor \cdot \right\rfloor$：位置除以 h 后取下整，保证坐标连续变化时格子编号整数跳跃。

$$\displaystyle \left\lfloor \frac{x}{h} \right\rfloor$$

3. 三维坐标分量：x、y、z 三个方向各算一次，组成完整的空间格子编号。

$$\displaystyle \frac{x}{h},\ \frac{y}{h},\ \frac{z}{h}$$

实际操作中，空间哈希算法会把三维格子编号压缩成一个一维索引（通过乘法和取模），塞进哈希表。这样每个粒子只需维护一串“邻居列表”，每次查询平均只需检查常数个格子。密度、压力、粘性力的计算全都复用这份邻居列表，整体复杂度从 O(N²) 降到接近 O(N)。选择格子边长时，一般取刚好等于光滑长度 h——既能保证不漏邻居（粒子影响范围不超过 h 个格子），又不会让格子大而无当。

就像大厨不会把整锅食材都倒进碗里来数，SPH 通过空间哈希把自己的“社交圈”锁定在相邻几格内。邻域搜索不是锦上添花，而是让 SPH 从玩具演示走向百万粒子工业模拟的关键一跃。

时间积分：稳稳地一步步推进

掌握了每个粒子受到的力，接下来就是“让汤真正动起来”。物理模拟的本质，是用一个个微小的时刻去逼近连续的时间。我们把时间切成固定的小段 Δt，每段做一次“力→加速度→速度→位移”的推进，这就是时间积分。半隐式欧拉是最常用的起步方法：先用当前力更新速度，再用新速度更新位置——先算速度、再算位置，每一步都踩在前一步的“新鲜”结果上：

$$\displaystyle \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + \mathbf{v}_i(t + \Delta t) \Delta t, \qquad \mathbf{v}_i(t + \Delta t) = \mathbf{v}_i(t) + \frac{\mathbf{F}_i(t)}{m_i} \Delta t$$

1. 位置更新 $\mathbf{x}_i(t+\Delta t)$：新位置 = 旧位置 + 新速度 × 时间步长，速度驱动粒子走向下一帧。

$$\displaystyle \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + \mathbf{v}_i(t + \Delta t) \Delta t$$

2. 速度更新 $\mathbf{v}_i(t+\Delta t)$：新速度 = 旧速度 + 加速度 × 时间步长，力决定粒子如何加速。

$$\displaystyle \mathbf{v}_i(t + \Delta t) = \mathbf{v}_i(t) + \frac{\mathbf{F}_i(t)}{m_i} \Delta t$$

3. 合力与质量 $\dfrac{\mathbf{F}_i(t)}{m_i}$：把上节算出的压力、粘性、重力等合力除以质量，得到加速度。

$$\displaystyle \frac{\mathbf{F}_i(t)}{m_i}$$

4. 时间步长 $\Delta t$：每步推进的时间，过大则数值失稳（粒子飞散），过小则计算量剧增，通常受 CFL 类似条件约束，与光滑长度和最大速度相关。

$$\displaystyle \Delta t$$

时间步长是模拟稳定的命门。直觉上，Δt 必须小到“粒子在一步之内不会冲过整个光滑长度 h”——否则它会直接穿过邻居的“影响圈”，力还没来得及作用就已经飞走了，这就是数值失稳（模拟爆炸）。实践中常把 Δt 限制为 h 除以最大速度的若干分之一。除此之外还有最后一个环节：边界处理。容器壁就像锅沿，不能放食材跑出去。常用做法是当粒子越过边界时立刻反弹回内部，并附加轻微阻尼。

把整锅汤的周期串起来——每帧先做邻域搜索找到“勺子里的食材”，再算密度判断浓稠，用状态方程算压力、叠加各种力求合力，最后按时间步长推进速度和位置。这个循环就像大厨的节拍：尝一口、搅一搅、等一等，反反复复，一锅活灵活现的流体就在屏幕上沸腾起来了。

AI 应用场景

Python 演示

"""A simplified 2D SPH fluid simulation for teaching purposes."""import mathimport random  class SPHSimulation:    """Use particles to simulate fluid: density, pressure, force, integrate."""     def __init__(self, count=50):        self.count = count        self.smooth_len = 1.2          # smoothing length h, also the neighbor radius        self.target_density = 2.0      # rest density rho0        self.pressure_scale = 1.5      # stiffness coefficient k        self.mass = 1.0                # mass of each particle        self.time_step = 0.016         # time step dt        self.particles = []     def reset(self):        self.particles = []        for _ in range(self.count):            self.particles.append({                "x": random.uniform(0.5, 6.5),                "y": random.uniform(0.5, 3.5),                "vx": 0.0,                "vy": 0.0,                "density": 0.0,                "pressure": 0.0,            })     def kernel(self, r, h):        """Poly6 kernel: larger weight closer to the center, zero beyond h."""        if 0 <= r <= h:            q = h * h - r * r            return q * q * q / (4.0 * h * h * h)        return 0.0     def kernel_gradient(self, dx, dy, r, h):        """Gradient of the kernel, pointing to the neighbor direction."""        if 0 < r <= h:            scale = -3.0 * (h * h - r * r) ** 2 / (4.0 * h * h * h)            return scale * dx / r, scale * dy / r        return 0.0, 0.0     def compute_density(self):        """Step 1: sum kernel-weighted masses from all neighbors."""        for p in self.particles:            p["density"] = 0.0        for i, p in enumerate(self.particles):            for j, q in enumerate(self.particles):                if i == j:                    continue                dx = p["x"] - q["x"]                dy = p["y"] - q["y"]                r = math.sqrt(dx * dx + dy * dy)                p["density"] += self.mass * self.kernel(r, self.smooth_len)            p["pressure"] = self.pressure_scale * (p["density"] - self.target_density)     def apply_forces(self):        """Step 2: pressure force and gravity, then update velocity and position."""        for i, p in enumerate(self.particles):            fx = 0.0            fy = -9.8            for j, q in enumerate(self.particles):                if i == j:                    continue                dx = p["x"] - q["x"]                dy = p["y"] - q["y"]                r = math.sqrt(dx * dx + dy * dy) + 1e-6                gx, gy = self.kernel_gradient(dx, dy, r, self.smooth_len)                scale = p["pressure"] / (p["density"] ** 2) + q["pressure"] / (q["density"] ** 2)                fx += -self.mass * scale * gx                fy += -self.mass * scale * gy            p["vx"] += fx * self.time_step            p["vy"] += fy * self.time_step            p["x"] += p["vx"] * self.time_step            p["y"] += p["vy"] * self.time_step            if p["x"] < 0.5 or p["x"] > 6.5:                p["x"] = max(0.5, min(6.5, p["x"]))                p["vx"] *= -0.5            if p["y"] < 0.5 or p["y"] > 3.5:                p["y"] = max(0.5, min(3.5, p["y"]))                p["vy"] *= -0.5     def run(self, steps=30):        """Main loop: density -> forces -> integration, repeated."""        self.reset()        for _ in range(steps):            self.compute_density()            self.apply_forces()        return self.particles  if __name__ == "__main__":    sim = SPHSimulation(count=40)    result = sim.run(steps=30)    avg = sum(p["density"] for p in result) / len(result)    print("SPH simulation done, particle count:", len(result))    print("Average density:", round(avg, 4))

1. def kernel(self, r, h):定义 Poly6 核函数，距离越近权重越大，超过 h 返回 0。
2. q = h * h - r * r计算核函数的自变量，这是 Poly6 幂次项的核心。
3. return q * q * q / (4.0 * h * h * h)核函数值，数值平滑且在边界处归零。
4. p["density"] += self.mass * self.kernel(r, self.smooth_len)密度是邻域粒子质量的核加权累加，对应文章第 3 节的公式。
5. p["pressure"] = self.pressure_scale * (p["density"] - self.target_density)状态方程：密度高于静水密度则压力为正。
6. scale = p["pressure"] / (p["density"] ** 2) + q["pressure"] / (q["density"] ** 2)对称压力项，保证牛顿第三定律。
7. fx += -self.mass * scale * gx压力梯度力累加，负号让力指向压力减小方向。
8. p["vx"] += fx * self.time_step半隐式欧拉：先更新速度。
9. p["x"] += p["vx"] * self.time_step再用新速度更新位置。
10. p["vx"] *= -0.5边界反弹并施加阻尼，防止粒子逃出计算区域。
11. self.compute_density(); self.apply_forces()主循环两步走，对应文章第 6 节的模拟节奏。

常见误区

面试题