武從海 李虎 劉旭亮 羅勇 張樹海

(中國空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,四川綿陽 621000)

(中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣研究所,四川綿陽 621000)

引言

可壓縮流動(dòng)中存在激波、接觸間斷等結(jié)構(gòu).為了捕捉這些流動(dòng)結(jié)構(gòu),激波捕捉格式應(yīng)運(yùn)而生.早期的激波捕捉格式以總變差減小(total variation diminishing,TVD)格式[1]、無振蕩、無自由參數(shù)耗散(non-oscillatory and non-free parameter dissipation,NND)格式[2]為主.TVD 方法是當(dāng)前工業(yè)界最常用的激波捕捉方法.經(jīng)典TVD 方法在極值點(diǎn)附近降為一階精度,整體最多只能達(dá)到二階精度,不利于獲得更細(xì)致的流場(chǎng)結(jié)構(gòu)信息.湍流的大渦模擬、直接數(shù)值模擬以及氣動(dòng)聲學(xué)問題的數(shù)值模擬中通常需要捕捉這些細(xì)微結(jié)構(gòu),這要求格式具有良好的分辨率性質(zhì).如果問題中存在間斷,則可采用具有良好分辨率性質(zhì)的高精度激波捕捉格式進(jìn)行模擬.

高精度激波捕捉格式最受關(guān)注的一個(gè)類別是加權(quán)本質(zhì)無振蕩 (weighted essentially non-oscillatory,WENO)格式.WENO 格式最初由Liu 等[3]提出,后由Jiang 等[4]改進(jìn)并達(dá)到了模板上的最優(yōu)精度.后續(xù)研究者們針對(duì)WENO 格式做了大量的改進(jìn)研究,如為了獲得更好流場(chǎng)分辨率的WENO-M[5],WENOZ[6-7],WENO-NS[8],WENO-CU[9]等,為了獲得更好收斂性或穩(wěn)定性的ZSWENO[10],ESWENO[11]等.加權(quán)緊致非線性格式(weighted compact nonlinear scheme,WCNS)[12-13]是另一類借鑒了WENO 思想的激波捕捉方法,張樹海[14]對(duì)兩類格式進(jìn)行了比較.而關(guān)于非結(jié)構(gòu)網(wǎng)格WENO 格式的構(gòu)造可參看文獻(xiàn)[15-16].

WENO 方法采用多個(gè)候選模板,每個(gè)候選模板上均有一個(gè)逼近值,對(duì)這些逼近值進(jìn)行加權(quán)得到最終的逼近值.候選模板上函數(shù)的光滑程度直接影響這些逼近值的權(quán)重.含間斷候選模板的權(quán)重幾乎為0,因此可以達(dá)到捕捉間斷的效果.在連續(xù)區(qū)域中,這些權(quán)值需盡量接近于最優(yōu)權(quán)值,從而減小數(shù)值離散誤差.實(shí)際計(jì)算中,變量在連續(xù)區(qū)域小尺度的波動(dòng)使得WENO 格式可能明顯偏離線性格式,給數(shù)值模擬帶來不利影響.

Wu 等[17]提出了光滑因子單頻波保常設(shè)計(jì)準(zhǔn)則,即光滑因子滿足對(duì)單頻波為常數(shù).這種光滑因子在一般極值點(diǎn)附近的變化幅度小,從而使WENO 格式更接近線性格式.按照這一要求,Wu 等[17-18]構(gòu)造了一類光滑因子,并設(shè)計(jì)了相應(yīng)的WENO-S 格式.由于光滑因子滿足單頻波保常準(zhǔn)則,因此這類格式的近似色散關(guān)系[19]與線性基底格式一致.同時(shí),該類光滑因子在小尺度波動(dòng)區(qū)域變化小,間斷附近變化大,因而WENO-S 格式能更有效地區(qū)分間斷與小尺度波動(dòng),在數(shù)值模擬中能獲得良好的結(jié)果.

WENO-S 格式的光滑因子比經(jīng)典形式的光滑因子[4]更簡(jiǎn)潔,所需浮點(diǎn)運(yùn)算少,因而具有較高的計(jì)算效率.由于其光滑因子僅與子模板上幾個(gè)函數(shù)值相關(guān),與子模板在大模板中的位置無關(guān).對(duì)于線性對(duì)流方程,計(jì)算相鄰數(shù)值通量時(shí)部分光滑因子相同,無需重復(fù)計(jì)算.因此,可以通過消除冗余的光滑因子計(jì)算來提高計(jì)算效率.

本文以7 階WENO-S 格式為例介紹其構(gòu)造方法,包括光滑因子的設(shè)計(jì)準(zhǔn)則,然后對(duì)格式特點(diǎn)進(jìn)行分析,接下來討論計(jì)算效率,針對(duì)線性對(duì)流方程給出基于消除冗余光滑因子計(jì)算的算法,然后討論該方法的使用條件,并將其推廣到其他問題的數(shù)值模擬中.最后通過數(shù)值實(shí)驗(yàn)證明該算法的有效性.

1 7 階WENO-S 格式

1.1 光滑因子單頻波保常設(shè)計(jì)準(zhǔn)則

常見WENO 格式對(duì)于小尺度波動(dòng)通常表現(xiàn)出過多的耗散.其原因可歸結(jié)于極值點(diǎn)附近其權(quán)值偏離線性權(quán)值較大.如最常用的Jiang-Shu 型r點(diǎn)子模板光滑因子[4]

