基于映射函數的新型五階WENO格式

劉博，李詩堯，陳嘉禹，程啟豪，時曉天,*

1.中國航天空氣動力技術研究院，北京 100074

2.中國科學院 力學研究所 高溫氣體動力學國家重點實驗室，北京 100190

3.中國科學院大學 工程科學學院，北京 100049

4.天津大學 水利工程仿真與安全國家重點實驗室，天津 300072

5.天津大學 建筑工程學院，天津 300350

6.天津大學 數學學院，天津 300350

隨著航空航天技術的發展，計算流體力學(CFD)逐漸成為研究各類復雜氣體流動問題的主要手段，而計算機科技的進步為數值模擬提供了更廣闊的前景。針對有強非線性效應的物理問題(諸如激波問題)，使用擁有良好頻譜性的線性格式計算會出現非物理振蕩，嚴重時甚至會使計算發散，即使加密網格也不會減弱振蕩[1]，為解決此類問題需要發展非線性格式。TVD(Total Variation Diminishing)格式[2]具有良好的分辨率和穩定性,但是精度較低。DG(Discontinous Galerkin)[3]、FR(Flux Reconstruction)[4]等有限元方法雖然可以適應非結構網格，但在捕捉間斷方面仍有待完善。Harten[5]和Shu[6]等把TVD格式要求降低為滿足總變差有界TVB(Total Variation Bounded)的條件，提出了ENO(Essentially Non-Oscillatory scheme)格式。在此基礎上，Liu等[7]首次提出WENO(Weighted Essentially Non-Oscillatory scheme)思想，Jiang和Shu[8]隨后提出光滑因子和加權系數具體的表達式，并提出了三階和五階的WENO格式，此格式適合處理含有間斷的問題，被廣泛使用[9-14]。然而WENO-JS格式在極值點處的精度會降低，Henrick等[15]分析了其原因并通過引入映射函數重構加權系數，提出了WENO-M格式，Borges[16]和Castro[17]等對光滑因子進行改進，形成WENO-Z格式，克服了極值點降階的缺陷。隨后有研究者提出一系列針對WENO-Z的改進方案[18-20]，但仍然改進不了耗散較大的缺陷。Wu等[21]通過數學推導指出WENO-Z格式在高階極值點處會降階。為了減少耗散，Fan[22]提出WENO-η格式，Ha等[23]提出WENO-NS格式。另有研究人員基于Henrick的思想，構造新的映射函數，例如Feng等[24]將WENO-M格式推廣到更高階情形，提出WENO-IM格式，隨后Feng等[25]又提出在端點處應接近ENO性質，改進得到WENO-PM格式。Li[26]和劉朋欣[27]等則認為在端點處應逼近恒等映射，提出WENO-PPM格式。Wang等[28]通過大量實驗發現映射函數在右端點處的導數影響微小，從而改進出振蕩更小的WENO-RM格式。Vevek等[29]提出基于變化因子λ的WENO-AIM格式。

基于映射函數思想重構加權因子，大多數已有的映射函數都是多項式型函數，為了提高函數在端點處的收斂速度與理想權重值處穩定性，以指數函數作為基函數構造分段映射函數WENO-Pe(Piecewise exponential mapping function)，使其匹配常用WENO格式，并以五階為例通過若干經典算例驗證，與其他格式對比來驗證新格式的性能。

1 問題描述與數值離散

1.1 標量方程

一維對流方程是最簡單的雙曲線偏微分方程，該方程可分為線性與非線性兩類：

(1)

(2)

1.2 矢量方程

無黏性流體動力學中最重要的基本方程，即Euler方程組，是指對無黏性流體微團應用牛頓第二定律得到的運動微分方程。以一維守恒的情形為例:

(3)

式中：

(4)

其中：ρ、u和p分別代表流體密度、速度與壓強；E為單位體積流體的總能量，其表達式為

(5)

式中：γ為氣體常數。

1.3 數值離散

在對Euler方程組離散時，采取Steger-Warming流通矢量分裂法，通量F分解為正負兩部分分別參與運算：

F=F++F-

(6)

(7)

再對時間導數項進行離散，為了與空間項精度匹配，所用的五階格式均采用三階精度的Runge-Kutta格式[30]：

(8)

2 通量重構

2.1 WENO格式思想

(9)

該函數是通過每個子模板上使用Lagrange插值多項式得到代數多項式，正負通量模板如圖1(a)和圖1(b)所示。

圖1 WENO5總模板與子模板

(10)

(11)

并要求非線性權在光滑區域充分接近理想權重，在間斷附近退化為ENO格式，以保證減小振蕩，因此構造合適的非線性權成為關鍵問題。

