【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak

术数古籍专卖疤 · 发表于 2020-12-5 20:07

开局一张图！
浅水波方程
当水波的波长远远大于水的深度时，水波就表现出来一些特别的性质，比如江河就是典型的浅水，可以用浅水波方程(Shallow Water Equation)来模拟，这个方程又叫Saint Venant Equation。
要推导一维的浅水波方程，同样要从连续性方程和动力方程入手。对于X轴上某个坐标点来说，这点此时刻水深度和上时刻水深度不同，这是因为与附近水的深度差以及水的速度。这里假设速度在z方向上没有变化，或者说速度在z方向上的变化非常小以至于可以被忽略。于是就写出了浅水波方程中的连续性方程。这里的速度和水深有关，所以也是一个非线性方程。

$\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} = 0\tag{1}$
然后复习一下初中物理。水中任何一点的压力都是向四面八方的。将一块方形石头放入静水中，它的前后两面所受的压力可以抵消，左右两面也是如此，但是上下两面的压力因为重力的原因而无法抵消。这就是静水压力(hydrostatic pressure)。因此进入深水区需要特殊的潜水装备甚至是潜水艇，否则人类身体无法承受如此大的压力差。下面是简单的推导。

$p = \frac{F}{S} = \frac{mg}{S} = \frac{\rho V g}{S} = \frac{\rho S h g}{S} = \rho g h\\$ p是压力，F是受到的力，S是面积，m是质量，g是重力加速度，V是体积，h是水的深度。上式说明密度越大或水越深，那么受到的压力越大。在我们的实验中，密度和重力加速度是常数，只有水的深度会随着x轴或y轴变化。也就是说：

$\frac{\partial p}{\partial x} = \rho g \frac{\partial h}{\partial x}\\$ 然后开始写动量方程。由于水的粘度很小，这里继续忽略掉粘性项。改写一下第四章提到的Euler方程

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = -g\frac{\partial h}{\partial x}\\$
然后上式两边同乘以水深，就得到了浅水波方程中的动量方程

$\frac{\partial}{\partial t}(hu) + \frac{\partial }{\partial x}(hu^2 + \frac{1}{2}gh^2) = 0\tag{2}$
1,2式有一些相同的地方，我们可以把它们写成下面这个式子，再加上坡度和摩擦力的影响，就得到了一维浅水波方程(Shallow Water Equation)。

$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = S\\ U = \begin{bmatrix} h \\ hu \end{bmatrix}, F = \begin{bmatrix} hu \\ h^2 u + gh^2/2 \\\end{bmatrix},S = \begin{bmatrix} 0 \\ gh(S_0 - S_f) \end{bmatrix}$
等式右侧S代表浅水底坡度和摩擦力造成的影响。其中S0是坡度，Sf是摩擦力，其中\eta是Manning resistance coefficient取决于不同的摩擦规则与情况。

$S_0 = \frac{\partial z}{\partial x},S_f = \frac{u^2\eta^2}{gh}\\$
浅水波方程和Burgers方程很相似，其求解方法就是上一篇提到的各种方法。部分方法的python代码我也放到了与上一篇相同的仓库里。clatterrr/CFDcodepython。要注意的就是不要把高度初始值设置为零，否则到后面就要除以零了。而且浅水波方程使用变化的时间步dt可以获得更合理的结果，以及可以使用不同的边界条件。本篇稍后讨论。二维浅水波方程如下：

$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = S\\ U = \begin{bmatrix} h \\ hu \\ hv \end{bmatrix}, F = \begin{bmatrix} hu \\ h^2 u + gh^2/2 \\ huv \end{bmatrix},G = \begin{bmatrix} hv \\ huv \\ h^2 v + gh^2/2 \end{bmatrix},S = \begin{bmatrix} 0 \\ gh(S_{0x} - S_{fx}) \\ gh(S_{0y} - S_{fx}) \end{bmatrix}\\ S_{0x} = \frac{\partial z}{\partial x},S_{fx} = \frac{u\eta^2 \sqrt{u^2 + v^2}}{gh},S_{0y} = \frac{\partial z}{\partial y},S_{fx} = \frac{v\eta^2\sqrt{u^2 + v^2}}{gh}$
Unity上的实现

久违的Unity终于又出现了。代码可在这个仓库的09章找到。
。=我分别用CPU上的C#计算了LaxF。我分别用CPU上的C#计算了LaxFriedchs方法和MacCormack方法，对应着ShallowWater.cs。以及用GPU上的ComputeShader计算了LaxFriedchs方法，对应着ComputeManager.cs。GameObject上只要勾选其中一个脚本就行了。代码写起来并没有什么特别的地方，和python差不多，照搬就行了。
ShallowWater.cs我设置了四种模拟情况，也就是watercondition。第一种就是Dam Break，不同地方的初始高度值不同
if (i < nmax / 2) h[i, j] = 2.0f; else h[i, j] = 1.0f;
上面代码在第60行。这种情况模拟的是左侧和右侧原本有一个无限高的大坝，现在大坝突然消失的情况。最终结果如下：
第二种就是下雨了，也就是封图，相应代码可在第214~234行找到，我觉得意外地好看，比波动方程更加真实自然，用python或matlab或openfoam可是看不到这么漂亮的景象哦。于是录了个150秒的放在了B站了https://www.bilibili.com/video/BV1Ry4y167MV求三连~
第三种情况就是带有坡度的DamBreak。最终迭代完的结果如下图，红色三角形就是我设置的有坡度的土堆的地方，相应代码在第42~53行。
第四种情况是无坡度，但有反射边界条件的Dam Break
水面的着色器用的是《Unity着色器入门与精要》第十章的Glass Refraction。不过这个场景还有三处可以改进的地方：
1.只实现了LaxFriedchs和MacCormack，我实现其它方法时遇到了一些问题暂时无法解决。
2.我在使用ComputeShader的时候发现居然比CPU还要慢，查了一圈发现大概是GetData要等待的时间太久了。
3.现在水面的美术效果只有静态环境映射，还有很多效果比如焦散和水下体积光都还没有实现。我说的就是这个AsehesL/UnityWaveEquation
勇敢的读者们，何不尝试挑战一下呢？
变化的时间步