式中,q(x) 為該子模板上的重構(gòu)函數(shù),xj為大模板的中點(diǎn),h為空間步長(zhǎng).在遠(yuǎn)離臨界點(diǎn)的區(qū)域,可由泰勒展開得到臨界點(diǎn)表示一階導(dǎo)數(shù)為0 的點(diǎn),極值點(diǎn)則是一類常見的臨界點(diǎn).在臨界點(diǎn)附近,則有顯然,該類光滑因子在極值點(diǎn)附近下降明顯,這使得相應(yīng)的模板權(quán)重變化較大,對(duì)數(shù)值模擬產(chǎn)生不利影響.

針對(duì)上述情況,需要減小連續(xù)區(qū)域光滑因子的變化幅度.為此,Wu 等[17]提出了光滑因子單頻波保常準(zhǔn)則.依此準(zhǔn)則設(shè)計(jì)的光滑因子對(duì)于單頻波為常數(shù),相應(yīng)WENO 格式非線性權(quán)保持不變,退化為線性權(quán).而對(duì)于其他一般連續(xù)波形,非線性權(quán)的變化也明顯減小.該文認(rèn)為,光滑因子應(yīng)該對(duì)單頻波為常數(shù)的原因如下:

(1)單頻波可視為一個(gè)單頻振蕩的函數(shù),每一部分都具有相同的頻率和最大振幅;

(2)如圖1 所示,幾何意義下圓周處處同樣光滑,而單頻波可作為表示圓周橫坐標(biāo)(或者縱坐標(biāo))關(guān)于角度的函數(shù),可視為函數(shù)意義下處處同樣光滑.

圖1 圓周和正弦函數(shù)的光滑程度示意圖Fig.1 Schematic diagram of the smoothness of a circle and a sine function

Wu 等[17-18]給出了一系列滿足該準(zhǔn)則的光滑因子.以4 點(diǎn)等距模板為例,對(duì)應(yīng)的光滑因子為

1.2 7 階WENO-S 格式

考慮一維雙曲守恒律方程

對(duì)于等距網(wǎng)格xj=jh,其半離散守恒型差分格式為

下文中公式均假設(shè)特征速度 ?f/?u為正.對(duì)于特征速度為負(fù)的情況,可通過對(duì)稱性獲得相應(yīng)的結(jié)果.若無法確定正負(fù),則需要采用通量分裂方法,如Lax-Friedrichs 分裂,時(shí)間格式可以采用高階Runge-Kutta 方法,相關(guān)處理可參見文獻(xiàn)[4]及其參考文獻(xiàn).

2 7 階WENO-S 格式的分析

2.1 收斂精度

7 階WENO-S 格式在連續(xù)區(qū)域通常為7 階精度.但是非線性格式在臨界點(diǎn)附近可能發(fā)生降階,這里臨界點(diǎn)指一階導(dǎo)數(shù)為0 的點(diǎn).臨界點(diǎn)的階數(shù)用ncp表示,如ncp=1 表示f′=0,f′′≠0,ncp=2表示f′=0,f′′=0,f′′′≠0.根據(jù)Wu 等[17]的分析,WENO7-S 在二階及以下臨界點(diǎn)可保持7 階精度,而在3 階及以上臨界點(diǎn)可能會(huì)發(fā)生降階.

對(duì)于不小于3 階的臨界點(diǎn),如果要保持7 階收斂精度,可在7 階WENO-S 格式非線性權(quán)的計(jì)算中采用與空間步長(zhǎng)相關(guān)的 ε值[20].

2.2 極值點(diǎn)附近的模擬

極值點(diǎn)附近波形的數(shù)值模擬是激波捕捉方法的難點(diǎn).正是由于極值點(diǎn)的存在,經(jīng)典TVD 格式只能達(dá)到二階精度第1.1 節(jié)中光滑因子單頻波保常設(shè)計(jì)準(zhǔn)則的出發(fā)點(diǎn)之一也是降低極值點(diǎn)附近光滑因子的變化幅度[17].

極值點(diǎn)一般位于連續(xù)區(qū)域,為獲得更好的模擬結(jié)果,通常需要WENO 格式的非線性權(quán)系數(shù)盡量接近于線性權(quán)系數(shù).絕大多數(shù)極值點(diǎn)屬于一階臨界點(diǎn).由Taylor 展開可知,經(jīng)典光滑因子 β(JS)[4]在一階臨界點(diǎn)附近和一般連續(xù)區(qū)域分別為本文中簡(jiǎn)略起見,光滑因子在不涉及到具體點(diǎn)時(shí)省去下標(biāo).而7 階WENO-S 的光滑因子 β(S)則保持為即光滑因子變化較小.那么 β(S)不僅對(duì)單頻波為常數(shù),而且在一般波形的極值點(diǎn)附近變化也較小.因此,7 階WENO-S 在極值點(diǎn)附近更接近線性格式,從而數(shù)值模擬誤差通常更小.

2.3 間斷穩(wěn)定性

為了消除或減小捕捉間斷時(shí)的振蕩現(xiàn)象,要求WENO格式中含間斷子模板的權(quán)系數(shù)盡量接近于0.光滑因子 β(S)在間斷區(qū)域?yàn)镺(1),連續(xù)區(qū)域?yàn)槎?jīng)典光滑因子 β(JS)在兩類區(qū)域分別為O(1)和那么光滑因子 β(S)在兩類區(qū)域之間的差距更大.若權(quán)系數(shù)公式其他部分相同,則7 階WENO-S格式中含間斷子模板的權(quán)系數(shù)更小.這使得該格式具有良好的間斷穩(wěn)定性.