2.2 各類非線性權

2.2.1 WENO-JS格式

(12)

并提出一種非線性權構造方案：

(13)

式中:p為≥2的正整數;ε為防止分母為0的充分小正數，通常取10-6。

2.2.2 WENO-M格式

(14)

再將新得到的函數值重構為權重系數：

(15)

注意到該映射函數在(0,1)區間內嚴格單調遞增且具有二階連續導數，同時還應當滿足如下性質：

(16)

通過將ωM,k表達式在Ck點做Taylor展開可以證明：

ωM,k=Ck+O(Δx3)

(17)

2.2.3 WENO-PPM格式

Li等[26]認為映射函數應經過映射以后的非線性權在定義域邊界處應盡快地收斂到0或者1，于是提出分段多項式函數作為映射函數：

(18)

其重構加權因子為

(19)

2.2.4 WENO-PM格式

Feng等[25]注意到，Henrick等提出的映射函數(14)在端點處導數值偏大:

(20)

這樣會使靠近端點處的權重被放大，出現較大誤差，為了修正此問題，Feng等在式(16)基礎上對映射函數補充了新要求，即端點處函數的趨近速度：

g′(0)=g′(1)=0

(21)

并設計出新的映射函數：

gPM(ω)=

(22)

使用五階WENO格式時推薦參數n=4，重構權重系數為

(23)

2.2.5 WENO-RM格式

Wang等[28]認為分段函數的全局光滑性不足，應使用全局高階連續函數克服來自不光滑模板的影響，故提出新的映射函數：

(24)

(25)

函數(24)仍然具備如下性質：

(26)

2.2.6 WENO-Pe格式

總結以上的研究，一個好的映射函數應具有如下特點：

1) 在(0,1)區間內具有良好的光滑性。

2) 在各階極值點處精度不下降。

3) 有較強的克服不光滑模板的能力。

4) 在g(0)處平緩收斂到0。

5) 在g(1)處接近恒等映射。

多項式函數在固定區間的增長速率依賴次數與系數，故考慮使用其他類型函數替代，而指數函數具有無窮次可導且增長率較大的良好性質，可考慮使用其構造映射函數，使該函數具有如下性質：

(27)

2.3 構造WENO-Pe格式的非線性權

以Ck為分界點，將[0,1]區間分為[0,Ck)和[Ck,1]，在左半端內設

(28)

(29)

在右半段內設

(30)

αR,mR+1tRn+mR+1)+Ck

(31)

式中：AL和AR為正實數；mL、mR、n是正整數(n≥m≥2)；αL,i和αR,i是由AL、AR、mL、mR和n確定的系數，對于五階WENO而言，取n=5，mL=mR=2(一般情況下取mL=mR，以下簡稱m)，并且AL=AR即可。易證

(32)

為使在左端點處表現為ENO性質，需滿足

(33)

則αL,i應滿足：

(34)

同樣，為使右端點滿足WENO性質，需滿足：

(35)

αR,i應滿足

(36)

式中：

(37)

因此整體有

(38)

(39)

圖2是AL=AR=15時，五階WENO-Pe格式關于3個理想權重的映射函數圖像。

圖2 WENO5-Pe的映射函數

圖3是映射函數gPe同其他映射函數(按其文獻中推薦的最佳參數)圖像對比，其中各函數的具體參數如下：

(40)

圖3 WENO5格式中第2個模板的Pe格式的映射函數與M、IM、PPM、RM和PM函數的對比

如果應用于(2r-1)階WENO格式，可取n=r+2,m=r-1，相應的系數為

P(n+m;n,i,m-i)

(41)

式中：i=0,1,…,m；αR,0和αR,mR+1表達式同式(36)；P(n;q1,q2,…,qk)為多重全排列：

(42)

3 格式性能分析

3.1 截斷誤差

(43)

將光滑函數IS在xj處做Taylor展開

(44)

Henrick等[16]指出在f的一階極值點x0處，即x→x0時有

f′(x)～O(Δx)

(45)

(46)

根據映射函數重構的權重系數，在ω=Ck處進行Taylor展開有

(47)

即有

(48)

只要n取不小于2的正整數, 該格式就能在極值點處達到空間上的五階精度。

3.2 頻譜性

WENO是非線性格式，其頻譜性分析與線性格式Fourier頻譜分析類似，可以使用Pirozzol提出的近似法色散關系(ADR)。令fj=eikxj，其中i為虛數單位，則

fj+n=eikxj+n=fjeikxn

(49)

由此可以推出修正波數：

K=Re(k)+iIm(k)

