一類高效率高分辨率加映射的WENO 格式及其在復雜流動問題數值模擬中的應用

鐘 巍 賈雷明 王澍霏 田 宙

* (西北核技術研究所,西安 710024)

? (北京大學數學科學學院,北京 100871)

引言

針對包含激波或精細流場結構的復雜流動問題,往往要求數值模擬方法具備高精度和高分辨率.加權本質無振蕩 (weighted essentially non-oscillatory,WENO) 格式[1]在間斷附近擁有本質無振蕩 (essentially non-oscillatory,ENO) 特性,同時在光滑區域能夠獲得高階精度,因此,在復雜流動問題的數值模擬中廣受歡迎,并已經在諸多領域得到大量成功應用[2-25].

Liu 等[26]提出了第一個WENO 格式,他們在r階精度的ENO 格式基礎上,通過對其所有候選模板的權系數進行一個凸組合,獲得了精度為(r+1)階的WENO 格式.Jiang 等[1]通過定義新的光滑因子,提出了經典的WENO-JS 格式,該格式在普通光滑區域能夠獲得(2r? 1)階收斂精度.

Henrick 等[27]指出,WENO-JS 格式在光滑區域極值點附近會出現收斂精度丟失的問題.通過系統的理論分析,他們給出了WENO 格式的權系數在光滑區域極值點附近獲得最佳收斂精度的充分必要條件,并據此構造了一個關于WENO-JS 格式權系數的映射函數使得映射后的權系數滿足上述充分必要條件,相應的格式被稱為WENO-M.特別地,WENO-M格式在構造映射函數時采用的準則被學者們廣泛采用[28-32].然而,由于引入了映射操作,WENO-M 格式的計算效率受到了嚴重的影響,相對WENO-JS 格式而言,其額外計算時間增加了25%左右[33].不僅如此,Feng 等[28]發現WENO-M 格式的映射函數在權系數為0 的附近會放大不光滑模板的作用,這被認為是造成WENO-M 格式在長時間模擬時出現精度丟失的關鍵原因.于是,他們新增了一個限制條件,即要求映射函數在0 附近的導數值為0.為此,他們設計了一個分片多項式映射函數,相應的格式被稱為WENO-PM6.WENO-PM6 格式具有非常低的數值耗散,因此,其在長時間模擬時仍然能夠獲得非常高的分辨率,且對于精細流場結構的捕捉能力遠遠高于WENO-JS 和WENO-M 格式.但是,它引入的分片多項式映射函數帶來了更加復雜的數學運算操作,這使得其計算效率比WENO-M 格式還要低.本文后面給出數值算例測試結果表明,WENO-PM6 格式相對WENO-JS 格式的額外計算時間增加了50%以上.

除了前面提到的WENO-M 和WENO-PM6,學者們還提出了許多各種形式的加映射的WENO 格式[29-32].顯然,這些加映射的WENO 格式都能夠很好的克服WENO-JS 格式在極值點附近計算精度丟失的問題,擁有比WENO-JS 格式更低的耗散從而獲得更高的分辨率.無一例外的,這些加映射的WENO格式都存在計算效率下降這個普遍的缺陷.

本文巧妙地設計了一個全新的映射函數,其滿足現有映射函數的所有設計要求,但在絕大部分區域僅需要一次簡單的賦值運算即可完成映射操作.換句話說,新映射函數不會帶來嚴重的額外計算時間消耗,故而具有非常高的計算效率.此外,相應的加映射的WENO 格式還擁有非常低的數值耗散,因此,與已有加映射的WENO 格式相比,計算分辨率得到極大的提高.

1 問題描述與WENO 格式回顧

1.1 問題描述

無黏、可壓縮流動的控制方程為歐拉方程組.三維歐拉方程組為

其中

式中,ρ,u,v,w,p,E分別表示流體密度,x,y,z方向的速度分量,壓強和單位體積流體的總能量.以下理想氣體狀態方程與方程(1)組成封閉方程組

其中,γ 為絕熱指數,通常取 γ=1.4.

1.2 WENO 格式回顧

通常情況下,重構都是在特征空間中進行然后將結果投影回物理空間[1].因此,為了方便描述且不失一般性,針對以下一維標量雙曲守恒律方程展開討論

下面簡單介紹5 階WENO-JS、WENO-M 和WENO-PM6 這3 種格式的具體實現方法.

5 階WENO 重構通過對上述結果做加權平均得到uj+1/2,即

其中,ωk是非線性權系數.

