尺度無關(guān)高階WCNS 格式1)

張子軒 董義道 黃梓全 孔令發(fā) 劉偉

(國防科技大學(xué)空天科學(xué)學(xué)院,長(zhǎng)沙 410073)

引言

計(jì)算流體力學(xué)在航空航天等諸多領(lǐng)域逐漸發(fā)揮越來越重要的作用[1-2],為更加真實(shí)準(zhǔn)確地刻畫復(fù)雜的流動(dòng)細(xì)節(jié),實(shí)現(xiàn)流動(dòng)的高保真模擬,高精度數(shù)值方法應(yīng)運(yùn)而生.相比于高精度有限體積格式,在結(jié)構(gòu)網(wǎng)格上,滿足自由流保持特性的高精度有限差分格式的優(yōu)勢(shì)更加明顯[3].光滑流場(chǎng)的數(shù)值模擬,通常使用線性格式,如緊致格式[4]等.然而,當(dāng)計(jì)算包含激波時(shí),需采用激波裝配法或激波捕捉法.由于激波裝配法目前還僅適用于簡(jiǎn)單構(gòu)型,難以推廣到包含復(fù)雜激波結(jié)構(gòu)的流場(chǎng)的計(jì)算[5],故目前主流的求解激波的方法仍然是激波捕捉法.

Harten[6]提出的總變差不增(total variation diminishing,TVD)格式是高分辨率的激波捕捉格式之一.隨后Van Leer 基于TVD 格式構(gòu)造了二階MUSCL(monotonic upstream schemes for conservation laws)格式,以減少數(shù)值耗散色散誤差[7-8].然而,TVD 格式在計(jì)算復(fù)雜流動(dòng)問題,如湍流時(shí),極值點(diǎn)附近存在降階現(xiàn)象.為在不降階的前提下提高精度,高階格式得以發(fā)展.高階格式可以在一定程度上減少耗散和色散誤差[9],并且在相同精度的條件下較低階格式更為高效[4].此外,相較于低階格式,高階格式具有更好的并行可擴(kuò)展性[10].

作為對(duì)高階基本無振蕩ENO(essentially nonoscillatory)格式[11-12]的改進(jìn),Liu 等[13]提出的加權(quán)基本無震蕩WENO 格式是目前較為理想的高階激波捕捉格式之一.WENO 格式較好地兼顧了間斷區(qū)域的數(shù)值震蕩和光滑區(qū)域的計(jì)算精度,被廣泛應(yīng)用于各類流動(dòng)的數(shù)值模擬[14–16].隨后,Jiang 等[17]引入了新的光滑度量因子,使WENO 格式達(dá)到最優(yōu)精度,該格式被命名為WENO-JS.針對(duì)經(jīng)典的WENO 格式所存在的極值點(diǎn)附近降階并且耗散過大等問題,典型的改進(jìn)方法主要包括: Henrick 等[18]提出映射權(quán)函數(shù)以擴(kuò)大線性權(quán)的范圍,使得較大梯度下也能恢復(fù)線性權(quán),克服了極值點(diǎn)附近降階的問題但計(jì)算花費(fèi)較大;Borges 等[19]通過引入高階全局光滑度量因子構(gòu)造的形式更為簡(jiǎn)單的Z 型權(quán).Z 型權(quán)可以以較小的計(jì)算量獲得與原始JS 權(quán)基本一致的精度且不存在降階問題.此外,一些其他常見的WENO 改進(jìn)格式也被陸續(xù)提出,如CWENO[20],WENO-Z+[21]以及HWENO[22]等.

另一種典型的高階激波捕捉格式是由Deng等[23-24]結(jié)合緊致格式[4]和非線性加權(quán)技術(shù),在緊致非線性格式CNS (compact nonlinear scheme)[25]基礎(chǔ)上構(gòu)造的加權(quán)緊致非線性格式WCNS.WCNS 格式已經(jīng)成功應(yīng)用于各種類型的數(shù)值模擬,包括邊界層轉(zhuǎn)捩及湍流[26]、氣動(dòng)聲學(xué)[27]、多組分流[28]以及激波-邊界層干擾[29]等.總的來說,基于變量插值的WCNS 格式相較于與Jiang 等[17]提出的基于通量函數(shù)重構(gòu)的WENO 格式相比有如下優(yōu)點(diǎn)[30-33]: (1)相同精度下分辨率更高;(2)數(shù)值通量選擇更加靈活[34];(3)自由流和渦流保持特性上表現(xiàn)更好[35].

然而,經(jīng)典的WCNS/WENO 格式仍存在一些明顯不足,如靈敏度參數(shù)和尺度因子的取值對(duì)格式的數(shù)值性能有很大影響.例如,由Deng 等[23-24,36]的研究可知,使用WCNS 格式離散Navier-Stokes 方程時(shí),非線性加權(quán)函數(shù)中靈敏度參數(shù)被假定為 ε=1.0×10?6.但使用該格式離散雷諾應(yīng)力方程時(shí)[37],由于雷諾應(yīng)力的量級(jí)很小約為 1.0×10?6,ε=1.0×10?18才能保證計(jì)算穩(wěn)定進(jìn)行,而 ε=1.0×10?6或 ε=1.0×10?12將會(huì)導(dǎo)致計(jì)算發(fā)散.此外,當(dāng)使用靈敏度參數(shù)ε=1.0×10?12的W E N O 格式重 構(gòu)小尺 度因子λ=1.0×10?7下的Lax 問題時(shí)(Lax 問題的初始條件為Q(x,t=0)=λQ0(x)且Q0(x)=(ρ,ρu,E)),密 度ρ將會(huì)在間斷附近出現(xiàn)明顯的數(shù)值振蕩[38].因?yàn)槭褂? 階WENO 格式對(duì)小尺度下的不連續(xù)函數(shù)進(jìn)行重構(gòu)的過程中,相對(duì)較大的靈敏度參數(shù)對(duì)相對(duì)較小的光滑度量因子占主導(dǎo)地位,使得子模板喪失了自適應(yīng)性.

