石國權(quán)，燕振國，朱華君，*，馬燕凱，賈斐然，2

(1.空氣動力學國家重點實驗室，綿陽 621000；2.西北工業(yè)大學 動力與能源學院，西安 710000)

0 引言

近年來，計算流體力學（Computational Fluid Dynamics，CFD）快速發(fā)展，不管是基礎(chǔ)理論研究還是工程實際應(yīng)用都需要進行精細的流場結(jié)構(gòu)模擬。對于旋渦、分離、湍流等復雜流動現(xiàn)象，低階格式數(shù)值耗散和色散較大，難以給出精細的流場結(jié)構(gòu)（計算代價遠遠大于高階精度格式），需采用耗散和色散小的高階精度計算格式（3 階及以上）[1]。另一方面，對復雜構(gòu)型有著良好適應(yīng)性的非結(jié)構(gòu)網(wǎng)格，生成簡單，且隨機的數(shù)據(jù)結(jié)構(gòu)非常利于網(wǎng)格自適應(yīng)，非結(jié)構(gòu)網(wǎng)格算法受到國內(nèi)外學者重視[2]。近年來發(fā)展的高階CPR方法是一種既適用于結(jié)構(gòu)網(wǎng)格[3]，也適用于非結(jié)構(gòu)網(wǎng)格[4]的緊致的差分型格式，具有很高的計算效率。

CPR 方法最早是由Huynh[3]在2007 年為求解結(jié)構(gòu)網(wǎng)格下的雙曲守恒律方程提出的，當時被稱為通量重構(gòu)（flux reconstruction，F(xiàn)R）方法。一維FR 方法可以直接利用張量計算推廣到四邊形網(wǎng)格和六面體網(wǎng)格上。2009 年，Wang 和Gao[4]將FR 的概念推廣到了二維三角形網(wǎng)格及混合網(wǎng)格上，提出了提升配點懲罰(lifting collocation penalty,LCP)方法。由于FR 和LCP方法緊密聯(lián)系，相關(guān)學者將其統(tǒng)稱為CPR 方法[5]。

CPR 方法在選擇特殊的求解點和修正函數(shù)時，在線性方程情況下可以等價于間斷Garlerkin（discontinuous Galerkin，DG）方法、譜差分（spectral difference，SD）方法等，同時呈現(xiàn)出更簡化的形式，計算量又相對較小[3]。具體來看，DG 方法在求解過程中涉及到通量函數(shù)的積分，計算量較大，而CPR 方法是基于求解點的差分型方法，沒有積分計算的過程，因而計算方便且復雜度較低。此外，由于CPR 方法在具體計算某個單元內(nèi)的函數(shù)值時，只需要利用與其直接相鄰單元的函數(shù)值，因此CPR 格式具有緊致的特點，方便進行并行計算，適合求解大規(guī)模復雜構(gòu)型流動問題。

但CPR 方法是高階線性格式，在捕捉較強間斷時，容易產(chǎn)生明顯的虛假振蕩，需要發(fā)展適合非結(jié)構(gòu)CPR 方法的激波捕捉策略。目前對間斷問題，常見的限制技術(shù)有添加人工黏性[6-8]、添加限制器[9-11]、采用混合格式[12-13]等方法。這些常見的限制技術(shù)已有學者將其推廣應(yīng)用到CPR 方法中，如MLP 方法[14]、人工黏性方法[15]、斜率限制器[16]、WENO 限制器方法[17]等。值得一提的是子單元限制技術(shù)的發(fā)展。Persson 和Peraire[7]提出，對高階多項式近似，分辨率為 δ～h/p（h是 單元尺度，p是多項式階數(shù)），激波通?？缭?～4 個單元，而劃分子單元的做法可以將激波限制在子單元層面。如高階DG 方法在單元內(nèi)構(gòu)造了高次多項式分布，具有高分辨特性，子單元限制技術(shù)可以在激波捕捉過程中盡可能地保持這種高分辨特性。Huerta 等[18]提出的用于DG 單元內(nèi)分段積分常數(shù)的激波捕捉方法、Dumbser 等[19-20]基于多維最優(yōu)階 偵 測（multi-dimensional optimal order detection，MOOD）方法提出的后驗子單元有限體積限制器、Kochi 和Ramakrishna[21]提出的DG 上的緊致子單元加權(quán)本質(zhì)無振蕩（compact subcell weighted essentially non-oscillatory，CSWENO）限制方法、Guo 等[22]提出的用于守恒律方程激波捕捉的三階WCNS-CPR 混合格式，都是基于子單元的思想進行子單元激波捕捉格式的構(gòu)造。但是這些方法中子單元等距分布，需要考慮主單元與子單元之間的數(shù)據(jù)交換。Hennemann 等[23]提出可壓縮Euler 方程的高階分裂形式DG 的熵穩(wěn)定子單元激波捕捉方法，基于Gauss-Legendre-Lobatto（GLL）求解點，設(shè)計了一種非等距的子單元限制方法。Sonntag 和Munz[24]也提出了一種DG 方法的FV子單元限制方法，子單元劃分方式選取了Gauss積分權(quán)值作為子單元的邊長。非等距剖分下的子單元限制過程中，子單元的值直接由主單元上求解點的值得到，避免了平均劃分子單元時所遇到的主單元與子單元之間的數(shù)據(jù)傳遞問題。這些子單元限制方法的基本思想是偵測到問題單元后，將其剖分為子單元，然后在子單元上采取具有激波捕捉能力的格式。從這個角度來看，子單元限制屬于混合格式的范疇。