(50)

具體可參照文獻[32-33]。其中,實部Re(k)代表色散誤差，虛部Im(k)代表耗散誤差。圖4是格式WENO-Pe與其他格式的頻譜性對比，系數同式(34)和式(35)，其中WENO-AIM系數選擇k=4，m=2，C=104。

從圖4可見，WENO-Pe格式的色散誤差與數值耗散性均小于其他格式。

4 數值實驗

無特殊說明，本文的算例均采用無量綱形式的方程。

4.1 精度驗證

4.1.1 定常問題

下面驗證WENO5-Pe格式捕捉高階極值點的性能，仿照文獻[34]，構造如下函數族：

u0(x)=ea(x+1)xn+1

(51)

式中：a= 0.5；n為非負整數，當n=0時，該函數在[-1,1]上無極值點；當n= 1時，在此區間上只有唯一的一階極值點0；當n> 1時，僅有唯一的一階、二階極值點0。設N為網格數，分別使用WENO5-Pe、WENO5-PPM、WENO5-RM和WENO5-IM格式，具體格式參數同式(40)，在不同網格下計算其一階導數，選取n= 0、1、2時的函數，統計數值解在x= 0處與真實值的誤差并計算誤差階，結果如表1所示，關于WENO5-Z在極值點處降階可參考文獻[21]。由表1可見，WENO5-PPM格式在無極值點或僅有一階極值點情況下能達到五階，但有二階極值點時，明顯下降到三階左右；WENO5-RM和WENO5-IM在一階極值點處就會降階到三階，有二階極值點時甚至降到更低的二階；而WENO5-Pe格式無論有無二階極值點均能保持在五階左右精度，且在相同網格下，WENO5-Pe格式的誤差更小，這表明提出的格式具備更高的精度。

表1 4種數值格式計算結果對比(定常問題)

4.1.2 非定常問題

本節驗證一維線性對流方程，該算例選自文獻[15]：

(52)

此問題存在精確解：

(53)

顯然，這是一個周期為2的函數，且存在2個一階極值點和一個二階極值點，無三階極值點。選取時間步長為Δt=h5/3(匹配精度)。分別使用WENO5-Pe、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Z格式計算到時間t=2.2和t=3(均無量綱，下文同)，統計L2誤差及精度階，如表2所示。

(54)

由表2可以看出，WENO5-Pe、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Z格式均能達到預期精度，但WENO5-Pe的精度要優于其他4種格式。

4.2 耗散驗證

為了驗證WENO5-Pe格式長時間計算的穩定性，選取如下一維對流方程：

(55)

式中:取C=10 000，a=-0.5;分別使用WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式計算到時間t=100，計算域為[-1,1]，選取網格N=100、N=200、CFL=0.5，得到計算結果如圖5所示。

由圖5可見，使用同階格式時，WENO5-PPM與WENO5-RM效果相似，其耗散均小于WENO5-IM，WENO5-IM在網格加密后會出現較大非物理振蕩，而WENO5-Pe更接近于真實值，比其他4個格式的耗散更小，能夠更精確地捕捉高階極值點。

表2 5種數值格式計算的結果對比(非定常問題)

圖5 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式長時間計算的穩定性

4.3 Sod激波管

該算例的初始條件為[35]

(56)

圖6 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式模擬Sod激波管

設置兩邊界為流邊界，網格數J=200，計算時間為t=0.18，CFL=0.4, 并將WENO5-JS在2 000個網格下的計算結果作為準精確解，分別使用WENO5-Z、WENO5-PPM、WENO5-IM、WENO5-RM和WENO5-Pe格式，得到密度曲線如圖6所示。由文獻[16,24,26,28]的數值實驗可知，WENO5-Z、WENO5-PPM、WENO5-IM和WENO5-RM均比WENO5-Z計算效果好，而從圖6可見，WENO5-Pe的捕捉間斷能力相較于以上格式有明顯提升。

4.4 Lax激波管

此算例[36]描述的流場中同時包含有激波、接觸間斷、膨脹波和平滑區域，其初始條件為

(ρ,u,p)=

(57)

網格數J=200、CFL=0.2，計算終止到t=0.16，并將WENO5-JS在2 000個網格下的計算結果作為準精確解，分別使用WENO-Z、WENO5-PPM、WENO5-IM、 WENO5-RM和WENO5-Pe格式，得到密度曲線如圖7所示。由圖7(a)～圖7(c)可見，WENO5-Pe更接近于真實解，產生的非物理振蕩更小，在間斷處的識別更加敏感，魯棒性更強。

圖7 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式模擬Lax激波管

