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

(1.空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000；2.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院，西安 710000)

0 引言

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

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

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

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

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

1 控制方程及離散

1.1 控制方程

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

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

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

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

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

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

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

式中，

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

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

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

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

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

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

式中，PK是Legendre 多項(xiàng)式。

時(shí)間離散格式采取三階Runge-Kutta 格式，用RHS表示式（1）的右端項(xiàng)：

則三階Runge-Kutta 格式為：

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

2 子單元限制策略

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

2.1 問題單元偵測(cè)因子

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

2.1.1 TVB 偵測(cè)因子

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

TVB 限制器的思想是：

式中，

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

2.1.2 MDHE 偵測(cè)因子

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

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

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

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

定義閾值為：

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

2.1.3 JST 偵測(cè)因子

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

一維情況下：

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

二維情況下考慮兩個(gè)方向：

式中，k、l代表單元內(nèi)部求解點(diǎn)的局部坐標(biāo)索引，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）給出了每個(gè)子單元或CPR 求解點(diǎn)處的指示值，使用體積加權(quán)，得到單元的單個(gè)指標(biāo)值：

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

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

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

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

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

圖1 子單元?jiǎng)澐址绞紽ig.1 Layout of subcell division

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

圖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)四邊形單元，每個(gè)方向上界面附近的兩個(gè)求解點(diǎn)u1和u5對(duì) 應(yīng)的子單元。圖3 是求解u1對(duì)應(yīng)子單元的情況，插值過程只在求解界面處時(shí)涉及到單元之間的數(shù)據(jù)處理。將的計(jì)算由計(jì)算空間改為物理空間上進(jìn)行（反距離加權(quán)代入物理空間距離），在物理空間上計(jì)算得到。

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

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

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

2.3 CPR-CNNW2 混合策略

2.3.1 格式切換

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

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

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

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

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

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

2.4 子單元時(shí)間步長限制

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

式中，d是空間維度，N是多項(xiàng)式階數(shù)，h代表網(wǎng)格尺寸，|λmax|是最大特征波速度。對(duì)二維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í)間步長限制下，等距子單元的最大劃分個(gè)數(shù)為Ns=2N+1，其中N是多項(xiàng)式次數(shù)。當(dāng)Ns＞2N+1時(shí)，子單元時(shí)間步長限制會(huì)反過來約束DG 的最大時(shí)間步長。本文的非等距劃分子單元方式，標(biāo)準(zhǔn)單元內(nèi)子單元尺寸符合Gauss 積分權(quán)分布，和求解點(diǎn)個(gè)數(shù)相關(guān)，最小的為第一積分權(quán) ω1代表的子單元長度。CPR 格式使用五階精度，Guass 積分使用五點(diǎn)格式時(shí)，ω1＞1/(2N+1)。故在此可以統(tǒng)一使用CPR 的時(shí)間步長限制。

3 算例測(cè)試

3.1 Sod 激波管問題

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

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

計(jì)算結(jié)果如圖5 所示。為了方便對(duì)比三種激波偵測(cè)因子下的問題單元分布，將圖5 右側(cè)縱坐標(biāo)軸設(shè)為問題單元指示器 β，并將三種方法錯(cuò)位分布。在接觸間斷x=0.7 和x=0.85 附近，使用三種問題單元偵測(cè)因子都無明顯振蕩。對(duì)比之下，無論采取何種偵測(cè)方式，混合CPR-CNNW2 格式都比CNNW2 格式對(duì)間斷的抹平小。雖然三種偵測(cè)因子最終的偵測(cè)結(jié)果中，問題單元只存在于間斷附近，但可以看到采取TVB方法對(duì)接觸間斷的抹平較大，而MDHE 和JST 方法對(duì)間斷的抹平較小，數(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（對(duì)應(yīng)DoFs=500），左右邊界分別為流入流出邊界條件，計(jì)算終止時(shí)間為T=0.18，CFL=0.9 。在此算例中TVB 偵測(cè)因子的自由參數(shù)取M=3.0。

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

另外，我們將本文的CPR-CNNW2 結(jié)果與文獻(xiàn)[31]中的p-weighted 限制器在取p=5次多項(xiàng)式時(shí)采用Nc=100（DoFs=600）網(wǎng)格下的計(jì)算結(jié)果進(jìn)行對(duì)比?？梢钥吹皆诟哳l振蕩區(qū)，MDHE 和JST 偵測(cè)因子下的計(jì)算結(jié)果的分辨率比p-weighted 限制器的計(jì)算結(jié)果稍好。

3.3 等熵渦測(cè)試

二維等熵渦問題的初始條件是在一個(gè)均勻流上添加一個(gè)各向異性的等熵旋渦的擾動(dòng)。初始條件設(shè)置見文獻(xiàn)[32]，取自由度為200×200，時(shí)間步長dt=0.000 1，計(jì)算到T=1.0?；诘褥販u問題，對(duì)CNNW2與CPR 的計(jì)算時(shí)間進(jìn)行對(duì)比測(cè)試。

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

表1 單元混合與分維混合的計(jì)算時(shí)間Table 1 Calculation time for hybrid by cell and by dimension

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

表2 不同偵測(cè)因子對(duì)計(jì)算時(shí)間的影響Table 2 Influence of different indicators on the calculation time

3.4 雙馬赫反射

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

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

圖8 CPR-CNNW2 在不同偵測(cè)因子下求解雙馬赫反射問題（從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 雙馬赫反射問題的各種偵測(cè)因子下的問題單元分布圖Fig.9 Distribution of problematic cells under different indicators for the double Mach reflection problem

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

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

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

3.5 激波-旋渦干擾

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

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

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

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

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

圖11 不同偵測(cè)因子下CPR-CNNW2 求解激波-旋渦干擾問題的計(jì)算結(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

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

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

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

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

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

圖16 中展現(xiàn)了三種偵測(cè)因子下分維混合CPRCNNW2 方法在截線處的密度分布結(jié)果。在當(dāng)前閾值下，MDHE 分維偵測(cè)下的計(jì)算結(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 計(jì)算激波-旋渦干擾問題的問題單元分布Fig.17 Distribution of problematic cells of CPR-CNNW2 with detecting by dimension for the shock-vortex interaction problem

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

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

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

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

3.5.3 CPR-CNNW2 格式與其他格式的對(duì)比

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

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

圖21 截線x=1.05 上CPR-CNNW2 和WENO 的密度分布對(duì)比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 方法的子單元限制方法，對(duì)比了不同問題單元偵測(cè)因子對(duì)計(jì)算結(jié)果的影響，得到以下結(jié)論：

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

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

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

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