一種基于幅值和波數(shù)的耗散控制方法

魏皇生，黃柱，席光

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

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

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

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

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

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

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

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

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

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

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

1 數(shù)值方法

1.1 理論基礎

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

式中：Q為守恒變量；F為守恒通量；t為時間；x為空間坐標。

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

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

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

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

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

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

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

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

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

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

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

對其進行傅里葉分析可以得到：

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

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

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

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

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

1.3 小尺度波動

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.4 波數(shù)控制

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

對κi進行求導可得

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

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

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

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

1.5 幅值控制

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

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

1.6 尺度可控混合格式

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

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

在WENO-Z 的基礎上，TENO 進行尺度分離：

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

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

對γk進行無量綱化可以得到

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

之后對χk進行截斷從而區(qū)分光滑和不光滑的模板：

式中：δk為截斷系數(shù)；CT為截斷閾值，本文采取推薦值10-7。

最后，可以得到TENO 格式的非線性權為

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

當流動較為光滑時，6 階TENO 格式的非線性權ωk與理想權dk相等，此時TENO 格式將會恢復到6 階線性中心格式：

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

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

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

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

式中：子模版的選取與6 點WENO 格式相同。

此外，文獻［14］推薦式（39）中ε=10-4。然而，文獻［32］指出當ε較大時，會起到對光滑指示器進行截斷的效果，這導致對于不同問題往往需要調(diào)整ε以獲得最佳結果。為了規(guī)避ε的影響，在本文中，選取ε=10-40，并基于小尺度的幅值提出一種新的方法來替代ε的截斷效果。

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

因此，當式（42）的邏輯判斷為真時，Case 1 判斷為光滑。

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

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

2 算例驗證

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

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

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

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

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

2.1 一維相對坐標系下的正激波

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

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

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

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

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

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

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

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

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

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

如圖4 所示，在密度波區(qū)域，中心-TENO 格式和尺度可控的混合格式均能更好地保持密度波的振幅。尺度可控的混合格式與中心-TENO格式的計算結果幾乎相同，說明沒有過多的引入數(shù)值黏性，可以有效改善對大尺度流動結構的分辨能力。而在間斷區(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］，計算域設置為[0，1]×[0，1]，使用400×400 個網(wǎng)格進行離散。邊界條件均為狄利克雷邊界條件，初始條件為

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

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

圖5 二維黎曼問題t=0.8 時的密度等值線（在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 時的渦量等值線（渦量=0）Fig. 6 Vorticity contours of 2-D Riemann problem at t=0.8（vorticity is equal to 0）

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

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

2.4 雙馬赫反射問題

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

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

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

圖8 雙馬赫反射問題t=0.2 時的密度等值線（在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時的渦量等值線（渦量=0）Fig. 9 Vorticity contours of double Mach reflection problem at t=0.2 （ vorticity is equal to 0）

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

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

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

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

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

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

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

2.6 非物理波動尺度

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

圖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 計算效率

2.1 節(jié)～2.5 節(jié)的計算結果已經(jīng)展示了尺度可控混合格式對間斷的識別能力。本節(jié)對比了3 種格式的計算效率。表2 展示了2 種混合格式在計算過程中使用的TENO 格式的比例。可以看到，對于激波主導的算例，尺度可控混合格式的TENO 使用比例在3.2%左右，而中心-TENO 格式的使用比例在6.5%左右，尺度可控混合格式略好。而在三維算例中，尺度可控混合格式的TENO 使用比例僅為0.1%，而中心-TENO 格式的使用比例高達31.7%。

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

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

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

表3 不同算例的計算時間Table 3 Calculation time of different cases

3 結 論

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

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

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

3）基于小尺度的幅值提出一種新的方法替代ε在開關函數(shù)中的截斷效果，改善了對光滑區(qū)域的識別能力。在整個計算過程中，尺度可控混合格式平均僅有3.2% 左右的情況下啟用TENO 格式，而中心格式有6.5%左右的情況啟用TENO 格式。相較于TENO 格式，尺度可控混合格式效率提升2～4 倍。

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