劉旭亮，武從海，李 虎，范召林

(1. 空氣動力學(xué)國家重點實驗室，綿陽 621000；2. 中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引 言

雙曲守恒方程是流體力學(xué)領(lǐng)域基本控制方程之一，在特定物理參數(shù)下，即使初始條件是光滑的，也可能出現(xiàn)激波和接觸間斷等不連續(xù)解，采用傳統(tǒng)的線性格式進(jìn)行數(shù)值求解會導(dǎo)致虛假振蕩。準(zhǔn)確模擬該類方程中的間斷及精細(xì)結(jié)構(gòu)需要開展非線性激波捕捉方法研究，為此，許多學(xué)者建立了一系列理論和方法。

為了消除間斷處的數(shù)值振蕩，Godunov[1]在1959 年首次提出迎風(fēng)格式的思想，這是一個里程碑性的近似解方法。后來，Harten[2-3]建立了總變差遞減（total variation diminishing, TVD）格式。Godunov 格式僅有一階精度，TVD 格式在數(shù)學(xué)上已經(jīng)被證明至多能達(dá)到二階精度[4]。為了抑制數(shù)值解在激波上、下游的振蕩，同時保證熵增條件，張涵信[5]提出了一種無波動、無自由參數(shù)的耗散差分算法，即二階精度的NND 格式。 由于 NND 格式具有TVD 性質(zhì)，物理概念明確且形式簡潔，在工程領(lǐng)域得到了廣泛的應(yīng)用。隨后，張涵信及其團隊進(jìn)一步發(fā)展了三階精度的ENN 格式[6]。為了進(jìn)一步提高精度，Harten 等[7-9]通過放寬TVD 條件，提出了有限體積框架的本質(zhì)無振蕩（essentially non-oscillatory, ENO）重構(gòu)格式。ENO 格式采用Newton 差商作為光滑因子，其基本思想是選擇Newton 差商較小的候選模板來進(jìn)行重構(gòu)，從而達(dá)到一致高階精度，并且在間斷和局部極值點附近具有本質(zhì)無振蕩的性質(zhì)。之后，Shu 和Osher[10-11]將ENO 格式推廣到了有限差分框架。

ENO 格式在重構(gòu)過程中計算了多個候選模板的光滑因子，但最后只選擇了最光滑的一個模板進(jìn)行重構(gòu)。在間斷區(qū)域，ENO 格式的思想是十分成功的；但是在光滑區(qū)域，從多個光滑候選模板中只選擇最光滑的模板則未充分利用其余候選模板的信息。為了把多個光滑模板都利用起來并且提高格式精度， Liu等[12]在1994 年基于ENO 格式發(fā)展了有限體積框架的加權(quán)ENO（weighted ENO, WENO）格式。WENO格式把所有候選模板以凸組合的方式進(jìn)行非線性加權(quán)，從而實現(xiàn)了一致高階精度，并且保持了間斷附近的本質(zhì)無振蕩性質(zhì)。非線性權(quán)值分配的基本策略是在光滑區(qū)域恢復(fù)到最優(yōu)精度的迎風(fēng)格式，在間斷模板分配較小的權(quán)值，實現(xiàn)本質(zhì)無振蕩的性質(zhì)。由于Liu 等[12]的WENO 格式采用Newton 差商的平方和作為光滑因子，導(dǎo)致WENO 格式在光滑區(qū)不能達(dá)到最優(yōu)精度[13]，例如五點WENO 格式在光滑區(qū)只能達(dá)到四階精度。1996 年，Jiang 和Shu[14]將WENO 格式推廣到了有限差分框架（簡稱為WENO-JS 格式），這項工作更重要的是通過所有候選模板的子重構(gòu)多項式導(dǎo)數(shù)的歸一化平方和（L2范數(shù)）構(gòu)造了Jiang-Shu 光滑因子。由于該光滑因子的良好特性，使得五點WENO-JS 格式在光滑區(qū)能達(dá)到五階最優(yōu)精度。隨后的研究者[15-17]基于Jiang-Shu 光滑因子得到了更高階的WENO-JS 格式，并將其應(yīng)用到求解間斷和復(fù)雜流動問題。

Henrick 等[18]對五階WENO-JS 格式進(jìn)行研究，結(jié)果表明當(dāng)存在高階臨界點（critical points）時，五階WENO-JS 格式的收斂精度會降至三階甚至更低。為了修正這一缺陷，Henrick 等[18]巧妙設(shè)計了一個映射函數(shù)，對WENO-JS 格式的非線性權(quán)進(jìn)行映射，使得非線性權(quán)在光滑區(qū)盡可能接近線性權(quán)。該格式被稱為映射WENO（mapped WENO）格式，簡稱WENOM 格式。相比于WENO-JS 格式，WENO-M 格式可以在光滑區(qū)域的臨界點達(dá)到最優(yōu)收斂精度，同時降低數(shù)值耗散。Borges 等[19]通過對子模板上的Jiang-Shu 光滑因子進(jìn)行線性組合，提出了全局光滑因子，他們將全局光滑因子與局部光滑因子之比引入到非線性權(quán)的設(shè)計中，稱為WENO-Z 格式。五階WENO-Z 格式在一階臨界點處能達(dá)到最優(yōu)收斂精度，并且由于WENO-Z 格式不需要計算映射函數(shù)，使得WENOZ 格式計算效率高于WENO-M 格式。隨后，Castro等[20]發(fā)展了高階WENO-Z 格式，給出了所有奇數(shù)階精度的全局光滑因子的公式。WENO-M 和WENOZ 格式在子模板上依然采用的是Jiang-Shu 光滑因子，本質(zhì)上都是導(dǎo)數(shù)項的L2范數(shù)。如果基于L1范數(shù)設(shè)計光滑因子，則可能會在光滑區(qū)域丟失精度[21]。為了克服這個缺陷，Ha 等[22]通過控制Taylor 展開余項來構(gòu)造L1范數(shù)意義下的局部和全局光滑因子，然后采用WENO-Z 型非線性權(quán)得到滿足最優(yōu)精度的WENO 格式，簡記為WENO-NS 格式。Ha 等[22]驗證了存在一階臨界點時，五階WENO-NS 格式依然能達(dá)到最優(yōu)精度。然而，WENO-NS 格式的光滑因子存在不對稱性，并且全局光滑因子的計算量比較大，隨后的研究工作[23-25]針對這些問題進(jìn)行了改進(jìn)。WENO格式計算量最大的部分在于光滑因子，為了降低計算成本，Baeza 等[26]構(gòu)造了一類高效計算的光滑因子，其形式比Jiang-Shu 光滑因子更簡潔，計算效率更高?；谠摴饣蜃?，通過WENO-Z 型非線性權(quán)得到了滿足最優(yōu)精度的快速WENO（fast WENO）格式，簡稱為WENO-F 格式。但是，七階以上的WENO-F格式的計算穩(wěn)定性略顯不足[27]。

傳統(tǒng)的光滑因子在光滑區(qū)被設(shè)計成盡量接近零的小量，最近Wu 等[27-28]認(rèn)識到如果減小光滑區(qū)中候選模板的光滑因子之間的差異，不必滿足光滑因子盡量接近零的條件，就能夠減小WENO 格式的非線性權(quán)和線性權(quán)之間的差異。為了實現(xiàn)這一目標(biāo)，Wu等[27-28]構(gòu)造了一類光滑因子，使得各個候選模板上的光滑因子對于正弦函數(shù)為同一常數(shù)。基于該光滑因子通過WENO-Z 型非線性權(quán)得到了滿足七階直至十七階精度的WENO 格式[27,28]，簡稱為WENO-S 格式。然而，由于算子函數(shù)定義的限制，WENO-S 格式至少為七階，不能得到在實際應(yīng)用中需求更廣泛的五階及以下的WENO 格式。

對于WENO 格式，構(gòu)造新的光滑因子至少需要滿足以下兩點要求：1）WENO 格式在光滑區(qū)域滿足精度條件[18-19]；2）在間斷區(qū)保持本質(zhì)無振蕩的計算穩(wěn)定性。

為了滿足上述光滑因子的設(shè)計要求，本文基于WENO-S 格式光滑因子對于特定函數(shù)為常數(shù)的準(zhǔn)則，設(shè)計了一種適用于五階WENO 格式的新型光滑因子。數(shù)值實驗表明新型WENO 格式不僅能夠準(zhǔn)確地捕捉強激波，而且能夠模擬剪切層和聲波等精細(xì)流場結(jié)構(gòu)。

1 控制方程與數(shù)值方法

1.1 控制方程

本文的數(shù)值方法主要應(yīng)用于Euler 方程和Navier-Stokes 方程。以二維可壓縮守恒變量形式的無量綱Navier-Stokes 方程為例進(jìn)行說明。

其中守恒變量為：

x、y方向的通量分別為：

其中總能為：

粘性應(yīng)力項為：

傳熱項為：

其中，Re和Ma分別為參考雷諾數(shù)和參考馬赫數(shù)，Pr=0.72，比熱比γ=1.4，黏性系數(shù)由Sutherland 公式計算：

其中Tref為參考溫度。

無量綱形式的完全氣體狀態(tài)方程和聲速方程分別為：

1.2 數(shù)值方法

本文研究的數(shù)值格式主要用于離散Euler 或者Navier-Stokes 方程的對流項，下面以x方向的離散為例進(jìn)行說明。

WENO 格式采用如下公式來離散導(dǎo)數(shù)項：

其中通過WENO 重構(gòu)的方法來計算。對于系統(tǒng)方程，例如Euler 或者Navier-Stokes 方程，需要先把物理通量F投影到特征空間，然后對每一維的特征變量進(jìn)行標(biāo)量形式的WENO 重構(gòu)過程。計算出標(biāo)量形式的數(shù)值通量之后，再通過逆特征投影得到物理空間中的數(shù)值通量。詳細(xì)計算步驟可參考文獻(xiàn)[9,13]。

WENO 格式的本質(zhì)是對標(biāo)量形式的特征變量進(jìn)行計算，因此下文重點以標(biāo)量形式來說明五階WENO 格式的詳細(xì)計算過程。

首先對通量f進(jìn)行分裂f(u)=f+(u)+f?(u)，滿足迎風(fēng)性條件：常用的Lax-Friedrichs 分 裂 為其 中

WENO 格式需要對和分別進(jìn)行計算，然后組合成數(shù)值通量

1.2.1 線性子重構(gòu)格式和線性權(quán)

在三個候選子模板S0={xj?2,···,x j}、S1={xj?1,···,xj+1}和S2={xj,···,xj+2}上的子重構(gòu)分別為：

從而可以計算出線性權(quán)值為：d0=1/10，d1=6/10，d2=3/10。

本文中所有的五階WENO 格式的子重構(gòu)和線性權(quán)取值都是一致的，下文不再贅述。

1.2.2 WENO-JS 格式的光滑因子和非線性權(quán)

五階WENO-JS 格式[14]在三個候選子模板上的光滑因子為具體為如下形式：

對其Taylor 展開，可得：

根據(jù)線性權(quán)和光滑因子可以計算出非線性權(quán)：

其中，ε用于防止分母為零，一般取值為ε=1×10?6。歸一化后的非線性權(quán)為：

根據(jù)非線性權(quán)和三個子重構(gòu)式（12）可以得到五階WENO-JS 格式：

1.2.3 WENO-Z 格式的光滑因子和非線性權(quán)

WENO-Z 格式[19,20]在三個候選子模板上的光滑因子與WENO-JS 格式保持一致，而WENO-Z 格式的重要創(chuàng)新之處是引入了全局光滑因子τ5，它表征了總模板的光滑性，其中五階全局光滑因子為

WENO-Z 格式的非線性權(quán)取為如下形式：

其中，ε用于防止分母為零，一般取值為ε=1×10?40。歸一化后的非線性權(quán)為：

Borges 等[19]推導(dǎo)出WENO-Z 型非線性權(quán)滿足五階精度的充分條件為：

根據(jù)非線性權(quán)和三個子重構(gòu)式（12）可以得到五階WENO-Z 格式：

1.2.4 WENO-NS 格式的光滑因子和非線性權(quán)

五階WENO-NS 格式[22]在三個候選子模板上的光滑因子取為如下形式：

其中參數(shù)取為ξ=0.4。

五階全局光滑因子為：

對其Taylor 展開，可得：

WENO-NS 格式的非線性權(quán)取為如下形式：

其中，ε用于防止分母為零，一般取值為ε=1×10?40。歸一化后的非線性權(quán)為：

根據(jù)非線性權(quán)和三個子重構(gòu)式（12）可以得到五階WENO-NS 格式：

1.2.5 WENO-F 格式的光滑因子和非線性權(quán)

五階WENO-F 格式[26]在三個候選子模板上的光滑因子取為如下形式：

五階全局光滑因子為：

對其Taylor 展開，可得：

WENO-F 格式的非線性權(quán)取為如下形式：

其中，ε用于防止分母為零，一般取值為ε=1×10?40。歸一化后的非線性權(quán)為：

根據(jù)非線性權(quán)和三個子重構(gòu)式（12）可以得到五階WENO-F 格式：

1.3 新型光滑因子和非線性權(quán)

Wu 等[27,28]針對描述波動屬性的基本函數(shù)f(x)=sin(x)，得到其導(dǎo)數(shù)關(guān)系：對于r≥4， (f(r?2))2?f(r?3)f(r?1)=1。由此發(fā)展了一類適用于多尺度流動問題的七階及以上WENO 格式的光滑因子。通過簡單的導(dǎo)數(shù)運算，可以發(fā)現(xiàn)當(dāng)r= 3 時， (f′)2?f f′′=1也是成立的。本文以此為基礎(chǔ)，發(fā)展了五階WENO 格式的光滑因子。

首先將導(dǎo)數(shù)關(guān)系推廣到更一般的正弦函數(shù)f(x)=Asin(ωx+φ)，經(jīng)過簡單的導(dǎo)數(shù)運算可以得到(f′)2?f f′′=A2ω2是常數(shù)。然后將導(dǎo)數(shù)關(guān)系推廣到更一般的情形(f′)2?μf f′′=C，其中μ和C為常數(shù)。求解該微分方程可以得到對應(yīng)的函數(shù)形式。

以下舉例說明。對于微分方程(f′)2?f f′′=1，還可以得到另一個特解為f(x)=sinh(x)。對于微分方程(f′)2?2f f′′=0，可以得到一個特解為f(x)=x2。對于微分方程(f′)2?μf f′′=C2，其中μ和C為任意常數(shù)，可以得到一個特解為f(x)=Cx。對于微分方程(f′)2?μf f′′=0，其中μ= 0.025，可以得到一個特解為

因此，當(dāng)實際計算遇到的函數(shù)為這些特解時，對于特定的參數(shù)μ和C，就能滿足導(dǎo)數(shù)關(guān)系(f′)2?μf f′′=C。

定義如下的算子函數(shù)[27,28]：

其中h是網(wǎng)格寬度。

通過Taylor 展開，可得：

因此，可以給出如下的算子函數(shù)關(guān)系：

其中μ為參數(shù)。

根據(jù)光滑因子設(shè)計的原理：當(dāng)模板位于光滑區(qū)域，則光滑因子應(yīng)該盡量接近零；當(dāng)模板位于間斷區(qū)域，則光滑因子應(yīng)該盡量接近無窮大。因此符合光滑因子的設(shè)計原理，本文基于該函數(shù)來設(shè)計一類新型光滑因子。

由于光滑因子必須為正數(shù)才能保證計算的穩(wěn)定性，而上式并不能滿足該要求。經(jīng)過數(shù)值試驗，如果直接取上式的絕對值作為光滑因子會導(dǎo)致計算嚴(yán)重不穩(wěn)定。借鑒WENO-NS 和WENO-S 格式的設(shè)計經(jīng)驗，我們?nèi)∪缦潞瘮?shù)作為光滑因子：

