求解雙曲守恒律方程的三階修正模板WENO格式*

王亞輝

（鄭州師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，鄭州 450044）

引 言

非線性雙曲守恒律的解通常包含間斷和光滑的小尺度結(jié)構(gòu).在不產(chǎn)生虛假振蕩的情況下捕捉間斷，并盡可能多地求解小尺度結(jié)構(gòu)，是數(shù)值格式的必要條件.在求解雙曲守恒律方程的眾多數(shù)值格式中，加權(quán)本質(zhì)無振蕩(WENO)格式已成為最常用的方法之一.WENO 是求解不連續(xù)雙曲守恒律的一種常用方法，其抑制虛假振蕩的成功歸功于一組動態(tài)候選模板的低階多項式的非線性組合方法.加權(quán)過程是通過根據(jù)每個候選模板的光滑度評估其局部通量的貢獻來執(zhí)行的，具體地說，WENO 格式中的光滑度指示子可以驅(qū)使格式向最優(yōu)方向發(fā)展，并通過分配間斷模板的基本零權(quán)重來避免跨間斷的插值.

由Liu 等首先介紹的WENO 格式[1]使用了所有候選模板的非線性權(quán)來凸組合ENO 格式[2-5]，以提高解在光滑區(qū)域的精度，而不丟失ENO 格式在間斷附近的無振蕩特性.該過程通過根據(jù)候選模板上解的光滑度指示子對局部通量的貢獻進行加權(quán)來執(zhí)行，使得包含間斷局部通量的模板權(quán)重基本上為零.Jiang 和Shu[6]設(shè)計了一個光滑度指標(biāo)子，即所有低階多項式的各階導(dǎo)數(shù)的L2范數(shù)的歸一化平方和，從而發(fā)展了經(jīng)典的有限差分WENO(WENO-JS)格式.Henrick 等[7]提出了WENO-M 格式，用映射函數(shù)修改了權(quán)，使之滿足格式最優(yōu)收斂階的充分必要條件.Borges 等[8]發(fā)展了WENO-Z 格式，該方法將全局高階光滑度指示子納入權(quán)重的構(gòu)造中.WENO-Z 格式可以達(dá)到與WENO-M 格式相同的精度.由于該格式賦予間斷模板較大的權(quán)重，因此發(fā)展了七階到十一階保單調(diào)WENO 格式.Gerolymos 等[9]將WENO-Z 格式推廣到十一階WENO 格式，發(fā)展了七階到十一階WENO-Z 格式.

與高階WENO 格式相比，三階WENO 格式具有許多優(yōu)點.例如，它對激波問題具有很強的魯棒性，使用較少的網(wǎng)格點，從而降低了邊界處理和求解的難度，很容易地推廣到非結(jié)構(gòu)網(wǎng)格上.盡管三階WENO 格式具有上述優(yōu)點，但三階WENO 格式很少受到關(guān)注.這是因為經(jīng)典的三階WENO 格式在光滑區(qū)域的精度在2 到3 階之間，而在臨界點退化為1 階，而三階WENO 格式由于成本較低，是理想的格式.利用WENO-Z 格式的加權(quán)過程，可以解決提高分辨率和恢復(fù)最優(yōu)收斂階的問題.通過這種方法，Wu 和Zhao[10]開發(fā)了一個改進的三階WENO(WENO-N3)格式.他們設(shè)計了一個四階參考光滑度指示子 τ4N，由候選模板的光滑度指示子與全局模板的光滑度指示子的線性組合得到的.然而 τ4N不能使格式在臨界點附近達(dá)到最優(yōu)三階收斂.為了彌補這個缺陷，Wu 等[11]用τ4Np的p冪形式改進了參考光滑度指示子，得到了WENO-Np3 格式.Xu 和Wu[12]提出了另一個四階參考光滑度指示子 τ4P，這樣得到的WENO-P3 格式耗散性更小.然而，該算法也只能在沒有臨界點的光滑區(qū)域?qū)崿F(xiàn)三階精度，在有臨界點的光滑區(qū)域降到一階準(zhǔn)確度.Wang 等[13]給出了三階WENO 格式的參考光滑度指示子 τ4Rp，理論分析表明，所得到的WENO-Rp3 格式(p=1.5)在光滑區(qū)域(包括臨界點)能達(dá)到三階精度，并且WENO-R3 格式(p=1)保持了ENO 性質(zhì).關(guān)于三階WENO 格式的其他改進格式可參見文獻[14-15].

