程曉晗，封建湖，聶玉峰

(1.西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西 西安 710072； 2.長安大學(xué)理學(xué)院，陜西 西安 710064)

許多物理模型都采用雙曲型守恒律方程描述，如氣體動力學(xué)中的Euler方程、淺水動力學(xué)中的淺水波方程和等離子物理學(xué)中的磁流體方程等。在一維情況下，有如下形式：

(1)

本文中，進一步研究熵相容格式，通過在單元交界面處進行高階重構(gòu)，得到一類高精度的加權(quán)基本無振蕩(weighted essentially non-oscillatory,WENO)型熵相容格式。最后，通過一些數(shù)值算例驗證所構(gòu)造的新格式的性能。

1 WENO型熵相容格式

1.1 一維標(biāo)量守恒律方程

1.1.1熵守恒格式

(2)

(3)

與F相容。該格式被稱為熵守恒格式，具有二階精度[2]。

1.1.2熵穩(wěn)定格式

根據(jù)熵穩(wěn)定條件，在激波等間斷位置，總熵應(yīng)該耗散(減少)。E.Tadmor[2]指出，一個三點格式只需含有比熵守恒格式更多的黏性則是熵穩(wěn)定的。因此，P.L.Roe[3]提出用Roe格式的數(shù)值黏性項作為對熵守恒格式的補充，數(shù)值通量為：

(4)

該格式保持了Roe格式對間斷的良好捕捉效果，并且無需進行熵修正。

1.1.3熵相容格式

考慮黏性形式的Roe格式數(shù)值通量[5]：

(5)

(6)

式中：Δai+1/2=f′ui+1-f′ui，取α=1/6。式(6)定義的數(shù)值通量跨過激波時具有足夠熵增且數(shù)值解保持單調(diào)，具有這種性質(zhì)的通量稱為熵相容通量。熵相容格式是目前估計熵的變化最精確的一種熵穩(wěn)定格式。

1.1.4WENO型熵相容格式

(7)

1.2 一維守恒律方程組

氣體動力學(xué)中的Euler方程為：

(8)

1.2.1熵守恒格式

對于守恒律系統(tǒng)，一般采用逐段分解的構(gòu)造方法[7]獲得熵守恒通量。在該方法中，保持總熵不變的數(shù)值通量可以表示為沿熵變量空間內(nèi)不同路徑的積分和的形式。雖然該方法適用于各種守恒律系統(tǒng)，但其計算量過大并且有時候會產(chǎn)生數(shù)值不穩(wěn)定[8]。對于Euler方程和淺水波方程，根據(jù)其特定的熵函數(shù)，他們的熵守恒通量有比較簡單的構(gòu)造方法。

其中，Euler方程的熵守恒通量的計算方法如下。定義參數(shù)變量z為：

(9)

則熵守恒通量為：

(10)

式中：

1.2.2熵穩(wěn)定格式

與標(biāo)量守恒律一樣，對熵守恒格式添加適當(dāng)?shù)酿ば?，可使之滿足離散熵不等式?？紤]用Roe格式的數(shù)值黏性項作為補充，得到關(guān)于Euler方程的熵穩(wěn)定格式，數(shù)值通量為：

(11)

1.2.3熵相容格式

類似標(biāo)量守恒律，為了抵消解在跨過激波時產(chǎn)生的熵增，需要對D進行修正，令：

D1=Λ+αΔΛS

(12)

式中：ΔΛ=diagui+1-ai+1-ui-ai,ui+1-ui,ui+1+ai+1-ui+ai，于是可以得到Euler方程的熵相容格式，數(shù)值通量為：

(13)

1.2.4WENO型熵相容格式

(14)

對于一維守恒律方程組，本文中僅以氣體動力學(xué)中的Euler方程為例，設(shè)計WENO型熵相容格式，其他方程可以參照Euler方程做一些相應(yīng)的修改，在此不再贅述。

2 數(shù)值模擬算例

采用半離散格式：

(15)

時間上的推進一律采用三階強穩(wěn)定的Runge-Kutta方法[9]。為了比較，也給出了相應(yīng)算例的熵相容格式(EC格式)的計算結(jié)果。

2.1 線性對流方程連續(xù)初值問題

在區(qū)域x∈[-1,1]上求解初值問題：

采用周期邊界條件，計算到t=10。

本算例可以用來測試格式的收斂性和精度，表1給出了相應(yīng)的誤差L1和L，充分說明了WENO型熵相容格式(EC-WENO格式)的精度能達到五階。