為解決上述問題,Don 等[38]引入去尺度函數(shù)構(gòu)建了尺度無關(guān)的WENO 格式.受Don 等[38]工作啟發(fā),基于WCNS/WENO 格式存在的問題,本文采用函數(shù)的平均值構(gòu)造去尺度函數(shù).通過在WCNS 格式插值過程中引入去尺度函數(shù)修正非線性權(quán)重中的靈敏度參數(shù),衍生出一種更有效的非線性權(quán)重,從而提出了一種不受靈敏度參數(shù)和尺度因子控制的尺度無關(guān)的WCNS 格式.研究?jī)?nèi)容主要包括兩方面: 首先,探討WCNS 格式與靈敏度參數(shù)和尺度因子的關(guān)系;此外,將去尺度函數(shù)引入經(jīng)典的WCNS7-JS/Z/D 格式的非線性權(quán)重中以消除尺度的依賴性,構(gòu)造尺度無關(guān)的WCNS 格式(WCNS7-JSm/Zm/Dm).文章的結(jié)構(gòu)如下: 第1 節(jié)主要介紹了WCNS 格式離散方法,推導(dǎo)了7 階D 權(quán)函數(shù)的形式并且給出了WCNS7-JSm/Zm/Dm 格式插值算法的步驟;第2 節(jié)在雙精度算法下驗(yàn)證了格式的基本無振蕩性質(zhì),并在4 精度算法下驗(yàn)證了格式的尺度無關(guān)(scale-invariant,Si)性質(zhì)和極值點(diǎn)(critical point,Cp)性質(zhì);在文章的第3 節(jié)進(jìn)行了精度測(cè)試,并通過幾個(gè)一維和二維算例測(cè)試了新格式在激波捕捉以及分辨小尺度渦結(jié)構(gòu)等方面的性能;第4 節(jié)給出了本文結(jié)論.

1 控制方程及其離散

1.1 WCNS 差分格式

考慮如下一維雙曲守恒律

其中u(x,t) 是守恒變量,f(u)是通量.考慮在均勻網(wǎng)格上離散,空間坐標(biāo)xj=jh(j=0,1,2,···,N),網(wǎng)格間距h=(b?a)/N,式(1)的半離散形式為

早期的WCNS 格式[39]常采用Lele[4]提出的單元中心型緊致格式計(jì)算一階導(dǎo)數(shù)

式中 κ,a,b,c,d表示系數(shù).如果 κ ≠0,式(3)將需要進(jìn)行三對(duì)角矩陣求逆,此時(shí)被稱為隱式格式.已有研究[24,40]發(fā)現(xiàn),差分格式的形式,即顯式或隱式,對(duì)計(jì)算結(jié)果影響不大,故這里選擇更加簡(jiǎn)單高效的8 階顯式中心差分格式(κ=0),形式如下

1.2 WCNS 插值格式

WCNS 采用非線性加權(quán)技術(shù),對(duì)上述4 個(gè)子模板 (r=4)插值進(jìn)行加權(quán)組合

其中,ωk為非線性權(quán).當(dāng) ωk=dk時(shí),上述非線性插值恢復(fù)成線性插值格式.dk稱為線性權(quán)或最優(yōu)權(quán)

1.3 非線性權(quán)函數(shù)

非線性加權(quán)技術(shù)是在流場(chǎng)光滑區(qū)令非線性權(quán)逼近線性權(quán),以保證設(shè)計(jì)精度,即 ωk=dk;在流場(chǎng)間斷區(qū)減小不光滑子模板的權(quán)重,以避免跨間斷插值引起數(shù)值震蕩和計(jì)算不收斂甚至發(fā)散等情況,此時(shí) ωk為非線性權(quán),可以通過不同的非線性權(quán)函數(shù)求解得到.本文使用到的非線性權(quán)函數(shù)包括經(jīng)典的JS 權(quán)函數(shù)[17]

耗散更低的Z 權(quán)函數(shù)[19]

以及廣義上與靈敏度參數(shù)無關(guān)的D 權(quán)函數(shù)[44]

權(quán)函數(shù)中指數(shù)p和q通過改變子模板光滑度量因子的偏離程度來影響數(shù)值耗散特性,本文默認(rèn)取p=2q=2. ε 稱為靈敏度參數(shù),是避免分母為零的小量,JS,Z 和D 權(quán)函數(shù)中 ε 通常取1.0×10?6,1.0×10?40和1.0×10?40.βk為子模板的光滑度量,具體形式為各子模板在xj點(diǎn)處1～3 階導(dǎo)數(shù)的平方和

式(11)包含一個(gè)線性項(xiàng)L和一個(gè)非線性項(xiàng) Γk

τ

為全模板的高階全局光滑度量,7 階WCNS 格式對(duì)應(yīng)的 τ7=|β0+3β1?3β2?β3|[45]在xj點(diǎn)處的泰勒展開式為

式(12)中 Φ 為修正函數(shù),相當(dāng)于激波探測(cè)器: 在流場(chǎng)光滑區(qū),?很小且 Φ=?;在間斷區(qū),?很大且 Φ=1,此時(shí),非線性權(quán)函數(shù)恢復(fù)成Z 權(quán)函數(shù).根據(jù)文獻(xiàn)[44]和式(15)可推導(dǎo)修正函數(shù) Φ 的表達(dá)式為