2.4 分辨率性質(zhì)

Lele[21]闡述了分辨率性質(zhì)對(duì)差分格式的意義,并采用Fourier 方法分析了緊致差分格式的分辨率特性.但是,Fourier 方法不適用于非線性格式.針對(duì)此問題,文獻(xiàn)中最常用的方法是近似色散關(guān)系(approximate dispersion relation,ADR)分析[18].該方法首先使用數(shù)值方法對(duì)不同無量綱波數(shù) κ的單頻波進(jìn)行短時(shí)間的推進(jìn),然后采用離散傅里葉變換結(jié)合初值計(jì)算得到有效波數(shù).有效波數(shù)的實(shí)部代表數(shù)值解的相速度,通常也用來分析色散性能,而虛部代表耗散率.Li 等[22]和毛枚良等[23]中將ADR 方法進(jìn)行了簡(jiǎn)化,去掉了時(shí)間推進(jìn)過程.

由于WENO-S 格式采用了對(duì)單頻波為常數(shù)的光滑因子,該類格式在計(jì)算單頻波時(shí)退化為線性格式,因而其近似色散關(guān)系與線性基底格式一致,優(yōu)于絕大多數(shù)同階WENO 格式.文獻(xiàn)[17]中圖2 給出了幾種7 點(diǎn)WENO 格式及其基底格式的近似色散關(guān)系,其中7 階WENO-S 與線性基底格式完全重合.

Zhao 等[24]和Li 等[25-26]中在ADR 方法的基礎(chǔ)上考慮了非線性響應(yīng),即單頻波在非線性格式下產(chǎn)生的其他波數(shù)的雜波.Cravero 等[27]中借用熱力學(xué)中的“溫度”來表示非線性格式的這種性質(zhì).由于WENO-S 格式對(duì)于單頻波退化為線性格式而不產(chǎn)生雜波,因此在這類分析方法中,該類格式不存在非線性響應(yīng),具有溫度0 的特征.

3 7 階WENO-S 格式的計(jì)算效率

7 階WENO-S 格式具有良好的計(jì)算效率.其光滑因子 β(S)比經(jīng)典光滑因子 β(JS)更簡(jiǎn)潔,所需計(jì)算量更小.但是該格式的計(jì)算效率還需考慮參數(shù) τ的計(jì)算.Wu 等[17]對(duì)比了幾種WENO 格式單個(gè)數(shù)值通量的計(jì)算量,并通過數(shù)值實(shí)驗(yàn)說明了其效率優(yōu)勢(shì).但是,對(duì)于一些特殊問題,還可進(jìn)一步提升該格式的計(jì)算效率.

3.1 線性對(duì)流方程的效率提升方法

采用經(jīng)典的WENO 格式(3 階WENO 除外)時(shí),各子模板光滑因子的計(jì)算公式形式不同,在相鄰數(shù)值通量的計(jì)算中不會(huì)出現(xiàn)相同的光滑因子.而WENO-S 的光滑因子在各子模板上的計(jì)算公式除下標(biāo)不同外形式一致.那么對(duì)于線性對(duì)流方程,計(jì)算相鄰數(shù)值通量時(shí),部分光滑因子相同.因而可以通過消除WENO-S 的冗余光滑因子計(jì)算來提高計(jì)算效率.

具體計(jì)算中,可將所有光滑因子先行計(jì)算并存儲(chǔ).另外,考慮到大模板上的光滑因子 τ(S)的計(jì)算,同時(shí)需要存儲(chǔ)的還有=?fj+3fj+1?3fj+2+fj+3組成的序列.那么,當(dāng)網(wǎng)格點(diǎn)數(shù)較多時(shí),對(duì)于7 階WENO-S格式,子模板光滑因子的計(jì)算量約為原來的1/4.

考慮一維線性對(duì)流方程

假設(shè)對(duì)流速度a>0.設(shè)計(jì)算網(wǎng)格點(diǎn)為 {x1,x2,···,xN}.采用虛擬網(wǎng)格點(diǎn)的方法處理邊界條件,對(duì)于7 階WENO 格式,左端需設(shè)置4 個(gè)虛擬點(diǎn),右端需設(shè)置3 個(gè)虛擬點(diǎn).計(jì)算點(diǎn)和虛擬點(diǎn)上的通量組成的序列為 {f?3,f?2,f?1,···,fN+3}.對(duì)于a<0 的情況,可通過對(duì)稱性得到.

簡(jiǎn)便起見,這里用WENO7-JS 表示經(jīng)典的7 階WENO 格式,而WENO7-Z 在WENO7-JS 的基礎(chǔ)上采用了Z 型權(quán)系數(shù)[28],其中冪因子取1,WENO7-S表示7 階WENO-S 格式.

表1 幾種WENO 格式計(jì)算N+1 個(gè)數(shù)值通量所需浮點(diǎn)運(yùn)算量Table 1 The number of floating point operations required to calculate N+1 numerical fluxes for several WENO schemes

表1 中可以看到,采用該快速算法可以有效降低WENO-S 格式的計(jì)算量,N很大時(shí)數(shù)值通量計(jì)算量降低 5/13≈ 38.46%.實(shí)際程序運(yùn)行時(shí),由于還有其他部分的運(yùn)算,其整體計(jì)算效率提高的程度應(yīng)小于該值.在后文的數(shù)值算例測(cè)試中,其計(jì)算時(shí)間約節(jié)省20%.

