一、前言

游戏中的场景交互一直是人们追求的，比如人物与水的交互、物体与烟雾交互、可交互雪地等。它们都有一个共性，就是这些物体都是具有流体效果的（雪地可以看做是流体粘度无限大）。因此我想尝试一下流体的模拟是如何做出来的。这个难度是远高于预期的，其中充斥着无数的数学方程。与此同时，将这些数学公式离散化转换为分步求解的代码也是一大挑战。

我看了网上很多的文献以及许多大神的做法，但大神们普遍言简意赅，而有些在求解过程中出现了不必要的求解方式。我则在这些基础上，根据自己的理解整理了一份模拟思路。

然后我简要介绍一下NS方程，这个方程是流体模拟的关键。NS方程全程是纳维-斯托克斯偏微分方程组（Naiver-Stokes），使用这个方程基本可以描述各种流体在各种情况下的运动方式。求解三维空间中的N-S方程组光滑解的存在性问题是千禧年十大问题之一，它的难度可想而知。下面展示的是不可压缩流体的NS方程

 

在前言的最后，我们要明确以下几个问题，首先NS方程的求解是舍弃了很多项，并且简化很多条件后的数值解；我们始终是在2D维度求解NS方程；我所解释的更为具体化，尽管我已经尽可能的求证，但也可能存在错误的解释；下文涉及到的公式我都尽可能以普通高数水平和较高高中数学水平进行解释。一定要注意区分方程中矢量场和标量场。

 

二、NS方程推导

矢量微积分算子

不可压缩流体的NS方程中，很多量都是矢量，因此对于矢量的计算，我们需要引入一些算子，用于简化表示矢量的计算。下述的都是2维的算子公式。

梯度算子 ▽ ：，物理含义是矢量场在各个方向上的偏导的矩阵。更通俗来说，就是一个向量场上的某个点，在此位置的某个方向上的变化率，以速度场求梯度▽u，就可以理解为某点在x方向的加速度和在y方向的加速度所组成的向量。

散度算子 ▽· ：，物理含义是矢量场在此点偏导的求和，得到的结果是一个标量。也就是说，散度是当前点的物质增加还是减少。

旋度算子 ▽× ：，物理含义是矢量场在此点偏导的差，得到的也是一个标量。通俗来说就是此点的环流情况，也就是是否会产生涡旋现象。

拉普拉斯算子 ∇² = ▽ · ▽：，拉普拉斯算子不太好进行直观的解释，大致来说它提供了标量场的曲率和平滑度的信息。

 

重要的一点，速度场 = 速度的散度场 + 速度的旋度场

 

 

基础物理公式

首先先列出几个最基础的物理公式：

动量方程：p = m*v

动量方程另一形式：F = p / t

牛二定理：F = ma

密度公式：ρ = m/V

这些基础公式有了，再非常快速的介绍一下欧拉方法和拉格朗日方法，两个方法都是假设空间中有一个微元i，欧拉法相对参考系是整个世界，拉格朗日法相对参考系是以微元自身。我们将使用欧拉方法进行NS方程的推导，下面开导：

动量守恒

现在假设一个方程，u = u(t, X)，它用于描述一个粒子在时间 t 和位置X时的速度方程。方程u对时间求导，物理意义为一个二维平面上的点元在t时的加速度。这个式子被称为拉格朗日变量。注意这里是点乘而不是乘积。

 

x是二维的，∂X/ ∂t 就是该点的速度 u，然后我们将∂u / ∂X = (∂u / ∂x , ∂u / ∂y) 记作▽u，得到：

         （1）

 

再换一个角度，考虑一个质点受到的力，F总 =  接触力+黏性力（剪应力），这里是根据牛顿流体进行推导的。（牛顿流体（英语：Newtonian fluid）指应力与应变率成正比的流体。此比例系数为流体的黏度。）

 

就可以写成 F总 = F(压强) + τ + F(额外力)，我们分别求出F(压强) 和 τ

F(压强) =  VdP / dX = (∂p / ∂x, ∂p / ∂y)。我们将∂P/ ∂x记作▽p，得到：