1.4 WCNS7-JSm/Zm/Dm 格式

在JS/Z/D 權(quán)函數(shù)基礎(chǔ)上引入去尺度函數(shù),探索了一種更加簡(jiǎn)單穩(wěn)定的非線性權(quán)函數(shù)的構(gòu)造方法,并將其應(yīng)用于7 階WCNS 格式,構(gòu)造尺度無關(guān)的WCNS7-JSm/Zm/Dm 格式,具體步驟如下.

(1)引入去尺度函數(shù) ξ和 μ.

去尺度函數(shù) ξ定義為每個(gè)網(wǎng)格點(diǎn)上函數(shù)值{fj,,j=0,1,2,···,N}的絕對(duì)值和的平均,即

去尺度函數(shù) μ定義為

式中,μ0是避免f(x)=0 且 ξ=0時(shí)導(dǎo)致分母為零的情況,在雙精 度和四 精度中 μ0分別取1.0×10?40和1.0×10?100.

(2)使用去尺度函數(shù) μ進(jìn)行修正.

修正Z/D 權(quán)函數(shù)非線性項(xiàng)中的靈敏度參數(shù)

修正D 權(quán)函數(shù)中的修正函數(shù) Φ

(3)得到修正后的格式.

基于修正后的JS 權(quán)構(gòu)造的WCNS7-JSm 格式為

基于修正后的Z/D 權(quán)函數(shù)構(gòu)造的WCNS7-Zm/Dm格式為

WCNS7-JSm/Zm/Dm 格式的非線性權(quán)為

(4)根據(jù)上述步驟求得非線性權(quán) ωk,插值后可得本文采用特征變量插值,半節(jié)點(diǎn)上對(duì)流通量的計(jì)算使用Rusanov 格式[43],時(shí)間格式采用3 階Runge-Kutta,庫朗數(shù)為CFL=0.4.

2 性質(zhì)驗(yàn)證

本節(jié)驗(yàn)證了不同靈敏度參數(shù) ε和尺度因子 λ下WCNS 格式的3 種(ENO-,Si-和Cp-)性質(zhì).采用雙精度算法驗(yàn)證格式的ENO 性質(zhì);采用四精度算法分析驗(yàn)證Si 性質(zhì)和Cp 性質(zhì).

2.1 ENO 性質(zhì)

參照文獻(xiàn)[44],這里使用WCNS7-JS/Z/D 格式在 ε=1.0×10?6和ε=1.0×10?40時(shí)分別模擬小尺度λ=1.0×10?7和正常尺度λ=1下的Lax 問題,以驗(yàn)證WCNS7-JS/Z/D格式在不同靈敏度參數(shù)和尺度因子下的ENO 性質(zhì).

Lax 問題的初始條件為Q(x,t=0)=λQ0(x),Q0(x)=(ρ,ρu,E)且E=p/γ ?1 +ρu2/2,即[46]

計(jì)算終止時(shí)間為t=1.3,網(wǎng)格數(shù)為N=200.

圖1和圖2 給出了 在 ε=1.0×10?6和ε=1.0×10?40時(shí)分別使用WCNS 格式求解Lax 問題的計(jì)算結(jié)果,其中紅色實(shí)線和藍(lán)色實(shí)線分別為 λ=1,λ=1.0×10?7時(shí)的計(jì)算結(jié)果,黑色虛線代表精確解.ε=1.0×10?6時(shí),3 種格式的密度分布曲線基本相同,因此圖1 僅給出WCNS7-JS 格式的密度分布曲線;ε=1.0×10?40時(shí),WCNS7-JS 與WCNS7-Z 格式密度分布曲線基本相同但與WCNS7-D 格式不同,故圖2 給出WCNS7-JS 和WCNS7-D 格式的計(jì)算結(jié)果.

圖1 ε=1.0×10?6 時(shí)Lax 問題的密度分布曲線Fig.1 Density distribution curves of the Lax problem with ε=1.0×10?6

圖2 ε=1.0×10?40 時(shí)Lax 問題的密度分布曲線Fig.2 Density distribution curves of the Lax problem with ε=1.0×10?40

從圖1 和圖2 中可以看出,當(dāng) λ=1時(shí),無論靈敏度參數(shù) ε如何變化,密度ρ都滿足ENO 性 質(zhì).當(dāng)λ=1.0×10?7時(shí),ε=1.0×10?6下的WCNS7-JS/Z 格式在不連續(xù)區(qū)域附近會(huì)出現(xiàn)明顯的數(shù)值振蕩;WCNS7-D格式則是不管靈敏度參數(shù)如何取值,在不連續(xù)區(qū)域附近都會(huì)出現(xiàn)數(shù)值振蕩.因此,WCNS7-D 格式在小尺度下對(duì)任意靈敏度參數(shù)都不滿足ENO 性質(zhì),而WCNS7-JS/Z 格式僅在小尺度和較大的靈敏度參數(shù)下不滿足ENO 性質(zhì).綜上,靈敏度參數(shù)和尺度因子的取值對(duì)WCNS 格式的ENO 性質(zhì)有顯著的影響.

圖3 展示的是 ε=1.0×10?6和ε=1.0×10?40時(shí)的WCNS7-Dm 格式(WCNS7-JSm/Zm/Dm 格式的密度分布曲線基本相同) 求解 λ=1 和λ=1.0×10?7下Lax 問題的密度分布圖.從圖3 反映的結(jié)果可以看出,不管靈敏度參數(shù)和尺度因子如何變化,WCNS7-JSm/Zm/Dm 格式都滿足ENO 性質(zhì),即某種程度上可以認(rèn)為修正后的格式與靈敏度參數(shù)和尺度因子的取值無關(guān).

