間斷Galerkin方法中的加權本質無振蕩限制器述評

邱建賢,朱 君

(1.廈門大學數學科學學院，福建 廈門 361005;2.南京航空航天大學理學院，江蘇 南京 211106)

間斷Galerkin(DG)方法具有堅實的數學理論和實現方法，易于處理各類復雜區域和邊界問題，可實現高效并行和各類網格，非常適合于計算流體力學(CFD)、計算氣動聲學、計算電磁學等相關工程問題的計算，DG方法已成為目前CFD高精度數值方法的研究熱點之一，也是下一代CFD解算器高精度仿真的首要選擇.非線性雙曲守恒律方程作為流體力學中的重要方程，即使初始值光滑，其解也可能產生間斷，如激波等.而對于存在強間斷解的數值模擬中，DG方法容易產生較為明顯的非物理振蕩，進而產生非線性不穩定，導致數值解爆炸.所以對于此類問題往往需要使用非線性限制器來抑制非物理振蕩，從而實現相關物理過程的準確、穩健的數值模擬.而常用的限制器(如：總變差減小(TVD)限制器、總變差有界(TVB)限制器等)雖然可以控制間斷附近的偽振蕩，但是它們不僅會過渡抹平激波(包含需要人為調整的參數)，而且會降低DG方法的預期設計精度.另一方面，由于限制器的引入，修改了原DG的解，會造成解在局部不匹配，從而影響定常問題的收斂性.因此，近十年來，本課題組一直致力于DG緊致格式的高精度非線性限制器的研究，提出并發展了一系列兼具高精度、強穩健及緊致特性的限制器，在促進基礎算法發展的同時已廣泛應用于相應工程領域.

間斷有限元方法最早由Reed等[1]為解決中子輸運方程(線性雙曲方程)而提出的，它的一個很重要的發展是Cockburn等[2-6]構建了發展型的非線性雙曲守恒律的Runge-Kutta DG(RKDG)方法.這種DG方法[7-13]屬于線方法范疇，在時間離散系統中采用非線性穩定的高精度Runge-Kutta方法，在空間的離散系統中，DG有限元方法，界面間的數值通量采用真正的或逼近的Riemann解，并且使用TVB的限制器，得到了即使在強激波情況下也無偽振蕩的性質.此外，Dumbser等[14-15]在每個時間步使用重構算子來增加高階DG方法的數值精度和穩定性.文獻[16-20]提出了拉格朗日DG方法.Tia等[21]應用泰勒基構造了DG譜有限元方法.Hex等[22]設計了 一種新的加權RKDG方法用于求解三維聲波和彈性波，并為解擴散方程提出了rDG方法[23-24].由于DG方法優越的計算性能，該方法也被廣泛用于流體力學等科學計算領域[25-31]，其他的高階經典DG方法參見文獻[32-34].由于有限元的特性，DG方法和古典的有限差分和有限體積法相比有以下的優點：1)DG精度的階僅依賴于精確解.DG的階可通過適當選取逼近解的多項式的次數得到.2)DG方法具有高度的可并行性.3)單元劃分及對邊界處理的靈活性高，適合于處理復雜的區域的問題.4)DG方法易于進行自適應處理，自適應算法在具有間斷的雙曲問題的求解中非常重要.目前，DG方法已廣泛應用流體動力學、湍流、天氣預報、粒子流、海洋學、油氣儲藏模擬、淺水模型、多孔介質中流體輸運、半導體的模擬、電磁場、圖象處理等問題中，是目前CFD高精度數值方法的主流算法之一.

在DG方法的構造中，限制器是一個很重要的組成部分，它可以控制格式在間斷問題(如激波)計算中產生偽振蕩，且必須采用限制器，保證格式的穩定性.通常限制器的過程可以分為兩個部分：首先確定“壞單元”(“Troubled-cell”)，即解含有間斷的單元，在該單元需要做限制過程；然后在“壞單元”中修正DG的多項式解，由于守恒的要求，需要保證單元平均值不變，同時保證與原DG的解具有相同的精度，且振蕩減小.

在第一部分中，使用“壞單元”或稱間斷指示器，如基于最小模型指示器[4]，基于力矩的指示器[35]，以及改進的力矩指示器[36]，基于DG超收斂性質的KXRCF指示器[37]，基于Harten子單元分辨方法的指示器[38]等.

