【游戏流体力学基础及Unity代码（四）】用欧拉方程模拟无粘性染料之公式推导

IT圈老男孩1 · 发表于 2021-4-30 12:58

先放一张动态图吊一下胃口~下面就是最终的效果
不可压缩的欧拉方程只比NS方程少一个粘性项。所以下面的内容是完全适合NS方程的。各位请准备好！
散度定理

模拟流体的时候会遇到许多数学公式。为了深刻理解这些数学公式，我们先从最简单的复习起，比如说二十以内的加减法。
如下图，ABCD四个区域原本有一些人，此时刻这些人的离开速度为3xi + yj，也就是说位于（Xa = 1，Ya = 1）的A区域的人将有3 * Xa个人从x轴正方向离开A区域进入B区域，将有1 * Ya个人往y轴正方向离开A区域进入C区域。那么经过一个时刻后，一共有多少人离开了ABCD四个小区域组成的大区域？请写出这道题的两种解法，其中第一种方法分别计算四个小区域的离开人数，第二种方法则把四个区域当成一个整体。
如果你正确写出了两种解法，那么赶紧给自己戴一朵小红花吧！如果没有也没关系，我们来看看这两种解法的不同之处。如下图，第一种解法，由十以内的加减法显而易见，ABCD四个小区域人数的净离开量都是4，相加得到总净离开人数是16。第二种解法，我们只计算从大区域边界上离开的人数，毕竟对于小区域之间的人互相移动并不影响大区域的人数离开量。最后结果也就是2+2+6+6=16。
于是我们得出一个一般的结论，对二维空间来说，如果将一个大区域分成很多小区域，那么
从大区域的边界上的总离开量 = 各个小区域的净离开量相加
上面这个公式还有一个听起来高大上的名字，散度定理(divergence theorem)，或者高斯定理(Gauss&#39;s theorem)。所以毕导说很多东西小学二年级我们就学过真的不是瞎吹的。
上面这个公式写成积分形式就是

$\oint_{C} \mathbf{F} \cdot \mathbf{n} \, dl = \iint_{S} \nabla{\mathbf{F}} \, dx \, dy \\$
其中F是我们的向量场，在我们这道题中就是速度。n就是边界的法向量，l是边界的长度，倒三角就是向量场F的梯度。C就是curve曲线，也就是大区域的边界，S是surface表面，也就是整个大区域的面积。至于那个积分符号，在计算机离散的世界里，把它看成加号就行了。
然后三维散度定理，说人话就是，将三维坐标系上的一个三维大空间分成许多三维小空间，那么
从大空间的边界上的总离开量 = 各个小空间的净离开量相加
这个公式的积分形式在各大教科书上长这样，右边那个V是指Volume体积的意思：
散度定理在本系列中主要用于推导动量方程。
速度分解

现在我们来研究固体和流体有什么区别。
比如说一块砖头，第一种运动是平移(Translation)，也就是从一个地方到另一个地方，但不改变朝向。第二种运动是旋转(Rotation)。这两种情况下砖头的形状都不发生改变，是不可压缩的(incompressible)，因此也叫刚体(Rigid Body)。
但对于一只猫来说，情况就不同了。它不仅可以平移旋转，还可以任意改变自己的形状，这也就是为什么有“一杯猫”，“一滩猫”，“一盒猫”的说法。
因此研究人员Marc-Antoine Fardin用流变函数证实了广为流传的说法“猫是液体”，或者更准确的“猫是可压缩流体”，也因此获得了2017年搞笑诺贝尔物理学奖。
因此为了研究猫，或者说是可压缩流体，我们先来看看它能做什么类型的运动。如下图了，除了平移和旋转外。还包括线变形(linear defomation)和角变形(Angular Deformation)。
如果一个力垂直于可压缩物体的表面的话，那么这个力就会造成线变形。比如用力将枕头压扁，如果这个力与表面的切线方向一致，那么这个力就会造成角变形。
比如对鼠标滑轮来说，如果用力向下按，而不让滑轮滑动，那么这个力算垂直于鼠标滑轮的切线方向。如果前后滑动滑轮，那么就算是与表面切线方向一致。
对于由流体粒子组成流体微团来说，各点的平移速度不应该和坐标有关系，也就是一个常数，也就是

$V_{平移} = \{C_{x},C_y\} \\$
而各个点的线变形率应该是一个一次方程，也就是

$V_{线变形速度} = \{\frac{\partial u}{\partial x},\frac{\partial v}{\partial y}\}\\$
对于角变形来说，角变形率用如下的方法求：
如上图，当由橙色正方形变成蓝色平行四边形的时候，取一个非常短的时间，那么点B变成了点B星，它在y轴上的速度为v，在x轴上的速度很小所以忽略不计。点D点在x轴上的速度为u，在y轴上的速度很小所以忽略不计。
所以角度B星AB上的的正切值在单位时间上的变化量应该为y轴上的速度除以AB之间的距离

$tan(\theta) = \frac{v}{L_{AB}} \\$
又因为AB之间的距离很小，夹角B星AB也很小，所以得到下面这个式子。点D同理。

$\theta_B = tan\theta_B =\frac{\partial v}{\partial x}, \theta_D = tan\theta_D =\frac{\partial u}{\partial y} \\$
所以我们要求的第三种速率，也就是角变形速率如下：

$\frac{1}{2}(\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}) = \frac{1}{2}(\theta_B + \theta_D) = \theta \\$
然后对于旋转来说，如下图：
B在y轴上的速度为v，在x轴上的速度很小所以忽略不计。点D点在x轴上的速度为u，在y轴上的速度很小所以忽略不计。由于线速度等于角速度乘以半径，那么B点和D点的角速度：