圖3 ε=1.0×10?6和 ε=1.0×10?40 時(shí)使用WCNS7-Dm 格式求解Lax 問題的密度分布曲線Fig.3 Density distribution curves of the Lax problem computed by the WCNS7-Dm scheme with ε=1.0×10?6 and ε=1.0×10?40

2.2 Si 性質(zhì)

接下來驗(yàn)證WCNS 格式的Si 性質(zhì).在繼續(xù)進(jìn)行數(shù)值模擬之前,需要給出Si 性質(zhì)和弱Si 性質(zhì)的定義.

定義1:如果WCNS 格式的非線性插值算子W[·] 對(duì)任意尺度因子 λ ∈R和任何給定函數(shù)f都滿足W[λf]=λW[f],稱此WCNS 格式是尺度無關(guān)的,即滿足Si 性質(zhì).

然而,Si 性質(zhì)要求的條件非常嚴(yán)格,故引入臨界尺度因子 λ?和弱Si 性質(zhì).

定義2:如果存在臨界尺度因子 λ?≥0 且 λ?∈R使得WCNS 算子W[·] 對(duì)任何給定函數(shù)f和任意兩個(gè)尺度因子 λ1,λ2且λ1,λ2∈S∞={λ|λ ∈R,|λ|≥λ?}都滿足 λ2W[λ1f]=λ1W[λ2f],則稱此WCNS 格式是弱尺度無關(guān)的,即滿足弱Si 性質(zhì).

計(jì)算L∞范數(shù)下Si 的誤差ESi(λ)來驗(yàn)證WCNS格式的(弱)Si 性質(zhì),ESi(λ)的形式如下

式中,?0稱為機(jī)器舍入誤差,λ∞是參考尺度因子,λ ∈[1.0×10?19,1.0×1019].因?yàn)閰⒖汲叨纫蜃?λ∞=1.0×1020足夠大,故參考WCNS 算子W[λ∞f]/λ∞可以假設(shè)其計(jì)算是與尺度無關(guān)的.

為驗(yàn)證WCNS 格式的Si 性質(zhì),本節(jié)計(jì)算了求解雙曲守恒律中經(jīng)常遇到的一些典型函數(shù)[38]: (1)C∞正弦函數(shù);(2)C0分段光滑ReLU;(3) Heaviside 函數(shù),函數(shù)的定義如下

如圖4～圖6 給出了6 種WCNS 格式分別求解正弦函數(shù)、ReLU 函數(shù)和Heaviside 函數(shù)在尺度因子λ ∈[1.0×10?19,1.0×1019] 范圍內(nèi)的Si 誤差ESi(λ),其中 ε=1.0×10?6(上行)、ε=1.0×10?40(下行),并且有3 套不同的網(wǎng)格數(shù).從圖中可以看出,對(duì)于任意 λ誤差ESi(λ)通常很小,在雙精度算法中,計(jì)算誤差與機(jī)器舍入誤差 ?0≈1.0×10?16將無法區(qū)分,故需要在四精度算法中進(jìn)行,機(jī)器的舍入誤差為 ?0≈1.0×10?34.

圖4 3 套網(wǎng)格單元數(shù) N 下使用WCNS 格式計(jì)算Sine 函數(shù) f1(x) 的Si 誤差 E Si(λ),λ ∈[1.0×10?19,1.0×1019],ε=1.0×10?6 (上行)和ε=1.0×10?40 (下行)Fig.4 The Si-error E Si(λ) of the Sine function f1(x) calculated by the WCNS schemes with ε=1.0×10?6 (top row),ε=1.0×10?40 (bottom row),and different scale factors λ ∈[1.0×10?19,1.0×1019] for three number of cells N

圖5 3 套網(wǎng)格單元數(shù) N 下使用WCNS 格式計(jì)算 R eLU 函數(shù) f2(x) 的Si 誤差 E Si(λ),λ ∈[1.0×10?19,1.0×1019],ε=1.0×10?6 (上行)和ε=1.0×10?40 (下行)Fig.5 The Si-error E Si(λ) of the R eLU function f2(x) calculated by the WCNS schemes with ε=1.0×10?6 (top row),ε=1.0×10?40 (bottom row),and different scale factors λ ∈[1.0×10?19,1.0×1019] for three number of cells N

圖6 3 套網(wǎng)格單元數(shù) N 下使用WCNS 格式計(jì)算Heaviside 函數(shù) f3(x) 的Si 誤差 E Si(λ),λ ∈[1.0×10?19,1.0×1019],ε=1.0×10?6 (上行)和ε=1.0×10?40 (下行)Fig.6 The Si-error E Si(λ) of the Heaviside function f3(x) calculated by the WCNS schemes with ε=1.0×10?6 (top row),ε=1.0×10?40 (bottom row),and different scale factors λ ∈[1.0×10?19,1.0×1019] for three number of cells N