為了降低經(jīng)典三階WENO 格式的數(shù)值耗散，提高數(shù)值精度，這里我們提出了一種新的WENO 格式的改進模板近似方法.這里我們所采取的主要策略與文獻[10-13]的不同，不是構(gòu)造更高階的參考光滑度指示子，而是改進經(jīng)典WENO-JS3 格式中各候選模板上數(shù)值通量的一階多項式逼近，通過加入二次項使模板逼近達(dá)到三階精度，并且計算了相應(yīng)的候選通量.它能使所得格式(WENO-MS)在包括一階臨界點在內(nèi)的光滑區(qū)域內(nèi)達(dá)到三階收斂點數(shù).為了抑制強激波和間斷區(qū)域的明顯數(shù)值振蕩，在不影響計算精度的前提下，我們在候選模板通量的修正項中加入可調(diào)函數(shù) φ，并給出了一系列一維和二維數(shù)值算例來驗證新格式的性能.數(shù)值結(jié)果表明，與WENO-JS3、WENO-Z3 格式相比，WENO-MS 格式對精細(xì)光滑結(jié)構(gòu)具有相當(dāng)或更高的分辨率.

本文的主要結(jié)構(gòu)如下：第1 節(jié)簡要回顧了一維標(biāo)量守恒律的三階WENO 格式；第2 節(jié)提出了新的修正模板WENO 近似方法，并給出了WENO-MS 格式；第3 節(jié)給出了數(shù)值算例來說明新格式的性能；第4 節(jié)對全文進行了總結(jié).

1 三階WENO 格式綜述

在這一節(jié)中，我們簡要介紹了一維標(biāo)量守恒律的三階WENO 有限差分格式：

設(shè)Ij是 給定域的一個分區(qū)，第j個Ij:=[xj?1/2，xj+1/2].Ij的 中心表示為函數(shù)f在位置xj的 值由下標(biāo)j表示，即fj=f(xj).我們假設(shè)網(wǎng)格點{xj+1/2}均 勻分布，單元大小Δx=xj+1/2?xj?1/2.

為了構(gòu)造方程(1)的守恒有限差分格式，可以通過隱函數(shù)h(x)來定義數(shù)值通量[5]：

對方程(2)求關(guān)于x的 導(dǎo)數(shù)得到

因此，式(1)可由常微分方程組近似，其中空間導(dǎo)數(shù)由守恒有限差分近似，產(chǎn)生半離散形式：

圖1 三階WENO 數(shù)值通量的模板Fig.1 Stencils for the 3rd-order WENO numerical flux

1.1 三階WENO 格式

它們結(jié)合起來定義了格式的數(shù)值通量：

這里的非線性權(quán)重 ωk定義為

其中，常系數(shù)dk被稱為理想權(quán)重，因為其與的 線性組合保持最優(yōu)三階收斂到h(xj+1/2)，即

dk的具體值如下所示(參見文獻[16])：

式(9)中的0

其具體形式為

Borges 等[8]提出了一個WENO-Z 權(quán)：

其中 ? =10?40，τ3=|β1?β0|是 一個三階參考光滑探測子，它驅(qū)使非線性權(quán)重 ωk(式(13)的WENO-Z)向最優(yōu)權(quán)重dk靠近，并且比J-S 權(quán)重(9)快.當(dāng)冪q=2時，WENO-Z 格式在一階臨界點處可以達(dá)到最優(yōu)的三階，但在間斷附近具有更多的耗散.一個不那么耗散的方法是設(shè)計一個足夠的高階參考光滑探測子，并且總是使用冪q=1.

為了要為WENO-Z 權(quán)設(shè)計一個合適的參考光滑探測子(13)，需要一個充分條件來恢復(fù)三階WENO 格式的最優(yōu)階精度.文獻[12]給出了三階WENO 格式達(dá)到三階精度的一個條件：

方程(14)僅僅是一個充分條件，但不是必要條件.然而，這一條件為設(shè)計高階參考光滑性指示子提供了有用的信息.

2 新的修正三階WENO 格式

2.1 修正模板近似

“MS”表示“修正模板”.

注1修改后的候選通量(25)依賴于整個三點模板，實際上，通過簡單的計算不難發(fā)現(xiàn)，它就是三階迎風(fēng)格式，因此新的格式(26)失去了ENO 性質(zhì).為了抑制強激波和不連續(xù)區(qū)域內(nèi)明顯的數(shù)值振蕩，在這里，應(yīng)用了一個可調(diào)函數(shù) φ來限制式(25)中最后一項的影響，使新發(fā)展的格式具有ENO 性質(zhì).修正后的候選通量(25)變?yōu)?/p>

這里φ 應(yīng)該滿足以下條件：

