求解二維雙曲守恒律方程的自適應(yīng)人工粘性熵穩(wěn)定格式

龔承啟

(長安大學(xué)理學(xué)院，陜西 西安 710000)

0 引言

雙曲守恒律方程是流體力學(xué)中的一類重要方程，如歐拉方程、淺水波方程等。數(shù)值求解這類方程的主要困難是，即使初值充分光滑，解也會在有限的時間內(nèi)產(chǎn)生間斷，形成激波。為研究間斷解，Lax于1954年提出了弱解的概念[1]，但弱解并不唯一。1974年Lax又提出熵穩(wěn)定條件[2]，并證明滿足熵穩(wěn)定條件的數(shù)值格式可以得到唯一且滿足物理意義的解。1987年Tadmor構(gòu)造了一類二階熵守恒格式，并證明一個三點格式只須包含比熵守恒格式更多的粘性則是熵穩(wěn)定的[3]。2006年Roe提出在熵守恒格式的基礎(chǔ)上增加Roe格式的數(shù)值粘性項，得到的格式稱為ERoe格式[4]。該格式是目前為止比較理想的熵穩(wěn)定格式，但是一階精度的。為此，羅力、封建湖等人提出在ERoe格式的基礎(chǔ)上加入由Yee構(gòu)造的滿足TVD條件的Combined Superbee限制器來使格式達到高分辨率的效果，即在光滑區(qū)域采用熵守恒格式計算，保持二階精度，在間斷處轉(zhuǎn)用熵穩(wěn)定ERoe格式計算以避免振蕩。該格式稱為EYee格式[4]。2012年，Alexander Kurganov提出了自適應(yīng)人工粘性的方法[5]，該方法通過弱局部剩余量來構(gòu)造人工粘性系數(shù)。這樣構(gòu)造出來的人工粘性系數(shù)在解的光滑處非常小，幾乎不影響格式的精度，在間斷處又足夠大以保持穩(wěn)定。基于此，本文針對二維問題，在熵守恒格式的基礎(chǔ)上加入自適應(yīng)人工粘性，得到一類具有高分辨率特性的熵穩(wěn)定格式。本文將以ERoe格式和EYee格式的在幾個二維問題上的數(shù)值結(jié)果作為參照來說明該格式的特點。

1 熵守恒格式

本文討論二維雙曲守恒律方程

ut+f(u)x+g(u)y=0

(1)

其中u=u(x,y,t)∈RN是守恒變量，(x,y,t)∈R×R×R+，f(u) 是x方向的通量函數(shù)，g(u) 是y方向的通量函數(shù)。

1.1 熵不等式

對于方程組(1)，假定存在凸函數(shù)U:RN→R和函數(shù)F,G:RN→R滿足：

?uF(u)=vT?uf(u),?uG(u)=vT?ug(u)

(2)

其中v為熵變量，定義為

v=Uu(u)

(3)

則方程組(2)的熵解在弱意義下滿足熵不等式

U(u)t+F(u)x+G(u)y≤0

(4)

其中U稱為熵函數(shù)，F(xiàn)和G分別是x方向和y方向的熵通量函數(shù)。

1.2 熵守恒格式

對于求解方程組(1)的半離散守恒型差分格式

(5)

v=Uu(u)

(6)

熵勢為

ψx=vTf-F,ψy=vTg-G

(7)

(8)

對于二維歐拉方程

ut+f(u)x+g(u)y=0

(9)

2 自適應(yīng)人工粘性熵穩(wěn)定格式

2.1 二維局部弱剩余量[5]

(10)

(11)

作為數(shù)值解的估計。并定義局部測試函數(shù)

(12)

(13)

將(11)和(13)式代入(10) 式，得到

(14)

其中

(15)

(16)

其中r是所選用的數(shù)值格式的精度，本文選用的是熵守恒格式，則r=2。

2.2 二維自適應(yīng)人工粘性熵穩(wěn)定格式

在方程(1)的右端增加人工粘性項得到

(17)

(18)

3 數(shù)值算例