我們前期針對二維結(jié)構(gòu)網(wǎng)格提出了一種高階CPR 方法的子單元限制方法[25]并推廣到非結(jié)構(gòu)四邊形網(wǎng)格[26]，在光滑單元上采用CPR 方法，問題單元上采用二階非等距非線性加權(quán)激波捕捉格式。該子單元限制方法使用擴展多項式的最高模態(tài)衰減（highest modal decay of the extended polynomials,MDHE）激波偵測因子[23]偵測問題單元。本文在文獻[25]的基礎(chǔ)上，基于非結(jié)構(gòu)四邊形網(wǎng)格CPR 方法的子單元限制方法，討論了單元混合策略和分維混合的可行性，并對比了不同問題單元偵測因子下的結(jié)果，為間斷有限元類高階格式（DG、CPR、SD 等）的子單元限制提供可靠的參考。

1 控制方程及離散

1.1 控制方程

考慮守恒形式的二維Euler 方程

式中，U是守恒變量，F(xiàn)、G是通量。

式中，ρ 為密度，u、v為速度，p為 壓力，E為總能。對理想氣體，比熱比 γ=1.4。

1.2 二維非結(jié)構(gòu)四邊形網(wǎng)格上的CPR 方法

首先，建立物理空間中的四邊形網(wǎng)格單元與計算空間中的標準單元I=[-1,1]×[-1,1]之間的坐標變換關(guān)系：

式中，L是用來定義四邊形的點的個數(shù)，在直邊四邊形中L=4對 應(yīng)4 個頂點，(xl,yl)是這些點的笛卡爾坐標，Ml是形狀函數(shù)。Jacobi 矩陣形式如下：

假設(shè)該矩陣非奇異，則其逆矩陣為：

式中，

將Euler 方程（1）變換到計算空間后，變?yōu)椋?/p>

每個單元上的通量值利用單元內(nèi)部K×K個求解點處的通量值構(gòu)造多項式得到。二維情況下，通量的插值多項式為：

式中，Lk是Lagrange 多項式。此時的通量多項式函數(shù)是間斷通量函數(shù)，單元界面通量不連續(xù)，會破壞守恒律，需要構(gòu)造連續(xù)通量多項式函數(shù)，即：

式中，gL、gR是K次多項式修正函數(shù)，滿足：

在CPR 的計算過程中，需要求解通量導數(shù)。Huynh[3]基于通量的Lagrange 插值多項式和Riemann通量對間斷通量函數(shù)進行修正，將界面處間斷通量函數(shù)與Riemann 數(shù)值通量的差量引入到單元內(nèi)部，構(gòu)造連續(xù)的通量函數(shù)，建立了單元間的聯(lián)系。求解點處的通量導數(shù)可由連續(xù)通量函數(shù)式（14）直接求導。最終得到計算空間中的空間離散方程：

本文采取五階CPR 格式，求解點選擇五點Gauss-Legendre 點，單元界面處的Riemann 通量選擇局部 Lax-Friedrichs 通量。Riemann 通量基于原始變量計算，原始變量在界面的左右值通過基于求解點的插值多項式得到，多項式的形式與式（12、13）相同。修正函數(shù)選擇文獻[3]中g(shù)DG修正函數(shù)：

式中，PK是Legendre 多項式。

時間離散格式采取三階Runge-Kutta 格式，用RHS表示式（1）的右端項：

則三階Runge-Kutta 格式為：

式中，j是單元序號，k、l是四邊形單元內(nèi)部求解點序號，dt是時間步長。

2 子單元限制策略

本文子單元限制策略的主要思想是，利用問題單元偵測因子標記流場中含有大梯度或是間斷的不光滑單元為問題單元，然后將這些問題單元按照求解點位置劃分為非等距子單元，在子單元上添加限制措施。

2.1 問題單元偵測因子

