五階高分辨率熵穩(wěn)定算法

呂夢(mèng)迪,鄭素佩，陳 芳

(長(zhǎng)安大學(xué) 理學(xué)院, 陜西 西安 710064)

0 引言

雙曲守恒律方程是計(jì)算流體力學(xué)中一類(lèi)反映物理現(xiàn)象及其規(guī)律的方程.一維情況下，它們都具有如下形式：

ut+f(u)x=0,

(1)

其中：x∈R，t>0，u∈Rn是守恒變量，f(u)是通量函數(shù).在求解方程(1)時(shí)，即使初始條件充分光滑，其解在時(shí)間推進(jìn)的某個(gè)時(shí)刻也可能出現(xiàn)間斷，產(chǎn)生激波.出現(xiàn)的間斷解違背了古典解理論，為此LAX于1954年提出了弱解概念[1]，但是弱解并不是唯一確定的，據(jù)此LAX證明了,如果弱解u滿足熵穩(wěn)定條件

E(u)t+F(u)x≤0,

(2)

其中E(u)是u的一個(gè)凸函數(shù)，即E″(u)>0，且F(u)滿足

F′(u)T=E′(u)Tf′(u),

(3)

則u是方程(1)唯一且有物理意義的解.稱(chēng)E為熵函數(shù)，F(xiàn)為熵通量，(E,F)為熵對(duì).方程不同，熵對(duì)的具體表達(dá)式也可能不同.為了滿足熵穩(wěn)定條件，人們常用的方法是添加數(shù)值黏性，但添加多少黏性則依具體情況而定.為了確定熵穩(wěn)定和數(shù)值黏性之間的關(guān)系，TADMOR[2]于1987年定義了一類(lèi)二階熵守恒格式，還證明了一個(gè)三點(diǎn)格式只需包含比熵守恒格式更多的黏性則是熵穩(wěn)定的.之后TADMOR又提出的逐段分解顯式構(gòu)造方法[3]適用于各種守恒型系統(tǒng).2006年，ROE將他本人提出的Roe格式數(shù)值黏性項(xiàng)添加到熵守恒格式上，從而得到熵穩(wěn)定格式[4]，該格式不僅保留了Roe格式對(duì)激波的良好捕捉效果，而且避免了間斷處激波的產(chǎn)生，但是該格式在空間方向上只能達(dá)到一階精度.RAY[5]、封建湖[6-8]也對(duì)此類(lèi)格式進(jìn)行深入探究，構(gòu)造了一系列熵穩(wěn)定格式，但其高階格式還不夠穩(wěn)定.簡(jiǎn)而言之，熵穩(wěn)定格式目前還存在著精度不夠高、高階格式穩(wěn)定性差的問(wèn)題.

據(jù)此，本文通過(guò)在空間方向采用五階CWENO重構(gòu)，得到一種高精度中心加權(quán)基本無(wú)振蕩(Central Weighted Essentially Non-Oscillatory, CWENO)型熵穩(wěn)定格式，時(shí)間上的推進(jìn)采用四階強(qiáng)穩(wěn)定的Runge-Kutta方法[9].最后，將該算法應(yīng)用于若干典型問(wèn)題的數(shù)值求解中，以此來(lái)驗(yàn)證新算法的性能.

1 空間離散

則守恒型半離散格式的一般表達(dá)式為：

(4)

式中G是數(shù)值通量函數(shù).其對(duì)應(yīng)的熵穩(wěn)定條件由此變?yōu)槭?2)的半離散形式.

1.1 熵守恒格式

(5)

(6)

與F相容.該格式被稱(chēng)為熵守恒格式，且具有二階精度.

1.1.1 一維無(wú)黏Burgers方程

對(duì)于無(wú)黏Burgers方程

ut+f(u)x=0，

其中f(u)=u2/2，可取熵函數(shù)E(u)=u2/2，相應(yīng)的熵變量v(u)=u，易得熵通量F(u)=u3/3和熵勢(shì)ψ(v)=u3/6.由式(5)可得無(wú)黏Burgers方程的熵守恒通量為：

(7)

1.1.2 一維Euler方程組

對(duì)于氣體動(dòng)力學(xué)Euler方程組

ut+f(u)x=0,