其中，參數(shù)μ越大則格式的數(shù)值耗散越大；參數(shù)μ過小則會由于耗散不足而導(dǎo)致計算激波出現(xiàn)不穩(wěn)定。經(jīng)過大量的數(shù)值試驗，本文中取μ=0.025，該參數(shù)能符合WENO 格式對激波捕捉和模擬精細(xì)結(jié)構(gòu)的需求。由這種新型光滑因子構(gòu)造的WENO 格式簡記為WENO-XS 格式。

類似于三點模板的光滑因子推導(dǎo)過程，可以得到五點模板的新型五階全局光滑因子：

WENO-XS 格式的非線性權(quán)取為如下形式：

其中，ε用于防止分母為零，本文取值為ε=1×10?40。歸一化后的非線性權(quán)為：

根據(jù)非線性權(quán)和三個子重構(gòu)式（12）可以得到新型五階WENO-XS 格式：

2 精度分析和驗證

2.1 新型光滑因子和非線性權(quán)的精度分析

對新型光滑因子式（39）進(jìn)行Taylor 展開，可以得到如下的精度估計：

通過對比分析發(fā)現(xiàn)，如果忽略常數(shù)系數(shù)，WENO-JS、 WENO-Z、 WENO-NS 和 WENO-F 格式的光滑因子的Taylor 展開主項都是而新型光滑因子的Taylor 展開主項為可見新型光滑因子對候選模板的光滑性判斷多出一項