在WENO-JS 格式中,ωk按下式計算

式中,k=0,1,2 (下同),參數 ε 是一個非常小的正數,用來防止分母為0,(C0,C1,C2)=(0.1,0.6,0.3) 是理想權系數,I Sk被稱為光滑因子,用來度量子模板的光滑程度,其具體表達式如下

Henrick 等[27]發現WENO-JS 格式在極值點附近無法獲得最佳收斂精度,于是他們提出了WENO-M格式,其權系數按下式計算

1.3 時間離散

令

則可以采用合適的時間推進方法求解半離散方程組

本文采用3 階顯式強穩定總變差不增的Runge-Kutta 格式[35-37]進行時間推進,即

2 加近似常值映射的WENO 格式

2.1 新映射函數

先定義標準符號函數的一個光滑近似函數sg(x,μ,T,m),其表達式為

其中,m為正整數,其值越大s g(x,μ,T,m) 越接近實際符號函數;μ 為大于0 的實數且 μ→0 ,其值越小sg(x,μ,T,m)越接近實際符號函數;參數T為用于調整函數形狀的一個正實數,其值越小sg(x,μ,T,m) 越接近實際符號函數.在后文中將會看到,sg(x,μ,T,m)越接近實際符號函數對應的格式的計算效率會越高.數學上容易證明,sg(x,μ,T,m) 是一個單調遞增的光滑函數.

利用sg(x,μ,T,m) 可以構造一個全局單調遞增的映射函數,形式如下

至此,可以定義一個全新的加全局單調遞增光滑映射的5 階WENO 格式,記為WENO-ACM,其權系數按下式計算

作為示例,圖1 中分別給出了取m=3,T=30,μk=0.13,ζk=Ck/3和m=2,T=20,μk=10?6,ζk=0.1Ck時,5 階WENO-ACM 格式對應C k=0.6 的映射函數變化曲線.作為比較,在圖中同時給出了WENO-M和WENO-PM6 格式相應的映(射)函數變化曲線.為了方便討論,不妨稱映射結果由0 變為Ck(或由C k變為1)對應的 ω 所在區間為過渡區間.顯然,結合圖1(a)很容易完成對式(10)中性質的驗證.實際上,圖1(b)中的WENO-ACM 的映射函數同樣滿足式(10)中的性質,只是其過渡區間非常小.不過,非常重要的是,通過合理設置參數,在滿足光滑性等要求的前提下,盡可能減小過渡區間是非常有必要的.由式(9)可以知道,在過渡區間以外的區間僅進行一次賦值運算,而在過渡區間則要進行復雜的數學運算操作.因此,過渡區間越小,的映射過程所需數學運算操作數也越少,相應的計算效率也會越高.當過渡區間足夠小的時候,在絕大部分情況下都只進行一次賦值操作即完成了整個映射過程.因此,可以將看作是一個近似常值映射函數.另外,可以看到,WENO-M 和WENO-PM6 格式的過渡區間都非常大.

圖1 WENO-ACM 格式的映射函數曲線Fig.1 Mapping functions of WENO-ACM

2.2 收斂精度分析

對于5 階WENO-JS 格式,其權系數滿足[27]

而在1 階極值點處,上述關系式進一步退化為

由式(10)可得

2.3 頻譜特性分析

利用Pirozzoli[38]提出的近似色散關系(ADR),可以對非線性格式進行頻譜特性分析.ADR 分析通過求解模型方程,對給定網格上的全部Fourier 模式進行求解,揭示格式在整個頻域的數值表現.對于波數為 φ 的正弦波,從數值結果換算得到的修正波數為 Φ (φ),Φ (φ) 的具體計算方法在文獻[38-39]中有非常詳細的介紹.Φ (φ) 是一個復數,其實部反映格式的色散性質,虛部則反映格式的耗散性質.

圖2 中給出了幾種數值格式的ADR 分析結果,其中UW5 是5 階線性迎風格式,WENO-ACM-(a)和WENO-ACM-(b) 格式的參數分別與圖1(a) 和圖1(b)中一致.容易看到,WENO-ACM-(a)的色散和耗散性能均明顯不如WENO-ACM-(b),也比WENO-M 和WENO-PM6 要差.因此,接下來的內容中本文只對WENO-ACM-(b)進行討論,且為了簡化書寫,將用WENO-ACM 表示WENO-ACM-(b).由圖2(a)可見,在中波數區間,WENO-ACM 格式相比WENO-JS,WENO-M 和WENO-PM6 格式的色散是有所降低的,但在波數大于2.3 以上的高波數區間,WENO-ACM 格式的色散性能得到了提升.這里重點關注格式的耗散性質,如圖2(b)所示,在中波數范圍內,WENO-ACM 格式的耗散性質與WENO-M和WENO-PM6 格式非常接近,且均優于WENO-JS格式,而在波數較大的區間,WENO-ACM 格式相對其他格式則表現出了更好的耗散性能.

