五階WENO格式求解一維Euler方程

1、五階WENO格式求解一維Euler方程摘 要:應(yīng)用加權(quán)本質(zhì)無(wú)震蕩（WENO）格式耦合低耗散E-CUSP格式，空間采用低耗散E-CUSP迎風(fēng)格式，時(shí)間采 用顯式Runge-Kutta推進(jìn)方法，分別對(duì)一維激波管問(wèn)題和一維低密度流問(wèn)題進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)WENO格式耦 合低耗散E-CUSP算法在接觸間斷和激波處的捕捉效果有了顯著地提高。結(jié)果表明，WENO格式精度高，穩(wěn)定性 能好，將WENO格式耦合低耗散E-CUSP算法應(yīng)用到一維Euler方程的數(shù)值解中，對(duì)接觸間斷和激波的捕捉能力 較強(qiáng)，尤其是對(duì)激波的捕捉能力，僅需要23個(gè)網(wǎng)格單元，為后續(xù)將其推廣到二維問(wèn)題上做了良好的鋪墊。 關(guān)鍵詞：計(jì)算流體力學(xué)；E

2、-CUSP格式;加權(quán)本質(zhì)無(wú)震蕩格式；一維EuDc?方程;低密度流；激波管Fifth Order WENO Scheme for Solving One- dimensional Euler EquationsAbstract: The performance of coupled low dissipation E-CUSP scheme with high precision weighted essentially non - oscillatory WENO scheme is investigated by solving the one- dimensional shock tube

3、 problem and the one - dimensional low density flow problem. With the space adopts the low - dissipation ECUSP upwind scheme and the time adopts the explicit Runge - Kutta propulsion method. It is found that WENO scheme coupled with low dissipation E- CUSP scheme in space. The capture effect of - CU

4、SP algorithm in contact discontinuities and shock waves has been significantly improved. The results show that the WENO scheme has high accuracy and good stability. And WENO scheme coupled with low dissipation E- CUSP algorithm is applied to the numerical solution of one- dimensional Euler equation.

5、 It has strong ability to capture contact discontinuities and shock waves, especially for shock waves, only two or three grid elements are needed. It lays a good foundation for extending it to two- dimensional problems.Key words: computational fluid dynamics； E - CUSP scheme； weighted essentially no

6、n - oscillatory scheme； one - dimensional Euler equation； low density flow； shock tube引言計(jì)算流體力學(xué)（CFD）是伴隨著計(jì)算機(jī)技術(shù)的發(fā)展 而發(fā)展的一門學(xué)科，近年來(lái)得到了快速的發(fā)展，成為流 體力學(xué)中最活躍、最富有生命力的學(xué)科領(lǐng)域。盡管 CFD在眾多領(lǐng)域得到了廣泛應(yīng)用，但是CFD的模擬精 度和計(jì)算效率一直是人們關(guān)注的焦點(diǎn)問(wèn)題。隨著CFD 的快速發(fā)展，出現(xiàn)了許多經(jīng)典的CFD計(jì)算格式和數(shù)值 模擬算法,其中由Roe提出的近似Riemann解格式是 最著名的迎風(fēng)差分格式,一度被認(rèn)為是最精確的差分 格式12，但是Roe格式

7、需要花費(fèi)大量的時(shí)間計(jì)算耗 散矩陣，為了克服Roe格式的這一缺陷,人們研究了很 多不需要矩陣運(yùn)算的迎風(fēng)格式，其中Zha和ShenT 的低耗散E-CUSP迎風(fēng)格式將無(wú)粘通量分解成為守恒 型通量和壓力通量，有效地避免了矩陣運(yùn)算，但是其對(duì) 于間斷尤其是接觸間斷的抹平效應(yīng)太大。因此1994 年Liu等提出了加權(quán)本質(zhì)無(wú)振蕩（WENO）方法，并 且構(gòu)造了一個(gè)三階的WENO格式，詳細(xì)闡述了 WENO 格式的理論原理和其構(gòu)造過(guò)程。1996年Jiang和Shu67又給出了一個(gè)WENO格式的構(gòu)造方法，他們 所給出的方法可以在多維空間上構(gòu)造任一階精度的 WENO格式，并且設(shè)計(jì)了一個(gè)五階的WENO格式。 WENO格式是

