Stable Fluid流体模拟

, ,

前一段时间研究了一下流体模拟，还看了一门工程流体力学的公开课。在渲染中的流体模拟还是基于基础的流体力学相关的公式的。刚刚看到keijiro实现了一遍stable fluid，而之前自己实现的不太正确，又重新学习了一遍，做个笔记整理一下。

由于之前是直接看Stable Fluids, Jos Stam实现的，时间较早。而Real-Time Fluid Dynamics for Games省略了一部分繁杂的公式推导会更加简洁易懂。

在Webgl2上实现的Simulation

github: webgl2-stablefluids

LiveDemo:

向量分析

向量分析主要关注向量场的微分和积分性质，研究的对象为矢量场与标量场。向量分析的最重要的微分运算是梯度、散度、旋度、拉普拉斯算子和向量拉普拉斯算子。向量算子符号Nabla算子 ▽\triangledown。

梯度

grad(f)=▽fgrad(f) = \triangledown f
梯度是一个标量场f中某点某点增加率最大的速率与方向，标量场的梯度是一个向量场。

散度

div(F⃗)=▽⋅F⃗div(\vec{F}) = \triangledown \cdot \vec{F}

散度是向量场F中某点发散或汇聚的程度，向量场的散度是标量场。

旋度

curl(F⃗)=▽×F⃗curl(\vec{F}) = \triangledown \times \vec{F}
旋度表示向量场F中某点附近的旋转程序，向量场的旋度还是向量场。

拉普拉斯算子Laplace

△f=▽2f=▽⋅▽f\triangle f = \triangledown^2f = \triangledown \cdot \triangledown f

拉普拉斯预算表示对变量场f作梯度运算后，再对其作散度运算，标量场的拉普拉斯运算还是标量场。

详细可以参看知乎的几个回答

https://www.zhihu.com/question/21912411

https://www.zhihu.com/question/24074028

流体模拟与Navier-Stokes方程

游戏中的流体模拟

在游戏领域有很多流体模拟的方式，主要的分类和流体力学中一致，欧拉方法与拉格朗日法。对应的方法就是SPH(Smoothed-particle hydrodynamics)以及基于Grid的方法，也就是是stable fluid这篇paper提出的方法。

stable fluid这篇paper将NS方程的求解简化到在计算机中实现，同时满足stable的特点。即给定当前流体速度场的状态，以及任意有限的时间间隔t,可以计算出t时间后的流体速度场的精确变化。

对于流体模拟来说我们可以使用particle来模拟流体粒子的运动，从而得到整个流体场的运动状态和变化。但是对于烟雾这类粒度很小的流体来说，我们可以使用密度场来描述。
Real-Time Fluid Dynamics for Games中提到了两个公式，第一个是简化版的Navier-Stokes方程，第二个是密度场的方程。

∂u∂t=−(u⋅▽)u+ν▽2u+F\large \frac{\partial u}{\partial t} = -(u \cdot \triangledown) u + \nu \triangledown^{2} u + F

∂ρ∂t=−(u⋅▽)ρ+k▽2ρ+S\large \frac{\partial \rho }{\partial t} = -(u \cdot \triangledown)\rho + k \triangledown^{2} \rho + S

欧拉法侧重使用场的概念来计算流体的物理属性。对流体分割成一个个grid cell.对每个cell进行物理量的模拟和求解。最终的渲染也是基于网格的渲染。

对于最终的渲染，我们需要计算在 tt 时刻时某个流体单元的密度，只要在 Δt\Delta t 的时间范围内，流体的密度变化是稳定的，就可以渲染出连续的流体运动效果。

接下按照这篇paper的思路来计算流体。

求解NS方程

NS方程是流体力学的基本方程，描述的是流体的物理量随着时间的变化。那么对于流体的速度u、密度ρ\rho都是满足的。所以上文提到的两个方程的形式基本一致。

密度NS方程相对于速度的NS方程来说相对容易，因为密度方程是线性的，我们可以设定速度场是一个固定值，先求解密度NS方程

计算密度场变化

∂ρ∂t=−(u⋅▽)ρ+k▽2ρ+S\large \frac{\partial \rho }{\partial t} = -(u \cdot \triangledown)\rho + k \triangledown^2 \rho + S

