柴得林，李雨潤(rùn)，孫中國，席　光

(西安交通大學(xué) 能源與動(dòng)力工程學(xué)院， 陜西 西安　710049)

?

WENO格式權(quán)重分析與改進(jìn)*

柴得林，李雨潤(rùn)，孫中國，席光

(西安交通大學(xué) 能源與動(dòng)力工程學(xué)院， 陜西 西安710049)

為了優(yōu)化WENO格式計(jì)算性能，在對(duì)Jiang和Shu的經(jīng)典WENO格式(記為WENO-JS)加權(quán)方法分析的基礎(chǔ)上，通過引入間接光滑指數(shù)，構(gòu)造出一種新的WENO格式——WENO-E格式，取得減小間斷區(qū)耗散的效果。理論分析表明，該格式與WENO-JS格式計(jì)算效率基本相同，可達(dá)到相同階的計(jì)算精度；但在相同網(wǎng)格下，較之WENO-JS格式，該格式對(duì)光滑區(qū)域的求解有更小的截?cái)嗾`差，對(duì)間斷的捕捉有更高的分辨率。與WENO-JS格式相比，采用WENO-E格式進(jìn)行線性遷移方程、非線性Burgers方程、歐拉方程等相關(guān)問題的數(shù)值實(shí)驗(yàn)，均能取得更好的數(shù)值結(jié)果。

加權(quán)本質(zhì)無震蕩格式；光滑因子；高精度；高分辨率

合理、高效、精確地模擬非定常、含間斷、多尺度復(fù)雜流場(chǎng)是計(jì)算流體動(dòng)力學(xué)研究領(lǐng)域的一大難題。高精度格式對(duì)網(wǎng)格規(guī)模要求低，對(duì)復(fù)雜流動(dòng)分辨率高，計(jì)算精度高，適合計(jì)算包含各類間斷和多種復(fù)雜流動(dòng)光滑結(jié)構(gòu)以及二者相互作用的流場(chǎng)。而加權(quán)本質(zhì)無震蕩格式(Weighted Essentially Non-Oscillatory schemes，WENO)正是研究熱點(diǎn)之一。作為一種性能優(yōu)越的激波捕捉格式，經(jīng)過二十余年的發(fā)展，WENO格式以及各種基于WENO格式思想的改進(jìn)格式已經(jīng)被廣泛應(yīng)用于計(jì)算流體的各個(gè)領(lǐng)域，得到了學(xué)者們的認(rèn)可和關(guān)注。

WENO格式是由Liu，Osher[1]等在本質(zhì)無震蕩格式(Essentially Non-Oscillatory schemes，ENO)構(gòu)造思想基礎(chǔ)上于1994 年提出的。而后，Jiang和Shu于1996年對(duì)WENO 格式作了進(jìn)一步的改進(jìn)，給出了構(gòu)造任意階有限差分WENO格式的一般公式(WENO-JS)，并可以有效地推廣到多維空間問題[2]。但是，WENO格式在極值點(diǎn)處精度會(huì)降低，例如5階精度WENO格式在一階極值點(diǎn)處會(huì)降為3階。因而，對(duì)WENO格式的分析改進(jìn)研究仍備受關(guān)注，在這一方面也已經(jīng)有了大量的研究成果。其中，代表性的有WENO-M，WENO-Z格式等。在對(duì)WENO-Z格式的構(gòu)造中，Borges等指出適當(dāng)增大WENO格式在間斷附近含間斷模板函數(shù)對(duì)應(yīng)的非線性權(quán)重，可以減小耗散，增加該格式在極值點(diǎn)處的精度，提高其對(duì)間斷的分辨率[3]。

在Jiang和Shu的經(jīng)典WENO格式基礎(chǔ)上，本文在構(gòu)造非線性權(quán)時(shí)引入間接光滑指數(shù)，增大間斷模板函數(shù)對(duì)應(yīng)的非線性權(quán)重，提出一種新的WENO格式，記為WENO-E，并進(jìn)行了相關(guān)數(shù)值驗(yàn)證。

1　WENO格式分析

ENO格式通過選取多個(gè)備選模板，分別構(gòu)造擬合函數(shù)，然后采用一個(gè)最光滑的插值模板對(duì)網(wǎng)格界面處數(shù)值通量進(jìn)行逼近，具體的ENO格式處理方法參見文獻(xiàn)[4-5]。

而WENO格式的基本思想是對(duì)每一個(gè)備選模板函數(shù)的網(wǎng)格界面數(shù)值通量逼近值取一定的權(quán)重，作加權(quán)平均，即得所求通量。

1.1WENO格式基本原理

考慮一維雙曲守恒律方程：

ut+f(u)x=0,x∈[a,b]

(1)

進(jìn)行如下均勻網(wǎng)格劃分：

a=x1/2

給定單元{xi}，單元尺度、單元、單元中心分別為Δx=xi+1/2-xx-1/2，Ii=[xi-1/2,xx+1/2]，xi=(xi-1/2+xx+1/2)/2。

方程(1)半離散型守恒差分形式為：

(2)

1.2經(jīng)典WENO(WENO-JS)格式加權(quán)方法

經(jīng)典WENO格式引入表征插值模板對(duì)應(yīng)的插值函數(shù)光滑程度的參數(shù)——光滑因子[2]：

(3)

則非線性權(quán)ωm構(gòu)造如下：

(4)

其中：dm為光滑情況下由高階與低階插值函數(shù)確定的線性權(quán)，亦稱之為理想權(quán)重；ε為一極小量；防止分母為0，一般取為10-8～10-6。

2　WENO格式加權(quán)方法改進(jìn)

經(jīng)典WENO加權(quán)方法中取權(quán)重時(shí)考慮了βm的二次方對(duì)相應(yīng)模板函數(shù)權(quán)重的影響，在光滑處使得非線性權(quán)ωm偏離理想權(quán)重dm較大，計(jì)算截?cái)嗾`差偏大[6]；在間斷附近使得間斷模板的權(quán)重取值接近于0，但這樣產(chǎn)生了過大的耗散，降低了對(duì)間斷和波動(dòng)的捕捉性能[3,7]。 因而本文采用一種新的思路對(duì)其進(jìn)行了改進(jìn)。 下文討論中將改進(jìn)格式記為WENO-E。