(a) 如果模板S3={xj?1,xj,xj+1}是 一個間斷模板，則φ 是一個很小的值；

(b) 如果模板S3={xj?1,xj,xj+1}是 一個光滑模板，則φ 是 接近1，并且φ 不會影響新格式的收斂精度.

所以φ 的具體形式定義如下：

因而在接下來的數(shù)值算例中我們?nèi)《╭=2.

3 數(shù)值測試

在本節(jié)中，我們將在幾個數(shù)值例子中驗證所提出的WENO-MS 格式(WENO-MS-JS 和WENO-MS-Z)的性能，并與WENO-JS 和WENO-Z 格式進行比較.數(shù)值算例從簡單標(biāo)量對流方程求解開始，然后是一維和二維Euler 方程的數(shù)值求解.本文所有的格式和算例均采用局部Lax-Friedrichs 通量分裂.時間推進使用三階TVD-Runge-Kutta 方法[5]：

其 中 σ是CFL(Courant-Friedrichs-Lewy)條件數(shù).

3.1 標(biāo)量測試問題

考慮具有周期邊界的標(biāo)量對流方程條件并用不同的初始數(shù)據(jù)測試任意初始剖面的傳播，包含跳躍間斷和角點.

例1 (線性方程)讓我們考慮如下線性對流方程：

初值條件為

u(x,0)=u0(x).

我們首先在兩組初始數(shù)據(jù)的線性方程(31)上檢驗了WENO-MS 格式的數(shù)值收斂階：

時間步長為Δt=0.6Δx， 保證與空間精度兼容.誤差的L1范 數(shù)是通過與t=2的精確解比較得出的，根據(jù)以下式子計算：

表1 和表2 分別顯示t=2時 初值(32a)的L1、L∞誤差和收斂階.我們發(fā)現(xiàn)誤差隨著WENO-JS、WENOZ3、WENO-MS-JS3 和WENO-MS-Z3 的順序而減小.新的WENO-MS 格式能達(dá)到最優(yōu)三階收斂精度.表3 和表4顯示初值(32b)的結(jié)果.盡管有臨界點的出現(xiàn)，但是新的WENO-MS 格式的數(shù)值收斂階仍然能達(dá)到最優(yōu)階.相比之下，新的WENO-MS 格式的性能優(yōu)于其他WENO 格式.

表1 線性對流方程(31)在初值(32a)下，不同格式在t =2.0的L1 誤差和收斂階Table 1 L1 errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32a)

表2 線性對流方程(31)在初值(32a)下，不同格式在t =2.0 的 L∞誤差和收斂階Table 2 L ∞ errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32a)

表3 線性對流方程(31)在初值(32b)下，不同格式在t =2.0的L1 誤差和收斂階Table 3 L1 errors and convergence rates at t =2.0 of different schemes for linear advection eq.(31) with initial data eq.(32b)

表4 線性對流方程(31)在初值(32b)下，不同格式在t =2.0的L∞誤差和收斂階Table 4 L∞ errors and convergence rates at t= 2.0 of different schemes for linear advection eq.(31) with initial data eq.(32b)

然后我們討論WENO-MS 格式在長時間計算中的性能.具體來說，在式(31)上用

這是一個分段正弦函數(shù)，在x=0處 有一個跳躍，其中CFL 數(shù)為0.5，我們將其解到t=40處，以觀察格式在跳躍間斷處的行為.圖2 顯示了使用網(wǎng)格點N=400的不同WENO 格式的數(shù)值解與解析解的比較.我們可以看到，WENO-MS-JS3 和WENO-MS-Z3 格式的耗散比WENO-JS3 和WENO-Z3 格式小.

圖2 線性對流方程(31)在初值(33)下數(shù)值解與解析解的比較，t =40 Fig.2 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (33), at t =40

我們進一步考慮初始條件：

我們求解方程(31)，初始條件(34)在計算域 ?1≤x≤1中 時間為t=41， 其中 Δx=0.005，CFL 數(shù)為0.5.數(shù)值結(jié)果如圖3 所示.易見在間斷附近，WENO-MS-JS3 和WENO-MS-Z3 的耗散比WENO-JS3 和WENO-Z3 小.

圖3 線性對流方程(31)在初值(34)下不同格式的數(shù)值解與解析解的比較，t =41 ，N=400Fig.3 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (34), at t=41, N=400

最后，我們求解線性對流方程(31)，初始條件包括Gauss 波、方波、三角波和半橢圓波：