3.2 其他適用情況

上述效率提升方法也適用于其他一些問題.其關(guān)鍵在于一條網(wǎng)格線上用于WENO 重構(gòu)或插值的量可以表示為一個(gè)序列,進(jìn)行網(wǎng)格線上不同位置的WENO 重構(gòu)或插值時(shí),選取序列中相應(yīng)連續(xù)幾個(gè)位置的值,然后進(jìn)行WENO 重構(gòu)或插值.

除7 階W E N O-S 格式外,還有一些其他WENO 格式也可采用這種方法提升計(jì)算效率.其要求是光滑因子在各子模板上的計(jì)算公式除下標(biāo)不同外形式一致.如更高階的WENO-S 格式[18]、WENO-XS格式[31]、WENO-LOC 格式[32]和FWENO 格式[30]等.

在稍復(fù)雜流動(dòng)問題的數(shù)值模擬中,一些額外的處理過程可能會(huì)破壞這一條件.

(1) 通量分裂

采用守恒型差分方法求解非線性方程時(shí),考慮到迎風(fēng)性,通常采用通量分裂的方法.如果要采用這種效率提升方法,則必須采用全局通量分裂.而局部特征分裂將破壞上述條件,從而無法采用該方法提升效率.

對(duì)于有限體積法和WCNS 型[12]的差分方法,可采用Godunov 類方法計(jì)算通量.該類方法中進(jìn)行非線性插值或重構(gòu)的量不是通量,而是原始變量或守恒變量,避免了上述全局通量分裂的要求,可采用該加速方法計(jì)算.

(2) 特征投影

對(duì)于雙曲型偏微分方程組,由于信號(hào)沿特征方向傳播的特點(diǎn),數(shù)值求解采用特征投影處理更符合方程的特點(diǎn),在含間斷問題中結(jié)果更準(zhǔn)確.對(duì)于非線性問題,投影時(shí)一般采用局部特征矩陣.

不采用特征投影而直接對(duì)各個(gè)方程進(jìn)行數(shù)值離散求解的方法稱為組元型方法.組元型方法一般僅用于較弱激波的捕捉,對(duì)于強(qiáng)激波,該處理可能會(huì)產(chǎn)生較強(qiáng)的波后非物理振蕩,對(duì)計(jì)算結(jié)果產(chǎn)生明顯不利影響.但是局部特征投影破壞了該效率提升方法的前提.那么,這種效率提升方法只能結(jié)合組元型求解方法使用.

某些特定問題中部分方程具有明確的特征方向,即使剩下的方程組成的方程組采用特征投影求解,這些方程也可以采用該方法減少計(jì)算量.如多組分問題中的組分方程[33]、多介質(zhì)流問題中表示界面的方程[34]等.

4 數(shù)值算例

本節(jié)先進(jìn)行極值點(diǎn)附近權(quán)系數(shù)偏離情況測(cè)試,然后采用多個(gè)算例對(duì)7 階WENO-S 的性能進(jìn)行測(cè)試驗(yàn)證,并且同時(shí)測(cè)試所提效率提升方法的有效性.所選算例包括一維對(duì)流問題、球面波傳播問題、二維旋轉(zhuǎn)問題、二維小擾動(dòng)傳播問題及一維和二維無黏流動(dòng)問題.空間離散采用WENO7-JS、WENO7-Z和WENO7-S 格式,時(shí)間推進(jìn)采用3 階TVD Runge-Kutta 方法[4].在進(jìn)行時(shí)間效率對(duì)比時(shí),WENO7-JS和WENO-7S 在格式的后面加0 或1 以分別代表經(jīng)典計(jì)算方法和快速計(jì)算方法.

4.1 極值點(diǎn)附近權(quán)系數(shù)偏離程度測(cè)試

WENO 格式的權(quán)系數(shù)在極值點(diǎn)附近常常明顯偏離線性權(quán)系數(shù),這可能導(dǎo)致數(shù)值模擬誤差的增大.本算例比較極值點(diǎn)附近幾種WENO 格式的權(quán)系數(shù)偏離程度.下面用表示某點(diǎn)通量的權(quán)系數(shù)偏離值,統(tǒng)計(jì)極值點(diǎn)附近若干通量點(diǎn)的最大偏離值以及平均偏離值.所選取的3 個(gè)函數(shù)為:(1)f1(x)=,(2)f2(x)=ex?x?1,(3)f3(x)=sin4(πx/4).極值點(diǎn)均取x=0.該點(diǎn)為f1(x)和f2(x)的一階臨界點(diǎn),f3(x)的3 階臨界點(diǎn).

計(jì)算中取該極值點(diǎn)作為網(wǎng)格點(diǎn)之一,統(tǒng)計(jì)該點(diǎn)左右各4 個(gè)通量點(diǎn)的情況.網(wǎng)格間距均取h=0.1,其結(jié)果放在表2 中.可以看到表中WENO7-S 的權(quán)系數(shù)偏離值在極值點(diǎn)相對(duì)更小.其中一階臨界點(diǎn)附近小一個(gè)量級(jí)甚至更多,這與第2.2 小節(jié)的結(jié)論相符.

表2 權(quán)系數(shù)偏離值的平均值與最大值Table 2 The average values and maximum values of the deviations of weights

