楊 強, 李偉鵬

( 上海交通大學 航空航天學院,上海 200240)

可壓縮流動廣泛存在于超聲速/高超聲速飛行器中,激波間斷和湍流共存且相互耦合作用給高可靠性數(shù)值模擬帶來了巨大挑戰(zhàn).近40年,人們發(fā)展了眾多高精度激波捕獲方法,例如總變差不增(TVD)方法[1]、ASUM[2]、基本無振蕩/加權基本無振蕩(ENO/WENO)方法[3-5]、加權非線性緊致格式(WCNS)方法[6]、基于目標本質無振蕩(TENO)方法[7]和AFTENO[8]等.這些格式在激波間斷附近增加局部耗散來抑制非物理數(shù)值振蕩,而在光滑湍流區(qū)域減少耗散以保持對多尺度湍流結構的分辨率.經典激波捕捉格式,如TVD方法和ENO方法,雖然可以解決激波間斷無數(shù)值振蕩,但對湍流光滑區(qū)域產生了過多的數(shù)值耗散.WENO方法在ENO方法的基礎上,引入子模板光滑因子,通過非線性加權方法來動態(tài)調整局部耗散[5],使得高精度數(shù)值格式得到進一步發(fā)展.Borges等[9]在非線性權值中引入更高階的光滑指標,提出計算成本更低、精度更高的WENO-Z格式.雖然WENO格式的精度可以設計為任意高階[10-13],但WENO格式對光譜區(qū)域的分辨率仍然不能令人滿意[14],它的加權機制和光滑度指標設定依舊對小尺度湍流結構中引入過多的數(shù)值耗散[15].1999年,Deng等[16]提出非線性緊致格式(CNS),后續(xù)引入加權技術,構造了 WCNS[6].近年來,Fu等[7]提出了另一類非線性加權框架,即TENO 方法,通過類ENO方法的模板選擇來控制非線性耗散,在無間斷區(qū)域(包括光滑臨界點)保證最優(yōu)的高階空間重構,在間斷區(qū)域通過移除跨越間斷的模板來避免數(shù)值振蕩,這種方法有利于候選模板的使用和高階格式的改進.當子模板中存在多個相鄰的不連續(xù)間斷,該方法可退化為三階重構.即使對于高階版本,數(shù)值魯棒性也得到了明顯改善,且避免了階躍退化問題.Peng等[17]在此基礎上進一步改進了TENO格式的CT自適應函數(shù)策略,發(fā)展了TENO5-LAD格式,實驗結果證明新方法能夠更好地抑制不連續(xù)點附近的數(shù)值振蕩,并能以較低的額外計算成本保持TENO的低耗散特性.而Ye等[8]結合AFWENO格式和TENO格式的優(yōu)點構建了可以實現(xiàn)任意高階精度的AFTENO格式,在激波捕獲和波數(shù)分辨率方面,相較于AFWENO-JS和AFWENO-Z有著明顯的優(yōu)勢.之后,在加權緊致格式方面,Zhang等[18]將WCNS與TENO格式組合發(fā)展了TCNS,實驗結果顯示TCNS能夠在模擬中捕獲更豐富的波結構,特別是在對瑞利-泰勒不穩(wěn)定性問題的模擬中,TCNS在較粗網(wǎng)格上取得了與WCNS-JS[5]相似甚至更好的結果.

將TENO格式的模板選擇和加權策略引入到WCNS的構造過程,發(fā)展了WCNS-T格式,并在此基礎上測驗7階精度下的數(shù)值表現(xiàn)情況.新格式在擴展模板寬度的同時改變了模板選擇過程,當全局模板存在多個間斷時,WCNS7-T格式將不會使用所有候選子模板,通過引入截斷函數(shù),將光滑和非光滑候選模板完全分離,并將非光滑候選模板的權重直接設為0,徹底剔除非光滑模板的影響.該模板策略能使格式在高階精度時相較原來能保持更高的數(shù)值穩(wěn)健性.同時,光滑度因子的新評測策略較之前能更準確判定子模板的光滑程度,降低了TENO格式中復雜光滑指標所帶來的不必要的計算量,提高了格式捕捉激波的能力.一維和二維數(shù)值算例結果表明,WCNS7-T格式相比于傳統(tǒng)的WCNS7-Z格式具有更高的精度和分辨率.

1 數(shù)值方法

1.1 控制方程

針對黎曼問題,差分格式的推導僅考慮一維守恒型Euler方程:

?Q/?t+?E/?x=0

(1)

守恒變量Q和E為

(2)

式中:ρ為密度;u為來流速度;p為壓力;γ為絕熱指數(shù).