$w_b = \frac{\partial v}{\partial x},w_d = -\frac{\partial u}{\partial y} \\$
因为B点和D点的角速度实际上是相等的，所以当我们考虑流体微团的x轴速度u和y轴速度v与流体微团整体角速度的关系时，会得到下面这个式子，也就是我们要求的第四种速度，旋转角速度：

$\frac{1}{2}(\frac{\partial v}{\partial x} -\frac{\partial u}{\partial y})= \frac{1}{2}(w_b + w_d) = w \\$
等式最左边括号里的式子，就是旋度(curl)的一部分。符号如下。由上式看出，旋度刚好等于角速度的两倍。如果流体的速度的旋度等于零，那么流体微团的角速度也是零，那么这个流体就是无旋的(irrotational)。比如说如果地球仅仅绕太阳公转，而自身不旋转的话，那么地球的运动就是无旋的。但如果地球不仅公转，还自转的话，那么地球的运动就是有旋的。旋度在后面我们讨论涡流的时候还会遇到。
研究可压缩的“一滩猫”非常麻烦，所以我们先从不可压缩的流体研究起。之前我们说过，我们的不可压缩流体只能有平移和旋转，而不能有线变形和角变形，所以我们把这两种速度分开，也就是
速度 = （平移和旋转） + (线变形和角变形)
写成二维微分形式就是下面这样。

$\begin{pmatrix} \frac{\partial u}{\partial x}&\frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y}&\frac{\partial v}{\partial y}\\ \end{pmatrix} = \Large \begin{pmatrix} 0&-\frac{1}{2}(\frac{ \partial u}{\partial y} - \frac{ \partial v}{\partial x})\\ \frac{1}{2}(\frac{ \partial u}{\partial y} - \frac{ \partial v}{\partial x})&0 \end{pmatrix} + \Large \begin{pmatrix} \frac{\partial u}{\partial x}& \frac{1}{2}(\frac{ \partial u}{\partial y} + \frac{ \partial v}{\partial x})\\ \frac{1}{2}(\frac{ \partial u}{\partial y} + \frac{ \partial v}{\partial x})& \frac{\partial v}{\partial y}\\ \end{pmatrix}$
这个也叫亥姆霍兹速度分解(Helmholtz decomposition)。在张量分析中，这叫把一个张量分解成一个反对称张量和一个对称张量。这个等式的性质可以看看《A Mathematical Introduction to Fluid Mechanics》. 3rd ed. Springer.Chorin, A.J., and J.E. Marsden. 1993.的1.3章。
等式右边第一项就是我们想要的，它的正对角线上两项都是0，代表平移，反对角线上就是它的旋转角速度。它们互为相反数，这是旋转导致的。它也符合无散的性质。
等式右边第二项是我们不想要的，它的正对角线是线变形，而反对角线上的两项是一样的，它是角变形导致的。第二项符合无旋的性质。我们之后会用某个方法把第二项去掉。
我们再回顾一下简单的加减法。例如下图左边16个小区域原本都有有一些人，他们的离开速度是V = {2y,4x}。那么单位时间里人员的流动可以分解可以如下分解：
可以看到中间那个无散度场，如果16个小区域组成的大区域的边界上，进入的人数等于离开的人数，因此是无散度的，但它可能是有旋度的，比如有一个绕圈的人正在沿着EIMNOKGFE走。按照亥姆霍兹速度分解，也就是下面这样 $\begin{pmatrix} \frac{\partial u}{\partial x}&\frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y}&\frac{\partial v}{\partial y}\\ \end{pmatrix} = \begin{pmatrix} 0&4\\ 2&0\\ \end{pmatrix} = \begin{pmatrix} 0&1\\ -1&0 \end{pmatrix} + \begin{pmatrix} 0&3\\ 3&0 \end{pmatrix}\\$ 关于四种速度的分解还有一些不错的资料，比如：
http://users.metu.edu.tr/csert/me305/ME%20305%20Part%206%20Differential%20Formulation%20of%20Fluid%20Flow.pdf
http://book.ucdrs.superlib.net/views/specific/2929/bookDetail.jsp?dxNumber=000006086061&d=F62E39880A1BF4A158F629BA10D59031&fenlei=18210203 的第3.4章，文献传递50页至70页即可。
http://book.ucdrs.superlib.net/views/specific/2929/bookDetail.jsp?dxNumber=000012787629&d=4D3FB1D67093435A9FFB7DD6A7C3926D&fenlei=13020807 《粘性流体力学》的第3.2章，大约是在89页到99页。而且我觉得这本书整体非常不错。
动量定理