8、近幾年來(lái)比較流行的一種高分辨率數(shù)值 方法6 1%，在以對(duì)流為主的對(duì)流擴(kuò)散方程，尤其是雙 曲守恒律方程上應(yīng)用尤為廣泛。本研究將WENO格 式應(yīng)用到一維激波管問(wèn)題的求解過(guò)程中，將計(jì)算結(jié)果 與Zha格式進(jìn)行比較,WENO格式的計(jì)算解與一維激 波管問(wèn)題713的理論解模擬效果較好，表明采用WENO 格式建立的高分辨率模型可以很好地模擬激波管 問(wèn)題。它能夠保證在解的光滑區(qū)精度更高，在解的間 斷區(qū)保持陡峭的間斷過(guò)渡和本質(zhì)無(wú)振蕩性質(zhì)。WENO 格式不同于ENO格式單一的模板選取，它使用所有候 選模板的凸組合，至今，五階精度4,14 17的WENO格 式是使用最為廣泛的。1控制方程控制方程為一維Euler方程，

9、其守恒形式為: + 竺=#(1)a(1)-P -P 1pu為守恒型變量,F =pu + pL pe L ( pe + p) u式中,Q =為通量,p = ( - 1) (p(-32)。其中p是密度，u是速度，( 是比內(nèi)能,3是壓強(qiáng),-是取1.4的比熱值。(2)2數(shù)值方法(2)空間采用低耗散E-CUSP格式，時(shí)間使用四步龍 格庫(kù)塔推進(jìn)方法,對(duì)方程(1)在單元上進(jìn)行積分，利用 有限體積法可得：!Q- = (F 1 F 1 )式中為空間步長(zhǎng)，-表示時(shí)間層，#表示空間節(jié)點(diǎn)； Q-和Q-+I表示Q在t = t-和$ = t-+1時(shí)刻的空間平均 值;F表示穿越單元邊界a+x或x,+x的通量;F# +x和

10、F#+x表示為與通量F相容的數(shù)值通量。Zha的E-CUSP格式的基本思想是通過(guò)將通量分 解為守恒型通量和壓力通量,以避免矩陣運(yùn)算,減少計(jì) 算量，提高計(jì)算效率。其中雅克比矩陣可表示為：FQ-3) u(3) uTAT 1(3)+ ( - 3) upe + p) ( 1) u-u式中:u式中:u1uF4與壓力通量Fp的和。守恒型通量可表示為：F4=； (pu) X () I pu lx (-&*) ( 5)11ua%式中&在交界面處質(zhì)量通量為:可以得到如下關(guān)系:F = TAT 1Q式中:(4)L e(pu) _L = ( PlUL + pRUR ) 其中，+uLu% - a% -F = T%u%T-

