非結構網格通量重構算法下三種緊致WENO 限制器對比研究

石京昶，嚴 紅，*

(1. 西北工業大學 長三角研究院，太倉 215400；2. 陜西省航空發動機內流動力學重點實驗室，西安 710072)

0 引 言

高精度數值格式因日益強烈的工程需求被關注和研究。其中代表性的方法有加權本質無振蕩（weighted essentially non-oscillatory, WENO）、間斷伽遼金（discontinuous Galerkin methods, DG）、通量重構（flux reconstruction, FR）等。WENO 方法通過非線性重構保證了激波附近“本質無振蕩”。DG 和FR 方法天然適用于非結構網格處理復雜幾何外形，然而DG、FR 方法本身不能處理激波，所以需要引入限制器。

DG、FR 方法框架下激波捕捉方法的一般思路是：首先采用間斷探測器找到需要進行重構的網格單元，然后對目標網格單元采用限制器重構出新的解多項式。如何檢測間斷有不同的思路，參見Qiu 等的對比研究[1]。其中KXRCF 間斷檢測器[2]被廣泛用于高階DG 格式的激波捕捉方法研究中。有限體積方法下發展出的TVD、TVB 限制器可以應用于高階DG、FR 方法，但會丟失高階方法的高階特性。WENO 的思想可以用來構造適用于DG 和FR 方法的限制器。Qiu 等[3]和Luo 等[4]提出了Hermite WENO 限制器，但是并不緊致，目標單元解的重構需要用到相鄰網格的相鄰網格中的解。這種非緊致特性造成DG、FR 方法并行擴展性損失。Zhong 等[5]首先提出了緊致WENO 限制器用于DG 方法，核心思想是將相鄰網格單元中的解插值到目標單元，并將單元平均替換為目標單元的單元平均，形成的新多項式與目標單元原本的多項式以WENO 的方式加權平均得到最終的重構多項式解。此WENO 限制器思路清晰，實現簡單，故稱為簡單WENO 限制器（簡稱SWENO 限制器）。Zhu 等[6-7]在此基礎上加以改進，將簡單WENO 限制器中目標單元的相鄰單元中的解直接插值到目標單元的方式改進為相鄰單元中的解通過最小二乘法投影到目標單元。最近，Zhu 等[8-9]又進一步提出了新的多精度WENO 限制器（簡稱MWENO限制器），其相對于改進的簡單WENO 限制器做出的核心改進是，目標單元中的原始高階多項式分解成從一階至高階的多個不同階次多項式，相鄰單元中的解多項式經最小二乘法投影到目標單元，形成一階多項式后，用于計算最低階多項式的光滑因子，最終多個不同階次多項式仍采用WENO 方式加權重構得到最終的解多項式。Li 等[10]也于最近提出了p階加權WENO 限制器（簡稱PWENO 限制器）。核心思想與Zhu 等最新的多精度WENO 限制器類似，不過相鄰網格單元中的解經最小二乘法投影到目標單元，形成一階多項式，參與WENO 重構得到最終的解多項式，而不是只用于計算最低階多項式的光滑因子。以上WENO 限制器在極端低密度低壓的局部區域會重構出包含負密度負壓的非物理解，因此需要使用保正限制器避免此情形出現。Zhang 等[11]提出的保正限制器能夠保持原始高階精度無損失，被廣泛應用于高階DG 格式激波捕捉方法研究中。但值得一提的是，此保正保精度限制器證明過程中假設單元解點必須是Gauss-Lobatto 點，相比于Gauss-Legendre 點精度較低。

高階格式計算含有激波的穩態問題時較難收斂到機器零，WENO 格式存在此問題。Zhang 等[12]提出迎風型重構減弱激波后偽振蕩。Zhu 等[13]構造了新的多精度WENO 格式。DG 方法結合WENO 限制器激波捕捉同樣難以收斂到機器零。最近，Zhu 等[14]提出了一種新的間斷探測器，結合其多精度WENO限制器能夠在標準算例中獲得良好的收斂特性。該間斷探測器核心思想是將KXRCF 中基于單元面上解計算差值改為基于單元內解計算差值，本文稱之為CKXRCF 間斷探測器。

