深入理解流体模拟(SPH/PBD)

可视化演示

两大流派概览

泳池派对里，自由散漫型玩家靠“第六感”自动避让，而纪律严明型玩家服从哨声调位。这两种风格，正好对应流体模拟界的两大门派：SPH 与 PBD。游戏物理引擎中，它们就像武术里的刚柔两派——一个讲究实时传感下的随性而为，一个讲究规则约束下的秩序井然。假如今天你要为《原神》里的一汪湖水写物理，你需要先选好门派，因为后续所有代码都将沿着这条路线生长。

从更深一层看，无论哪一派，都是在求解同一组流体运动方程——纳维-斯托克斯方程。SPH 用“拉格朗日视角”把连续水体切成大量携带质量、速度的粒子，让流体性质附着在粒子上随水流动；而 PBD 则绕开显式积分，把“密度不能超限”翻译成约束条件，通过反复修正粒子位置来逼近真实流动。前者更像生物课上的“跟随模型”，后者更像工程课上的“公差配合”。

两者的差异最终体现在数学上：SPH 直接离散动量方程，PBD 离散约束投影。以下是流体运动方程的简化形式：

表示单位体积流体的动量变化率，本质是牛顿第二定律在连续介质上的体现——力越大，流速变化越剧烈。

方向从高压指向低压，驱动流体从拥挤区域涌向空旷区域，就像泳池里人太挤会自然向人少的地方移动。

描述内摩擦力对速度场的平滑作用，它抑制涡流和抖动，让水流看起来“柔顺”而不是“炸裂”。

代表重力、浮力、玩家交互力等外部作用，游戏里你踩水时的反作用力就挂在这一项上。

回到泳池比喻：SPH 玩家手臂一伸，靠自己感知的压力差挪动身体；PBD 玩家则听着哨声，每次哨响都对照规则挪到合规位置。方法不同，但都能让派对不失控——选谁，取决于你的派对需要精细的水花还是稳定的秩序。

核函数与密度

第六感不是玄学，它是有作用范围的。就像你只能察觉手臂长度内的同伴，SPH 里的核函数 $W(\mathbf{r}, h)$ 也定义了一个“感知半径” $h$，距离超过 $h$ 的邻居完全不参与计算。这就是 SPH 中“光滑”二字的来历：每个粒子的物理量不再是一个孤立的点，而是被核函数抹开成一片柔和的“影响力场”。

在设计核函数时，工程师要求它满足三个条件：归一化（全场积分等于 1）、紧致支持域（超过 $h$ 就归零）、以及随距离平滑衰减。满足这些条件的函数可以有很多种，游戏界最常用的是 Poly6（用于计算密度）、Spiky（用于压力梯度）和 Viscosity（用于粘性力）。它们的分工就像派对里的三样工具：卷尺量距离、安全员推人、润滑剂防磕碰。

有了核函数，密度就顺理成章地变成加权求和。粒子 $i$ 感受到的局部密度，等于所有邻居粒子质量乘以核函数值再累加：

遍历所有邻居粒子 $j$，但核函数会自然截断距离超过 $h$ 的粒子，不必担心遍历到全场景。

每个粒子携带的质量，游戏中通常恒定；它是加权求和的基础权重，代表每个同伴“有多重”。

以两个粒子的距离为输入，距离越近输出越大，且满足归一化与紧致支持域性质。

最终得到的是粒子 $i$ 私人视角的局部密度，而不是整个流体的平均密度——每个人感受到的拥挤各不相同。

泳池里，每个人都在心中默默估算：我这 1 平方米站了几个人？核函数就是那把“估算尺”，密度则是算出来的“重量感”。尺子的刻度（$h$）选多大，直接决定了你是只注意贴脸的朋友，还是能感受到半个泳池的拥挤。

邻居搜索策略

当你想知道周围有多少人时，最笨的办法是全场 100 人挨个问一遍。对于 10 万粒子的流体模拟，这种“全场咨询”意味着每秒数十亿次距离判断，再强的 GPU 也会被拖垮。聪明的做法是：先把泳池划分成边长等于 $h$ 的小方格，你只需要低头看看自己站在哪个格子里，再向相邻的 8 个格子张望一眼就够。

这个“抬头看格子”的操作，在算法里叫空间哈希网格。每个粒子根据位置可以快速算出自己所在格子的编号：

将粒子 $x$ 坐标除以格子边长并向下取整，得到它所在的列索引。

同理得到行索引，与列号共同构成二维网格坐标。

行号乘以总列数后与列号相加，把二维坐标压成一维数组下标，方便存入连续内存。

当前格加周围 8 格就是全部可能邻居的候选集，将复杂度从 $O(n^2)$ 降到 $O(n)$。

泳池管理员还会更进一步：把每个格子里的玩家按顺序排好队（桶排序），这样当你喊“第三格的人举手”时，GPU 可以一次性连续读到一整块名单，而不是满世界乱翻。这正是空间网格 + 桶排序这对组合在 GPU 上高效的原因——世界一下子从 100 人缩小到 7、8 人，而且名单就在手边。

