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

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

引言

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

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

作為對高階基本無振蕩ENO(essentially nonoscillatory)格式[11-12]的改進(jìn),Liu 等[13]提出的加權(quán)基本無震蕩WENO 格式是目前較為理想的高階激波捕捉格式之一.WENO 格式較好地兼顧了間斷區(qū)域的數(shù)值震蕩和光滑區(qū)域的計算精度,被廣泛應(yīng)用于各類流動的數(shù)值模擬[14–16].隨后,Jiang 等[17]引入了新的光滑度量因子,使WENO 格式達(dá)到最優(yōu)精度,該格式被命名為WENO-JS.針對經(jīng)典的WENO 格式所存在的極值點附近降階并且耗散過大等問題,典型的改進(jìn)方法主要包括: Henrick 等[18]提出映射權(quán)函數(shù)以擴大線性權(quán)的范圍,使得較大梯度下也能恢復(fù)線性權(quán),克服了極值點附近降階的問題但計算花費較大;Borges 等[19]通過引入高階全局光滑度量因子構(gòu)造的形式更為簡單的Z 型權(quán).Z 型權(quán)可以以較小的計算量獲得與原始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]、氣動聲學(xué)[27]、多組分流[28]以及激波-邊界層干擾[29]等.總的來說,基于變量插值的WCNS 格式相較于與Jiang 等[17]提出的基于通量函數(shù)重構(gòu)的WENO 格式相比有如下優(yōu)點[30-33]: (1)相同精度下分辨率更高;(2)數(shù)值通量選擇更加靈活[34];(3)自由流和渦流保持特性上表現(xiàn)更好[35].

然而,經(jīng)典的WCNS/WENO 格式仍存在一些明顯不足,如靈敏度參數(shù)和尺度因子的取值對格式的數(shù)值性能有很大影響.例如,由Deng 等[23-24,36]的研究可知,使用WCNS 格式離散Navier-Stokes 方程時,非線性加權(quán)函數(shù)中靈敏度參數(shù)被假定為 ε=1.0×10?6.但使用該格式離散雷諾應(yīng)力方程時[37],由于雷諾應(yīng)力的量級很小約為 1.0×10?6,ε=1.0×10?18才能保證計算穩(wěn)定進(jìn)行,而 ε=1.0×10?6或 ε=1.0×10?12將會導(dǎo)致計算發(fā)散.此外,當(dāng)使用靈敏度參數(shù)ε=1.0×10?12的W E N O 格式重 構(gòu)小尺 度因子λ=1.0×10?7下的Lax 問題時(Lax 問題的初始條件為Q(x,t=0)=λQ0(x)且Q0(x)=(ρ,ρu,E)),密 度ρ將會在間斷附近出現(xiàn)明顯的數(shù)值振蕩[38].因為使用5 階WENO 格式對小尺度下的不連續(xù)函數(shù)進(jìn)行重構(gòu)的過程中,相對較大的靈敏度參數(shù)對相對較小的光滑度量因子占主導(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 格式.研究內(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é)在雙精度算法下驗證了格式的基本無振蕩性質(zhì),并在4 精度算法下驗證了格式的尺度無關(guān)(scale-invariant,Si)性質(zhì)和極值點(critical point,Cp)性質(zhì);在文章的第3 節(jié)進(jìn)行了精度測試,并通過幾個一維和二維算例測試了新格式在激波捕捉以及分辨小尺度渦結(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]提出的單元中心型緊致格式計算一階導(dǎo)數(shù)

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

1.2 WCNS 插值格式

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

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

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

非線性加權(quán)技術(shù)是在流場光滑區(qū)令非線性權(quán)逼近線性權(quán),以保證設(shè)計精度,即 ωk=dk;在流場間斷區(qū)減小不光滑子模板的權(quán)重,以避免跨間斷插值引起數(shù)值震蕩和計算不收斂甚至發(fā)散等情況,此時 ω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點處1～3 階導(dǎo)數(shù)的平方和

式(11)包含一個線性項L和一個非線性項 Γk

τ

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

式(12)中 Φ 為修正函數(shù),相當(dāng)于激波探測器: 在流場光滑區(qū),?很小且 Φ=?;在間斷區(qū),?很大且 Φ=1,此時,非線性權(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ù),探索了一種更加簡單穩(wěn)定的非線性權(quán)函數(shù)的構(gòu)造方法,并將其應(yīng)用于7 階WCNS 格式,構(gòu)造尺度無關(guān)的WCNS7-JSm/Zm/Dm 格式,具體步驟如下.

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

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

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

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

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

修正Z/D 權(quán)函數(shù)非線性項中的靈敏度參數(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é)點上對流通量的計算使用Rusanov 格式[43],時間格式采用3 階Runge-Kutta,庫朗數(shù)為CFL=0.4.

2 性質(zhì)驗證

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

2.1 ENO 性質(zhì)

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

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

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

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

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

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

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

圖3 展示的是 ε=1.0×10?6和ε=1.0×10?40時的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 時使用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ì)

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

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

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

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

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

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