圖2 不同WENO 格式的頻譜特性曲線比較Fig.2 Comparison of the spectral properties of different WENO schemes

上述ADR 分析結果表明,WENO-ACM 格式是一種低耗散的數值格式.在接下來的一節中,將通過大量的數值算例來證實這一點.

3 數值算例與結果分析

本節通過豐富的數值算例進一步表明本文所提出的新格式在計算效率和低耗散性方面顯著的優勢.為了更好的比較,將同時給出WENO-JS,WENO-M,WENO-PM6 這3 種非常有代表性格式的計算結果.

在經過大量的數值實驗后,本文發現選取如圖1(b)中所示參數的時候,WENO-ACM 格式的數值表現是最好的.因此,接下來的部分,WENOACM 格式均指采用如圖1(b)中所示參數的情況.同時,根據數值經驗,本文建議在實際使用WENOACM 格式時選取這組參數,即在選擇參數時確保映射函數的過渡區盡可能小且理想權區間盡量大.

3.1 精度測試

利用一維線性對流方程

進行精度測試.采用周期邊界條件,輸出時間為t=2.分別考慮下面兩種不同的初始條件

為了確保時間收斂精度與空間收斂精度一致,取CFL 數為 ?x2/3.利用特征線法,容易求得上述兩種初始條件下問題的精確解分別為

于是,可得數值格式的L2和L∞誤差分別為

表1 給出了不同格式對初始條件1 的L2和L∞誤差及其收斂精度.由表1 可知,所有格式均能夠獲得最佳收斂精度.整體上,WENO-M,WENOPM6 和WENO-ACM 的絕對計算誤差是相當的,且要明顯小于WENO-JS.不同格式對初始條件2 的L2和L∞誤差及其收斂精度計算結果如表2 所示.可以看到,WENO-JS 的L∞收斂精度只有3 階左右,其L2收斂精度為下降到4 階以下,而WENO-M,WENOPM6 和WENO-ACM 則仍然能夠獲得最佳收斂精度,且它們的絕對計算誤差明顯高于WENO-JS 格式.

表1 不同WENO 格式的 L2 和 L ∞ 誤差及其收斂精度 (初始條件1)Table 1 L 2 and L ∞ errors and the convergence properties of various schemes for initial condition 1

表2 不同WENO 格式的 L2 和 L ∞ 誤差及其收斂精度 (初始條件2)Table 2 L 2 and L ∞ errors and the convergence properties of various schemes for initial condition 2

為了表明WENO-ACM 格式在計算效率方面的優勢,圖3 和圖4 分別給出了不同格式針對初始條件1 和初始條件2 時的L2和L∞計算誤差與CPU 計算時間關系曲線.由圖可見,在獲得相同計算精度的情況下,WENO-ACM 所消耗的CPU 計算時間是最小的,這表明WENO-ACM 的計算效率是最高的.另外,可以發現,在獲得相同的計算精度時,WENOPM6 的CPU 計算時間要大于WENO-M,這表明WENO-PM6 的映射操作帶來的額外計算時間要高于WENO-M,從而導致其計算效率出現了嚴重的下降.在本文最后一小節通過二維數值算例專門討論上述格式的計算效率問題時,將會再次證實這一結論.

圖3 不同WENO 格式對初始條件1 的計算誤差與CPU 計算時間關系曲線Fig.3 Comparison of various WENO schemes for initial condition 1 in CPU time and numerical computing errors

圖4 不同WENO 格式對初始條件2 的計算誤差與CPU 計算時間關系曲線Fig.4 Comparison of various WENO schemes for initial condition 2 in CPU time and numerical computing errors

3.2 長時間模擬測試

為了驗證WENO-ACM 格式長時間計算時在穩定性和低耗散性方面的優勢,選取一維線性對流方程進行測試.采用周期邊界條件,初始條件如下

取輸出時間為t=1200,CFL 數為 ?x2/3.利用特征線法,容易求得問題的精確解為

