邱建賢,朱 君

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

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

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

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

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

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

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

1 DG方法簡述

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

(1)

(2)

(3)

(4)

(5)

(6)

2 高精度限制器

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

定義

從式(2)可以看出

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

其中m定義為

(7)

或者定義為TVB修正的minmod函數(shù)[71]

(8)

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

2.2 KXRCF間斷指示器[39]

(9)

2.3 WENO限制器

2.3.1 重構(gòu)步驟1

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

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

對應的一組線性權(quán)為

對應的一組線性權(quán)為:

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

(10)

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

4)根據(jù)光滑指示器計算非線性權(quán)

(11)

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

(12)

2.3.2 重構(gòu)步驟2

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

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

2.4 HWENO限制器

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

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

得到：

2)計算線性權(quán)γ0,γ1和γ2.由

得到

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

(13)

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

得到：

2)計算線性權(quán)γ0,γ1和γ2:

得到

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

(14)

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

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)根據(jù)式(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)根據(jù)式(11)計算非線性權(quán)叫ωn,n=1,2,3.并通過高斯-洛巴托求積點xG處的重構(gòu)點值u(xG)和數(shù)值積分得到

(18)

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

2.6 簡單緊湊型HWENO限制器

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

(19)

(20)

(21)

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

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

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

2.7 Multi-resolution WENO限制器

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

(22)

(23)

(24)

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

(25)

非線性權(quán)為

(26)

(27)

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

3 總結(jié)及展望

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