在第二部分中，包括經典的最小模型限制器[2-4,6]，基于力矩的限制器[37]，以及改進的力矩限制器[36]等.這些限制器屬于斜率型限制器的范疇，它們的優點是可以在強間斷附近有效抑制偽振蕩的出現，但付出的代價是在解的光滑極值點處有可能降低格式的數值精度.另一種類型的限制器是基于本質無振蕩(ENO)和加權ENO(WENO)方法[39-46]的思想構造的，具有很強的激波穿透能力和保持格式的高精度.這種類型的限制器通過將問題單元上的DG方法數值解的自由度進行重構，可以在解的光滑區域實現高精度并在強間斷附近保持本質無振蕩的性質.這些高精度WENO限制器基本上均為根據有限體積WENO格式的構造方法設計的，隨著DG方法格式精度的提高該限制器需要更寬的空間模板，有時會破壞原始DG方法本身具有的空間模板緊湊的性質.此外，眾多的WENO限制器[47-54]，包括新型WENO限制器[51,55-57]，自適應階WENO限制器[58]，中心型WENO(CWENO)限制器[59]，和HermiteWENO(HWENO)限制器[48,50,60]等都是第二類型的限制器.此外，由于CWENO格式[61-64]比經典的WENO格式[65-67]在數值計算上的代價更低，它們可以作為DG格式的后驗子單元限制器[30].

基于上述研究，本課題組近十年來一直致力于DG方法的高精度WENO限制器的研究工作，提出并發展了一系列兼具高精度、強穩健及緊致特性的限制器，在促進基礎算法發展的同時已廣泛應用于相應工程領域.主要包含如下工作：針對已有限制器不能保持DG方法高精度或需要選取參數的問題，結合DG方法和有限體積WENO格式的優良特性設計了WENO限制器，有效地解決了限制器中含有參數和在部分光滑區域中降低DG方法數值精度的問題，實現了對非線性雙曲守恒律的數值模擬，并通過數值結果驗證了該類WENO限制器的一致高精度、無振蕩和高分辨率等特性，解決了經典WENO限制器的模板比DG方法模板更寬的問題；將Hermite插值理論用于限制器構造過程設計了HWENO限制器，有效地緩解了WENO限制器寬模板問題，進而更好地發揮DG方法的作用，實現了對結構網格和非結構網格上非線性雙曲守恒律的數值模擬；為了更好地模擬具有高頻振蕩的物理問題，設計了一類基于三角基函數空間的TWENO格式，并將其發展為DG方法的限制器.該限制器在相同條件下，能夠更好地數值求解復雜波形或高頻振蕩問題；為了減少傳統WENO限制器構造過程過于復雜繁瑣以及解決寬模板問題，設計了一類緊湊簡單的新型HWENO限制器，該限制器解決了經典WENO限制器的構造破壞了DG方法緊致性的弊端，且避免了線性權在復雜網格體系下的計算，從而大大提高了計算效率.

1 DG方法簡述

首先考慮一維標量守恒律方程

(1)

(2)

(3)

(4)

(5)

(6)

2 高精度限制器

2.1 TVB問題單元間斷指示器[4]

定義

從式(2)可以看出

這是修改的標準minmod限制器[70]

其中m定義為

(7)

或者定義為TVB修正的minmod函數[71]

(8)

使用上述TVB問題單元指示器識別問題單元.其中M>0是常數且M的選擇取決于問題的解.對于標量情形，可以通過文獻[4]中的初始條件估計M(M與初始條件在光滑極值處的二階導數的絕對值成正比)；然而，在方程組情形下較難估計M.如果M選擇得太小，DG方法的數值精度可能在解的光滑極值點處降低；如果M選擇得太大，DG方法會在解的非光滑區域產生偽振蕩.由于下面介紹的高精度WENO限制器可以在問題單元上對DG方法數值解中的自由度進行重構時保證格式的高階精度不被破壞，所以是否選擇一個精確的M就顯得不太重要.

2.2 KXRCF間斷指示器[39]

(9)

2.3 WENO限制器

2.3.1 重構步驟1

在高斯或高斯-洛巴托積分點上重構u的點值.對于Pk的DG方法(k+1階精度)，需要使用至少達到2k+2階精度的高斯或高斯-洛巴托求積公式，WENO重構的數值精度至少必須達到2k+1階.為此，本課題組需要用2k+1個相鄰單元Ii-k，…，Ii+k的單元平均值構造多項式并在高斯或高斯-洛巴托求積節點上得到u的高階逼近值.然后執行如下重構步驟[41-42,72]:

2)計算線性權γ0,…，γk.在不同的高斯或高斯-洛巴托求積節點xG上滿足

對應的一組線性權為

對應的一組線性權為:

3)計算光滑指示器βj，用于衡量問題單元Ii中代數多項式pj(x)的光滑程度.光滑指示器βj越小，問題單元中的代數多項式pj(x)越光滑[41,72]:

(10)

光滑指示器βj常被寫成u的單元平均值的二次方形式，詳情參見[41,43,72].

4)根據光滑指示器計算非線性權

(11)

其中γj是上述步驟2)中確定的線性權，ε取一個避免分母為零的小正數.最終的近似值為

