unity 高斯模糊_【游戏流体力学基础及Unity代码（八）】激波捕捉法和RiemannSolver...

本篇介绍各种解决burgers方程的数值方法，相应的代码可在clatterrr/CFDcodepython找到。

上一章介绍了最简单的非线性方程，即一维无粘性burgers方程，但我们将它稍微改写一下

微分项用泰勒展开微分实际上得到了

上式是精确解。但是把右侧的值全算出来不现实。如果只取等式右侧第一项，就是有限差分(finite difference)形式的迎风格式(Upwind Scheme)。

一阶upwind格式的缺点就是舍弃了太多泰勒展开项，所以很不精确。而有限差分法规定每个网格由一个离散的值表示，在处理有关激波的不连续问题时并不那么自然，方便。所以有限体积(finite volume)法便被提了出来，这种方法将一个格子里的值看成连续的，从前到后积分后再除以格子长度就是格子的值，这种方法在处理不连续问题时有很大优势。比如最简单的有限体积法下的Lax-Friedrichs方法：

这里的f就是flux的意思，它的位置处在格子之间，它所代表的意思就是两个格子之间的流量。假如有n个网格，那么就应该有n+1个f。代码可以在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxFriedchs.py

然而这种方法也只有一阶精度，并且有一些人工粘性(artifial viscosity)所导致的耗散(disspative)。它看起来像是这样的。虚线是准确的解，而其它颜色的线是不同情况下的Lax-Friedrhs。人工粘性就像高斯模糊样，让尖锐的部分变得平滑了，然而我们要模拟的激波也消失了，这显然不是一种很好的方法。

所以我们继续尝试新方法，即Lax-Wendroff方法:

Lax-Wendroff方法有二阶精度，并且人工粘性的影响非常小，二阶精度在大多数时候都够用了。其代码可在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxWendroff.py。与Lax-Wendroff方法很相似的MacCormack方法也经常被使用：

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。如下图的左下，浅蓝是初始位置，深蓝是半个时间步后的位置，深蓝画得有点偏为了方便比较。

最终的五种情况如下：

这样极大值附近就不会出现比极大值还大的值了，极小值也是如此，振荡自然也就消失了。这就是一种Total Variation Diminishing方法，简称TVD方法，可以用于消除振荡。另一种分情况讨论的方法是Roe方法。标准的Roe方法要先线性化，算出一堆特征值和矩阵，然后再计算f，但本篇并不打算讨论这个，所以仅仅用一种简单的方法代替这个，也就是

上面的a其实就是上篇所说的激波的速度，因此也被称为Roe Speed。还有一种是HLL方法。

这就是早期人类用数值模拟驯服激波的珍贵方法。【误】

然而Godunov，Roe和HLL方法的准确度很差，都要先判断激波两侧的速度然后决定F的值。速度仅仅被划分为两个区间，即大于零或小于零，F也只能从一开始就设定好的几个结果中选一个，因此只有一阶精度。能不能不要这么一刀切，而是根据向前速度与向后速度的比例来确定最终的F呢？

我们要继续尝试新方法，也就是通量限制器(flux limiter)，也叫slop limiter。比如对于之前的LaxWendroff方法，要计算半步长的U，实际上应该用网格中心的U，加上自己斜率乘以网格长的二分之一，向前的斜率由前一格的值减去这一格的值得到，也就是下式

那个phi，就是flux limiter。对于LaxWendroff来说，它永远是一，因此也导致了振荡问题，这样不好。所以要继续尝试不同的方法。比如一个叫flux limiter的方法。因此我们先定义向前的斜率与向后的斜率之比为r，用r来决定phi。简便起见，下面的数学公式用u来代替U - dtdx*F。

继续分情况讨论，如果r < 0，也就是前后斜率的符号不一样，那就说明出现了局部最大值或最小值，我们不能让半步的值比局部最大值还大，也不能比局部最小值还小，那么就定义此时的斜率为零，也就是 phi = 0。如果向前的斜率与向后的斜率符号相同，比如都是正数，但前向的斜率绝对值略小，那么就要注意向前半步的值不能比前一格的值还要大，因为此时前一格可能是局部最大值。所以此时phi = r。如果向前斜率的绝对值比向后的更大，那么此时我们就使用laxWendroff所使用的phi = 1就行了，这就是minmod limiter。下图红线代表半步长的值。

表达式可以写成下面两种形式：

但实际写代码的时候更常使用下面这种函数的形式

除了minmod之外，还可以选择其它的flux limiter，最常用的如下：

上面的改进是基于LaxWendroff。但也可以使用另一种方法，也就是Monotonic Upstream-centered Scheme for Conservation Laws，简称MUSCL。它的就是在计算半步F的时候，不像之前那样使用此处的半步U，而且同时将而是把相邻网格所计算的此处的半步值也计算进来。

如下图，i+1/2处的值，即可以由左边U_i计算得到U^L，也可以由U_{i+1}计算得到U^R，而MUSCL方法将这两种值的影响都考虑了进来。

而将两个值都考虑进来后半步长F的具体算法，仍然可以由godunov或其它的flux limiter算法计算。比如由MUSCL + minmod来计算burgers方程，那么公式如下：

用上面介绍的所有方法解1dBurgers问题的代码都可以clatterrr/CFDcodepython找到。除了上面这些显式方法外，还有隐式方法也可以用，比如Runge-Kutta法，本篇不再介绍了。这些方法主要用来解决本系列第四篇提到的Euler方程，第七篇提到的Burgers方程，以及下一篇第九篇提到的浅水波方程。了解这些求解非线性方程的方法之后，就可以开始在unity上模拟了。

另外，我对一些概念理解还不够深入，所以本篇和代码可能有一些错误。还请大佬们发现后指导我一下。

上一篇：光影帽子：【游戏流体力学基础及Unity代码（七）】车流量问题，非线性水波以及burgers方程

下一篇：光影帽子：【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak

参考

书籍

[1]《Computational Fluid Mechanics and Heat Transfer》 by Anderson, Dale Pletcher, Richard H. Tannehill, John C

[2]《Finite Volume Methods for Hyperbolic Problems》by Randall J. LeVeque

[3]《Handbook of Numerical Methods for Hyperbolic Problems Basic and Fundamental Issues》by Rémi Abgrall and Chi-Wang Shu

[4]《Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction》 by Eleuterio F. Toro非常经典的书

[5]《Numerical Methods for Conservation Laws From Analysis to Algorithm》 by Jan S. Hesthaven

[6]《Solving Hyperbolic Equations with Finite Volume Methods》 by M. Elena Vázquez-Cendón 这书后面有一些代码可以参考

[7]《Numerical Solution of Hyperbolic Conservation Laws》by John A. Trangenstein

[8] https://www.ita.uni-heidelberg.de/~dullemond/lectures/num_fluid_2011/ 这里的讲义都挺不错的

[9]《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