7 階WENO-S 格式的計算效率研究1)

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

(中國空氣動力研究與發展中心空氣動力學國家重點實驗室,四川綿陽 621000)

(中國空氣動力研究與發展中心計算空氣研究所,四川綿陽 621000)

引言

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

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

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

Wu 等[17]提出了光滑因子單頻波保常設計準則,即光滑因子滿足對單頻波為常數.這種光滑因子在一般極值點附近的變化幅度小,從而使WENO 格式更接近線性格式.按照這一要求,Wu 等[17-18]構造了一類光滑因子,并設計了相應的WENO-S 格式.由于光滑因子滿足單頻波保常準則,因此這類格式的近似色散關系[19]與線性基底格式一致.同時,該類光滑因子在小尺度波動區域變化小,間斷附近變化大,因而WENO-S 格式能更有效地區分間斷與小尺度波動,在數值模擬中能獲得良好的結果.

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

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

1 7 階WENO-S 格式

1.1 光滑因子單頻波保常設計準則

常見WENO 格式對于小尺度波動通常表現出過多的耗散.其原因可歸結于極值點附近其權值偏離線性權值較大.如最常用的Jiang-Shu 型r點子模板光滑因子[4]

式中,q(x) 為該子模板上的重構函數,xj為大模板的中點,h為空間步長.在遠離臨界點的區域,可由泰勒展開得到臨界點表示一階導數為0 的點,極值點則是一類常見的臨界點.在臨界點附近,則有顯然,該類光滑因子在極值點附近下降明顯,這使得相應的模板權重變化較大,對數值模擬產生不利影響.

針對上述情況,需要減小連續區域光滑因子的變化幅度.為此,Wu 等[17]提出了光滑因子單頻波保常準則.依此準則設計的光滑因子對于單頻波為常數,相應WENO 格式非線性權保持不變,退化為線性權.而對于其他一般連續波形,非線性權的變化也明顯減小.該文認為,光滑因子應該對單頻波為常數的原因如下:

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

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

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

Wu 等[17-18]給出了一系列滿足該準則的光滑因子.以4 點等距模板為例,對應的光滑因子為

1.2 7 階WENO-S 格式

考慮一維雙曲守恒律方程

對于等距網格xj=jh,其半離散守恒型差分格式為

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

2 7 階WENO-S 格式的分析

2.1 收斂精度

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

對于不小于3 階的臨界點,如果要保持7 階收斂精度,可在7 階WENO-S 格式非線性權的計算中采用與空間步長相關的 ε值[20].

2.2 極值點附近的模擬

極值點附近波形的數值模擬是激波捕捉方法的難點.正是由于極值點的存在,經典TVD 格式只能達到二階精度第1.1 節中光滑因子單頻波保常設計準則的出發點之一也是降低極值點附近光滑因子的變化幅度[17].

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

2.3 間斷穩定性

為了消除或減小捕捉間斷時的振蕩現象,要求WENO格式中含間斷子模板的權系數盡量接近于0.光滑因子 β(S)在間斷區域為O(1),連續區域為而經典光滑因子 β(JS)在兩類區域分別為O(1)和那么光滑因子 β(S)在兩類區域之間的差距更大.若權系數公式其他部分相同,則7 階WENO-S格式中含間斷子模板的權系數更小.這使得該格式具有良好的間斷穩定性.

2.4 分辨率性質

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

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

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

3 7 階WENO-S 格式的計算效率

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

3.1 線性對流方程的效率提升方法

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

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

考慮一維線性對流方程

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

簡便起見,這里用WENO7-JS 表示經典的7 階WENO 格式,而WENO7-Z 在WENO7-JS 的基礎上采用了Z 型權系數[28],其中冪因子取1,WENO7-S表示7 階WENO-S 格式.

表1 幾種WENO 格式計算N+1 個數值通量所需浮點運算量Table 1 The number of floating point operations required to calculate N+1 numerical fluxes for several WENO schemes

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

3.2 其他適用情況

上述效率提升方法也適用于其他一些問題.其關鍵在于一條網格線上用于WENO 重構或插值的量可以表示為一個序列,進行網格線上不同位置的WENO 重構或插值時,選取序列中相應連續幾個位置的值,然后進行WENO 重構或插值.

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

在稍復雜流動問題的數值模擬中,一些額外的處理過程可能會破壞這一條件.

(1) 通量分裂

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

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

(2) 特征投影

對于雙曲型偏微分方程組,由于信號沿特征方向傳播的特點,數值求解采用特征投影處理更符合方程的特點,在含間斷問題中結果更準確.對于非線性問題,投影時一般采用局部特征矩陣.

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

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

4 數值算例

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

4.1 極值點附近權系數偏離程度測試

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

計算中取該極值點作為網格點之一,統計該點左右各4 個通量點的情況.網格間距均取h=0.1,其結果放在表2 中.可以看到表中WENO7-S 的權系數偏離值在極值點相對更小.其中一階臨界點附近小一個量級甚至更多,這與第2.2 小節的結論相符.

表2 權系數偏離值的平均值與最大值Table 2 The average values and maximum values of the deviations of weights

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

4.2 一維對流問題

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

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