在每個網(wǎng)格節(jié)點xi處,對式(1)進行有限差分格式離散得:

(3)

(4)

1.2 WCNS7-Z格式構造

WCNS是在CNS基礎上構造的,在解決CNS三對角線反演效率低于顯式格式的問題同時,使用類似WENO格式加權的思想,充分利用計算模板的信息,對光滑模板和包含間斷的模板分配對應的權值系數(shù),從而使格式更高效.WCNS的構造方法為3個步驟:① 特征變量節(jié)點到半節(jié)點的非線性加權插值;② 半節(jié)點的通量求解;③ 節(jié)點-半節(jié)點到節(jié)點的混合差分求導.重點介紹第1步驟中的非線性加權插值過程.

與5階格式的模板不同,在7階WCNS[19]的構造中使用4個4點子模板:(xi-3,xi-2,xi-1,xi),(xi-2,xi-1,xi,xi+1),(xi,xi+1,xi+2,xi+3),(xi,xi+1,xi+2,xi+3),如圖1所示.圖中:L為全局模板;K1、K2、K3、K4分別是4個子模板.每個子模板上的4階緊致格式設定為特定偏置形式:

圖1 WCNS7-Z格式的模板構造Fig.1 Candidate stencils of WCNS7-Z scheme

(5)

設定子模板的線性權值為

(6)

組合可得到7階緊致插值格式:

(7)

經過非線性權w0、w1、w2、w3再組合后可得:

(8)

對于非線性權wk的使用策略,Jiang等[5]設計的經典權函數(shù)方案對去除跨越間斷模板帶來的振蕩現(xiàn)象有很好的效果,但也引入了較大的非線性誤差,且在極值點附近格式的精度不能得到保證,Borges等[9]引入更高階的光滑指標,并因此提出了計算成本更低、精度更高的Z型權函數(shù),顯著地降低了由非線性權帶來的格式耗散.權重方案采用形式簡潔的Z型權函數(shù):

(9)

式中:k=0,1,2,3;ε=10-12為防止分母出現(xiàn)零除引入的小量;τ7為全局模板光滑度指標,

τ7=|β0+3β1-3β2-β3|=O(h7)

(10)

式中:O(h7)是數(shù)學上表示“與h7同階的項”的一種記法,并且可明顯地表示出截斷誤差的量級.

非線性權重中局部平滑度指標βk的大小決定了每一個子模板的光滑度,通式為

(11)

對于k=0的子模板:

(12)

對于k=1的子模板:

(13)

對于k=2的子模板:

(14)

對于k=3的子模板:

(15)

1.3 WCNS7-T格式構造

不同于傳統(tǒng)的權函數(shù)方案策略,TENO格式的模板加權策略在優(yōu)化WENO格式模板選擇的同時,更為準確地評估了模板光滑度,起到了間斷檢測器的作用,并在重構中以完全去除振蕩模板來取代WENO格式通過非線性加權策略減小振蕩模板影響的方式,使格式具有更加健壯的間斷捕捉能力.在高階的TENO方案中,如Hiejima[23]所采取的策略一樣,使用了遞增寬度的模板序列.圖2展示了WCNS7-T格式的模板構造,為了獲得7階空間精度,使用了3個3點模板和2個4點模板,從而使模板的選擇更加緊致.

圖2 WCNS7-T格式的模板構造Fig.2 Candidate stencils of WCNS7-T scheme

傳統(tǒng)WCNS中插值過程計算量大且復雜,對光滑區(qū)域來說這種插值沒有必要.利用評估策略發(fā)現(xiàn)間斷點,對于格式計算效率和精度都有非常重要的意義.將TENO格式的模板加權和寬度遞增思想引進WCNS構造,在構造間斷檢測器的同時,更快地實現(xiàn)線性格式和非線性格式之間的重建切換.因此,在WCNS7-T格式中將表征插值模板光滑度的權重分配方案更替為TENO的光滑度評測策略,有利于增強格式的激波捕捉能力.優(yōu)化后WCNS7-T格式的權重方案為

(16)

式中:k=0, 1, …,K-3;γk、χk、δk、q皆為光滑因子指示參數(shù),經截斷函數(shù)CT過濾得到最終非線性權值wk.

光滑因子計算式與WCNS的不同,具體形式為

對于非光滑的模板,根據(jù)類ENO模板的選擇,δk=0,這樣振蕩模板將會被剔除;而對于光滑模板,δk=1,格式將恢復為線性方案.顯然,截斷參數(shù)CT的取值至關重要,而控制耗散特性的CT取值依賴于不同的問題,且在實際應用中,參數(shù)調整通常非常耗時.Fu等[7]給出了CT的自適應函數(shù)取值公式.對于本文中的算例,控制耗散特性的CT值在幾組實驗下對比測算后,選定為10-7,其中K階精度τK[7]對應各個階數(shù)的計算量過大,尤其在高階時計算成本過高,本文將其值依舊限定為1,ε=10-12,其余光滑度指標參數(shù)與TENO給定的值保持一致.

