高精度激波捕捉格式的性能分析

屈 峰 閻 超 于 劍 陳嘉陽

(北京航空航天大學 國家計算流體實驗室，北京100191)

無論是在學術界還是工程界，湍流與激波都是普遍存在并且至關重要的流動現象.鑒于湍流的極度復雜性，短時間求出其解析解的可能性極小[1]，目前人們只有通過DNS(直接數值模擬)、LES(大渦模擬)[2]等數值方法來對其進行近似模擬.但由于數值格式本身具有數值耗散，其計算得到的數值解會耗散掉較小尺度的渦.也就是說，這種數值上的耗散無論在計算精度方面還是在計算效率方面都對湍流的模擬起著負面的作用.與此相對地，為了不因為激波這種強間斷的存在而使計算發散，數值格式又必須具有一定程度的數值耗散以使計算穩定、無震蕩[3].因此，湍流結構的模擬與激波的模擬對數值格式性能上的要求是相互矛盾的，即:如果在包含激波的流動中直接應用數值耗散較小的格式會造成所謂的Gibbs(非物理的)震蕩[4－5]，從而導致計算的失真或者發散;如果在湍流的模擬中，不加任何處理地使用數值耗散較大的低階迎風格式，會導致精度較差，從而無法合理地模擬出湍流脈動結構.本文通過對目前工程中應用較為廣泛的三階MUSCL插值、五階 WENO(加權本質無震蕩)重構[6－8](基于特征型以及基于非特征型)以及五階WENO格式改進版Map-WENO[9]等高階格式系統地研究對比，得到了一些有關這些常用高階格式性能的結論，從而為今后格式的構造與工程化應用提供思路.

1 控制方程與計算方法

1.1 NS方程有限體積的一般形式

考慮如下三維雙曲守恒律組[10]:

定義空間單元Iijk為

在單元Iijk上對式(1)進行積分，整理得

其中

1.2 計算方法

1.2.1 MUSCL 格式

以界面i+1/2為例，MUSCL插值可以表示為[11]

其中Δqi為單元i中流場變量的坡度.

其表達式為

其中k為任意常數，且k∈[－1，1].

1.2.2 基于變量插值的WENO格式

為了獲得更高階的計算精度(大于等于三階精度)，本文沿著MUSCL格式的思路，采用提高變量插值精度的方法，即將MUSCL格式中的限制器用WENO方法來代替.上述改進方法只增加很小一部分計算量，卻實現了高階格式的優異性能[12].一些著名的 CFD(Computational Fluid Dynamics)代碼(如 CFL3D，OVERFLOW)也都采用了類似的思想來提高數值方法的精度.

其中

q'0，q'1，q'2代表3個三階精度的線性插值模板，通過非線性權將它們組合成1個五階精度的非線性插值模板.非線性權ωk的表達式如下:

其中，Ck稱為最優權，用于保證格式在光滑區的插值精度;ε為一小量，用來防止分母為0;IS為光滑指示因子.

1.2.3 Map-WENO 格式

實際應用中，盡管上文所述的五階WENO格式在求解大部分問題中都表現出了其“五階”的精度，但在許多問題特別是含有強間斷時其精度要小于理論精度.文獻[9]證明在通量一階導數消失而其三階導數不同時隨之消失的極點處，上述的WENO格式會損失部分精度以致理論上只能達到三階精度.

針對上述問題，Henrick等[9]通過構造滿足一定條件的映射函數對原WENO格式非線性系數的選取進行優化.其具體表達形式如下:

2 算例及結果分析

2.1 Taylor-Green 渦

通過理論分析可知，在預先給定初始值[3]的情況下，該問題所描述的無黏Taylor-Green渦將隨著時間推移不斷破裂形成更小尺度的渦，從而形成了渦尺度沒有最小邊界的問題.因此，這個算例的計算可以在一定程度上反映格式數值耗散的大小.

計算區域 xi∈[0，2π]，網格間距為 Δxi=2π/64.采取周期性邊界條件.計算初始條件為

圖1給出了幾種格式在不同時刻計算所得動能與初始時刻動能之比(MUSCL_Three-order為三階 MUSCL 插值[7]，WENO_Original代表基于原始變量即非特征型的五階 WENO格式，WENO_Character代表基于特征變量的五階WENO格式，Map-WENO表示前文所述的Map-WENO格式).其中，K為對應t(單位:s)時刻的動能值，K0為初始時刻的動能值.理論上，該動能比值隨著時間的推移保持不變.

