程曉晗,封建湖,聶玉峰
(1.西北工業大學應用數學系,陜西 西安 710072; 2.長安大學理學院,陜西 西安 710064)
許多物理模型都采用雙曲型守恒律方程描述,如氣體動力學中的Euler方程、淺水動力學中的淺水波方程和等離子物理學中的磁流體方程等。在一維情況下,有如下形式:
(1)

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

(2)
(3)
與F相容。該格式被稱為熵守恒格式,具有二階精度[2]。
1.1.2熵穩定格式
根據熵穩定條件,在激波等間斷位置,總熵應該耗散(減少)。E.Tadmor[2]指出,一個三點格式只需含有比熵守恒格式更多的黏性則是熵穩定的。因此,P.L.Roe[3]提出用Roe格式的數值黏性項作為對熵守恒格式的補充,數值通量為:
(4)
該格式保持了Roe格式對間斷的良好捕捉效果,并且無需進行熵修正。
1.1.3熵相容格式
考慮黏性形式的Roe格式數值通量[5]:
(5)
(6)
式中:Δai+1/2=f′ui+1-f′ui,取α=1/6。式(6)定義的數值通量跨過激波時具有足夠熵增且數值解保持單調,具有這種性質的通量稱為熵相容通量。熵相容格式是目前估計熵的變化最精確的一種熵穩定格式。
1.1.4WENO型熵相容格式

(7)
氣體動力學中的Euler方程為:
(8)


1.2.1熵守恒格式
對于守恒律系統,一般采用逐段分解的構造方法[7]獲得熵守恒通量。在該方法中,保持總熵不變的數值通量可以表示為沿熵變量空間內不同路徑的積分和的形式。雖然該方法適用于各種守恒律系統,但其計算量過大并且有時候會產生數值不穩定[8]。對于Euler方程和淺水波方程,根據其特定的熵函數,他們的熵守恒通量有比較簡單的構造方法。
其中,Euler方程的熵守恒通量的計算方法如下。定義參數變量z為:
(9)
則熵守恒通量為:
(10)
式中:

1.2.2熵穩定格式
與標量守恒律一樣,對熵守恒格式添加適當的黏性,可使之滿足離散熵不等式。考慮用Roe格式的數值黏性項作為補充,得到關于Euler方程的熵穩定格式,數值通量為:
(11)

1.2.3熵相容格式
類似標量守恒律,為了抵消解在跨過激波時產生的熵增,需要對D進行修正,令:
D1=Λ+αΔΛS
(12)
式中:ΔΛ=diagui+1-ai+1-ui-ai,ui+1-ui,ui+1+ai+1-ui+ai,于是可以得到Euler方程的熵相容格式,數值通量為:
(13)
1.2.4WENO型熵相容格式

(14)
對于一維守恒律方程組,本文中僅以氣體動力學中的Euler方程為例,設計WENO型熵相容格式,其他方程可以參照Euler方程做一些相應的修改,在此不再贅述。
采用半離散格式:
(15)
時間上的推進一律采用三階強穩定的Runge-Kutta方法[9]。為了比較,也給出了相應算例的熵相容格式(EC格式)的計算結果。
在區域x∈[-1,1]上求解初值問題:
采用周期邊界條件,計算到t=10。
本算例可以用來測試格式的收斂性和精度,表1給出了相應的誤差L1和L,充分說明了WENO型熵相容格式(EC-WENO格式)的精度能達到五階。

表1 EC-WENO格式的數值精度Table 1 Numerical accuracy of EC-WENO scheme
在區域x∈[-1,1]上求解初值問題:
采用周期邊界條件,取CFL系數為0.8,空間網格數為40,計算到t=0.3。
本算例的解主要包括2部分:在左邊x=-1/3處由-1到1的跳躍產生了一個稀疏波,在右邊x=1/3處由1到-1的跳躍產生了一個定常激波。計算結果如圖1所示。由圖可知,EC-WENO格式算得的解不僅準確地捕捉到稀疏波,而且在間斷位置也沒有產生偽振蕩,符合“高分辨率”的特性。

圖1 無黏Burgers方程間斷初值問題Fig.1 Discontinuous initial value problem of Burgers equation
在區域[-0.5,0.5]上求解初值問題:

采用Neumann邊界條件,取CFL系數為0.4,空間網格數為100,計算到t=0.1。
精確解包括稀疏波、接觸間斷和激波。密度的計算結果如圖2所示。由圖可知,EC-WENO格式比EC格式能對解的結構有更準確地捕捉,尤其是對激波的捕捉非常銳利。

圖2 一維Euler方程Sod激波管問題Fig.2 Sod shock tube problem of 1D Euler equation
在區域[-0.5,0.5]上求解初值問題:

采用Neumann邊界條件,取CFL系數為0.4,空間網格數為100,計算到t=0.05。
精確解包括2個強稀疏波和1個平凡的接觸間斷。由于中間部分的壓力非常小,接近于真空狀態,很多數值方法容易在該位置出現壓力為負的情況(如Roe格式),導致計算失敗。壓力的計算結果如圖3所示。由圖可知,EC-WENO格式能成功地計算低密度流問題,且效果較好,能有效避免非物理現象(如負壓力)。

圖3 一維Euler方程低密度流問題Fig.3 Low density problem of 1D Euler equation
在區域[-10,15]上求解初值問題:

采用Neumann邊界條件,取CFL系數為0.4,空間網格數為100,計算到t=0.01。
精確解包括強稀疏波、接觸間斷和激波。密度的計算結果如圖4所示。由圖可知,EC-WENO格式取得了非常良好的效果,解在稀疏波有一個平滑的過渡,且在波頭波尾位置都非常貼近精確解。

圖4 一維Euler方程強稀疏波問題Fig.4 Density of strong expansion problem of 1D Euler equation
以熵相容格式為基礎,通過在單元交界面處進行WENO重構,構造了一類求解雙曲守恒律方程的WENO型熵相容格式,有效地提高了熵相容格式的精度。數值模擬結果表明,新格式具有高精度、基本無振蕩性等特點。
[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.