【游戏流体力学基础及Unity代码（八）】激波捕捉法和RiemannSolver
本篇介绍各种解决burgers方程的数值方法，相应的代码可在clatterrr/CFDcodepython找到。上一章介绍了最简单的非线性方程，即一维无粘性burgers方程，但我们将它稍微改写一下
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+%5Cfrac%7B%5Cpartial+%7D%7B%5Cpartial+t%7D%28u%29+%2B+%5Cfrac%7B%5Cpartial+%7D%7B%5Cpartial+x%7D%28%5Cfrac%7B1%7D%7B2%7Du%5E2%29+%3D+%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+u+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5C%2C+%5Cspace+%5Cspace+%5Cspace+%5Cspace+F+%3D+%5Cfrac%7B1%7D%7B2%7Du%5E2 微分项用泰勒展开微分实际上得到了
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%3D+%5Cfrac%7Bu_%7Bi%2B1%7D+-+u_i%7D%7B%5CDelta+x%7D+-+%5Cfrac%7B%5CDelta+x%7D%7B2%7D%28%5Cfrac%7B%5Cpartial%5E2+u%7D%7B%5Cpartial+x%5E2%7D%29+-+%5Cfrac%7B%5CDelta+x%7D%7B6%7D%28%5Cfrac%7B%5Cpartial%5E3+u%7D%7B%5Cpartial+x%5E3%7D%29+-+...%5C%5C+ 上式是精确解。但是把右侧的值全算出来不现实。如果只取等式右侧第一项，就是有限差分(finite difference)形式的迎风格式(Upwind Scheme)。
https://www.zhihu.com/equation?tex=%5Cfrac%7BU_i%5E%7Bt%2B1%7D+-+U_i%5Et%7D%7B%5CDelta+t%7D+%2B+%5Cfrac%7BF_%7Bi%7D%5Et+-+F_%7Bi-1%7D%5Et%7D%7B%5CDelta+x%7D+%3D+0%5C%5C 一阶upwind格式的缺点就是舍弃了太多泰勒展开项，所以很不精确。而有限差分法规定每个网格由一个离散的值表示，在处理有关激波的不连续问题时并不那么自然，方便。所以有限体积(finite volume)法便被提了出来，这种方法将一个格子里的值看成连续的，从前到后积分后再除以格子长度就是格子的值，这种方法在处理不连续问题时有很大优势。比如最简单的有限体积法下的Lax-Friedrichs方法：
https://www.zhihu.com/equation?tex=U_i%5E%7Bt%2B1%7D+%3D+U_i%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D%28f_%7Bi%2B1%2F2%7D%5Et+-+f_%7Bi-1%2F2%7D%5Et%29%5C%5C+f_%7Bi-1%2F2%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28F_%7Bi%7D%5Et+%2B+F_%7Bi-1%7D%5Et%29+-+%5Cfrac%7B%5CDelta+x%7D%7B2%5CDelta+t%7D%28U_%7Bi%7D%5Et+-+U_%7Bi-1%7D%5Et%29
这里的f就是flux的意思，它的位置处在格子之间，它所代表的意思就是两个格子之间的流量。假如有n个网格，那么就应该有n+1个f。代码可以在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxFriedchs.py
然而这种方法也只有一阶精度，并且有一些人工粘性(artifial viscosity)所导致的耗散(disspative)。它看起来像是这样的。虚线是准确的解，而其它颜色的线是不同情况下的Lax-Friedrhs。人工粘性就像高斯模糊样，让尖锐的部分变得平滑了，然而我们要模拟的激波也消失了，这显然不是一种很好的方法。
所以我们继续尝试新方法，即Lax-Wendroff方法:
https://www.zhihu.com/equation?tex=U_i%5E%7Bt%2B1%7D+%3D+U_i%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D%28F_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+-+F_%7Bi-1%2F2%7D%5E%7Bt%2B1%2F2%7D%29%5C%5C+F_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28U_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D%29%5E2%28burgers%E6%96%B9%E7%A8%8B%29%5C%5C+U_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28U_%7Bi%2B1%7D%5Et+%2B+U_%7Bi%7D%5Et%29+-+%5Cfrac%7B%5CDelta+t%7D%7B2%5CDelta+x%7D%28F_%7Bi%2B1%7D%5Et+-+F_%7Bi%7D%5Et%29
Lax-Wendroff方法有二阶精度，并且人工粘性的影响非常小，二阶精度在大多数时候都够用了。其代码可在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxWendroff.py。与Lax-Wendroff方法很相似的MacCormack方法也经常被使用：
https://www.zhihu.com/equation?tex=U_i%5E%7Bt%2B1%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28U_i%5Et+%2B+U_i%5E%7Bt%2B1%2F2%7D%29+-+%5Cfrac%7B%5CDelta+t%7D%7B2%5CDelta+x%7D%28F_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+-+F_%7Bi-1%2F2%7D%5E%7Bt%2B1%2F2%7D%29%5C%5C+U_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+%3D+U_%7Bi%7D%5E%7Bt%7D+-+%5Cfrac%7B%5CDelta+t%7D%7B2%5CDelta+x%7D%28F_%7Bi%2B1%7D%5Et+-+F_%7Bi%7D%5Et%29 LaxWendroff和MacCormack都是二阶精度，看起来很方便好用，然而有着隐藏的问题，那就是振荡(oscillation)。它看起来像是这样的，虚线是准确的解，其它颜色则是在不同情况下的LaxWendroff方法，在临近激波的时候，虽然人工粘性几乎没了，但解的跳动非常严重，忽上忽下，偏离准确结果非常多。
振荡出现的原因也非常简单。比如我们要更新U_i的值，就要考虑它的U_i+1/2和U_i-1/2的值，这两个值又和U_i-1和U_i+1有关。所以U_i实际上的增加量，就是它自己的斜率，而斜率通过U_i-1和U_i+1计算出来。这就是问题所在。如下图，黑线是精确的解析解，而蓝线是近似的数值解。当U_i是全局最高值，理论上不能再有值比这个值高，而由于U_i-1小于U_i+1导致U_i的斜率为正的话，那么U_i+1/2会比U_i的值还要高，如下图，红线就是斜率，红点就是不该出现的高点，这就造成了振荡。
如何消除这种振荡？Godunov直接采用一种非常简单的方法，也就是我们小学就学过的分情况讨论。比如对于下图，黑色箭头所指的位置，半个时间步后的值应该是多少？如果uL和uR都大于零，那么它们都会向右移动，那么结果就是uL。如下图的左下，浅蓝是初始位置，深蓝是半个时间步后的位置，深蓝画得有点偏为了方便比较。
最终的五种情况如下：
https://www.zhihu.com/equation?tex=U_%7Bi%2B1%2F2%7D+%3D+uL%2C%E5%BD%93+uR+%3E+0%2CuL+%3E+0%5C%5C+U_%7Bi%2B1%2F2%7D+%3D+uR%2C%E5%BD%93+uR+%3C+0%2CuL+%3C+0%5C%5C+U_%7Bi%2B1%2F2%7D+%3D+0%2C%E5%BD%93+uL+%3C+0+%3C+uR%5C%5C+U_%7Bi%2B1%2F2%7D+%3D+uL%2C%E5%BD%93+uL+%3E+0+%3E+uR+%5Cspace%E4%B8%94+%5Cfrac%7BuL+%2B+uR%7D%7B2%7D+%3E+0%5C%5C+U_%7Bi%2B1%2F2%7D+%3D+uR%2C%E5%BD%93+uL+%3E+0+%3E+uR+%5Cspace+%E4%B8%94+%5Cfrac%7BuL+%2B+uR%7D%7B2%7D+%3C+0%5C%5C 这样极大值附近就不会出现比极大值还大的值了，极小值也是如此，振荡自然也就消失了。这就是一种Total Variation Diminishing方法，简称TVD方法，可以用于消除振荡。另一种分情况讨论的方法是Roe方法。标准的Roe方法要先线性化，算出一堆特征值和矩阵，然后再计算f，但本篇并不打算讨论这个，所以仅仅用一种简单的方法代替这个，也就是
https://www.zhihu.com/equation?tex=f_%7Bi%2B1%2F2%7D%3D+%5Cbegin%7Bcases%7D+F%28U_i%29%2C%E5%BD%93a_%7Bi%2B1%2F2%7D+%3E%3D+0%5C%5C+F%28U_%7Bi%2B1%7D%29%2C%E5%BD%93a_%7Bi%2B1%2F2%7D+%3C+0%5C%5C+%5Cend%7Bcases%7D%5C%5C+a_%7Bi%2B1%2F2%7D+%3D+%5Cfrac%7BF%28U_%7Bi%2B1%7D%29+-+F%28U_%7Bi%7D%29%7D%7BU_%7Bi%2B1%7D+-+U_%7Bi%7D%7D
上面的a其实就是上篇所说的激波的速度，因此也被称为Roe Speed。还有一种是HLL方法。
https://www.zhihu.com/equation?tex=f_%7Bi%2B1%2F2%7D%3D+%5Cbegin%7Bcases%7D+F_%7Bi%7D%2C%E5%BD%93sl+%3E%3D+0%5C%5C+F_%7Bhll%7D%2C%E5%BD%93sl+%3C%3D+0+%3C%3D+sr%5C%5C+F_%7Bi%2B1%7D%2C%E5%BD%93+sr+%3C%3D+0+%5Cend%7Bcases%7D%5C%5C+F_%7Bhll%7D+%3D+%5Cfrac%7Bsr%2AF_i+-+sl+%2A+F_%7Bi%2B1%7D+%2B+sr+%2A+sl+%2A+%28u_%7Bi%2B1%7D+-+u_%7Bi%7D%29%7D%7Bsr+-+sl%7D%5C%5C+%E5%AF%B9%E4%BA%8EBurgers%E6%96%B9%E7%A8%8B%EF%BC%9Asl+%3D+U_%7Bi%7D+%5Cspace+%5Cspace+%5Cspace+sr+%3D+U_%7Bi%2B1%7D
这就是早期人类用数值模拟驯服激波的珍贵方法。【误】
然而Godunov，Roe和HLL方法的准确度很差，都要先判断激波两侧的速度然后决定F的值。速度仅仅被划分为两个区间，即大于零或小于零，F也只能从一开始就设定好的几个结果中选一个，因此只有一阶精度。能不能不要这么一刀切，而是根据向前速度与向后速度的比例来确定最终的F呢？
我们要继续尝试新方法，也就是通量限制器(flux limiter)，也叫slop limiter。比如对于之前的LaxWendroff方法，要计算半步长的U，实际上应该用网格中心的U，加上自己斜率乘以网格长的二分之一，向前的斜率由前一格的值减去这一格的值得到，也就是下式
https://www.zhihu.com/equation?tex=U_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28U_%7Bi%2B1%7D%5Et+%2B+U_%7Bi%7D%5Et%29+-+%5Cfrac%7B%5CDelta+t%7D%7B2%5CDelta+x%7D%28F_%7Bi%2B1%7D%5Et+-+F_%7Bi%7D%5Et%29+%5C%5C+%5Clarge++%3D+%5Cfrac%7B1%7D%7B2%7D%28%28U_%7Bi%2B1%7D%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7DF_%7Bi%2B1%7D%5Et%29+-+%28U_%7Bi%7D%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7DF_%7Bi%7D%5Et%29%29%5Cphi++%2B+U_%7Bi%7D%5Et
那个phi，就是flux limiter。对于LaxWendroff来说，它永远是一，因此也导致了振荡问题，这样不好。所以要继续尝试不同的方法。比如一个叫flux limiter的方法。因此我们先定义向前的斜率与向后的斜率之比为r，用r来决定phi。简便起见，下面的数学公式用u来代替U - dt\dx*F。
https://www.zhihu.com/equation?tex=u_%7Bi%2B1%2F2%7D%5E%7Bt%2B1%2F2%7D++%3D+%5Cfrac%7B1%7D%7B2%7D%28u_%7Bi%2B1%7D%5Et+-+u_%7Bi%7D%5Et%29%5Cphi%28r%29+%2B+u_%7Bi%7D%5Et%5C%5C+%5Cphi%28r%29+%3D+%5Cphi%28%5Cfrac%7Bu_%7Bi%2B1%7D+-+u_%7Bi%7D%7D%7Bu_%7Bi%7D+-+u_%7Bi-1%7D%7D%29 继续分情况讨论，如果r < 0，也就是前后斜率的符号不一样，那就说明出现了局部最大值或最小值，我们不能让半步的值比局部最大值还大，也不能比局部最小值还小，那么就定义此时的斜率为零，也就是 phi = 0。如果向前的斜率与向后的斜率符号相同，比如都是正数，但前向的斜率绝对值略小，那么就要注意向前半步的值不能比前一格的值还要大，因为此时前一格可能是局部最大值。所以此时phi = r。如果向前斜率的绝对值比向后的更大，那么此时我们就使用laxWendroff所使用的phi = 1就行了，这就是minmod limiter。下图红线代表半步长的值。
表达式可以写成下面两种形式：
https://www.zhihu.com/equation?tex=minmode+%3A%5Cphi+%3D++%5Cbegin%7Bcases%7D++++++0+%2C%26+%5Ctext%7B%E5%A6%82%E6%9E%9C%7D+r+%3C+0%5C%5C++++++1+%2C%26+%5Ctext%7B%E5%A6%82%E6%9E%9C%7D+r+%3E+1%5C%5C++++++++++r+%2C%26+%5Ctext%7B%E5%85%B6%E5%AE%83%7D++%5C%5C++++%5Cend%7Bcases%7D%5C%5C+++minmode%3A%5Cphi+%3D+max%5B0%2Cmin%281%2Cr%29%5D+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Clim_%7Br-%3E%5Cinf%7D%5Cphi%28r%29+%3D+1 但实际写代码的时候更常使用下面这种函数的形式
https://www.zhihu.com/equation?tex=u_%7Bi%2B1%2F2%7D+%3D+u_i+%2B+%5Cfrac%7B1%7D%7B2%7D%2A+slope%5C%5C+++slope+%3D+minmod%28%5Cfrac%7Bu_%7Bi%2B1%7D+-+u_%7Bi%7D%7D%7B%5CDelta+x%7D%2C%5Cfrac%7Bu_%7Bi%7D+-+u_%7Bi-1%7D%7D%7B%5CDelta+x%7D%29%5C%5C+++minmod%28a%2Cb%29+%3D+%5Cbegin%7Bcases%7D+0%2C%E5%BD%93ab+%3C%3D+0%5C%5C+a%2C%E5%BD%93b+%3D+0+%E6%88%96+ab+%3E+0%E4%B8%94abs%28a%29+%3E%3D+abs%28b%29%5C%5C+a%2Fb%2C%E5%BD%93ab+%3E+0%E4%B8%94abs%28a%29+%3C+abs%28b%29%5C%5C+%5Cend%7Bcases%7D 除了minmod之外，还可以选择其它的flux limiter，最常用的如下：
https://www.zhihu.com/equation?tex=monotonized+central+%28MC%29%EF%BC%9A%5Cphi+%3D+max%5B0%2Cmin%282r%2C0.5%281%2Br%29%2C2%29%5D+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Clim_%7Br-%3E%5Cinf%7D%5Cphi%28r%29+%3D+2%5C%5C+superbee%EF%BC%9A%5Cphi+%3D+max%5B0%2Cmin%282r%2C1%29%2Cmin%28r%2C2%29%5D+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Clim_%7Br-%3E%5Cinf%7D%5Cphi%28r%29+%3D+2%5C%5C+van+Leer%EF%BC%9A%5Cphi+%3D+%5Cfrac%7Br+%2B+%7Cr%7C%7D%7B1+%2B+%7Cr%7C%7D+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Cspace+%5Clim_%7Br-%3E%5Cinf%7D%5Cphi%28r%29+%3D+2
上面的改进是基于LaxWendroff。但也可以使用另一种方法，也就是Monotonic Upstream-centered Scheme for Conservation Laws，简称MUSCL。它的就是在计算半步F的时候，不像之前那样使用此处的半步U，而且同时将而是把相邻网格所计算的此处的半步值也计算进来。
https://www.zhihu.com/equation?tex=U_i%5E%7Bt%2B1%7D+%3D+U_i%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D%28F_%7Bi%2B1%2F2%7D%5Et+-+F_%7Bi-1%2F2%7D%5Et%29%5C%5C+MUSCL%E6%96%B9%E6%B3%95+%EF%BC%9AF_%7Bi%2B1%2F2%7D%5E%7Bt%7D+%3D+F%28U_%7Bi%2B1%2F2%7D%5EL%2CU_%7Bi%2B1%2F2%7D%5ER%29%5C%5C+%E4%B9%8B%E5%89%8D%E7%9A%84%E6%96%B9%E6%B3%95+%EF%BC%9AF_%7Bi%2B1%2F2%7D%5E%7Bt%7D+%3D+F%28U_%7Bi%2B1%2F2%7D%5Et%29 如下图，i+1/2处的值，即可以由左边U_i计算得到U^L，也可以由U_{i+1}计算得到U^R，而MUSCL方法将这两种值的影响都考虑了进来。
而将两个值都考虑进来后半步长F的具体算法，仍然可以由godunov或其它的flux limiter算法计算。比如由MUSCL + minmod来计算burgers方程，那么公式如下：
https://www.zhihu.com/equation?tex=U_i%5E%7Bt%2B1%7D+%3D+U_i%5Et+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D%28F_%7Bi%2B1%2F2%7D%5Et+-+F_%7Bi-1%2F2%7D%5Et%29%5C%5C+F_%7Bi%2B1%2F2%7D%5Et+%3D+%5Cfrac%7B1%7D%7B2%7D%28F%28%5BU_%7Bi%2B1%2F2%7D%5E%7Bt%7D%5D%5EL%29+%2B+F%28%5BU_%7Bi%2B1%2F2%7D%5E%7Bt%7D%5D%5ER%29+-+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D%28%5BU_%7Bi%2B1%2F2%7D%5E%7Bt%7D%5D%5ER+-+%5BU_%7Bi%2B1%2F2%7D%5E%7Bt%7D%5D%5EL%29%29%5C%5C+%5BU_%7Bi%2B1%2F2%7D%5E%7Bt%7D%5D%5EL+%3D+U_%7Bi%7D%5Et+%2B+%5Cfrac%7B%5CDelta+x%7D%7B2%7D%2Aminmod%28%5Cfrac%7BU_%7Bi%2B1%7D+-+U_%7Bi%7D%7D%7B%5CDelta+x%7D%2C%5Cfrac%7BU_%7Bi%7D+-+U_%7Bi-1%7D%7D%7B%5CDelta+x%7D%29%5C%5C+%5BU_%7Bi-1%2F2%7D%5E%7Bt%7D%5D%5ER+%3D+U_%7Bi%7D%5Et+-+%5Cfrac%7B%5CDelta+x%7D%7B2%7D%2Aminmod%28%5Cfrac%7BU_%7Bi%2B1%7D+-+U_%7Bi%7D%7D%7B%5CDelta+x%7D%2C%5Cfrac%7BU_%7Bi%7D+-+U_%7Bi-1%7D%7D%7B%5CDelta+x%7D%29
用上面介绍的所有方法解1dBurgers问题的代码都可以clatterrr/CFDcodepython找到。除了上面这些显式方法外，还有隐式方法也可以用，比如Runge-Kutta法，本篇不再介绍了。这些方法主要用来解决本系列第四篇提到的Euler方程，第七篇提到的Burgers方程，以及下一篇第九篇提到的浅水波方程。了解这些求解非线性方程的方法之后，就可以开始在unity上模拟了。
另外，我对一些概念理解还不够深入，所以本篇和代码可能有一些错误。还请大佬们发现后指导我一下。
上一篇：光影帽子：【游戏流体力学基础及Unity代码（七）】车流量问题，非线性水波以及burgers方程
下一篇：光影帽子：【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak
参考
书籍
《Computational Fluid Mechanics and Heat Transfer》 by Anderson, Dale Pletcher, Richard H. Tannehill, John C
《Finite Volume Methods for Hyperbolic Problems》by Randall J. LeVeque
《Handbook of Numerical Methods for Hyperbolic ProblemsBasic and Fundamental Issues》by Rémi Abgrall and Chi-Wang Shu
《Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction》 by Eleuterio F. Toro非常经典的书
《Numerical Methods for Conservation Laws From Analysis to Algorithm》 by Jan S. Hesthaven
《Solving Hyperbolic Equations with Finite Volume Methods》 by M. Elena Vázquez-Cendón 这书后面有一些代码可以参考
《Numerical Solution of Hyperbolic Conservation Laws》by John A. Trangenstein
https://www.ita.uni-heidelberg.de/~dullemond/lectures/num_fluid_2011/ 这里的讲义都挺不错的
《Exact solution of the Riemann problem for the shallow water equations》
推荐代码，包括了下一篇的浅水波方程。每一份都是经过精心挑选，简短而又清晰的适合初学者的代码~
https://www.mathworks.com/matlabcentral/fileexchange/46475-1d-shallow-water-equations-dam-break
https://math.la.asu.edu/~gardner/526.html
http://folk.ntnu.no/leifh/teaching/tkt4140/
https://github.com/zerothphase/DamBreakSimulation
https://github.com/sagarbhatt0904/Dam-Break-Problem
https://github.com/jostbr/shallow-water
https://github.com/cangyu/Riemann-Solvers
https://github.com/JMFaundez/Shallow_waterFVM
https://github.com/mehradans92/Shallow-water-dynamics
https://github.com/sammaster/sw-lax-friedrichs-python
https://github.com/iCFD/NonlinearScalarAdvectionLaws
https://github.com/shivams15/DamBreak1D
https://github.com/celinezou/AM205-Numerical-Method-Project
页:
[1]