本文的重點在于介紹不同偵測因子對結(jié)果的影響，選擇了目前較為流行的TVB[12,27]、MDH[23,25-26]、JST[24,28]三種問題單元偵測因子進行對比。下面簡單介紹這三種偵測因子的基本方法。

2.1.1 TVB 偵測因子

總變差有界（total variation bounded,TVB）是有限體積中經(jīng)典的minmod 型限制器[27]。Qiu 等[12]基于TVB 構(gòu)造了問題單元偵測因子，每當minmod 限制器試圖改變斜率時，單元就被標記為問題單元。

TVB 限制器的思想是：

式中，

若(a1,a2,···,ak)沒 有返回a1值，則說明此單元是問題單元，定義問題單元指示器 β=1，否則為光滑單元，定義 β=0。式（24）中的M是人工定義的常量，通過改變M的大小可以調(diào)節(jié)問題單元判斷條件的嚴格程度。

2.1.2 MDHE 偵測因子

正如文獻[7]中提到的那樣，解的表達多項式的最高模態(tài)能量系數(shù)在光滑區(qū)域會迅速衰減，而在非光滑區(qū)域則衰減較慢，表征了單元的光滑程度。因此，可以針對最高模態(tài)定義閾值，來判斷單元內(nèi)是否存在間斷。該方法被稱作最高模態(tài)衰減（highest modal decay,MDH）方法[23]。Hennemann 等使用最高和次高模態(tài)系數(shù)改進了此方法以避免奇偶效應(yīng)[29]，但是只利用了單元內(nèi)部求解點，對于存在于單元界面處的間斷可能無法識別。我們前期通過引入左右界面處Roe 平均值構(gòu)造新的多項式[25]，構(gòu)造新的多項式所用模板為 {u0,u1,···,uK+1}，其中u0和uK+1為左右界面處Roe 平均值。Roe 平均值基于界面通量點附近的兩個求解點值計算得到[30]。由于它是通過擴展MDH 偵測因子的多項式階來構(gòu)造的，所以它被命名為MDHE[26]。一維多項式的模態(tài)能量定義為：

式中，K是單元內(nèi)求解點個數(shù)，是模態(tài)系數(shù)，P j（j=0,···,K+1）是 Legendre 基函數(shù)，取 ?=ρp為偵測變量。解的多項式可以表達為Lagrange 多項式展開式：

也可以表達為Legendre 基函數(shù)展開式：

基于最高和次高模態(tài)能量的占比定義為：

定義閾值為：

按照大量算例測試的經(jīng)驗[23]，式（30）中取a=0.5、c=1.8 。如果IE＞IT，則單元被標記為問題單元，并定義 β=1；否則判定為光滑單元，β=0。

2.1.3 JST 偵測因子

JST 偵測因子源于Jameson 等[28]提出的經(jīng)典的Jameson-Schmidt-Turkel（JST）方法。該方法根據(jù)基于壓力梯度的開關(guān)函數(shù)將人工二階和四階耗散項添加到歐拉方程以抑制振蕩。Sonntag 等[24]在研究DG 方法的子單元限制策略中，改進了JST 偵測因子。本文采用該改進形式的JST 偵測因子。

一維情況下：

式中，k代表單元內(nèi)部求解點的局部坐標索引，pmin,k=min(pk±d,d=1,2,3)是求解點附近節(jié)點值的最小壓力，pmax,k=max(pk±d,d=1,2,3)是求解點附近節(jié)點值的最大壓力。

二維情況下考慮兩個方向：

式中，k、l代表單元內(nèi)部求解點的局部坐標索引，pmin,k,l=min(pk±d,l±d,d=1,2,3)，pmax,k,l=max(pk±d,l±d,d=1,2,3)。在非結(jié)構(gòu)四邊形網(wǎng)格中，模板由本單元及共邊的相鄰單元提供。式（32）給出了每個子單元或CPR 求解點處的指示值，使用體積加權(quán)，得到單元的單個指標值：

式中，K是單元內(nèi) ξ 或 η方 向的求解點個數(shù)，VElement和Vk,l在標準單元上計算，分別是標準單元和子單元的面積。

如不特殊說明，我們?nèi)￠撝礗T=0.01來判斷間斷的存在。當IJST＞IT時，認為該單元是問題單元，記為β=1。另外，偵測變量的選取與文獻[24]一致，使用壓力p作為偵測變量。

2.2 基于非等距非線性插值的子單元激波捕捉格式（CNNW2）

在子單元上采取常值分布函數(shù)可以構(gòu)造類似于Godunov 格式的一階迎風格式，但是一階格式的耗散太大，分辨率會嚴重下降。為提高子單元上的分辨率，在前期的工作中[26]，我們將Zhu 和Liu 等[24]提出的用作五階CPR 方法的子單元限制器的二階緊致非等距非線性加權(quán)格式（CNNW2）推廣到了非結(jié)構(gòu)網(wǎng)格。由于該子單元限制方法也是一種局部混合格式，故稱作CPR-CNNW2 格式。

