魏皇生，黃柱，席光

西安交通大學(xué) 能源與動力工程學(xué)院，西安 710049

大渦模擬（Large Eddy Simulation， LES）可以用于研究如湍流的復(fù)雜流動問題，近年來受到越來越多的關(guān)注［1-2］。由于只解析大尺度結(jié)構(gòu)，LES 的最高分辨率集中在對湍流演化有重要影響的湍流慣性子區(qū)域。然而，可壓縮流中不僅存在如湍流的復(fù)雜流動結(jié)構(gòu)，還存在大量間斷。為了模擬間斷，數(shù)值格式往往具有較大的耗散和很強(qiáng)的非線性。但是，過大的耗散會抹平流動細(xì)節(jié)，且強(qiáng)非線性會誘發(fā)如“數(shù)值湍流”的小尺度非物理波動［3］。高質(zhì)量的LES 既要分辨大尺度流動結(jié)構(gòu)，又要抑制小尺度非物理波動，這對耗散控制方法提出了挑戰(zhàn)。

在面向超聲速流動的高精度方法中，加權(quán)本質(zhì)無振蕩格式（Weighted Essentially Non-Oscillatory Schemes， WENO）［4-5］可以根據(jù)局部流場的光滑程度動態(tài)地調(diào)整耗散，是具有代表性的激波捕捉格式［6］。在此基礎(chǔ)上，為了更好地分辨大尺度流動結(jié)構(gòu)和抑制小尺度非物理波動，人們對耗散的控制方法進(jìn)行了大量研究。

為了更好地分辨大尺度流動結(jié)構(gòu)，基于保單調(diào)技術(shù)，Balsara 和Shu［7］構(gòu)造了具有更小色散誤差與更高精度的WENO 格 式。Henrick［8］和Borges［9］等通過改進(jìn)光滑指示器，改善了WENO格式在奇點(diǎn)處精度下降的問題。為了降低耗散，Martín 等［10］通過引入順風(fēng)模板構(gòu)造了具有中心型背景模板的WENO-SYMOO 格式。在此基礎(chǔ)上，Hu 等［11］通過修改光滑因子，構(gòu)造了自適應(yīng)中心/迎風(fēng)格式的WENO-CU 格式，顯著降低了耗散。最近，F(xiàn)u 等［12-13］基于新的模板選擇和加權(quán)策略，提出了定向本質(zhì)無振蕩格式（Targeted Essentially Non-Oscillatory schemes， TENO），進(jìn)一步降低了格式的耗散。Hill 和Pullin［14］通過在光滑區(qū)域使用針對LES 進(jìn)行優(yōu)化的中心差分格式（TCD），在間斷附近使用WENO 格式，構(gòu)造了具有更高分辨率的混合格式。