本文在FR 框架下，結合間斷探測器和保正保精度限制器，對比研究簡單WENO 限制器、p階加權WENO 限制器和多精度WENO 限制器在多個經典算例中的性能，并討論其高精度的優勢和收斂性問題。常見文獻中WENO 限制器驗證測試算例集中在無黏Euler 方程，本文將其擴展到激波與層流邊界層相互干擾的算例中，針對N-S 方程測試以上WENO 限制器的性能，為相關研究提供初步結果。

1 FR 方法簡介

非定?？蓧嚎sNavier-Stokes 方程的微分形式可寫成如下雙曲型守恒律的形式：

式中，通量F(Q)=Fc(Q)?Fv(Q)，守恒變量Q、無黏通量和黏性通量分別如下所示：

本文簡介FR 方法二維情形下的基本思想，詳細闡述見Wang 等[15]和Huynh 等[16]的綜述。計算域離散成N個網格單元{Vi}iN=1。方程（1）的加權弱形式可改寫為以下形式：

式中，Qi是網格單元Vi上的近似解，且屬于k階多項式空間，即Qi∈Pk(Vi)。Π[?·F(Qi)]是?·F(Qi)在多項式空間Pk上的投影。δi∈Pk(Vi)是Vi上的修正項，由式（3）定義：

式中，W是權函數，[Fn]=Fcnom?Fn(Qi)是網格界面上黎曼通量與單元內部局部通量的差值。將自由度定義在網格單元內部的解點。為計算修正項，在單元界面處定義通量點，則解點上的修正項由式（4）計算：

式中，αj,f,l是獨立于解的常數，Sf是單元界面面積，j是解點在當前網格單元內的編號，f是當前網格單元面的編號，l是通量點在當前面的編號。綜上所述，FR 方法由以下離散方程近似求解原雙曲律方程：

本文所有算例均使用自主開發的N-S 方程非結構網格并行求解器NFR，時間推進方法均采用三階強穩定性顯式SSP-RK3 方法，無黏通量采用Roe 通量或Lax–Friedrichs 通量，黏性通量均采用BR2 格式[17]。本文所有算例的結果均將原始網格單元多項式階次均等剖分為子單元，如一個四邊形網格單元的P2 剖分為4 個四邊形子單元。

2 WENO 限制器簡介

本文采用的WENO 限制器的基本思路是使用間斷探測器找出需要限制解的問題網格單元，然后對其應用WENO 限制器重構得到新解。對于歐拉方程系統，為了更好地避免激波附近的振蕩，使用通量雅可比矩陣將守恒變量轉換為特征變量，經WENO 限制器重構得到新的特征變量，再經通量雅可比矩陣轉換回守恒變量。下述3 種WENO 限制器均只簡介其針對一維標量方程的限制過程。

2.1 簡單WENO 限制器

簡單WENO 限制器的核心是將問題網格單元Ij及其相鄰網格單元的解多項式加權組合得到重構解，保證其與原始解有相同的單元平均和k階精度。權重由解多項式的光滑因子確定。

記網格單元Ij?1、Ij、Ij+1的解為pL(x)、p0(x)、pR(x)，簡單WENO 限制器重構的新解為：

式中，ωL、ωR是歸一化的非線性權重，(x)、(x)是相鄰網格單元內解多項式pL(x)、pR(x)投影到Ij中的新多項式。為保證守恒性，相鄰網格單元內的解向Ij中的投影形式為：

式中，單元平均定義為：

注意以上積分均為目標單元Ij內的積分。

上述歸一化非線性權重ωl計 算式如下：

式中，

其中，常數ε=1×10?6，r=2。線性權重γL=0.001，γ0=0.998,γR=0.001。光滑因子βl采用經典的WENO方式，定義如下：

2.2 p 階加權WENO 限制器

記網格單元Ij?1、Ij、Ij+1的解為pL(x)、p0(x)、pR(x)。重構問題網格單元為Ij，其內各階解多項式p0,s可通過將原始高階解多項式p0轉換為模態多項式獲得。為保證守恒性且將相鄰網格單元內的線性解多項式投影到問題網格單元，在問題網格單元上可定義來自相鄰網格單元Ij?1的線性多項式如下：

