柴得林,王 強(qiáng),2,易 賢,2*,劉 宇

(1. 中國空氣動力研究與發(fā)展中心 結(jié)冰與防除冰重點實驗室, 四川 綿陽 621000;2. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000)

隨著電子計算機(jī)技術(shù)的飛速發(fā)展,計算流體動力學(xué)在流體力學(xué)領(lǐng)域的理論研究與工程應(yīng)用中起到越來越大的作用。作為計算流體動力學(xué)的基礎(chǔ)之一,離散格式的性能對流場數(shù)值模擬具有重要影響。特別地,加權(quán)本質(zhì)無振蕩(weighted essentially non-oscillatory,WENO)格式的提出極大地推進(jìn)了含激波等復(fù)雜流動結(jié)構(gòu)的流場的精確數(shù)值模擬。

Liu等[1]在本質(zhì)無振蕩(essentially non-oscillatory, ENO)[2-3]格式的基礎(chǔ)上創(chuàng)造性地提出了WENO 格式,采用子模板上的低階格式的非線性凸組合使得格式兼具高精度與本質(zhì)無振蕩特性,并設(shè)計了有限體積形式的3階、4階精度 WENO 格式。Jiang等[4]對 WENO 格式進(jìn)行理論分析,將WENO格式拓展至有限差分形式,系統(tǒng)地設(shè)計了任意階有限差分形式格式的光滑因子與非線性權(quán)計算方法,他們提出的 5 階WENO 格式成為最經(jīng)典的 WENO 格式之一,一般記為WENO-JS格式。

盡管WENO-JS格式具有優(yōu)越的激波捕捉性能,但仍存在耗散過大,極值點處精度降階等不足。圍繞WENO-JS格式,學(xué)者們開展了大量的性能優(yōu)化研究。在非線性權(quán)計算方面,Henrick等[5]指出5階WENO-JS格式在求解雙曲守恒律時并未完全滿足 5 階精度,在極值點附近發(fā)生降階。為此,他們在非線性權(quán)的計算中引入了一個非線性權(quán)映射函數(shù),設(shè)計了一種完全滿足 5 階精度的改進(jìn)WENO格式,記為 WENO-M 格式。根據(jù) WENO-M 格式的思想,多種新型映射函數(shù)被提出并應(yīng)用于 WENO 格式的優(yōu)化[6-8]。Borges等[9]從增大間斷子模板的權(quán)重分配的角度開展研究,指出增大非線性加權(quán)時間斷子模板的權(quán)重,可降低格式耗散,優(yōu)化格式性能;在這一理論上,他們設(shè)計了一個高階光滑因子,構(gòu)造了新的非線性權(quán),提出了耗散更低、分辨率更高的 5 階 WENO-Z 格式。

在WENO-Z格式的基礎(chǔ)上,Liu等[10]改進(jìn)了5階WENO-Z格式的高階光滑因子及其在非線性權(quán)中的應(yīng)用,使得格式既滿足5階精度充分條件,又具有較低耗散。Castro等[11]給出了高階光滑因子的通用公式,將 WENO-Z 格式拓展至任意奇數(shù)階。Acker等[12]在WENO-Z格式的非線性權(quán)公式中增加了光滑因子比值相關(guān)項,提出WENO-Z+格式,進(jìn)一步提高了格式中間斷子模板的權(quán)重,改善了格式對高頻波的分辨率,并指出在通常的應(yīng)用中,間斷子模板上的權(quán)重對WENO格式實際計算性能的影響起主要因素。文獻(xiàn)[13-16]均引入光滑因子比值優(yōu)化了格式權(quán)重, 文獻(xiàn)[17-20]對光滑因子進(jìn)行了重新設(shè)計與構(gòu)造。徐維錚等[21]則對3階WENO-Z格式的光滑因子進(jìn)行了多種設(shè)計與系統(tǒng)研究,研究表明格式在連續(xù)解非極值點處的理論精度對實際計算性能起決定性的作用,極值點處的精度影響則較小。上述研究表明改進(jìn)非線性權(quán)計算方法可有效實現(xiàn)WENO-Z格式耗散降低,性能提升。

