屈 峰 閻 超 于 劍 仲崇巖

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

隨著計算流體力學(xué)(CFD，Computational Fluid Dynamics)需求的日益發(fā)展，人們對通量格式計算精度的要求也日益提高.1959年，俄羅斯科學(xué)家Godunov提出了基于Riemann精確解的一階計算格式[1].理論與實踐表明，它的解對于守恒型標(biāo)量方程而言，滿足熵條件且具有保單調(diào)性質(zhì)和總變差不增的總差不減(TVD，Total Variation Diminishing)性質(zhì)，并一定是物理解.因此該格式引起了廣泛的注意.自此之后，許多構(gòu)思精妙的通量格式雨后春筍般接連問世并極大促進了CFD的發(fā)展.而在這眾多格式中，目前世界上公認性能最為優(yōu)秀并被廣泛應(yīng)用的有Jameson的中心格式[2]、Roe的通量差分(FDS，F(xiàn)lux Difference Splitting)格式[3]、van Leer的矢通量分裂(FVS，F(xiàn)lux Vector Splitting)格式[4]、Liou 的對流上風(fēng)分裂(AUSM，Advection Upstream Splitting Method)格式[5]及在其構(gòu)造基礎(chǔ)上衍生的一系列格式如AUSM+[6]，SLAU2[7]，AUSM+UP[8]等.盡管在實際應(yīng)用中上述格式均表現(xiàn)出較為良好的特性，其依然因為各自構(gòu)造思想中的“bug”而距離人們理想中的“完美”格式相距甚遠.如van Leer的FVS格式無法精確捕捉接觸間斷，黏性分辨率較差;Roe的FDS格式在計算高速流動時會出現(xiàn)“粉刺”現(xiàn)象且對定常無黏通量的計算不滿足焓守恒條件;Jameson的中心格式和AUSM+UP需要根據(jù)實際問題人工設(shè)置經(jīng)驗參數(shù)等[9].2010年，Kitamura等人通過一系列的數(shù)值實驗得出結(jié)論[10]:目前已有的格式在高超聲速計算中依然無法得到準(zhǔn)確可靠的氣動熱結(jié)果.而影響其計算精度的因素有格式的激波穩(wěn)定性、總焓守恒性、黏性分辨率.Kim等人通過對奇偶失聯(lián)問題的理論分析，將Roe格式進行改進得到了RoeM類格式[11]，計算結(jié)果表明，RoeM類格式在保持Roe格式黏性分辨率高、魯棒性強等優(yōu)點的基礎(chǔ)上避免了“粉刺”現(xiàn)象的出現(xiàn)并保持了定常無黏通量計算的總焓守恒.但是格式構(gòu)造過程中，Kim等人僅簡單地通過計算相鄰網(wǎng)格格心處的壓強比來判斷該處是否有激波.這在較為簡單的流動計算中尚不會造成影響，但在較為復(fù)雜的流動模擬中會因非物理情況的出現(xiàn)而導(dǎo)致其可靠性降低.為了克服這一問題，本文利用Jiang提出的光滑指示因子法來判斷激波位置[12]，在 RoeM 格式的基礎(chǔ)上構(gòu)造了 RoeMW1，RoeMW2格式.下文的相關(guān)分析將表明，這兩種格式均滿足以下特點:

1)無需通過設(shè)置經(jīng)驗性參數(shù)就可保證其通量函數(shù)的激波穩(wěn)定性;

2)黏性分辨率高，接觸間斷捕捉精確;

3)在高速定常流的計算中保證總焓守恒;

4)膨脹區(qū)計算時魯棒性較強，無膨脹激波出現(xiàn).

1 控制方程與計算方法

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

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

定義空間單元為

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

其中

1.2 RoeMW類格式的構(gòu)造

迎風(fēng)類通量格式的數(shù)值質(zhì)量通量均可由式(3)表示.

2000 年，Liou[13]發(fā)現(xiàn)式(3)中 D(p)的存在會導(dǎo)致Roe格式在計算超聲速流動時容易出現(xiàn)“粉刺”現(xiàn)象.文獻[14]也得到類似的結(jié)論.而文獻[11]通過理論分析得出D(ρ)的大小亦影響著此格式的激波穩(wěn)定性.為了平衡質(zhì)量通量在數(shù)值計算中壓強梯度對密度耗散的影響，D(p)可取為以下形式:

其中Bi采用文獻[12]提出的形式，即

式中，s(x)為插值函數(shù);k表示精度.實際計算表明，k的取值大小對最終計算結(jié)果影響不大.因此，考慮到計算效率的問題，在這里取k=2.即

后面的算例結(jié)果表明，在大部分情況下，該處理足以抑制“粉刺”現(xiàn)象的發(fā)生.但在某些情況如高速非定常流動的計算中，它可能因為壓強擾動與密度擾動之間的平衡難以保證而出現(xiàn)激波不穩(wěn)定現(xiàn)象.為了解決這一問題，本文在上述處理的基礎(chǔ)上再利用函數(shù)f來加強對密度擾動的耗散.其具體形式如下:

眾所周知，無黏定常流計算中Roe格式不能保證總焓守恒的特性.根據(jù)Jameson在文獻[15]中給出的結(jié)論，焓守恒格式的連續(xù)性方程和能量方程需滿足如下形式:

同時，為了消除Roe格式在實際計算時經(jīng)常出現(xiàn)的“膨脹激波”現(xiàn)象，本文借鑒HLLEM[16]格式以及AUSM+[6]的構(gòu)造思想將上式稍作修正，其具體形式如下:

綜上所述，本文所得格式的最終形式如下所述:

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

2.1 一維接觸間斷問題

本問題的初始條件為

庫蘭特數(shù)(CFL，Courant，F(xiàn)riedrichs，and Lewy)取0.85.為了實現(xiàn)對上述通量格式進行“純粹”地分析，下文所有的算例結(jié)果均由一階格式計算得到.

從圖1～圖3可以清晰地看出，3種格式在捕捉到的接觸間斷基本一致，且速度壓強分布曲線無振蕩.這表明RoeMW1和RoeMW2格式同Roe格式一樣具有較高的黏性分辨率.

圖1 一維接觸間斷的密度分布曲線

圖2 一維接觸間斷的壓強分布曲線

圖3 一維接觸間斷的速度分布曲線

2.2 無黏圓柱超音速繞流

1988年，文獻[17]研究發(fā)現(xiàn)Roe格式在計算鈍頭體激波時會產(chǎn)生“粉刺”現(xiàn)象(圖4a).而FVS格式則無此現(xiàn)象發(fā)生.因此，本算例是判斷格式激波穩(wěn)定性的經(jīng)典算例.在本文中，取來流馬赫數(shù)為8，網(wǎng)格量為102×171，CFL數(shù)為5.時間格式統(tǒng)一選用無條件穩(wěn)定的LUSGS方法.

圖4 無黏圓柱超音速繞流等壓強線圖

從圖4可以看出，RoeMW1與RoeMW2格式在保證了激波分辨率高的同時均可有效地避免“粉刺”現(xiàn)象的發(fā)生.圖5則表明RoeMW1與RoeMW2格式在計算該問題時保持總焓守恒.而這一點在預(yù)測氣動熱等氣動性能參數(shù)時非常重要.

圖5 對稱線上激波參數(shù)顯示圖

2.3 雙馬赫反射

該算例最初由Woodward和colella[18]給出，并已成為格式研究的標(biāo)準(zhǔn)算例.其具體描述如下:計算域為[0，4]×[0，1]，其中在下邊界[1/6，4]的位置為滑移壁面條件.初始時刻，在計算域中給定一個與x軸成60°夾角的馬赫數(shù)為10的激波，該激波從x=1/6，y=0的位置一直延伸到上邊界，并向右側(cè)快速移動.下邊界[0，1/6]的區(qū)域以及左邊界均給定激波的波后條件，右邊界給定出流條件，上邊界根據(jù)激波每一時刻所處的位置分別給定激波的波前和波后條件.計算時間推進到t=0.2.采用的計算網(wǎng)格為960×240.