來自相鄰網格單元Ij+1的線性多項式可通過相似方式得到。

最終的WENO 重構解多項式定義為目標單元內的解與相鄰網格單元解的加權組合，保證了和原始解p0有相同的單元平均和k階精度。

上述歸一化非線性權重ωl計算同公式（9），式中ωl計算式如下：

線性權重γl和參數εl與簡單WENO 限制器中的不同，不再是常數，而是因為l對應的階次不同而有區別，設定如下：

式中，參數Kε取0.01。

顯然，線性權重γl和參數εl對 階次的依賴關系弱化了相鄰網格單元Ij±1內的解投影到目標網格單元Ij中的線性多項式，強化了Ij中的高階模態。另外，當低階模態的振蕩程度比高階模態更嚴重時，令低階模態的線性權重為0，去掉振蕩程度更嚴重的部分模態。按照如下平均方式比較兩者的光滑因子：

式中，參數Ktrunc常取1，τs,r是多項式p0,s的光滑因子βs的分量，即因此，計算光滑因子βs時存儲τs,r留待此處使用。方程（19）中右端高階模態的光滑因子平均值對每個s階模態是變化的。令s從最高階k降序至2，對每個s檢查；如果低階模態的光滑因子平均值更大，則令高階模態之外的模態多項式對應的線性權重為0。

2.3 多精度WENO 限制器

記網格單元Ij?1、Ij、Ij+1的解為pL(x)、p0(x)、pR(x)。重構問題網格單元為Ij，其內各階解多項式p0,s可通過將原始高階解多項式p0轉換為模態多項式獲得?；诟麟A解多項式重構得到新的線性解多項式p和非線性解多項式p?。

上述非線性權重ωl1,l2計算如下：

式中，光滑因子β0,1的計算并不依賴于零階多項式p0,0，而是依賴相鄰網格單元中重構得到的多項式：

式中，σL、σR為相鄰網格單元重構解的光滑因子。為保證守恒性，且將相鄰網格單元內的線性解多項式投影到目標單元，在目標單元上可定義來自相鄰網格單元Ij?1的線性多項式如下：

來自相鄰網格單元Ij+1的線性多項式用相似方式得到。

最終的WENO 重構解多項式定義為最高k階非線性解多項式，保證了和原始解p0有相同的單元平均和k階精度。

2.4 保正保精度限制器

盡管前述限制器抑制了FR 方法在激波附近的振蕩，但無法完全避免出現負密度和壓力，限制器對解多項式的重構并不考慮是否會重構出非物理的解。因此還需要保正保精度限制器配合激波捕捉方法。本文實現的保正保精度限制器[11]要求使用Legendre-Gauss-Lobatto 點，以保證采用顯式RK 時間推進格式和一定的CFL（Courant-Friedrichs-Lewy）數下，可以從數學上證明該限制器保正保精度且守恒。記目標網格單元為Ij，k階多項式解分布在N(k)個Legendre-Gauss-Lobatto 點上，記守恒標量qj,i、密度ρ j,i、壓力pj,i，i∈[1,N(k)]，ε=1×10?15表征極小的正值。保正保精度限制器的算法流程如下：

1）限制密度。找出當前網格單元所有解點上密度ρ j,i的最小值，若該值小于允許的最小正值ε，則將該網格單元中所有解點上的密度ρ j,i等比例縮小。

本文所有算例均使用了保正保精度限制器。

3 間斷探測器簡介

3.1 KXRCF 間斷探測器

記網格單元Ij?1、Ij、Ij+1的解為pL(x)、p0(x)、pR(x)。將網格單元Ij的所有面分為兩組：流入面組?I?j和流出面組?I+j，則可由如下條件判斷間斷：

其中，m=1，Ck=1，hj是網格單元Ij的特征半徑。公式（33）的物理含義是，如果該網格單元與相鄰網格單元在“流入面組”上的間斷差值占比超過了Ck值，就判定其為問題網格單元。