其中：u=(ρ,ρu,ρE)T是守恒型向量，f(u)=(ρu,ρu2+p,ρuH)T是通量函數(shù).ρ、u、p、E分別為氣體密度、速度、壓強(qiáng)、總能，且E=ρ(1/2u2+e)，比內(nèi)能e=p/{ρ(γ-1)}(其中γ為理想氣體常數(shù)，通常取γ=1.4)，總焓H=E+p/ρ.一般情況下，取一維Euler方程組的一組熵對(duì)為

E(u)=-ρS/(γ-1)，F(xiàn)(u)=-ρuS/(γ-1)，

S是物理熵，且S=lnp-γlnρ,則對(duì)應(yīng)的熵變量v=[(γ-S)/(γ-1)-ρu2/(2p),ρu/p,-ρ/p]T，熵勢(shì)ψ(u)=ρu.

熵守恒通量函數(shù)滿足式(5)，這里的守恒通量函數(shù)有三個(gè)分量需確定，然而僅有一個(gè)方程，因此就需要先確定一些量.本文采用Roe的熵守恒數(shù)值通量[10]，進(jìn)而得到另外一些量的表達(dá)式.其熵守恒通量表達(dá)式如下：

(8)

其中

變量左右狀態(tài)的算術(shù)平均值表達(dá)式為：

對(duì)數(shù)平均值表達(dá)式為：

L(x,y)=(x-y)/{ln(x)-ln(y)}.

1.2 熵穩(wěn)定格式

由熵穩(wěn)定條件可知，在激波等間斷處，總熵會(huì)耗散.因此，若對(duì)熵守恒格式加入適當(dāng)?shù)酿ば?，使之滿足離散熵不等式，即能避免偽振蕩的產(chǎn)生.此格式被稱(chēng)為熵穩(wěn)定格式.

1.2.1 一維無(wú)黏Burgers方程

2006年，ROE將其自己提出的Roe格式黏性項(xiàng)添加到熵守恒格式上，得到熵穩(wěn)定通量函數(shù)為：

(9)

因此，Burgers方程的熵穩(wěn)定通量函數(shù)為：

(10)

式(10)保持了Roe格式對(duì)激波的良好捕捉效果.

1.2.2 一維Euler方程組

根據(jù)Roe和Ismail的熵增定義，ROE[4]于2006年得到一維Euler方程組的熵穩(wěn)定通量函數(shù)為：

(11)

式中R是Euler方程組相對(duì)應(yīng)的雅克比矩陣的右特征向量所構(gòu)成的均值矩陣，

1.3 五階熵穩(wěn)定格式

(12)

(13)

2 數(shù)值算例

算例1無(wú)黏Burgers方程連續(xù)初值問(wèn)題

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

采用周期性邊界條件，取空間網(wǎng)格數(shù)為40，CFL=0.2，其中振幅u1=0.5，速度u0=0，分別計(jì)算到t=0.32和t=0.96.圖1中的ROE[16]算法在空間、時(shí)間方向的離散方法均與本文相同.

(a) t=0.32時(shí)刻EC、Roe和ES通量結(jié)果;(b) t=0.32時(shí)刻熵穩(wěn)定通量結(jié)果;(c) t=0.96時(shí)刻EC、Roe和ES通量結(jié)果;(d) t=0.96時(shí)刻熵穩(wěn)定通量結(jié)果

在激波產(chǎn)生前，熵穩(wěn)定格式的數(shù)值解比熵守恒格式和Roe格式的數(shù)值解更接近精確解；三階、四階、五階CWENO型熵穩(wěn)定格式都可以捕捉到精確解(見(jiàn)圖1(a)、(b)).但在激波產(chǎn)生后，熵守恒格式產(chǎn)生了明顯的偽振蕩；Roe格式和熵穩(wěn)定格式?jīng)]有產(chǎn)生振蕩，且熵穩(wěn)定格式的數(shù)值解比Roe格式的數(shù)值解更貼近精確解(見(jiàn)圖1(c)).由圖1(d)發(fā)現(xiàn)：五階CWENO型熵穩(wěn)定格式比其他格式對(duì)解的捕捉更銳利.本算例通過(guò)激波產(chǎn)生前后解的變化情況充分說(shuō)明了新算法的數(shù)值解收斂于物理解.

算例2無(wú)黏Burgers方程間斷初值問(wèn)題

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