圖6給出了Roe，RoeM，RoeMW1以及RoeMW2這4種格式計算得到的密度等值線圖.從圖上可清晰地看到Roe格式在計算正激波時會發(fā)生激波不穩(wěn)定現(xiàn)象.RoeM格式雖然較之于Roe格式有所改進，但依然出現(xiàn)比較明顯的激波不穩(wěn)定現(xiàn)象.相比于RoeM格式，RoeMW1對Roe格式有了更大的改善，但依然無法完全抑制激波不穩(wěn)定現(xiàn)象的發(fā)生.這也與上文所述的僅僅利用函數(shù)f來控制壓強耗散系數(shù)有時并不能完全抑制激波不穩(wěn)定現(xiàn)象發(fā)生的結(jié)論一致.同時，從圖6中RoeMW2的計算結(jié)果可看出，在RoeMW1格式基礎(chǔ)上又增加了對密度擾動耗散的RoeMW2格式，則成功地抑制了激波不穩(wěn)定現(xiàn)象的發(fā)生.

圖6 雙馬赫反射密度等值線圖

3 結(jié)論

本文通過利用WENO格式中的光滑指示因子來判斷流場間斷，在RoeM格式的基礎(chǔ)上得到了RoeMW1，RoeMW2格式.計算結(jié)果表明:一般情況下，構(gòu)造相對簡單的RoeMW1格式在擁有Roe格式黏性分辨率高的優(yōu)點同時，能夠抑制激波不穩(wěn)定現(xiàn)象的發(fā)生并且保證無黏定常計算時的總焓守恒.但在模擬較為復(fù)雜的流動結(jié)構(gòu)時，雖然比RoeM格式更加穩(wěn)定，但RoeMW1格式僅對壓強擾動采取控制的方式依然無法完全抑制激波不穩(wěn)定現(xiàn)象的發(fā)生.而同時控制壓強擾動與密度擾動的RoeMW2格式則可以在保持RoeMW1格式上述所有優(yōu)點的基礎(chǔ)上解決此問題.因此在實際復(fù)雜流動模擬中，它比現(xiàn)有的通量格式(如Roe，RoeM等格式)具有更廣闊的應(yīng)用前景.

References)

[1]Toro E F.Riemann solvers and numerical methods for fluid dynamics[M].Berlin:Springer-Verlag，2009

[2]Jameson A，Schmidt W，Turkel E.Numerical solutions of the Euler equations by finite volume methods using Runge-Kuttta timestepping schemes[R].AIAA 81-1259，1981

[3]Roe P L.Approximate Riemann solvers，parameter vectors and difference schemes[J].J Comp Phys，1981，43:357－372

[4]van Leer B.Flux vector splitting for the Euler equations[J].Lecture Notes in Physics，1982，170:507－512

[5]Liou M S，Steffen Jr C J.A new flux splitting scheme[J].Journal of Computational Physics，1993，107(1):23－39

[6]Liou M S.Progress toward an improved CFD method:AUSM+[R].AIAA 95-1701-CP，1995

[7]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

[8]Liou M S.A further development of the AUSM_scheme，towards robust and accurate solutions for all speeds[R].AIAA 2003-4116，2003

[9]閻超.計算流體力學(xué)方法及應(yīng)用[M].2版.北京:北京航空航天大學(xué)出版社，2006

Yan Chao.Computational methods and applications[M].2nd ed.Beijing:Beihang University Press，2006(in Chinese)

[10]Kitamura K，Shima E，Nakamura Y，et al.Evaluation of Euler fluxes for hypersonic heating computations[J].AIAA Journal，2010，48(4):763－776

[11]Kim K H，Kim C，Rho O H.Cures for the shock instability，development of a shock-stable Roe scheme[J].Journal of Computational Physics，2003，185(2):342－374

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

[13]Liou M S.Mass flux schemes and connection to shock instability[J].Journal of Computational Physics，2000，160(2):623－648

[14]Kitamura K，Shima E.Improvements of simple low-dissipation AUSM against shock instabilities in consideration of interfacial speed of sound[C]//Proceedings of ECCOMAS CFD.Libson:[s.n.]，2010:14－17

[15]Jameson A.Analysis and design of numerical schemes for gas dynamics，2:artificial diffusion and discrete shock structure[J].International Journal of Computational Fluid Dynamics，1995，5(1/2):1－38

[16]Einfeldt B，Munz C D，Roe P L，et al.On Godunov-type methods near low densities[J].J Comput Phys，1991，92(2):273－295

[17]Peery K M，Imlay S T.Blunt-body flow simulations[R].AIAA88-2904，1988

[18]Woodward P，Colella P.The numerical simulation of two-dimensional fluid flow with strong shocks[J].Journal of Computational Physics，1984，54(1):115－173