圖5 給出了300 個均勻網格下不同WENO 格式的計算結果.由圖可見,WENO-JS 和WENO-M 均出現了非常嚴重的分辨率下降,這是由于它們的數值耗散太大引起的.WENO-ACM 和WENO-PM6 則在長時間模擬時仍然能夠獲得非常高的分辨率,這表明它們的數值耗散非常小.此外,圖6 給出了不同WENO 格式對該問題的L2和L∞計算誤差與CPU計算時間關系曲線.圖中結果表明,在獲得相同計算精度的情況下,WENO-ACM 所消耗的CPU 計算時間是最少的,這表明其計算效率是最高的.可以看到,WENO-PM6 雖然計算精度高,但其消耗的CPU 計算時間非常大,導致其計算效率甚至比WENO-M還低.

圖5 不同WENO 格式對包含高階臨界點光滑問題的長時間模擬結果,N=300, t=1200Fig.5 Results of various WENO schemes for the smooth problem with high-order critical points at long output time,N=300, t=1200

圖6 不同WENO 格式對高階臨界點光滑問題的計算誤差與CPU 計算時間關系曲線Fig.6 Comparison of various WENO schemes for the smooth problem with high-order critical points in CPU time and numerical computing errors

3.3 一維歐拉方程組

本節對一維歐拉方程組進行模擬,方程組的具體表達形式由式(1)和式(2)很容易得到.在進行有限體積框架下WENO 重構的時候,采用局部特征分解技術,具體實現過程可以參考文獻[1,40].所有算例,CFL 數均取0.5,初始條件為U0=(ρ,u,p)(x,0).

3.3.1 激波管問題

先考慮兩個標準的激波管問題,即Sod 和Lax激波管問題,其初始條件如下

采用透射邊界條件,取200 個均勻網格進行模擬,圖7 給出了WENO-JS,WENO-M,WENOPM6 和WENO-ACM 格式將Sod 和Lax 激波管問題分別計算到t=0.25 和t=1.3 的結果(密度).由圖可見,WENO-ACM 和WENO-M,WENO-PM6 均給出了比WENO-JS 更高的分辨率.更仔細的觀察可以發現,WENO-ACM 的模擬效果要輕微地優于WENO-M 和WENO-PM6,這表明WENO-ACM 的數值耗散是最小的.

圖7 不同WENO 格式對Sod 和Lax 激波管問題的模擬結果Fig.7 Results of various WENO schemes for the Sod and Lax shocktube problems

3.3.2 爆炸波碰撞問題

計算由Woodward 等[43]提出的一維爆炸波碰撞問題,其初始條件如下

采用反射邊界條件,取400 個均勻網格進行模擬,圖8 給出了WENO-JS,WENO-M,WENOPM6 和WENO-ACM 格式將該問題計算到t=0.038的結果(密度),其中參考解是取10 000 個網格時WENO-JS 的計算結果.由圖可見,WENO-JS 的分辨率明顯低于其他幾種格式,而WENO-ACM 的分辨率則是所有格式中最高的,雖然只是輕微的好于WENO-M 和WENO-PM6 格式.同樣地,上述結果證明WENO-ACM 具有非常低的數值耗散.

圖8 不同WENO 格式對Woodward-Colella 爆炸波碰撞問題的模擬結果Fig.8 Results of various WENO schemes for the Woodward-Colella blast wave interacting problem

3.3.3 激波?熵波相互作用問題

這里模擬兩個經典的激波與熵波相互作用問題,其初始條件如下

采用透射邊界條件,分別取300 和1500 個均勻網格進行模擬,圖9 和圖10 依次給出了WENO-JS,WENO-M,WENO-PM6 和WENO-ACM 格式將Shu-Osher 和Titarev-Toro 激波?熵波相互作用問題分別計算到t=1.8 和t=5.0 的結果(密度),其中參考解是取10 000 個網格時WENO-JS 的計算結果.由圖9(b)可見,在高頻光滑波區域,由于數值耗散太大,WENO-JS 的分辨率是最低的,而得益于低耗散的優點,WENO-ACM 的分辨率是最高的,盡管此時其分辨率只是略微的高于WENO-M 和WENOPM6.然而,由圖10(b) 可見,當面對更為苛刻的Titarev-Toro 問題時,WENO-ACM 的分辨率要顯著的高于所有其他格式.此時,WENO-JS 的分辨率仍然是最低的,而WENO-PM6 的分辨率也要明顯的高于WENO-M.