對(duì)于f3(x)=sin4(πx/4) 的3 階臨界點(diǎn)x=0,幾種格式的權(quán)系數(shù)偏離值均較大,這表明WENO 格式在這種臨界點(diǎn)附近與線性基底格式有明顯差異,最終導(dǎo)致它們可能在此問題中發(fā)生降階.注意當(dāng)網(wǎng)格間距進(jìn)一步減小時(shí),WENO7-JS 由于更大的ε值,其權(quán)系數(shù)偏離值可能會(huì)小于其他兩個(gè)格式,使得其更容易在密網(wǎng)格下恢復(fù)7 階精度.

4.2 一維對(duì)流問題

考慮一組包括高斯波、方波、尖三角波和橢圓波在內(nèi)的復(fù)雜組合包以速度1 向右傳播的問題[4],其控制方程為

其中的常數(shù)a=0.5,z=?0.7,δ=0.005,α=10,β=.兩端采用周期邊界條件.

這里考慮波形的遠(yuǎn)距傳輸,計(jì)算終止時(shí)間取T=2000,網(wǎng)格點(diǎn)數(shù)取400.為減小時(shí)間離散誤差的影響,CFL 數(shù)取0.01.計(jì)算結(jié)果顯示在圖2 中.對(duì)于這種遠(yuǎn)距傳輸問題,WENO7-JS 耗散過大,而WENO7-Z 則在間斷附近出現(xiàn)了明顯的數(shù)值振蕩,但是在三角波峰值附近表現(xiàn)最好.WENO7-S 格式對(duì)于其他幾種波形均獲得了優(yōu)良的結(jié)果.WENO7-JS,WENO7-Z 和WENO7-S 計(jì)算結(jié)果的L1誤差依次為2.51 × 10?1,4.67 × 10?2和4.79 × 10?2.此問題中WENO7-Z 的誤差最小,這可能是因?yàn)槠淙遣ǜ浇膬?yōu)良表現(xiàn)和更小的間斷抹平程度.

圖2 組合波的遠(yuǎn)距傳輸問題,網(wǎng)格點(diǎn)數(shù)N=400,計(jì)算終止時(shí)間T=2000Fig.2 Long-distance advection of combined waves with the number of grid points N=400 and the end time T=2000

表3 給出了幾種格式的計(jì)算時(shí)間,結(jié)果表明WENO7-S0 比WENO7-JS1 略快,后續(xù)算例也大多如此,而表1 中的浮點(diǎn)計(jì)算統(tǒng)計(jì)則剛好相反.WENO7-JS0 和WENO7-Z 之間也有類似表現(xiàn).在實(shí)際計(jì)算中,運(yùn)行時(shí)間與計(jì)算平臺(tái)、浮點(diǎn)計(jì)算單元使用率和緩存命中率等均有較大關(guān)系[35].本文計(jì)算均采用Fortran 程序在CPU 為AMD?Ryzen?R5 4500 U 的計(jì)算機(jī)上運(yùn)行.WENO7-S0 的光滑因子具有相同的形式,編譯器可能將其處理為向量計(jì)算,這會(huì)提升其計(jì)算速度.WENO7-Z 相比WENO7-JS0 在光滑因子部分的計(jì)算時(shí)間減少了,而τ及權(quán)系數(shù)的計(jì)算時(shí)間增大了.可能在程序運(yùn)行中,增大的這一部分時(shí)間大于減小的時(shí)間.

表3 組合波遠(yuǎn)距傳輸問題的計(jì)算時(shí)間Table 3 Computational time of long-distance advection of combined waves

當(dāng)消除WENO7-S 的冗余光滑因子計(jì)算后,其效率得到進(jìn)一步提高,該問題中WENO7-S1 比WENO7-S0 提高24.64%.

4.3 球面波傳播問題

球面波傳播問題是一個(gè)計(jì)算氣動(dòng)聲學(xué)的標(biāo)準(zhǔn)測(cè)試算例[36],其控制方程為

計(jì)算區(qū)域?yàn)閇5,450],初始條件為u=0,在r=5處的邊界條件為u=sin(ωt).

取 ω=π/4,計(jì)算終止時(shí)間為T=400,空間步長(zhǎng)為 dr=1.為減小時(shí)間推進(jìn)誤差的影響,時(shí)間步長(zhǎng)取dt=0.1.計(jì)算結(jié)果顯示在圖3 中.圖中幾種格式計(jì)算結(jié)果的波幅有明顯差距,WENO7-JS 耗散最大,而WENO7-S 由于其良好的分辨率性質(zhì),獲得了最優(yōu)的計(jì)算結(jié)果.

圖3 球面波傳播問題,計(jì)算終止時(shí)間為T=400,空間步長(zhǎng)為dr=1.時(shí)間步長(zhǎng)為dt=0.1Fig.3 The spherical wave propagation problem with the end time T=400,the grid length dr=1,and the time step dt=0.1

由于該問題計(jì)算時(shí)間過短,計(jì)算時(shí)間差距過小.因此,表3 給出了時(shí)間步長(zhǎng) dt取0.01 時(shí)的計(jì)算時(shí)間統(tǒng)計(jì).結(jié)果表明,WENO7-S 計(jì)算效率較高,而消除冗余計(jì)算后,效率進(jìn)一步提升21.82%.

表4 球面波傳播問題的計(jì)算時(shí)間,dt=0.01Table 4 Computational time of spherical wave propagation problem with dt=0.01