F(压强) = V▽p                 （2）

然后请看牛顿流体特性的最后一句话，因为这里的速度场是矢量场，因此需要先对速度场求梯度，得到方向速率变化，得到的值点乘梯度才能得到速度场在各方向的梯度，得到：

τ = μ*▽·(▽u)= μ*∇²u       （3）

下面，我们综合（1） = （2）+ （3）+F(额外力)

m(∂u / ∂t + u·▽u) = V▽p +  μ*∇²u + f ， 两边同时除m 并移项

=>  ∂u / ∂t = -u·▽u - 1/ρ* ▽p + v*∇²u + F   (常数项v = μ/m， F = f/m)

 

质量守恒：

m0 - m1 = ρ*△V = u*△t*ρ =ρ* ▽ · u = 0

=> ▽ · u = 0  

 

到此为止，令人惊叹的，我们推导出来了不可压缩流体的NS方程！

 

三、数学方法求解方程

现在我们假设一个情况，我们有一个当前时刻的速度场，那么在接下来的一段时间内，这个场将如何变化呢？也就是说，我们已知u，如何去求解上面的NS方程，得到加速度，进而得到下一刻的速度场。以u = at为例，这是连续的方程，由于是在电脑中模拟，我们还需要离散化。

很明显这是一个欧拉视角的求解。

亥姆霍兹-霍奇分解

用公式表示为

这个在讲解算子的时候说到了，任何向量场都可以分解为两个其他向量场的和，

速度场 = 速度的散度场 + 速度的旋度场。

回看一下动量守恒项，我们简化一下NS方程，假设这个流体是无粘度的理想液体，公式的黏性项v*∇²u 就可以舍去了：

 

也就是说力只剩下了F = -1/ ρ ▽p

又有：F = ma = m * (u(n+1)-un) / △t,带入得：

      （1）

又因为（之前说到的质量守恒公式）：

（2）

 

现在令△x = △y = △t = 1，将（1）带入（2）得：

 

 

这个式子被称为压力泊松方程(Pressure Possion Equation)

−Δt/ρ∇⋅∇p=−Δ⋅⃗u

 

四、转换为代码

上面主要是说了如何将▽p转换为关于u的离散化方程，速度场也是要迭代的，这个迭代就是平流项。那么大致来说呢，先初始化一下速度场，求个平流项（平流项是速度场点乘速度场的梯度），然后加一个压力场（压力场是通过速度场迭代求解出来的），最后再保证一下速度场的散度为0。就这几步就ok了！

1、初始化一个速度场和密度场（也可以叫染料场）。

两个场的的区别只有颜色是不同的。下图是我自己提前画了一个速度场存储在了RenderTexture中，（这一篇我希望更偏学术向，我不再介绍什么是RT，什么是shader等等了。不懂的可以简单理解为RT就是图片，shader就是如何计算像素（体素）的。）为什么颜色是这样的呢，画图解释：颜色是rgb，我们只取rg通道作为速度向量的xy，黑色是由于负数，由于0也会显示黑色，造成了颜色变成了这样。

 

模拟后可以发现确实如我们所想的，鼠标初始化速度场也就是这么做的。密度场相当于对有速度场的区域上色，比如空气流动是看不到的，但我们喷一些烟上去，也就看到了整个速度场的形状。换成鼠标上色的方式。

2、计算平流项（对流项）。

平流项也就是计算：-（u·▽）u = -（▽u）·u，速度场刚刚初始化了，对它求梯度再相乘就好了，这里能转化是因为梯度算子满足分配律、分解率，但不满足交换律。那么对速度场求梯度，我们只需要代入之前的离散化方程就可以了。

看代码，我们在shader中对每一个像素点进行平流，简化代码只剩下平流项，画圆圈可以发现颜料都好像是沿着切线方向移动了。看shader代码，先来一个简化版本的，看公式这不就是求平均值嘛！。采样速度场，得到的是当前点的速度向量，x = vt, 采用前向欧拉法，i.uv=> i.uv-△x(i - (i-1))。运行一下，没问题！But！注意步长，我们选的很小，如果稍微大一点呢？

 