圖9 不同WENO 格式對Shu-Osher 激波-熵波相互作用問題的模擬結果Fig.9 Results of various WENO schemes for the Shu-Osher shockentropy interaction problem

圖10 不同WENO 格式對Titarev-Toro 激波-熵波相互作用問題的模擬結果Fig.10 Results of various WENO schemes for the Titarev-Toro shockentropy interaction problem

3.4 二維歐拉方程組

接下來,對二維歐拉方程組進行模擬.同樣的,方程組的具體表達形式由式(1)和式(2)很容易得到.仍然在有限體積框架下進行求解,逐維進行WENO 重構,并采用局部特征分解技術.對所有算例,CFL 數均取0.5,在沒有特別說明的情況下,絕熱指數取 γ=1.4,初始條件為U0=(ρ,u,v,p)(x,y,0).

3.4.1 二維黎曼問題

首先模擬非常經典的二維黎曼問題[48-50],其計算區域為[0,1] × [0,1].已有研究表明,二維黎曼問題有很多不同的解結構[49],這里選取其中一種結構進行模擬,其初始條件如下

4 條邊界均取透射邊界條件,使用400 × 400 個均勻網格,輸出時間為t=0.8.

圖11 中給出了WENO-JS,WENO-M,WENOPM6 和WENO-ACM 計算得到的密度等值線圖.由圖可見,WENO-JS 的分辨率最低,其次是WENO-M.WENO-PM6 的分辨率要略高于WENO-M,但要明顯低于WENO-ACM.這表明,WENO-ACM 的數值耗散在所有格式中是最小的.

圖11 不同WENO 格式對二維黎曼問題的模擬結果Fig.11 Results of various WENO schemes for the 2D Riemann problem

同時,還可以觀察到,WENO-M,WENOPM6 和WENO-ACM 均出現了比較明顯的數值振蕩——如圖中粉色虛線框標注所示.分析原因如下(為方便描述,稱映射函數在 ωk=Ck附近用理想權系數代替原權系數對應的 ωk所在的區間為理想權區間): 文獻[51]中已經指出,更大的理想權區間是獲得更高分辨率的關鍵,但增大理想權區間同時也提高了放大非光滑子模板上權系數的風險進行導致產生數值振蕩.由圖1(b)容易發現,這里所考慮的格式的理想權區間大小滿足關系: WENO-ACM >WENO-PM6 >WENO-M >WENO-JS,而圖11 中數值分辨率高低以及數值振蕩的大小同樣滿足上述關系.顯然,這與文獻[51]的結論是一致的.盡管如此,這些數值振蕩并未影響流場主要波系結構的捕捉,因此,當以追求高分辨率為主要目標時,它們是可以容忍的.這種情況在后面的算例(特別是前臺階流動算例)中同樣非常顯著,為了節省篇幅,后文中將不再對此進行贅述.

3.4.2 雙馬赫反射問題

雙馬赫反射問題[43]被廣泛用于測試高分辨率格式的性能.該問題的計算區域為[0,4] × [0,1],初始時刻計算區域的流場分布如下

其中,x0=1/6.在左側入口處采用入流邊界條件,流場狀態與初始條件一致,在右側出口處采用出流邊界條件.在底部 [ 0,x0] 區域采用初始條件,[x0,4] 區域采用反射邊界條件.在頂部,其邊界條件設定如下

圖12 中給出了WENO-JS,WENO-M,WENO-PM6和WENO-ACM 計算得到的密度等值線圖,其中粉色盒子標注的滑移線附近的不穩定渦結構及主激波后方的伴隨結構通常用來衡量數值格式的耗散性能.由圖可見,所有格式均能夠很好的捕捉到主激波后方的伴隨結構.WENO-JS 對滑移線附近的不穩定結構捕捉效果最差,因此其數值耗散最大,然后是WENO-M.WENO-PM6 和WENO-ACM 對滑移線附近不穩定結構的捕捉效果相差不大,均明顯優于WENO-JS 和WENO-M.顯然,本算例表明WENO-ACM具有非常低的數值耗散.

3.4.3 前臺階流動問題

前臺階流動問題[43]在高分辨率格式的數值測試中同樣很受歡迎,特別是對發源于馬赫桿處的物理不穩定現象即漩渦面卷曲現象的捕捉,對數值格式的低耗散性能提出了非常高的要求[13,52-55].該問題的計算區域為 ?=[0,0.6]×[0,1]∪[0.6,3]×[0.2,1],即在一個長為3、高為1 的小型風洞內距離風洞左側入口0.6 的位置放置一個長為2.4、高為0.2 的臺階.整個計算區域在初始時刻的流場狀態如下