2.1WENO-E格式的加權(quán)處理

(5)

(6)

這樣使得新的非線性權(quán)中更多地考慮了各個(gè)模板函數(shù)光滑度的相互影響。其他計(jì)算與經(jīng)典WENO格式相同。

值得指出的是，式(6)中ε仍為一極小量，與式(4)中意義相同。在下文相關(guān)數(shù)值計(jì)算中，如無特別說明，均取ε=10-6。相關(guān)計(jì)算表明其取值對(duì)格式計(jì)算結(jié)果影響很小。

2.2權(quán)重分析

通過上文加權(quán)方法的闡述，現(xiàn)對(duì)5階的改進(jìn)前后兩種WENO格式中非線性權(quán)進(jìn)行分析?？紤]一般情況，以下分析取ε=0，省去其影響。

在光滑區(qū)，具體考察第一個(gè)插值模板的非線性權(quán)，有：

(7)

(8)

(9)

接下來討論間斷處二者的權(quán)重比較。假設(shè)5階WENO格式中第一個(gè)模板包含間斷，第二、三個(gè)模板光滑，即β0遠(yuǎn)大于β1，β2，則：

(10)

考慮如下間斷問題，數(shù)值上進(jìn)行驗(yàn)證：

ut+ux=0

(11)

對(duì)該問題，邊界條件取周期性邊界條件，網(wǎng)格數(shù)取200，第一次計(jì)算所得權(quán)重分布(局部區(qū)域)如圖1所示。

圖1表示非線性權(quán)在間斷附近的分布情況，縱坐標(biāo)取以10為底數(shù)的對(duì)數(shù)坐標(biāo)，其中虛線代表各權(quán)重對(duì)應(yīng)的理想權(quán)重，可以看出，在遠(yuǎn)離間斷處兩種格式的非線性權(quán)均與理想權(quán)重符合得很好。而在間斷附近兩者分布差異較大，較之WENO-JS格式，WENO-E格式對(duì)含有間斷的模板函數(shù)的權(quán)重分配都有所增大，尤其是在間斷左右兩個(gè)單元處，相應(yīng)權(quán)重從10-8～10-6增加至10-5～10-3，而這可以減小耗散，增加格式在極值點(diǎn)等處的精度。