表1 EC-WENO格式的數(shù)值精度Table 1 Numerical accuracy of EC-WENO scheme

2.2 無黏Burgers方程間斷初值問題

在區(qū)域x∈[-1,1]上求解初值問題：

采用周期邊界條件，取CFL系數(shù)為0.8，空間網(wǎng)格數(shù)為40，計算到t=0.3。

本算例的解主要包括2部分：在左邊x=-1/3處由-1到1的跳躍產(chǎn)生了一個稀疏波，在右邊x=1/3處由1到-1的跳躍產(chǎn)生了一個定常激波。計算結(jié)果如圖1所示。由圖可知，EC-WENO格式算得的解不僅準(zhǔn)確地捕捉到稀疏波，而且在間斷位置也沒有產(chǎn)生偽振蕩，符合“高分辨率”的特性。

圖1 無黏Burgers方程間斷初值問題Fig.1 Discontinuous initial value problem of Burgers equation

2.3 一維Euler方程Sod激波管問題

在區(qū)域[-0.5,0.5]上求解初值問題：

采用Neumann邊界條件，取CFL系數(shù)為0.4,空間網(wǎng)格數(shù)為100，計算到t=0.1。

精確解包括稀疏波、接觸間斷和激波。密度的計算結(jié)果如圖2所示。由圖可知，EC-WENO格式比EC格式能對解的結(jié)構(gòu)有更準(zhǔn)確地捕捉，尤其是對激波的捕捉非常銳利。

圖2 一維Euler方程Sod激波管問題Fig.2 Sod shock tube problem of 1D Euler equation

2.4 一維Euler方程低密度流問題

在區(qū)域[-0.5,0.5]上求解初值問題：

采用Neumann邊界條件，取CFL系數(shù)為0.4,空間網(wǎng)格數(shù)為100，計算到t=0.05。

精確解包括2個強稀疏波和1個平凡的接觸間斷。由于中間部分的壓力非常小，接近于真空狀態(tài)，很多數(shù)值方法容易在該位置出現(xiàn)壓力為負的情況(如Roe格式)，導(dǎo)致計算失敗。壓力的計算結(jié)果如圖3所示。由圖可知，EC-WENO格式能成功地計算低密度流問題，且效果較好，能有效避免非物理現(xiàn)象(如負壓力)。

圖3 一維Euler方程低密度流問題Fig.3 Low density problem of 1D Euler equation

2.5 一維Euler方程強稀疏波問題

在區(qū)域[-10,15]上求解初值問題：

采用Neumann邊界條件，取CFL系數(shù)為0.4,空間網(wǎng)格數(shù)為100，計算到t=0.01。

精確解包括強稀疏波、接觸間斷和激波。密度的計算結(jié)果如圖4所示。由圖可知，EC-WENO格式取得了非常良好的效果，解在稀疏波有一個平滑的過渡，且在波頭波尾位置都非常貼近精確解。

圖4 一維Euler方程強稀疏波問題Fig.4 Density of strong expansion problem of 1D Euler equation

3 結(jié) 論

以熵相容格式為基礎(chǔ)，通過在單元交界面處進行WENO重構(gòu)，構(gòu)造了一類求解雙曲守恒律方程的WENO型熵相容格式，有效地提高了熵相容格式的精度。數(shù)值模擬結(jié)果表明，新格式具有高精度、基本無振蕩性等特點。

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

[2] Tadmor E. The numerical viscosity of entropy stable schemes for systems of conservation laws, Ⅰ[J]. Mathematics of Computation, 1987,49(179):91-103.

[3] Roe P L. Affordable, entropy-consistent, flux functions[C]∥Oral Talk at Eleventh International Conference on Hyperbolic Problems: Theory, Numerics, Applications. Lyon, France, 2006.

[4] Ismail F, Roe P L. Affordable, entropy-consistent Euler flux functions, Ⅱ: Entropy production at shocks[J]. Journal of Computational Physics, 2009,228(15):5410-5436.

[5] Tadmor E. Numerical viscosity and the entropy conditions for conservative difference schemes[J]. Mathematics of Computation, 1984,43(168):369-381.

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

[7] Tadmor E. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[J]. Acta Numerica, 2003,12:451-512.

[8] Fjordholm U S, Mishra S, Tadmor E. Energy preserving and energy stable schemes for the shallow water equations[R]. Hong Kong: Foundations of Computational Mathematics, 2008.

[9] Gottlieb S, Shu C W, Tadmor E. High order time discretizations with strong stability properties[J]. SIAM Review, 2001,43(1):89-112.