邱建賢,朱 君

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

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

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

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

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

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

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

1 DG方法簡(jiǎn)述

首先考慮一維標(biāo)量守恒律方程

(1)

(2)

(3)

(4)

(5)

(6)

2 高精度限制器

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

定義

從式(2)可以看出

這是修改的標(biāo)準(zhǔn)minmod限制器[70]

其中m定義為

(7)

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

(8)

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

2.2 KXRCF間斷指示器[39]

(9)

2.3 WENO限制器

2.3.1 重構(gòu)步驟1

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

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

對(duì)應(yīng)的一組線性權(quán)為

對(duì)應(yīng)的一組線性權(quán)為:

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

(10)

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

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

(11)

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

(12)

2.3.2 重構(gòu)步驟2

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

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

2.4 HWENO限制器

在本小節(jié)中，本課題組將HWENO重構(gòu)過程應(yīng)用于2.3.2節(jié).該重構(gòu)過程的優(yōu)點(diǎn)是可在保持DG方法相同數(shù)值精度的同時(shí)有效減少所用空間模板的個(gè)數(shù)和降低模板的空間尺度，使之具有更好的空間緊湊性.在此，基于一維標(biāo)量守恒律方程(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)多項(xiàng)式p0(x),p1(x),p2(x)和四次重構(gòu)多項(xiàng)式q(x):

得到：

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

得到

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

(13)

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

得到：

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

得到

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

(14)

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

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)計(jì)算光滑指示器：

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

(18)

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

2.6 簡(jiǎn)單緊湊型HWENO限制器

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

(19)

(20)

(21)

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

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

二維結(jié)構(gòu)網(wǎng)格和三角形網(wǎng)格上簡(jiǎn)單緊湊型HWENO限制器的構(gòu)造過程參見文獻(xiàn)[54，75].該簡(jiǎn)單緊湊型HWENO限制器的主要?jiǎng)?chuàng)新點(diǎn)是重構(gòu)過程僅使用來自問題單元及其近鄰單元的DG方法數(shù)值解的信息，空間模板非常緊湊.在重構(gòu)過程中使用的線性權(quán)不再需要根據(jù)計(jì)算網(wǎng)格的空間拓?fù)浣Y(jié)構(gòu)通過繁冗的數(shù)值計(jì)算獲得，可以是任意和為1的正數(shù).因此，該限制器能在結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格、雜交網(wǎng)格、運(yùn)動(dòng)網(wǎng)格、自適應(yīng)網(wǎng)格等情況下靈活使用且對(duì)網(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)格上針對(duì)一維標(biāo)量雙曲守恒律方程(1)構(gòu)造multi-resolution WENO限制器[89].具體過程如下所示：

(22)

(23)

(24)

4)計(jì)算非線性權(quán).根據(jù)文獻(xiàn)[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)造過程參見文獻(xiàn)[89,96].DG方法結(jié)合這種multi-resolution WENO限制器在求解雙曲守恒律方程組時(shí)，可以在光滑區(qū)域保持DG方法的數(shù)值精度且在強(qiáng)間斷附近有效抑制偽振蕩的出現(xiàn).這種新型multi-resolution WENO限制器使用的線性權(quán)可以是任意和為1的正數(shù).這種新的WENO限制器構(gòu)造非常簡(jiǎn)單且適用于任意高階DG方法及應(yīng)用于二維和三維問題的高精度數(shù)值求解.

3 總結(jié)及展望

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