首先觀察尺度無關(guān)的WCNS 格式,如圖4～圖6(左),對(duì)于任意的尺度因子 λ和靈敏度參數(shù) ε,測(cè)試函數(shù) {fl(x),l=1,2,3}的 誤差ESi(λ)基 本上為 1.0×10?34,這說明WCNS7-JSm/Zm/Dm 格式滿足Si 性質(zhì).其次,觀察經(jīng)典的WCNS 格式,如圖4～圖6(中和右),誤差ESi(λ)明 顯取決于 λ 并 且在較小程度上也取決于 ε,例如,對(duì)于分段光滑ReLU 函數(shù)f2(x)和Heaviside 函數(shù)f3(x),臨界尺度因子 λ?隨著靈敏度參數(shù) ε的減小而減小.具體而言,對(duì)于Heaviside 函數(shù)f3(x),當(dāng)靈敏度參數(shù)從 1.0×10?6減小到1.0×10?40時(shí),WCNS7-JS 和WCNS7-Z格式的臨界尺度因子 λ?從O(1014)減少 到O(10?3);同理,ε從1.0×10?6減 小到1.0×10?40時(shí),WCNS7-D 格式的臨界尺度因子 λ?從O(1014)減少到O(100),這說明經(jīng)典的WCNS 格式只滿足弱Si 性質(zhì).最后,結(jié)合文獻(xiàn)[38]可以得出如下結(jié)論:

(1) 對(duì)于任意給定的靈敏度參數(shù) ε,WCNS7-JS 和WCNS7-Z 格式在臨界尺度因子時(shí)滿足弱Si 性質(zhì);

(2) 對(duì)于任意給定的靈敏度參數(shù) ε,WCNS7-D格式在臨界尺度因子時(shí)滿足弱Si 性質(zhì);

(3)WCNS7-JSm/Zm/Dm 格式滿足Si 性質(zhì)且與 ε和 λ無關(guān),故被稱為尺度無關(guān)的WCNS 格式.

2.3 Cp 性質(zhì)

為討論WCNS 格式的Cp 性質(zhì),結(jié)合文獻(xiàn)[21,44]給出極值點(diǎn)、Cp 性質(zhì)的定義.

定義3:若f′(xc)=...=且ncp>0,稱f(x)在xc處有一 個(gè)ncp階的極值點(diǎn).若ncp=0,則xc不是極值點(diǎn).

定義4:對(duì)于給定指數(shù)p和靈敏度參數(shù) ε,如果一個(gè) (2r?1) 階的WCNS 格式直到第r階的極值點(diǎn)處逼近光滑函數(shù)的一階導(dǎo)數(shù)時(shí)都保持其最優(yōu)精度,稱此WCNS 格式滿足Cp 性質(zhì).

為研究新格式在不同尺度因子 λ下對(duì)存在高階極值點(diǎn)的光滑函數(shù)所能達(dá)到的收斂精度,使用如下測(cè)試函數(shù)

式中 λ >0,函數(shù)在x=0且ncp=n處有一個(gè)n階極值點(diǎn).圖7 給出了當(dāng) λ=1,ncp=1,2,3,4時(shí)WCNS 格式在 ε=1.0×10?6和 ε=1.0×10?40時(shí)的L∞誤差和收斂階,線條顏色(藍(lán)/粉/橙/紅/綠/黑)分別代表WCNS7-(JS/Z/D/JSm/Zm/Dm)格式.表1 給出了 λ=1.0×10?7和 λ=100 時(shí),ε=1.0×10?40且ncp=3 下的L∞誤差和WCNS 格式的收斂階,由于WCNS7-JS 和WCNS7-JSm 格式的誤差值和收斂階大致相同,WCNS7-Z和WCNS7-Zm 格式的誤差值和收斂階也基本相同,因此,表1 僅給出了WCNS-JSm/Zm/D/Dm 格式的結(jié)果.根據(jù)圖7、表1 以及參考文獻(xiàn)[38,44]可以得到如下結(jié)論.

圖7 在極值點(diǎn) n cp=n=1,2,3,4 下WCNS 格式的收斂精度和 L∞ 誤差,且 λ=1,ε=1.0×10?6 (上行),ε=1.0×10?40 (下行)Fig.7 L∞ error of the WCNS schemes with ε=1.0×10?6 (top row),ε=1.0×10?40 (bottom row),and the scale factor λ=1 for ncp=n=1,2,3,4

表1 WCNS 格式在 ε=1.0×10?40,n cp=3 時(shí)的 L∞ 誤差 E N 和收斂階 O(?xs)Table 1 L∞ error E N,and order of accuracy O (?xs) of WCNS schemes with ε=1.0×10?40 and ncp=3

(1)當(dāng)靈敏度參數(shù)較大時(shí) ε=1.0×10?6,6 種WCNS格式都能達(dá)到設(shè)計(jì)精度O(?x7).

(2) 當(dāng)靈敏度參數(shù)較小時(shí) ε=1.0×10?40,僅有WCNS7-D/Dm 格式能夠達(dá)到設(shè)計(jì)精度O(?x7),WCNS7-JS/JSm 和WCNS7-Z/Zm 格式只有在ncp≥3時(shí)精度才能達(dá)到O(?xncp).

(3)對(duì)于小靈敏度參數(shù) ε=1.0×10?40,WCNS7-Z/Zm 格式只有在ncp=1時(shí)能夠達(dá)到設(shè)計(jì)精度;WCNS7-JS/JSm 格式在ncp=1,2,3,4時(shí)都無法達(dá)到設(shè)計(jì)精度.

(4)對(duì)于任意靈敏度參數(shù) ε,WCNS7-D/Dm 格式都能滿足Cp 性質(zhì)并且能保持最優(yōu)精度而WCNS7-JS/JSm 和WCNS7-Z/Zm 格式則不能.

綜上,WCNS7-D/Dm 格式比WCNS7-JS/JSm和WCNS7-Z/Zm 格式表現(xiàn)出更好的精度特性.這些結(jié)論與Don 等[38]的觀察是一致的,盡管基于通量重構(gòu)的WENO 格式與基于變量插值的WCNS 格式的離散化過程存在明顯差異,但這兩種格式在Cp 性質(zhì)方面沒有明顯的差異.