其中G(x,β,z)=e?β(x?z)2，F(xiàn)(x,α,a)=a=0.5,z=?0.7，δ=0.005，α=10和 β=(ln2)/(36δ2).這種情況的解包括接觸間斷，角點奇異和光滑區(qū)域.我們用時間步長 Δt=CCFLΔx，CFL 數(shù)CCFL為0.5，計算到最終時間為t=6.數(shù)值結(jié)果如圖4 所示.結(jié)果表明，WENO-MS 格式的性能優(yōu)于WENO-JS3 和WENO-Z3 格式，特別是在間斷面附近.所有改進的三階WENO 格式都比經(jīng)典的WENO-JS3 格式具有更好的分辨率，其中WENO-MS 格式的耗散最小.

圖4 線性對流方程(31)在初值(35)下不同格式的數(shù)值解與解析解的比較，t =6 ，N=400Fig.4 Comparison of the analytical solution with the numerical solutions of linear advection eq.(31) with initial value (35), at t=6, N=400

3.2 一維Euler 系統(tǒng)

在這一小節(jié)中，我們將給出雙曲守恒律方程的數(shù)值結(jié)果.考慮一維Euler 方程組：

其中

U=(ρ,ρu,E)T,F(U)=(ρu,ρu2+p,u(E+p))T.

狀態(tài)方程為

ρ,u,p和E分別是密度、速度、壓強和總能， γ 是比熱比.相應(yīng)的Jacobi 矩陣A(U)=?F/?U的特征值為

λ1=u?c, λ2=u, λ3=u+c,

這里c=(γp/ρ)1/2為聲速.采用特征場的特征分解與局部Lax-Friedrichs 分裂，并將WENO 格式推廣到一維Euler 系統(tǒng).

例2 (一維 Riemann 問題)我們考慮具有如下形式的三個經(jīng)典初值問題：

① Sod 激波管[17-18]

其中γ =1.4,CCFL=0.6， 網(wǎng)格尺度為Δx=0.005， 計算時間為t=0.2.密度場的數(shù)值結(jié)果與精確的Riemann 解進行了比較，如圖5 所示.可以看出，WENO-MS-Z3 格式的性能優(yōu)于WENO-MS-JS3、WENO-Z3 和WENO-JS3 格式.與其他WENO 格式相比，WENO-MS-Z3 格式在很好地捕捉激波的同時，具有最小的耗散性.

圖5 Sod 激波管問題[16]的數(shù)值結(jié)果，t =0.2 ，N =200，CCFL=0.6Fig.5 Numerical results of the Sod problem[16], at t=0.2 ， N=200，CCFL=0.6

② 沖擊波相互作用[17-18]

其中γ =1.4,CCFL=0.6.我們用網(wǎng)格尺度為 Δx=0.001 25計 算到t=0.038.密度場的數(shù)值結(jié)果如圖6 所示.由于精確解是未知的，參考解是通過使用五階WENO-JS 格式[6]在3 201 個網(wǎng)格點上獲得的.易知, 在接觸間斷處的耗散的順序是WENO-MS-Z3

圖6 沖擊波的相互作用的數(shù)值結(jié)果，t =0.038 ，N =800，CCFL=0.6Fig.6 Numerical results of interacting blast waves, at t =0.038 ， N =800，CCFL = 0.6

③ 一維激波等熵波作用(Shu-Osher 問題)[5]我們計算[ ?5,5]區(qū)域上的近似解，這里采用周期邊界條件.初始條件為

其中 γ =1.4,CCFL=0.6， ? 和k分 別是熵波的振幅和波數(shù).這里，我們設(shè)置? =0.2，k=5.在這個問題中，一個右移的超音速(Ma=3)激波與密度擾動中的正弦波相互作用，產(chǎn)生一個具有光滑結(jié)構(gòu)和不連續(xù)性的流場.這股氣流在波數(shù)高于初始密度變化波數(shù)k的右行激波后面誘導(dǎo)波跡.由于精確解是未知的，參考解是通過使用五階WENO-JS 格式[6]在3201 個網(wǎng)格點上獲得的.圖7 將密度剖面的數(shù)值結(jié)果與參考解進行了比較.所有的WENO 格式都比WENO-JS3 格式的耗散性小.WENO-MS 格式比WENO-Z3 和WENO-JS3 格式更好地捕捉了右行激波后面的高頻波，而且WENO-MS 格式的分辨率最接近WENO-JS5 格式.

圖7 激波等熵波相互作用(Shu-Osher)[5]的密度分布，t =1.8 ，N =401，CCFL=0.6Fig.7 Density profiles of the shock entropy interacting of Shu-Osher[5] at t =1.8 with 401 points and CCFL=0.6

3.3 二維 Euler 系統(tǒng)

