一種求解Euler方程的高階格式

任琴琴, 鄭秋亞, 梁益華

(1. 長安大學(xué) 理學(xué)院, 西安 710064； 2. 中國航空計算技術(shù)研究所 航空氣動力數(shù)值模擬重點實驗室, 西安 710068)

在計算流體力學(xué)中, 許多經(jīng)典的數(shù)值方法已被廣泛應(yīng)用. 一個好的數(shù)值方法應(yīng)該在光滑區(qū)域獲得較高的精度, 同時在間斷附近避免非物理意義的偽振蕩. Liu等[1]提出了加權(quán)本質(zhì)無振蕩(WENO)方法, 并構(gòu)造了三階WENO格式. 目前, WENO格式在數(shù)值模擬中已成為求解非線性雙曲守恒律方程常用的數(shù)值方法之一, 但WENO格式計算量大、 數(shù)值耗散偏大. Jiang等[2]在多維空間上建立了具有一般平滑因子和非線性權(quán)重的五階WENO-JS格式; Henrick等[3]對WENO-JS進(jìn)行分析表明, 在臨界點處并未達(dá)到最優(yōu)解, 利用映射函數(shù)提出了WENO-M格式; 文獻(xiàn)[4-6]進(jìn)一步分析改進(jìn)并提出了WENO-Z格式, 這兩種格式主要改進(jìn)了WENO-JS格式中的非線性權(quán)重, 從而在一定程度上降低了耗散, 提高了分辨率; 文獻(xiàn)[7-8]通過引入一種基于Lagrange插值的新局部光滑因子η, 并給出η的顯式表達(dá)式, 構(gòu)造了WENO-η格式.

為提高計算效率, 在間斷處達(dá)到高分辨率, 研究者們又提出了對流迎風(fēng)分裂(AUSM)格式[9]與對流迎風(fēng)和分壓(CUSP)格式[10]. 根據(jù)其物理性質(zhì), AUSM格式將無黏通量分裂為對流項和壓力項, CUSP格式根據(jù)其對流項的物理特征, 分為總焓對流迎風(fēng)和分壓(H-CUSP)格式[11]與總能對流迎風(fēng)和分壓(E-CUSP)格式[12-15]. 兩種格式計算過程中都避免了矩陣運算, 簡化格式的同時提高了計算效率. E-CUSP格式從Jacobi矩陣的特征值出發(fā), 分裂得到對流項和壓力項. 通過將無黏通量分解為守恒型通量和壓力型通量構(gòu)造通量. 該格式擴(kuò)散小, 能清晰捕捉激波和接觸間斷. 本文在空間上采用E-CUSP格式與WENO-η格式耦合得到一種新格式, 時間上采用4階TVD Runge-Kutta方法[16].

1 控制方程

1.1 Euler方程

考慮一維標(biāo)量雙曲守恒律方程:

(1)

其中：U=(ρ,ρu,E)T為守恒變量,F=(ρu,ρu2+p,(E+p)u)T為通量,ρ表示密度,u表示速度,E表示單位體積的內(nèi)能,p=(γ-1)(E-ρu2/2)表示壓強(qiáng),γ=1.4表示比熱比;t表示時間變量；x表示空間變量.

1.2 E-CUSP格式

在區(qū)間單元上利用有限體積法, 可得方程(1)的數(shù)值離散形式為

(2)

E-CUSP格式[15]的主要思想是通量分裂數(shù)值耗散項的系數(shù)不是矩陣而是標(biāo)量, 避免了矩陣運算, 提高了計算效率.

式(1)中F的Jacobi矩陣為

其中，

a為音速,H=(E+p)/ρ為總焓.

F與A之間的齊次性關(guān)系式為

F=TΛT-1U.

(3)

根據(jù)Jacobi矩陣特征值的特性, 將式(3)拆分可得

其中：Fc=u(ρ,ρu,E)T為守恒型通量；Fp=(0,p,pu)T為壓力型通量. 因此，通量F可表示為守恒型通量Fc與壓力型通量Fp之和. 通量Fc可表示為

其中Uc=(1,u,E/ρ)T. 在交界面處質(zhì)量通量的表達(dá)式為(ρu)1/2=(ρLuL+ρRuR), 其中速度uL和uR定義如下：