3 數(shù)值實(shí)驗(yàn)

接下來的數(shù)值實(shí)驗(yàn)將進(jìn)一步研究新格式的性質(zhì)以及新格式在典型算例中的表現(xiàn).數(shù)值實(shí)驗(yàn)包括一維線性對(duì)流方程的精度測(cè)試,一維Euler 方程,包括Lax 問題[46]、Sod 問題[47]、激波密度波相互作用問題[19]以及爆炸波相互作用問題[48],以及二維Euler方程,如二維Riemann 問題[49]、雙馬赫反射[48]以及Richtmyer-Meshko 不穩(wěn)定性問題[50].所有算例均在雙精度算法下進(jìn)行并采用特征變量插值和比熱比γ=1.4的量熱完全氣體模型.

3.1 精度測(cè)試

求解3 種尺度因子(λ=1.0×10?7,1,100) 下的線性對(duì)流方程以測(cè)試 ε=1.0×10?6和ε=1.0×10?40時(shí)WCNS 格式的收斂精度,由于兩種靈敏度參數(shù)下的結(jié)果相似,故省略了 ε=1.0×10?6時(shí)的結(jié)果.測(cè)試用到高斯波函數(shù)

計(jì)算一個(gè)周期t=2,計(jì)算域x∈[?1,1],左右邊界設(shè)為周期性邊界且?guī)炖蕯?shù)為CFL=0.1.簡(jiǎn)便起見,這里僅給出WCNS7-D 和WCNS7-Dm 格式的計(jì)算結(jié)果,如表2 所示.數(shù)據(jù)結(jié)果表明,WCNS7-Dm 格式的收斂精度與經(jīng)典的WCNS7-D 格式的收斂精度基本一致,并且隨著網(wǎng)格加密都可以達(dá)到設(shè)計(jì)(7 階)精度.WCNS7-JSm 和WCNS7-Zm 格式有相同的結(jié)論.

表2 WCNS7-D/Dm 格式的誤差值和收斂階 (ε=1.0 × 10?40)Table 2 Error and accuracy statistics for WCNS7-D/Dm schemes (ε=1.0 × 10?40)

3.2 一維Euler 方程

3.2.1 Lax 和Sod 問題

Lax 問題的初始條件如式(25).Sod 問題的初始條件[47]為

計(jì)算終止時(shí)間為t=0.2,網(wǎng)格數(shù)為N=200.

圖8 給出了使用靈敏度參數(shù)為 1.0×10?6和1.0×10?40的WCNS 格式求解3 種尺度(λ=1.0×10?7,1,100)下Lax 問題的計(jì)算結(jié)果.由于Lax 和Sod 問題的結(jié)論相同,為簡(jiǎn)便起見,圖9 僅給出了 ε=1.0×10?40時(shí)使用WCNS格式計(jì)算 λ=1.0×10?7下的Sod 問題的結(jié)果.觀察圖8 和圖9,對(duì)于 λ=1和λ=102,6 種WCNS 格式的表現(xiàn)基本相同,均沒有出現(xiàn)數(shù)值振蕩并滿足ENO 性質(zhì);當(dāng) λ=1.0×10?7時(shí),WCNS7-J Sm/Zm/Dm 格式?jīng)]有出現(xiàn)數(shù)值振蕩,依舊保持ENO性質(zhì),而WCNS7-JS/Z/D 格式在不連續(xù)區(qū)域附近出現(xiàn)了明顯的數(shù)值振蕩,失去了ENO 性質(zhì),這就驗(yàn)證2.1 節(jié)的分析結(jié)果.

圖8 ε=1.0×10?6和 ε=1.0×10?40 時(shí)分別使用WCNS 格式計(jì)算 λ=1.0×10?7,1,100 下的Lax 問題Fig.8 Density of the Lax problem computed by the WCNS schemes with ε=1.0×10?6 and ε=1.0×10?40 for λ=1.0×10?7,1,100

圖9 ε=1.0×10?40 時(shí)使用WCNS 格式計(jì)算 λ=1.0×10?7 的Sod 問題Fig.9 Density of the Sod problem computed by the WCNS schemes with ε=1.0×10?40 for λ=1.0×10?7

3.2.2 激波密度波相互作用問題

該問題描述一道馬赫數(shù)為3 的右行激波與一系列密度波相互作用,可用來測(cè)試格式的激波捕捉能力、耗散性、數(shù)值穩(wěn)定性和高頻波分辨率,初始條件為[19]

其中,a=0.2,x0=?4 且k=5.計(jì)算終 止時(shí)間t=5,網(wǎng)格數(shù)N=600.

簡(jiǎn)便起見,圖10 僅給出了 λ=1.0×10?7且 ε=1.0×10?40時(shí)6 種WCNS 的計(jì)算結(jié)果.觀察圖10(b),注意到WCNS7-D 格式(綠色)曲線在區(qū)間[0.73,0.8]上存在嚴(yán)重的過沖,這可能與WCNS7-D 格式低耗散性質(zhì)有關(guān).其他5 種格式?jīng)]有明顯的過沖或下沖現(xiàn)象,這說明WCNS7-Dm 格式相較于WCNS7-D 格式實(shí)現(xiàn)了對(duì)過沖較好的抑制作用.曲線過沖或下沖可能會(huì)增大誤差,影響計(jì)算的穩(wěn)定性,故除存在過沖現(xiàn)象的WCNS7-D 格式(綠色),WCNS7-Dm 格式(紅色)相較于其他格式表現(xiàn)略好且分辨率略高.

