非結(jié)構(gòu)網(wǎng)格通量重構(gòu)算法下三種緊致WENO 限制器對比研究

石京昶，嚴 紅，*

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

0 引 言

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

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

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

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

1 FR 方法簡介

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

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

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

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

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

式中，αj,f,l是獨立于解的常數(shù)，Sf是單元界面面積，j是解點在當前網(wǎng)格單元內(nèi)的編號，f是當前網(wǎng)格單元面的編號，l是通量點在當前面的編號。綜上所述，F(xiàn)R 方法由以下離散方程近似求解原雙曲律方程：

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

2 WENO 限制器簡介

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

2.1 簡單WENO 限制器

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

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

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

式中，單元平均定義為：

注意以上積分均為目標單元Ij內(nèi)的積分。

上述歸一化非線性權(quán)重ωl計 算式如下：

式中，

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

2.2 p 階加權(quán)WENO 限制器

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

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

最終的WENO 重構(gòu)解多項式定義為目標單元內(nèi)的解與相鄰網(wǎng)格單元解的加權(quán)組合，保證了和原始解p0有相同的單元平均和k階精度。

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

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

式中，參數(shù)Kε取0.01。

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

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

2.3 多精度WENO 限制器

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

上述非線性權(quán)重ωl1,l2計算如下：

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

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

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

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

2.4 保正保精度限制器

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

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

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

3 間斷探測器簡介

3.1 KXRCF 間斷探測器

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

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

3.2 CKXRCF 間斷探測器

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

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

4 算例結(jié)果及討論

4.1 雙馬赫反射

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

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

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

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

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

4.2 激波與渦相互作用

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

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

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

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

4.3 激波反射

激波反射算例[14]用來測試激波捕捉格式對穩(wěn)態(tài)激波問題的收斂能力。該二維問題的計算域定義為[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)。二維數(shù)值模擬到收斂，檢查數(shù)值解與無黏激波關(guān)系式導出的解析解之間的誤差。模擬基于三套連續(xù)加密的四邊形網(wǎng)格，采用Lax–Friedrichs 通量；基于四邊形網(wǎng)格中Legendre-Gauss-Lobatto 積分點；PWENO 限制器參數(shù)Ktrunc=10、Kε=0.01；時間推進從初場模擬到密度殘差收斂到1×10?12。粗、中、細3 套網(wǎng)格單元數(shù)目依次為：120×30、240×60、480×120。

圖6 （a）顯示了細網(wǎng)格P2 采用PWENO 限制器得到的壓力等值線和CKXRCF 間斷探測器標記的紅色問題網(wǎng)格單元，從中可以看到CKXRCF 間斷探測器精準探測了激波所在的網(wǎng)格單元，沒有標記非激波區(qū)域。密度殘差收斂后，提取位于y= 0.5 上的壓力。3 套連續(xù)加密網(wǎng)格上P2 的壓力分布如圖6（b）所示，精確解根據(jù)正激波關(guān)系式計算得到。對比連續(xù)加密網(wǎng)格上的壓力分布發(fā)現(xiàn)，激波下游附近存在振蕩，但隨著網(wǎng)格加密，振蕩迅速變?nèi)酰壹げㄏ掠未嬖诿黠@振蕩的網(wǎng)格數(shù)目明顯減少。

圖6 PWENO 限制器P2 在細網(wǎng)格上的壓力等值線圖、y =0.5 處壓力分布和CKXRCF 間斷探測器標記的問題網(wǎng)格單元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 限制器在不同的間斷探測器下表現(xiàn)不同：采用的KXRCF 間斷探測器無法收斂；采用的CKXRCF 間斷探測器可以收斂；MWENO 限制器配合CKXRCF 限制器也可以收斂；但即使應用CKXRCF 限制器，SWENO 限制器也無法收斂。在定常激波反射這一算例中，圖7（a）中無法收斂的特性沒有對圖7（b）中的y= 0.5 處壓力分布造成明顯影響。SWENO、PWENO、MWENO 限制器在激波附近均存在一定程度的數(shù)值振蕩。

圖7 SWENO、PWENO 和MWENO 限制器P2 在中網(wǎng)格上使用不同間斷探測器的密度殘差收斂曲線及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 黏性弓形激波

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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