4.5 Osher-Shu問題

該問題[37]描述一道馬赫數3 的右行激波與熵波的相互作用，熵波在激波作用下被壓縮，并向下游傳播方向生成一系列聲波。此算例可以檢驗數值格式的間斷識別能力、穩定性和分辨率，其初始條件為

(58)

選取網格點J= 240、CFL = 0.4，計算終止時間為t= 1.8，分別使用WENO5-Z、WENO5-PPM、WENO5-IM、WENO5-RM和WENO5-Pe格式計算該問題, 選取WENO5-JS格式在2 000個網格下計算結果為準精確解，得到密度曲線如圖8所示，統計計算所用時間如表3所示。

由圖8可見，在同樣網格下WENO5-Pe格式具備更好的分辨率，而且比其他4個格式計算出的函數極值點偏移量明顯的小，而WENO5-Z的分辨率要略差于其他4種基于函數映射的WENO格式，說明WENO5-Pe的色散誤差更小。

圖8 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式模擬Shu-Osher問題

表3 指定網格下各格式計算耗時

另一方面，使用網格N= 360的WENO5-Z格式，與N= 240的WENO5-Pe格式計算結果對比，可見在CPU耗時接近的情況下，240網格數的WENO5-Pe計算效果依舊優于360網格數的WENO5-Z計算效果。

從以上一維的算例中不難看出，WENO5-IM與WENO5-PPM在精度與分辨率的表現上都優于WENO5--RM和WENO5-Z，而WENO5-Pe對極值點捕捉位置、識別間斷的敏感性、數值精度以及計算的魯棒性都要優于其他4個格式。

4.6 二維格式性能驗證

二維歐拉方程的控制方程為

(59)

式中：

(60)

通量F與G的離散與分裂形式與一維情形一致。下面分別選取對流問題與定常問題分別驗證格式精度與收斂性。

4.6.1 精度驗證

二維對流問題：

(61)

該問題存在精確解：

u(x,y,t)=sin(x+y-2t)+1.200

(62)

且不需要通量分裂，可以更好地檢驗格式精度，分別使用WENO5-Z、WENO5-PPM、WENO5-IM、WENO5-RM和WENO5-Pe在不同網格下計算到t= 0.2和t= 1.2時刻，統計L2誤差并計算誤差階，結果如表4所示。通過表4可見，5種格式在達到理論精度階的同時，WENO5-Pe的精度要略優于WENO5-PPM、WENO5-IM和WENO5-RM，而它們均優于WENO5-Z。

4.6.2 平板激波反射問題

有一道與x方向成29°角的激波射向平面上一壁面，并與壁面碰撞進而發生反射，形成入射波和反射波[38]，如圖9(a)所示。當流動達到穩定后，該問題存在精確解，如圖9(b)所示，該問題屬于無黏流動，可以用Euler方程描述該過程。該問題的計算域為[0,4]×[0,1]的矩形區域，當t= 0時，使激波剛好與底部壁面相遇，具體初始條件可以根據激波前后密度、壓強等相關公式計算，其左側初始條件：

(63)

頂部初始條件：

[1.699 97,2.619 34,-0.506 32,1.528 19]

(64)

表4 5種數值格式在二維算例下精度對比

圖9 二維平板激波反射示意圖

式中：右側邊界為自由流出條件，底部邊界為反射條件。分別使用WENO5-Z、WENO5-PPM、WENO5-IM、WENO5-RM和WENO5-Pe在200×50的網格下計算該流動，終止時間為t=4.5，CFL=0.5。將計算獲得的結果與精確解對比，結果如圖10所示，并統計L2誤差如表5所示。

另外，比較相同網格下不同格式的L2誤差隨迭代步數發展的關系曲線如圖10(e)和圖10(f)所示，從而對比其收斂速度，計算網格仍然為200×50。在相同網格，相同物理時間下，從圖10(a)～圖10(c)可見，WENO5-Pe的分辨率更高，穩定性更好，捕捉激波的位置更準確，非物理振蕩更小，而WENO5-Z對激波初始位置的捕捉誤差較大，WENO5-IM出現了較明顯的非物理振蕩；通過表5可見，計算到相同時刻WENO5-Pe的精度高于其他五階格式；通過圖10(e)和圖10(f)可見，相同網格和耗時的情況下，WENO5-Pe相較于其他格式計算到收斂的耗時更短，而最終收斂時的精度也更高。

圖10 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式模擬平板激波反射問題

表5 5種WENO格式計算二維激波反射到t=4.5時刻得到的結果對比

4.7 二維Riemann問題