以下從光滑函數(shù)問題和包含間斷的問題兩方面來分別說明：

1）光滑函數(shù)問題。

取f(x)=sin(πx)，其計算域為?1 ≤x≤1，計算網(wǎng)格取為N= 20，網(wǎng)格寬度為h= 0.1。

圖1 給出了WENO5-XS 格式光滑因子的兩個分量在計算域內(nèi)的分布，從圖中可以發(fā)現(xiàn)在光滑函數(shù)的極值附近迅速減小，在其余區(qū)域的值則顯著大于分量基本保持常數(shù)，這是由于對于函數(shù)滿足導(dǎo)數(shù)關(guān)系式為常數(shù)，使得其中C為常數(shù)。 采用這兩個分量作為光滑因子，分別計算得到的非線性權(quán)記為其分布如圖2 所示。對于光滑函數(shù)，非線性權(quán)偏離線性權(quán)（d1=0.6）嚴(yán)重，顯然是不符合設(shè)計要求的?？偟姆蔷€性權(quán)ω1與線性權(quán)比較接近，是滿足光滑問題的高精度計算需求的。

圖1 對光滑函數(shù)WENO5-XS 格式光滑因子的兩個分量Fig. 1 Smoothness indicator and its components of WENO5-XS scheme in smooth region

圖2 對光滑函數(shù)WENO5-XS 格式光滑因子的兩個分量計算的非線性權(quán)Fig. 2 Nonlinear weights of WENO5-XS scheme in smooth region