因此, E-CUSP格式中通量F可表示為

2 數(shù)值方法

2.1 WENO-JS格式數(shù)值重構(gòu)

WENO格式用于評估變量UL和UR. 采用qk的凸組合近似變量UL的WENO格式, 可寫成

(5)

其中:qk(k=0,1,…,r-1)是不同模板中變量的重構(gòu)多項式;γk(k=0,1,…,r-1)是權(quán)重.

(6)

對七階(r=4)WENO-JS格式, 有

2.2 WENO-η格式數(shù)值重構(gòu)

類似上述光滑因子的構(gòu)造過程, Fan等[7]基于Lagrange插值多項式重新定義了一種新的局部光滑因子, 該光滑因子η的簡潔公式為

(7)

(8)

由ηk重構(gòu)的新WENO格式, 記為WENO-η.

因此, 根據(jù)式(7)和式(8)計算出對于七階(r=4)的WENO-η格式, 新光滑因子ηk可表示為

(9)

顯然, 式(9)中ηk的顯式表達(dá)式比七階的WENO-JS格式中ISk的表達(dá)式更緊湊. 通過Lagrange插值多項式Pi,r(x)的Taylor展開式, 可觀察到Lagrange插值多項式有以下遞歸關(guān)系：

(10)

根據(jù)式(10), 易從低階的插值多項式遞歸導(dǎo)出更高階的插值多項式, 利用該性質(zhì)可很大程度簡化更高階光滑因子ηk顯式表達(dá)式的推導(dǎo)過程.

類似文獻(xiàn)[4-6]中WENO-Z格式非線性權(quán)重定義的思想： 利用經(jīng)典的局部光滑因子IS構(gòu)造高階全局光滑指標(biāo)τ2r-1, 使得WENO-Z權(quán)值的凸組合比WENO-JS更接近于中心格式, 從而減少耗散. WENO-η也可構(gòu)造一系列的全局光滑因子.

而且基于η的對稱性, 分析ηk的Taylor級數(shù)展開式可得τ2r-1的截斷誤差為τ2r-1=ο(Δxr+2).

2.3 高階全局光滑指標(biāo)τ2r-1的改進(jìn)

除給出全局光滑性因子τ2r-1的規(guī)則形式外,τ2r-1也可用如下方法[8]表示:

(11)

(12)

其中Fk(k=1,2,…,r-1)的表達(dá)式為

2.4 Runge-Kutta方法

對Euler方程在時間方向上的離散, 本文利用如下4階TVD Runge-Kutta方法[16]:

求解

(13)

其中： Δt表示計算的時間步長；L(·)表示差分算子.

3 數(shù)值結(jié)果與分析

其中:y(xj)為方程的數(shù)值解;y(x)為方程的精確解. 對比不同格式的數(shù)值模擬結(jié)果, 從而分析其性能.

3.1 一維Euler方程激波管問題

例1在區(qū)域[0,1]上求解初值問題:

3.2 一維Euler方程Lax激波管問題

例2在區(qū)域[0,1]上求解初值問題:

圖1 一維Euler方程激波管問題的數(shù)值解Fig.1 Numerical solution of shock tube problem for one-dimensional Euler equation

圖2 一維Euler方程Lax激波管問題的數(shù)值解Fig.2 Numerical solution of Lax shock tube problem for one-dimensional Euler equation

綜上, 首先本文基于Riemann問題解求解一維Euler方程, 在已有WENO格式的基礎(chǔ)上通過重新構(gòu)造基于Lagrange插值多項式的局部光滑因子和全局光滑因子, 改變非線性權(quán)值構(gòu)造過程, 構(gòu)造出高階WENO-η格式. 與其他經(jīng)典WENO格式相比, WENO-η格式收斂階數(shù)更高. 其次, 本文考察了低耗散 E-CUSP格式與高階精度WENO-η格式的耦合性能. 在相同網(wǎng)格條件下, 理論分析和數(shù)值實驗結(jié)果表明, 新格式對接觸間斷和激波的捕捉更精確, 分辨率和計算精度均更高, 數(shù)值耗散更小. 因此, 新格式具有更好的穩(wěn)定性.

表1 不同格式求解例1和例2的相對誤差對比