為驗證WCNS 格式的Si 性質(zhì),本節(jié)計算了求解雙曲守恒律中經(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ù).從圖中可以看出,對于任意 λ誤差ESi(λ)通常很小,在雙精度算法中,計算誤差與機器舍入誤差 ?0≈1.0×10?16將無法區(qū)分,故需要在四精度算法中進(jìn)行,機器的舍入誤差為 ?0≈1.0×10?34.

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

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

(2) 對于任意給定的靈敏度參數(shù) ε,WCNS7-D格式在臨界尺度因子時滿足弱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]給出極值點、Cp 性質(zhì)的定義.

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

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

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

式中 λ >0,函數(shù)在x=0且ncp=n處有一個n階極值點.圖7 給出了當(dāng) λ=1,ncp=1,2,3,4時WCNS 格式在 ε=1.0×10?6和 ε=1.0×10?40時的L∞誤差和收斂階,線條顏色(藍(lán)/粉/橙/紅/綠/黑)分別代表WCNS7-(JS/Z/D/JSm/Zm/Dm)格式.表1 給出了 λ=1.0×10?7和 λ=100 時,ε=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 在極值點 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 時的 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ù)較大時 ε=1.0×10?6,6 種WCNS格式都能達(dá)到設(shè)計精度O(?x7).

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

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

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

3.1 精度測試

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

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

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

圖8 ε=1.0×10?6和 ε=1.0×10?40 時分別使用WCNS 格式計算 λ=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 時使用WCNS 格式計算 λ=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 的右行激波與一系列密度波相互作用,可用來測試格式的激波捕捉能力、耗散性、數(shù)值穩(wěn)定性和高頻波分辨率,初始條件為[19]

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

簡便起見,圖10 僅給出了 λ=1.0×10?7且 ε=1.0×10?40時6 種WCNS 的計算結(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 格式實現(xiàn)了對過沖較好的抑制作用.曲線過沖或下沖可能會增大誤差,影響計算的穩(wěn)定性,故除存在過沖現(xiàn)象的WCNS7-D 格式(綠色),WCNS7-Dm 格式(紅色)相較于其他格式表現(xiàn)略好且分辨率略高.

圖10 ε=1.0×10?40 時使用WCNS 格式計算 λ=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 爆炸波相互作用問題

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

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

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

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

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

3.3 二維Euler 方程

3.3.1 二維Riemann 問題

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

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

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

圖13 ε=1.0×10?40 時使用WCNS 格式計算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 時使用WCNS 格式計算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)

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

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

圖14 ε=1.0×10?40 時使用WCNS7-D/Dm 格式計算 λ=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 時使用WCNS7-D/Dm 格式計算 λ=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)生的流場中包含豐富的旋渦結(jié)構(gòu)和十分復(fù)雜的激波反射,可用來測試數(shù)值格式捕捉激波和回流區(qū)小尺度結(jié)構(gòu)等能力.計算域為 [ 0,4]×[0,1],初始條件為

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

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

圖15 ε=1.0×10?40 時使用WCNS7-Dm 格式計算雙馬赫反射問題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]實際上是激波加速兩種密度不同的流體所導(dǎo)致的界面不穩(wěn)定現(xiàn)象.計算域為 [ 0,4]×[0,1],初始條件為

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

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

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

其次,從設(shè)計思路來看,經(jīng)典的WCNS 格式不滿足Si 性質(zhì),因為小尺度下函數(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ì),為解決極值點附近出現(xiàn)降階現(xiàn)象,基于WCNS7-Z/Zm 格式,又提出了WCNS7-D/Dm 格式.

此外,從計算結(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ù)值振蕩的問題并且在激波捕捉能力和對復(fù)雜流場結(jié)構(gòu)的分辨率上表現(xiàn)良好,以期為WCNS 格式的改進(jìn)和應(yīng)用提供新的思路.接下來的工作將從兩個方面開展,首先從應(yīng)用層面,將進(jìn)一步測試尺度無關(guān)的WCNS 格式在具體算例中的數(shù)值表現(xiàn).其次從數(shù)值分析層面,去尺度函數(shù)的優(yōu)化工作也是下一步研究的重點.