孫迪，閻超，于劍，屈峰，華俊

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

伴隨著計算流體力學(xué)(CFD)技術(shù)的日益發(fā)展，人們對高精度格式的需求愈發(fā)強烈.譬如研究具有多尺度的湍流問題以及對一些涉及氣動聲學(xué)、熱流、摩阻及非定常流動問題的計算，需要采用具有多維分辨能力的高階格式才能得到較為滿意的結(jié)果［1-2］.

目前，CFD界對高階格式的研究方興未艾，ENO (Essentially Non-oscillatory)、 WENO［3］(Weighted Essentially Non-Oscillatory)、ENN(Essentially Non-oscillatory，No-free-parameter)、緊致格式及譜方法等格式雨后春筍般應(yīng)運而生.但它們都因為計算量過大或是魯棒性較差而無法被廣泛應(yīng)用于復(fù)雜飛行器的流動模擬當(dāng)中.此外，上述格式的構(gòu)造思想都是基于一維插值，而實際流動是多維的，因此在計算中它們無法判斷多維流動現(xiàn)象，特別是多維流動間斷(如與網(wǎng)格斜交的斜激波)，以致分辨率降低.為了克服這一缺陷，近十年來，國內(nèi)外的學(xué)者發(fā)展了一系列基于多維流動的計算格式，比如:模化波速法、旋轉(zhuǎn)通量法等［4-6］.然而實踐表明，這些方法雖然在計算精度方面比傳統(tǒng)方法有顯著改善，但魯棒性較差，在多維間斷處的保單調(diào)性不足，特別是在模擬高超聲速時流場時，強激波附近常常產(chǎn)生劇烈的非物理振蕩［7］.另外，這些方法的求解效率較差，因而未被廣泛使用.

與上述方法相比，Kim和Noh等［8-9］提出的多維限制器MLP［10］有較好的魯棒性和收斂性，將其與傳統(tǒng)的TVD［11］高階插值格式相結(jié)合可以更好地抑制多維振蕩.由于 Roe［12］格式是被廣泛認可的空間離散方法，因此本文將結(jié)合Roe格式的MLP限制器方法與一些常用數(shù)值計算方法的計算結(jié)果進行對比，發(fā)現(xiàn)該多維限制器有如下優(yōu)點:①相較于傳統(tǒng)二階TVD格式，在效率相當(dāng)?shù)那闆r下有更高的分辨率.②在多維間斷處能夠有效抑制振蕩.③保證強魯棒性及收斂性的前提下有較高的求解精度.④算法簡單，容易實現(xiàn).

1 計算方法

1.1 控制方程及空間離散方法

二維N-S方程(Navier-Stokes equations)的守恒形式［13］為

式中流通矢量為

1.2 MUSCL(TVD)限制器的一般形式

MUSCL 格式的具體形式［14-15］為

式中，φ(rL)φ(rR)為限制器;Φ表示通量矢量.

1.3 高階多維限制器的構(gòu)造

高精度多維限制器的主要構(gòu)造方法如下［16］.

首先構(gòu)造高階TVD插值:

因此，為了抑制在間斷處的振蕩，高階插值也應(yīng)滿足TVD限制條件:

式中3階限制器的表達形式為

5階限制器的表達形式為

其次，構(gòu)造多維限制函數(shù)，由單調(diào)性條件，界面通量函數(shù)值應(yīng)在空間相鄰的8個函數(shù)值范圍之間，即

在充分考慮了4個方向(上下左右)的θj(如圖1所示)后，得出多維限制條件如下:

綜上所述，得到高階限制器的表述形式:

圖1 格點值與格心值分布示意圖Fig.1 Distribution schematic of cell-vertex value and cell-center value of grid

2 算例及結(jié)果分析

為進一步分析驗證該多維限制器的性能，本文做了一系列數(shù)值實驗，包含無黏非定常問題及有黏的激波邊界層干擾問題.簡便起見，下文算例中3階多維限制器用MLP3表示，3階WENO格式用WENO3表示，5階WENO格式用WENO5表示.時間格式方面，非定常問題采用3階龍格庫塔法，定常問題采用無條件穩(wěn)定的LUSGS方法.

2.1 Sod 問題

計算域為［0，1］，網(wǎng)格點數(shù)為 100，初始條件為