牛顿爷爷的第二定律告诉我们，外力是导致一个物体速度改变的原因，也就是

$ma = F \\$ 上面F就是外力，m是物体的质量，a是加速度。但如果这个力持续作用一个时间，那么这个物体就会从这个时间开始时的速度V0增加到时间结束后的速度V1。 $mat = Ft = mV_1 - mV_0 \\$ 这里的mv，就是动量(momentum)吧。动量的概念在高中物理上就有。而流体中动量守恒的概念为

$单位时间内控制体动量的增加量 = \\ 通过控制面流入控制体的动量总和 - 通过流出控制面流出控制体内的动量总和 + 外界作用的力 \tag{1}$
单位时间内控制体动量的增加量通过控制面流入控制体的动量总和通过流出控制面流出控制体内的动量总和外界作用的力
其中控制体(control volume)就是空间中一个固定的地方，比如一个立方体。而控制面(control suface)就是这个立方体的六个面。流体微团可以自由从控制面流进流出。控制体和控制面实际上并不存在，只是为了方便研究而已。
我们先看看把((流入动量) - (流出动量))简化为(-(净流出动量))。我们假设控制面的法向量指向控制面外面。那么流出动量的速度方向就和控制面的方向相同。
如果控制体很小，那么速度的变化就可以忽略，因此动量的变化主要是质量的变化。而质量的变化量等于密度乘以离开控制体的流体体积。密度在这里是常数，而离开控制体的体积等速度乘以控制面的面积。也就是

$m = \rho VS \\$
上式中，m是质量，rho是密度，V是速度，S是表面的面积。那么最终动量的变化是

$mV = \rho V^2S \\$
现在让表面S细分为无穷小的区域，然后将其结果相加，也就是积分起来，就得到了净流出动量，而(-(净流出动量))还要在下面这个式子前加个负号。
然后是流体受到的力。比如说一个水箱破了一个洞，那么这个洞附近的水就会受到外面空气压力的作用而流向外面。但压力只能让在流体表面的水受到影响，所以压力也被称为表面力(Surface Force)。我们现在需要考虑的表面力只有压力。
由于压力的方向由外面指向控制体里面，与控制面的方向相反，所以这个项的前面应该加个负号。然后让压力在表面上积分，得到下面这个式子：
除了表面力，流体整体也会受到力的作用。比如重力，电磁力。这种力作用在所有流体微团上，不管它是不是在控制体的表面，又叫做体积力(body force)。不过，我们这里模拟的流体会忽略掉重力的影响，而加上鼠标拖动的影响。这时候体积力项就是：

$\iiint_{Volume} \rho f \, d(Volume) \tag{4}$
为了将体积与速度的符号区分开，上面这个式子我用Volume来代表体积。
最后，1式等式左边动量可以用下面这个式子表示：

$mV = \rho (Volume)V \\$
单位时间内动量的变化就是

$\iiint_{Volume} \frac{\partial (\rho V)}{\partial t} \, d(Volume) \tag{5}$
集齐2~5这四个式子，就召唤出了流体力学三大基本守恒方程之一——动量方程。
它的意思就和1式差不多，不过得稍微修改一下：
单位时间内控制体动量的增加量 = 通过流出控制面净流出控制体内的动量总和 + 外界作用的表面力 + 外界作用的体积力
不过上面这个式子有个缺点，就是它的积分符号多得吓人。因此我们用散度定理把面积分换成体积分，然后换一下流出动量的位置。于是式子中多了两个梯度符号，也就是倒三角形：

