【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak
开局一张图！浅水波方程
当水波的波长远远大于水的深度时，水波就表现出来一些特别的性质，比如江河就是典型的浅水，可以用浅水波方程(Shallow Water Equation)来模拟，这个方程又叫Saint Venant Equation。
要推导一维的浅水波方程，同样要从连续性方程和动力方程入手。对于X轴上某个坐标点来说，这点此时刻水深度和上时刻水深度不同，这是因为与附近水的深度差以及水的速度。这里假设速度在z方向上没有变化，或者说速度在z方向上的变化非常小以至于可以被忽略。于是就写出了浅水波方程中的连续性方程。这里的速度和水深有关，所以也是一个非线性方程。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+h%7D%7B%5Cpartial+t%7D+%2B+%5Cfrac%7B%5Cpartial+%28hu%29%7D%7B%5Cpartial+x%7D+%3D+0%5Ctag%7B1%7D
然后复习一下初中物理。水中任何一点的压力都是向四面八方的。将一块方形石头放入静水中，它的前后两面所受的压力可以抵消，左右两面也是如此，但是上下两面的压力因为重力的原因而无法抵消。这就是静水压力(hydrostatic pressure)。因此进入深水区需要特殊的潜水装备甚至是潜水艇，否则人类身体无法承受如此大的压力差。下面是简单的推导。
https://www.zhihu.com/equation?tex=p+%3D+%5Cfrac%7BF%7D%7BS%7D+%3D+%5Cfrac%7Bmg%7D%7BS%7D+%3D+%5Cfrac%7B%5Crho+V+g%7D%7BS%7D+%3D+%5Cfrac%7B%5Crho+S+h+g%7D%7BS%7D+%3D+%5Crho+g+h%5C%5C p是压力，F是受到的力，S是面积，m是质量，g是重力加速度，V是体积，h是水的深度。上式说明密度越大或水越深，那么受到的压力越大。在我们的实验中，密度和重力加速度是常数，只有水的深度会随着x轴或y轴变化。也就是说：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+x%7D+%3D+%5Crho+g+%5Cfrac%7B%5Cpartial+h%7D%7B%5Cpartial+x%7D%5C%5C 然后开始写动量方程。由于水的粘度很小，这里继续忽略掉粘性项。改写一下第四章提到的Euler方程
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+t%7D+%2B+u%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%3D+-g%5Cfrac%7B%5Cpartial+h%7D%7B%5Cpartial+x%7D%5C%5C
然后上式两边同乘以水深，就得到了浅水波方程中的动量方程
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial+t%7D%28hu%29+%2B+%5Cfrac%7B%5Cpartial+%7D%7B%5Cpartial+x%7D%28hu%5E2+%2B+%5Cfrac%7B1%7D%7B2%7Dgh%5E2%29++%3D+0%5Ctag%7B2%7D
1,2式有一些相同的地方，我们可以把它们写成下面这个式子，再加上坡度和摩擦力的影响，就得到了一维浅水波方程(Shallow Water Equation)。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+t%7D+%2B+%5Cfrac%7B%5Cpartial+F%7D%7B%5Cpartial+x%7D++%3D+S%5C%5C+U+%3D+%5Cbegin%7Bbmatrix%7D+h+%5C%5C+hu++%5Cend%7Bbmatrix%7D%2C+F+%3D+%5Cbegin%7Bbmatrix%7D+hu+%5C%5C+h%5E2+u+%2B+gh%5E2%2F2+%5C%5C%5Cend%7Bbmatrix%7D%2CS+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+gh%28S_0+-+S_f%29++%5Cend%7Bbmatrix%7D
等式右侧S代表浅水底坡度和摩擦力造成的影响。其中S0是坡度，Sf是摩擦力，其中\eta是Manning resistance coefficient取决于不同的摩擦规则与情况。
https://www.zhihu.com/equation?tex=S_0+%3D+%5Cfrac%7B%5Cpartial+z%7D%7B%5Cpartial+x%7D%2CS_f+%3D+%5Cfrac%7Bu%5E2%5Ceta%5E2%7D%7Bgh%7D%5C%5C
浅水波方程和Burgers方程很相似，其求解方法就是上一篇提到的各种方法。部分方法的python代码我也放到了与上一篇相同的仓库里。clatterrr/CFDcodepython。要注意的就是不要把高度初始值设置为零，否则到后面就要除以零了。而且浅水波方程使用变化的时间步dt可以获得更合理的结果，以及可以使用不同的边界条件。本篇稍后讨论。二维浅水波方程如下：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+t%7D+%2B+%5Cfrac%7B%5Cpartial+F%7D%7B%5Cpartial+x%7D+%2B+%5Cfrac%7B%5Cpartial+G%7D%7B%5Cpartial+y%7D+%3D+S%5C%5C+U+%3D+%5Cbegin%7Bbmatrix%7D+h+%5C%5C+hu+%5C%5C+hv+%5Cend%7Bbmatrix%7D%2C+F+%3D+%5Cbegin%7Bbmatrix%7D+hu+%5C%5C+h%5E2+u+%2B+gh%5E2%2F2+%5C%5C+huv+%5Cend%7Bbmatrix%7D%2CG+%3D+%5Cbegin%7Bbmatrix%7D+hv+%5C%5C+huv+%5C%5C+h%5E2+v+%2B+gh%5E2%2F2+%5Cend%7Bbmatrix%7D%2CS+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+gh%28S_%7B0x%7D+-+S_%7Bfx%7D%29+%5C%5C+gh%28S_%7B0y%7D+-+S_%7Bfx%7D%29+%5Cend%7Bbmatrix%7D%5C%5C+S_%7B0x%7D+%3D+%5Cfrac%7B%5Cpartial+z%7D%7B%5Cpartial+x%7D%2CS_%7Bfx%7D+%3D+%5Cfrac%7Bu%5Ceta%5E2+%5Csqrt%7Bu%5E2+%2B+v%5E2%7D%7D%7Bgh%7D%2CS_%7B0y%7D+%3D+%5Cfrac%7B%5Cpartial+z%7D%7B%5Cpartial+y%7D%2CS_%7Bfx%7D+%3D+%5Cfrac%7Bv%5Ceta%5E2%5Csqrt%7Bu%5E2+%2B+v%5E2%7D%7D%7Bgh%7D
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 = 2.0f; else h = 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
勇敢的读者们，何不尝试挑战一下呢？
变化的时间步
对于如下不考虑坡度的非线性浅水波方程：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+t%7D+%2B+%5Cfrac%7B%5Cpartial+F%7D%7B%5Cpartial+x%7D+%3D+0%5C%5C+U+%3D+%5Cbegin%7Bbmatrix%7D+h+%5C%5C+hu++%5Cend%7Bbmatrix%7D%2C+F+%3D+%5Cbegin%7Bbmatrix%7D+hu+%5C%5C+h%5E2+u+%2B+gh%5E2%2F2+%5C%5C%5Cend%7Bbmatrix%7D
可以将其拟线性化(quasi-linear)
https://www.zhihu.com/equation?tex=AU+%3D+F%EF%BC%8CA+%3D+%5Cbegin%7Bbmatrix%7D+0+%26+1%5C%5C+-u%5E2h+%2B+gh+%26+2u++%5Cend%7Bbmatrix%7D%5C%5C 矩阵A的特征值为
https://www.zhihu.com/equation?tex=%5Clambda_1+%3D+u+-+%5Csqrt%7Bgh%7D+%5Cqquad+%5Clambda_2+%3D+u+%2B+%5Csqrt%7Bgh%7D%5C%5C 对于浅水波方程这样的双曲线(hyberbolic)方程来说，这个特征值就是它们波中不同组分的传播速度，也叫特征速度( characteristic speeds)。其中u就是形成的水波的速度，而sqrt(gh)就是水流本身的速度。关于这两者请参考本篇的超临界流和亚临界流部分。正如Fluid Mechanics 101在他的视频https://www.youtube.com/watch?v=WBWY46ynRk0&t=1450s中说的那样，我们不希望CFL数大于一，也就是希望数值传播一次不要超过一个格子，也就是速度乘以单位时间不要超过一个格子长度
https://www.zhihu.com/equation?tex=CFL+%3D+%5Cfrac%7B%5Clambda_%7Bmax%7D+%5CDelta+t%7D%7B%5CDelta+x%7D+%3C+1+%5Cquad+%5CRightarrow+%5Cquad+%5CDelta+t+%3D+CFL+%2A+%5CDelta+x+%2A+%5Clambda_%7Bmax%7D%28%E4%BF%9D%E8%AF%81CFL+%3C+1%29%5C%5C 这里CFL数和dx都是一开始就固定好的，变化的只有特征值和dt。写成代码就是
dt = cfl*dx/max(abs(u) + np.sqrt(abs(10*h)))到以后模拟更复杂的流体时，网格长度也需要发生变化，所以变化的时间步非常方便。
边界条件
不同边界条件的代码在ShallowWater.cs的第363~413行，主要有三类边界条件。
u = -u;//速度的反射边界条件
u = u;//速度的自由边界条件，梯度为零即黎曼边界条件
u = u;//速度的周期边界条件
h = h;//无论哪种边界条件，边界上高度h一直保持梯度为零
超临界流和亚临界流
水波本身是可以携带信息的。如果在水流中产生一个水波，可以根据这个水波是否能传递到上游而区分不同的水流。因此我们定义一个Froude数，其中u是水波的速度，而sqrt(gh)是水流的速度。
https://www.zhihu.com/equation?tex=froude+%3D+%5Cfrac%7Bu%7D%7B%5Csqrt%7Bgh%7D%7D%5C%5C+froude+%3C+1+%3A+supercritial+%5Cspace+flow%5C%5C+froude+%3E+1+%3A+subcritial+%5Cspace+flow%5C%5C 如果波的速度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
页:
[1]