2 數(shù)值測試算例

為了驗證WCNS7-T格式計算表現(xiàn)和性能,開展了一維和二維問題的算例驗證,計算結果與傳統(tǒng)的WCNS7-Z格式進行對比.求解過程中采用了HLLC通量分裂方法[24],時間推進采用3階TVD Runge-Kutta格式[25],實驗的參考精確解由WCNS5-JS通過網(wǎng)格加密獲得.

2.1 數(shù)值精度測試

為驗證WCNS7-T格式的計算性能,首先以文獻[8]中的數(shù)值精度實驗為算例,給定初始條件為

(17)

計算域[x,y]同樣設定為[0, π]×[0, π],u、v分別為x、y方向上的速度分量.各個方向上使用周期邊界條件,精確的密度解為ρ(x,y,t)=1+0.2sin(x+y-2t),計算時間直到t=2 s,表1中給出了WCNS7-Z和WCNS7-T格式L1和L∞的誤差(E1和E2)和精度(P1和P2)比較,網(wǎng)格單元劃分Nx從102到3202.可知,相較于WCNS7-Z格式,WCNS7-T格式體現(xiàn)出更好的網(wǎng)格收斂性,且格式的設計精度也得到了較好的數(shù)值驗證.

表1 二維歐拉方程的WCNS7-Z 和 WCNS7-T 格式的誤差和精度Tab.1 Errors and accuracy with two dimensional Euler equations for WCNS7-Z and WCNS7-T schemes

2.2 Sod問題

一維Sod激波管問題是檢驗差分格式性質的經典算例之一.該問題包含一道激波、一個接觸間斷和一個膨脹波,問題的具體描述參考文獻[26].給定初始條件為

計算域為[-0.5, 0.5],網(wǎng)格均勻劃分為101個網(wǎng)格點,其中Δx=0.01.在邊界處設置零梯度邊界條件,以收斂條件判斷數(shù)CFL為0.4和定時間步長Δt=0.01 s,計算到t=0.2 s.計算結果如圖3所示,WCNS7-Z格式和WCNS7-T格式都能很好地捕捉間斷,不會產生虛假的數(shù)值振蕩,表明在間斷處兩種格式都引入了足夠且穩(wěn)定的耗散.然而,WCNS7-T格式要比WCNS7-Z格式在激波和間斷處斜率更陡峭,對激波間斷位置的分辨率更準確.計算結果表明,WCNS7-T格式引入的數(shù)值耗散相對較小,對激波捕獲的精度更高.

圖3 Sod激波管問題的的密度分布結果Fig.3 Density distribution results for the Sod problem

2.3 Lax問題

Lax問題是左稀疏波、右激波類型的黎曼問題,同樣能夠用于檢驗數(shù)值算法捕捉間斷能力.設定初始條件為

(ρ,u,p)=

計算域為[0, 1],均勻劃分為101個網(wǎng)格點,其中Δx= 0.01,CFL為0.4,Δt=0.01 s,計算到t=0.14 s,密度分布計算結果如圖4所示.由放大視圖可見,數(shù)值解未觀察到過沖或振蕩問題,與精確解較為吻合.WCNS7-T通過加權組合多個低階緊致格式達到了高階緊致格式,在間斷附近采用TENO格式完全剔除振蕩模板的方法,獲得了比WCNS7-Z更加陡峭的間斷捕獲,表明WCNS7-T格式具有更弱的局部耗散特性使其能夠更高效地捕捉間斷.

圖4 Lax問題的密度分布結果Fig.4 Density distribution results for the Lax problem

2.4 Shu-Osher問題

Shu-Osher問題涉及向右移動的超音速激波與正弦分布密度場相互作用, 構造一個既有平滑結構又有不連續(xù)結構的流場,用于測試格式捕獲流場間斷和解析平滑多尺度結構的能力.設定初始條件為

(ρ,u,p)=

接觸間斷由兩個強激波的碰撞產生,以 4 001 網(wǎng)格點計算結果為參考,終止時間為t=1.8 s,兩種格式計算網(wǎng)格取為400,CFL數(shù)設為0.5,測試結果如圖5所示.從WCNS7-Z和WCNS7-T格式的計算對比結果可見,WCNS7-T格式的分辨率更高,對高頻波動的求解有著更好的波形和幅度,對可分辨頻域的波峰波谷區(qū)域抹平也更少,表明WCNS7-T格式比經典的WCNS7-Z格式具有更優(yōu)的間斷捕捉能力和脈動分辨率.

