【游戏流体力学基础及Unity代码（三）】用波动方程模拟三维落雨池塘，连续性方程

skyangel · 发表于 2020-11-29 11:23

散度

最近要进行全国第七次人口普查啦，为了全社会的繁荣稳定，我们决定研究一下和水流很相似的人流的性质。假如有一个无限长的大厅，我们将它划分为很多方形区域。每个区域都有一些人，人可能会到处移动：如下图
现在我们想知道在某个区域的人，离开这个区域的程度，怎么算呢？比如区域E，原来有一些人，现在单位时间内有u个人从左右边界离开，有v个人从上下边界离开，现在区域E少了多少人？
答案很简单，就是u+v，我们可以叫它&#34;某个东西离开这个区域的程度&#34;或”分散程度“，或者简称叫做“散度”，英文Divergence。

$div V=\nabla \cdot V = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \\$
V就是一个速度向量，包括x轴的速度u和y轴的速度的v。倒三角即为散度符号。假如进入这个区域的人数和离开这个区域的人数相等，那么这个区域的散度就为零。如果离开的速度大于进入的速度，那么散度就为正，比如源（Source）就是这样的。如果进入的速度大于离开的速度，那么散度就为负，比如汇（Sink）就是这样的。
比如对于上面的大厅，有大厅入口的区域A就是源，因为并没有别的区域的人进入区域A，而一直都是区域A的人离开区域A而进入别的区域。有大厅出口的区域I就是汇。对于一个雨中池塘来说，雨滴落到水面之处就是池塘的源。如果有个抽水管正将水抽离池塘，那么抽水管之处就是池塘的汇。
我们再来看看隔壁电磁场是怎么玩的。电场和流场和一样也有源和汇，不过电场的源是正电荷，汇是负电荷。
电场强度的散度和什么有关呢？当然就是正电荷和负电荷的分布情况啦！

$\nabla \cdot E = \frac{\rho}{\varepsilon_0} \tag{1}$
上面E就是电场强度，p就是单位体积电荷密度。
而磁场有一点不一样。虽然我们可以在电场中找出正电荷和负电荷，但却找不到磁场单独的源和汇——磁单极子。我们现在只能观察到成对出现的磁场正负极，而没法将它们分开。既然没有单独的源和汇，那么磁场中的散度就到处都是0啦。

$\nabla \cdot B = 0 \tag{2}$
公式1和2就是麦克斯韦方程组中的两个方程。是不是看起来有点不对称？巧了，不仅是你和我，很多科学家也这么想，他们想找出磁单极子，来稍微完善一下大一统理论。而且很多理论研究借助磁单极子的思想将更加方便。但是，无奈现在就是找不到磁单极子。但是无论是证明了磁单极子的存在还是不存在，都将影响很多科学研究和理论。
或许你也可以试着找找看？磁单极子是否存在的答案，也许就隐藏在某本书，某篇论文或某个公式之中。
到那时候就可以批量生成可以产生磁场的皮卡丘或御坂了。【误】
欧拉形式连续性方程

我们再来看一个简单的问题。在之前提到的那些区域中，假设区域E一开始有8个人，现在这个时刻区域D有3个人离开了区域D，穿过了区域D与区域E的交界处，来到了区域E。那么下一个时刻区域E有多少人？
“太简单了”，你一边想一边算出了结果“8 + 3 = 11”。
就在你用0.00001秒心算出结果的时候，你不仅严格遵循了牛顿第一定律“物质不会凭空产生，也不会凭空消灭”，还顺便将其扩展到了流体力学，推导出了流体力学三大基础守恒方程之一——连续性方程。
我们再来分析一遍你思考的过程，对区域E来说：
区域E增长的人数  = （从区域DE交界处进来的人数  + 从区域BE交界处进来的人数 + 从区域FE交界处进来的人数 + 从区域HE交界处进来的人数）
这么长一串文字写起来太麻烦了，我们化简一下。首先：
从区域DE交界处进入区域E的人数 = -从区域DE交界离开区域E的人数
注意等式右边有个负号，也就是如果进来3个人，相当于离开了（-3）个人。并且DE和FE都是x轴方向，剩下两个都是y轴方向。那么对区域E来说：
区域E增长的人数  + （从x轴方向离开的人数  + 从y轴方向离开的人数）=  0
然后把时间也考虑上：
区域E在单位时间增长的人数  + （区域E的人数乘以从x轴方向离开的速度  + 区域E的人数乘以从y轴方向离开的速度）=  0
写成微分形式就是