PBD约束求解

PBD 的玩家不关心你推我的力有多大，只关心一条铁律：每平方米人数不能超过上限。哨声一响，所有人按规则自动调整站位。物理学家把这条“上限”写成密度约束：每个粒子的局部密度必须回归到静水密度。于是，流体问题被转化成一个“位置修正”问题。

PBD 的求解过程分三步走：先按当前速度预测新位置，再检查约束是否被违反，最后沿约束梯度方向修正位置。其中最关键的一步是，密度约束函数被定义为：

由第 2 节的核函数密度估计得到，描述粒子 $i$ 感受到的实际拥挤程度。

水的静止参考密度，可以理解为派对场地的“额定容量”，超过它就代表超员。

实际密度与额定密度的比值再减 1，得到一个无量纲的“超员指标”。

PBD 只惩罚“过于拥挤”的情况，不强行限制“过于稀疏”。

找到违反约束的粒子后，还需要知道往哪个方向修正。数学上，这需要计算约束的梯度，然后沿负梯度方向投影：

指向密度增长最快的方向，也就是“最容易造成拥挤的方向”。

沿梯度的反方向移动，才能让密度下降，把粒子从拥挤区“推”向空旷区。

把修正量缩放到安全尺度，避免一步跨过头导致新的违反。

最终得到的位置偏移，将它加到粒子上，即完成一次约束投影。

PBD 的妙处在于：约束求解可以只迭代 2~4 轮就收手，允许大时间步且几乎不会爆炸。派对管理员不需要精确计算每人的受力，只需喊几轮“那边的人往右挪两步”，全场就能快速恢复秩序。这正是 PBD 在移动端游戏和 VR 交互场景中广受欢迎的原因——稳、快、可控。

GPU并行实现

假如派对规模从 100 人变成 100 万人，你不可能再靠几个管理员来回跑。你得给每个人配备一个“分身管理员”，让他们同时工作——这正是 GPU 并行计算的本质。每个 GPU 线程负责一个粒子，所有线程同时计算密度、压力和位置更新，理论上速度可以提升几十上百倍。

但并行计算有个绕不开的麻烦：写冲突。当两个线程同时往同一个密度累加器里写入数据时，结果会互相覆盖。GPU 为此提供了原子操作（atomicAdd），让多个线程排队写入，保证结果正确。同时，为了减少全局内存的反复访问，聪明的工程师会把邻居粒子索引缓存到共享内存中，让同一线程块内的线程快速复用。

GPU 上并行计算密度的过程，可以看作先分块求和、再跨块汇总：

将粒子 $i$ 的邻居集合按 GPU 线程块划分成 $B$ 个互不相交的子集。

每个线程块内的线程并行计算各自邻居的贡献，再通过块内归约得到子集和。

将各线程块的部分和汇总，得到完整的密度值——这一步通常需要少量原子操作或全局同步。

同时维护新旧两份位置数据，所有线程读旧数据、写新数据，避免读到半更新的中间状态。

派对策划的智慧在于：让管理员们先在各自小区域内统计人数（共享内存），再汇总到总指挥；新旧场地分开（双缓冲），旧场地的信息不会干扰新场地的布置。并发高、排队少，派对效率自然直线上升。

性能参数调优

派对开久了，总有人嗨过头，突然加速乱跑。如果管理员的响应节奏跟不上，就会撞成一团。在流体模拟里，这个“响应节奏”就是时间步长 $\Delta t$。步长取得太小，模拟一帧要跑几十次，性能告急；步长取得太大，粒子一步跨出老远，数值爆炸，画面瞬间变成“雪花屏”。

幸运的是，数值分析早就给出了一个经验判据：CFL 条件。它指出时间步长的上界与核函数半径和粒子最大速度相关：

物理仿真每次推进的时间间隔，直接影响模拟精度与性能的平衡。

粒子的感知半径，也是空间网格的格子边长，决定了空间分辨率。

当前所有粒子中速度的最大值，代表“跑得最快的那个人”。

经验常数，给预测留足余量，防止最坏情况下的数值失稳。

除了控制时间步长，工程上还有三招常用。其一是固定物理帧率：物理系统固定在 60Hz 更新，渲染层插值插帧，避免步长突变；其二是 SoA 内存布局：把粒子位置、速度分别放进独立数组，而不是塞进一个结构体数组，让 GPU 缓存命中率大增；其三是粒子简并 LOD：远处的水面用稀疏大粒子，近处才细化小粒子，像派对里普通区少派管理员、VIP 区多派人手一样，预算有限但体验不减。

好的派对一定有时间表：每隔几分钟清点一次人数（固定 tick），场地按区域划好动线（空间排序），重点区域多配管理员（LOD）。流程越规整，意外越少，这就是调优的全部秘密。