這里我們簡要介紹一下CNNW2，詳細可見文獻[26]。非等距子單元的劃分如圖1 所示，子單元的求解點與CPR 單元的求解點重合[23]（紅色圓點是Gauss 求解點，視作子單元“中心”；藍色方塊是子單元通量點，視作子單元邊界），避免了子單元等距剖分下，主單元與子單元之間的數(shù)據(jù)交換的麻煩。每個子單元的長度可以表示為ξfpi+1-ξfpi=ωi，（i=1,2,…,K），其中：ωi為 Gauss 積分權(quán)；ξfpi為通量點位置，ξf p1=-1。

圖1 子單元劃分方式Fig.1 Layout of subcell division

對于五階CPR，考慮單元五個Gauss 求解點的中間三個求解點對應(yīng)的子單元，以第二個求解點為例，如圖2 所示，CNNW2 插值利用標準單元內(nèi)部三個求解點得到通量點處的值。

圖2 二階非等距非線性插值模板Fig.2 Stencil for the second-order nonuniform nonlinear interpolation

具體求解步驟如下：

4）添加限制器控制數(shù)值振蕩，得到最終的子單元界面處物理量的逼近值。

式中限制器為：

式中，m=min(u1,u2,u3),M=max(u1,u2,u3)。

再考慮非結(jié)構(gòu)四邊形單元，每個方向上界面附近的兩個求解點u1和u5對 應(yīng)的子單元。圖3 是求解u1對應(yīng)子單元的情況，插值過程只在求解界面處時涉及到單元之間的數(shù)據(jù)處理。將的計算由計算空間改為物理空間上進行（反距離加權(quán)代入物理空間距離），在物理空間上計算得到。

圖3 物理空間上插值模板Fig.3 Interpolation stencil in the physical space

插值過程使用特征變量，根據(jù)以上插值方法得到子單元界面的左右值，特征投影的過程參考文獻[30]。最終得到子單元界面的原始變量，進而計算子單元界面處的Riemann 通量。半離散形式使用二階差分算子：

問題單元內(nèi)部的子單元分布是結(jié)構(gòu)化的，由一維空間離散格式可以直接推廣得到二維空間離散格式。

2.3 CPR-CNNW2 混合策略

2.3.1 格式切換

一個單元所采用的格式由問題單元指示器 β決定，當 β=1時 采用CNNW2，當 β=0時采用CPR。問題單元指示器的值發(fā)生變化，格式對應(yīng)地進行切換。單元與單元的交界面處需要計算共用的Riemann 通量。Riemann 通量的左值和右值的求解方式根據(jù)單元是否為光滑單元而定，對于CPR 單元和CNNW2 單元的界面，假設(shè)左側(cè)單元為CPR 單元，右側(cè)為CNNW2 單元。左值由光滑單元提供，使用CPR 中的Lagrange 插值獲得；右值由問題單元提供，使用二階非等距非線性插值獲得，然后得到Riemann 通量為：

也就是說，從求解點到界面通量點的插值方法由單元的類型進行控制。

2.3.2 單元/分維混合CPR-CNNW2 格式

四邊形網(wǎng)格中，問題單元的偵測是逐維進行的，當單元內(nèi)任一行/列被標記為問題單元，則整個單元都被標記為問題單元。按照問題單元分布，以整個單元為基本單元進行CPR-CNNW2 格式的切換計算，稱作單元混合CPR-CNNW2 格式。

單元混合CPR-CNNW2 格式是常見的一種混合方法，但在四邊形網(wǎng)格CPR 方法中，計算是逐維進行的。在單元內(nèi)部只標記出現(xiàn)問題的行或列，如二維問題中激波附近區(qū)域的求解點的空間離散，在一個方向使用CNNW2 格式，另一個方向使用CPR 的模式是可行的，我們稱之為分維混合CPR-CNNW2 格式，如圖4。特定情況下，采用分維偵測可以減少偵測單元，提高格式分辨率，同時提高計算效率。我們在算例測試中，測試了分維混合CPR-CNNW2 格式下不同偵測因子的表現(xiàn)。

圖4 分維混合CPR-CNNW2 示意圖Fig.4 Schematic of the hybrid CPR-CNNW2 scheme by dimension

2.4 子單元時間步長限制

CPR 和DG 方法的CFL 數(shù)是和多項式階數(shù)有關(guān)的，參照Dumbser 在文獻[19]中提到的時間步長限制：