$\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} = \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho V) = 0 \tag{3}$
p就是人口密度，t时间，u就是x方向上的速度，v是y方向上的速度，倒三角形就是散度符号，
拉格朗日形式连续性方程

在刚才我们提到的大厅中，假如有一个旅行团（流体微团），由一些旅客和导游（流体粒子）组成，他们会尽量聚集在一起，但是整个旅行团不同的时间会前往不同的地方。
这个旅行团成员会占据一定区域。这个区域的面积都会随着时间和空间变化。因此这个区域的人口密度也会变化，这就是物质导数。

$\frac{D \rho}{D t} = \frac{\partial \rho}{\partial t} + u\frac{\partial \rho}{\partial x} + v\frac{\partial \rho}{\partial y} = \frac{\partial \rho}{\partial t} + V(\nabla \cdot \rho) \\$
旅行团的人数遵循牛顿第一定律，不会无缘无故减少。所以旅行团占据区域的人口密度变化的原因就是成员在空间上分散了，如下图，原来四名成员占据四个方格，现在占据六格，导致人口密度变低。
或者说，人口密度之所以变低，是因为他们有一个互相离开的速度，写成数学公式就是如下，也就是拉格朗日形式的连续性方程：

$\frac{D \rho}{D t} + \rho \nabla \cdot V = 0\tag{4}$
两种形式的联系

可以看到我们在推导欧拉形式的连续性方程时，我们关注的是某个区域，这个区域位置和大小都不会变化，流体微团可以自由进出。而推导拉格朗日形式的连续性的连续性方程时，关注的是某个封闭系统，其位置和大小都会变化，但是流体微团不能自由进出。这两种方法各有优劣。虽然现在我们使用欧拉方法，但前者使用也很多。比如涡方法(vortex method)和光滑粒子流体动力学(smoothed particle hydrodynamics)使用的是拉格朗日方法，
这两种形式的连续性方程是如何联系到一起的呢？根据《高等数学》中导数概念有：

$\frac{\partial (\rho u)}{\partial x} =\rho \frac{\partial ( u)}{\partial x} + u \frac{\partial (\rho)}{\partial x}\\$
因此3式中的散度符号可以换成

$\nabla \cdot (\rho V) = V(\nabla \rho) + \rho(\nabla \cdot V) \\$
代入3式得，再和4式比较：

$\frac{\partial \rho}{\partial t} + V(\nabla \cdot \rho) + \rho(\nabla \cdot V) = \frac{D \rho}{D t} + \rho \nabla \cdot V = 0 \\$
不可压缩流体中的不可压缩是指，流体微团的密度不会改变。如下图：
因此得到如下结果：

$\frac{D \rho}{D t} = \nabla \cdot V = 0 \\$
第二项写成微分形式就是下面这样，我们以后还会经常见到这个公式的：

$\nabla \cdot V = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \\$
好了，我们旅游团的成员现在只会跟随旅游团一起移动，而不会独自到处乱跑啦，这能让导游的工作量大大减轻。把流体看成不可压缩的，也能让计算大大简化。
波动方程

前面扯了一堆连续性方程，其实和波动方程关系并不大，完全只是为了减少下一篇的篇幅~
回顾第二章的一维对流方程，用这个方程去模拟波的话，这个波只能向一个方向传播。

$\frac{\partial \rho}{\partial t} + u\frac{\partial \rho}{\partial x} = 0\tag{5}$ &#39;如果我们想要它往反方向传播，就把u改成负数：

$\frac{\partial \rho}{\partial t} - u\frac{\partial \rho}{\partial x} = 0\tag{6}$
但上式仍然只能让波往一个方向传播。如果想让它同时向正反两个方向传播怎么办？很简单，将5式和6式乘起来，就得到了一维的波动方程，它能让波同时向两个方向传播。思考一下为什么是相乘而不是相加。