3.2 CKXRCF 間斷探測器

記網格單元Ij?1、Ij、Ij+1的解為pL(x)、p0(x)、pR(x)。先計算目標網格單元及其相鄰網格單元內解的積分，再由其差值判斷間斷。具體如下：

其中，Ck=1，hj是網格單元Ij的特征半徑。CKXRCF間斷探測器與KXRCF 間斷探測器的不同在于前者使用了單元內解的積分，而非網格單元交界面通量點處解的積分。一般來說，流場中存在間斷時，高階解多項式在相鄰網格單元中存在振蕩，而且在網格單元內積分的振蕩比網格單元界面處積分的振蕩小，有利于保持穩定地探測包含間斷的網格單元，有利于穩定收斂。這一特性也將在后續算例討論中得到驗證。

4 算例結果及討論

4.1 雙馬赫反射

雙馬赫反射算例是驗證高精度激波捕捉格式的經典算例[18]。控制方程為歐拉方程，計算域為[0,4]×[0,1]。初始條件為下表面x=1/6處一道右行激波Ma=10，與x軸夾角為60°。x=1/6之后下邊界為滑移壁面，上邊界變量隨時間變化與右行激波匹配，右邊界為超聲速出流邊界。計算時間為0～0.2。模擬采用三階空間精度P2 和Roe 通量，基于均勻分布的四邊形網格中Legendre-Gauss-Lobatto 積分點。

圖1～圖3 顯示了SWENO、PWENO 和MWENO限制器P2 在三組網格上的密度云圖結果，圖中藍色線是[1.75, 21.5]范圍內均勻分布的30 條密度等值線，紅色表示t= 0.2 時當地網格單元被KXRCF 間斷探測器標記為問題網格單元。對比顯示PWENO 限制器和MWENO 限制器均能在較細的網格上獲得清晰的渦結構，SWENO 限制器則分辨率較差且數值噪聲較大。SWENO 限制器從相鄰網格單元中直接取多項式替換單元平均，這顯然是其存在較強數值振蕩的直接原因。PWENO 限制器與MWENO 限制器表現類似。PWENO 限制器對渦的分辨率更高，但在上游區域的數值噪聲相對較大。PWENO 限制器與MWENO限制器的基本思路都是將目標網格單元內的高階多項式與相鄰網格單元的低階多項式加權組合，所以表現相似。PWENO 限制器存在可調的參數用于控制數值耗散，而MWENO 限制器基本沒有可調參數。PWENO 限制器在當前參數Ktrunc=1、Kε=0.01下，數值耗散較小，因此對渦的分辨率相對更高，但同時造成數值噪聲更明顯。

圖1 SWENO 限制器P2 在三組網格上的密度云圖Fig. 1 Contour of density on 3 grids using SWENO limiter P2

圖2 PWENO 限制器P2 在三組網格上的密度云圖Fig. 2 Contour of density on 3 grids using PWENO limiter P2

圖3 MWENO 限制器P2 在三組網格上的密度云圖Fig. 3 Contour of density on 3 grids using MWENO limiter P2

4.2 激波與渦相互作用

激波與渦相互作用算例在第五屆高階CFD 國際研討會上用來檢驗高階激波捕捉格式預測強渦與激波相互作用的復雜現象[19]。這是一個二維非定常無黏流內含激波造成的間斷。渦傳播過程中遇到激波后，激波結構會出現顯著畸變，隨之產生多道向下游傳播的波。在這個過程中存在兩種物理現象：第一個流動特征是初始的渦通過激波時由于壓縮效應分裂成兩個渦結構，激波下游的渦結構依賴于激波與渦的相對強度；第二個流動特征是激波下游的波結構反射傳播。本文采用的計算域及設定來自第五屆高階CFD 國際研討會，網格為研討會提供的網格RQ300，細節見第五屆高階CFD 國際研討會[19]。模擬采用四階空間精度P3 和Roe 通量，基于均勻分布的四邊形網格中的Legendre-Gauss-Lobatto 積分點，PWENO 限制器參數Ktrunc=10、Kε=0.01。