(a) WENO-JS權(quán)重(a) Weights of the WENO-JS

(b) WENO-E權(quán)重(b) Weights of the WENO-E圖1　含間斷問題間斷附近權(quán)重分布Fig.1　Distribution of the weights ωkclose to the discontinuity

圖2給出網(wǎng)格數(shù)為200，t=8時(shí)，式(11)表示的問題在間斷附近的計(jì)算結(jié)果，圖中虛線代表解析解，實(shí)線代表WENO-JS的解。可以看出，在對(duì)間斷的捕捉上，WENO-E的計(jì)算結(jié)果更為陡峭，結(jié)果優(yōu)于WENO-JS。

圖2　線性遷移方程含有階躍的初始間斷條件數(shù)值計(jì)算結(jié)果Fig.2　Numerical solution of the advection equation with the discontinuous initial condition containing a step

2.3收斂精度分析

(12)

(13)

由以上分析可得WENO-E格式也具有5階精度。

3　數(shù)值實(shí)驗(yàn)

3.1線性遷移方程

對(duì)線性遷移方程ut+ux=0，分別取光滑與間斷初始條件進(jìn)行數(shù)值模擬，時(shí)間離散采用3階Runge-Kutta格式離散，邊界條件取周期性邊界條件。

3.1.1光滑問題

表1　初始條件為u(x,0)=sin(πx)(-1≤x≤1)時(shí)線性遷移方程計(jì)算誤差

3.1.2間斷問題

根據(jù)文獻(xiàn)[2]，對(duì)ut+ux=0取初始條件為如下間斷問題(從左至右依次為高斯分布曲線、矩形波、三角形波、半橢圓)：

(14)

取網(wǎng)格數(shù)為200，在t=8時(shí)，WENO-JS格式與WENO-E格式計(jì)算結(jié)果如圖3(a)所示，圖3(a)中虛線代表解析解，實(shí)線代表WENO-JS的解。

(a) 計(jì)算結(jié)果(a) Numerical solution

(b) 誤差(b) Error圖3　線性遷移方程取含有多種波形的初始間斷條件計(jì)算結(jié)果與誤差Fig.3　Numerical solution and absolute pointwise error of the advection equation with the discontinuous initial condition consisting of various waves

從圖出可以看出兩種格式定性上都可以較好地捕捉各間斷的形狀。從局部放大圖圖4中可得在間斷附近WENO-E的計(jì)算結(jié)果更為陡峭，尤其是在對(duì)高斯曲線的捕捉上，但是在矩形波的捕捉上WENO-E格式會(huì)引入輕微的震蕩。從誤差分布圖3(b)上可以看出這一震蕩影響很小，基本可以忽略不計(jì)；而在間斷附近的大多數(shù)網(wǎng)格節(jié)點(diǎn)處WENO-E的誤差小于WENO-JS。

圖4　圖3(a)局部放大圖Fig.4　Locally enlarged plot of Fig.3(a)

3.2一維非線性Burgers方程

一維非線性Burgers方程為：

(15)

對(duì)方程(15)取初始條件u(x,0)=0.5+sin(πx)，取周期為2的邊界條件，時(shí)間離散采用3階Runge-Kutta格式離散，進(jìn)行數(shù)值計(jì)算。

3.2.1光滑問題

時(shí)間為t=0.5/π時(shí)，計(jì)算結(jié)果仍是光滑的。表2給出了不同網(wǎng)格下兩種格式的誤差分析。由表2可以看出兩種格式整體均可以達(dá)到5階精度。但是，在相同網(wǎng)格下，WENO-E格式的計(jì)算誤差低于WENO-JS。

表2　初始條件為u(x,0)=0.5+sin(πx)時(shí)非線性Burgers方程計(jì)算誤差

表2(續(xù))

3.2.2間斷問題

時(shí)間為t=1.5/π時(shí)，解會(huì)產(chǎn)生激波。圖5給出網(wǎng)格數(shù)為40時(shí)WENO-JS格式與WENO-E格式的計(jì)算結(jié)果。

