【游戏流体力学基础及Unity代码（八）】激波捕捉法和RiemannSolver

六翼天使494 · 发表于 2020-12-6 08:34

本篇介绍各种解决burgers方程的数值方法，相应的代码可在clatterrr/CFDcodepython找到。
上一章介绍了最简单的非线性方程，即一维无粘性burgers方程，但我们将它稍微改写一下

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = \frac{\partial }{\partial t}(u) + \frac{\partial }{\partial x}(\frac{1}{2}u^2) = \frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0\\ U = u \space \space \space \space \, \space \space \space \space F = \frac{1}{2}u^2$ 微分项用泰勒展开微分实际上得到了

$\frac{\partial u}{\partial x} = \frac{u_{i+1} - u_i}{\Delta x} - \frac{\Delta x}{2}(\frac{\partial^2 u}{\partial x^2}) - \frac{\Delta x}{6}(\frac{\partial^3 u}{\partial x^3}) - ...\\$ 上式是精确解。但是把右侧的值全算出来不现实。如果只取等式右侧第一项，就是有限差分(finite difference)形式的迎风格式(Upwind Scheme)。

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

$U_i^{t+1} = U_i^t - \frac{\Delta t}{\Delta x}(f_{i+1/2}^t - f_{i-1/2}^t)\\ f_{i-1/2} = \frac{1}{2}(F_{i}^t + F_{i-1}^t) - \frac{\Delta x}{2\Delta t}(U_{i}^t - U_{i-1}^t)$
这里的f就是flux的意思，它的位置处在格子之间，它所代表的意思就是两个格子之间的流量。假如有n个网格，那么就应该有n+1个f。代码可以在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxFriedchs.py
然而这种方法也只有一阶精度，并且有一些人工粘性(artifial viscosity)所导致的耗散(disspative)。它看起来像是这样的。虚线是准确的解，而其它颜色的线是不同情况下的Lax-Friedrhs。人工粘性就像高斯模糊样，让尖锐的部分变得平滑了，然而我们要模拟的激波也消失了，这显然不是一种很好的方法。
所以我们继续尝试新方法，即Lax-Wendroff方法:

$U_i^{t+1} = U_i^t - \frac{\Delta t}{\Delta x}(F_{i+1/2}^{t+1/2} - F_{i-1/2}^{t+1/2})\\ F_{i+1/2}^{t+1/2} = \frac{1}{2}(U_{i+1/2}^{t+1/2})^2(burgers方程)\\ U_{i+1/2}^{t+1/2} = \frac{1}{2}(U_{i+1}^t + U_{i}^t) - \frac{\Delta t}{2\Delta x}(F_{i+1}^t - F_{i}^t)$
Lax-Wendroff方法有二阶精度，并且人工粘性的影响非常小，二阶精度在大多数时候都够用了。其代码可在这里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxWendroff.py。与Lax-Wendroff方法很相似的MacCormack方法也经常被使用：

$U_i^{t+1} = \frac{1}{2}(U_i^t + U_i^{t+1/2}) - \frac{\Delta t}{2\Delta x}(F_{i+1/2}^{t+1/2} - F_{i-1/2}^{t+1/2})\\ U_{i+1/2}^{t+1/2} = U_{i}^{t} - \frac{\Delta t}{2\Delta x}(F_{i+1}^t - F_{i}^t)$ 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。如下图的左下，浅蓝是初始位置，深蓝是半个时间步后的位置，深蓝画得有点偏为了方便比较。
最终的五种情况如下：

$U_{i+1/2} = uL,当 uR > 0,uL > 0\\ U_{i+1/2} = uR,当 uR < 0,uL < 0\\ U_{i+1/2} = 0,当 uL < 0 < uR\\ U_{i+1/2} = uL,当 uL > 0 > uR \space且 \frac{uL + uR}{2} > 0\\ U_{i+1/2} = uR,当 uL > 0 > uR \space 且 \frac{uL + uR}{2} < 0\\$ 这样极大值附近就不会出现比极大值还大的值了，极小值也是如此，振荡自然也就消失了。这就是一种Total Variation Diminishing方法，简称TVD方法，可以用于消除振荡。另一种分情况讨论的方法是Roe方法。标准的Roe方法要先线性化，算出一堆特征值和矩阵，然后再计算f，但本篇并不打算讨论这个，所以仅仅用一种简单的方法代替这个，也就是

$f_{i+1/2}= \begin{cases} F(U_i),当a_{i+1/2} >= 0\\ F(U_{i+1}),当a_{i+1/2} < 0\\ \end{cases}\\ a_{i+1/2} = \frac{F(U_{i+1}) - F(U_{i})}{U_{i+1} - U_{i}}$
上面的a其实就是上篇所说的激波的速度，因此也被称为Roe Speed。还有一种是HLL方法。

$f_{i+1/2}= \begin{cases} F_{i},当sl >= 0\\ F_{hll},当sl <= 0 <= sr\\ F_{i+1},当 sr <= 0 \end{cases}\\ F_{hll} = \frac{sr*F_i - sl * F_{i+1} + sr * sl * (u_{i+1} - u_{i})}{sr - sl}\\ 对于Burgers方程：sl = U_{i} \space \space \space sr = U_{i+1}$
这就是早期人类用数值模拟驯服激波的珍贵方法。【误】
然而Godunov，Roe和HLL方法的准确度很差，都要先判断激波两侧的速度然后决定F的值。速度仅仅被划分为两个区间，即大于零或小于零，F也只能从一开始就设定好的几个结果中选一个，因此只有一阶精度。能不能不要这么一刀切，而是根据向前速度与向后速度的比例来确定最终的F呢？
我们要继续尝试新方法，也就是通量限制器(flux limiter)，也叫slop limiter。比如对于之前的LaxWendroff方法，要计算半步长的U，实际上应该用网格中心的U，加上自己斜率乘以网格长的二分之一，向前的斜率由前一格的值减去这一格的值得到，也就是下式

