一種求解歐拉方程的高精度數(shù)值格式

任琴琴,白曉雅,鄭秋亞,梁益華

(1.長(zhǎng)安大學(xué) 理學(xué)院,陜西 西安 710064; 2.中國(guó)航空計(jì)算技術(shù)研究所 航空氣動(dòng)力數(shù)值模擬重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710068)

計(jì)算流體力學(xué)(computational fluid dynamics,CFD)需要借助計(jì)算機(jī)完成,隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展,CFD的理論基礎(chǔ)和實(shí)際應(yīng)用得到了快速發(fā)展。雙曲守恒律方程經(jīng)常被應(yīng)用在流體力學(xué)等實(shí)際問題中。眾所周知,即使守恒律方程的初始條件是光滑的,方程的解也可能產(chǎn)生間斷,如激波、接觸間斷波等。激波捕捉的計(jì)算格式是計(jì)算流體力學(xué)的常用工具,這些格式常用于模擬流體力學(xué)中出現(xiàn)的各種現(xiàn)象;對(duì)大量格式研究發(fā)現(xiàn),高階激波捕捉格式具有較好的計(jì)算效率。目前,加權(quán)本質(zhì)無振蕩(weighted essentially non-oscillatory,WENO)格式是一類應(yīng)用廣泛的激波捕捉格式,是近幾年求解雙曲守恒律方程常用的數(shù)值方法之一。

文獻(xiàn)[1-2]提出了總變差遞減(total variation decreasing Runge-Kutta,TVD-Runge-Kutta)的概念;文獻(xiàn)[3]證明了TVD-Runge-Kutta格式在接近光滑極值的情況下最多具有一階精度;文獻(xiàn)[4-5]導(dǎo)出了具有放寬TVD-Runge-Kutta條件和允許在數(shù)值格式中出現(xiàn)偽振蕩性質(zhì)的高階格式,被稱為本質(zhì)無振蕩(essentially non-oscillatory,ENO)格式;文獻(xiàn)[6]在有限差分框架下構(gòu)造ENO格式;文獻(xiàn)[7]提出了有限體積框架下的WENO格式;文獻(xiàn)[8]在有限差分框架下提出了經(jīng)典的WENO格式,簡(jiǎn)稱WENO-JS,該格式在多維空間上建立了局部光滑因子和非線性權(quán)重;文獻(xiàn)[9]對(duì)WENO-JS進(jìn)行分析發(fā)現(xiàn),在臨界點(diǎn)處并未達(dá)到最優(yōu)精確階,利用映射加權(quán)公式提出了一種改進(jìn)的格式,即WENO-M格式;文獻(xiàn)[10]進(jìn)一步改進(jìn),提出了一種含有全局光滑因子的非線性權(quán)重,構(gòu)造WENO-Z格式,數(shù)值實(shí)驗(yàn)表明這種格式主要是對(duì)傳統(tǒng)的WENO-JS格式中非線性權(quán)重的改進(jìn),減小了數(shù)值耗散,同時(shí)提高了分辨率,避免了像WENO-M格式在臨界點(diǎn)附近精度的損失;文獻(xiàn)[11]將這種思想推廣到高階格式,提出了一種全局光滑因子的廣義公式;文獻(xiàn)[12-13]通過引入一種新的基于拉格朗日插值的局部光滑因子η,構(gòu)造WENO-η格式;文獻(xiàn)[14]基于L1范數(shù)提出了用于五階WENO格式的新光滑因子,稱為WENO-NS格式,該格式將一階導(dǎo)數(shù)的高階近似引入到全局光滑因子,相對(duì)于其他WENO格式,WENO-NS格式的全局光滑因子是由五點(diǎn)模板和中間三點(diǎn)模板的光滑因子的平均得到的;文獻(xiàn)[15]導(dǎo)出了七階WENO-NS格式,并用廣義不可分差分算子得到了光滑因子,通過引入?yún)?shù)ζ來平衡光滑區(qū)域與間斷區(qū)域之間精度,由此得到的全局光滑度指標(biāo)滿足在存在臨界點(diǎn)的情況下得到最優(yōu)收斂速度階的充分條件。

迎風(fēng)格式因其物理意義清晰、分辨率高,至今仍是空間離散的主要格式之一。為了達(dá)到較高的計(jì)算效率,并且在間斷處達(dá)到高分辨率,文獻(xiàn)[16-18]提出對(duì)流迎風(fēng)分裂(advection upstream splitting method,AUSM)格式;文獻(xiàn)[19-20]提出對(duì)流迎風(fēng)和分壓(convective upwind and split pressure schemes,CUSP)格式;文獻(xiàn)[21-24]表明,CUSP格式根據(jù)其對(duì)流項(xiàng)的總能量是總焓還是總能，分為總焓對(duì)流迎風(fēng)和分壓(total enthalpy convective upwind and split pressure schemes,H-CUSP)格式、總能對(duì)流迎風(fēng)和分壓(total energy convective upwind and split pressure schemes,E-CUSP)格式。2種格式計(jì)算過程中均避免了矩陣運(yùn)算,提高了計(jì)算效率。E-CUSP格式數(shù)值耗散小,分辨率高,能清晰地捕捉激波剖面和精確的接觸間斷。