式中，d是空間維度，N是多項式階數(shù)，h代表網(wǎng)格尺寸，|λmax|是最大特征波速度。對二維Euler 方程，其特征矩陣A(U)=?F/?U的特征值λA=[u-c,u,u,u+c]，特征矩陣B(U)=?G/?U的特征值 λB=[v-c,v,v,v+c]。令

式中，Δdj是相鄰單元重心距離最小值。Dumbser[19]利用CFL 條件得到DG 的時間步長限制下，等距子單元的最大劃分個數(shù)為Ns=2N+1，其中N是多項式次數(shù)。當Ns＞2N+1時，子單元時間步長限制會反過來約束DG 的最大時間步長。本文的非等距劃分子單元方式，標準單元內(nèi)子單元尺寸符合Gauss 積分權(quán)分布，和求解點個數(shù)相關(guān)，最小的為第一積分權(quán) ω1代表的子單元長度。CPR 格式使用五階精度，Guass 積分使用五點格式時，ω1＞1/(2N+1)。故在此可以統(tǒng)一使用CPR 的時間步長限制。

3 算例測試

3.1 Sod 激波管問題

Sod 激波管問題的初始條件為：

采用單元數(shù)為40，左右邊界為Dirichlet 邊界條件，計算終止時間為T=0.2，CFL=0.5，對三種激波偵測因子下的CPR-CNNW2 的計算結(jié)果進行對比。在此算例中TVB 偵測因子的自由參數(shù)取M=100.0。

計算結(jié)果如圖5 所示。為了方便對比三種激波偵測因子下的問題單元分布，將圖5 右側(cè)縱坐標軸設(shè)為問題單元指示器 β，并將三種方法錯位分布。在接觸間斷x=0.7 和x=0.85 附近，使用三種問題單元偵測因子都無明顯振蕩。對比之下，無論采取何種偵測方式，混合CPR-CNNW2 格式都比CNNW2 格式對間斷的抹平小。雖然三種偵測因子最終的偵測結(jié)果中，問題單元只存在于間斷附近，但可以看到采取TVB方法對接觸間斷的抹平較大，而MDHE 和JST 方法對間斷的抹平較小，數(shù)值表現(xiàn)更好。

圖5 Sod 激波管的密度分布和問題單元分布圖Fig.5 Distributions of the density and problematic cells for the Sod shock tube problem

3.2 Shu-Osher 問題

Shu-Osher 問題的初始條件為：

采用單元數(shù)為100（對應(yīng)DoFs=500），左右邊界分別為流入流出邊界條件，計算終止時間為T=0.18，CFL=0.9 。在此算例中TVB 偵測因子的自由參數(shù)取M=3.0。

計算結(jié)果見圖6，其中參考結(jié)果為混合CPR-CNNW2格式在x方向DoFs=2 000 下的計算結(jié)果?；旌细袷皆趩栴}單元采取低階格式計算，在流場的大部分光滑區(qū)域采取高階CPR 方法計算，能夠顯著提高算法的分辨率。通過比較得到，MDHE 偵測因子下的分辨率最高，JST 次之，TVB 最低，相對應(yīng)的問題單元占總單元數(shù)分別為5%、9%、23%。這反映了問題單元偵測區(qū)域與分辨率之間的關(guān)系：問題單元越多，即采取二階格式計算的單元越多，引入的耗散越大，分辨率越低。問題區(qū)域的分布隨時間變化見圖7，TVB 偵測標記的問題單元較多，尤其是在高頻振蕩區(qū)，無法有效識別問題單元，這也導致了在高頻振蕩區(qū)TVB 偵測因子下混合CPR-CNNW2 格式的結(jié)果幾乎和CNNW2的結(jié)果重合，而MDHE 和JST 的結(jié)果具有更高的分辨率，表現(xiàn)較好。

圖6 Shu-Osher 問題的密度分布和問題單元分布Fig.6 Distributions of the density and prolematic cells for the Shu-Osher problem

圖7 問題單元隨時間變化的分布圖Fig.7 Distributions of problematic cells over time

另外，我們將本文的CPR-CNNW2 結(jié)果與文獻[31]中的p-weighted 限制器在取p=5次多項式時采用Nc=100（DoFs=600）網(wǎng)格下的計算結(jié)果進行對比。可以看到在高頻振蕩區(qū)，MDHE 和JST 偵測因子下的計算結(jié)果的分辨率比p-weighted 限制器的計算結(jié)果稍好。

3.3 等熵渦測試

二維等熵渦問題的初始條件是在一個均勻流上添加一個各向異性的等熵旋渦的擾動。初始條件設(shè)置見文獻[32]，取自由度為200×200，時間步長dt=0.000 1，計算到T=1.0。基于等熵渦問題，對CNNW2與CPR 的計算時間進行對比測試。