本文有三個數(shù)值算例，每個算例都以ERoe格式和EYee格式作為參照。對于半離散形式的格式(18)，在時間上采用三階Runge-Kutta法[8]。

3.1 算例1

在空間區(qū)域[0,1]×[0,1] 上取初值

取0階外推邊界條件，CFL=0.4，Δx=Δy=0.01計算到時間T=0.3。該問題由四個一維激波構(gòu)成[9]。圖1(a)、(b)、(c)分別是ERoe、EYee、EA格式的計算結(jié)果的密度的等高線圖(取30條等高線)。可以看到EYee、EA格式對激波的捕獲明顯好于ERoe格式。在四個激波的交界處，EA格式對細節(jié)的捕獲略差于EYee格式。

圖1

3.2 算例2

取初值

取0階外推邊界條件，CFL=0.4，Δx=Δy=0.01，計算到時間T=0.3。該問題包括一個激波，兩個接觸間斷和一個稀疏波[9]。圖2(a)、(b)、(c)分別是ERoe、EYee、EA格式的計算結(jié)果的密度的等高線圖(取30條等高線)。

3.3 算例3：超音速激波衍射問題[10]

該問題的計算區(qū)域如圖3所示其中黑色區(qū)域(0, 0.05)×(0, 0.5)為擋板，灰色區(qū)域(0, 0.05)×(0.5, 0.9)的初值為

(p1,ρ1,u1,v1)=(30.0594,7.0411,4.0799,0)

其余白色區(qū)域的初值為

(p2,ρ2,u2,v2)=(1,1.4,0,0)

擋板處取反射邊界條件，灰色區(qū)域左側(cè)的邊界值固定為(p1,ρ1,u1,v1)，其余邊界處取0階外推邊界條件，CFL=0.4，Δx=Δy=0.01，計算到時間T=0.1561。圖4(a)、(b)、(c)分別是ERoe、EYee、EA格式的計算結(jié)果的密度的等高線圖(取30條等高線)。

圖2

圖3 超音速激波衍射問題的計算區(qū)域

4 結(jié)語

本文針對二維雙曲守恒律方程，在熵守恒格式的基礎(chǔ)上添加自適應(yīng)人工粘性，使新格式在保持熵穩(wěn)定的同時又具有高分辨率的特征。在捕獲激波方面，該格式明顯優(yōu)于ERoe格式，相對于EYee格式又避免了限制器的使用。但該格式中包含可調(diào)整的參數(shù)C，需要針對不同問題具體調(diào)整，這是該格式的主要缺陷。該參數(shù)通用的選取方法有待于以后進一步探討。

圖4

[1] LAX P D. Weak solutions of non-linear hyperbolic equations and their numerical computations[J]. Communication on Pure and Applied Mathematics,1954, 7(1):159-193.

[2] LAX P D. Hyperbolic systems of conservation laws and the mathematical theory of shock waves[A] //Volume 11 of SIAM Regional Conferences Lectures in Applied Mathematics[C]. Philadelphia: Society for Industrial and Applied Mathematics, 1973:1-48.

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

[4] 羅力,封建湖,唐小娟,等.求解雙曲守恒律方程的高分辨率熵穩(wěn)定格式[J]. 計算物理, 2010, 27(5): 671-678.

[5] KURGANOV A, LIU Y. New adaptive artificial viscosity method for hyperbolic systems of conservation laws[J]. Journal of Computational Physics, 2012, 231(24):8114-8132.

[6] FJORDHOLM U, MISHRA S, TADMOR E. Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws[J]. SIAM Journal on Numerical Analysis, 2012, 50(2):544-573.

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

[8] GOTTLIEB S, SHU C W, TADMOR E. Strong stability-preserving high-order time discretizations methods[J]. Society for Industrial and Applied Mathematics Review, 2001, 43(1):89-112.

[9] KURGANOV A, TADMOR E. Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[J]. Numerical Methods for Partial Differential Equations, 2002, 18(5):584-608.

[10] NISHIKAWA H, KITAMURA K. Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers[J]. Journal of Computional Physics, 2008, 227(4):2560-2851.