與上述僅改進(jìn)非線性權(quán)計算方法的研究不同,模板優(yōu)化是WENO格式改進(jìn)的另一重要方法。Martín等[22]和HU等[23]通過在WENO構(gòu)造模板中引入下迎風(fēng)模板,優(yōu)化權(quán)重,分別提出了WENO-SYMBO和WENO-CU6格式;Zhu等[24]則創(chuàng)造性地設(shè)計了由一個5點模板和兩個2點模板加權(quán)得到的5 階有限差分 WENO格式,該格式對線性權(quán)的選擇更為靈活,實現(xiàn)更為簡單。這些研究在優(yōu)化模板的基礎(chǔ)上,對非線性權(quán)計算方法,包括線性權(quán)的選取、光滑因子的計算等,進(jìn)行了一定改進(jìn),最終實現(xiàn)了格式性能的提升。

盡管上述研究已經(jīng)實現(xiàn)了WENO格式性能的大幅改善,模板優(yōu)化方法也往往伴隨著非線性權(quán)計算方法的重新設(shè)計,但鮮有研究對非線性權(quán)計算方法改進(jìn)與模板優(yōu)化之間的關(guān)系進(jìn)行研究,鮮有研究提出可以轉(zhuǎn)化為改進(jìn)非線性權(quán)計算方法的模板優(yōu)化方法。

本文以5階WENO-Z格式為研究對象,借鑒文獻(xiàn)[22-24]的優(yōu)化模板思路,在WENO-Z格式構(gòu)造中引入一個由3點模板重新組合形成的4點模板,取其上重構(gòu)格式為對應(yīng)3點模板格式的線性組合,通過這種模板重組方法實現(xiàn)了格式的非線性權(quán)調(diào)節(jié)與性能提升;同時借助子模板的線性組合特性將所提模板優(yōu)化方法等效轉(zhuǎn)化為格式的改進(jìn)非線性權(quán)計算方法;采用一系列基準(zhǔn)問題對改進(jìn)格式性能提升進(jìn)行數(shù)值驗證。

1 WENO-Z格式

以一維雙曲守恒律為研究對象

(1)

(2)

式中,h(x)為數(shù)值通量函數(shù),其隱式定義為

(3)

(4)

5階WENO-JS格式或WENO-Z格式的重構(gòu)模板如圖1所示。

圖1 5階WENO-JS/WENO-Z格式重構(gòu)模板Fig.1 Reconstruction stencils of fifth-order WENO-JS/WENO-Z

S(5)={Ii-2,Ii-1,Ii,Ii+1,Ii+2}為構(gòu)造大模板,可劃分為3個互有重疊的3點子模板,即S0={Ii-2,Ii-1,Ii}、S1={Ii-1,Ii,Ii+1}和S2={Ii,Ii+1,Ii+2}。在各子模板上計算數(shù)值通量對應(yīng)的線性格式,可得

(5)

而大模板上數(shù)值通量對應(yīng)的5階線性格式為

(6)

式中:dm為線性權(quán)重,d0=0.1,d1=0.6,d2=0.3;m=0,1,2。引入與其相對應(yīng)的非線性權(quán)ωm,則可得WENO格式的一般形式

(7)

5階WENO-JS格式的非線性權(quán)計算公式為

(8)

式中,βm為光滑因子,表征相應(yīng)模板上變量的光滑程度,ε為一正數(shù)小量,防止分母為0。文獻(xiàn)[4]給出的k階WENO格式的光滑因子βm通用計算公式為

(9)

對于5階WENO-JS格式各光滑因子為

(10)

而ε則取經(jīng)驗值10-6。

5階WENO-Z格式的非線性權(quán)計算公式為

(11)

式中:βm為式(10)所示的WENO-JS光滑因子;τ5=|β0-β2|為高階光滑因子;p為度量光滑因子對非線性權(quán)的影響的指數(shù),一般取1;εZ為一正數(shù)小量,防止分母為0。參考文獻(xiàn)[9],本文εZ取一極小量10-40使其僅起到防止分母為0的作用。較之WENO-JS格式,WENO-Z格式增大了間斷模板的非線性權(quán),降低了格式耗散,提高了格式分辨率。

2 重組模板WENO-Z格式

2.1 模板重組

圖1中WENO-Z格式計算xi+1/2處的值時所用整體模板是上迎風(fēng)的,為了降低格式耗散,應(yīng)使得模板更接近中心模板,因而本文在WENO格式的構(gòu)造中引入中心模板S3={Ii-1,Ii,Ii+1,Ii+2},如圖2所示。則改進(jìn)格式含有4個子模板,記改進(jìn)格式為WENO-ZF,F為單純記號,則新模板組合下,對應(yīng)加權(quán)公式為

(12)

2.2 線性權(quán)與子模板S3參數(shù)計算

由圖2可知子模板S3為S1、S2的組合,假設(shè)S3上線性格式為S1、S2對應(yīng)線性格式的線性組合,即