密度方程描述了密度的变化量受到三个部分的影响：

密度发散

k▽2ρ k \triangledown^2 \rho 密度随着密度场自身的散度发散。对密度ρ\rho拉普拉斯运算表示密度场梯度的散度。并且扩散以参数kk的程度进行。

1
2
3
4
5
6
7
8
9
10
11
float k = dt * diff * N*N;

//unstable
for(var i=1;i<=N;i++){
    for(var j=1;j<=N;j++){
        q[i,j] = q0[i,j] + k *(q0[i-1,j] + q0[i+1,j] + q0[i,j-1]+ q0[i,j+1] - 4*q0[i,j]);
    }
}

//stable
q0[i,j] = q[i,j] - k *(q[i-1,j] + q[i+1,j] + q[i,j-1]+ q[i,j+1] - 4* q[i,j]);

对于上述方程 我们需要求得q[i,j]。可以看成是一个q[i,j]的线性方程组，我们可以使用高斯-赛德尔迭代或者是雅可比迭代进行求解。

雅可比迭代与高斯赛德尔迭代的区别

Gauss-Seidel迭代k+1次的值用到了第k+1次迭代其他成员的值与k次迭代的值，收敛速度更快。而雅可比迭代，每个成员的计算不需要依赖其他成员只需要第k此迭代的结果，可以进行并行运算，更加适合GPU实现。

速度梯度

−(u⋅▽)ρ-( u \cdot \triangledown)\rho

密度的变化量是随着速度梯度的方向减少的，直观的表述就是速度场将物质搬运到梯度的反方向的位置。

物理量的随速度梯度的变化也称作平流Advection。对于物理量q，我们可以使用以下方程直接计算。

x=x0+u(x0,t)x = x0+ u(x0,t)

q(x0+u(x0,t),t+δt)=q0(x0,t)q(x0 + u(x0,t),t+ \delta t) = q0(x0,t)

然而这个方程,最终的位置取决于速度u，而要使得其稳定则要满足 x=x0+u(x,t)x = x0 + u(x,t) ,而不是 x=x0+u(x0,t)x = x0+ u(x0,t) 这样这个方程更加难解。
采用backtrace的思路得到公式：

q(x,t+δt)=q0(x−u(x,t)δt,t)q(x,t+\delta t) = q0(x-u(x,t)\delta t, t)

对于当前位置x，我们根据x处的速度u逆推出前一步的位置，然后通过该位置的周围四个grid的值进行插值得到密度的平流。

密度源

SS 表示密度的源，可以设定为一个固定初始状态下的密度分布。

最终计算密度场的代码大致如下：

1
2
3
4
5
AddSource() //密度源
while(true){
    Difffuse()  //发散
    Move()      //按照速度梯度移动
}

计算速度场变化

由速度的NS方程可知，速度的变化量也由三个部分组成:

  • 速度梯度的平流
  • 粘滞发散
  • 外力作用

速度梯度的平流

对于速度平流的计算和上文密度的方式相同也是使用backtrace的方式来求解。

粘滞发散

由于流体的粘度会导致流体的速度发散，其发散的速率与流体的粘性系数 vv 有关。

在计算黏性发散时也需要用雅可比迭代计算速度的拉普拉斯运算。

稳定流体与亥姆霍兹分解

不过需要注意的是我们计算稳定流体时要满足NS方程中对不可压缩流体，速度的散度为0，也就是满足方程

▽⋅u=0\triangledown \cdot u = 0

对于经过平流以及发散计算后的速度场，其散度已经不一定为0。只有当速度散度为0，NS方程的迭代才会收敛稳定。

所有我们需要一个project的方法将速度场投影到一个无散的速度场。

亥姆霍兹定理描述了，一个足够光滑的三维矢量场可以分解为一个无旋矢量场与一个螺旋线矢量场的和。

螺旋线矢量场即一个散度为0的矢量场。所以我们需要迭代的速度场为只有旋度的场。同时速度的梯度是一个旋度为0的矢量场。根据亥姆霍兹定理，只需要将计算出的速度场减去速度梯度就能过够得到稳定流体的速度场。

(未完)

References

github StableFluids, keijiro

Stable Fluids, Jos Stam

Real-Time Fluid Dynamics for Games

github 2DFluidSim, candycat1992