(12)

2.3.2 重構步驟2

基于2.3.1重構的點值u(xG)和數值積分獲得自由度的重構值.

對于方程組情形，WENO限制器總是與局部特征分解一起使用，詳情參見文獻[43].三角形和四面體網格上使用DG方法結合WENO限制器求解雙曲守恒律方程組的詳細過程如下[57,73]：首先在非結構網格上使用TVB問題單元指示器識別出問題單元，通過使用相鄰三角形(四面體)的單元平均值進行WENO重構，在保持問題單元上物理量的單元平均值不變的情況下得到重構多項式并對其上DG方法數值解的自由度進行修正.

2.4 HWENO限制器

在本小節中，本課題組將HWENO重構過程應用于2.3.2節.該重構過程的優點是可在保持DG方法相同數值精度的同時有效減少所用空間模板的個數和降低模板的空間尺度，使之具有更好的空間緊湊性.在此，基于一維標量守恒律方程(1)來介紹三階精度(k=2)HWENO限制器的構造情況[73].

1)給定小模版S0={Ii-1,Ii}，S1={Ii，Ii+1}，S2={Ii-1,Ii,Ii+1}和大模板T={S0,S1}，構造Hermite二次重構多項式p0(x),p1(x),p2(x)和四次重構多項式q(x):

得到：

2)計算線性權γ0,γ1和γ2.由

得到

3)根據式(10)計算光滑指示器βj，然后根據式(11)計算非線性權.重構多項式的第一個自由度為

(13)

1)給定小模版S0={Ii-1,Ii}，S1={Ii,Ii+1}，S2={Ii-1,Ii,Ii+1}和大模板T={S0，S1}，構造Hermite三次重構多項式p0(x)，p1(x)，p2(x)和五次重構多項式q(x):

得到：

2)計算線性權γ0,γ1和γ2:

得到

3)根據式(10)計算光滑指示器βj，然后根據式(11)計算非線性權.重構多項式的第二個自由度為

(14)

本小節只考慮了一維HWENO限制器的構造流程，該限制器可以推廣到多維情形.文獻[74]詳細地介紹了使用DG方法結合HWENO限制器在非結構網格上求解三維非線性雙曲守恒律方程組.同樣該文首先在三維四面體網格上使用TVB問題單元指示器找出問題單元，然后通過使用相鄰四面體的單元平均值或導數平均值，在保持問題單元上物理量的單元平均值不變的情況下進行HWENO重構，得到問題單元內的多項式并對其上DG方法數值解的自由度進行修正.

2.5 TWENO限制器

……

(15)

(16)

和

(17)

γ1=(cos(h(-1+α))-cos(h(2+α))+

cos(h(2+3α))-cos(h+3haα)+4hsin(hα)+

4hαsin(hα)-2hsin(2hα)-2hαsin(2hα)-

4hαsin(h(1+α))+2hαsin(2h(1+α)))/

γ2=1-γ1-γ3,

γ3=((-cos(h(-1+α))+cos(h(2+α))-

cos(h(2+3α))+cos(h+3hα)+2hsin(h)+

4hαsin(h)+hsin(h(-1+α))+hαsin(h(-1+

α))-hαsin(h(2+α))-hαsin(h(2+3α))-

2hsin(h+2hα)+hsin(h+3hα)+hαsin(h+

γ1=(-cos(h(-1+α))+cos(h(2+α))-

cos(h(2+3α))+cos(h+3hα)+2hsin(h)+

4hαsin(h)+hsin(h(-1+α))+hαsin(h(-1+

α))-hαsin(h(2+α))-hαsin(h(2+3α))-

2hsin(h+2hα)+hsin(h+3hα)+hαsin(h+

γ2=1-γ1-γ3,

γ3=(cos(h(-1+α))-cos(h(2+α))+

cos(h(2+3α))-cos(h+3hα)+4hsin(hα)+

4hαsin(hα)-2hsin(2hα)-2hαsin(2hα)-

4hαsin(h(1+α))+2hαsin(2h(1+α)))/

3)根據式(10)計算光滑指示器：

4hα(ui-1-ui-2)(ui-1-ui)cos(hα)+2(ui-1-

ui-2)(ui-1-ui)sin(hα)-(ui-2-ui)(-2ui-1+

ui-2+ui)sin(2hα)-2(ui-1-ui-2)(ui-1-

ui)sin(3hα)-(ui-1-ui)2sin(4hα)+

4hα(ui-1-ui-2)(ui-1-ui)cos(hα)-2(ui-1-

ui-2)(ui-1-ui)sin(hα)+(ui-2-ui)(-2ui-1+

ui-2+ui)sin(2hα)+2(ui-1-ui-2)(ui-1-

ui)sin(3hα)+(ui-1-ui)2sin(4hα)))/256,

cos(hα)-4(-1+h2α2)(ui-1-ui)(ui-ui+1)

sin(2hα)))/256,

