高精度多維限制器的性能分析

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

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

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

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

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

1 計(jì)算方法

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階限制器的表達(dá)形式為

5階限制器的表達(dá)形式為

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

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

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

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

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

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

2.1 Sod 問題

計(jì)算域?yàn)椋?，1］，網(wǎng)格點(diǎn)數(shù)為 100，初始條件為

計(jì)算推進(jìn)到t=0.2 s中止，此時(shí)，流場(chǎng)中包含一個(gè)激波、一個(gè)接觸間斷和一個(gè)膨脹波.

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

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

2.2 二維渦流動(dòng)問題

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

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

理論上該無黏渦的渦核壓強(qiáng)隨著時(shí)間的推移是保持不變的.

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

圖3 初始流場(chǎng)密度分布Fig.3 Density distribution of initial vortex flow field

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

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

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

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

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

2.3 激波邊界層干擾

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

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

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

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

圖8 壓強(qiáng)分布曲線Fig.8 Pressure distribution curves

圖9 壁面壓強(qiáng)分布曲線Fig.9 Pressure distribution curves along wall surface

3 結(jié)論

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

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

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

3)激波邊界層干擾算例表明:超聲速有黏流動(dòng)計(jì)算時(shí)，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］ 楊建龍，劉猛.限制器對(duì)高超聲速表面熱流數(shù)值模擬的影響［J］.北京航空航天大學(xué)學(xué)報(bào)，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é)報(bào)，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］ 閻超.計(jì)算流體力學(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.