這里考慮波形的遠距傳輸,計算終止時間取T=2000,網格點數取400.為減小時間離散誤差的影響,CFL 數取0.01.計算結果顯示在圖2 中.對于這種遠距傳輸問題,WENO7-JS 耗散過大,而WENO7-Z 則在間斷附近出現了明顯的數值振蕩,但是在三角波峰值附近表現最好.WENO7-S 格式對于其他幾種波形均獲得了優良的結果.WENO7-JS,WENO7-Z 和WENO7-S 計算結果的L1誤差依次為2.51 × 10?1,4.67 × 10?2和4.79 × 10?2.此問題中WENO7-Z 的誤差最小,這可能是因為其三角波附近的優良表現和更小的間斷抹平程度.

圖2 組合波的遠距傳輸問題,網格點數N=400,計算終止時間T=2000Fig.2 Long-distance advection of combined waves with the number of grid points N=400 and the end time T=2000

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

表3 組合波遠距傳輸問題的計算時間Table 3 Computational time of long-distance advection of combined waves

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

4.3 球面波傳播問題

球面波傳播問題是一個計算氣動聲學的標準測試算例[36],其控制方程為

計算區域為[5,450],初始條件為u=0,在r=5處的邊界條件為u=sin(ωt).

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

圖3 球面波傳播問題,計算終止時間為T=400,空間步長為dr=1.時間步長為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

由于該問題計算時間過短,計算時間差距過小.因此,表3 給出了時間步長 dt取0.01 時的計算時間統計.結果表明,WENO7-S 計算效率較高,而消除冗余計算后,效率進一步提升21.82%.

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

4.4 二維旋轉問題

這個問題來自于Zalesak[37]的文章.該問題中,一個割開的圓柱繞著附近的軸旋轉.該運動的控制方程為

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

該問題中兩個分量速度既可能為正又可能為負,通常同時包含正負速度時需要進行通量分裂.但是該問題每條網格線上特征速度可確定為正或負,無需通量分裂處理.

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

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

圖4 二維旋轉問題計算結果與精確解Fig.4 The computational results and the exact solution of the twodimensional rotation problem

圖5 二維旋轉問題兩條線上的結果對比圖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 更容易出現上沖/下沖現象.而3 種格式的L1誤差依次為2.298,1.774 和1.959.說明盡管WENO7-Z 有較明顯的上沖/下沖,但是其誤差相對最小.從圖5(b)和圖5(c)中也可看出,WENO7-Z 的間斷抹平程度最小,而間斷附近貢獻了大部分誤差量,這使得其誤差小于WENO7-S.

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

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

表6 給出了200 × 200 網格旋轉5 個周期的計算時間對比.結果表明WENO7-S 計算效率較高,新的計算方法進一步提升了其計算效率,該問題中提升23.48%.

表6 二維旋轉問題的計算時間Table 6 Computational time for the two-dimensional rotation problem

4.5 二維小擾動傳播問題

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

Mx=0.5,My=0是x和y方向的平均流馬赫數.U′表示密度、速度及壓強的擾動量.計算區域為(x,y)∈[?100,100]×[?100,100].初值為

該問題為計算氣動聲學的一個標準算例[36,38].

根據迎風格式的要求,采用了每條網格線上的全局通量分裂.x方向計算時正負通量分別為方向計算時正負通量分別為

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

表7 二維小擾動傳播問題的計算時間Table 7 Computational time of two-dimensional small disturbance propagation problem

圖6 二維小擾動傳播問題計算結果,網格200 × 200,終止時間T=30Fig.6 Computational results of two-dimensional small disturbance propagation problem,with grid 200 × 200 and end time T=30

4.6 一維無黏流動問題

控制方程為一維Euler 方程

式中 ρ,u,p,E表示密度、速度、壓強和內能.氣體狀態方程為

其中比熱比 γ=1.4.

這里計算Shu-Osher 問題,該問題中馬赫數3 的激波與正弦波相互干擾,誘發高頻振動,其初值條件為

計算終止時間T=1.8,這里取網格點為200.

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

圖7 Shu-Osher 問題計算結果的密度對比Fig.7 The density distributions of computational results of Shu-Osher problem

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

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

表8 Shu-Osher 問題計算時間,N=2000Table 8 Computational time of Shu-Osher problem,N=2000

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

4.7 二維無黏流動問題

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

式 中 ρ,u,v,p,E表示密 度、x向速度、y向速度、壓強和內能.氣體狀態方程為

其中比熱比 γ=1.4.

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

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

采用組元型方法計算時,x方向正負通量分別為其 中λx為 該x向網格 線上所有矩陣 ?F/?U特征值絕對值的最大值.y方向計算時正負通量分別為其中λy為該y向網格線上所有矩陣 ?G/?U特征值絕對值的最大值.

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

圖8 激波旋渦相互作用問題,壓強等值線范圍1.19-1.37,共90 條,網格200 × 100,終止時間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 給出了幾種格式的計算時間.結果表明WENO7-S 格式效率相對較高,WENO7-S1 的效率提升方法進一步減少了21.66%的計算時間.

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

5 結論

WENO-S 格式具有許多優良的性質.特別是其光滑因子對單頻波為常數,使得其是一類滿足近似色散關系與線性基底格式一致的激波捕捉格式.另外,該格式能有效區分間斷與小尺度波動,在連續區域誤差小,間斷附近也能保證良好的穩定性,在數值模擬中能獲得良好的結果.

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

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

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

數值實驗結果表明7 階WENO-S 格式同時具有良好的小尺度結構分辨率和激波捕捉能力.在計算效率方面,本文所提方法能有效減少7 階WENOS 格式約20%的計算時間.