本篇介紹各種解決burgers方程的數(shù)值方法，相應(yīng)的代碼可在clatterrr/CFDcodepython找到。

上一章介紹了最簡單的非線性方程，即一維無粘性burgers方程，但我們將它稍微改寫一下

微分項(xiàng)用泰勒展開微分實(shí)際上得到了

上式是精確解。但是把右側(cè)的值全算出來不現(xiàn)實(shí)。如果只取等式右側(cè)第一項(xiàng)，就是有限差分(finite difference)形式的迎風(fēng)格式(Upwind Scheme)。

一階upwind格式的缺點(diǎn)就是舍棄了太多泰勒展開項(xiàng)，所以很不精確。而有限差分法規(guī)定每個(gè)網(wǎng)格由一個(gè)離散的值表示，在處理有關(guān)激波的不連續(xù)問題時(shí)并不那么自然，方便。所以有限體積(finite volume)法便被提了出來，這種方法將一個(gè)格子里的值看成連續(xù)的，從前到后積分后再除以格子長度就是格子的值，這種方法在處理不連續(xù)問題時(shí)有很大優(yōu)勢。比如最簡單的有限體積法下的Lax-Friedrichs方法：

然而這種方法也只有一階精度，并且有一些人工粘性(artifial viscosity)所導(dǎo)致的耗散(disspative)。它看起來像是這樣的。虛線是準(zhǔn)確的解，而其它顏色的線是不同情況下的Lax-Friedrhs。人工粘性就像高斯模糊樣，讓尖銳的部分變得平滑了，然而我們要模擬的激波也消失了，這顯然不是一種很好的方法。

所以我們繼續(xù)嘗試新方法，即Lax-Wendroff方法:

Lax-Wendroff方法有二階精度，并且人工粘性的影響非常小，二階精度在大多數(shù)時(shí)候都夠用了。其代碼可在這里找到https://github.com/clatterrr/CFDcodepython/blob/main/RiemannSolver/Burgers1dLaxWendroff.py。與Lax-Wendroff方法很相似的MacCormack方法也經(jīng)常被使用：

LaxWendroff和MacCormack都是二階精度，看起來很方便好用，然而有著隱藏的問題，那就是振蕩(oscillation)。它看起來像是這樣的，虛線是準(zhǔn)確的解，其它顏色則是在不同情況下的LaxWendroff方法，在臨近激波的時(shí)候，雖然人工粘性幾乎沒了，但解的跳動(dòng)非常嚴(yán)重，忽上忽下，偏離準(zhǔn)確結(jié)果非常多。

振蕩出現(xiàn)的原因也非常簡單。比如我們要更新U_i的值，就要考慮它的U_i+1/2和U_i-1/2的值，這兩個(gè)值又和U_i-1和U_i+1有關(guān)。所以U_i實(shí)際上的增加量，就是它自己的斜率，而斜率通過U_i-1和U_i+1計(jì)算出來。這就是問題所在。如下圖，黑線是精確的解析解，而藍(lán)線是近似的數(shù)值解。當(dāng)U_i是全局最高值，理論上不能再有值比這個(gè)值高，而由于U_i-1小于U_i+1導(dǎo)致U_i的斜率為正的話，那么U_i+1/2會(huì)比U_i的值還要高，如下圖，紅線就是斜率，紅點(diǎn)就是不該出現(xiàn)的高點(diǎn)，這就造成了振蕩。

關(guān)于這個(gè)問題，油管上有三個(gè)非常不錯(cuò)的視頻Need For Flux Limiter, Flux Limiters , How Does Flux Limiter Work，這三個(gè)視頻的博主是同一人，似乎是MIT的教師，他的有關(guān)計(jì)算流體力學(xué)的其它視頻質(zhì)量也非常高，推薦關(guān)注。

如何消除這種振蕩？Godunov直接采用一種非常簡單的方法，也就是我們小學(xué)就學(xué)過的分情況討論。比如對于下圖，黑色箭頭所指的位置，半個(gè)時(shí)間步后的值應(yīng)該是多少？如果uL和uR都大于零，那么它們都會(huì)向右移動(dòng)，那么結(jié)果就是uL。如下圖的左下，淺藍(lán)是初始位置，深藍(lán)是半個(gè)時(shí)間步后的位置，深藍(lán)畫得有點(diǎn)偏為了方便比較。

最終的五種情況如下：

這樣極大值附近就不會(huì)出現(xiàn)比極大值還大的值了，極小值也是如此，振蕩自然也就消失了。這就是一種Total Variation Diminishing方法，簡稱TVD方法，可以用于消除振蕩。另一種分情況討論的方法是Roe方法。標(biāo)準(zhǔn)的Roe方法要先線性化，算出一堆特征值和矩陣，然后再計(jì)算f，但本篇并不打算討論這個(gè)，所以僅僅用一種簡單的方法代替這個(gè)，也就是

上面的a其實(shí)就是上篇所說的激波的速度，因此也被稱為Roe Speed。還有一種是HLL方法。