圖10 ε=1.0×10?40 時(shí)使用WCNS 格式計(jì)算 λ=1.0×10?7 激波密度波問題Fig.10 Density of the shock-density wave interaction problem computed by the WCNS schemes with ε=1.0×10?40 for λ=1.0×10?7

3.2.3 爆炸波相互作用問題

接下來計(jì)算爆炸波相互作用問題[48],初始條件為

計(jì)算終止時(shí)間為t=0.038,網(wǎng)格數(shù)N=400.

圖11 給出 λ=1.0×10?7和 ε=1.0×10?40下的計(jì)算結(jié)果.觀察圖11(b),注意到WCNS7-D 格式(綠色)曲線在區(qū)間 [ 0.73,0.8]上存在微弱的下沖和過沖,其他5 種格式表現(xiàn)良好均不存在過沖或下沖.同樣可以看出,WCNS-Dm 格式相較于WCNS-D 格式實(shí)現(xiàn)了對(duì)過沖和下沖行為較好的抑制作用.除存在過沖和下沖現(xiàn)象的WCNS7-D 格式(綠色),WCNS7-Dm 格式(紅色)相較于其他格式表現(xiàn)略好且分辨率略高.此外,在相同條件下使用網(wǎng)格數(shù)N=800重新計(jì)算該問題,如圖12 所示.觀察圖11 和圖12,可以看出兩種網(wǎng)格分辨率下的現(xiàn)象基本相同.總的來說,WCNS7-Dm 格式性能表現(xiàn)最優(yōu).

圖11 N=400 時(shí)使用WCNS 格式計(jì)算爆炸波相互作用問題Fig.11 Density of the blast-waves interaction problem computed by the WCNS schemes with N=400

圖12 N=800 時(shí)使用WCNS 格式計(jì)算爆炸波相互作用問題Fig.12 Density of the blast-waves interaction problem computed by the WCNS schemes with N=800

3.3 二維Euler 方程

3.3.1 二維Riemann 問題

該問題包含多尺度的流場(chǎng)結(jié)構(gòu),可測(cè)試數(shù)值格式的多尺度分辨率.使用Riemann 問題[49]中的構(gòu)型3 和構(gòu)型6,這兩種構(gòu)型的流場(chǎng)中包含豐富的小尺度結(jié)構(gòu),可以更好地測(cè)試數(shù)值格式分辨小尺度結(jié)構(gòu)的能力.計(jì)算域?yàn)?[ 0,1]×[0,1],構(gòu)型3 的初始條件為

計(jì)算終止時(shí)間為t=0.8,網(wǎng)格數(shù)量為 4 00×400.

圖13 給出了 ε=1.0×10?40時(shí)WCNS7-JSm/Zm/Dm格式的計(jì)算結(jié)果,密度等值線是 [ 0.1,1.8]×λ等間距的30 條.觀察圖13,WCNS7-JSm/Zm/Dm 格式都可以基本無振蕩地捕捉不連續(xù),并且在兩種不同的尺度因子下,λ=1 (左) 和λ=1.0×10?7(右),密度等值線分布基本相同.ε=1.0×10?6時(shí)的WCNS7-JSm/Zm/Dm 格式也有類似的結(jié)果,簡(jiǎn)便起見,在此省略.結(jié)果表明,WCNS7-JSm/Zm/Dm 格式同時(shí)滿足Si 性質(zhì)和ENO 性質(zhì).構(gòu)型6 的初始條件為

圖13 ε=1.0×10?40 時(shí)使用WCNS 格式計(jì)算Riemann 問題中的構(gòu)型3Fig.13 Density of the configuration 3 of 2-D Riemann problems computed by the WCNS7-JSm/Zm/Dm schemes with ε=1.0×10?40

圖13 ε=1.0×10?40 時(shí)使用WCNS 格式計(jì)算Riemann 問題中的構(gòu)型3 (續(xù))Fig.13 Density of the configuration 3 of 2-D Riemann problems computed by the WCNS7-JSm/Zm/Dm schemes with ε=1.0×10?40 (continued)

計(jì)算終止時(shí)間為t=0.3,網(wǎng)格數(shù)量為 4 00×400.

圖14 比較了小尺度因子 λ=1.0×10?7下靈敏度參數(shù)分別為 1.0×10?6和 1.0×10?40時(shí)的WCNS7-D和WCNS-Dm 格式的計(jì)算結(jié)果,密度等值線是[0.1,2.9]×λ等間距的40 條.觀察圖14,WCNS7-D格式圖14(a) 沿滑移線附近的振蕩更加明顯,而WCNS7-Dm 格式圖14(b)則依舊能夠基本無振蕩地捕捉不連續(xù).

圖14 ε=1.0×10?40 時(shí)使用WCNS7-D/Dm 格式計(jì)算 λ=1.0×10?7 下的二維 Riemann 問題中的構(gòu)型6Fig.14 Density of the configuration 6 of 2-D Riemann problems computed by the WCNS7-D/Dm schemes with λ=1.0×10?7 and ε=1.0×10?40

圖14 ε=1.0×10?40 時(shí)使用WCNS7-D/Dm 格式計(jì)算 λ=1.0×10?7 下的二維 Riemann 問題中的構(gòu)型6 (續(xù))Fig.14 Density of the configuration 6 of 2-D Riemann problems computed by the WCNS7-D/Dm schemes with λ=1.0×10?7 and ε=1.0×10?40(continued)

3.3.2 雙馬赫反射