在這一算例中，SWENO 中途崩潰，故未顯示結果。圖4 顯示了MWENO 和PWENO 限制器P3 的數值紋影云圖。圖4（a）中紅色顯示了KXRCF 間斷探測器中的參數Ck=0.01時標記的問題網格單元，即KXRCF 間斷探測器只標記了激波，沒有標記激波下游存在數值偽振蕩的區域為問題網格單元。對比圖4（a、b）中的結果顯示，對整個計算域應用PWENO限制器后，激波下游的數值偽振蕩明顯減弱，即PWENO 限制器抑制了FR 格式本身在激波附近的數值偽振蕩。同時也應看到，如圖5 所示，渦核分裂核心區域的流場結構沒有受到數值偽振蕩的明顯影響。圖4（c）顯示了在所有網格單元應用MWENO 限制器后得到的數值紋影。對比圖4（b、c）的結果顯示，相比PWENO 限制器，MWENO 限制器在激波下游區域的數值偽振蕩更加明顯。

圖4 PWENO 和MWENO 限制器P3 的數值紋影云圖Fig. 4 Contour of numerical schlieren using PWENO and MWENO limiter P3

圖5 核心渦核分裂區PWENO 限制器P3 數值紋影云圖Fig. 5 Contour of numerical schlieren in the core region usingPWENO limiter P3

4.3 激波反射

激波反射算例[14]用來測試激波捕捉格式對穩態激波問題的收斂能力。該二維問題的計算域定義為[0,4]×[0,1]。初場設定為(ρ,u,v,p)=(1.0,2.9,0,1.0/1.4)。邊界條件設定如下：計算域下邊界為無黏壁面；右邊界為超聲速出流邊界；左邊界設為Dirichlet 邊界條件，(ρ,u,v,p)=(1.0,2.9,0,1.0/1.4)；上邊界設定為Dirichlet邊界條件，有(ρ,u,v,p)=(1.699 97, 2.619 34, ?0.506 32,1.528 19)。二維數值模擬到收斂，檢查數值解與無黏激波關系式導出的解析解之間的誤差。模擬基于三套連續加密的四邊形網格，采用Lax–Friedrichs 通量；基于四邊形網格中Legendre-Gauss-Lobatto 積分點；PWENO 限制器參數Ktrunc=10、Kε=0.01；時間推進從初場模擬到密度殘差收斂到1×10?12。粗、中、細3 套網格單元數目依次為：120×30、240×60、480×120。

圖6 （a）顯示了細網格P2 采用PWENO 限制器得到的壓力等值線和CKXRCF 間斷探測器標記的紅色問題網格單元，從中可以看到CKXRCF 間斷探測器精準探測了激波所在的網格單元，沒有標記非激波區域。密度殘差收斂后，提取位于y= 0.5 上的壓力。3 套連續加密網格上P2 的壓力分布如圖6（b）所示，精確解根據正激波關系式計算得到。對比連續加密網格上的壓力分布發現，激波下游附近存在振蕩，但隨著網格加密，振蕩迅速變弱，且激波下游存在明顯振蕩的網格數目明顯減少。

圖6 PWENO 限制器P2 在細網格上的壓力等值線圖、y =0.5 處壓力分布和CKXRCF 間斷探測器標記的問題網格單元Fig. 6 Pressure contour and pressure distribution at y = 0.5 using PWENO limiter P2 with contour of troubled cell by CKXRCF shock sensor

圖7 顯示了SWENO、PWENO 和MWENO 限制器在CKXRCF 間斷探測器下的收斂特性及其在y=0.5 處的壓力分布。圖7（a）中PWENO 限制器在不同的間斷探測器下表現不同：采用的KXRCF 間斷探測器無法收斂；采用的CKXRCF 間斷探測器可以收斂；MWENO 限制器配合CKXRCF 限制器也可以收斂；但即使應用CKXRCF 限制器，SWENO 限制器也無法收斂。在定常激波反射這一算例中，圖7（a）中無法收斂的特性沒有對圖7（b）中的y= 0.5 處壓力分布造成明顯影響。SWENO、PWENO、MWENO 限制器在激波附近均存在一定程度的數值振蕩。