針對小尺度非物理波動，Pirozzoli［15］指出，數(shù)值格式中有必要包含一定的耗散以抑制高波數(shù)的數(shù)值振蕩。為了對數(shù)值格式的計(jì)算誤差進(jìn)行分析，Pirozzoli［16］提出了近似色散關(guān)系方法（Approximate Dispersion Relation， ADR），得到了激波捕捉格式的耗散和波數(shù)的關(guān)系?；诤纳⒄`差大小應(yīng)與色散誤差相匹配的假設(shè)，Hu 等［17］提出了一種耗散與色散的關(guān)系，可以用于衡量計(jì)算誤差的傳播距離。Sun 等［18］根據(jù)Hu 等的耗散色散關(guān)系確定了格式中的耗散系數(shù)。最近，李妍慧等［19］提出了一種波數(shù)識別器，并基于此得到了一種根據(jù)局部流場的波數(shù)控制耗散的方法。

上述數(shù)值方法基于局部流場的光滑程度或者波數(shù)控制耗散，可以在提升大尺度流動結(jié)構(gòu)分辨能力的同時(shí)有效地抑制小尺度非物理波動。然而，劉君等［20］發(fā)現(xiàn)5 階WENO 格式在計(jì)算激波或接觸間斷時(shí)存在放大計(jì)算結(jié)果誤差的風(fēng)險(xiǎn)，揭示了上述耗散控制方法的局限性。因此，不能簡單地將如間斷區(qū)域引入的非物理波動視為一種高波數(shù)的數(shù)值振蕩。如何在穩(wěn)定計(jì)算激波的同時(shí)降低耗散仍然是一個亟需研究的問題［21］，這要求更加合理的耗散控制方法。

小尺度非物理波動不僅具有波數(shù)高的特點(diǎn)，同時(shí)具有幅值小的特點(diǎn)。為了更好地識別小尺度非物理波動，進(jìn)而更精確地控制耗散，本文提出了一種采用局部流場的幅值和波數(shù)控制耗散的方法。

首先，針對如激波主導(dǎo)的或各向同性湍流的具有強(qiáng)烈非定常性的問題，根據(jù)一維非定常歐拉方程，確定小尺度下不同物理量之間的關(guān)系，并根據(jù)數(shù)值實(shí)驗(yàn)或Kolmogorov 尺度理論［22］確定小尺度波動幅值的閾值。根據(jù)所采用的人工黏性格式，推導(dǎo)了小尺度所對應(yīng)的人工黏性通量幅值的閾值。通過比較人工黏性通量幅值的閾值與實(shí)際計(jì)算值，給出一種根據(jù)局部流場的幅值控制耗散的方法。

其次，基于傅里葉分析，構(gòu)造了6 階人工黏性的對稱黏性。這個人工黏性僅有2 階，在低波數(shù)區(qū)域耗散誤差增長迅速而在高波數(shù)區(qū)域耗散誤差增長緩慢。由于耗散誤差與人工黏性通量的幅值相關(guān)，因此6 階人工黏性通量和其對稱黏性通量幅值的比值是一個4 階小量。使用這個比值去控制6 階人工黏性，就可以得到一個非線性的10 階人工黏性。該人工黏性在中低波數(shù)區(qū)域耗散小，在高波數(shù)區(qū)域耗散迅速增大。

最后，同時(shí)采用上述2 種耗散控制方法，就可以針對小幅值或高波數(shù)的區(qū)域，采用較大的耗散來抑制小尺度的非物理波動。為了獲得激波捕捉能力，將上述人工黏性格式與中心格式和TENO 格式進(jìn)行混合，從而構(gòu)造尺度可控的混合格式。

本文的具體內(nèi)容如下，1.1 節(jié)介紹有限差分格式的理論基礎(chǔ)；1.2 節(jié)討論了非物理波動的幅值和波數(shù)；1.3 節(jié)推導(dǎo)并討論了小尺度波動；1.4 節(jié)和1.5 節(jié)分別介紹波數(shù)控制與幅值控制；1.6 節(jié)介紹尺度可控的混合格式的實(shí)施方法。第2 節(jié)進(jìn)行數(shù)值測試，驗(yàn)證所提出的格式具有很好的分辨率并可以有效抑制小尺度非物理波動。第3 節(jié)進(jìn)行總結(jié)。

1 數(shù)值方法

1.1 理論基礎(chǔ)

以一維雙曲型方程組為例：

式中：Q為守恒變量；F為守恒通量；t為時(shí)間；x為空間坐標(biāo)。

假定對式（1）在空間方向上使用均勻一致的網(wǎng)格Δx=xi-1/2-xi+1/2進(jìn)行離散，則有

式中：L為守恒形式的空間導(dǎo)數(shù)算子；F為對數(shù)值通量h的高階逼近。h可以隱式地定義為

對于歐拉方程，為了改善高階格式的穩(wěn)定性，通常使用左特征矩陣R-1將守恒變量和守恒通量投影到特征空間進(jìn)行求解：

式中：q和f分別為特征變量與特征通量。為簡單起見，在下文中用標(biāo)量形式進(jìn)行分析，記q和f分別為標(biāo)量形式的變量與通量。

Kim 和Kwon［23］指出，迎風(fēng)型 的WENO 格 式可以拆分成一個中心格式和一個非線性耗散項(xiàng)；之后，Sun 等［18］通過釋放數(shù)值格式中的自由度，從而可以獨(dú)立地優(yōu)化色散和控制耗散；最近，F(xiàn)u等［24］也提出了一套可以獨(dú)立優(yōu)化色散和耗散的框架，這些研究為獨(dú)立地控制耗散提供了理論基礎(chǔ)。在下文中，將數(shù)值通量拆分成中心通量和人工黏性通量：

1.2 非物理波動的幅值和波數(shù)

針對光滑區(qū)域中產(chǎn)生的非物理波動，Hu等［17］使用一種耗散色散關(guān)系去估計(jì)必要的數(shù)值耗散大小，其可以表示為

式中：ζ定性表示由于色散誤差導(dǎo)致的數(shù)值偽波包可傳播的網(wǎng)格個數(shù)；κ為波數(shù)；κr為修正波數(shù)的實(shí)部；κi為修正波數(shù)的虛部；ε=10-3為一個可容許的誤差界限。式（6）中分子部分表示數(shù)值偽波包的群速度誤差尺度；分母部分表示耗散數(shù)值偽波包所需要的時(shí)間尺度。當(dāng)|dκr/dκ-1|～ε時(shí)，ζ～O(1) ，即數(shù)值偽波包的傳播距離僅為個位數(shù)網(wǎng)格間距。注意到，色散誤差主要集中在高波數(shù)區(qū)域，因此數(shù)值格式可以在中低波數(shù)區(qū)域保持較低的耗散，而在高波數(shù)區(qū)域引入較多的耗散，這就是基于波數(shù)控制耗散的思想。

然而，上述分析不適用于間斷區(qū)域產(chǎn)生的非物理波動。為了說明這個問題，考慮一個僅存在一道間斷的一維問題。當(dāng)網(wǎng)格較密時(shí)，光滑區(qū)域的數(shù)值通量逼近精確值，即

式中：h(x)和(x)分別為標(biāo)量形式的數(shù)值通量及其逼近；xs為間斷所在位置。Casper 和Carpenter［25-26］發(fā)現(xiàn)，除非間斷正好處于網(wǎng)格邊界上，激波捕捉方法僅有一階精度，與數(shù)值方法的設(shè)計(jì)精度無關(guān)。因此，在間斷附近的數(shù)值通量的誤差會遠(yuǎn)大于光滑區(qū)域，即

為簡單起見，不妨假定數(shù)值通量的誤差函數(shù)為一脈沖函數(shù)，其滿足式（9）和式（10）。

對其進(jìn)行傅里葉分析可以得到：

式中：F（·）表示傅里葉變換。式（11）的能譜為

其中：F*（·）為F（·）的共軛。因此，間斷區(qū)域產(chǎn)生的非物理波動服從均勻一致的能譜，這些非物理波動會傳播到并污染光滑區(qū)域的流動結(jié)構(gòu)。這使得一系列中心型的激波捕捉格式或混合格式不再合理，因?yàn)樗鼈冊谥械筒〝?shù)區(qū)域會完全恢復(fù)為中心型格式，缺乏對這部分非物理波動的抑制能力，最終形成大量的“數(shù)值湍流”。而迎風(fēng)型激波捕捉格式盡管能更好地抑制間斷區(qū)域產(chǎn)生的非物理波動，但容易抹平流動細(xì)節(jié)。上述分析僅僅是一個簡單的包含激波的算例，對于更加復(fù)雜的問題，上述現(xiàn)象會更加嚴(yán)重，這說明了基于波數(shù)控制耗散的局限性。

另一方面，間斷附近的計(jì)算誤差主要由網(wǎng)格與間斷不匹配或通量分裂格式所引起，可以視為一個小量。因此，可以定性認(rèn)為數(shù)值計(jì)算中產(chǎn)生的非物理波動是小幅值或高波數(shù)的小尺度波動。

為了抑制小尺度非物理波動，本文采取幅值和波數(shù)分別進(jìn)行控制。首先，在小幅值區(qū)域添加人工黏性，從而抑制間斷區(qū)域產(chǎn)生的小幅值非物理波動；其次，在高波數(shù)區(qū)域增強(qiáng)人工黏性的耗散，從而抑制光滑區(qū)域產(chǎn)生的高波數(shù)非物理波動。

在此強(qiáng)調(diào)，歐拉網(wǎng)格下的激波捕捉方法必然會在間斷附近產(chǎn)生非物理波動，上述耗散控制的目的在于減少非物理波動對光滑區(qū)域的影響，無法徹底消除該非物理波動。若要減少間斷區(qū)域產(chǎn)生的非物理波動，需要加密網(wǎng)格或改用其他方法。

1.3 小尺度波動

如1.2 節(jié)所述，非物理波動往往是具有小幅值或高波數(shù)的小尺度波動。因此，當(dāng)某局部區(qū)域物理量的波動幅值小于某一閾值時(shí)，可以認(rèn)為該區(qū)域存在小尺度的非物理波動。為了得到小尺度非波動幅值的閾值，首先需要推導(dǎo)小尺度下物理量之間的關(guān)系?？紤]一維非定常歐拉方程，將密度方程代入到動量方程中，可以得到：

式中：ρ為密度；u為速度；p為壓力。在可壓縮流動中，往往會出現(xiàn)如激波和湍流等復(fù)雜情況。針對這些具有強(qiáng)烈非定常性的問題，可以忽略式（13）中的定常部分ρu(?u/?x)的影響，得到小尺度下密度與速度波動之間的關(guān)系式：

式中：ρΔ、uΔ和pΔ分別為小尺度下的密度波動、速度波動和壓力波動；τ和η為小尺度下的時(shí)間和空間；Ma為馬赫數(shù)。

根據(jù)壓力和密度之間的關(guān)系，可以得到：

式中：γ=1.4 為空氣比熱比；c為聲速。若給定小尺度波動幅值的閾值為σ=uΔ/u，則只需根據(jù)所計(jì)算問題的給定其他特征量ρ、p、u、Ma，就可以估計(jì)出小尺度下的物理量波動ρΔ、pΔ、uΔ。易得，小尺度下的守恒變量為

其中：e為總能量。類似的，可以得到小尺度下的守恒通量FΔ以及特征空間的小尺度變量qΔ和通量fΔ：

根據(jù)這些小尺度下的物理量，可以識別出小幅值的非物理波動，從而更合理地控制格式的耗散并區(qū)分光滑區(qū)域和間斷區(qū)域。具體實(shí)施過程見1.5 節(jié)和1.6 節(jié)。

本文重點(diǎn)考慮激波主導(dǎo)的問題和具有經(jīng)典理論的各向同性湍流問題。對于不同問題，小尺度波動幅值的閾值可能會差別很大，因此需要采用不同的方法進(jìn)行確定。

針對激波主導(dǎo)的問題，本文在2.1 節(jié)中通過正激波算例研究了閾值σ的敏感性。當(dāng)閾值σ較大時(shí)，會導(dǎo)致對間斷的誤判；當(dāng)閾值σ較小時(shí)，無法有效抑制小幅值的非物理波動。 考慮σ=0.1，其含義為小尺度的波動要比宏觀量小一個量級。注意到，當(dāng)σ在0.1 附近時(shí)，計(jì)算結(jié)果對σ并不敏感，因此對于所有激波主導(dǎo)的問題，均選取σ=0.1。

對于具有成熟理論的各向同性湍流問題，可以類比Kolmogorov 尺度與積分尺度的關(guān)系，根據(jù)所求解問題的特征參數(shù)得到網(wǎng)格尺度波動幅值的閾值σΔx，使用該閾值可以得到更精確的結(jié)果。

根據(jù)能量的耗散率和輸運(yùn)率應(yīng)當(dāng)相等的條件，基于量綱分析，Kolmogorov 尺度與積分尺度的物理量存在如下關(guān)系：

式中：ηk為Kolmogorov 尺度；? 為積分尺度；Re為雷諾 數(shù)；uk為Kolmogorov 尺度的速度波動；ν為動力黏性系數(shù)；u為積分尺度的速度波動。為了將網(wǎng)格無法解析的小尺度非物理波動耗散掉，在網(wǎng)格尺度，數(shù)值的耗散率與輸運(yùn)率也應(yīng)當(dāng)相等。因此，對于網(wǎng)格尺度遠(yuǎn)大于Kolmogorov 尺度的情況，通過用Δx替代式（20）中的ηk，結(jié)合式（21），可以類比得到網(wǎng)格尺度下的速度波動uΔx與積分尺度的速度波動u之間的關(guān)系為

值得注意的是，式（22）并不代表數(shù)值計(jì)算中某波數(shù)區(qū)域內(nèi)的速度遵從1/3 冪律，僅用于估計(jì)網(wǎng)格尺度下的速度波動幅值大小。此外，與波動能量的-5/3 冪律不同，上述推導(dǎo)沒有用到慣性子區(qū)中流體黏性與耗散率相比是次要的這一假設(shè)［27-28］。唯一的假設(shè)為在網(wǎng)格尺度Δx附近，數(shù)值的耗散率和輸運(yùn)率應(yīng)當(dāng)匹配。

對于曲線坐標(biāo)系下的問題，可以通過坐標(biāo)變換得到式（16）～式（19）在計(jì)算坐標(biāo)系下的對應(yīng)物理量。然而，仍需進(jìn)一步研究上述耗散控制方法對曲線坐標(biāo)系下產(chǎn)生的坐標(biāo)變換誘導(dǎo)誤差［29］的影響。此外，對于高各向異性的湍流、強(qiáng)分離或轉(zhuǎn)捩等復(fù)雜情形，仍需要根據(jù)數(shù)值實(shí)驗(yàn)或經(jīng)典理論確定閾值σ，這將是接下來研究的重點(diǎn)。

1.4 波數(shù)控制

如1.2 節(jié)所述，非物理波動具有高波數(shù)的特點(diǎn)。因此數(shù)值計(jì)算中的高波數(shù)流動細(xì)節(jié)往往是非物理的，需要加以控制，本節(jié)采用波數(shù)對人工黏性進(jìn)行控制。在Jameson 等的人工黏性方法［30］中，人工黏性的大小依賴于流動的光滑程度。然而這種方法容易對極值點(diǎn)產(chǎn)生誤判。Sun［18］和Fu［24］等主要通過采用線性的人工黏性合理地控制黏性。對于6 點(diǎn)模板的通量格式，線性的人工黏性最高僅有6 階，因此無法在對高波數(shù)區(qū)域產(chǎn)生足夠耗散的同時(shí)較好地分辨中低波數(shù)流動細(xì)節(jié)。而基于傅里葉分析，可以構(gòu)造出更高階的人工黏性，從而在對高波數(shù)區(qū)域保持足夠耗散的同時(shí)減少對低波數(shù)區(qū)域的耗散。首先，6 個模板點(diǎn)的線性人工黏性通量可以表示為

對κi進(jìn)行求導(dǎo)可得

注意到式（25）中等號右邊第2 項(xiàng)在κ∈(0，π]區(qū)間的積分為0，而第1 項(xiàng)和第3 項(xiàng)共同的系數(shù)共同決定波數(shù)κ=π 時(shí)的耗散誤差。特別地，第2 項(xiàng)的系數(shù)可以用于調(diào)整耗散誤差隨著波數(shù)的增長速率在低波數(shù)和高波數(shù)的比例關(guān)系。對于5 階迎風(fēng)格式所對應(yīng)的6 階人工黏性，其系數(shù)為c1=1/60，c2-c1=-6/60，c3-c2=15/60 ，其耗散誤差在低波數(shù)增長緩慢，而在高波數(shù)增長迅速。將式（23）中第2 項(xiàng)取反數(shù)，可以得到c1=1/60，c2-c1=6/60，c3-c2=15/60，其耗散誤差在低波數(shù)增長迅速，在高波數(shù)增長緩慢。這種人工黏性作為6 階人工黏性的對稱黏性項(xiàng)，是一個耗散非常大的2 階人工黏性。由于耗散誤差與人工黏性通量的幅值相關(guān)，因此6 階人工黏性通量和其對稱黏性通量幅值的比值是一個4 階小量。使用這個比值去控制6 階人工黏性的通量，就可以得到10 階精度的非線性人工黏性：

式 中：為10階非線性人工黏性；為6 階線性人工黏性為6 階線性人工黏性的對稱黏性；ε=10-40。

為了展示非線性人工黏性的特點(diǎn)，對各種格式進(jìn)行傅里葉分析，其結(jié)果如圖1 所示。與預(yù)期的一樣，6 階人工黏性的對稱黏性和5 階迎風(fēng)格式的耗散誤差具有對稱性，標(biāo)定過的1 階迎風(fēng)格式位于兩者中間。10 階非線性人工黏性的耗散誤差在波數(shù)2.1 之前小于9 階迎風(fēng)格式，在高波數(shù)逐漸恢復(fù)到5 階迎風(fēng)格式。盡管傅里葉分析不能完全地展示非線性的影響，但上述分析仍有一定參考意義，后續(xù)的數(shù)值實(shí)驗(yàn)將展示10 階非線性人工黏性對中低波數(shù)流動細(xì)節(jié)的保持和高波數(shù)非物理波動的控制能力。

圖1 各種格式的耗散特性Fig. 1 Dissipation properties for various schemes

1.5 幅值控制

如1.2 節(jié)所述，非物理波動具有小幅值的特點(diǎn)。因此數(shù)值計(jì)算中幅值小于閾值σ的波動是非物理的，需要加以控制。注意到，6 階人工黏性的耗散誤差在波數(shù)＜2π/5 時(shí)耗散誤差比較小。因此假定一道波數(shù)為2π/5，幅值為A的小尺度波動正弦波。定義在不同相位下，6 階人工黏性對于該正弦波計(jì)算出的最大值為過濾系數(shù)εf，易得

在保證光滑區(qū)域數(shù)值格式連續(xù)［31］的前提下，對10 階非線性人工黏性進(jìn)行幅值控制：

1.6 尺度可控混合格式

1.4 和1.5 節(jié)所提出的方法可以根據(jù)幅值和波數(shù)控制耗散，但是無法分辨出如間斷的不光滑流動結(jié)構(gòu)。為了獲得激波捕捉能力，需要將上述方法與激波捕捉格式進(jìn)行混合。在激波捕捉格式中，TENO 格式克服了WENO 格式在奇點(diǎn)處降階的問題，往往具有更好的分辨率。因此，本文選用6 階TENO 格式［12］進(jìn)行混合，構(gòu)造了尺度可控混合格式。TENO 格式繼承WENO 格式的思想，其通量可以通過如下方法得到：

表1 子模板數(shù)值通量的系數(shù)Table 1 Coefficients of numerical flux of sub-stencils

在WENO-Z 的基礎(chǔ)上，TENO 進(jìn)行尺度分離：

式中：γk為尺度分離的光滑因子；τ6為以6 點(diǎn)模板計(jì)算的參考光滑因子；βk為每個子模版的光滑因子；ε=10-40是一個足夠小的數(shù)，其作用僅為防止分母為0。根據(jù)Jiang 和Shu 的定義［5］，光滑因子βk的計(jì)算式為

參考光滑因子τ6可以定義為

對γk進(jìn)行無量綱化可以得到

式中：χk為無量綱化的光滑因子。

之后對χk進(jìn)行截?cái)鄰亩鴧^(qū)分光滑和不光滑的模板：

式中：δk為截?cái)嘞禂?shù)；CT為截?cái)嚅撝担疚牟扇⊥扑]值10-7。

最后，可以得到TENO 格式的非線性權(quán)為

式中：dk為理想權(quán)，具體為d0=9/20，d1=6/20，d2=1/20，d3=4/20。

當(dāng)流動較為光滑時(shí)，6 階TENO 格式的非線性權(quán)ωk與理想權(quán)dk相等，此時(shí)TENO 格式將會恢復(fù)到6 階線性中心格式：

為了在更合理地控制耗散的同時(shí)獲得激波捕捉能力，將基于幅值和波數(shù)的耗散控制方法與TENO 格式進(jìn)行混合，構(gòu)造了尺度可控混合格式。其數(shù)值通量為

式中：為混合格式的通量；下標(biāo)“comp”和“char”指該通量是通過組份方向或特征方向重構(gòu)得到的；Case 1 指在坐標(biāo)方向判斷為光滑的情況；Case 2 指在坐標(biāo)方向判斷為光滑、特征方向判斷為不光滑的情況；Case 3 指在坐標(biāo)方向和特征方向均判斷為不光滑的情況。在這套混合格式結(jié)構(gòu)中，坐標(biāo)方向的混合可以免除特征變換和非線性格式的計(jì)算，從而提升計(jì)算效率；特征方向的混合可以最大程度上減少非線性機(jī)制的啟動。特別地，坐標(biāo)方向上僅使用密度進(jìn)行判斷，這進(jìn)一步地減少了計(jì)算量。

是否光滑可以通過光滑指示器SI 和一個人為給定的閾值TH 判斷。當(dāng)SI＜TH 時(shí)判斷為不光滑，否則為光滑。本文所采用的光滑指示器是由Hill 和Pullin［14］提出的WENO形式的開關(guān)函數(shù)進(jìn)行推廣得到的，定義為

式（39）含義為計(jì)算任意相鄰的3 個子模版的最小值和最大值的比值，這樣做的目的在于防止模板增大時(shí)對光滑區(qū)域的誤判。特別地，為了改善格式的數(shù)值對稱性，將Jiang 和Shu 對βk的定義中所使用的積分區(qū)間改為[xi，xi+1]：

式中：子模版的選取與6 點(diǎn)WENO 格式相同。

此外，文獻(xiàn)［14］推薦式（39）中ε=10-4。然而，文獻(xiàn)［32］指出當(dāng)ε較大時(shí)，會起到對光滑指示器進(jìn)行截?cái)嗟男Ч?，這導(dǎo)致對于不同問題往往需要調(diào)整ε以獲得最佳結(jié)果。為了規(guī)避ε的影響，在本文中，選取ε=10-40，并基于小尺度的幅值提出一種新的方法來替代ε的截?cái)嘈Ч?/p>

該方法的核心思想為將所有波動幅值小于閾值的小尺度非物理波動均視為光滑區(qū)域，并通過幅值控制耗散掉。因此假定一道波數(shù)為π、幅值為A的小尺度波動正弦波。定義在不同相位下，光滑因子βk的最大值為過濾系數(shù)εf，易得

因此，當(dāng)式（42）的邏輯判斷為真時(shí)，Case 1 判斷為光滑。

最后，在坐標(biāo)方向，閾值THcomp=1/4 設(shè)為一個非常保守的數(shù)；在特征方向，使用特征通量進(jìn)行判斷，閾值THcomp=1/14 使用Zhao 等［33］的推薦值（Hill和Pullin 的推薦值為1/150～1/120［14］）。

為了更好地顯示本文所提出的耗散控制方法的優(yōu)越性，在此定義一種中心-TENO 混合格式。其與尺度可控混合格式在實(shí)施上的區(qū)別在于，不施加耗散和基于幅值的光滑判斷。此外，ε取為文獻(xiàn)［14］的推薦值ε=10-4，并隨問題調(diào)整。

2 算例驗(yàn)證

為了展示基于幅值和波數(shù)的耗散控制方法的優(yōu)越性，展示了4 個激波主導(dǎo)的標(biāo)準(zhǔn)算例和三維可壓縮各向同性衰減湍流的計(jì)算結(jié)果。

其中，激波主導(dǎo)的無黏算例存在奇異性。即黏性越小，雷諾數(shù)越高，流動細(xì)節(jié)越豐富；一旦黏性為0，卻“搓不出渦”了。由于激波捕捉過程需要一定的數(shù)值黏性，計(jì)算結(jié)果與數(shù)值方法的修正方程有關(guān)。因此數(shù)值黏性越小，修正方程的雷諾數(shù)越大，流動細(xì)節(jié)越多。這些算例的流動細(xì)節(jié)的豐富程度可以用于評估數(shù)值方法引入的數(shù)值黏性大小。

盡管修正方程不再無黏，但對于波后未受到影響的區(qū)域，渦量仍應(yīng)當(dāng)?shù)扔?。渦量的強(qiáng)度越小、旋渦的尺度越大，說明能更好地抑制小幅值或高波數(shù)的小尺度非物理波動。

針對激波主導(dǎo)的問題，特征物理量ρ、p、u、Ma選為波前波后參數(shù)的平均值。針對可壓縮各向同性衰減湍流算例，密度和壓力均取初始值，速度與積分尺度均取最大含能渦所在波數(shù)值。

所有格式均采用3 階龍格-庫塔法進(jìn)行顯式時(shí)間推進(jìn)。坐標(biāo)方向使用局部Lax-Friedrichs 矢流通量分裂格式。按照文獻(xiàn)［13］中推薦的，特征方向使用Rusanov 矢流通量分裂格式。黏性項(xiàng)使用6 階中心格式離散。

2.1 一維相對坐標(biāo)系下的正激波

參考文獻(xiàn)［20］中的正激波算例，計(jì)算域設(shè)置為[0，1]，使用200個網(wǎng)格進(jìn)行離散。激波運(yùn)動速度為3.0。本算例固結(jié)于激波的相對坐標(biāo)系下進(jìn)行計(jì)算，激波位置固定在x=0.5處。波前波后參數(shù)為

計(jì)算時(shí)長為t=100 。由于在相對坐標(biāo)系下，可以統(tǒng)計(jì)出所有相對位置物理量的最大值和最小值，從而可以更直觀地比較激波在從初始間斷演化成數(shù)值激波的過程中產(chǎn)生的非物理波動。

圖2 研究了激波主導(dǎo)的問題對閾值σ的敏感性。圖中，黑色、紅色、藍(lán)色線代表不同閾值σ的計(jì)算結(jié)果。實(shí)線為t=100 時(shí)刻的瞬時(shí)值；長虛線為整個計(jì)算過程中各點(diǎn)處物理量的最大值；點(diǎn)線為最小值。好的數(shù)值方法應(yīng)當(dāng)僅在初始間斷附近形成一道數(shù)值激波，在遠(yuǎn)離間斷的位置，物理量的最大值和最小值應(yīng)始終與初始值保持一致。因此，圖2 中的長虛線和點(diǎn)線可以代表非物理波動的傳播與衰減過程?？紤]σ=0.1，其含義為小尺度的波動要比宏觀量小一個量級。數(shù)值實(shí)驗(yàn)的結(jié)果發(fā)現(xiàn)，σ?0.1 時(shí)，對小尺度非物理波動的抑制能力下降；σ?0.1 時(shí)，會在間斷附近出現(xiàn)振蕩。進(jìn)一步實(shí)驗(yàn)發(fā)現(xiàn)，σ在0.1 附近時(shí)對數(shù)值實(shí)驗(yàn)的影響不明顯。因此，對于激波主導(dǎo)的算例，σ=0.1 可以良好地區(qū)分光滑區(qū)域和間斷區(qū)域并抑制小尺度的非物理波動。

圖2σ 對激波主導(dǎo)問題的敏感性Fig. 2 Sensitivity of σ to shock dominant problem

圖3展示了各種格式產(chǎn)生的非物理波動。3 種格式均能得到穩(wěn)定的激波結(jié)構(gòu)，尺度可控混合格式的小尺度非物理波動略小于TENO 格式，遠(yuǎn)小于中心-TENO 格式。在遠(yuǎn)離激波處，尺度可控混合格式的計(jì)算誤差依然存在衰減趨勢；TENO 格式和中心-TENO 的計(jì)算誤差在一定距離之后保持不變。這是因?yàn)門ENO 和中心-TENO 格式僅依據(jù)光滑程度控制耗散，當(dāng)非物理波動較為光滑時(shí)會退化為無耗散的中心格式，無法進(jìn)一步加以抑制。而根據(jù)幅值和波數(shù)控制耗散不會受到非物理波動光滑程度的影響。

圖3 計(jì)算過程中的最大值和最小值及t=100 時(shí)的數(shù)值激波Fig. 3 Maximum and minimum values in calculation process and numerical shock wave at t=100

最后，在激波附近區(qū)域，不同格式的非物理波動完全相同，這是因?yàn)楦鞣N格式均采用了相同的格式對間斷進(jìn)行捕捉。特別強(qiáng)調(diào)，本文無意于消除激波捕捉過程中產(chǎn)生的非物理波動，本文所提出的幅值和波數(shù)控制僅用于抑制已經(jīng)產(chǎn)生的小尺度非物理波動的進(jìn)一步傳播。

2.2 激波密度波相互作用問題

現(xiàn)在考慮Shu 和Osher 提出的激波密度波相互作用問題［34］，計(jì)算域設(shè)置為[0，10]，使用200個網(wǎng)格進(jìn)行離散，初始條件為

選取計(jì)算時(shí)長為t=1.8。一道馬赫數(shù)為3的激波與密度波相互作用，并產(chǎn)生豐富的流動細(xì)節(jié)和間斷。因此，該問題可用于檢驗(yàn)對流動細(xì)節(jié)的分辨能力。其中，“精確解”由經(jīng)典的5 階WENO 格式使用4 000 個網(wǎng)格點(diǎn)得到。

如圖4 所示，在密度波區(qū)域，中心-TENO 格式和尺度可控的混合格式均能更好地保持密度波的振幅。尺度可控的混合格式與中心-TENO格式的計(jì)算結(jié)果幾乎相同，說明沒有過多的引入數(shù)值黏性，可以有效改善對大尺度流動結(jié)構(gòu)的分辨能力。而在間斷區(qū)域，中心-TENO 格式在x=2.8 和x=3.6 均產(chǎn)生了數(shù)值振蕩。而尺度可控的混合格式則完全避免了數(shù)值偽波。

圖4 激波密度波相互作用問題的密度分布（t=1.8）Fig. 4 Density distribution of shock-density interaction problem （t=1.8）

2.3 二維黎曼問題

現(xiàn)在考慮Lax 和Liu 提出的二維黎曼問題［35］，計(jì)算域設(shè)置為[0，1]×[0，1]，使用400×400 個網(wǎng)格進(jìn)行離散。邊界條件均為狄利克雷邊界條件，初始條件為

選取計(jì)算時(shí)長為t=0.8。每個域之間的差異導(dǎo)致許多復(fù)雜的現(xiàn)象，包括激波、稀疏波、滑移線。這些由激波相互作用引起的渦流非常適合測試數(shù)值方案的分辨率。理想的數(shù)值格式應(yīng)該捕獲更多由滑移線引起的開爾文-亥姆霍茲不穩(wěn)定性以及它們在射流頭部的復(fù)雜相互作用引起的渦流。除此之外，在未受到激波相互作用問題干擾的波后區(qū)域，流場物理量應(yīng)當(dāng)均勻一致，也就是說這些區(qū)域的渦量應(yīng)當(dāng)?shù)扔?。

由圖5 展示t=0.8 時(shí)的密度等值線可知，TENO格式，中心-TENO格式，尺度可控的混合格式均能捕捉到大量流動結(jié)構(gòu)。相比TENO格式，中心-TENO 格式和尺度可控混合格式均能捕捉到更多的流動細(xì)節(jié)，且尺度可控混合格式捕捉到的流動細(xì)節(jié)更加豐富。結(jié)合圖6中不同格式的渦量等值線可知，尺度可控混合格式在波后區(qū)域渦量等值線的密度更小，更符合物理規(guī)律。相比起中心-TENO格式，尺度可控混合格式減少了小尺度非物理波動對主要流動結(jié)構(gòu)的干擾，增強(qiáng)了對光滑區(qū)域的識別能力。根據(jù)圖7中2個混合格式的間斷探測結(jié)果，可以發(fā)現(xiàn)混合格式在主要流動結(jié)構(gòu)附近更少地啟用TENO格式，因此更好地保留了流動細(xì)節(jié)。

圖5 二維黎曼問題t=0.8 時(shí)的密度等值線（在0.2～1.6 之 間33 等 分）Fig. 5 Density contours of 2-D Riemann problem at t=0.8（33 equally spaced from 0.2 to 1.6）

圖6 二維黎曼問題t=0.8 時(shí)的渦量等值線（渦量=0）Fig. 6 Vorticity contours of 2-D Riemann problem at t=0.8（vorticity is equal to 0）

圖7 二維黎曼問題t=0.8 時(shí)的間斷探測結(jié)果Fig. 7 Discontinuous detection result of 2-D Riemann problem at t=0.8

另一方面，2 個混合格式在數(shù)值不對稱性上也有一定的改進(jìn)。混合格式在光滑區(qū)域減少了強(qiáng)非線性機(jī)制的啟動，因此提升了數(shù)值對稱性。最后，尺度可控混合格式給出非常干凈的間斷識別結(jié)果，這歸功于小尺度波動的引入。

2.4 雙馬赫反射問題

現(xiàn)在考慮雙馬赫反射問題［36］，計(jì)算域設(shè)置為[0，4]×[0，1]，使用960×240 個網(wǎng)格進(jìn)行離散。在本例中，一道馬赫數(shù)為10 的激波以與x軸夾角為60°的方向向右移動，并與下壁面發(fā)生反射。入口和出口分別采用進(jìn)口與出口邊界條件。上邊界采用精確的波前波后參數(shù)，并隨著激波的移動而變化。在下邊界，x∈[0，1.666 7]被設(shè)定為波后邊界條件，其余部分被設(shè)定為反射邊界條件。初始條件為

選取計(jì)算時(shí)長為t=0.2。同樣的，在未受到激波相互作用問題干擾的波后區(qū)域，流場物理量應(yīng)當(dāng)均勻一致，即這些區(qū)域的渦量應(yīng)當(dāng)?shù)扔?。

由圖8 展示t=0.2 時(shí)的密度等值線可知，在滑移線附近，尺度可控混合格式能捕捉到更加細(xì)致的流動結(jié)構(gòu)。結(jié)合圖9 中不同格式的渦量等值線可知，尺度可控混合格式很好地控制了激波演化過程中產(chǎn)生的小尺度非物理波動。

圖8 雙馬赫反射問題t=0.2 時(shí)的密度等值線（在2～22之間33 等分）Fig. 8 Density contours of double Mach reflection problem at t=0.2（33 equally spaced from 2 to 22）

圖9 雙馬赫反射問題t=0.2時(shí)的渦量等值線（渦量=0）Fig. 9 Vorticity contours of double Mach reflection problem at t=0.2 （ vorticity is equal to 0）

根據(jù)圖10 中的間斷探測結(jié)果，尺度可控混合格式再一次給出非常干凈的間斷識別結(jié)果，說明尺度可控混合格式可以在不降低分辨率的前提下控制小尺度非物理波動并識別光滑區(qū)域。

圖10 雙馬赫反射問題t=0.2 時(shí)的間斷探測結(jié)果Fig. 10 Discontinuous detection result of double Mach reflection problem at t=0.2

2.5 可壓縮各向同性湍流衰減問題

考慮可壓縮各向同性衰減湍流算例［37］，計(jì)算域設(shè)置為[0，2π]×[0，2π]×[0，2π]，網(wǎng)格分辨率為64×64×64。初始泰勒馬赫數(shù)為Mat=0.6，初始雷諾數(shù)Reλ=100，與文獻(xiàn)［37］中相同。采用周期性邊界條件。設(shè)定初始湍能譜為

式中：urms為速度均方根；κ0為峰值波數(shù)，本文取4。選取計(jì)算時(shí)長為t/τ=4，τ=2/(κ0urms)。在本算例中，存在最大含能渦，因此將積分尺度設(shè)定為最大含能渦尺度?=2π/κ0。

圖11展示了計(jì)算過程中速度、渦量、溫度的無量綱2 范數(shù)的變化曲線與DNS（Direct Numerical Simulation）結(jié)果的對比?？梢钥吹剑谟?jì)算的前期，如t/τ＜1 時(shí)，尺度可控混合格式的耗散與其他格式的相當(dāng)。當(dāng)t/τ＞1 時(shí)，隨著湍流衰減的演化，流場中高波數(shù)流動增加，TENO 和中心-TENO 格式的耗散顯著增加，而尺度可控混合格式仍然保持了較低的耗散。尺度可控混合格式與DNS 的計(jì)算結(jié)果吻合良好，說明網(wǎng)格尺度的推導(dǎo)可靠，在數(shù)值計(jì)算中起到湍流模型的作用。

圖11 可壓縮各向同性衰減湍流的計(jì)算結(jié)果Fig. 11 Calculation results of compressible isotropic decaying turbulence

2.6 非物理波動尺度

為了進(jìn)一步驗(yàn)證本格式對非物理波動的控制效果，研究非物理波動的幅值和長度尺度。因?yàn)? 個二維算例的計(jì)算結(jié)果相似，選擇二維黎曼問題中t=0.8 時(shí)刻在區(qū)域[0.85，0，95]×[0.85，0.95]處的渦量進(jìn)行研究。這個區(qū)域作為波后區(qū)域，理想情況下的渦量應(yīng)當(dāng)為0。圖12給出渦量幅值的云圖?？梢钥吹?，在該區(qū)域，TENO 格式的計(jì)算誤差不到0.1，中心-TENO格式的計(jì)算誤差不到0.7，尺度可控混合格式計(jì)算誤差不到0.004。圖13展示渦量為0 的等值線和網(wǎng)格的放大圖?？梢钥吹?，TENO 格式和中心-TENO 格式均呈現(xiàn)隨機(jī)混沌狀的小尺度非物理波動，而尺度可控混合格式的非物理波動呈條帶結(jié)構(gòu)。其次，TENO 格式的非物理波動尺度大約在4個網(wǎng)格長度左右，而中心-TENO 格式在2個網(wǎng)格長度左右，尺度可控混合格式的波動尺度在8～10 個網(wǎng)格左右。事實(shí)上，這是因?yàn)橹行?TENO 格式和TENO 格式僅依據(jù)光滑程度控制耗散，因此會把小幅值的非物理波動視為光滑區(qū)域。此時(shí)，2個格式均會恢復(fù)為無耗散的中心格式，因此難以抑制具有小幅值和高波數(shù)特點(diǎn)的小尺度非物理波動。

圖12 波后區(qū)域渦量幅值云圖的放大圖Fig. 12 Enlarged view of vorticity magnitude contours in post-shock region

圖13 波后區(qū)域的渦量等值線放大圖和網(wǎng)格Fig. 13 Enlarged view of vorticity contours and mesh in post-shock region

圖12和圖13 表示，簡單的基于光滑程度控制耗散無法有效抑制小尺度非物理波動；而基于幅值和波數(shù)，可以更好地抑制小尺度非物理波動。

2.7 計(jì)算效率

2.1 節(jié)～2.5 節(jié)的計(jì)算結(jié)果已經(jīng)展示了尺度可控混合格式對間斷的識別能力。本節(jié)對比了3 種格式的計(jì)算效率。表2 展示了2 種混合格式在計(jì)算過程中使用的TENO 格式的比例?？梢钥吹剑瑢τ诩げㄖ鲗?dǎo)的算例，尺度可控混合格式的TENO 使用比例在3.2%左右，而中心-TENO 格式的使用比例在6.5%左右，尺度可控混合格式略好。而在三維算例中，尺度可控混合格式的TENO 使用比例僅為0.1%，而中心-TENO 格式的使用比例高達(dá)31.7%。

表2 TENO 格式的使用比例Table 2 Usage ratio of TENO scheme

除了基于Kolmogorov 尺度推導(dǎo)的網(wǎng)格尺度波動幅值的閾值更加合理之外，導(dǎo)致上述現(xiàn)象的還有以下原因。首先，可壓縮衰減湍流不是激波主導(dǎo)的，同時(shí)還是一個三維算例，因此激波占比更低。其次，基于幅值和波數(shù)可以有效抑制小尺度非物理波動，從而確保大多數(shù)區(qū)域被識別為光滑的。最后，隨著時(shí)間的推進(jìn)，流場中高波數(shù)的流動結(jié)構(gòu)增加，傳統(tǒng)的開關(guān)函數(shù)容易對高波數(shù)的流動結(jié)構(gòu)產(chǎn)生誤判，而使用小尺度幅值輔助的開關(guān)函數(shù)不受流動結(jié)構(gòu)波數(shù)的影響。

表3 統(tǒng)計(jì)了不同算例的計(jì)算時(shí)間。盡管由于編程手段的差異，計(jì)算時(shí)間不能完全表示算法的計(jì)算效率，但仍有一定借鑒意義。首先，尺度可控混合格式與中心-TENO 格式的計(jì)算效率接近。盡管尺度可控的混合格式更少地使用了TENO 格式，但是增加了程序復(fù)雜程度。此外，尺度可控混合格式計(jì)算效率為TENO 格式的2～4 倍。隨著計(jì)算問題維數(shù)的增加，程序復(fù)雜性的增強(qiáng)，混合格式所帶來的效率提升被程序中的其他模塊逐漸稀釋。通過對程序進(jìn)行優(yōu)化，尺度可控混合格式的計(jì)算效率提升幅度會更大。

表3 不同算例的計(jì)算時(shí)間Table 3 Calculation time of different cases

3 結(jié) 論

本文提出一種基于幅值和波數(shù)控制耗散的方法。為了獲得激波捕捉能力，與TENO 格式進(jìn)行混合構(gòu)造了尺度可控混合格式。主要結(jié)論如下：

1）對于激波主導(dǎo)的或各向同性湍流的具有強(qiáng)烈非定常性的問題，根據(jù)一維非定常歐拉方程，推導(dǎo)了小尺度下不同物理量之間的關(guān)系，并通過數(shù)值實(shí)驗(yàn)或Kolmogorov 尺度理論確定小尺度波動幅值的閾值。基于該閾值構(gòu)造了耗散大小和局部流場幅值的關(guān)系，從而有效抑制小幅值的非物理波動。一維正激波的結(jié)果表明，非物理波動的幅值顯著減小。激波密度波相互作用算例中，在間斷區(qū)域避免了數(shù)值振蕩。二維算例（以二維黎曼問題為例）在波后區(qū)域的渦量分布表明，本格式在波后區(qū)域的渦量計(jì)算誤差不到0.004，遠(yuǎn)遠(yuǎn)小于TENO 格式的0.1 和中心-TENO 格式的0.7。

2）基于傅里葉分析構(gòu)造了10 階非線性人工黏性，其耗散誤差在波數(shù)2.1 之前小于9 階迎風(fēng)格式，在高波數(shù)逐漸恢復(fù)到5 階迎風(fēng)格式。二維算例在波后區(qū)域的渦量分布表明，尺度可控混合格式減少了非物理的數(shù)值湍流，其非物理波動的尺度在8～10 個網(wǎng)格左右，遠(yuǎn)遠(yuǎn)優(yōu)于TENO 格式的4 個網(wǎng)格和中心-TENO 格式的2 個網(wǎng)格。

3）基于小尺度的幅值提出一種新的方法替代ε在開關(guān)函數(shù)中的截?cái)嘈Ч纳屏藢饣瑓^(qū)域的識別能力。在整個計(jì)算過程中，尺度可控混合格式平均僅有3.2% 左右的情況下啟用TENO 格式，而中心格式有6.5%左右的情況啟用TENO 格式。相較于TENO 格式，尺度可控混合格式效率提升2～4 倍。

4）對于本文中的算例，尺度可控混合格式可以捕捉到更多的大尺度流動細(xì)節(jié)。激波密度波算例中更好地保持了密度波的幅值?？蓧嚎s各項(xiàng)同性衰減湍流算例與DNS 結(jié)果對比良好。即使相對于中心-TENO 格式，尺度可控混合格式分辨率也有所提升。