圖1 各格式在不同時刻得到的動能與初始時刻動能之比隨時間變化的分布曲線Fig.1 Comparison of schemes’turbulent kinetic energy decays

從圖1可以清晰看出，相較于其他格式，三階MUSCL格式在起始的一段時間內具有較大的數值耗散.但隨著時間的推移，其耗散性弱于五階WENO格式.其原因應該是隨著時間的推移，流場中形成各種尺度相差很大的渦結構，而這使得流場中出現很多流動參數間斷.在捕捉這些間斷時，五階WENO格式中的非線性加權系數會逐漸遠離最優權以增強其在計算過程中的魯棒性，最終導致其數值耗散大于三階MUSCL格式.基于特征變量的五階WENO格式是為了計算穩定性而在特征場中進行插值的五階WENO格式，因此理論上其耗散性大于基于原始變量的五階WENO格式.Map-WENO格式是對原WENO格式非線性系數的優化，因此理論上其耗散性應弱于原始WENO格式，而圖1中也很準確地表明了相關特性.

上述結論亦可從圖2(4種格式在不同時刻求解得到的擬渦能與初始時刻擬渦能的比值分布曲線)以及圖3(4種格式在t=5 s時求解得到的能譜分布曲線)中分析得到.關于擬渦能以及能譜的具體求解方法，可參見文獻[13].

圖2 各格式在不同時刻得到的擬渦能隨時間的分布曲線Fig.2 Comparison of schemes’turbulent enstrophy decays

圖3 各格式在t=5 s時求解得到的能譜分布曲線Fig.3 Comparison of schemes’energy spectrum distributions at t=5 s

2.2 Shu-Osher問題

Shu-Osher問題是一維的激波湍流相互干擾問題.這個算例主要是用來驗證各格式的激波捕捉能力.

計算區域為 x∈[－5，5]，網格間距 Δx=0.05.計算初始條件為

圖4給出了通過3種格式計算得到的t=1.8 s時流場密度分布的局部放大圖.通過與精確解的對比可以發現，3類WENO格式捕捉到的激波較之于三階MUSCL格式更加銳利.其中Map-WENO格式捕捉到的激波最銳利，與精確解也最為接近.而基于原始變量的五階WENO格式較之于基于特征變量的五階WENO格式計算得到的激波更銳利一些，但在捕捉某些密度波時其解發生震蕩.這表明基于特征變量形式的五階WENO格式在捕捉間斷時更加穩定.

圖4 t=1.8 s時流場密度分布的局部放大圖Fig.4 Distribution of density at t=1.8 s

2.3 可壓縮各向均勻同性湍流[13－15]

計算區域 xi∈[0，2π]為周期性的，網格間距為Δxi=2π/64，采取周期性邊界條件.初始流動參數為:Reλ=50;MaT=1.2.湍流能譜[4]E(k)=參數為:k0=4，合適的A以保證初始的 urms等于1.

圖5給出了幾種格式計算得到的動能與初始動能之比隨時間變化的分布曲線.從圖中可以看出，三階MUSCL格式的數值耗散要小于2類原始的五階WENO格式，而Map-WENO格式的耗散性最弱.此結論與算例1中所得到的結論略有不同.其原因應該是此問題中在初始時刻流場中就有局部激波，為了穩定地捕捉這些小激波，五階WENO格式加大了自身在計算過程中的耗散性.此外，基于一般變量的五階WENO格式的耗散性要略小于基于特征變量的五階WENO格式.此結論與算例1中所得到的結論一致.上述結論亦可從圖6(4種格式在不同時刻得到的擬渦能分布曲線)中得到.但從圖6中可以得到，2類原始五階WENO格式所能得到的最大擬渦能遠大于三階MUSCL格式.這表明三階MUSCL格式對小尺度流動的耗散性大于五階WENO格式.此結論亦可從圖7(能譜圖)中得出.同時從圖7中可以看到，三階MUSCL格式在慣性子區的數值耗散要明顯小于2類原始的五階WENO格式而略大于Map-WENO格式.這應該是因為流場在慣性子區這個尺度范圍中有較大的局部間斷，2類原始的五階WENO格式為了穩定地捕捉間斷增加了數值耗散以致其超過了三階MUSCL格式的數值黏性.