本文在空間上采用基于E-CUSP格式與七階的WENO格式耦合得到新的格式,捕捉間斷,實(shí)現(xiàn)高精度,時(shí)間上采用文獻(xiàn)[25]提出的四階TVD-Runge-Kutta格式來推進(jìn)。

1 控制方程

1.1 歐拉方程

考慮一維雙曲守恒律方程:

(1)

1.2 E-CUSP格式

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

(2)

E-CUSP格式的主要思想就是將通量分解為守恒型通量和壓力型通量?jī)身?xiàng)。該格式值耗散項(xiàng)的系數(shù)不是矩陣而是標(biāo)量,避免了矩陣運(yùn)算所帶來的時(shí)間消耗,可以減少計(jì)算量,提高計(jì)算效率。(1)式中的F雅克比矩陣可表示為:

其中

由F和A之間的齊次性可以得到如下關(guān)系:

F=TΛT-1U

(3)

根據(jù)雅克比矩陣的特征值的特性,將(3)式拆分為:

其中,Fc=u[ρ,ρu,E]T、Fp=[0,p,pu]T分別為守恒型通量和壓力型通量。

守恒型通量Fc可表示為:

(4)

在壓力型通量Fp中,動(dòng)量方程中壓力的表達(dá)式為:

(5)

能量方程中的壓力可分解為:

(6)

最終E-CUSP格式中通量F可表示為:

(7)

其中:a為音速;M為馬赫數(shù)。

2 數(shù)值方法

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

WENO格式用于評(píng)估變量UL和UR。采用qk的凸組合近似變量UL的WENO格式可以寫成:

(8)

其中：qk(k=0,1,…,r)為不同模板中變量的重構(gòu)多項(xiàng)式；ωk(k=0,1,…,r)為權(quán)重。

(8)式給出的重構(gòu)形式既適合光滑流場(chǎng)也適合含間斷流場(chǎng)。對(duì)于含激波的間斷流場(chǎng),(8)式中非線性權(quán)函數(shù)ωk的定義如下:

(9)

文獻(xiàn)[8]提出的光滑因子ISk為:

(10)

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

文獻(xiàn)[10]通過引入全局光滑因子τ,重新定義了WENO-JS格式的非線性權(quán)值,并得到接近理想權(quán)值Ck的非線性權(quán)值:

(11)

其中

對(duì)于七階WENO-Z7格式,全局光滑因子定義為τ7=|IS0+3 IS1-3 IS2-IS3|。

2.3 WENO-NS格式數(shù)值重構(gòu)

WENO格式的基本要素是光滑因子的計(jì)算,文獻(xiàn)[15]基于L1范數(shù)光滑因子建立了七階高精度的WENO格式,稱為WENO-NS7格式。因?yàn)榛贚1范數(shù)的光滑因子在光滑區(qū)域可能會(huì)導(dǎo)致精度損失,所以該格式利用導(dǎo)數(shù)的高階近似。

光滑因子的構(gòu)造如下。

定義運(yùn)算符Ls,kU,通過估計(jì)導(dǎo)數(shù)的近似值大小,測(cè)量每一個(gè)四點(diǎn)模板

Sk(j)={xj+k-3,xj+k-2,xj+k-1,xj+k},

k=0,1,2,3,

中解的規(guī)律性。得到這些算符后,光滑因子ISk的計(jì)算公式為:

ISk=ξ|L1,kU|+|L2,kU|+|L3,kU|

(12)

其中:ξ∈(0,1]為平衡光滑區(qū)域和間斷區(qū)域精度之間的平衡的參數(shù);算子Ls,kf為f的廣義不可分差分,公式如下:

(13)

運(yùn)算符L3,kU與WENO-JS7格式中的相同,但是L3,kU在WENO-NS7格式中使用絕對(duì)值,在WENO-JS7格式中使用平方值。

算子Ls,kU的優(yōu)點(diǎn)是導(dǎo)數(shù)ΔxsU(s)在點(diǎn)xj+1/2處的逼近具有更高的精度:

WENO-NS7格式的非線性權(quán)重定義為:

(14)

其中