圖7 SWENO、PWENO 和MWENO 限制器P2 在中網格上使用不同間斷探測器的密度殘差收斂曲線及y = 0.5 處壓力分布Fig. 7 History of density residual and pressure distribution at y = 0.5 of SWENO, PWENO, MWENO limiter P2 on medium grid with different shock sensors

4.4 黏性弓形激波

弓形激波是驗證激波捕捉格式的經典算例，脫體激波可以盡量隔絕邊界條件的影響來驗證激波捕捉格式, 同時穩態激波用來檢驗高階激波捕捉格式的收斂特性。黏性弓形激波算例在檢驗激波捕捉格式的同時還對數值格式解析高超聲速下很薄的邊界層提出了較高要求。該問題中的幾何構型是一個半圓柱，自由來流Ma=8.03，基于圓柱半徑的雷諾數Re=1.835×105。雖然幾何構型對稱，但仍然保留整體構型以檢驗高階格式是否會產生偽波動。計算域不包括鈍體下游部分，避免出現非定常尾跡區。初場設定為自由來流。邊界條件設定如下：來流邊界為自由來流的 Dirichlet 邊界條件，出流邊界為外插邊界條件，圓柱壁面為黏性絕熱壁面。模擬基于3 套連續加密的四邊形網格，采用三階空間精度P2 和Lax-Friedrichs 通量，基于四邊形網格中Legendre-Gauss-Lobatto 積分點，PWENO 限制器參數Ktrunc=10、Kε=0.01，時間推進從初場模擬到密度殘差收斂到1×10?10。三套網格均設定為：沿壁面均勻分布，在壁面法向方向上自壁面向自由來流邊界拉伸，靠近圓柱壁面附近的拉伸比為1.05，靠近自由來流邊界附近的拉伸比為2。粗、中、細3 套網格單元數目依次為：34×45、68×65、136×95。

在這一算例中，SWENO 中途崩潰，故未顯示結果。圖8 顯示了PWENO 限制器P2 的收斂結果。圖8（a）顯示了在細網格的所有網格單元中應用PWENO 限制器P2 的壓力云圖。壓力進行了歸一化處理：p/p02，其中p02是根據激波關系式計算得到的滯止點處的壓力。3 套連續加密網格中，圓柱壁面上PWENO 限制器P2 的壓力分布如圖8（b）所示。實驗數據提取自 Wieting[20]的實驗，參考數值模擬結果提取自Jiang 等[18]的DG/FV 混合格式結果。結果表明，網格加密后PWENO 限制器P2 在圓柱壁面上的壓力分布收斂到了參考的數值結果。

圖8 PWENO 限制器P2 在細網格上的壓力云圖和3 套網格上的壁面壓力分布Fig. 8 Pressure contour and the converged wall pressure distribution using PWENO limiter P2

圖9 顯示了PWENO 和MWENO 限制器的收斂特性及壁面壓力分布。圖9（a）中PWENO 限制器在不同的間斷探測器下表現不同：采用的KXRCF 間斷探測器和采用的CKXRCF 間斷探測器均無法收斂；全場應用PWENO 限制器可以收斂；但即使全場應用限制器，MWENO 限制器也無法收斂。圖9（a）中無法收斂的特征與圖9（b）中的壁面壓力分布是對應的。PWENO 限制器在KXRCF 和CKXRCF 間斷探測器下無法收斂是中軸線處數值振蕩導致的，MWENO 限制器無法收斂則是出口邊界附近的數值振蕩導致的。

圖9 PWENO 和MWENO 限制器P2 在粗網格上使用不同間斷探測器的密度殘差收斂曲線及壁面壓力分布Fig. 9 History of density residual and wall pressure distribution of PWENO, MWENO limiter P2 on coarse grid with different shock sensors

4.5 激波與層流邊界層相互作用