這就是早期人類用數(shù)值模擬馴服激波的珍貴方法。【誤】

然而Godunov，Roe和HLL方法的準(zhǔn)確度很差，都要先判斷激波兩側(cè)的速度然后決定F的值。速度僅僅被劃分為兩個(gè)區(qū)間，即大于零或小于零，F也只能從一開始就設(shè)定好的幾個(gè)結(jié)果中選一個(gè)，因此只有一階精度。能不能不要這么一刀切，而是根據(jù)向前速度與向后速度的比例來確定最終的F呢？

我們要繼續(xù)嘗試新方法，也就是通量限制器(flux limiter)，也叫slop limiter。比如對于之前的LaxWendroff方法，要計(jì)算半步長的U，實(shí)際上應(yīng)該用網(wǎng)格中心的U，加上自己斜率乘以網(wǎng)格長的二分之一，向前的斜率由前一格的值減去這一格的值得到，也就是下式

那個(gè)phi，就是flux limiter。對于LaxWendroff來說，它永遠(yuǎn)是一，因此也導(dǎo)致了振蕩問題，這樣不好。所以要繼續(xù)嘗試不同的方法。比如一個(gè)叫flux limiter的方法。因此我們先定義向前的斜率與向后的斜率之比為r，用r來決定phi。簡便起見，下面的數(shù)學(xué)公式用u來代替U - dt\dx*F。

繼續(xù)分情況討論，如果r < 0，也就是前后斜率的符號不一樣，那就說明出現(xiàn)了局部最大值或最小值，我們不能讓半步的值比局部最大值還大，也不能比局部最小值還小，那么就定義此時(shí)的斜率為零，也就是 phi = 0。如果向前的斜率與向后的斜率符號相同，比如都是正數(shù)，但前向的斜率絕對值略小，那么就要注意向前半步的值不能比前一格的值還要大，因?yàn)榇藭r(shí)前一格可能是局部最大值。所以此時(shí)phi = r。如果向前斜率的絕對值比向后的更大，那么此時(shí)我們就使用laxWendroff所使用的phi = 1就行了，這就是minmod limiter。下圖紅線代表半步長的值。

表達(dá)式可以寫成下面兩種形式：

但實(shí)際寫代碼的時(shí)候更常使用下面這種函數(shù)的形式

除了minmod之外，還可以選擇其它的flux limiter，最常用的如下：

上面的改進(jìn)是基于LaxWendroff。但也可以使用另一種方法，也就是Monotonic Upstream-centered Scheme for Conservation Laws，簡稱MUSCL。它的就是在計(jì)算半步F的時(shí)候，不像之前那樣使用此處的半步U，而且同時(shí)將而是把相鄰網(wǎng)格所計(jì)算的此處的半步值也計(jì)算進(jìn)來。

如下圖，i+1/2處的值，即可以由左邊U_i計(jì)算得到U^L，也可以由U_{i+1}計(jì)算得到U^R，而MUSCL方法將這兩種值的影響都考慮了進(jìn)來。

而將兩個(gè)值都考慮進(jìn)來后半步長F的具體算法，仍然可以由godunov或其它的flux limiter算法計(jì)算。比如由MUSCL + minmod來計(jì)算burgers方程，那么公式如下：

用上面介紹的所有方法解1dBurgers問題的代碼都可以clatterrr/CFDcodepython找到。除了上面這些顯式方法外，還有隱式方法也可以用，比如Runge-Kutta法，本篇不再介紹了。這些方法主要用來解決本系列第四篇提到的Euler方程，第七篇提到的Burgers方程，以及下一篇第九篇提到的淺水波方程。了解這些求解非線性方程的方法之后，就可以開始在unity上模擬了。

另外，我對一些概念理解還不夠深入，所以本篇和代碼可能有一些錯(cuò)誤。還請大佬們發(fā)現(xiàn)后指導(dǎo)我一下。

宣傳：我創(chuàng)建了一個(gè)模擬流體交流Q群，歡迎各位喜歡自己編寫代碼模擬流體的小伙伴們?nèi)肴夯ハ嘟涣鲗W(xué)習(xí)！群號：1001290801

參考

僅僅看我寫的文章肯定是不夠的，所以我選了一些不錯(cuò)的書籍列在下面

[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非常經(jīng)典的書，超級推薦Riemann Solvers and Numerical Methods for Fluid Dynamics

[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

[9]《Exact solution of the Riemann problem for the shallow water equations》

[10]《計(jì)算淺水動(dòng)力學(xué)——有限體積法的應(yīng)用》譚維炎 著

推薦代碼，包括了下一篇的淺水波方程。每一份都是經(jīng)過精心挑選，簡短而又清晰的適合初學(xué)者的代碼~

總結(jié)

以上是生活随笔為你收集整理的flux unity 流体_【游戏流体力学基础及Unity代码（八）】激波捕捉法和RiemannSolver...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。