對(duì)于WENO-NS7格式,全局光滑因子定義為τ7=(IS0-3 IS1+3 IS2-IS3)2。

2.4 Runge-Kutta方法

對(duì)于歐拉方程在時(shí)間方向上的離散,本文使用文獻(xiàn)[25]的四階TVD-Runge-Kutta方法來解決等式:

(15)

其中,L(u)為導(dǎo)數(shù)-f(u)x的近似。

四階TVD-Runge-Kutta方法公式為:

(16)

其中:Δt為時(shí)間步長(zhǎng);L(·)為差分算子。

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

為了驗(yàn)證本文E-CUSP-WENO-NS7格式數(shù)值模擬效果,針對(duì)具有Riemann初值問題的Euler方程,應(yīng)用E-CUSP格式、E-CUSP-WENO7格式、E-CUSP-WENO-Z7格式和E-CUSP-WENO-NS7格式進(jìn)行數(shù)值模擬,對(duì)比它們的結(jié)果,并將數(shù)值模擬結(jié)果與理論解進(jìn)行比較,從而對(duì)其性能進(jìn)行分析。

3.1 一維Euler方程修正的Sod激波管問題

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

采用Neumann邊界條件:

|λi(G(Δx,k))|≤1+cΔx

(17)

其中:λi(G(Δx,k))為增長(zhǎng)矩陣G(Δx,k)的特征值;Δx為空間步長(zhǎng);c為常數(shù)。取CFL=0.1,空間取等距的200個(gè)網(wǎng)格點(diǎn)。

該問題的精確解主要由稀疏波、接觸間斷和激波構(gòu)成,密度、速度和壓力在t=0.2時(shí)用E-CUSP格式、E-CUSP-WENO7格式、E-CUSP-WENO-Z7格式和E-CUSP-WENO-NS7格式計(jì)算一維Euler方程激波管問題的數(shù)值模擬結(jié)果如圖1所示。通過對(duì)圖1及其局部放大圖對(duì)比分析,可以看出,E-CUSP-WENO-NS7格式很好地捕捉到了激波和接觸間斷波。

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

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

采用Neumann邊界條件(17),取CFL=0.1,空間取等距的200個(gè)網(wǎng)格點(diǎn)。密度、速度和壓力在t=0.16時(shí)用E-CUSP格式、E-CUSP-WENO7格式、E-CUSP-WENO-Z7格式和E-CUSP-WENO-NS7格式計(jì)算一維Euler方程Lax激波管問題的數(shù)值模擬結(jié)果如圖2所示。通過對(duì)圖2及其局部放大圖對(duì)比分析,可以看出,E-CUSP-WENO-NS7格式在激波處的數(shù)值模擬結(jié)果優(yōu)于其他格式,對(duì)激波的捕捉更為陡峭,具有較小的數(shù)值耗散。

圖2 一維Euler方程Lax激波管問題的數(shù)值解

3.3 一維Euler方程激波管問題

本文在區(qū)域[0, 1]上求解初值問題：

(ρ,u,p)=

采用Neumann邊界條件(17),取CFL=0.1, 空間取等距的200個(gè)網(wǎng)格點(diǎn)。

密度、速度和壓力在t=0.25時(shí)用E-CUSP格式、E-CUSP-WENO7格式、E-CUSP-WENO-Z7格式和E-CUSP-WENO-NS7格式計(jì)算一維Euler方程Lax激波管問題的數(shù)值模擬結(jié)果如圖3所示。

圖3 一維Euler方程激波管問題的數(shù)值解

通過對(duì)圖3中各局部放大圖對(duì)比分析,可以看出,E-CUSP-WENO-NS7格式在激波處過渡帶窄,而且對(duì)解的捕捉比其他格式更加精確,計(jì)算結(jié)果更接近理論值。

4 結(jié) 論

本文提出E-CUSP與WENO-NS7耦合格式:E-CUSP-WENO-NS7格式,逼近非線性雙曲守恒定律的解,基于模板上插值多項(xiàng)式導(dǎo)數(shù)的不可分差分格式,構(gòu)造局部和全局光滑因子。理論分析和數(shù)值實(shí)驗(yàn)結(jié)果表明,在相同網(wǎng)格情況下,與其他格式相比,E-CUSP-WENO-NS7格式對(duì)接觸間斷和激波的捕捉更為精確,分辨率更高,計(jì)算精度更高,數(shù)值耗散更小,新格式具有更好的穩(wěn)定性。

接下來研究的方向是:將耦合的新格式與自適應(yīng)階結(jié)合提高計(jì)算效率,并應(yīng)用到文獻(xiàn)[26]所提到的海底輸油管道石油泄漏和文獻(xiàn)[27]提到的通風(fēng)空調(diào)等實(shí)際問題中。