圖2 5階WENO-ZF格式重構(gòu)模板Fig.2 Reconstruction stencils of fifth-order WENO-ZF

(13)

(14)

同時有

(15)

由數(shù)值通量線性格式精度式(4)得

(16)

(17)

在單元Ii=[xi-1/2,xi+1/2]上采用與式(15)相同的線性組合,則有

將式(18)代入光滑因子的定義式(9)可得S3上光滑因子計算式為

2.3 非線性加權(quán)與簡化

參考WENO-Z格式的加權(quán)方法,仍將WENO-ZF格式的非線性權(quán)取為

(20)

其中:m=0,1,2,3;參考文獻(xiàn)[9],取τ5=|β0-β2|,p=1,εZF=10-40。

(21)

式(21)中,τ5、βm、εZF與WENO-ZF格式中式(20)中的定義及取值相同。對比式(20)和式(21)可知,重組模板的模板優(yōu)化方案最終等效為非線性權(quán)計算方法的優(yōu)化。為減少計算量,后續(xù)數(shù)值實驗中采用式(21)進(jìn)行WENO-ZF格式的計算。

(22)

將其代入式(18)得

(23)

將式(22)、式(23)代入光滑因子的定義式(9),可得

(24)

因而,計算非線性權(quán)時可取

(25)

這樣可以減少計算量。為便于表示,下文中式(19)和式(25)計算β3的改進(jìn)格式分別記為WENO-ZF1、WENO-ZF2。

2.4 近似色散關(guān)系分析

近似色散關(guān)系(approximate dispersion relation,ADR)分析方法[25]可有效分析非線性格式的譜特性,即色散耗散特性,廣泛應(yīng)用于非線性格式的設(shè)計與優(yōu)化。采用該方法計算所得WENO-ZF1、WENO-ZF2格式的色散耗散特性如圖3所示。圖中Φ(φ)為波數(shù)φ對應(yīng)的修正波數(shù),其實部和虛部分別反映格式的色散與耗散。

比較圖3中各格式色散與耗散特性:低頻波范圍內(nèi),所有格式均有理想的色散耗散;中高頻波范圍內(nèi),WENO-ZF1與WENO-ZF2格式的色散與耗散均小于WENO-Z格式,而WENO-JS格式的色散和耗散則明顯大于這三種格式,這表明采用模板重組引入模板S3改進(jìn)WENO格式的方法可以有效降低格式色散與耗散。對比WENO-ZF1與WENO-ZF2格式的色散與耗散特性,WENO-ZF2格式除對部分中頻波的色散略小于WENO-ZF1格式外,其他情況下對中高頻波的色散和耗散均大于WENO-ZF1格式。

(a) 色散(a) Dispersion

(b) 耗散(b) Dissipation

3 數(shù)值實驗

本節(jié)采用5階WENO-JS、WENO-Z、WENO-ZF1和WENO-ZF2格式對一系列經(jīng)典的數(shù)值算例進(jìn)行模擬,以驗證格式的精度及其色散耗散特性、分辨率等方面的性能。需說明的是,算例中變量如無特別說明,均為無量綱量;各算例的時間離散均采用3階TVD Runge-Kutta 格式計算。

3.1 線性對流方程

對線性對流方程

(26)

取初始條件u(x,0)=sin(πx),在網(wǎng)格數(shù)為20至320的均勻網(wǎng)格上進(jìn)行計算,結(jié)束時間取t=2,計算所得L1、L∞誤差與精度如表1所示。

由表1中數(shù)據(jù)可得,四種WENO格式均可達(dá)到設(shè)計的5階精度;相同網(wǎng)格數(shù)下,WENO-Z、WENO-ZF1、WENO-ZF2三種格式誤差較為接近,較WENO-JS格式小一個數(shù)量級。比較WENO-Z、WENO-ZF1、WENO-ZF2三種格式誤差,對于L∞誤差,相同網(wǎng)格下,WENO-ZF1格式與WENO-ZF2格式誤差相當(dāng),均小于WENO-Z格式。對于L1誤差,相同粗網(wǎng)格下,WENO-ZF1格式與WENO-ZF2格式誤差相當(dāng),均小于WENO-Z格式;相同細(xì)網(wǎng)格下,三種格式誤差基本相同。綜上,該算例中WENO-ZF1和WENO-ZF2格式可達(dá)到設(shè)計收斂精度,且有較WENO-Z格式更高的分辨率。

表1 t=2時線性對流方程各WENO格式的誤差與精度

3.2 一維歐拉方程