不同五階WENO 格式的光滑因子和非線性權(quán)的比較如圖3、圖4 所示，可以發(fā)現(xiàn)WENO5-JS 格式的非線性權(quán)ω1與線性權(quán)d1的偏離較大，其余格式非線性權(quán)與線性權(quán)的差異較小。

圖3 對光滑函數(shù)不同五階WENO 格式的光滑因子比較Fig. 3 Smoothness indicators of different WENO5 schemes in smooth region

圖4 對光滑函數(shù)不同五階WENO 格式的非線性權(quán)比較Fig. 4 Nonlinear weights of different WENO5 schemes in smooth region

2）包含間斷的問題。

間斷函數(shù)取為如下形式： 當(dāng)?1

圖5 給出了WENO5-XS 格式光滑因子的兩個分量在計算域內(nèi)的分布，從圖中可以發(fā)現(xiàn)：對于光滑區(qū)域，IS1(1)在函數(shù)極值附近迅速減小，在其余光滑區(qū)域IS1(1)的值則顯著大于IS1(2)。對于間斷區(qū)域，IS1(1)和IS1(2)都能準(zhǔn)確地識別間斷點，并且IS1(1)的值比IS1(2)更大。采用IS1(1)和IS1(2)這兩個分量作為光滑因子，分別計算得到非線性權(quán)記為ω1(1)和ω1(2)，其分布如圖6 所示。在間斷處，非線性權(quán)ω1(2)的數(shù)值比ω(11)更小，說明IS1(2)對間斷的識別起著重要作用，因此光滑因子兩個分量之和IS1能更準(zhǔn)確地識別間斷。