我們通過調(diào)節(jié)TVB 偵測因子中的參數(shù)M減小閾值，將一些光滑單元設(shè)置為問題單元，對比單元混合策略和分維混合策略的計算時間，如表1 所示。在相同的偵測閾值和計算條件下，分維混合策略的計算時間比單元混合花費時間略少，計算效率更高。但是在單元混合下，四邊形單元出現(xiàn)不滿足條件的維度時便可定義為問題單元，分維混合下需要計算所有維度。故在偵測時間上，分維混合所花費的時間較多。

表1 單元混合與分維混合的計算時間Table 1 Calculation time for hybrid by cell and by dimension

另外，調(diào)整不同偵測因子的自由參數(shù)使偵測到的問題單元占比為0，以測試不同偵測因子對計算時間的影響。由表2 中計算時間可知，在相同計算網(wǎng)格下，CPR 計算效率比CNNW2 高。另外，對比分維混合CPR-CNNW2 格式和單個CPR 格式，偵測因子會帶來額外的計算量。表2 中單個CPR 的計算時間和偵測時間的總和與混合格式的總計算時間大致相同（相對誤差小于1%），觀察偵測時間的影響，MDHE偵測因子占總計算時間的6.71%，TVB 和JST 偵測因子計算開銷更多。問題單元的占比對計算時間有一定影響，問題單元越多，采取CNNW2 格式計算的單元越多，計算效率越低。同時，偵測因子本身的計算時間也會對總計算時間產(chǎn)生影響，計算時間最小的是MDHE 偵測因子。

表2 不同偵測因子對計算時間的影響Table 2 Influence of different indicators on the calculation time

3.4 雙馬赫反射

雙馬赫反射的初始條件為：

要提到的是，TVB 偵測因子魯棒性比較差，為解決這個問題，我們引入了過渡單元[22]，把問題單元的共邊相鄰單元也標記為問題單元，增加了偵測區(qū)域。雙馬赫反射問題計算結(jié)果見圖8、圖9。從圖8 可以看出，三種偵測因子下CPR-CNNW2 都能準確地捕捉到馬赫桿。JST 偵測下的計算結(jié)果中剪切渦更明顯，分辨率更高，TVB 偵測與MDHE 偵測下的計算結(jié)果相差不大。對比問題單元的分布，TVB、MDHE 和JST 偵測因子標記的問題單元占總單元數(shù)分別為22.66%、3.84%、3.71%。TVB 偵測因子標記的問題單元最多，MDHE 明顯少了很多，但是密度等值線結(jié)果的分辨率相差不多。JST 偵測因子下剪切渦附近的問題單元最少，所以對剪切渦的影響最小，分辨率最高。

圖8 CPR-CNNW2 在不同偵測因子下求解雙馬赫反射問題（從1.5 到21.7 的30 條密度等值線）Fig.8 Density contours with 30 levels ranging from 1.5 to 21.7 for CPR-CNNW2 with different indicators in the double Mach reflection problem

圖9 雙馬赫反射問題的各種偵測因子下的問題單元分布圖Fig.9 Distribution of problematic cells under different indicators for the double Mach reflection problem

從對比結(jié)果中可以看出，不同的問題單元偵測因子，偵測問題單元的區(qū)域不同，對計算結(jié)果的影響不同。問題單元越多，計算時間越長，分辨率越低。

由于TVB 偵測因子在分維混合計算雙馬赫問題，即使加上過渡單元且參數(shù)M取到1 ×10-6，依然計算失敗，故以MDHE 和JST 偵測因子測試分維混合下的CPR-CNNW2 方法。雙馬赫問題計算結(jié)果密度等值線見圖10（a、b）。圖10（c、d）中紅色區(qū)域內(nèi)的求解點 ξ、η兩個方向都采取CNNW2 格式，綠色部分只有一個方向采取CNNW2 格式計算。分維混合計算會使魯棒性下降，我們按照一維形式增加了過渡單元。其中JST 偵測器在閾值不變的情況下計算穩(wěn)定，而MDHE 偵測器需要減小閾值才能穩(wěn)定計算。從圖10（c、d）可以看到，過渡單元的添加和閾值的降低有可能會造成分維偵測的問題單元多于單元偵測的問題單元。在魯棒性上，JST 偵測器的表現(xiàn)要好于MDHE 偵測器，但是JST 偵測因子下的計算結(jié)果仍需關(guān)注邊界處出現(xiàn)局部振蕩的現(xiàn)象。分維混合計算需要結(jié)合適合的問題單元偵測因子，才能更好地發(fā)揮減少問題單元區(qū)域、提高計算結(jié)果分辨率的優(yōu)勢。

圖10 雙馬赫反射的分維混合CPR-CNNW2的計算結(jié)果Fig.10 Numerical results of CPR-CNNW2 with detecting by dimension for the double Mach reflection problem

