高精度WCNS格式加權(quán)策略改進(jìn)與數(shù)值驗(yàn)證

楊 強(qiáng), 李偉鵬

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

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

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

1 數(shù)值方法

1.1 控制方程

針對(duì)黎曼問(wèn)題,差分格式的推導(dǎo)僅考慮一維守恒型Euler方程:

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

(1)

守恒變量Q和E為

(2)

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

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

(3)

(4)

1.2 WCNS7-Z格式構(gòu)造

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

與5階格式的模板不同,在7階WCNS[19]的構(gòu)造中使用4個(gè)4點(diǎn)子模板:(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個(gè)子模板.每個(gè)子模板上的4階緊致格式設(shè)定為特定偏置形式:

圖1 WCNS7-Z格式的模板構(gòu)造Fig.1 Candidate stencils of WCNS7-Z scheme

(5)

設(shè)定子模板的線性權(quán)值為

(6)

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

(7)

經(jīng)過(guò)非線性權(quán)w0、w1、w2、w3再組合后可得:

(8)

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

(9)

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

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

(10)

式中:O(h7)是數(shù)學(xué)上表示“與h7同階的項(xiàng)”的一種記法,并且可明顯地表示出截?cái)嗾`差的量級(jí).

非線性權(quán)重中局部平滑度指標(biāo)βk的大小決定了每一個(gè)子模板的光滑度,通式為

(11)

對(duì)于k=0的子模板:

(12)

對(duì)于k=1的子模板:

(13)

對(duì)于k=2的子模板:

(14)

對(duì)于k=3的子模板:

(15)

1.3 WCNS7-T格式構(gòu)造

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

圖2 WCNS7-T格式的模板構(gòu)造Fig.2 Candidate stencils of WCNS7-T scheme

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

(16)

式中:k=0, 1, …,K-3;γk、χk、δk、q皆為光滑因子指示參數(shù),經(jīng)截?cái)嗪瘮?shù)CT過(guò)濾得到最終非線性權(quán)值wk.

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

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

2 數(shù)值測(cè)試算例

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

2.1 數(shù)值精度測(cè)試

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

(17)

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

表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問(wèn)題

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

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

圖3 Sod激波管問(wèn)題的的密度分布結(jié)果Fig.3 Density distribution results for the Sod problem

2.3 Lax問(wèn)題

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

(ρ,u,p)=

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

圖4 Lax問(wèn)題的密度分布結(jié)果Fig.4 Density distribution results for the Lax problem

2.4 Shu-Osher問(wèn)題

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

(ρ,u,p)=

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

圖5 Shu-Osher問(wèn)題的密度分布結(jié)果Fig.5 Density distribution results for the Shu-Osher problem

2.5 Titarev-Toro問(wèn)題

Titarev-Toro問(wèn)題是對(duì)Shu-Osher問(wèn)題的強(qiáng)擴(kuò)展,初始條件給定為

(ρ,u,p)=

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

圖6 Titarev-Toro問(wèn)題的密度分布結(jié)果Fig.6 Density distribution results for the Titarev-Toro problem

2.6 二維黎曼問(wèn)題

Lax等[27]使用的19個(gè)二維黎曼問(wèn)題中的第3個(gè)問(wèn)題和第6個(gè)問(wèn)題所產(chǎn)生的流場(chǎng)都具有豐富的小尺度特征,非常適合用于測(cè)驗(yàn)數(shù)值格式求解精細(xì)流場(chǎng)結(jié)構(gòu)的能力,因此選擇其中的這兩個(gè)問(wèn)題來(lái)評(píng)估目前WCNS7-T格式的性能.

2.6.1算例3初始條件給定為

(ρ,u,v,p)=

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

圖7 黎曼問(wèn)題3密度輪廓結(jié)果Fig.7 Density contour results for Riemann problem 3

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

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

圖8 黎曼問(wèn)題6密度輪廓結(jié)果Fig.8 Density contour results for Riemann problem 6

3 結(jié)論

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

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

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