简介：这份基于 Visual Studio 2010 与 OpenSceneGraph 3.4.1 的 SPH（平滑粒子流体动力学）流体仿真项目，面向图形学、物理模拟方向的开发者和学生，用于学习无网格流体的数值计算与三维可视化。它通过将流体离散为有质量粒子，借助加权函数求解密度、压力、动量与能量方程，可模拟液面飞溅、容器内水体等动态效果。压缩包共 40 个文件，约 43.37MB，包括 10 个 h 头文件与 8 个 cpp 源文件（实现粒子系统、时间积分、碰撞边界等），另有 VS2010 解决方案（sln/vcxproj）、OSG 场景文件、BMP 纹理和说明文档，整体目录清晰，便于直接编译学习。已有 430 人学习下载。项目不仅给出 SPH 核心算法与完整工程配置，还涉及密度估计、光滑核函数、粘性与表面张力、自适应时间步长等关键知识点，帮助读者理解从理论推导到 C++ 实现的完整链路，并利用 OpenSceneGraph 实时观察粒子运动与流体形态，适合希望深入流体模拟或开展课程设计的开发者参考。 前几年我在游戏项目里接到过一个需求：让一桶水泼到地上以后，地面能看到真实的水花飞溅，而不是一个提前做好的透明度动画。调研了一大圈之后，我把目标锁定在了一个词上——SPH，全称是 Smoothed Particle Hydrodynamics，中文一般翻译成“平滑粒子流体动力学”。这是一个典型用粒子做流体模拟的方法，把流体拆成成千上万个带着物理属性的小球，然后通过这些小球之间的相互作用，算出水花、波浪、烟雾、岩浆这些动态效果。
它最早出现在 1977 年，当时天文学家 Lucy、Gingold 和 Monaghan 提出这个思路，想解决的是恒星形成、星系碰撞这类没有规则边界的演化问题。后来到了九十年代，计算机图形学的研究者把它搬进了视觉特效领域，SPH 才成了游戏引擎和影视特效里绕不开的经典方法。不管你是做实时渲染、离线 CG、工程仿真，还是单纯研究数值方法，这套思路都值得完整过一遍。这篇文章我就把 SPH 的核心原理、一版能跑通的最小模拟器，以及我在调参过程中踩过的坑一起讲清楚。
1. SPH 到底是个啥：从“用网格还是用粒子”说起
如果你去问一个做流体模拟的人，最基础的问题一定是“你算的是网格还是粒子”。这个选择决定了整个技术路线的走向。欧拉网格方法习惯于把空间切成固定的小格子，在每个格子上记录速度、压力、密度，然后不断更新这些场。它的强项是处理大范围水体、稳定流动，数据天然有拓扑结构，方便做并行和矩阵迭代。但遇到强烈飞溅、液体分离、自由表面破碎这些场景，固定网格会被大量空气混入导致数值发散，处理起来非常痛苦。
SPH 走的是另一条路：不切网格，直接给流体“贴标签”。每个粒子代表一小团流体微元，携带质量、速度、密度、压力这些属性，跟着流体一起运动。宏观上的水花、涡旋、飞沫，本质上都是大量粒子共同运动的结果。这种视角叫拉格朗日视角，它天然适合处理自由表面和剧烈形变的问题，因为粒子想怎么飞就怎么飞，不会被网格束缚住。
打个比方，网格方法就像一个交警站在路口，只知道某块区域的车流量和平均速度；粒子方法则像是给每台车都装上了定位器，你能完整追踪每一台车的轨迹。应用到水花问题上，显然是后者更适合——你需要看到水珠散开、空中分离、落回水面后再弹起的小细节，这些细节恰恰是粒子方法最容易表现的。
至于适合谁看这篇文章，我建议三类人重点关注：一类是做游戏特效或者 TA 的，想在引擎里实现真实水体交互；一类是搞工程仿真的，想评估 SPH 是否能处理自己领域里的自由表面问题；还有一类纯粹是学数值计算的学生，拿 SPH 当个入门的无网格方法案例来学习理解。下面的内容我尽量从原理讲到实现，再讲到调参经验，让你读完就能动手搭一个自己的流体模拟器。
2. 核心原理拆解：粒子、核函数和三条关键公式
2.1 一切从核函数开始
SPH 的底层逻辑非常朴素：既然流体是连续的，那任意位置上的物理量就应该由它周围一小片区域里的物质共同贡献。怎么把“周围一小片”的贡献算进来？靠的就是核函数 W。核函数是一个只和距离有关的函数，作用有点像加权平均里的权重：离目标点越近的粒子，权重越大；超过某个半径之后，权重直接归零。这个半径通常记作 h，也就是平滑长度。
核函数必须满足几个基本条件：归一性（积分和为 1）、紧凑支撑（在 h 之外为零）、还有足够的平滑性。归一性保证插值不会引入系统偏差，紧凑支撑则保证了计算量可控，不会每个粒子都要和全场所有粒子发生关系。实际做模拟的时候，常用的核函数有三个：Poly6 核适合算密度，因为它处处光滑，形式简单；Spiky 核适合算压力梯度，因为它在粒子靠近时会产生比较大的排斥力，能有效防止粒子互相穿透；粘性核则专门用来算 Laplace 算子，让粘度在近距离处更明显。从我个人的测试经验来看，密度用 Poly6、压力用 Spiky，基本是手工实现 SPH 的标配组合，很少有需要改的时候。
2.2 密度和压力：两个最核心的计算量
SPH 里没有显式的连续性方程去演算密度，而是直接通过粒子分布来“统计”密度。对于粒子 i，它的密度可以写成：
ρ_i = Σ_j m_j · W(|r_i - r_j|, h)这个公式的含义是：把周围所有邻居粒子的质量，用一个距离相关的权重累加起来，就得到了粒子 i 所在位置的流体密度。你不需要额外求解质量守恒方程，只要粒子数量不变，体系总质量就是守恒的，这是 SPH 非常讨喜的一点。
有了密度，压力就好算了。实际工程里很少去求解一个完整的不可压缩压力方程，大部分实现用的是状态方程，也叫弱可压缩模型：
P_i = k · (ρ_i - ρ_0)其中 ρ0 是初始参考密度，k 是刚度系数。这个公式非常简洁：密度偏离初始值越多，压力越大，粒子就会拼命往密度低的地方跑，从而恢复体积。k 的取值很敏感，取小了流体像气体一样软绵绵，取大了时间步长必须降得很小，否则数值直接爆炸。这个部分的具体调参经验，我在后面专门用一节来讲。
2.3 三股力：压力、粘性和重力
粒子在流体里的受力可以拆成三部分。首先是压力梯度力，它让流体从高压区域往低压区域流动，但在 SPH 里直接对压力场求梯度并不方便，所以一般用对称形式来保证动量守恒：
a_pressure = -Σ_j m_j · (P_i / ρ_i² + P_j / ρ_j²) · ∇W_ij其次是粘性力，它让相邻粒子的速度趋于一致，也就是流体天然有“搅匀速度差”的倾向：
a_viscosity = μ · Σ_j m_j · ((v_j - v_i) / ρ_j) · ∇²W_ij粘性系数 μ 越大，流体看起来越稠，比如蜂蜜和泥浆；μ 越小，流体越活泼，比如清水和酒精。最后就是重力，直接给所有粒子加同一个加速度向量就行。把这些力加起来，用牛二定律算出每个粒子的加速度，再积分更新速度和位置，整个模拟就在时间轴上一帧一帧地走下去了。
3. 手把手写一个 2D 水花模拟器
3.1 数据结构与初始化
先声明一下，我下面给的是一份最小可行的 CPU 实现，语言不重要，重点是流程。粒子结构体大概长这样：
struct Particle { float x, y; // 位置 float vx, vy; // 速度 float rho; // 密度 float p; // 压力 float ax, ay; // 加速度 };初始化的时候，把流体区域按规则的网格排列粒子，而不是随机撒点。随机撒点会带来密度场的大幅波动，模拟一开始就会出现局部高压、粒子乱飞的现象。我自己第一次尝试就是图省事用随机分布，结果前 20 帧就有大量粒子飞到了十万八千里外，排查了很久才意识到是初始化密度不均匀的问题。
粒子之间的初始间距记作 dp，每个粒子质量通常设成相同值，比如 m = ρ0 · dp²（在 2D 情况下）。这样整个水体的密度从一开始就接近 ρ0，压力项不会产生虚假的巨大梯度。
3.2 邻居搜索：从暴力到空间哈希
如果你直接两两配对判断邻居，算法复杂度是 O(N²)，几千个粒子还能忍，几万几十万个粒子就直接卡死。所以正规的 SPH 模拟里一定会做邻居搜索。最常见的方案是空间哈希，基本思路是把整个模拟区域分成格子，格子边长取核函数的支撑半径 h，这样每个粒子只需要查它自己所在的格子以及周围相邻的格子，就能找到所有可能的邻居。
在 2D 情况下，每个粒子只需要检测周围 3×3 个格子；3D 情况下是 3×3×3 个格子。对于粒子规模在十万以内的模拟，这是性价比最高的方案。实现空间哈希时要注意一个细节：格子边长不要取得比 h 小太多，否则同一个粒子的邻居会分布在很多个格子里，索引起来反而麻烦；也不要取得比 2h 大太多，否则每个格子里的粒子过多，依然退化成暴力搜索。我的经验是格子边长取 h 或者 1.1h 最稳。
3.3 完整的更新伪代码
下面这一段是整个模拟器的骨架，一共六步，每帧（每个时间步）循环执行：
for step in range(total_steps): // 1. 用空间哈希重建网格，收集每个粒子的邻居 build_hash_grid() find_neighbors() // 2. 计算每个粒子的密度 for each particle i: rho_i = 0 for each neighbor j: rho_i += m_j * W_poly6(x_i - x_j, h) // 3. 更新压力 for each particle i: p_i = k * (rho_i - rho_0) // 4. 计算加速度：压力项 + 粘性项 + 重力 for each particle i: a_i = (0, -g) for each neighbor j: a_i += -m_j * (p_i/rho_i² + p_j/rho_j²) * gradW_spiky(x_i-x_j, h) a_i += mu * m_j * (v_j - v_i) / rho_j * lapW_viscosity(x_i-x_j, h) // 5. 半隐式欧拉积分更新速度和位置 for each particle i: v_i += dt * a_i x_i += dt * v_i // 6. 处理边界，防止粒子飞出模拟区域 handle_boundaries()这个流程看着简单，但每一步都有隐藏的难度。密度计算用的是 Poly6 核，因为它公式简洁、数值平滑；压力梯度用的是 Spiky 核，因为它随距离递减更快，排斥力更“硬”，粒子不容易互相穿透。积分方式我建议先用半隐式欧拉，也就是先更新速度再用新速度更新位置，稳定性比显式欧拉好不少，实现成本却几乎为零。
3.4 边界处理的经验之谈
模拟区域边界如果不处理，粒子会在重力作用下直接落到地面以下，然后逐渐堆积、穿透、形成不可控的数值膨胀。最简单的边界方案是惩罚力：当粒子离墙太近时，施加一个沿法线方向的排斥力，距离越近，力越大。比如对地面 y = 0，可以写成：
if y_i < r0: a_y += spring * (r0 - y_i) - damping * vy_i这里的 r0 是粒子半径或者一个很小的阈值，spring 是两个系数。弹簧系数越大，粒子被推回越快，但过大又会让粒子在边界反弹得过于剧烈，像撞到蹦床一样。阻尼系数的作用是吸收法向速度，让粒子落在边界上时能安静下来。实际调参时，我习惯先让 water 落到地面上不动，再逐步加大扰动，一步一步观察边界的表现。
另外，还有一种更物理的做法是用一层“虚拟 ghost 粒子”填充边界区域，让内部粒子自然感受到墙壁的排斥作用。这个方法稳定性更好，但需要额外生成和处理虚拟粒子，初次实现时不建议上来就搞，惩罚力已经能解决大部分场景的问题了。
4. 调参避坑：一个稳定水花背后藏着的那几个数
4.1 平滑长度 h 是第一个要小心的数
SPH 的很多行为都由 h 决定。h 如果太大，每个粒子的邻居会非常多，流体被“抹”得过平，水花细节全丢，计算量也会明显增大；h 如果太小，邻居数量不足，密度场会出现严重振荡，直接导致模拟崩溃。工程里常用的经验是让 h 约为初始粒子间距的 1.2 到 1.5 倍：
h ≈ 1.3 · dp这个值既保证每个粒子周围有一定数量的邻居，又不至于让影响范围过大。在 2D 情况下，1.3dp 大概能带来 10 到 15 个有效邻居粒子，足够让密度场平滑稳定。如果你的场景需要更细的飞沫，应该是去减小 dp（也就是增加粒子数），而不是把 h 调小，两者不是同一个粒度上的事情。
4.2 时间步长：为什么“看着还行”突然就崩了
弱可压缩 SPH 的时间步长限制比一般想象中苛刻得多。因为它引入了人工声速，信息在粒子之间的传播速度变得很快，时间步长必须满足 CFL 条件，经验公式是：
dt < 0.4 · h / c_max其中 c_max 是波速上限，通常和压力刚度系数 k 有关。你如果把 k 调大来让流体更“硬”，就必须同步把 dt 缩小，否则模拟会在几十个时间步内直接爆掉。此外还有一个基于最大加速度的经验约束：
dt < sqrt(0.1 · h / max_accel)很多初学朋友遇到“粒子飞散”第一反应就是把 h 或者 k 乱调，其实大概率是 dt 超标了。我自己的做法是先在每个时间步里扫一遍所有粒子的速度和加速度，计算出推荐 dt，再乘以一个 0.8 的安全系数，上线之后基本不会因为时间步长出问题。
4.3 刚度和粘性：流体“手感”的平衡杆
刚度系数 k 决定了流体抵抗压缩的能力。k 太小，水柱落地后会像面团一样摊开甚至不反弹；k 太大，模拟会变得像弹球游戏，水面轻轻一碰就剧烈震荡，甚至粒子炸开。工程上，我通常会先用最小可接受的水花效果去标定 k，再逐步增加，直到找到“水落地后弹起但不碎成雾”的临界值。这个过程有点玄学，但对最终效果影响非常大。
粘性系数 μ 则是另一种手感控制器。清水需要较小的粘性才能有那种活泼的水花，但太小的话自由表面会产生微小的高频抖动，画面会显得很脏。我的习惯是：先用稍大的粘性跑稳，再逐步减小看效果哪里先出现抖动，把 μ 停在“抖动刚要出现还没出现”的位置。特别提醒一下，μ 过大会让流体看起来像糖浆，哪怕物理上很稳定，视觉上也完全不对，所以这个值最后一定要放到真实渲染里去确认，而不是只看模拟时的粒子位置。
5. 常见问题排查速查表
很多事情光讲原理不够，真到自己写代码的时候，犯的错误往往是类型完全不同的。这里把我在开发中最常遇到的几类问题直接整理出来，碰到症状直接对照找原因就行。
| 现象 | 最可能的原因 | 处理办法 |
|---|---|---|
| 模拟刚开始粒子就炸飞 | 初始化密度不均匀 / dt 太大 | 规则网格排布粒子；缩小 dt 到 0.2 倍再试 |
| 水体像棉花一样软塌塌 | 压力刚度 k 太小 | 增大 k，同时相应缩小 dt |
| 粒子互相穿透严重 | Spiky 核使用不当或 h 太小 | 确认压力梯度用的是 Spiky 核；提高 h 到 1.2dp 以上 |
| 水面高频抖动、画面脏 | 粘性太小 | 增大 μ 观察抖动是否消失 |
| 粒子堆积在边界上不反弹 | 边界惩罚力 spring 太小 | 增大 spring，并加适当的阻尼 |
| 粒子数量一多就非常卡 | 邻居搜索用了暴力法 | 换成空间哈希，检查格子尺寸是否约等于 h |
| 水花飞得太整齐，没有随机感 | 初始配置过于规则 | 在初始化位置上加很小的随机扰动 |
其中“粒子刚开始就炸飞”是新手最容易遇见的状况，也是最劝退的。排查顺序应该先看 dt，再看初始化密度，最后检查邻居搜索是否正确。如果 dt 已经缩小到理论安全的十分之一还崩溃，那大概率不是时间步长的问题，而是邻居列表本身找错了——请注意空间哈希在查询邻居时要算上周围所有格子，漏掉任意一个格子都会导致密度少算，形成巨大的虚假压力梯度。
6. 扩展方向：SPH 不只用来“溅水花”
6.1 游戏与影视里的降级与增强策略
在游戏引擎里，SPH 通常不会直接跑几十万个完整粒子，因为性能撑不住。常用的做法是分两套系统：一套是数量较少的大粒子，负责模拟宏观水体的走势；另一套是渲染层的飞沫粒子，从大粒子运动轨迹中抽取速度、位置和生命周期，只负责视觉效果。这种做法能在保证水花观感的同时，把粒子数量压到几千，手机和低端 PC 也可以流畅运行。
影视 CG 则更倾向于先跑高精度 SPH 模拟，再用 meshing 算法把粒子云提取成三角形网格，最后在渲染器里加材质和贴图。这一步叫 surface reconstruction，常见方案是 marching cubes 或者更平滑的起泡算法。如果你是从游戏转到影视方向，反而会觉得渲染阶段比模拟阶段还复杂，需要补不少以“液体形状到网格”的相关知识。
6.2 工程与科研领域：经典仍然能长寿
工程仿真领域里，SPH 特别擅长处理大变形和自由表面问题，被广泛应用在金属加工时的飞溅、液滴碰撞、船体与波浪相互作用、血液灌注等场景。原因很简单：网格方法遇到大变形通常需要不断重剖分网格，代价极高；SPH 不需要重剖分，流体跑到哪里，粒子就追随到哪里。
而在早期天体物理模拟中，SPH 的位置也非常稳定。因为星际尘埃、分子云碰撞、双星并合这些问题完全没有天然边界，且往往伴随极端的密度对比和大范围流动，SPH 的粒子表达天然适合这种稀疏环境。当然它也有缺点：捕捉激波需要额外引入人工粘性，流体混合界面的锐利度一般不如高精度网格方法。了解这些局限性其实比了解优势更重要，我在实际项目里经常遇到有人拿着 SPH 去做它不擅长的事，比如高马赫数可压缩流，结果效果不佳又回头怪方法不行，这其实是选型问题。
最后再分享一个我吃了很多次亏之后的经验教训：第一次写 SPH，千万别一上来就追求几万粒子的大水花。先把场景缩小到一千个粒子左右，用一个最简单的 2D 水箱塌落实验去验证密度、压力、粘性这几项的数值行为，确认粒子不会炸、水不会穿透边界之后，再逐步加粒子、加场景复杂度。我的第一个能跑起来的 3D 版本，就是从这种“小水花”里一步步爬出来才最终遮住一片水池的。如果你的目标不只是还原一个视觉效果，还想理解背后的计算原理，这套从小到大的递进路径会帮你省下大量排查问题的力气。
本文还有配套的精品资源，点击获取