采用周期性邊界條件，取空間網(wǎng)格數(shù)為80，CFL=0.2，計(jì)算到t=0.32.

(a) ES通量結(jié)果;(b) Roe通量結(jié)果; (c) ES-CWENO5通量結(jié)果;(d) 總熵變化;(e) 熵穩(wěn)定通量結(jié)果

由圖2(a)、(b)、(c)可以看出：熵穩(wěn)定算法在稀疏波處產(chǎn)生跳躍；五階CWENO型Roe算法產(chǎn)生了膨脹激波；而五階CWENO型熵穩(wěn)定算法不僅能夠準(zhǔn)確地捕捉到稀疏波，且在間斷處也沒(méi)有出現(xiàn)振蕩；圖2(d)表明：三階、四階和五階CWENO型熵穩(wěn)定算法總熵持續(xù)耗散，其變化均滿足總熵耗散規(guī)律; 由圖2(e)可知：三階、四階、五階CWENO型熵穩(wěn)定算法是在保持Roe格式對(duì)激波良好捕捉效果的同時(shí)避免了間斷處膨脹激波的產(chǎn)生，都得到了數(shù)值解與精確解相吻合的結(jié)果；由兩個(gè)局部放大圖可以明顯地看出五階算法的結(jié)果優(yōu)于三階和四階算法的結(jié)果，即五階CWENO型熵穩(wěn)定算法的分辨率高于三階和四階算法的分辨率.

算例3 Sod激波管問(wèn)題

在區(qū)域x∈[0,1]上求解初值問(wèn)題

采用紐曼邊界條件，取空間網(wǎng)格數(shù)為100，CFL=0.2，計(jì)算到t=0.25.壓強(qiáng)的計(jì)算結(jié)果如圖3所示.由圖3(a)、(b)可知：五階CWENO型熵穩(wěn)定算法的數(shù)值結(jié)果比熵穩(wěn)定算法更貼近精確解，也沒(méi)有出現(xiàn)大的抹平現(xiàn)象和膨脹激波；圖3(c)表明了五階CWENO型熵穩(wěn)定算法總熵的耗散量比三階算法和四階算法少;圖3(d)顯示：五階算法的分辨率比三階和四階算法的分辨率高.

算例4強(qiáng)稀疏波問(wèn)題

在區(qū)域x∈[-10,15]上求解初值問(wèn)題

采用紐曼邊界條件，取空間網(wǎng)格數(shù)為100，CFL=0.2，計(jì)算到t=0.01.壓強(qiáng)的計(jì)算結(jié)果如圖4所示.由圖4(a)、(b)可知：新算法結(jié)果比熵穩(wěn)定算法更接近精確解；圖4(c)顯示了算法總熵變化均滿足總熵耗散規(guī)律;圖4(d)表明：五階算法的分辨率優(yōu)于三階和四階算法的分辨率.

(a) ES通量結(jié)果;(b) ES-CWENO5通量結(jié)果;(c) 總熵變化;(d) 熵穩(wěn)定通量結(jié)果

(a) ES通量結(jié)果; (b) ES-CWENO5通量結(jié)果; (c)總熵變化 ;(d) 熵穩(wěn)定通量結(jié)果

3 結(jié)語(yǔ)

本文針對(duì)一維雙曲守恒律方程(組)的數(shù)值求解問(wèn)題，構(gòu)造五階CWENO型熵穩(wěn)定算法，并把新算法應(yīng)用于Burgers方程和Euler方程組的數(shù)值求解中.通過(guò)對(duì)所得結(jié)果的比較和分析，得出如下結(jié)論：1)新算法所得數(shù)值結(jié)果與精確解吻合效果非常好；2)新算法的分辨率比三階、四階CWENO型熵穩(wěn)定算法的分辨率高；3)新算法沒(méi)有出現(xiàn)大的抹平現(xiàn)象和膨脹激波，且能有效抑制偽振蕩的產(chǎn)生.

新算法在單元交界面進(jìn)行五階CWENO重構(gòu)，時(shí)間方向采用四階強(qiáng)穩(wěn)定的Runge-Kutta方法，具有精度高、穩(wěn)定性強(qiáng)的特點(diǎn).同時(shí)也易于計(jì)算機(jī)編程實(shí)現(xiàn)及向高維雙曲守恒律方程推廣.