$(\frac{\partial \rho}{\partial t} + u\frac{\partial \rho}{\partial x})(\frac{\partial \rho}{\partial t} - u\frac{\partial \rho}{\partial x}) = \frac{\partial^2 \rho}{\partial t^2} - u^2\frac{\partial^2 \rho}{\partial x^2} = 0 \\$
这样一来，二维波动方程也很容易得出：

$\frac{\partial^2 \rho}{\partial t^2} = u^2\frac{\partial^2 \rho}{\partial x^2} + v^2\frac{\partial^2 \rho}{\partial y^2} + f\\$
一维的波的”所有方向“只包括两个方向，而二维的”所有方向“组合起来成了一个圆。三维当然就是一个球啦。
注意等式右边还有一个f，这个代表的是源，与初始条件不同的是，源每时每秒都在发挥作用。离散化的公式直接截图吧：
解得
于是得到2D波动方程的着色器主要代码如下：

float2 uv = float2(uvx, uvy);
float TimeMinus2Center = tex2D(_MainTex, uv).g;
float centerT = tex2D(_TimeTexM1, uv).g;
float FSource = 0.0f;
if (flag == 2)
{
FSource = sin(_Time.y * 10.0f) * 50.0f;
}
uv = float2(uvxPlus, uvy);
float rightT = tex2D(_TimeTexM1, uv).g;
uv = float2(uvxMinus, uvy);
float leftT = tex2D(_TimeTexM1, uv).g;
uv = float2(uvx, uvyPlus);
float topT = tex2D(_TimeTexM1, uv).g;
uv = float2(uvx, uvyMinus);
float BottomT = tex2D(_TimeTexM1, uv).g;
_Speed = 0.6f;
_Height = 2 * centerT - TimeMinus2Center + (leftT + rightT + BottomT + topT - 4 * centerT) * _Speed * _Speed + FSource;

复制代码

完整代码请看
效果如下。不过这阴间配色看起来不像是有雨滴落在池塘里，而是有人试图在岩浆里游泳。
嗯，整天玩全屏着色器玩腻了，我们试试三维的吧，画三角形时间到！虽然说是三维的，但实际上还是用二维波动方程，不过是把颜色的变化改成了三角形高度的变化。wave3d.cs中主要代码如下：

for (int i = 0; i < Dimension; i++)
{
for (int j = 0; j < Dimension; j++)
{
int VerIdx = j * Dimension + i;
if (i > 0) LeftHeight = verts[VerIdx - 1].y;
<span class="p">(i < Dimension - 1) RightHeight = verts[VerIdx + 1].y;
if (j > 0) BottomHeight = verts[VerIdx - Dimension].y;
if (j < Dimension - 1) TopHeight = verts[VerIdx + Dimension].y;
CenterHeight = verts[VerIdx].y;
float Fsource = 0.0f;
if (i == Dimension / 2 && j == Dimension / 2) Fsource = Mathf.Sin(Time.time * 2.0f) * 2;
TimeMius2CenterHeight = TimeMius2Verts[VerIdx].y;
float _Speed = 0.6f;
float _Height = 2 * CenterHeight - TimeMius2CenterHeight + (LeftHeight + RightHeight + BottomHeight + TopHeight - 4 * CenterHeight) * _Speed * _Speed + Fsource;
NewVerts[j * Dimension + i] = new Vector3(i / 10.0f, _Height, j / 10.0f);
}
}

复制代码

完整代码请看：
最终效果如下：
虽然能用了，但是效果非常之丑。所以，多试一下不同的源，也就是FSource，以及不同的传播速度，换上不同的材质，用波动方程把场景打造成一个真正的意境唯美的三维落雨池塘，就是本篇的挑战啦。
此时我们离NS方程已经很接近了。不过在此之前，我们还需要一点开胃小菜，也就是下一篇光影帽子：【游戏流体力学基础及Unity代码（四）】用欧拉方程模拟无粘性染料之公式推导
上一篇：光影帽子：【游戏流体力学基础及Unity代码（二）】用平流方程模拟染料流动

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

[笔记] 【游戏流体力学基础及Unity代码（三）】用波动方程模拟三维落雨池塘，连续性方程

本帖子中包含更多资源