求解歐拉方程時,為了計算的穩(wěn)定性,將方程投影至特征空間,對特征變量進(jìn)行求解,并采用Lax-Friedrichs通量分裂方法對通量進(jìn)行分裂,采用Roe平均方法計算網(wǎng)格單元界面上的變量。

對一維歐拉方程

(27)

進(jìn)行計算,式(27)中守恒變量U=[ρ,ρu,ρE]T,通量F=[ρu,ρu2+p,(ρE+p)u]T。

3.2.1 Sod激波管

計算Sod激波管問題,其初始條件為

(28)

左右邊界設(shè)置為零梯度邊界條件,取200個均勻網(wǎng)格進(jìn)行計算,模擬該問題至t=0.4,計算結(jié)果密度分布及激波與接觸間斷附近的密度分布放大圖如圖4所示。

(a) 密度分布(a) Density distributions

(b) 接觸間斷附近局部放大圖(b) Locally enlarged plot near the contact discontinuity

(c) 激波附近局部放大圖(c) Locally enlarged plot near the shock

圖4表明各格式均可無振蕩計算得到激波與接觸間斷。但比較圖4(b)和圖4(c)中各格式在接觸間斷與激波附近的結(jié)果,WENO-ZF1、WENO-ZF2格式結(jié)果接近,略優(yōu)于WENO-Z格式,而WENO-JS格式抹平較大。這表明,基于模板重組的WENO-ZF1、WENO-ZF2格式耗散低于WENO-Z格式。

3.2.2 激波-密度波相互作用問題

計算激波-密度波相互作用問題,即Shu-Osher問題,該問題的初始條件為

(29)

計算中左右邊界設(shè)置為零梯度邊界條件,取200個均勻網(wǎng)格進(jìn)行計算,模擬該問題至t=1.8。由于該問題無解析解,故將在網(wǎng)格數(shù)為2 000的均勻網(wǎng)格上采用WENO-JS格式計算所得結(jié)果作為基準(zhǔn)解。圖5為網(wǎng)格數(shù)為200時各格式計算結(jié)果的密度分布。

圖5中各格式對鋸齒狀的聲波形成的一系列小激波及x=2.4附近的激波都有較好的捕捉效果,但對高頻振蕩的折射熵波計算效果較差。對比激波附近放大圖,就所得激波陡峭程度而言,WENO-ZF1、WENO-ZF2格式結(jié)果接近,優(yōu)于WENO-Z格式,且優(yōu)于WENO-JS格式。對比折射熵波附近放大圖,在對幅值的計算上,WENO-ZF2結(jié)果接近甚至優(yōu)于WENO-ZF1格式,兩者均優(yōu)于WENO-Z格式,WENO-JS格式則最差。可見,改進(jìn)格式的色散耗散特性得到了優(yōu)化。

(a) 密度分布(a) Density distributions

(b) 圖(a)局部放大圖(b) Locally enlarged plot of graph (a)

3.3 二維歐拉方程

對二維歐拉方程

(30)

做計算,式(30)中守恒變量U=[ρ,ρu,ρv,ρE]T,各方向通量F=[ρu,ρu2+p,ρuv,(ρE+p)u]T,G=[ρv,ρuv,ρv2+p,(ρE+p)v]T。特征投影、Lax-Friedrichs通量分裂等處理方法與求解一維歐拉方程相同。

3.3.1 二維黎曼問題

該問題中取計算區(qū)域為[0,1]×[0,1],初場為

(31)

計算區(qū)域劃分為400×400的均勻網(wǎng)格,各邊界均采用零梯度邊界條件,計算至t=0.8。計算結(jié)果中含一道射流和四道激波等復(fù)雜結(jié)構(gòu),可以有效校驗數(shù)值格式的色散耗散特性和激波捕捉特性。各格式計算結(jié)果密度分布如圖6所示,圖中展示0.14到 1.70均勻分布的40條密度等值線。

由圖6可得,四種WENO格式計算結(jié)果均可正確反映流場分布,實現(xiàn)對激波的有效捕捉。但對比滑移線附近的流場,在對復(fù)雜離散渦的計算上,WENO-ZF1和WENO-ZF2格式計算得到了更多細(xì)節(jié)結(jié)構(gòu),WENO-Z格式次之,WENO-JS更次之。這表明,四種格式中,WENO-ZF1與WENO-ZF2格式耗散最小,WENO-Z格式次之,WENO-JS格式耗散最大。

(a) WENO-JS

(b) WENO-Z

(c) WENO-ZF1

(d) WENO-ZF2

3.3.2 激波-渦相互作用