(ui+1-ui+2)cos(hα)+2(-1+h2α2)(ui-ui+1)

2h2α2uiui+1sin(2hα)+2ui+1ui+2sin(2hα)-

2h2α2uiui+2sin(3hα)+2ui+1ui+2sin(3hα)-

4)根據式(11)計算非線性權叫ωn,n=1,2,3.并通過高斯-洛巴托求積點xG處的重構點值u(xG)和數值積分得到

(18)

本小節介紹了三角函數多項式空間中TWENO限制器的構造過程.其構造思想與前面類似，首先使用TVB問題單元指示器找到問題單元，然后使用相鄰單元的單元平均值通過TWENO限制器在問題單元內重構三角函數多項式.文獻[74]中的數值結果表明，當基于三角函數多項式空間的格式用于模擬波類和高頻振蕩問題時能得到更好的數值結果.

2.6 簡單緊湊型HWENO限制器

首先使用KXRCF問題單元指示器[37]來檢測問題單元，接著在一維標量情況下闡述簡單緊湊型HWENO限制器的構造過程.具體過程如下所示：

(19)

(20)

(21)

2)選擇線性權γ0,γ1，γ2.如文獻[51]所述，本課題組在3個單元Ii-1，Ii+1，Ii中使用了3個多項式p0(x)，p1(x)，p2(x)的完整信息，因此線性權只需任意選擇和為1的正數即可保持高階數值逼近.

3)參照式(10)計算光滑指示器，參照式(11)計算非線性權.

二維結構網格和三角形網格上簡單緊湊型HWENO限制器的構造過程參見文獻[54，75].該簡單緊湊型HWENO限制器的主要創新點是重構過程僅使用來自問題單元及其近鄰單元的DG方法數值解的信息，空間模板非常緊湊.在重構過程中使用的線性權不再需要根據計算網格的空間拓撲結構通過繁冗的數值計算獲得，可以是任意和為1的正數.因此，該限制器能在結構網格、非結構網格、雜交網格、運動網格、自適應網格等情況下靈活使用且對網格質量要求較低.

2.7 Multi-resolution WENO限制器

本小節借鑒multi-resolution方法[76-86]構造了一種新型multi-resolution WENO限制器[87-88].本課題組仍在結構網格上針對一維標量雙曲守恒律方程(1)構造multi-resolution WENO限制器[89].具體過程如下所示：

(22)

(23)

(24)

4)計算非線性權.根據文獻[94-95]中提出的類似思想，定義：

(25)

非線性權為

(26)

(27)

本小節只介紹了一維結構網格上multi-resolution WENO限制器的構造過程，二維結構網格和三角形網格上multi-resolution WENO限制器的構造過程參見文獻[89,96].DG方法結合這種multi-resolution WENO限制器在求解雙曲守恒律方程組時，可以在光滑區域保持DG方法的數值精度且在強間斷附近有效抑制偽振蕩的出現.這種新型multi-resolution WENO限制器使用的線性權可以是任意和為1的正數.這種新的WENO限制器構造非常簡單且適用于任意高階DG方法及應用于二維和三維問題的高精度數值求解.

3 總結及展望

本文綜述了高精度DG方法和高精度WENO限制器的一般構造過程及其在雙曲守恒律中的應用.與現有高精度限制器相比，新型WENO限制器在以下幾個方面獨具創新：可以在多維(一維、二維以及三維)情況下保持物理量的守恒性和計算格式的魯棒性，在解的光滑區域不會因限制器的原因導致計算精度下降，同時在強激波或接觸間斷區域保持本質無振蕩性質；在解的光滑區域依然保持高階數值精度的同時間斷過渡區域更為狹窄，能顯著提高激波分辨率；很好地保持了DG方法的空間緊致特性，對DG方法在求解定常可壓縮流體問題的收斂性和收斂速度有較大改進；易于在非結構網格、混合網格、運動網格、自適應網格等復雜網格體系下構造和編程實現，隨著計算問題的維數增加，限制器的構造和編程依然非常簡單，顯著降低計算機內存資源的占用進而顯著提高計算效率.截至目前，新型WENO限制器是已發表相關文獻中空間最緊致、魯棒性最高、編程實現最為簡便的一類高精度非線性限制器.該類限制器具有較強的通用性，能在結構網格、非結構網格、雜交網格、運動網格、自適應網格等情況下靈活使用且對網格質量要求較低，且能夠推廣至任一緊致類格式，如重構校正方法(CPR)/通量重構(FR)格式、譜差分(SD)/譜體積(SV)方法等，因此具有較廣泛的軍事和商業工程應用前景.