SPH in a weekend

做物理不懂流体，就好比打 dota 不玩卡尔，缺少一种品味的美。于是花了几天时间补了一下，虽然玩的烂但是好歹能放出民工三连的水平。

读完以后感觉粒子这个东西是真的好用啊。凡是解析方法要搞一堆积分的东西，换成粒子瞬间就简单好多，梯度也好算了。

已经有非常多介绍 sph 算法推导过程的文章了，最原始的论文链接也放在这，就只简单讲下算法的思路了，网页上敲公式也是累死人。
Matthias Muller, SIGGRAPH 2003

简而言之，原先你要模拟流体，就需要求解大名鼎鼎的 Navier-Stokes 方程，长这样：
$\rho(\frac{\partial \pmb{v}}{\partial t} + \pmb{v} \cdot \nabla \pmb{v})=-\nabla p + \rho \pmb{g} + \mu \nabla^2 \pmb{v}$

这是一个 force-based 模型，右侧三项分别是 $f_{pressure}$ ,, $f_{external}$ , $f_{viscosity}$ 。左侧是动量对 t, x, y, z 的全导数拆解，等价于动量的变化。简单来说就是 $mv=ft$ 的流体版。

现在我们换粒子解决这个方程的求解，简单方便。

模拟上存在两个视角，Euler 和 Lagrangian。Euler methods 求解连续的场，例如速度场 v(x,y,z,t)。而 Lagrangian 方法求解离散的粒子。这篇文章比较清楚地总结了两种方法的优缺点。SPH_2017 对于一片规则的水域，可以采用网格法求解场，但是对于经常需要产生碰撞、形变的水体，粒子方法比网格法有着明显优势。

SPH 全称是 Smoothed Particle Hydrodynamics，Smoothed 是一个重要的特性，这使得每个粒子的特性都可以由临近的粒子计算得出。从而产生了邻域核函数的应用。N-S 方程右侧的两个梯度都可以用核函数的方式进行重写，从而得到简单的求梯度公式。论文中列的很清晰，这里不再赘述。

既然要 in a weekend 那么我们就撸起袖子开干，发挥面向 google 编程的优秀传统。

在 github 上看了几个 sph 相关的代码库。SplishSplash 是工业化做的比较好的，但是不能满足我们 in a weekend 快速阅读实现的需求。所以又找了一个 [SPH-3D-Fluid-Simulation] (https://github.com/saeedmahani/SPH-3D-Fluid-Simulation)。这个项目写的很 easy，课程作业的难度。ok 就是它了，拿来改改。Visual studio，启动！

sph_demo.gif

能运行，效果还凑合，那我们直接进行一个代码的读。

首先是基础流程，所有物理模拟都是类似这个。这里建了一个网格主要用来做邻域的圈定，每个粒子只需要计算邻域粒子的影响。

PARTICLE_SYSTEM::stepVerlet(double dt)
calculateAcceleration()  // 更新所有粒子的加速度

//用加速度更新速度 位置 (Verlet)
VEC3D newPosition = particle.position() + particle.velocity()*dt + particle.acceleration()*dt*dt;
VEC3D newVelocity = (newPosition - particle.position()) / dt;

updateGrid()  //粒子移动后更新所有网格中保存的粒子

求 density , pressure，参考 Matthias 论文式(3) ，（12）

for particle 的所有邻域粒子 neighborParticle:
    VEC3D diffPosition = particle.position() - neighborParticle.position();
    double radiusSquared = diffPosition.dot(diffPosition);
    particle.density() += Wpoly6(radiusSquared);

particle.density() *= PARTICLE_MASS;
particle.pressure() = GAS_STIFFNESS * (particle.density() - REST_DENSITY);

各种核函数及其梯度：

Wpoly6, Wpoly6Gradient, Wpoly6Laplacian : Matthias 论文 式（20） 
WspikyGradient : Matthias 论文 式（21） 
WviscosityLaplacian: Matthias 论文 式（22）

然后计算 N-S 右侧。对于 p0 粒子有：

for p0 的所有邻域粒子 p1:
    f_pressure += p0.pressure/pow(p0.density,2)  +  
                    p1.pressure/pow(p1.density,2) * spikyGradient;
    f_viscosity += (p1.v - p0.v) * WviscosityLaplacian(rsqd) / p1.density;

对于表面张力要特殊处理。计算表面粒子与表面张力如下。

for p0 的所有邻域粒子 p1:
    colorFieldNormal += poly6Gradient / p1.density;
    colorFieldLaplacian += Wpoly6Laplacian(rsqd) / p1.density;

colorFieldNormalMagnitude = colorFieldNormal.magnitude();
if (colorFieldNormalMagnitude > SURFACE_THRESHOLD) {
    particle.flag() = true;
    f_surface = -SURFACE_TENSION * colorFieldNormal / colorFieldNormalMagnitude * colorFieldLaplacian;

汇总计算加速度

f_gravity(0.0, particle.density() * -9.80665, 0.0);
particle.acceleration() = 
(f_pressure + f_viscosity + f_surface + f_gravity) / particle.density();

上面跑完就是一个 SPH 的完整流程了。公式很复杂，但流程其实很简单。

至于渲染，其实粒子可以视为一种点云，而点云在点的数量比较少的情况下，可以通过 marching cubes 方法通过查找等值面来将点云重建为 mesh。然后就可以走正常的渲染了。这也是很经典的方法了，至于有没有其他的特技，就交给更专业的人吧。

SPH in a weekend，跑通收工。

SPH in a weekend

推荐阅读更多精彩内容