在風洞的左側入口、右側出口處分別采用入流邊界條件、出流邊界條件,在風洞及臺階的壁面都采用反射邊界條件,計算輸出時間為t=4.0.

圖13 和圖14 分別給出了網格數取900 ×300 和1800 × 600 時,WENO-JS,WENO-M,WENOPM6 和WENO-ACM 計算得到的密度等值線圖.首先,對于上述兩種不同網格尺寸,所有格式均能夠很好的捕捉到主要的波系結構.由圖13 可見,網格數取900 × 300 時,只有WENO-ACM 能夠模擬出發源于馬赫桿處的漩渦面卷曲現象.當將網格加密到1800 × 600 時,所有格式均能夠模擬出上述漩渦面卷曲現象,但WENO-JS 的模擬結果最不顯著.WENO-M 和WENO-PM6 都能夠觀察到非常清晰的漩渦面卷曲結構,但從尺寸上看,WENO-ACM 的模擬效果是最佳的.本算例非常有力的證明了WENOACM 在低數值耗散方面的優勢.

圖13 不同WENO 格式對前臺階流問題的模擬結果,900 × 300個網格Fig.13 Results of various WENO schemes for the forward facing step problem,900 × 300 cells

圖14 不同WENO 格式對前臺階流問題的模擬結果,1800 × 600個網格Fig.14 Results of various WENO schemes for the forward facing step problem,1800 × 600 cells

3.4.4 瑞利?泰勒不穩定性問題

本算例刻畫初始流場中兩種不同密度的流體在重力作用下上方密度較大的流體加速進入下方密度較小的流體流動不穩定性現象[56].計算區域為[0,0.25] × [0,1],初始時刻計算區域的流場分布如下

左、右側均取反射邊界條件,上、下則采用常值邊界條件,即

需要注意的是,該問題的絕熱指數取為 γ=5/3,且在二維歐拉方程組的右側要添加源項S=(0,0,ρ,ρv)T來體現重力的影響.

圖15 給出了網格數取240 × 960 時,WENO-JS,WENO-M,WENO-PM6 和WENO-ACM 對該問題模擬到t=1.95 時的密度等值線圖.可以非常清晰的觀察到,WENO-ACM 的分辨率要顯著的高于其他幾種格式,它能夠識別出最為豐富的渦結構,這表明它的數值耗散是最小的.而數值耗散最大的WENOJS 模擬結果的分辨率是最低的,隨后是WENO-M.WENO-PM6 雖然分辨率要高于WENO-JS 和WENO-M,但仍遠遠低于WENO-ACM.

圖15 不同WENO 格式對瑞利-泰勒不穩定性問題的模擬結果Fig.15 Results of various WENO schemes for the Rayleigh-Taylor instability

3.4.5 開爾文?亥姆霍茲不穩定性問題

最后,對開爾文?亥姆霍茲不穩定性問題[57-60]進行模擬.該問題的計算區域為[0,1] × [0,1],初始時刻計算區域的流場分布如下

其中

四周均取透射邊界條件,輸出時間為t=0.8.

圖16 給出了網格數取512 × 512 時,WENO-JS,WENO-M,WENO-PM6 和WENO-ACM 模擬得到的密度等值線圖.由圖可見,WENO-JS 模擬得到的渦的尺寸是最小的,其次是WENO-M.WENOPM6 模擬得到的渦的尺寸要明顯大于WENO-M,但仔細觀察發現要略微的小于WENO-ACM.

圖16 不同WENO 格式對開爾文-亥姆霍茲不穩定性問題的模擬結果Fig.16 Results of various WENO schemes for the Kelvin-Helmholtz instability

3.5 計算效率測試

前面的算例已經充分證明WENO-ACM 在低耗散性上與其他幾種WENO 格式相比有著顯著的優勢.現在,仍然利用上一小節的算例,對WENOACM 的計算效率進行專門測試.為了更好的比較,這里對WENO-JS,WENO-M 和WENO-PM6 也進行效率測試.測試時的基本計算條件與上一小節一致,但新增了一組網格條件.測試程序采用 C++代碼編寫,并在個人計算機上運行,計算機CPU 為 Intel(R)Core(TM)i9-9880 H.為了盡量排除其他因素的影響,測試時僅考慮單個Runge-Kutta 迭代步的CPU 時間消耗.同時,為了消除隨機性的影響,每次測試均在相同條件下重復了3 次.