3.5 激波-旋渦干擾

激波-旋渦干擾問題初始條件設(shè)置可以參考文獻[33]。在該問題中，運動的旋渦和靜止激波相互干擾，發(fā)展出具有平滑特征和間斷的復雜流動結(jié)構(gòu)，常用來測試格式捕捉激波和維持分辨率的能力。計算域[0,2]×[0,1]，初始條件由位于x=0.5的靜止激波Ms和 中心位于 (xc,yc)=(0.25,0.5)的 旋渦Mv給出。其中旋渦Mv的切向速度表達式為：

式中，r2=(x-xc)2+(y-yc)2，vm是最大速度。

采用 120×60 (DoFs=600×300)的均勻網(wǎng)格計算該問題。左右邊界為Dirichlet 邊界條件，上下邊界設(shè)置為壁面滑移邊界條件，計算終止時間為T=0.7，CFL=0.9 。本算例中，TVB 偵測因子參數(shù)M=100.0，JST 偵測因子閾值IT=0.02。

3.5.1 單元混合CPR-CNNW2 格式計算結(jié)果

計算結(jié)果如圖11、圖12。從圖11 中，明顯觀察到TVB 偵測因子下的結(jié)果在旋渦附近的分辨率明顯不如MDHE 和JST 偵測下的計算結(jié)果，但相對于單一CNNW2 格式計算結(jié)果，分辨率改善很多。圖12 中的TVB、MDHE、JST 三種偵測因子標記的問題單元區(qū)域占比依次為3.01%、1.72%、1.67%，且TVB 在旋渦處標記了較多問題單元，對耗散的引入比較大。

圖11 不同偵測因子下CPR-CNNW2 求解激波-旋渦干擾問題的計算結(jié)果Fig.11 Numerical results of CPR-CNNW2 with different indicators for the shock-vortex interaction problem

圖12 激波-旋渦干擾問題的問題單元分布Fig.12 Distribution of problematic cells for the shock-vortex interaction problem

為了更精細地對比，在不同結(jié)構(gòu)特征處，通過截線進行定量分析[34]，兩條截線分別為截線①（y=0.40）和截線②（x=1.05），如圖13 所示。CNNW2 格式和三種偵測因子下的單元混合CPR-CNNW2 格式在截線處的密度分布結(jié)果見圖14。從圖中看出，TVB 偵測因子下旋渦處的耗散遠大于JST 和MDHE 偵測因子下的結(jié)果。其次，從問題單元占比隨時間變化的情況（見圖15）可以看出，整個求解過程中TVB 的問題單元最多，MDHE 次之，JST 最少。問題單元偵測區(qū)域的大小不僅和偵測因子類型相關(guān)，也和閾值的設(shè)置相關(guān)，如圖15 中黑色實線為JST 偵測因子閾值IT=0.01時問題單元的時間歷史，旋渦區(qū)域也被標記為問題單元，引入耗散過大。MDHE 和JST 表現(xiàn)的微小差距和閾值的設(shè)置也是有關(guān)的。

圖13 激波-旋渦干擾問題的參考線位置Fig.13 Reference lines for the shock-vortex interaction problem

圖14 不同偵測因子下單元混合CPR-CNNW2 格式在兩條截線上的密度分布圖Fig.14 Density distribution along the two sliced lines of CPRCNNW2 with detecting by cell for different cell indicators

圖15 單元混合下CPR-CNNW2 格式問題單元占比隨時間變化圖Fig.15 Time variation of the problematic cell proportion for CPR-CNNW2 with detecting by cell

3.5.2 分維混合CPR-CNNW2 格式計算結(jié)果

圖16 中展現(xiàn)了三種偵測因子下分維混合CPRCNNW2 方法在截線處的密度分布結(jié)果。在當前閾值下，MDHE 分維偵測下的計算結(jié)果分辨率在三者之中最高。另外從圖17 可以看到，分維混合CPR-CNNW2方法下的問題單元分布和激波的分布相關(guān)。

圖16 分維混合CPR-CNNW2 格式在兩條截線上的密度分布Fig.16 Density distribution along two sliced lines for CPR-CNNW2 with detecting by dimension

圖17 分維混合CPR-CNNW2 計算激波-旋渦干擾問題的問題單元分布Fig.17 Distribution of problematic cells of CPR-CNNW2 with detecting by dimension for the shock-vortex interaction problem

觀察圖18 中分維混合CPR-CNNW2 格式分別采取不同偵測因子的問題單元占比隨時間變化情況，整個計算過程中，ξ 方向和 η方向的問題單元的占比相差很大，說明單元偵測的方式會標記“多余”部分。另外，最終的分維混合格式計算結(jié)果，MDHE 比JST 分辨率略高，可能是由于在T=0.1到T=0.4激波和旋渦干擾期間，MDHE 分維偵測比JST 分維偵測標記的問題單元少。