对于如下不考虑坡度的非线性浅水波方程：

$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0\\ U = \begin{bmatrix} h \\ hu \end{bmatrix}, F = \begin{bmatrix} hu \\ h^2 u + gh^2/2 \\\end{bmatrix}$
可以将其拟线性化(quasi-linear)

$AU = F，A = \begin{bmatrix} 0 & 1\\ -u^2h + gh & 2u \end{bmatrix}\\$ 矩阵A的特征值为

$\lambda_1 = u - \sqrt{gh} \qquad \lambda_2 = u + \sqrt{gh}\\$ 对于浅水波方程这样的双曲线(hyberbolic)方程来说，这个特征值就是它们波中不同组分的传播速度，也叫特征速度( characteristic speeds)。其中u就是形成的水波的速度，而sqrt(gh)就是水流本身的速度。关于这两者请参考本篇的超临界流和亚临界流部分。正如Fluid Mechanics 101在他的视频https://www.youtube.com/watch?v=WBWY46ynRk0&t=1450s中说的那样，我们不希望CFL数大于一，也就是希望数值传播一次不要超过一个格子，也就是速度乘以单位时间不要超过一个格子长度

$CFL = \frac{\lambda_{max} \Delta t}{\Delta x} < 1 \quad \Rightarrow \quad \Delta t = CFL * \Delta x * \lambda_{max}(保证CFL < 1)\\$ 这里CFL数和dx都是一开始就固定好的，变化的只有特征值和dt。写成代码就是
dt = cfl*dx/max(abs(u[k,:]) + np.sqrt(abs(10*h[k,:])))到以后模拟更复杂的流体时，网格长度也需要发生变化，所以变化的时间步非常方便。
边界条件

不同边界条件的代码在ShallowWater.cs的第363~413行，主要有三类边界条件。
u[0] = -u[1];//速度的反射边界条件
u[0] = u[1];//速度的自由边界条件，梯度为零即黎曼边界条件
u[0] = u[N-1];//速度的周期边界条件
h[0] = h[1];//无论哪种边界条件，边界上高度h一直保持梯度为零
超临界流和亚临界流

水波本身是可以携带信息的。如果在水流中产生一个水波，可以根据这个水波是否能传递到上游而区分不同的水流。因此我们定义一个Froude数，其中u是水波的速度，而sqrt(gh)是水流的速度。

$froude = \frac{u}{\sqrt{gh}}\\ froude < 1 : supercritial \space flow\\ froude > 1 : subcritial \space flow\\$ 如果波的速度u大于水流速度sqrt(gh)，因此下游形成的波可以影响到上游，这就是亚临界流(subcritical flow)，特点是水深而水流速度慢。而如果波的速度比水流速度小，因此下游形成的波无法影响到上游，就是超临界流(supercritical flow)，特点是水浅而水流速度快。
仅仅这么介绍非常无聊，来看一点好玩的。我们可以做一个水箱，左侧保持较深的水，并且在中间将左侧与右侧隔开，隔板比左侧水的高度稍微低一点，让左侧的水能流到右侧。这样左侧的水就是亚临界流，而右侧的水就是速度很快的超临界流。
关于超/亚临界流推荐这个视频https://www.youtube.com/watch?v=7AKNR881BFE，上面两张图片也是来源于此。于是我们就可以用浅水波方程来模拟了。在ShallowWater.cs中，第三种情况就是专门模拟超/亚临界流的。如果将grad变量改为1.0，就模拟了一个在水底的非常陡峭的山坡，水流最终的高度也会受到影响，如下图，红线就是模拟的山坡处。
运行clatterrr/CFDcodepython 里面的Shallow1dSlopeLaxFriedrichs.py，也可以发现，迭代了几十次多次后，速度场在有坡度的地方发生了剧烈变化。
如果将坡度改为0.12让它不那么陡峭，可以看到最终水面的高度在山坡附近仍然有一些变化。
事实上水中的障碍物一般都会让水流的旋度发生变化，造成湍流，不过模拟湍流并不是本篇的内容。与超/亚临界流联系很紧密的就是Open Channel Flow了，但那也是另一个话题了。
介绍非线性方程的第七，八，九篇至此结束。下一篇，我计划尝试新的模拟方法，也就是光滑粒子动力学，并且重点介绍网上大部分教程都没怎么说到的核函数。
上一篇：光影帽子：
参考

大部分参考都已经在前两篇给出了，所以新的参考很少
《ON NUMERICAL SOLUTIONS OF THE SHALLOW WATER WAVE EQUATIONS》
《A kinetic scheme for the Saint-Venant system with a source term》
《Shallow Water Hydraulics》by Oscar Castro-Orgaz Willi H. Hager这本书后有一些关于超临界流和亚临界流模拟的内容
《Unsteady Flow in Open Channels》 by Jurjen A. Battjes, Robert Jan Labeur

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

[笔记] 【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak

本帖子中包含更多资源