圖5 對間斷函數(shù)WENO5-XS 格式光滑因子的兩個分量Fig. 5 Smoothness indicator and its components of WENO5-XS scheme in discontinuous region

圖6 對間斷函數(shù)WENO5-XS 格式光滑因子的兩個分量計算的非線性權(quán)Fig. 6 Nonlinear weights of WENO5-XS scheme in discontinuous region

不同五階WENO 格式的光滑因子和非線性權(quán)的比較如圖7、圖8 所示，可以發(fā)現(xiàn)間斷點處WENO5-JS 格式的非線性權(quán)值ω1相對較小，對間斷的抹平最嚴(yán)重；在間斷點處WENO5-XS 格式的非線性權(quán)值ω1相對較大，對間斷的分辨率更好。

圖7 對間斷函數(shù)不同五階WENO 格式的光滑因子比較Fig. 7 Smoothness indicators of different WENO5 schemes in discontinuous region

圖8 對間斷函數(shù)不同五階WENO 格式的非線性權(quán)比較Fig. 8 Nonlinear weights of different WENO5 schemes in discontinuous region

對新型五階全局光滑因子進(jìn)行Taylor 展開，可以得到：

把以上的精度估計代入非線性權(quán)，可以得到如下的關(guān)系：

歸一化后的非線性權(quán)為：

當(dāng)不存在臨界點，即當(dāng)fj′≠0時，式（47）滿足五階精度的充分條件式（20），此時WENO-XS 格式達(dá)到五階最優(yōu)精度要求。

當(dāng)存在一階臨界點，即當(dāng)fj′=0時，

同樣的推導(dǎo)方法可以得到：ωXrS=dr+Oh4，因此仍然滿足五階精度的充分條件式（20）。

2.2 數(shù)值精度驗證和計算效率

為了驗證數(shù)值方法的精度，我們選取如下的格式：時間離散采用三階TVD Runge-Kutta 格式[10]，通量分裂采用Lax-Friedrichs 分裂，空間離散采用五階WENO-XS 差分格式。

以一維Euler 方程的對流密度波問題為驗證標(biāo)準(zhǔn)[29]，其精確解為：