圖18 分維混合CPR-CNNW2 格式的問題單元占比隨時間變化圖Fig.18 Time variation of the problematic cell proportion for CPR-CNNW2 with detecting by dimension

基于閾值對計算結(jié)果的影響，我們對比了三種偵測因子在不同閾值下的表現(xiàn)，其中TVB 偵測因子中分別設(shè)置M=100.0和M=10.0，MDHE 偵測因子中分別 設(shè) 置a=0.5和a=0.05，JST 偵 測 因 子 中 分 別 設(shè) 置IT=0.02和IT=0.01。如圖19，相同偵測閾值下，分維混合策略比單元混合策略的分辨率高。需要注意的是，相比于其他兩種偵測因子，在TVB 偵測因子下，分維的混合策略對結(jié)果的影響明顯大于調(diào)整閾值的影響。而MDHE 偵測下，分維混合CPR-CNNW2 方法比單元混合分辨率更高，而且閾值越大，分辨率越高。JST 偵測下，分維偵測的影響較小，閾值的影響較大，閾值越大分辨率越高。

圖19 截線y=0.40 上三種偵測因子在不同閾值下的密度分布Fig.19 Density distributions along the sliced line at y =0.40 for three indicators under different thresholds

3.5.3 CPR-CNNW2 格式與其他格式的對比

基于MDHE 偵測因子的分維混合CPR-CNNW2格式，在 400×200（DoFs=2 000×1 000）的矩形網(wǎng)格下的計算結(jié)果數(shù)值紋影（密度梯度云圖）如圖20 所示。在相同自由度下，與采用改進的五階WENO-Z+格式[35]得到的計算結(jié)果相比，混合CPR-CNNW2 格式的計算結(jié)果流場細節(jié)更豐富。從圖21 的截線處的密度分布結(jié)果來看，分維混合CPR-CNNW2 格式在旋渦處的分辨率更高。

圖20 全局/局部數(shù)值紋影對比Fig.20 Global/local numerical schlieren comparison

圖21 截線x=1.05 上CPR-CNNW2 和WENO 的密度分布對比Fig.21 Density distributions along the sliced line at x=1.05 compared between CPR-CNNW2 and WENO

4 結(jié)論

本文基于非結(jié)構(gòu)四邊形網(wǎng)格CPR 方法的子單元限制方法，對比了不同問題單元偵測因子對計算結(jié)果的影響，得到以下結(jié)論：

1）一維算例測試結(jié)果顯示，不同偵測因子對結(jié)果的分辨率產(chǎn)生影響。TVB、MDHE、JST 三種偵測因子中，TVB 對結(jié)果的抹平最大，MDHE 和JST 的分辨率較高。

2）二維算例結(jié)果顯示，基于多項式模態(tài)衰減的MDHE 偵測因子在偵測到很少一部分問題單元的情況下，計算穩(wěn)定，而且偵測本身花費的時間較小?；趬毫μ荻鹊腏ST 偵測因子偵測到的問題單元也較少，但是偵測本身花費時間較多。而TVB 偵測因子對流場的適應(yīng)性不強，對高馬赫數(shù)問題需要采取一些措施提高魯棒性，比如調(diào)整自由參數(shù)M或者添加過渡單元。

3）基于四邊形單元分維計算的特點，測試了分維混合CPR-CNNW2 策略。和單元混合相比，分維混合可以提高分辨率。在含有明顯方向特征的流場中，劃分均勻網(wǎng)格并考慮分維偵測可以減少問題單元的標記，減小耗散的引入。激波旋渦干擾的算例顯示，采取TVB 和MDHE 偵測因子時，分維混合策略計算穩(wěn)定且結(jié)果分辨率更高；采取JST 偵測因子時，分辨率相差不大。同時，閾值的調(diào)整也會對計算結(jié)果的分辨率產(chǎn)生影響。

綜合來看，MDHE 偵測因子計算過程花費時間少，對問題單元偵測比較準確，且在分維混合措施下對計算結(jié)果的分辨率有所提高，是一種適合于四邊形網(wǎng)格的問題單元偵測因子。這三種偵測因子都需要人工設(shè)置參數(shù)，若不能有效準確地確定閾值，會造成計算的浪費。下一步我們將研究基于流場自適應(yīng)選擇閾值，進而減少人工干預，進一步提高計算效率。另外，通過將三角形單元分解為三個四邊形單元的方式，該子單元限制方法可以推廣到三角形網(wǎng)格或者混合網(wǎng)格，但這對網(wǎng)格質(zhì)量的要求可能更高。