圖5中虛線代表解析解，實(shí)線代表WENO-JS的解。從圖5中可以看出兩種格式都可以較好地捕捉激波，且都沒有震蕩。由局部放大圖可以看出在激波附近WENO-E格式計(jì)算結(jié)果的點(diǎn)值在WENO-JS結(jié)果與解析解之間；計(jì)算激波左右點(diǎn)值(橫坐標(biāo)分別為1.2，1.25)連線的斜率，WENO-E和WENO-JS對(duì)應(yīng)的斜率分別為-24.3和-23.4，可見WENO-E格式計(jì)算結(jié)果更為陡峭。

圖5　初始條件為u(x,0)=0.5+sin(πx)時(shí)非線性Burgers方程計(jì)算結(jié)果Fig.5　Numerical solution of the nonlinear Burgers equation with the initial condition u(x,0)=0.5+sin(πx)

3.3一維歐拉方程

本文對(duì)一維歐拉方程含間斷的典型相關(guān)問題進(jìn)行了計(jì)算，考察WENO-E格式的計(jì)算性能。求解過程中使用Lax-Friedrichs通量分裂[5,9]，時(shí)間離散仍采用3階Runge-Kutta格式離散[5，10]。

歐拉方程控制方程組為：

(16)

3.3.1Sod問題

初始條件為：

(17)

這一問題，即經(jīng)典的黎曼激波管問題，表征被薄膜分開的兩種狀態(tài)不同的理想氣體在薄膜破裂后的運(yùn)動(dòng)狀態(tài)。取網(wǎng)格數(shù)為50，分別采用兩種WENO格式計(jì)算該問題，t=0.4時(shí)刻所得速度分布如圖6所示，圖中虛線表示解析解，實(shí)線代表WENO-JS格式的解。由圖6可得，兩種格式都可以較好地捕捉間斷，但從速度局部放大圖中可以看到WENO-E格式對(duì)激波的捕捉更為陡峭，更為準(zhǔn)確。

圖6　Sod問題速度分布Fig.6　Distribution of the velocity of the Sod problem

3.3.2LeBlanc激波管

初始條件為：

(18)

這是初始條件為式(18)的黎曼激波管問題。取網(wǎng)格數(shù)為2400，分別采用兩種WENO格式計(jì)算該問題，t=0.001時(shí)刻所得密度與壓力分布如圖7所示，圖中虛線表示解析解，實(shí)線代表WENO-JS格式的解。由圖7可得，兩種格式對(duì)激波位置的捕捉誤差較大，但從局部放大圖中可以看到WENO-E格式的計(jì)算結(jié)果點(diǎn)值在WENO-JS結(jié)果左側(cè)，離解析解激波位置稍近一些，對(duì)激波的捕捉有一定改善。此外，這一問題具有很強(qiáng)的間斷性，WENO-E格式可以取得較好的計(jì)算結(jié)果，足以證明WENO-E格式的穩(wěn)定性。

(a) 密度(a) Density

(b) 壓力(b) Pressure圖7　LeBlanc激波管問題數(shù)值結(jié)果Fig.7　Numerical solution of the LeBlanc shock tube problem

3.3.3Shu-Osher問題

初始條件為：

(19)

該問題表征一維情況下馬赫數(shù)為3的激波與正弦熵波的相互作用[5，11]。取網(wǎng)格數(shù)為200，t=1.8時(shí)刻所得計(jì)算結(jié)果密度分布如圖8所示。圖8中虛線表示基準(zhǔn)解(基準(zhǔn)解取網(wǎng)格數(shù)為2000時(shí)采用WENO-JS格式求得的數(shù)值結(jié)果)。由圖8可得，WENO-E格式相對(duì)于WENO-JS格式可分辨出更多波動(dòng)細(xì)節(jié)，對(duì)波峰波谷的捕捉更為精確，尤其是在激波后的高頻波動(dòng)分布區(qū)。而WENO-JS格式的計(jì)算結(jié)果則對(duì)波動(dòng)細(xì)節(jié)抹平更為嚴(yán)重。

圖8　Shu-Osher問題密度分布Fig.8　Distribution of the density of the Shu-Osher problem

此外，還對(duì)兩種格式的計(jì)算時(shí)間進(jìn)行了比較。數(shù)值結(jié)果表明WENO-E格式計(jì)算所用時(shí)間與WENO-JS格式接近，計(jì)算效率基本相同。