計算推進到t=0.2 s中止，此時，流場中包含一個激波、一個接觸間斷和一個膨脹波.

比較密度分布曲線(圖2)可得，minmod限制器耗散較大，WENO3格式次之，WENO5格式及MLP3限制器的分辨率較高，而局部放大圖更清晰地表明了這一結(jié)論.

圖2 密度分布曲線Fig.2 Density distribution curves

2.2 二維渦流動問題

渦流動問題是一個典型的非定常多維流動問題.初始流場是在均勻流場上添加了一個等熵渦，該渦的強度為Γ=5.計算區(qū)域為10×10的方形區(qū)域，網(wǎng)格點數(shù)為80×80.擾動的形式如下.

式中，γ為比熱比;(x0，y0)為渦核的坐標(biāo).

理論上該無黏渦的渦核壓強隨著時間的推移是保持不變的.

初始流場密度分布如圖3所示.圖4給出了漩渦轉(zhuǎn)動100個無量綱時間后，3種不同插值方法得到的沿AB線(圖3)的密度分布曲線.結(jié)果表明，在多維光滑的情況下，minmod限制器有較大的耗散;與WENO3格式相比，MLP3限制器有較好的多維分辨能力.

圖3 初始流場密度分布Fig.3 Density distribution of initial vortex flow field

圖4 沿AB線的壓力分布Fig.4 Pressure distribution along line AB

圖5 熵隨時間變化的曲線Fig.5 Variation curves of entropy with change of time

分析熵隨時間的變化曲線(圖5)可知:superbee限制器計算得到的熵值一直非物理地減小，渦量變強且無收斂趨勢.而其他格式的計算結(jié)果表明:minmod限制器、WENO3及MLP3限制器的熵增依此減小，精度依此增高.

取2套不同密度的網(wǎng)格(40×40，100×100)并分別采用3種限制器(minmod，WENO3，MLP3)計算進行對比.圖6給出了各格式在不同密度網(wǎng)格中計算時，渦核壓強隨時間推進的變化情況.從中可以看出，3種限制器的耗散性均隨網(wǎng)格量減小有所增大，但MLP3限制器的數(shù)值耗散增加最小，網(wǎng)格收斂性最佳.

圖6 不同格式的網(wǎng)格收斂性比較Fig.6 Comparison of grid convergence with different schemes

2.3 激波邊界層干擾

本算例描述的物理問題為一斜激波入射到平板邊界層上，使得邊界層局部產(chǎn)生分離，并在平板下游再附［17］.流動參數(shù)如下:Ma=2，T=117 K，Re=296000，激波與平板之間的夾角為32.598°.網(wǎng)格單元數(shù)是200×200.計算時CFL數(shù)取5，壁面條件為無滑移絕熱壁.

圖7給出了各限制器計算得到的壓力等值線圖，圖8則給出了圖7(a)所示虛線處的壓強分布曲線.從中可以看出:minmod限制器及WENO3格式耗散較大;相比之下，MLP3限制器和WENO5格式可以更為清晰地分辨出入射激波及通過分離區(qū)的反射激波、膨脹扇區(qū)和再附激波等流動結(jié)構(gòu).但在圖8局部放大圖中，WENO5格式由于耗散過小，再附激波后出現(xiàn)了非物理的虛假振蕩，魯棒性較差.

圖9表明MLP3限制器與WENO5格式計算所得的壁面壓強分布相近，均與實驗值吻合較好，較為正確地預(yù)測了分離區(qū)的大小，而其他格式的耗散較大，求得的分離區(qū)過小.

圖7 壓力等值線圖Fig.7 Pressure contours

圖8 壓強分布曲線Fig.8 Pressure distribution curves

圖9 壁面壓強分布曲線Fig.9 Pressure distribution curves along wall surface

3 結(jié)論

本文基于無黏渦算例及激波邊界層干擾等算例，研究分析了高精度多維限制器MLP3的相關(guān)特性.通過將其與傳統(tǒng)TVD限制器及高階WENO格式進行比較，得到主要結(jié)論如下:

1)MLP3限制器與常見二階精度TVD限制器相比具有明顯的優(yōu)勢:與minmod限制器相比，MLP3有較高的精度，且通過多維限制函數(shù)，它避免了過多的耗散，可以較為精確地捕捉到激波等間斷;與superbee限制器相比，MLP3具有良好的保單調(diào)性以避免非物理解的產(chǎn)生.