該問題選自Liu等的數值實驗[39]，初始時刻4個計算子區域內流體有不同的初值，當瞬間除去4個計算子區域內膜后，在計算區域形成激波、渦和接觸間斷相互作用的復雜流動。該問題求解域為[0,2]×[0,2]，初始條件為

(65)

4條邊界均為自由輸出邊界，即在每一條邊界上，都有

(66)

使用400×400均勻網格，計算終止時間t=0.8,條件數CFL=0.5,分別使用五階的WENO5-IM、WENO5-PPM、WENO5-RM、WENO5-Z和WENO5-Pe計算其流場變化。圖11(a)～圖11(e)是從0.2～1.7的30條密度等值線。

觀察圖11(a)～圖11(e)可見，WENO5-RM和WENO5-PPM比WENO5-Z的分辨率更高，同時WENO5-IM比WENO5-RM和WENO5-PPM的分辨率稍高，而WENO5-Pe比WENO5-IM的分辨率更明顯，可以觀察到更多的渦結構細節，密度的分界線輪廓也更精準，由此可見其捕捉流場細節能力更強。

4.8 Rayleigh-Taylor不穩定性問題

該算例描述初始流場中存在兩種密度不同的流體，在重力場作用下，位于上方的密度較大的流體加速進入下方的密度較小的流體的失穩過程，最終形成復雜的流場[40]。其計算域為[0,0.25]×[0,1]，計算初始條件為在控制方程式(59)和式(60)右側添加源項以模擬重力影響。

(67)

S=[0,0,ρ,ρv]T

(68)

氣體指數取γ=5/3，設置左右邊界為反射壁面，上下邊界設為常值。

(69)

計算網格為240×960，計算終止時間t= 1.95，CFL = 0.5, 圖12(a)～圖12(e)給出的是5種WENO格式的計算結果，密度等值線為0.965 0～2.1481之間15等份。

圖12 WENO5-Z、WENO5-PPM、WENO5-RM、WENO5-IM和WENO5-Pe格式計算Rayleigh-Taylor問題得到的密度等值線

由圖12可見，同樣網格下WENO5-Z分辨率最低，甚至在y= 0.5附無法識別渦結構，WENO5-RM的分辨率略高于WENO5-PPM和WENO5-IM，而WENO5-Pe的分辨率更好，能夠在擴散界面捕捉到更多流場結構。

4.9 雙Mach反射問題

強激波雙Mach反射問題[41]是測試數值格式分辨率的經典算例之一。本算例的計算域為[0,4]×[0,1], 初始條件為一道與x軸正方向成60°，Ma= 10的右行斜激波在x= 1/6處與下壁面相遇，激波一直延伸至上壁面，其右側參數[ρ,u,v,p]為[1.4,0,0,1.0]，左側參數為[8.0,7.144 7,-4.125,116.5]。上邊界條件按激波傳播精確解給出，左邊界與x≤1/6的下邊界按流入條件給出，右邊界為流出條件，x>1/6的下邊界部分為壁面條件。計算網格為960×240，終止時間為t=0.2，CFL = 0.5，圖13(a)～圖13(e)為最終時刻尾部區域5種WENO格式計算的密度2～22的40條等值線圖。

圖13 WENO5-PPM、WENO5-IM、WENO5-RM、WENO5-Z和WENO5-Pe格式計算雙馬赫反射問題得到的密度等值線

由圖13可見，WENO5-Z、WENO5-PPM和WENO5-RM格式的計算結果相似，三者的分辨率略低于WENO5-IM，而WENO5-Pe的分辨率又高于WENO5-IM，其計算結果顯示的激波更清晰，滑移線上卷起的渦結構也更豐富。

5 結 論

針對傳統WENO格式在極值點處精度降低的情況，本文基于對權重系數重構的思想，設計出一族以WENO-JS格式的權重系數為自變量的映射函數，從而得到新格式WENO-Pe，克服了一階和高階極值點處降階的缺陷，且可以推廣到高階格式，并通過數值實驗與其他格式對比驗證了WENO-Pe格式的性能：

1) WENO-Pe格式的色散誤差與數值耗散均小于WENO-JS、WENO-Z、WENO-M以及列舉的其他映射函數型格式。

2) 對比WENO-IM和WENO-RM格式，WENO-Pe在一階極值點處能夠保持理論精度。對比WENO-PPM格式，WENO-Pe在二階極值點處也基本保持精度不下降。

3) 在相同精度情況下，本文WENO-Pe格式擁有更良好的分辨率，在Riemann問題和雙馬赫反射問題中，明顯可見其計算結果能捕捉更多的間斷、渦結構等流場細節。