4.4 二維旋轉(zhuǎn)問題

這個(gè)問題來自于Zalesak[37]的文章.該問題中,一個(gè)割開的圓柱繞著附近的軸旋轉(zhuǎn).該運(yùn)動(dòng)的控制方程為

這里u=??(y?y0),v=?(x?x0),其中 ?(=2π/360)是旋轉(zhuǎn)角速度.這里計(jì)算區(qū)域?yàn)?(x,y)∈[0,10]×[0,10],(x0,y0)=(5,5)是旋轉(zhuǎn)軸的坐標(biāo).初始時(shí)刻,割開圓柱的中心坐標(biāo)為 (xc,yc)=(5,7.5),圓柱半徑為1.5,兩條割線的橫坐標(biāo)分別為4.75 和5.25,割開截止位置的縱坐標(biāo)為8.5;割開圓柱上的 ?值為3.0,其他區(qū)域則是1.0.根據(jù)這個(gè)控制方程,旋轉(zhuǎn)過程中割開圓柱應(yīng)該保持原始的形狀.本問題將該割開圓柱旋轉(zhuǎn)5 個(gè)周期.

該問題中兩個(gè)分量速度既可能為正又可能為負(fù),通常同時(shí)包含正負(fù)速度時(shí)需要進(jìn)行通量分裂.但是該問題每條網(wǎng)格線上特征速度可確定為正或負(fù),無需通量分裂處理.

采用200 × 200 均勻網(wǎng)格進(jìn)行計(jì)算,網(wǎng)格點(diǎn)坐標(biāo)為xi=(i?0.5)·?x,yj=(j?0.5)·?y,其中Δx和Δy分別為兩個(gè)方向的網(wǎng)格間距.圖4 給出了計(jì)算結(jié)果,其中高度代表 ?值.圖中可以看到,WENO7-JS 對(duì)圓柱邊緣處的抹平最為嚴(yán)重,頂部區(qū)域有塌陷現(xiàn)象.WENO7-Z 的邊緣更銳利,但是附近有輕微的過沖現(xiàn)象.WENO7-S 則避免了過度的抹平和邊緣處過沖現(xiàn)象.

為了更清晰地顯示計(jì)算結(jié)果的差異,截取直線y=y150=7.475和x=x100=4.975上的結(jié)果進(jìn)行對(duì)比.圖5(a)為兩條直線的位置示意圖,圖5(b)和圖5(c)給出了兩條線上的結(jié)果.從圖5(b) 中可以看到WENO7-JS在割開部分表現(xiàn)出明顯的耗散,WENO7-Z 和WENO7-S 則出現(xiàn) 了下沖,其 中WENO7-Z 的下沖程度相對(duì)較大,這也與圖5(c)中區(qū)間[6,8.5]的結(jié)果相符.圖5(b)中WENO7-Z 在x=6.5附近出現(xiàn)了微弱的上沖.圖5(c)中區(qū)間[8.5,9]處WENO7-JS 的峰值明顯小于精確值3,這也對(duì)應(yīng)圖4(b)中的頂部塌陷現(xiàn)象.

圖4 二維旋轉(zhuǎn)問題計(jì)算結(jié)果與精確解Fig.4 The computational results and the exact solution of the twodimensional rotation problem

圖5 二維旋轉(zhuǎn)問題兩條線上的結(jié)果對(duì)比圖Fig.5 Comparison of results on two lines of the two-dimensional rotation problem

采用

來表示上沖/下沖幅度.WENO7-JS,WENO7-Z和WENO7-S 的 δ值分別為1.158 × 10?4,0.156和5.944 × 10?2.這說明WENO7-Z 更容易出現(xiàn)上沖/下沖現(xiàn)象.而3 種格式的L1誤差依次為2.298,1.774 和1.959.說明盡管WENO7-Z 有較明顯的上沖/下沖,但是其誤差相對(duì)最小.從圖5(b)和圖5(c)中也可看出,WENO7-Z 的間斷抹平程度最小,而間斷附近貢獻(xiàn)了大部分誤差量,這使得其誤差小于WENO7-S.

為了考察不同網(wǎng)格下的計(jì)算情況,采用200 ×200,400 × 400 和800 × 800 這3 套均勻網(wǎng)格對(duì)該問題進(jìn)行模擬,其時(shí)間步長(zhǎng)依次設(shè)為0.1,0.05 和0.025.為減少模擬時(shí)間,這里僅將圓柱旋轉(zhuǎn)一個(gè)周期.計(jì)算結(jié)果的L1誤差和上沖/下沖幅度由表5 給出.顯然WENO7-Z 的L1誤差依然最小.而由于間斷的存在,3 種格式的L1誤差收斂精度均約為1 階.對(duì)于上沖/下沖幅度 δ,200 × 200 和400 × 400 網(wǎng)格時(shí)WENO7-JS 最小,而800 × 800 時(shí)WENO7-S 最小.隨著網(wǎng)格的加密,WENO7-JS 和WENO7-Z 沒有表現(xiàn)出明顯的變化趨勢(shì),而WENO7-S 的 δ值迅速減小,體現(xiàn)了其良好的間斷穩(wěn)定性.

表5 幾種格式不同網(wǎng)格下計(jì)算結(jié)果的L1 誤差和上沖/下沖幅度Table 5 The L1 error and up/down overshooting amplitude of the computational results with different schemes and grids