圖5 Shu-Osher問題的密度分布結果Fig.5 Density distribution results for the Shu-Osher problem

2.5 Titarev-Toro問題

Titarev-Toro問題是對Shu-Osher問題的強擴展,初始條件給定為

(ρ,u,p)=

計算域為[-5, 5],均分為1 001個網(wǎng)格點,其中Δx=0.01,邊界設置為零梯度邊界條件.給定時間步長Δt=0.5 s,CFL數(shù)設為0.5,終止時間為t=5 s.由于沒有理論上的精確解,所以參考結果由WCNS5-JS格式在 5 001 個網(wǎng)格點計算獲得.圖6給出了兩種格式計算的密度分布, 與Shu-Osher問題測試算例相似,WCNS7-T格式在激波后解析出了更好的密度波形,而WCNS7-Z格式結果中小尺度波結構較為模糊,數(shù)值耗散也較強,分辨率也較低,表明WCNS7-T格式比WCNS7-Z格式能更優(yōu)、更好地解析局部大振幅密度波動問題.

圖6 Titarev-Toro問題的密度分布結果Fig.6 Density distribution results for the Titarev-Toro problem

2.6 二維黎曼問題

Lax等[27]使用的19個二維黎曼問題中的第3個問題和第6個問題所產生的流場都具有豐富的小尺度特征,非常適合用于測驗數(shù)值格式求解精細流場結構的能力,因此選擇其中的這兩個問題來評估目前WCNS7-T格式的性能.

2.6.1算例3初始條件給定為

(ρ,u,v,p)=

設定計算域為[0, 1]×[0, 1],網(wǎng)格尺寸劃分為 1 001×1 001,終止時間為t=0.3 s,邊界條件都設定為無反射邊界條件.Kelvin-Helmholtz不穩(wěn)定性在滑移線上產生的小尺度復雜結構情況可用來評估數(shù)值格式的耗散表現(xiàn).圖7所示為黎曼問題3的密度輪廓結果.從圖7可以看出,對比于WCNS7-Z格式,WCNS7-T格式在滑移線上捕獲的渦量更多,能夠更好地捕捉到激波間斷導致的流動不穩(wěn)定性.此外,在這種數(shù)值情況下,WCNS7-T的結果與文獻[28-30]中一樣在引用TENO非線性加權后出現(xiàn)了對稱性不足的現(xiàn)象,本文推測可能是控制格式耗散表現(xiàn)的CT值的取值問題所致.

圖7 黎曼問題3密度輪廓結果Fig.7 Density contour results for Riemann problem 3

2.6.2算例6第2個黎曼的計算域、網(wǎng)格劃分和邊界條件不變,終止時間設定為t=0.3 s.初始條件為

圖8所示為黎曼問題6的密度輪廓結果.由圖可見,兩種格式都解析出了接觸線上豐富的小尺度流動結構,而WCNS7-T格式的結果有更豐富的小尺度結構,流場分辨率更高,產生的耗散也小于傳統(tǒng)的WCNS7-Z格式.表明WCNS7-T格式經過TENO非線性加權后數(shù)值耗散相比之前有所減小,且模板選擇策略的改變也進一步降低了整體的數(shù)值耗散,使格式在流場計算時能更準確模擬流場細節(jié)結構.

圖8 黎曼問題6密度輪廓結果Fig.8 Density contour results for Riemann problem 6

3 結論

基于TENO格式的模板選擇和加權策略,將其間斷檢測和光滑度因子引入WCNS構造中,發(fā)展了一種7階精度的WCNS7-T格式,在抑制數(shù)值耗散的同時增強了激波捕捉能力,通過一維和二維問題開展算例測試,驗證了WCNS7-T格式的良好性能,重要結論總結如下:

(1) WCNS7-T格式體現(xiàn)出了流場光滑區(qū)域耗散小、激波間斷區(qū)域分辨率高、計算結果更準確的特點,表明經權重優(yōu)化后格式不僅保持了原WCNS的穩(wěn)定,也強化了小擾動分辨率和激波捕獲能力.

(2) WCNS7-T格式以完全剔除振蕩模板來取代傳統(tǒng)的加權分配,而在光滑區(qū)域恢復最優(yōu)非線性插值,不僅降低了格式在間斷處需要更復雜的平滑指標所帶來的較高計算量,保證了格式的數(shù)值精度,也提高了格式的間斷檢測能力,增強了加權方法抑制數(shù)值振蕩的能力.