简介这份基于 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 版本就是从这种“小水花”里一步步爬出来才最终遮住一片水池的。如果你的目标不只是还原一个视觉效果还想理解背后的计算原理这套从小到大的递进路径会帮你省下大量排查问题的力气。本文还有配套的精品资源点击获取