Stable Fluids 笔记 （2022更新）
First Added: Dec, 2009. Updated: Jan, 2021. Oct, 2022
Stable Fluids and Real-Time Fluid Dynamics for Games (Jos Stam 1999)
这两篇文章是我在本科期间开始研究图形学看的第一篇文章，前前后后看了很多遍，每一次都有新的感悟。
1、简介
首先Stable Fluids的基础是Navier-Stokes方程组（下称NS方程组）。这个方程组在流体动力系统的计算中举足轻重，可以说是作为流体模拟的首选。而Jos Stam的Stable Fluids一文在近些年的研究中具有里程碑式意义，很多人都是以此为基础进行其研究。
（1）（2）的NS方程组的形式是众多表达形式中的一种.这里的流体是不可压缩(incompressible)流体，也就是密度、温度不变。同时,质量守恒(1式)和动量守恒(2式)需要严格保证.
（1）式中计算散度(Divergence), 速度场u的散度为0。具体到场内某一点，各个方向变化量之和为零。(不会积蓄和释放能量的情况).
（2）式中各项的物理意义,P是压力场(标量场)，希腊字母p表示密度，希腊字母v表示粘滞系数(viscosity)，最后一项f是外力。 对（2）式的把握是个关键，现在我就对等式右边逐项讨论：第一项叫对流项（advection），速度会沿着自身的方向运动。可以看出，速度快的地方就会越来越快，漩涡阿、龙卷风就是这个原理。因为u出现了两次，所以该项非线性。第二项的意思是，速度受压强梯度场变化，压强大速度小，反之速度快压强小(苏联驶出站台的火车把月台的送客的人吸进去这道理?)。第三项为扩散项,它告诉我们，速度在自身二阶导下扩散。扩散越强烈，场就越平滑。最后一项最容易，是外力项。
（2）式左边是加速度（速度关于时间求导），其实这个式子是符合基本形式f=ma的。
现在看上面这个式子，形式上和（2）式有很多相似的地方。事实上，（2）式可以用来建模很多东西。在Real-Time Fluid Dynamics for games里,Jos就用了密度来产生效果。
2、离散化 Discretization
计算机模拟是要把公式进行离散化，加以编程技巧和数值方法.
i)对导数的离散化理解。二维三维中，对于各种差分（forward,back ,central）的实施，我们一般会用网格铺分，一般而言，我们是均匀的铺分成小格子，体素，欧拉栅格, 也就是说dx=dy=dz.这样在梯度、散度、拉普拉斯算子什么的我们分母就统一了，减少麻烦。
ii)边界条件（Boundary Condition, 下称BC）。Jos 此文提到了两个，Perodic BC和Fixed BC. 简而言之，Perodic BC就是说空间中重复着无限个我们的主模拟空间central box, 从central box 左边出去的粒子，会从central box 右边进入。Fixed BC就是说，边界值取一个固定的值。为什么会有边界问题啊？因为我们处理导数的离散问题时候，比如第（3，8）号点用到用到了（3，9）（3，7）（2，8）（2，9）四个点，但是边界（0，1）就没有四邻域了，（0，0）更没有了，我们需要给他们个交代。（以上描述的仅仅是最最简单的边界情况）
iii)合适的dt。很多前人的算法比不上Jos Stam就是dt使模拟不稳定，数据跳跃。这里指出一个经常用到的条件，Courant-Friedrichs-Lewy Condition，下称CFL。CFL的说白了就是，dt不能太大，dt*u<dx,只有这样才能保证信息能够随着栅格传播。实施中dt的取值有多种方法。
iv)数值方法。因为离散计算不得不考虑精度，不同阶的精度在时间累积和错误累积下的效果是迥异的。Upwind Differencing, ENO, WENO, TVD Runge-Kutta。膜拜龙格库塔、切比雪夫、雅克比、欧拉这些大牛。
3、Implementation
根据Helmholtz-Hodge Decomposition，任意矢量场可以被唯一的分解为一个散度为0的场和一个梯度场的叠加。我们要求解的正是散度为零的，所以我们要从总场中Cast掉那个梯度场,投射操作P. 上面是我的一页幻灯片。Jos的求解速度场就是按照（5）式的思路，一步步进行（i）加力、（ii）扩散、(iii)对流、(iv)投射。
i) 加力/源：对每个单元 x+=f*dt
ii）扩散。也就是 拉普拉斯算子(u)， 每个单元和周围的邻域交换。
x0: previous state, x: current state, a: diffuse rate
num : number of its neighbors
sum_p: the sum of its neighbor (previous moment)
sum_c: the sum of its neighbor (current moment)
一个简单的实施
For Each Cell:
x = x0 +a*(sum_p – num* x0 )*dt
(第二项可能为负，产生摇摆和发散，所以不这么做)
Jos’ implement: trace BACKWARD in time
For every cell:
x0 = x - a*(sum_c – num*x)*dt
x0= x – a*sum_c*dt+a*num*x*dt
(1+a*num*dt)x = x0+a*sum_c*dt
x = (x0+a*sum_c*dt)/(1+a*num*dt)
上式的分母衡大于1，所以呢比较稳定
Jos 认为对这个的求解没必要进行求逆（matrix inversion routine）他用了一个比较实用的方法就是Gauss-Seidel relaxation, 迭带若干次，松散化一下。我的体验就是，其实没必要带很多次，稍微几次效果可以很不错了：）
（iii）对流
这里要求速度沿自身方向运动了，前面说了该项已经不是线性的了，这会用上了半拉个朗日方法。我们把每个网格的中心建模为一个粒子，trace BACKWARD in time的思想，沿着速度场回溯，看dt时刻之前位于哪个位置的粒子在会运动到当前的位置（s=so+v*dt, so=s-v*dt），用线性插值的方法求出dt时间前那个粒子，再把该粒子的属性值赋给当前粒子。
（iv）投射
对于矢量场，比如速度场，
第一步，求“负”的散度场，这里负号非常关键，因为会用到和扩散项同样的solver，而其实在解泊松方程，可以少一次乘以-1的操作，
-div = -0.5*(ui+1 – ui-1+ vi+1 – vi-1 + wi+1 – wi-1)/dx
第二步，高斯赛德尔迭代法，解泊松方程 ——也就是什么场p的二阶导会成造成这个散度场（标量场）？这个和之前的扩散项用了同样的solver。
第三步，解出了p，我们再减去p的梯度场，就满足了所谓投射，也是确保了稳定性。这一步是本方法的关键！
u -=(p[i+1,,]-p[i-1,,])/2dx
v -=(p[,j+1,]-p[,j-1,])/2dx
w -=(p[,,k+1]-p[,,k-1])/2dx
4、总结
要从根本层理解求解思路，或者发明新的方法，数学还是最强武器！实践种，离散化的思想很重要，这些是《数据结构》和《算法》之外计算机科学系没有系统交的，《数值分析》是非常好的课程，不过是门选修课。基于物理的渲染流体，说白是个计算物理的（近似）过程，数学是重头戏，图形是个输出接口。