4　結(jié)論

在Jiang和Shu的經(jīng)典WENO格式的基礎(chǔ)上，通過引入間接光滑指數(shù)，改變非線性權(quán)構(gòu)造過程，構(gòu)造出WENO-E格式。

理論分析和數(shù)值結(jié)果表明，該格式保持了與WENO-JS格式基本相同的計(jì)算效率，取得了與之同階的計(jì)算精度；但在相同網(wǎng)格下，較WENO-JS格式，該格式對(duì)光滑區(qū)域的計(jì)算有更小的截?cái)嗾`差，計(jì)算更為準(zhǔn)確，對(duì)間斷的捕捉更為陡峭，分辨率更高。研究結(jié)果為WENO格式的改進(jìn)提出一種新的思路，即引入間接光滑指數(shù)，考慮了各個(gè)模板函數(shù)光滑度的相互影響，對(duì)非線性權(quán)優(yōu)化構(gòu)造具有一定的指導(dǎo)意義。

References)

[1]Liu X D, Osher S, Chan T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200-212.

[2]Shan G J, Wang C S . Efficient implementation of weighted ENO schemes[R]. Efficient Implementation of Weighted ENO Schemes, 1995.

[3]Borges R, Carmona M, Costa B, et al. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[J]. Journal of Computational Physics, 2008, 227(6): 3191-3211.

[4]Shu C W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes[J]. Journal of Computational Physics, 1988, 77(2): 439-471.

[5]Shu C W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes Ⅱ[J]. Journal of Computational Physics, 1989, 83(1): 32-78.

[6]Shen Y, Zha G. Improvement of the WENO scheme smoothness estimator[J]. International Journal for Numerical Methods in Fluids, 2010, 64(6): 653-675.

[7]Arshed G M, Hoffmann K A. Minimizing errors from linear and nonlinear weights of WENO scheme for broadband applications with shock waves[J]. Journal of Computational Physics, 2013, 246(8): 58-77.

[8]Roe P L. Approximate riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[9]Henrick A K, Aslam T D, Powers J M. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points [J]. Journal of Computational Physics, 2005, 207(2): 542-567.

[10]Castro M, Costa B, Don W S. High order weighted essentially non-oscillatory WENO[J]. Journal of Computational Physics, 2011, 230(5): 1766-1792.

[11]侯中喜, 易仕和, 李樺. 高精度高分辨率WEN0格式分析與改進(jìn)[J]. 國防科技大學(xué)學(xué)報(bào), 2003, 25(1): 17-20.

HOU Zhongxi, YI Shihe, LI Hua. Analysis and improvement of high precision, high resolution WENO schemes[J]. Journal of National University of Defense Technology, 2003, 25(1): 17-20. (in Chinese)

Analysis and improvement of weights in WENO schemes

CHAI Delin， LI Yurun， SUN Zhongguo， XI Guang

(School of Energy and Power Engineering， Xi′an Jiaotong University， Xi′an 710049， China)

In order to improve the computational performance of the WENO scheme, a new WENO scheme, namely WENO-E scheme was constructed, which reduces dissipation close to discontinuities. Based on the analysis of the algorithm for weighted factors in the classical WENO scheme (namely WENO-JS) proposed by Jiang and Shu, the new scheme was constructed by introducing indirect smooth indicator. Theoretical analysis shows that the WENO-E scheme can reach the same convergence order of WENO-JS with the same computational efficiency; while it can obtain smaller truncation errors at smooth parts of the solution and higher resolution close to the discontinuities with the same grids than the WENO-JS. Subsequently, compared with the classical WENO scheme, when numerical experiments with the linear transport equation, the nonlinear Burgers equation and the one dimensional Euler system of equations are conducted, the WENO-E scheme achieves better numerical solutions.

weighted essentially non-oscillatory schemes； smoothness indicator； high precision; high resolution

10.11887/j.cn.201604007http://journal.nudt.edu.cn

2016-01-15

國家自然科學(xué)基金資助項(xiàng)目(51236006)

柴得林(1991—)，男，山西臨汾人，博士研究生，E-mail：cunzhang.chai@stu.xjtu.edu.cn；席光(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：xiguang@mail.xjtu.edu.cn

V211.1; O241

A

1001-2486(2016)04-039-07