爆炸原因是：零空间问题，有兴趣可以自己搜索一下，简单解释就是这种方法是有偏的，并且精度是O(△x)。短步长会让离散解较为接近，增加步长后由于没有考虑中间时刻的变化，导致数值解出现高频振荡。反映在屏幕上，远处还好一点，鼠标位置爆炸最严重。这就是在鼠标地方速度有一个突变。

看来需要一种好一些的方法来计算步长。这里我们引入一种网格结构：MAC网格，速度分量记录在网格边界上，然后我们采用一种半拉格朗日法（拉格朗日视角解决欧拉视角的问题），保证了函数解最终会收敛。GPUGems图中展示了一个向量如何反向求解，得到的是质点在△t后的新的速度。为什么要双线性差值，向量很大概率不指向uv中心，而双线差实际上就是对周围取均值，精度也提高到了O(△x^2)。并且好就好在，tex2D函数，或者UE中的Sample2D采样Texture都是先二线差再返回颜色，因此从代码角度来说简化了我们的编程。

 

coord对应实际应该偏移的距离（目标位置的uv坐标），将这个点对应到MAC网格上，求MAC网格附近线性插值后的颜色。

 

3、求压力项▽p

亥姆霍兹-霍奇分解，得到的公式，这个公式很明显是一个迭代的形式，这个就被称为雅可比（Jacobi）迭代。

这一项不就是▽·u嘛！因此我们先求出速度场的散度，再代入散度值进行迭代。乘0.9999是浮点数精度问题。

 

还有一种迭代法叫高斯-赛德尔迭代法，也很不错，而且更为易懂，我就不详解了。

4、组合起来

上面的步骤我们得到了迭代后的压力场，还有速度场，我们再应用一下NS方程的第二项就是那个不可压缩项（u无散度 = u-p），最终得到的是无散度速度场，将这个场应用到密度场上，就可以了。

举一反三，如果把速度扩散项设置为比较小的值，那么是不是可以做一个雪地效果，或者是书法笔迹？如果在最后不对密度场平流，而对压力场平流，好像水波效果有了哦！（可以看一下games103 P10课程，思路是一样的。）再想想风场呢……这些只是一个大概思路，有兴趣的可以继续探索哦！

5、扩展计算    粘度项（涡度）

把他放在最后是因为其实目前的模拟效果已经可以接受了，而从原理上来讲，我们之前计算压力项时，实际是舍弃了这一项，现在我们将这一项加回来。

还记得这个公式吧，剪应力的表达式，仔细看看，你会发现为什么这么巧，这不就是对速度场求旋度嘛

我们令 ：

w = ▽×u

则 f = ε (N×w)， 其中N = ▽|w| / （| ▽|w| |），w就表示涡量，对涡量归一化再乘w的大小，实质就是让涡量方向统一。求得的f = ma = mvt，m = 1，则△v = f*△t。让速度场叠加上去就可以了。

 

五、总结

用物理推导公式 -》从公式中可以看出来我们要求的是平流项、压力项、扩散项、质量守恒项 -》离散化后发现很多项和提到的算子息息相关 -》 离散化求解平流项用半拉格朗日法求 -》 压力项先求散度再用雅各比迭代求解 -》 粘度项先求旋度再归一化求解 -》 得到的结果还需要满足质量守恒 -》最后将速度场应用到密度场就行了。

 

六、参考文献

https://yangwc.com/2019/05/01/fluidSimulation/

https://zhuanlan.zhihu.com/p/374280832

https://zhuanlan.zhihu.com/p/270531017 https://developer.nvidia.com/gpugems/gpugems/part-vi-beyond-triangles/chapter-38-fast-fluid-dynamics-simulation-gpu https://www.bilibili.com/video/BV1i54y1F73W/?spm_id_from=333.788&vd_source=e883ef8d644cde54ca3470fb7bdc548d

 

源码地址： https://github.com/clatterrr/FluidSimulationTutorialsUnity

 

b站 邱丑丑啊