表3～表7 給出了效率測試的結果,其中第一至四行給出了每個Runge-Kutta 迭代步的CPU 時間消耗,并在括號內給出了相對WENO-JS 的額外時間消耗.同時,為了更好的比較,在每個表格的最后兩行還以百分比的形式分別給出了WENO-ACM 相對WENO-M 和WENO-PM6 的“額外計算時間消耗減少量”.由表可見: WENO-M 和WENO-PM6 相對WENO-JS 的額外時間消耗分別大于 24%和54%,而WENO-ACM 相對WENO-JS 的額外時間消耗則不超過5%;WENO-ACM 相對于WENO-M 和WENO-PM6的“額外計算消耗減少量”分別大于80%和90%?遠遠高于WENO-PDM[61]相對WENO-M 的52%的“額外計算消耗減少量”.這表明,WENO-ACM 能夠極大的提高加映射的WENO格式的計算效率.

表3 每個Runge-Kutta 步的CPU 計算時間(s)及相對WENO-JS 的額外計算時間(二維黎曼問題)Table 3 CPU time (s) and the extra computational cost per Runge-Kutta step for the 2D Riemann problem

表4 每個Runge-Kutta 步的CPU 計算時間(s)及相對WENO-JS 的額外計算時間(雙馬赫反射問題)Table 4 CPU time (s) and the extra computational cost per Runge-Kutta step for the double Mach reflection problem

表5 每個Runge-Kutta 步的CPU 計算時間(s)及相對WENO-JS 的額外計算時間(前臺階流問題)Table 5 CPU time (s) and the extra computational cost per Runge-Kutta step for the forward facing step problem

表6 每個Runge-Kutta 步的CPU 計算時間(s)及相對WENO-JS 的額外計算時間(瑞利-泰勒不穩定性問題)Table 6 CPU time (s) and the extra computational cost per Runge-Kutta step for the Rayleigh-Taylor instability

為了獲得更清晰的計算時間消耗分布,這里進一步詳細統計程序計算所有各部分的CPU 時間,即分別統計特征分解以及特征空間與物理空間來回投影部分(Tch)、重構部分(Trec)、映射部分(Tmap)、其他部分(Toth)的CPU 時間.為了節省篇幅但又不失一般性,以400 × 400 網格的二維黎曼問題和1000 ×250 網格的雙馬赫反射問題為例,表8 和表9 中分別給出了一個Runge-Kutta 迭代步中不同格式上述各部分CPU 時間的統計結果,并在括號中給出了各部分CPU 時間所占百分比.

表8 每個Runge-Kutta 步的CPU 計算時間(s)分布情況(二維黎曼問題)Table 8 CPU time (s) allocation per Runge-Kutta step for the 2D Riemann problem

表9 每個Runge-Kutta 步的CPU 計算時間(s)分布情況(雙馬赫反射問題)Table 9 CPU time (s) allocation per Runge-Kutta step for the double Mach reflection problem

由表可見: 顯然,重構部分消耗的CPU 時間占比是最大的;WENO-M 的映射消耗時間(Tmap) 在25%以上,而WENO-PM6 則在40%以上且與重構部分消耗的CPU 時間Trec差不多;和預期的一樣,WENO-ACM 的映射消耗時間(Tmap)非常低,僅為3%～ 5%.上述結果表明,傳統加映射的WENO 格式由于映射操作消耗的CPU 時間是非常驚人的,這一結論與文獻[33]中給出的測試結果是一致的.此外,還可以看到: 其他部分(Toth)占用的時間是最少的且幾乎是恒定的,約0.063 s;二維黎曼問題的特征分解以及特征空間與物理空間來回投影部分(Tch)占用的時間在0.2 s 左右,而雙馬赫反射問題則在0.24～0.27 s 左右.

3.6 具體復雜問題驗證

最后,對WENO-ACM 模擬具體的復雜問題進行驗證,這里選擇豎直平面激波與熱層作用[15]問題進行模擬.

問題描述: 如圖17 所示,在[?5 m,5 m] × [?5 m,5 m]的區域布滿了理想完全氣體.流場中存在物質界面P1OP2(圖中虛線所示),點O與坐標原點重合,點P2固定在(0 m,?5 m)處,點P1位于x=5 m 上,坐標記為(5 m,y0).初始時刻,整個流場內的介質處于靜止狀態,且各處壓力均為100 kPa,比熱比均為1.4.界面左側氣體密度為1.00 kg/m3,右側的氣體被加熱形成形狀不規則的熱層,因此其密度降低為0.25 kg/m3.馬赫數2.0 且強度恒定的平面激波沿x軸正向從左向右傳播,假設其到達界面OP2處時的時刻為t=0.已知平面激波后流場壓力為450 kPa,密度為2.667 kg/m3,速度467.7 m/s,比熱比為1.4.除左側為入流邊界外,其他三側均為出流邊界條件.現在取y0=2.886 8 m,模擬t=4.0 ms 時的流場分布狀態.