$\iiint_{Volume} \frac{\partial (\rho V)}{\partial t} \, d(Volume) + \iiint_{Volume} \nabla \cdot (\rho V V) \, d(Volume) \\= -\iiint_{Volume}\nabla p \, d(Volume) + \iiint_{Volume} \rho f \, d(Volume)$
然后默念咒语，让体积分这些妖魔鬼怪快点一起离开，那么就简洁多了

$\frac{\partial (\rho V)}{\partial t} +\nabla \cdot (\rho V V) = -\nabla p + \rho f \\$
当然也可以同除以密度，就变成了下面这种更常见的形式：

$\frac{\partial V}{\partial t} +(\nabla \cdot V) V = -\frac{1}{\rho}\nabla p + f \\$
以上的动量定理的推导在各大教科书上都有，很容易查到。但是我目前觉得推导得最好的是《Fundamentals of Aerodynamics》by John D. Anderson, Jr的第2.5节，这本书本身也非常不错。不过好像只有译注版的。
无粘性不可压缩流体的欧拉方程

世界万物都遵循牛顿定律，流体也不例外。不过牛顿第一定律在流体力学中就是连续性方程，对不可压缩流体而言，上篇我们已经推导过了，也就是速度的散度为0。牛顿第二定律则变成了动量方程。联合起来，就得了欧拉方程:

$\nabla \cdot V = 0 \tag{5}$

$\frac{\partial V}{\partial t} + V(\nabla \cdot V) = -\frac{1}{\rho}\nabla p + f\tag{6}$
欧拉方程与NS方程唯一的区别就是，欧拉方程没有计算流体的粘性，而NS方程计算了流体的粘性。也就是
欧拉方程 = 重力 + 压力
NS方程 = 重力 + 压力 + 粘性力
但这里我们把重力省略，而加上了来自鼠标拖动产生的力。这里我们不能直接把5式代入6式，因为6式里的那项实际上是个非线性方程

$(\nabla \cdot V) V = u\frac{\partial V}{\partial x}+v\frac{\partial V}{\partial y} = u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} + u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} \\$
直接求解欧拉方程太过困难，所以我们将其分成两步，第一步，由上一时刻的速度，以及已知的力，在这里就是鼠标的拖动力，求解出中间速度(intermediate velocity)，也就是V*。

$\frac{V^{*} - V^n}{\Delta t} = F^n\\$
第二步，由中间速度，和压力，求解出下一时刻的速度。

$\frac{V^{n+1} -V^{*}}{\Delta t} =-\frac{1}{\rho}\nabla p \\$
这种方法也被称为预测校正步(predictor - corrector)。第一步并不困难，就是本系列第二篇介绍的平流方程。但这时候求解出来的中间速度，它的散度并不是零，因此第二步我们需要一点技巧将其的散度变为零，让下一时刻的流体继续符合不可压缩的性质。
这篇字实在有点太多了，剩下的公式推导和代码就放在下一篇了。
下一篇：光影帽子：【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现
上一篇光影帽子：【游戏流体力学基础及Unity代码（三）】用波动方程模拟三维落雨池塘，连续性方程
宣传：我创建了一个模拟流体交流Q群，欢迎各位喜欢自己编写代码模拟流体的小伙伴们入群互相交流学习！群号：1001290801

kyuskoj · 发表于 2021-4-30 13:06

好像shadertoy上的大佬把这些都玩出花了

NoiseFloor · 发表于 2021-4-30 13:11

请问公式3，压力积分那个，由面积积分换体积积分的时候为什么不用散度而是用的梯度

xiaozongpeng · 发表于 2021-4-30 13:19

那个下三角，对标量场来说就是求梯度，对矢量场就是求散度，不过公式下三角和矢量之间还要加个点。压力场是标量场，所以是梯度。这个可能帮到你为什么空间二阶导（拉普拉斯算子）这么重要？

APSchmidt · 发表于 2021-4-30 13:24

谢谢，受教了

rustum · 发表于 2021-4-30 13:32

m＝rho* V*s那个公式，V是速度，你笔误写成体积了吧~

量子计算9 · 发表于 2021-4-30 13:39

已修正。谢啦

Zephus · 发表于 2021-4-30 13:44

讲的真棒[赞同]

fwalker · 发表于 2021-4-30 13:49

流体力学老师没你讲的明白Orz 大神

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（四）】用欧拉方程模拟无粘性染料之公式推导

本帖子中包含更多资源

浏览过的版块