圖5 各格式在不同時刻得到的動能與初始時刻動能之比隨時間變化的分布曲線Fig.5 Comparison of schemes’turbulent kinetic energy decays

圖6 各格式在不同時刻得到的擬渦能分布曲線Fig.6 Comparison of schemes’turbulent enstrophy decays

圖8和圖9分別給出了通過幾種格式計算所得到的t=1 s時密度以及速度梯度在流場中的分布.從圖8和圖9中可以看到幾種格式均可以準確地捕捉到激波.整體來說，線性MUSCL插值的耗散性弱于2類原始的WENO格式.Map-WENO格式的耗散性最弱.同時，相比較于整體耗散較小的MUSCL格式，2類原始WENO格式捕捉到的激波更加銳利.

圖7 各格式在t=1 s時求解得到的能譜分布曲線Fig.7 Comparison of schemes’energy spectrum distributions at t=1 s

圖8 t=1 s時流場的密度分布圖Fig.8 Distribution of density at t=1 s

圖9 t=1 s時流場的速度散度分布圖Fig.9 Distribution of velocity divergence at t=1s

3 結論

本文主要結論如下:①一般狀況下，三階MUSCL格式的數值耗散大于五階原始WENO格式.但當流場中出現間斷時，五階WENO格式為了穩定地捕捉間斷，其非線性系數會遠離最優權系數以致其耗散性大于理論上精度更低的三階MUSCL格式.②與基于原始變量的五階WENO格式相比，基于特征變量的五階WENO格式具有更大的數值耗散.③4種常用的高精度激波捕捉格式中，Map-WENO格式具有最小的數值耗散.④在激波-湍流干涉的問題中，目前廣受好評并被廣泛應用的4種高精度格式數值黏性都明顯偏大.因此，對于高精度高分辨率格式需要更深一步的研究.

References)

[1]Yan C.Computational fluid dynamic’s methods and applications[M].Beijing:Beihang University Press，2006

[2]Li X S，Xu J Z，Gu C W.Preconditioning method and engineering application of large eddy simulation[J].Science in China Series G:Physics，Mechanics and Astronomy，2008，51(6):667－677

[3]Johnsen E，Larsson J，Bhagatwala A V，et al.Assessment of highresolution methods for numerical simulations of compressible turbulence with shock waves[J].Journal of Computational Physics，2010，229(4):1213 －1237

[4]Toro E F.Riemann solvers and numerical methods for fluid dynamics[M].British:Springer，2009

[5]Kitamura K，Shima E.Towards shock-stable and accurate hypersonic heating computations:a new pressure flux for AUSM-family schemes[J].Journal of Computational Physics，2013，245:62 －83

[6]Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics，1996，126(1):202－228

[7]Zhang S H，Jiang S F，Shu C W.Improvement of convergence to steady state solutions of Euler equations with the WENO schemes[J].Journal of scientific Computing，2011，47(2):216 －238

[8]Qiu J X，Shu C W.On the construction，comparison and local characteristic decomposition forhigh ordercentralWENO schemes[J].Journal of Computational Physics，2002，183(1):187－209

[9]Henrick A K，Aslam T D，Powers J M.Mapped weighted essentially non-oscillatory schemes:achieving optimal order near critical points[J].Journal of Computational Physics，2005，207(2):542－567

[10]Kitamura K，Shima E.Towards shock-stable and accurate hypersonic heating computations:a new pressure flux for AUSM-family schemes[J].Journal of Computational Physics，2013，245:62－83

[11]Van Leer B.Towardsthe ultimate conservative difference scheme.V.A second-order sequel to Godunov’s method[J].Journal of Computational Physics，1997，135(2):229 －248

[12]Li X S，Gu C W.Mechanism of Roe-type schemes for all-speed flows and its application[J].Computers and Fluids，2013，86:56－70

[13]Fu D X，Ma Y W，Li X L，et al.Direct numerical simulation of the compressible turbulences[M].Beijing:Scientific Press，2010

[14]Zhang Z S，Cui G X.The principles and applications of the turbulent flows[M].Beijing:Tsinghua University Press，2008

[15]Zhang Z X，Dong Z N.The viscous fluid dynamics[M].Beijing:Tsinghua University Press，2011