在雙馬赫反射問題[48]中,馬赫數(shù)為10 的斜激波以60°角入射到靜止平板上并繼續(xù)傳播,所產(chǎn)生的流場(chǎng)中包含豐富的旋渦結(jié)構(gòu)和十分復(fù)雜的激波反射,可用來測(cè)試數(shù)值格式捕捉激波和回流區(qū)小尺度結(jié)構(gòu)等能力.計(jì)算域?yàn)?[ 0,4]×[0,1],初始條件為

計(jì)算終止時(shí)間為t=0.2,網(wǎng)格數(shù)量為 9 61×241.

圖15 給出了 ε=1.0×10?40時(shí)的WCNS7-Dm 格式的計(jì)算結(jié)果,密度等值線是[2,23]× λ等間距的39 條.觀察圖15,WCNS7-Dm 格式在 λ=1.0×10?7和 λ=1的情況下都可以基本無振蕩地捕捉不連續(xù),并且兩種尺度因子下的密度等值線分布基本相同.WCNS7-JSm 和WCNS7-Zm 格式均有類似結(jié)果,簡(jiǎn)便起見,在此省略.這說明尺度無關(guān)的WCNS7-JSm/Zm/Dm 格式同時(shí)滿足Si 性質(zhì)和ENO 性質(zhì).

圖15 ε=1.0×10?40 時(shí)使用WCNS7-Dm 格式計(jì)算雙馬赫反射問題Fig.15 Density of the double Mach reflection problem computed by the WCNS7-Dm scheme with ε=1.0×10?40

3.3.3 Richtmyer-Meshko 不穩(wěn)定性問題

Richtmyer-Meshko 不穩(wěn)定性問題[50]實(shí)際上是激波加速兩種密度不同的流體所導(dǎo)致的界面不穩(wěn)定現(xiàn)象.計(jì)算域?yàn)?[ 0,4]×[0,1],初始條件為

計(jì)算終止時(shí)間為t=9,網(wǎng)格數(shù)量為 3 21×81.

圖16 僅給出了 ε=1.0×10?40時(shí)WCNS7-Dm 格式在 λ=1.0×10?7和 λ=1下的計(jì)算結(jié)果,密度等值線是 [ 2,7.5]×λ等間距的24 條.結(jié)論與3.3.2 節(jié)相同.

圖16 ε=1.0×10?40 時(shí)使用WCNS7-Dm 格式計(jì)算Richtmyer-Meshko 不穩(wěn)定性問題Fig.16 Density of the Richtmyer-Meshkov instability problem computed by the WCNS7-Dm scheme with ε=1.0×10?40

4 結(jié)論

本文通過一系列的數(shù)值實(shí)驗(yàn)驗(yàn)證了尺度無關(guān)的WCNS7-JSm/Zm/Dm 格式的ENO 性質(zhì)、Si 性質(zhì)和Cp 性質(zhì),對(duì)比了新格式與經(jīng)典的WCNS 格式在激波捕捉能力和對(duì)復(fù)雜流場(chǎng)結(jié)構(gòu)的分辨率上的表現(xiàn).首先,從操作層面而言,新格式只需要在WCNS格式基礎(chǔ)上,使用函數(shù)的平均值構(gòu)建去尺函數(shù),再引入去尺度函數(shù)來修正非線性權(quán)重中的靈敏度參數(shù),這個(gè)過程并未引入任何復(fù)雜性.

其次,從設(shè)計(jì)思路來看,經(jīng)典的WCNS 格式不滿足Si 性質(zhì),因?yàn)樾〕叨认潞瘮?shù)的插值更容易在不連續(xù)區(qū)域附近產(chǎn)生數(shù)值振蕩,而尺度無關(guān)的WCNS 格式則是使格式性能與尺度因子和靈敏度參數(shù)無關(guān),首先將去尺度函數(shù)引入WCNS7-JS 和WCNS7-Z 格式的權(quán)重中以消除尺度的依賴性,衍生出尺度無關(guān)的WCNS7-JSm 和WCNS7-Zm 格式,但WCNS7-JSm 和WCNS7-Zm 格式僅滿足ENO 性質(zhì)和Si 性質(zhì),不滿足Cp 性質(zhì),為解決極值點(diǎn)附近出現(xiàn)降階現(xiàn)象,基于WCNS7-Z/Zm 格式,又提出了WCNS7-D/Dm 格式.

此外,從計(jì)算結(jié)果來看,新格式隨著網(wǎng)格加密能夠達(dá)到最優(yōu)精度;WCNS7-Dm 格式滿足所有3 種(ENO-,Si-和Cp-)性質(zhì)而不受靈敏度參數(shù)和尺度因子的影響;WCNS7-JSm/Zm 格式滿足ENO 性質(zhì)、Si 性質(zhì)和弱Cp 性質(zhì);經(jīng)典的WCNS7-JS/Z/D 格式只滿足上述3 種性質(zhì)中的某些性質(zhì)或較弱的情況.

綜上,尺度無關(guān)的WCNS 格式性能與靈敏度參數(shù)和尺度因子無關(guān),解決了經(jīng)典的WCNS 格式在小尺度下易產(chǎn)生數(shù)值振蕩的問題并且在激波捕捉能力和對(duì)復(fù)雜流場(chǎng)結(jié)構(gòu)的分辨率上表現(xiàn)良好,以期為WCNS 格式的改進(jìn)和應(yīng)用提供新的思路.接下來的工作將從兩個(gè)方面開展,首先從應(yīng)用層面,將進(jìn)一步測(cè)試尺度無關(guān)的WCNS 格式在具體算例中的數(shù)值表現(xiàn).其次從數(shù)值分析層面,去尺度函數(shù)的優(yōu)化工作也是下一步研究的重點(diǎn).