計算域取為[0,2]，邊界取為周期性邊界條件。計算的最終時刻取為100 個周期t=200，初始的計算網(wǎng)格為80 個點，CFL = 0.5。隨著網(wǎng)格點數(shù)增大，CFL 數(shù)取為其中下標(biāo)“1”表示上一個時間層的值，下標(biāo)“2”表示下一個時間層的值，N表示網(wǎng)格點數(shù)，n= 5 表示空間離散精度。計算得到的數(shù)值結(jié)果如表1 所示。該問題包含一階臨界點，數(shù)值結(jié)果驗證了WENO-XS 格式在一階臨界點處也達(dá)到了五階的設(shè)計精度。

表1 五階WENO-XS 格式的數(shù)值精度驗證Table 1 Numerical accuracy of WENO5-XS scheme

采用CPU 計算時間和L1平均誤差作為綜合因素來評價數(shù)值格式的計算效率。計算結(jié)果如圖9 所示，可以發(fā)現(xiàn)計算效率從高到低依次是 WENO5-F、WENO5-XS、WENO5-Z、WENO5-NS 和WENO5-JS 格式。其中WENO5-F 格式的計算效率最高，主要得益于WENO5-F 的光滑因子的形式最簡潔。WENO5-JS 格式計算效率最低，主要是由于在相同網(wǎng)格下其計算誤差較大。WENO5-XS 格式的計算效率比較高，這是因為其光滑因子的形式比WENO5-JS 和WENO5-NS 的形式更簡潔，并且WENO5-XS 格式的三個光滑因子在計算時只需平移一個坐標(biāo)點，此優(yōu)點更利于編程和數(shù)值計算。

圖9 不同五階WENO 格式計算效率的比較Fig. 9 The comparison of computational efficiency for different fifth-order WENO schemes

3 數(shù)值實驗與分析

為了驗證新型WENO-XS 格式在雙曲守恒方程中的計算準(zhǔn)確度和穩(wěn)定性，本文計算了一維Euler 方程、二維Euler 方程和二維Navier-Stokes 方程。時間離散采用三階TVD Runge-Kutta 格式[10]，通量分裂采用Lax-Friedrichs 分裂，空間離散采用多種不同的五階WENO 格式用于比較計算結(jié)果。

3.1 Titarev–Toro 問題

本問題的控制方程是一維Euler 方程。Titarev–Toro 問題[30]流場中含有豐富的密度波結(jié)構(gòu)，能夠考察計算格式對流場細(xì)節(jié)的分辨能力，是驗證數(shù)值格式的標(biāo)準(zhǔn)算例。該問題的初始條件為：

計算區(qū)域取為[?5, 5]，初始間斷位于x=?4.5處，最終計算時刻取為t= 5。邊界處采用初始邊界值。采用10000 個網(wǎng)格點的五階WENO-JS 格式的計算結(jié)果來作為參考解。

本問題分別采用五階WENO-JS 格式、WENOZ 格式、 WENO-F 格式、 WENO-NS 格式和新型WENO-XS 格式來進(jìn)行計算，密度波的分布如圖10所示。為了便于比較，計算網(wǎng)格采用1000 個點，使得密度波的一個周期內(nèi)有10 個點。由于高頻的密度波區(qū)域是光滑區(qū)域，非常適合考察格式的短波分辨率和耗散性。從計算結(jié)果可以看出，新型WENO-XS 格式在本問題中的計算結(jié)果與參考解最為接近，這說明WENO-XS 格式在光滑區(qū)的非線性權(quán)更接近線性權(quán)，分辨率更高，耗散更小。WENO-NS 格式的計算結(jié)果次之， WENO-F 格式比WENO-Z 格式計算結(jié)果略好，而WENO-JS 格式在本問題中對高頻密度波的分辨率結(jié)果相對較差。

圖10 不同五階WENO 格式求解Titarev-Toro 問題Fig. 10 Solutions of the Titarev-Toro problem obtained by different fifth-order WENO schemes

3.2 雙沖擊波問題

本問題的控制方程是一維Euler 方程。雙沖擊波問題[31]包括了強激波、接觸間斷和膨脹波之間的相互作用，對計算格式的穩(wěn)定性要求比較高。其初值為：

計算域兩端采用反射邊界條件，計算終止時間為t= 0.038，采用10000 個網(wǎng)格點的五階WENOJS 格式的計算結(jié)果來作為該問題的參考解。