在本小節(jié)中，考慮二維可壓縮Euler 系統(tǒng)：

其中

這里 ρ ,u,v,p和E分別是密度，x和y坐 標(biāo)方向上的速度分量，壓力和總能量.另外，U是守恒變量的向量，F(xiàn)(U)是x方向的通量分量，G(U)是y方向的通量分量.對二維可壓縮Euler 方程進行了逐維求解.

例3 (Rayleigh-Taylor 不穩(wěn)定)當(dāng)加速度從較重流體定向到較輕流體時，Rayleigh-Taylor 不穩(wěn)定性發(fā)生在兩種密度流體之間的界面上.這一問題在文獻中得到了廣泛的模擬[19-20].計算域為[ 0,1/4]×[0,1]，初始條件為

比熱比γ =5/3.引力效應(yīng)是通過將源項S=(0,0,ρ,ρv)T添加到二維Euler 方程(37)的右側(cè)來引入的.左右邊界為反射邊界條件，頂部和底部邊界設(shè)置為 (ρ,u,v,p)=(1,0,0,2.5)和 (ρ,u,v,p)=(2,0,0,1)，CFL 數(shù)為0.6，時間為t=1.95.圖8 顯示W(wǎng)ENO-MS-Z3 格式比WENO-JS3、WENO-Z3 和WENO-MS-JS3 格式結(jié)構(gòu)更為豐富，尤其是在界面上產(chǎn)生的渦流較多，表明其耗散性比其他格式小.

圖8 不同格式關(guān)于Rayleigh-Taylor 不穩(wěn)定性問題的密度等值線，Δ x=Δy=1/240，C CFL=0.6，t=1.95Fig.8 Density contours of the Rayleigh-Taylor instability at t =1.95 with Δ x=Δy=1/240，CCFL=0.6

例4 (雙Mach 反射問題)[19]激波在斜面上的二維雙Mach 反射描述了平面Mach 激波在空氣中撞擊楔形物的反射.這種試驗被廣泛用于驗證數(shù)值方法的性能.我們在[ 0,4]×[0,1]上計算這個測試問題，并像往常一樣在[0,3]×[0,1]上 顯示結(jié)果.一開始施加一個右移Mach 數(shù)為10 的激波在x=1/6處 ，激波前緣與x軸成 60?角.在沿底部邊界x=0至x=1/6的區(qū)域內(nèi)，為初始激波后流動賦值，其余區(qū)域采用反射邊界條件.左邊界為初始激波后流動賦值.對于x=4處的右側(cè)邊界都設(shè)置為零.上邊界是用來描述Mach 數(shù)為10 的激波的精確運動.有關(guān)此問題的詳細(xì)描述，請參見文獻[19].我們采用的網(wǎng)格尺度為 Δx=Δy=1/480， 比熱比 γ =1.4，CCFL=0.6，計算時間t=0.2.圖9 比較了WENO-MS 格式與WENO-JS3 和WENO-Z3 格式的數(shù)值結(jié)果.可見，WENO-MS-JS3 格式的分辨率與WENO-Z3 格式大致相當(dāng)，WENO-MS-Z3 格式分辨率最高，WENO-MS-Z3 較好地解決了Mach桿附近的不穩(wěn)定性問題.

圖9 雙Mach 反射問題在t =0.2時 的密度等值線，(網(wǎng)格點為1 920×480)Fig.9 Density contours of the double Mach reflection problem at t =0.2 with 1 920×480 grid points

4 結(jié) 論

為了減少傳統(tǒng)的三階加權(quán)本質(zhì)無振蕩格式的數(shù)值耗散以及提高經(jīng)典三階WENO 格式的數(shù)值精度，本文發(fā)展了非線性雙曲守恒律數(shù)值解WENO 三階WENO-MS 格式的一種新的修正模板逼近方法.通過改進經(jīng)典WENO-JS3 格式中各候選模板上數(shù)值通量的一階多項式逼近，并加入二次項使模板逼近達(dá)到三階精度.計算出了相應(yīng)的候選通量.它可以新發(fā)展的WENO-MS 格式在光滑區(qū)域(包括一階臨界點)實現(xiàn)三階收斂.在不影響計算精度的前提下，為了使新的修正模板WENO-MS 格式具有ENO 特性，我們在候選模板通量的修正項中加入可調(diào)函數(shù) φ，并且給出了一系列一維和二維數(shù)值算例，驗證了新格式的有效性.數(shù)值結(jié)果表明，與WENOJS3、WENO-Z3 格式相比，WENO-MS 格式對精細(xì)光滑結(jié)構(gòu)具有相當(dāng)或更高的分辨率.