表6 給出了200 × 200 網(wǎng)格旋轉(zhuǎn)5 個(gè)周期的計(jì)算時(shí)間對(duì)比.結(jié)果表明WENO7-S 計(jì)算效率較高,新的計(jì)算方法進(jìn)一步提升了其計(jì)算效率,該問題中提升23.48%.

表6 二維旋轉(zhuǎn)問題的計(jì)算時(shí)間Table 6 Computational time for the two-dimensional rotation problem

4.5 二維小擾動(dòng)傳播問題

控制方程采用二維線化Euler 方程

Mx=0.5,My=0是x和y方向的平均流馬赫數(shù).U′表示密度、速度及壓強(qiáng)的擾動(dòng)量.計(jì)算區(qū)域?yàn)?x,y)∈[?100,100]×[?100,100].初值為

該問題為計(jì)算氣動(dòng)聲學(xué)的一個(gè)標(biāo)準(zhǔn)算例[36,38].

根據(jù)迎風(fēng)格式的要求,采用了每條網(wǎng)格線上的全局通量分裂.x方向計(jì)算時(shí)正負(fù)通量分別為方向計(jì)算時(shí)正負(fù)通量分別為

計(jì)算采用的網(wǎng)格為200 × 200,終止時(shí)間為T=30.圖6 給出了水平中心線y=0上的密度分布圖.可以看到幾種WENO 格式均取得了與精確解非常接近的結(jié)果.相比而言,WENO7-S 在峰值附近取得了最佳的結(jié)果.表7 給出了幾種格式的計(jì)算時(shí)間.本算例中WENO7-S 計(jì)算效率最高,且還可受益于本文計(jì)算效率提升方法,該問題中提升17.84%.

表7 二維小擾動(dòng)傳播問題的計(jì)算時(shí)間Table 7 Computational time of two-dimensional small disturbance propagation problem

圖6 二維小擾動(dòng)傳播問題計(jì)算結(jié)果,網(wǎng)格200 × 200,終止時(shí)間T=30Fig.6 Computational results of two-dimensional small disturbance propagation problem,with grid 200 × 200 and end time T=30

4.6 一維無黏流動(dòng)問題

控制方程為一維Euler 方程

式中 ρ,u,p,E表示密度、速度、壓強(qiáng)和內(nèi)能.氣體狀態(tài)方程為

其中比熱比 γ=1.4.

這里計(jì)算Shu-Osher 問題,該問題中馬赫數(shù)3 的激波與正弦波相互干擾,誘發(fā)高頻振動(dòng),其初值條件為

計(jì)算終止時(shí)間T=1.8,這里取網(wǎng)格點(diǎn)為200.

分別采用組元型和特征型方法結(jié)合幾種幾種WENO 格式進(jìn)行計(jì)算.采用5 階WENO-JS 在4000 個(gè)均勻網(wǎng)格點(diǎn)的計(jì)算結(jié)果作為參考解,圖7(a)為參考解的密度分布圖.圖7(b)對(duì)比了幾種WENO 格式的計(jì)算結(jié)果.圖例中幾種WENO 格式省略了名稱前綴WENO7,并用后綴Comp 表示組元型方法,Char 表示特征型方法.

圖7 Shu-Osher 問題計(jì)算結(jié)果的密度對(duì)比Fig.7 The density distributions of computational results of Shu-Osher problem

從圖7(b)中可以看出,在高頻震蕩區(qū)域,特征型的結(jié)果更好,這在WENO7-JS 格式的結(jié)果中體現(xiàn)得更為明顯.在圖中高頻震蕩左側(cè),WENO7-Z 和WENO7-S 的組元型結(jié)果相對(duì)更接近參考解.幾種WENO 格式中,WENO7-S 的結(jié)果最好,其中特征型的結(jié)果比組元型的結(jié)果略好.

該問題計(jì)算時(shí)間較短,為此將網(wǎng)格數(shù)設(shè)為2000,此時(shí)特征型WENO7-JS 的計(jì)算時(shí)間為5.317 s.對(duì)于組元型計(jì)算方法,表8 給出了幾種方法的計(jì)算時(shí)間.WENO7-JS 特征型的計(jì)算時(shí)間約為組元型的1.6 倍.組元型方法中WENO7-S 格式計(jì)算時(shí)間最短,其中本文加速方法提升了其25.26%的計(jì)算速度.

表8 Shu-Osher 問題計(jì)算時(shí)間,N=2000Table 8 Computational time of Shu-Osher problem,N=2000

該問題中激波馬赫數(shù)為3,并非弱激波.此時(shí)組元型和特征型的結(jié)果之間的差距可能相對(duì)較明顯,如WENO7-JS 格式.但是用WENO7-S 格式計(jì)算時(shí)兩種方法的差距較小,可考慮采用組元型方法并結(jié)合本文加速方法.

4.7 二維無黏流動(dòng)問題

控制方程為守恒形式的二維Euler 方程

式 中 ρ,u,v,p,E表示密 度、x向速度、y向速度、壓強(qiáng)和內(nèi)能.氣體狀態(tài)方程為

其中比熱比 γ=1.4.