本問題分別采用五階WENO-JS 格式、WENOZ 格式、 WENO-F 格式、 WENO-NS 格式和新型WENO-XS 格式來進(jìn)行計算，密度的分布如圖11 所示。計算網(wǎng)格采用250 個點。本問題中由雙沖擊波形成的間斷區(qū)域非常適合考察格式的激波捕捉特性。從計算結(jié)果可以看出，新型WENO-XS 格式在本問題中的計算結(jié)果與參考解最為接近， 對間斷的識別較為準(zhǔn)確。WENO-F 格式比WENO-NS 格式計算結(jié)果略好，WENO-Z 格式計算結(jié)果比WENO-JS 格式好。

圖11 不同五階WENO 格式求解雙沖擊波問題Fig. 11 Solutions of two interacting blast waves obtained bydifferent fifth-order WENO schemes

由于WENO-Z 格式、WENO-F 格式、WENO-NS格式和新型WENO-XS 格式采用的都是WENO-Z 類型的非線性權(quán)，對于高頻波的分辨率和激波捕捉優(yōu)于傳統(tǒng)的WENO-JS 類型的非線性權(quán)。因此，對于二維問題，不再給出WENO-JS 格式的計算結(jié)果。

3.3 激波雙馬赫反射問題

本問題的控制方程是二維Euler 方程。雙馬赫反射問題[31]包含強激波和滑移線，非常適合于考察格式的激波捕捉能力和流場精細(xì)結(jié)構(gòu)的分辨率。本問題描述的是馬赫數(shù)為10 的強運動斜激波以與x軸方向呈60°角的方向入射，入射點在（1/6,0），計算區(qū)域取[0,4] × [0,1]。激波波前靜止空氣的密度為1.4，壓力為1，波后按激波關(guān)系給定初始條件。初始物理參數(shù)為：

下邊界在[1/6,4]的區(qū)域給定壁面反射邊界條件，其它邊界按照激波運動所在的位置分別給定波前或波后的值。本文采用1500×375 的均勻網(wǎng)格計算到無量綱時間t= 0.2。

本問題分別采用五階WENO-Z 格式、WENOF 格式、WENO-NS 格式和新型WENO-XS 格式來進(jìn)行計算，結(jié)果如圖12 所示，取密度為1.6～21 共30 條等值線作圖。本問題比較不同格式主要關(guān)注馬赫桿和滑移線附近的分辨率，因此為了更好地比較不同的計算結(jié)果，在圖中僅給出渦卷起區(qū)域（blow-up region）。從計算結(jié)果可以看出，新型WENO-XS 格式得到的旋渦結(jié)構(gòu)是最豐富的，對流場精細(xì)結(jié)構(gòu)的分辨率最高， WENO-NS 格式比WENO-F 格式計算結(jié)果略好，WENO-Z 格式的計算結(jié)果相對較差。

圖12 不同五階WENO 格式求解雙馬赫反射問題Fig. 12 Solutions of the double Mach reflection problem obtained by different fifth-order WENO schemes

本問題的計算結(jié)果可以驗證新型WENO-XS 格式不僅能夠模擬滑移線等精細(xì)結(jié)構(gòu)，同時能夠準(zhǔn)確地捕捉激波。

3.4 前臺階問題

本問題的控制方程是二維Euler 方程。前臺階問題[31]的初始條件為馬赫數(shù)為3 的強運動激波沿著帶有臺階的風(fēng)洞向右傳播。臺階高度為0.2，前緣位置在0.6，計算域為[0,3] × [0,1]。初始物理參數(shù)為：

左、右邊界為開邊界，其余邊界取壁面反射條件。本文采用720×240 的均勻網(wǎng)格，計算截止到無量綱時間t= 4 的時刻。

本問題分別采用五階WENO-Z 格式、WENOF 格式、WENO-NS 格式和新型WENO-XS 格式來進(jìn)行計算，結(jié)果如圖13 所示，取密度為0.2568～6.607共60 條等值線作圖。本問題比較不同格式主要關(guān)注滑移線附近的分辨率。從計算結(jié)果可以看出，新型WENO-XS 格式識別的旋渦結(jié)構(gòu)是最豐富的，對流場精細(xì)結(jié)構(gòu)的分辨率最高，WENO-NS 格式的計算結(jié)果次之，WENO-F 格式與WENO-Z 格式的計算結(jié)果相當(dāng)。

圖13 不同五階WENO 格式求解前臺階問題Fig. 13 Flow fields around a forward facing step obtained by different fifth-order WENO schemes

本問題的計算結(jié)果可以驗證新型WENO-XS 格式能夠準(zhǔn)確地捕捉強激波，同時能夠分辨出滑移線附近更多的旋渦精細(xì)結(jié)構(gòu)，說明WENO-XS 格式的耗散是相對最小的。