2)二維渦流動算例表明:較之于高階WENO格式及傳統(tǒng)TVD限制器，MLP3限制器更容易實現(xiàn)，且在花費更小計算量的同時保持強魯棒性及高精度.

3)激波邊界層干擾算例表明:超聲速有黏流動計算時，MLP3限制器的求解精度高于同階WENO格式，與5階WENO格式結(jié)果相似.因此MLP3具有較高的黏性分辨率和保單調(diào)特性.

References)

［1］ Zhang H X.On problems to develop physical analysis in CFD［C］//Proceedings of the Fourth Asian Computational Fluid Dynamics Conference.Chengdu:IEEE，2000:3-19.

［2］ 楊建龍，劉猛.限制器對高超聲速表面熱流數(shù)值模擬的影響［J］.北京航空航天大學(xué)學(xué)報，2014，40(3):417-421.Yang J L，Liu M.Influence of limiters on numerical simulation of heating distributions for hypersonic bodies［J］.Journal of Beijing University of Aeronautics and Astronautics，2014，40(3):417-421(in Chinese).

［3］ Shu C W，Osher S.Efficient implementation of essentially nonoscillatory shock-capturing schemes［J］.Journal of Computational Physics，1988，77(2):439-471.

［4］ Roe P L.Discrete models for the numerical analysis of time-dependent multidimensional gas dynamics［J］.Journal of Computational Physics，1986，63(2):458-476.

［5］ Lacor C，Hirsch C.Genuinely upwind algorithms for multidimensional Euler equations［J］.AIAA Journal，1992，30(1):56-63.

［6］ Deconinck H，Paillere H，Struijs R，et al.Multidimensional upwind schemes based on fuctuation-splitting for systems of conservation laws［J］.Computational Mechanics，1993，11(5-6):323-340.

［7］ 屈峰，閻超，于劍，等，高精度激波捕捉格式的性能分析［J］.北京航空航天大學(xué)學(xué)報，2014，40(8):1085-1089.Qu F，Yan C，Yu J，et al.Assessment of shock capturing methods for numerical simulations of compressible turbulence with shock waves［J］.Journal of Beijing University of Aeronautics and Astronautics，2014，40(8):1085-1089(in Chinese).

［8］ Kim K H，Kim C.Accurate，efficient and monotonic numerical methods for multi-dimensional compressible flows，Part II:multidimensional limiting process［J］.Journal of Computational Physics，2005，208(2):570-615.

［9］ Noh S J，Lee K R，Park J H O，et al.An accurate and efficient calculation of high enthalpy flows using a high order new limiting process［J］.Journal of the Korean Society for Industrial and Applied Mathematics，2011，15(1):67-82.

［10］ Kang H M，Kim K H，Lee D H.A new approach of a limiting process for multi-dimensional flows［J］.Journal of Computational Physics，2010，229(19):7102-7128.

［11］ Harten A.High resolution schemes for hyperbolic conservation laws［J］.Journal of Computational Physics，1983，49(3):357-393.

［12］ Roe P L.Approximate Riemann solvers，parameter vectors and difference schemes［J］.Journal of Computational Physics，1981，43(2):357-372.

［13］ 閻超.計算流體力學(xué)方法及應(yīng)用［M］.北京:北京航空航天大學(xué)出版社，2006:15-25.Yan C.Computational fluid dynamic’s methods and applications［M］.Beijing:Beihang University Press，2006:15-25(in Chinese).

［14］ van Leer B.Toward the ultimate conservative difference scheme［J］.Journal of Computational Physics，1997，135(2):229-248.

［15］ Hirsch C.Numerical computation of international and external flows:Volume 2［M］.Hoboken，NJ:John Wiley＆ Sons Publish，1990:121-156.

［16］ Park J S，Chang T K，Kim C.Higher-order multi-dimensional limiting strategy for correction procedure via reconstruction［C］//52nd Aerospace Sciences Meeting.Maryland:AIAA，2014:2014-0772.

［17］ Knight D.RTO WG 10-Test cases for CFD validation of hypersonic flight，AIAA-2002-0433［R］.Reston:AIAA，2002.