激波-渦相互作用問題描述二維空間下運動的渦穿過靜止的馬赫數(shù)為1.1的激波時,二者相互作用的問題。計算區(qū)域為[0,2]×[0,1]。初始條件為:在x=0.5處給定一道垂直于x軸的靜止的激波,激波左側(cè)物理量為

(32)

激波右側(cè)物理量由Rankine-Hugoniot條件計算得出。同時,激波左側(cè)疊加一個以(0.25,0.5)為中心的渦,即在激波左側(cè)速度、溫度和熵給定值上施加渦對應(yīng)的脈動量,各脈動量為

(33)

(a) WENO-ZF1

(b) WENO-ZF2

圖8 t=0.6時激波-渦相互作用y=0.5壓力分布Fig.8 Pressure distributions of the shock-vortex interaction problem on central line y=0.5 at t=0.6

圖7與圖8中計算結(jié)果表明WENO-ZF1、WENO-ZF2格式均可較好地模擬這一問題,得到基本的流場分布,得到正確的激波結(jié)構(gòu)和渦核處的低壓區(qū)。比較圖9中激波附近和渦核附近的壓力分布局部放大圖,在對激波的捕捉和渦核壓力的計算上,WENO-ZF1格式與WENO-ZF2格式計算結(jié)果相當(dāng),優(yōu)于WENO-Z格式,WENO-JS格式計算結(jié)果的偏差最大。綜上,改進(jìn)的WENO-ZF1格式、WENO-ZF2格式對激波的捕捉特性和小尺度流場結(jié)構(gòu)的分辨率均優(yōu)于WENO-Z格式與WENO-JS格式。

(a) 激波附近壓力分布(a) Pressure distributions near the shock wave

(b) 渦核附近壓力分布(b) Pressure distributions near the vortex core

表2給出采用不同WENO格式計算上述算例中線性對流方程(初始分布為正弦波,網(wǎng)格數(shù)為320)、二維黎曼問題、激波-渦相互作用問題的相對用時,其他算例由于耗時較短,難以精確統(tǒng)計,未予列出。表2中分別以各算例對應(yīng)WENO-Z格式所用時間作為單位時間對其他格式所用時間進(jìn)行無量綱化,其中WENO-ZF格式表示采用最初的4個子模板加權(quán)公式計算的WENO-ZF格式(該格式與WENO-ZF1格式完全等價,計算結(jié)果完全相同,故在各算例結(jié)果中沒有列出)。由表中數(shù)據(jù)可得,最初的WENO-ZF格式較WENO-Z格式用時長7%～19%,WENO-ZF1格式較WENO-Z格式用時長4%～16%,WENO-ZF2格式較WENO-Z格式用時長1%～4%?？梢?WENO-ZF1格式的等效轉(zhuǎn)化較WENO-ZF格式提高了格式效率,WENO-ZF2格式的β3計算簡化則進(jìn)一步提高了格式效率。

表2 計算不同算例時各WENO格式所用相對時間

4 結(jié)論

本文在5階有限差分WENO-Z格式的基礎(chǔ)上,提出一種模板重組技術(shù):將原格式中的子模板組合得到一中心子模板S3,S3上逼近函數(shù)取對應(yīng)原始子模板的線性組合;在包括S3的4個子模板上采用WENO-Z格式加權(quán)方法得到WENO-ZF格式。通過數(shù)學(xué)運算將WENO-ZF格式轉(zhuǎn)換為原始子模板上的改進(jìn)非線性權(quán)計算方法的WENO-ZF1格式。WENO-ZF和WENO-ZF1格式的構(gòu)建首次實現(xiàn)了格式模板優(yōu)化與非線性權(quán)改進(jìn)二者的等價轉(zhuǎn)換。

為簡化計算,提高格式效率,本文又將WENO-ZF1中S3上光滑因子β3簡化為原始子模板光滑因子的線性組合,得到了計算更為簡潔的WENO-ZF2格式。

ADR分析表明,WENO-ZF1格式和WENO-ZF2格式色散耗散特性得到不同程度的優(yōu)化。數(shù)值實驗表明,WENO-ZF1格式和WENO-ZF2格式可達(dá)到理想的設(shè)計精度;較WENO-Z格式、WENO-JS格式,具有更低的耗散,對激波有更優(yōu)的捕捉特性,對小尺度流場結(jié)構(gòu)有更高的分辨率。格式效率上,WENO-ZF2格式效率高于WENO-ZF1格式,又高于采用4個子模板加權(quán)公式計算的最初WENO-ZF格式。