WENO 格式的一種新型光滑因子及其應用

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

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

0 引 言

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

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

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

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

傳統的光滑因子在光滑區被設計成盡量接近零的小量，最近Wu 等[27-28]認識到如果減小光滑區中候選模板的光滑因子之間的差異，不必滿足光滑因子盡量接近零的條件，就能夠減小WENO 格式的非線性權和線性權之間的差異。為了實現這一目標，Wu等[27-28]構造了一類光滑因子，使得各個候選模板上的光滑因子對于正弦函數為同一常數?；谠摴饣蜃油ㄟ^WENO-Z 型非線性權得到了滿足七階直至十七階精度的WENO 格式[27,28]，簡稱為WENO-S 格式。然而，由于算子函數定義的限制，WENO-S 格式至少為七階，不能得到在實際應用中需求更廣泛的五階及以下的WENO 格式。

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

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

1 控制方程與數值方法

1.1 控制方程

本文的數值方法主要應用于Euler 方程和Navier-Stokes 方程。以二維可壓縮守恒變量形式的無量綱Navier-Stokes 方程為例進行說明。

其中守恒變量為：

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

其中總能為：

粘性應力項為：

傳熱項為：

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

其中Tref為參考溫度。

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

1.2 數值方法

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

WENO 格式采用如下公式來離散導數項：

其中通過WENO 重構的方法來計算。對于系統方程，例如Euler 或者Navier-Stokes 方程，需要先把物理通量F投影到特征空間，然后對每一維的特征變量進行標量形式的WENO 重構過程。計算出標量形式的數值通量之后，再通過逆特征投影得到物理空間中的數值通量。詳細計算步驟可參考文獻[9,13]。

WENO 格式的本質是對標量形式的特征變量進行計算，因此下文重點以標量形式來說明五階WENO 格式的詳細計算過程。

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

WENO 格式需要對和分別進行計算，然后組合成數值通量

1.2.1 線性子重構格式和線性權

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

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

本文中所有的五階WENO 格式的子重構和線性權取值都是一致的，下文不再贅述。

1.2.2 WENO-JS 格式的光滑因子和非線性權

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

對其Taylor 展開，可得：

根據線性權和光滑因子可以計算出非線性權：

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

根據非線性權和三個子重構式（12）可以得到五階WENO-JS 格式：

1.2.3 WENO-Z 格式的光滑因子和非線性權

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

WENO-Z 格式的非線性權取為如下形式：

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

Borges 等[19]推導出WENO-Z 型非線性權滿足五階精度的充分條件為：

根據非線性權和三個子重構式（12）可以得到五階WENO-Z 格式：

1.2.4 WENO-NS 格式的光滑因子和非線性權

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

其中參數取為ξ=0.4。

五階全局光滑因子為：

對其Taylor 展開，可得：

WENO-NS 格式的非線性權取為如下形式：

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

根據非線性權和三個子重構式（12）可以得到五階WENO-NS 格式：

1.2.5 WENO-F 格式的光滑因子和非線性權

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

五階全局光滑因子為：

對其Taylor 展開，可得：

WENO-F 格式的非線性權取為如下形式：

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

根據非線性權和三個子重構式（12）可以得到五階WENO-F 格式：

1.3 新型光滑因子和非線性權

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

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

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

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

定義如下的算子函數[27,28]：

其中h是網格寬度。

通過Taylor 展開，可得：

因此，可以給出如下的算子函數關系：

其中μ為參數。

根據光滑因子設計的原理：當模板位于光滑區域，則光滑因子應該盡量接近零；當模板位于間斷區域，則光滑因子應該盡量接近無窮大。因此符合光滑因子的設計原理，本文基于該函數來設計一類新型光滑因子。

由于光滑因子必須為正數才能保證計算的穩定性，而上式并不能滿足該要求。經過數值試驗，如果直接取上式的絕對值作為光滑因子會導致計算嚴重不穩定。借鑒WENO-NS 和WENO-S 格式的設計經驗，我們取如下函數作為光滑因子：

其中，參數μ越大則格式的數值耗散越大；參數μ過小則會由于耗散不足而導致計算激波出現不穩定。經過大量的數值試驗，本文中取μ=0.025，該參數能符合WENO 格式對激波捕捉和模擬精細結構的需求。由這種新型光滑因子構造的WENO 格式簡記為WENO-XS 格式。

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

WENO-XS 格式的非線性權取為如下形式：

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

根據非線性權和三個子重構式（12）可以得到新型五階WENO-XS 格式：

2 精度分析和驗證

2.1 新型光滑因子和非線性權的精度分析

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

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

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

1）光滑函數問題。

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

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

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

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

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

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

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

2）包含間斷的問題。

間斷函數取為如下形式： 當?1

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

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

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

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

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

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

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

把以上的精度估計代入非線性權，可以得到如下的關系：

歸一化后的非線性權為：

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

當存在一階臨界點，即當fj′=0時，

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

2.2 數值精度驗證和計算效率

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

以一維Euler 方程的對流密度波問題為驗證標準[29]，其精確解為：

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

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

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

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

3 數值實驗與分析

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

3.1 Titarev–Toro 問題

本問題的控制方程是一維Euler 方程。Titarev–Toro 問題[30]流場中含有豐富的密度波結構，能夠考察計算格式對流場細節的分辨能力，是驗證數值格式的標準算例。該問題的初始條件為：

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

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

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

3.2 雙沖擊波問題

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

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

本問題分別采用五階WENO-JS 格式、WENOZ 格式、 WENO-F 格式、 WENO-NS 格式和新型WENO-XS 格式來進行計算，密度的分布如圖11 所示。計算網格采用250 個點。本問題中由雙沖擊波形成的間斷區域非常適合考察格式的激波捕捉特性。從計算結果可以看出，新型WENO-XS 格式在本問題中的計算結果與參考解最為接近， 對間斷的識別較為準確。WENO-F 格式比WENO-NS 格式計算結果略好，WENO-Z 格式計算結果比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 類型的非線性權，對于高頻波的分辨率和激波捕捉優于傳統的WENO-JS 類型的非線性權。因此，對于二維問題，不再給出WENO-JS 格式的計算結果。

3.3 激波雙馬赫反射問題

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

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

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

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

本問題的計算結果可以驗證新型WENO-XS 格式不僅能夠模擬滑移線等精細結構，同時能夠準確地捕捉激波。

3.4 前臺階問題

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

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

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

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

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

3.5 激波與剪切層相互作用

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

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

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

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

上邊界條件取激波波后的值，即取Ⅲ區的值。下邊界條件取為滑移壁面邊界條件。出口處用遠場特征邊界條件。

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

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

圖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

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

4 結 論

本文針對廣泛應用的五階WENO 格式，采用算子函數近似導數關系的方法，設計了一種新型光滑因子。在間斷區域，與傳統的光滑因子相比，新型光滑因子對模板的光滑性判斷更好，對間斷的識別更準確。在光滑區域，新型光滑因子使得新型WENO 格式的非線性權更接近線性權。

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

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

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

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