3.5 激波與剪切層相互作用

本問題的控制方程是二維可壓縮Navier-Stokes方程。本問題是典型的激波噪聲問題，在論文[32，33]中有詳細(xì)的描述和分析。以下簡要說明計算的初始條件和邊界條件。整個剪切層在初始時刻是等壓的，并且總的焓值也處處相等。計算域取為[0,200]×[?20,20]。初始時刻可以把計算域分成如圖14 所示的5 個部分，表2 給出了各個區(qū)域具體的流動參數(shù)。

表2 初始時刻流場的物理參數(shù)Table 2 Parameters of the initial flow field

圖14 初始流場計算域的劃分Fig. 14 Decomposition of the computational domain at the initial state

網(wǎng)格在x方向取均勻網(wǎng)格，在y方向取拉伸網(wǎng)格，拉伸變換關(guān)系為y方向的計算域長度為Ly=40，拉伸因子by=1，計算網(wǎng)格η∈[?1,1]為均勻網(wǎng)格。

上邊界條件取激波波后的值，即?、髤^(qū)的值。下邊界條件取為滑移壁面邊界條件。出口處用遠(yuǎn)場特征邊界條件。

本問題采用1000×200 的網(wǎng)格規(guī)模，計算截止到無量綱時間t= 120 的時刻。

采用3800×900 網(wǎng)格點的五階WENO-JS 格式的計算結(jié)果作為該問題的參考解。本問題分別采用五階WENO-Z 格式、WENO-F 格式、WENO-NS 格式和新型WENO-XS 格式來進(jìn)行計算，結(jié)果如圖15 所示，取密度陰影為?3.5～0.5 共30 條等值線作圖。本問題比較不同格式主要關(guān)注剪切層和聲波的分辨率。新型WENO-XS 格式對超聲速剪切層外部的小激波結(jié)構(gòu)識別相對較為清晰，對聲波的分辨率相對較好。為了能夠更清晰地顯示不同格式之間的差別，取t= 120 時刻橫向中心線位置（y= 0）的密度剖面進(jìn)行比較（圖16）。從密度剖面的計算結(jié)果可以看出，新型WENO-XS 格式在本問題中的表現(xiàn)是最好的，其計算結(jié)果與參考解最為接近， WENONS 格式的計算結(jié)果比WENO-F 格式略好， WENOZ 格式的計算結(jié)果相對較差。

圖15 不同五階WENO 格式求解激波剪切層相互作用問題Fig. 15 The interaction of shock wave and shear layer obtained by different fifth-order WENO schemes

圖16 不同五階WENO 格式沿y = 0 的密度分布Fig. 16 The distribution of density along y = 0 obtained by different fifth-order WENO schemes

計算結(jié)果可以驗證新型WENO-XS 格式不僅能精細(xì)模擬剪切層和聲波，同時能準(zhǔn)確捕捉強激波。

4 結(jié) 論

本文針對廣泛應(yīng)用的五階WENO 格式，采用算子函數(shù)近似導(dǎo)數(shù)關(guān)系的方法，設(shè)計了一種新型光滑因子。在間斷區(qū)域，與傳統(tǒng)的光滑因子相比，新型光滑因子對模板的光滑性判斷更好，對間斷的識別更準(zhǔn)確。在光滑區(qū)域，新型光滑因子使得新型WENO 格式的非線性權(quán)更接近線性權(quán)。

理論分析和數(shù)值驗證表明新型WENO 格式具有五階精度，即使在一階臨界點也能保持一致五階精度。

對于一維激波管問題，通過與WENO-JS、WENOZ、WENO-F 和WENO-NS 格式的比較，從計算結(jié)果可以發(fā)現(xiàn)新型WENO-XS 格式與參考解最為接近，對短波的分辨率最高，格式的耗散性最小，同時對激波的間斷識別最準(zhǔn)確，具有良好的激波捕捉特性。

對于二維包含激波的流動問題，通過與其他WENO 格式的比較，計算結(jié)果表明新型WENO-XS 格式不僅能夠準(zhǔn)確地捕捉強激波，而且能夠模擬剪切層和聲波等精細(xì)流場結(jié)構(gòu)，得到的旋渦結(jié)構(gòu)是最豐富的，對旋渦和聲波的分辨率最高。

本文提出的五階WENO-XS 格式在求解包含激波的復(fù)雜流動問題中具有潛在的應(yīng)用價值。提出的新型光滑因子構(gòu)造思路不局限于五階，可以推廣到更高階的WENO 格式，相關(guān)的工作還需進(jìn)一步研究探索。