二維激波與層流邊界層相互作用[21]用以檢驗激波捕捉格式模擬入射激波與平板層流邊界層相互作用。入射激波與層流邊界層相互作用，導致流動分離并產生分離泡。計算域為[?0.2,2]×[0,1]，對層流邊界層的模擬包含前緣。自由來流Ma∞=2.15，基于激波入射點位置的雷諾數Re=ρ∞U∞xsh/μ∞=1×105，激波入射點位于xsh=1處。入射激波與x軸之間的夾角θ= 30.8°。激波入射通過設定邊界條件實現。在y=1.2tanθ處將入口邊界分為上、下兩段。下段y＜1.2tanθ設定為無攻角的自由來流，Ma∞=2.15；上段y＞1.2tanθ設定為帶攻角的來流，其流動參數與自由來流之間滿足Rankine-Hugoniot 關系。其余邊界條件設為：上邊界為遠場自由來流邊界，參考狀態同入口邊界上段；右邊界為外插邊界條件；下邊界上游x＜0部分為對稱面，下游為平板無黏絕熱壁面。初場設定為自由來流Ma∞=2.15的狀態。模擬基于三套連續加密的四邊形網格，采用三階空間精度P2 和Lax–Friedrichs 通量，基于四邊形網格中Legendre-Gauss-Lobatto 積分點，PWENO 限制器參數Ktrunc=10、Kε=0.01，時間推進從初場模擬到密度殘差收斂到1×10?12。三套網格均設定為在壁面法向方向上自壁面向遠場邊界拉伸，靠近壁面附近的拉伸比為1.05，靠近遠場邊界附近的拉伸比為2；在壁面切向方向上自激波入射點向上下游拉伸。粗、中、細3 套網格單元數目依次為：76×32、130×47、154×54。

在這一算例中，SWENO 中途崩潰，故未顯示結果。圖10 顯示了PWENO 限制器P2 的收斂結果。圖10（a）顯示了在細網格的所有網格單元中應用PWENO 限制器P2 的密度云圖，密度進行了歸一化處理：ρ/ρ∞，其中ρ∞是自由來流的密度。三套連續加密網格中，壁面上PWENO 限制器P2 的Cf分布如圖10 （b）所示。參考數值模擬結果提取自第四屆高階CFD 國際研討會上的總結報告[22]，其采用DG 在11041 個網格單元上進行P6 模擬得到收斂解。結果表明，網格加密后PWENO 限制器P2 在壁面上的Cf分布收斂到了參考的數值結果。

圖10 PWENO 限制器P2 在細網格上的密度云圖和3 個網格上的 Cf曲線Fig. 10 Density contour on grid 3 and Cf distribution of PWENO limiter P2 on refined grid

圖11 顯示了PWENO 和MWENO 限制器的收斂特性及其壁面分布。圖11（a）中PWENO 限制器在不同的間斷探測器下表現不同：采用的KXRCF 間斷探測器無法收斂；采用的CKXRCF 間斷探測器逐漸趨向收斂，但存在振蕩；全場應用PWENO 限制器可以收斂；全場應用限制器，MWENO 限制器也可以收斂。圖11（a）中無法收斂的特征對圖11（b）中的壁面分布影響不大，只有CKXRCF 間斷探測器的“振蕩式”收斂對應壁面前半部分的較小偏差。

圖11 PWENO 和MWENO 限制器P2 在中網格上使用不同間斷探測器的密度殘差收斂曲線及壁面壓力分布Fig. 11 History of density residual and wall friction distribution of PWENO, MWENO limiter P2 on medium grid with different shock sensors

5 結 論

本文通過在FR 方法框架下對比研究簡單WENO限制器、p階加權WENO 限制器和多精度WENO 限制器在多個二維無黏及黏性算例中的性能表現。確認簡單WENO 限制器性能較差，多精度WENO 限制器與p階加權WENO 限制器性能相似，捕捉激波的解析精度較高，而且能夠一定程度地抑制FR 格式本身在激波附近的數值偽振蕩。與多精度WENO 限制器相比，p階加權WENO 限制器的穩態收斂性相對更好。雖然本文研究的多個算例中全場應用WENO限制器可以得到收斂的結果，但多精度WENO 限制器和p 階加權WENO 限制器均亟需設計新型的間斷探測器避免限制光滑區域，以減少計算量并保證收斂。