应用与AI融合

掌握了两种派对管理术，我们终于可以讨论它们如何为水上乐园、电影特效和数字孪生城市服务。在《原神》这类开放世界游戏中，玩家踩水、施法、扔炸弹，水面必须在毫秒级内做出反应，这背后是成千上万个 SPH 粒子在 GPU 上疯狂运算；而在移动端切水果游戏里，水果被切开瞬间的果汁飞溅，则更多依赖 PBD 用极低成本“装出”液体的样子。

更高阶的应用来自 AI 与物理模拟的联姻。传统的 SPH 每一步都要遍历邻居、算核函数、求压力梯度，计算量巨大。研究者开始尝试用神经网络预测粒子的粗粒度速度场——让网络承担 80% 的“大方向”预测，只在需要精细细节的局部区域跑完整 SPH。这种“先预测、后修正”的思想，与 PBD 的约束投影异曲同工：

粒子在第 $n$ 帧的位置，是预测的起点。

用当前速度乘以时间步长，估算出粒子在下一帧的大致位置。

未经修正的“猜测位置”，已经包含惯性运动趋势。

把预测位置拉回合法状态，使其满足密度、边界等物理约束。

修正完毕的合法位置，进入下一帧循环。

未来的派对策划也许不用挨个讲规则：让一个看过无数派对录像（训练数据）的智能体提前预判哪里会挤爆，只派管理员去关键区域现场协调。SPH 与 PBD 依然是舞台上的演员，而 AI 正在成为那个运筹帷幄的总导演。

AI 应用场景

Python 演示

import mathimport random H = 0.3        # 核函数支持半径M = 1.0        # 粒子质量RHO0 = 100.0   # 静水密度DT = 0.0005    # 时间步长GRAVITY = 9.8  # 重力加速度  def kernel_w(r):    # 2D Poly6 核函数, 用于密度估计    if r >= H:        return 0.0    diff = H * H - r * r    return 4.0 * diff ** 3 / (math.pi * H ** 8)  def kernel_grad(r):    # 2D Spiky 核函数梯度, 用于压力计算    if r >= H or r == 0.0:        return 0.0    return -60.0 * (H - r) ** 2 / (math.pi * H ** 6)  class SPHParticle:    def __init__(self, x, y):        self.x = x        self.y = y        self.vx = random.uniform(-0.5, 0.5)        self.vy = random.uniform(-0.5, 0.5)        self.density = RHO0  def simulate_step(particles):    # 第一遍: 计算每个粒子的密度    for p in particles:        total = 0.0        for q in particles:            dx = p.x - q.x            dy = p.y - q.y            r = math.sqrt(dx * dx + dy * dy)            total += M * kernel_w(r)        p.density = total     # 第二遍: 计算压力梯度并更新位置    for p in particles:        fx = 0.0        fy = -GRAVITY * M  # 重力作为外力        for q in particles:            if p is q:                continue            dx = p.x - q.x            dy = p.y - q.y            r = math.sqrt(dx * dx + dy * dy)            if r > 0 and r < H:                press = (p.density - RHO0) + (q.density - RHO0)                grad = kernel_grad(r)                fx += -press * grad * dx / r                fy += -press * grad * dy / r        p.vx += fx * DT        p.vy += fy * DT        p.x += p.vx * DT        p.y += p.vy * DT         # 边界反弹        if p.x < 0.0 or p.x > 1.0:            p.vx *= -0.5            p.x = max(0.0, min(1.0, p.x))        if p.y < 0.0:            p.vy *= -0.5            p.y = 0.0        if p.y > 1.0:            p.vy *= -0.5            p.y = 1.0  def main():    particles = []    for row in range(6):        for col in range(8):            x = 0.2 + col * 0.08            y = 0.7 + row * 0.05            particles.append(SPHParticle(x, y))     print("开始 2D SPH 简化模拟, 粒子数 = %d" % len(particles))    for step in range(300):        simulate_step(particles)        avg_density = sum(p.density for p in particles) / len(particles)        max_density = max(p.density for p in particles)        min_y = min(p.y for p in particles)        if step % 20 == 0:            print("step=%03d avg_density=%.1f max_density=%.1f min_y=%.2f"                  % (step, avg_density, max_density, min_y))    print("仿真完成, 粒子在重力下自然沉降并因压力互相推开")  if __name__ == "__main__":    main()

1. p.density = total用 Poly6 核函数累加所有邻居粒子的质量贡献，得到局部密度
2. press = (p.density - RHO0) + (q.density - RHO0)密度差之和近似压力源项，密度越高排斥越强
3. grad = kernel_grad(r)计算 Spiky 核函数梯度，它决定了压力推力的方向与强度
4. p.x += p.vx * DT用显式欧拉法更新粒子位置，时间步长受 CFL 条件约束
5. p.vx *= -0.5边界反弹系数，模拟粒子与容器壁的碰撞耗散

常见误区

面试题