圖17 豎直平面激波與熱層作用示意圖Fig.17 Illustration of the interaction of a vertical planar shock with a thermal layer

針對這一典型的涉及復雜流場的力學問題,采用成熟商業軟件Autodyn 及本文提出的WENOACM 格式進行了模擬,其中Autodyn 和WENOACM 分別采用了400 萬和25 萬的均勻網格,圖18中給出了密度的模擬結果.作為對比,圖中還給出了相同計算條件下WENO-M 和WENO-PM6 的模擬結果.由圖18(a) 和圖18(d) 可見,盡管WENOACM 的網格數量只有Autodyn 的1/16,但其模擬結果與Autodyn 是相當的,甚至在局部位置(如圖18(a)中紅色虛線框區域) WENO-ACM 的模擬分辨率要高于Autodyn.可以看到,在圖18 中WENO-M 和WENO-PM6 的模擬結果與WENO-ACM 非常接近,幾乎難以直接觀察到它們之間的區別.

圖18 Autodyn 和WENO-ACM 計算所得流場密度分布Fig.18 Density distribution computed by Autodyn and WENO-ACM

此外,如圖18(a)中粉色點劃線所圍區域所示,Autodyn 的模擬結果在激波后方出現了明顯的數值擾動現象,而圖18(d)所示的WENO-ACM 的模擬結果中似乎不存在這樣的擾動.為了進一步驗證這一點,在圖19 中畫出了Y=4 截面上Autodyn,WENOACM 以及WENO-M 和WENO-PM6 的密度變化曲線,其中參考解是采用1000 × 1000 網格的WENOJS 格式的模擬結果.由圖19 可以清楚的看到,在激波后方,Autodyn 出現了顯著的數值擾動,而WENO-ACM 則成功的避免了這些數值擾動.

由圖19 可知,WENO-M 和WENO-PM6 同樣能夠避免上述數值擾動,但仔細觀察可以看到,WENO-M 的分辨率要低于WENO-ACM,WENOPM6 的分辨率非常接近但略微低于WENO-ACM.此外,如前所述,WENO-ACM 相比于WENO-M 和WENO-PM6 的另一個顯著的優勢是其具備更高的計算效率.為了驗證這一點,針對本算例同樣進行了計算效率測試,結果如表10 所示.由表10 可見,WENO-ACM 相對WENO-JS 的額外計算時間控制在5% 以內,而WENO-M 和WENO-PM6 相對WENO-JS 的額外計算時間則分別達到了30% 及70%,遠遠高于WENO-ACM 的結果.顯然,這一結果與前面的結論是一致的.

圖19 Y=4 截面上的密度變化曲線比較Fig.19 The cross-sectional slices of density plot along the plane Y =4

表10 每個Runge-Kutta 步的CPU 計算時間(s)及相對WENO-JS 的額外計算時間(豎直平面激波與熱層作用問題)Table 10 CPU time (s) and the extra computational cost per Runge-Kutta step for the interaction of a vertical planar shock with a thermal layer

以上結果成功的驗證了WENO-ACM 格式在模擬具體的復雜問題時,具有很好的魯棒性以及非常高的分辨率.同時,相對文獻中已有的加映射的WENO 格式,WENO-ACM 格式在計算效率方面有明顯的優勢.

4 結論

本文給出了一種全新的加近似常值映射的WENO 格式,稱為WENO-ACM.借助符號函數的一種近似函數,構造出了一種全新的映射函數.該映射函數完全滿足傳統映射函數構造過程中關鍵的約束條件,因此,其相應的WENO-ACM 格式能夠獲得最佳收斂精度.最重要的是,由于在大部分情況下,新的映射函數只需要一次賦值運算操作,因此可以極大的降低計算消耗,進而提高計算效率.數值測試結果表明,WENO-ACM 能夠極大的減少相對WENOJS 格式的額外計算消耗.此外,本文推薦的WENOACM 還能夠顯著降低數值耗散,進而提高計算結果的分辨率.