$U_{i+1/2}^{t+1/2} = \frac{1}{2}(U_{i+1}^t + U_{i}^t) - \frac{\Delta t}{2\Delta x}(F_{i+1}^t - F_{i}^t) \\ \large = \frac{1}{2}((U_{i+1}^t - \frac{\Delta t}{\Delta x}F_{i+1}^t) - (U_{i}^t - \frac{\Delta t}{\Delta x}F_{i}^t))\phi + U_{i}^t$
那个phi，就是flux limiter。对于LaxWendroff来说，它永远是一，因此也导致了振荡问题，这样不好。所以要继续尝试不同的方法。比如一个叫flux limiter的方法。因此我们先定义向前的斜率与向后的斜率之比为r，用r来决定phi。简便起见，下面的数学公式用u来代替U - dt\dx*F。

$u_{i+1/2}^{t+1/2} = \frac{1}{2}(u_{i+1}^t - u_{i}^t)\phi(r) + u_{i}^t\\ \phi(r) = \phi(\frac{u_{i+1} - u_{i}}{u_{i} - u_{i-1}})$ 继续分情况讨论，如果r < 0，也就是前后斜率的符号不一样，那就说明出现了局部最大值或最小值，我们不能让半步的值比局部最大值还大，也不能比局部最小值还小，那么就定义此时的斜率为零，也就是 phi = 0。如果向前的斜率与向后的斜率符号相同，比如都是正数，但前向的斜率绝对值略小，那么就要注意向前半步的值不能比前一格的值还要大，因为此时前一格可能是局部最大值。所以此时phi = r。如果向前斜率的绝对值比向后的更大，那么此时我们就使用laxWendroff所使用的phi = 1就行了，这就是minmod limiter。下图红线代表半步长的值。
表达式可以写成下面两种形式：

$minmode :\phi = \begin{cases} 0 ,& \text{如果} r < 0\\ 1 ,& \text{如果} r > 1\\ r ,& \text{其它} \\ \end{cases}\\ minmode:\phi = max[0,min(1,r)] \space \space \space \space \space \space \space \lim_{r->\inf}\phi(r) = 1$ 但实际写代码的时候更常使用下面这种函数的形式

$u_{i+1/2} = u_i + \frac{1}{2}* slope\\ slope = minmod(\frac{u_{i+1} - u_{i}}{\Delta x},\frac{u_{i} - u_{i-1}}{\Delta x})\\ minmod(a,b) = \begin{cases} 0,当ab <= 0\\ a,当b = 0 或 ab > 0且abs(a) >= abs(b)\\ a/b,当ab > 0且abs(a) < abs(b)\\ \end{cases}$ 除了minmod之外，还可以选择其它的flux limiter，最常用的如下：

$monotonized central (MC)：\phi = max[0,min(2r,0.5(1+r),2)] \space \space \space \space \space \space \space \lim_{r->\inf}\phi(r) = 2\\ superbee：\phi = max[0,min(2r,1),min(r,2)] \space \space \space \space \space \space \space \lim_{r->\inf}\phi(r) = 2\\ van Leer：\phi = \frac{r + |r|}{1 + |r|} \space \space \space \space \space \space \space \lim_{r->\inf}\phi(r) = 2$
上面的改进是基于LaxWendroff。但也可以使用另一种方法，也就是Monotonic Upstream-centered Scheme for Conservation Laws，简称MUSCL。它的就是在计算半步F的时候，不像之前那样使用此处的半步U，而且同时将而是把相邻网格所计算的此处的半步值也计算进来。

$U_i^{t+1} = U_i^t - \frac{\Delta t}{\Delta x}(F_{i+1/2}^t - F_{i-1/2}^t)\\ MUSCL方法：F_{i+1/2}^{t} = F(U_{i+1/2}^L,U_{i+1/2}^R)\\ 之前的方法：F_{i+1/2}^{t} = F(U_{i+1/2}^t)$ 如下图，i+1/2处的值，即可以由左边U_i计算得到U^L，也可以由U_{i+1}计算得到U^R，而MUSCL方法将这两种值的影响都考虑了进来。
而将两个值都考虑进来后半步长F的具体算法，仍然可以由godunov或其它的flux limiter算法计算。比如由MUSCL + minmod来计算burgers方程，那么公式如下：

$U_i^{t+1} = U_i^t - \frac{\Delta t}{\Delta x}(F_{i+1/2}^t - F_{i-1/2}^t)\\ F_{i+1/2}^t = \frac{1}{2}(F([U_{i+1/2}^{t}]^L) + F([U_{i+1/2}^{t}]^R) - \frac{\Delta t}{\Delta x}([U_{i+1/2}^{t}]^R - [U_{i+1/2}^{t}]^L))\\ [U_{i+1/2}^{t}]^L = U_{i}^t + \frac{\Delta x}{2}*minmod(\frac{U_{i+1} - U_{i}}{\Delta x},\frac{U_{i} - U_{i-1}}{\Delta x})\\ [U_{i-1/2}^{t}]^R = U_{i}^t - \frac{\Delta x}{2}*minmod(\frac{U_{i+1} - U_{i}}{\Delta x},\frac{U_{i} - U_{i-1}}{\Delta x})$
用上面介绍的所有方法解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

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

[笔记] 【游戏流体力学基础及Unity代码（八）】激波捕捉法和RiemannSolver

本帖子中包含更多资源