11、1 Q + T%L%uL%aT-1Q=al.+ I ml i12+ I:* hl/ + 彳(M* +1)2-里廣-P -P -% -其中F4 = upu,Fp =PL peL pH _=Fc + 32(p/P) (PCP)* +(pCPp) p/p) +(pCPp) (6) 在壓力通量Xp中，動(dòng)量方程中的壓力可表示為：F4表示守恒型通量，特征值為u,F3表示壓力型通量, 特征值為a這樣通量F就被分解成為守恒型通量嚴(yán) + r-： (Mr 1) 2 本文定義如下：2( p/p)lRP( P，P& * + ( P_P)( 7)式中* =：( : 1( 2 ： aM( M2 -1) ,a =+(o能量

12、方程中的壓力可分解為:(PN)+ =； P (N + -* +； P (N - -(8)式中,-=；( -* +)，M* = j,M(8)22最終Zha的E-CUSP格式中X多的的可表示為: 當(dāng) I u I #-,(9)廣 ( pu) 1 ( &* +&(9)-%P+P%(u + - 多)當(dāng) u Z - a,Fx = X;而當(dāng) u a,Fx = X* ,式中 - 是 音速,M是馬赫數(shù)。23 WENO格式2WENO格式的思想是:不僅僅是選擇ENO重構(gòu)模 板中的某一個(gè)模板,而是所有的模板的凸組合。設(shè)其 中 S 個(gè)模板為 S,( #) = A, - r ,A# - r + 1，,A, - r + S

13、-1 ,J =%, 1,#,S-1,可求得！*+土不同的S個(gè)重構(gòu)方式：S-1= ! Cj/U.-f ,r = %,1 ,# ,S - 1( 1%)J = %其中可能會(huì)有以下三個(gè)插值模板：，%= U.-2 ,U.-1 ,U.,1= U.-1 ,U,U+1，2= U,u. + 1 ,u# +2模板如圖1所示：圖1 WENO格式模板選擇圖UJ的表達(dá)如下所示：7-7( %)1 7,11U#+i=T U#-2-6U#-1+6 U#U( +)1 土 U. 1 U. + U. +1(11)U(3) 1 = U- + U-，- U- i1+t 3 i 6. +16(11)WENO格式的思想是通過(guò)多個(gè)模板的線性

14、組合低 階通量得到高階通量。通過(guò)構(gòu)造合適的加權(quán)系數(shù)以及 選取適當(dāng)?shù)墓饣蜃?，可以使得WENO格式的數(shù)值解 在間斷附近保持基本無(wú)振蕩特性，而且能夠在光滑區(qū) 域獲得一致高精度。圖2 VFOA算法流程圖三階ENO格式需要在所有模板中選擇一個(gè)最為 光滑的模板;而五階WENO格式則是將所有模板進(jìn)行 了加權(quán)組合。對(duì)于三階ENO格式，有3個(gè)備選模板, 每個(gè)模板包含3個(gè)格點(diǎn)信息。將它們組合在一起可以 形成包含5個(gè)格點(diǎn)信息的插值區(qū)域。通過(guò)插值，可以 獲得五階精度的WENO格式，插值表達(dá)式可以寫為：2fl+ i! H?3(f.+S-2 J.+S-1，,M)( 12)S=%式中,HS為相對(duì)應(yīng)模板的權(quán)重，與每個(gè)模板的

15、光滑性 有關(guān);尤為模板的狀態(tài)變量;qk為模板的插值多項(xiàng)式, 一般采取Newton插值可得到。將式(12)用于狀態(tài)變 量的插值計(jì)算上，以.方向?yàn)槔瑢?duì)于.+；處的左狀 態(tài)量,插值模板示意圖如圖3所示,?% ,?1和?2分別為 3個(gè)模板。知9102i-2i-1 i /+1i+2圖3左狀態(tài)量插值模板示意圖其加權(quán)重構(gòu)公式為：U+H%?% + H1 ?1 + H2?2( 13)所使用的3個(gè)重構(gòu)模板的表達(dá)式如下(U為原始 變量) ：% S-2U-1 U ?1= -： U.-1+； U. +； U. + 1( 14)?2 =y U. +； U. + 1 -_( U.+2Jiang和Shu的權(quán)因子的定義為:S

16、HSSv M #(S 0,1,2) (15)_ O-(#S + .) 2式中，權(quán)系數(shù)分別為:0% =0. 3，S =0. 6,0 =0. 1,參數(shù) .的作用是避免分母為0,這里取.= 10r6；#s為模板 的光滑度量系數(shù)，用于度量第S個(gè)模板的光滑性。采 用文獻(xiàn)計(jì)算方式,其表達(dá)式為：# = !3 (撰T ()dA(16)I = 1 A-2Oa /具體計(jì)算如下：#2 = 12( !, FT1 + !,T)2 +：(3U# -4!+1+ !T)2#1=13 (U-1 -2U, + U+1)2 +： (U-1- U+1)2#0 = 1!( U#-2 -2!-1+ U#)2 +土( U#-2 -4U#+

17、1 +3!)2時(shí)間離散將式(14)、式(15)和式(16)代入 式(13)中,就可計(jì)算得到具有五階精度的左狀態(tài)量。 對(duì)于處右狀態(tài)量，計(jì)算方法與上文相同，根據(jù)對(duì)稱性便 可計(jì)算得到。4結(jié)果分析算例1: 一維Euler方程低密度流問(wèn)題為了檢驗(yàn)本文格式的穩(wěn)定性，選用一維Euler方 程低密度流問(wèn)題。該問(wèn)題在中間部分的壓強(qiáng)和密度接 近于零，即接近真空，以往一些數(shù)值格式在計(jì)算時(shí)中間 部分的密度時(shí)常會(huì)出現(xiàn)負(fù)值，造成計(jì)算失敗，給數(shù)值模 擬也帶來(lái)一定的難度。一維Euler方程低密度流問(wèn)題 描述如下：(1,-2,0.4),a 0. 5本文在空間區(qū)域0,1上，網(wǎng)格為與文獻(xiàn)11具 有相同網(wǎng)格點(diǎn)數(shù)(200個(gè))的均勻分布

18、網(wǎng)格，計(jì)算中條 件數(shù)CFL =0.05。圖4展示了 $ =0.1 s時(shí)刻本文壓 力、密度和速度與理論值與文獻(xiàn)11中采用Roe格式 計(jì)算結(jié)果的比較。從圖4可以看出，本文格式就壓力、 密度和速度來(lái)說(shuō),均得到了與理論值更接近的計(jì)算結(jié) 果,對(duì)接觸間斷和稀疏波的捕捉效果優(yōu)于文獻(xiàn)11中 的結(jié)果，說(shuō)明本文算法具有較高的數(shù)值模擬精度和良 好的穩(wěn)定性。另外，本文E-CUSP方法所用計(jì)算時(shí)間 比文獻(xiàn)11中Roe格式所用時(shí)間縮短了 2. 7%，在本 算例中,E-CUSP方法相對(duì)于Roe格式具有較高的計(jì) 算效率這一特性并沒(méi)有得到好的體現(xiàn),這是因?yàn)樗?慮問(wèn)題是簡(jiǎn)單的一維無(wú)粘流動(dòng)問(wèn)題，網(wǎng)格規(guī)模小。隨 著問(wèn)題的復(fù)雜度和網(wǎng)

19、格規(guī)模的增大,E-CUSP格式的 計(jì)算高效性才能得到較好的體現(xiàn)。圖4 $ =0.1 s時(shí)的壓力圖算例2: 一維Euler方程激波管問(wèn)題由于一維激波管非定常流動(dòng)包含有激波、接觸間 斷等流動(dòng)特征，網(wǎng)格生成容易且質(zhì)量能夠得到保證，加 上計(jì)算中初始和邊界條件容易處理，常常被人們選作 考察格式性能的理想模型。本文在區(qū)域0,1上求解 初值問(wèn)題(1, -0,1) ,a 0. 5采用Neumann邊界條件，取空間網(wǎng)格點(diǎn)數(shù)400 個(gè),計(jì)算中條件數(shù)CFL =0.3。圖5是$ =0.25 s時(shí)刻 本文密度、馬赫數(shù)、壓力和速度與理論值以及Zha格式 計(jì)算結(jié)果的比較。從圖5可以看出，本文計(jì)算結(jié)果更 接近理論值，新格式捕捉到的接觸間斷和激波的質(zhì)量 明顯優(yōu)于Zha格式的計(jì)算結(jié)果結(jié)果，尤其