這里采用二維激波旋渦相互作用問題作為測(cè)試算例[4],該問題描述了一個(gè)靜止的激波和旋渦的相互作用.計(jì)算區(qū)域?yàn)?[ 0,2]×[0,1].一個(gè)馬赫數(shù)為1.1 的靜止激波放置在x=0.5 處且與x軸垂直.左狀態(tài)為 (ρ,u,v,p)=右狀態(tài)由Rankine-Hugoniot 關(guān)系給出.一個(gè)中心在 (xc,yc)=(0.25,0.5)的小旋渦放在激波左側(cè).這個(gè)旋渦可被視為速度 (u,v)、熵溫度T=p/ρ的小擾動(dòng),即

由于特征投影方法采用了流場(chǎng)當(dāng)?shù)氐耐队熬仃?WENO7-S1 的效率提升方法無法使用.這里采用組元型方法,該方法計(jì)算量比特征投影方法小,但是不宜用于強(qiáng)激波問題.該問題來流馬赫數(shù)為1.1,說明激波較弱,計(jì)算結(jié)果對(duì)比也表明兩種方法差距很小.另外,特征型WENO 格式中特征投影及逆投影所需計(jì)算量很大.本算例采用WENO7-JS1 結(jié)合特征型方法時(shí)計(jì)算時(shí)間為41.893 s,約為該格式結(jié)合組元型方法計(jì)算時(shí)間13.195 s 的3 倍.這高于一維問題中的1.6 倍.原因是特征投影矩陣在二維問題中是四維,其乘法計(jì)算量顯著高于一維問題中的三維矩陣.可以預(yù)計(jì),在三維問題中,特征型算法的計(jì)算時(shí)間相比組元型多3 倍以上.

采用組元型方法計(jì)算時(shí),x方向正負(fù)通量分別為其 中λx為 該x向網(wǎng)格 線上所有矩陣 ?F/?U特征值絕對(duì)值的最大值.y方向計(jì)算時(shí)正負(fù)通量分別為其中λy為該y向網(wǎng)格線上所有矩陣 ?G/?U特征值絕對(duì)值的最大值.

圖8 中給出了計(jì)算結(jié)果的壓強(qiáng)云圖,其等值線范圍為1.19～ 1.37,共90 條.其中圖8(a)為特征型WENO7-JS 的結(jié)果,可以看到在較弱的激波下,組元型和特征型結(jié)果相差較小.組元型方法的計(jì)算結(jié)果中,WENO7-JS 和WENO7-S 較為相近,而WENO7-Z在激波附近聚集了更多等值線.這可能是由于WENO7-Z 在間斷附近耗散更小,激波波后非物理振蕩更為明顯,而這種振蕩會(huì)導(dǎo)致相應(yīng)值的等值線多次出現(xiàn).

圖8 激波旋渦相互作用問題,壓強(qiáng)等值線范圍1.19-1.37,共90 條,網(wǎng)格200 × 100,終止時(shí)間T=0.6Fig.8 Shock vortex interaction problem.90 pressure isolines ranging from 1.19 to 1.37.Component-wise seventh-order WENO schemes with grid 200 × 100 and end time T=0.6

表9 給出了幾種格式的計(jì)算時(shí)間.結(jié)果表明WENO7-S 格式效率相對(duì)較高,WENO7-S1 的效率提升方法進(jìn)一步減少了21.66%的計(jì)算時(shí)間.

表9 激波旋渦相互作用問題的計(jì)算時(shí)間Table 9 Computational time of shock vortex interaction problem

5 結(jié)論

WENO-S 格式具有許多優(yōu)良的性質(zhì).特別是其光滑因子對(duì)單頻波為常數(shù),使得其是一類滿足近似色散關(guān)系與線性基底格式一致的激波捕捉格式.另外,該格式能有效區(qū)分間斷與小尺度波動(dòng),在連續(xù)區(qū)域誤差小,間斷附近也能保證良好的穩(wěn)定性,在數(shù)值模擬中能獲得良好的結(jié)果.

由于WENO-S 格式的光滑因子在各子模板上的計(jì)算公式除下標(biāo)不同外形式一致,在計(jì)算線性對(duì)流方程的相鄰數(shù)值通量時(shí),部分光滑因子相同.本文基于7 階WENO-S 格式,介紹了如何通過避免冗余光滑因子計(jì)算提高格式計(jì)算效率.除WENO-S 格式外,還有一些其他WENO 格式也可采用這種方法提升計(jì)算效率.

這種方法可推廣至多種問題,其應(yīng)用條件是一條網(wǎng)格線上用于WENO 重構(gòu)或插值的量可以表示為一個(gè)序列.因此在非線性方程中使用通量分裂時(shí)只能采用全局通量分裂.非線性方程的另一類處理方法是對(duì)原始變量或守恒變量直接進(jìn)行WENO 重構(gòu)或插值,然后采用Godunov 類方法求得數(shù)值通量.這時(shí)也可以采用避免冗余光滑因子計(jì)算的算法.

對(duì)于可壓縮流動(dòng)問題,局部特征投影求解會(huì)破壞這種消除冗余計(jì)算的方法.因此,在可壓縮流動(dòng)模擬中,該方法只能結(jié)合組元型方法求解.某些流動(dòng)問題中部分方程的求解不受該限制,仍可采用該方法減小計(jì)算量,如多組分問題中的組分方程、多介質(zhì)流問題中表示界面的方程等.

數(shù)值實(shí)驗(yàn)結(jié)果表明7 階WENO-S 格式同時(shí)具有良好的小尺度結(jié)構(gòu)分辨率和激波捕捉能力.在計(jì)算效率方面,本文所提方法能有效減少7 階WENOS 格式約20%的計(jì)算時(shí)間.