崔樹鑫 高 琳 高 歌

(北京航空航天大學 能源與動力工程學院，北京100191)

精確捕捉激波是計算流體動力學(CFD，Computational Fluid Dynamics)研究的重要內(nèi)容.從20世紀60年代CFD的逐漸形成，到后來各種格式的不斷興起及其在計算機上的成功應用，幾乎大多數(shù)研究都是圍繞激波、間斷等內(nèi)容展開的.無論哪種格式，其主要任務都是解決單元交界面無粘數(shù)值通量的確定問題.作為迎風格式的兩類典型代表——矢通量分裂格式(FVS，F(xiàn)lux Vector Splitting)［1］和通量差分分裂格式 (FDS，F(xiàn)lux Difference Splitting)［2］雖然已能較好地捕捉激波和間斷，并且得到了廣泛應用，但是其通量確定的復雜性，向三維推廣的非精確性以及會出現(xiàn)非物理解等一系列問題，仍困擾著眾多CFD工作者.

CE/SE方法源于20世紀90年代，由參考文獻［3］提出.它將空間量與時間量統(tǒng)一處理，通過交錯的時空網(wǎng)格分布，巧妙地解決數(shù)值通量的確定難題.即引入包含交界面的解元，交界面的數(shù)值通量僅由與已知求解點相關(guān)的泰勒展開式得到，既保證了界面數(shù)值通量的唯一性，又維持了時空單元的守恒性.因此，它無需方向分裂，也不需要求得黎曼解.更為突出的優(yōu)點是，基于加權(quán)平均思想的空間導數(shù)求解，可以靈活地控制數(shù)值耗散.近來，新發(fā)展的 CNIS格式［4］，克服了克朗數(shù)(CFL，Courant-Friedrichs-Lewy)遠遠小于1時，增加的數(shù)值耗散對物理解的污染.迄今為止，CE/SE方法已成功應用于氣動聲學、化學反應和磁流體［5－8］等領(lǐng)域.而關(guān)于翼型流場的激波捕捉問題，國內(nèi)外學者研究甚少.

文獻［9］采用CE/SE方法對無激波的翼型繞流問題進行了數(shù)值研究，并取得了很好的結(jié)果.但是，由于守恒元和解元布置方式的特殊性，在每個時間層內(nèi)只有一半空間網(wǎng)格點參與計算.若要分辨諸如激波等間斷現(xiàn)象，必須提供兩倍的空間網(wǎng)格點，大大降低了CE/SE的高分辨率特性，不適合在翼型流場中捕捉激波.為此，本文對守恒元和解元重新定義，將計算點放在網(wǎng)格節(jié)點和單元中心上，在一個時間步內(nèi)交錯計算單元、節(jié)點上的流場變量.以此為基礎(chǔ)，對激波翼型流場進行數(shù)值模擬.通過對數(shù)值結(jié)果的討論與分析，說明新的CE/SE方法捕捉激波的有效性和精度.

1 控制方程及新的CE/SE方法

二維無量綱歐拉方程各分量表達式為

式中，［u1，u2，u3，u4］= ［ρ，ρu，ρv，ρE］;［f1，f2，f3，f4］= ［ρu，ρu2+p，ρuv，(ρE+p)u］;［g1，g2，g3，g4］=［ρv，ρvu，ρv2+p，(ρE+p)v］;m=1，2，3，4分別對應連續(xù)方程、x方向動量方程、y方向動量方程和能量方程;ρ是氣體密度;u是x方向速度分量;v是y方向速度分量;p是壓力;E是單位比能量.

圖1 空間網(wǎng)格

因為 (fm)Q*，(gm)Q*，(fmx)Q*，(gmx)Q*，(fmy)Q*，(gmy)Q*，(fmt)Q*，(gmt)Q*是(um)Q*，(umx)Q*，(umy)Q*，(umt)Q*的函數(shù)，由微分形式的歐拉方程(1)可知，(umt)Q*是(fmx)Q*和(gmy)Q*的函數(shù).故而，對于任意求解點Q*，其獨立的變量僅為(um)Q*，(umx)Q*和(umy)Q*.

令x1=x，x2=y，x3=t作為三維歐幾里德空間E3的坐標.對方程(1)在守恒元CE(Q)上積分，并運用高斯散度定理，可得

式中，S是面 A1B1A2B2A3B3A4B4的面積;S1，S2，S3，S4，(x1，y1)，(x2，y2)，(x3，y3)，(x4，y4)分別為面 A1B1QB4，面 A2B2QB1，面 A3B3QB2和 面A4B4Q B3的面積和質(zhì)心坐標.(x11，y11)，(x21，y21)，(x12，y12)，(x22，y22)，(x13，y13)，(x23，y23)，(x14，y14)，(x24，y24)，λ11，λ21，λ12，λ22，λ13，λ23，λ14，λ24，(n11x，n11y)，(n21x，n21y)，(n12x，n12y)，(n22x，n22y)，(n13x，n13y)，(n23x，n23y)，(n14x，n14y)，(n24x，n24y)分別為側(cè)邊 B4A1，A1B1，B1A2，A2B2，B2A3，A3B3，B3A4，A4B4的中心坐標、邊長和外法向矢量分量.

式(3)就是新的CE/SE求解守恒變量的表達式.相應地，可以推導出以單元質(zhì)心為存儲點的計算公式.此外，CE/SE方法還需要計算求解點的空間導數(shù)，該計算方法與文獻［10］一致，這里不再累述.至此，形成了兩套分別基于單元中心和網(wǎng)格節(jié)點對應求解點的計算公式.在計算過程中，將單元中心、網(wǎng)格節(jié)點上的求解點分別放置相鄰的半時間層上，通過交錯的方式相互求解，從而，組成一套完整的CE/SE方法計算體系.

新方法增加了儲存變量節(jié)點，相對于原CE/SE方法［10］，在一定程度上增加了計算復雜度.另外，在本文中還要用到CNIS格式，關(guān)于該格式的推導過程以及克朗數(shù)的計算方法參見文獻［11］.

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

本文采用兩個算例對新CE/SE方法在翼型流場中捕捉激波的能力進行分析和驗證，同時，與原方法的計算結(jié)果進行了比較.根據(jù)文獻［8］的結(jié)論，在以下算例中，均采用CNIS格式和當?shù)貢r間步長法(CFL=0.5)，加權(quán)平均因子取為2.

2.1 NACA0012有攻角跨音流動

采用C型網(wǎng)格，網(wǎng)格數(shù)為265X45(翼型周向布置147個節(jié)點)，計算域的上下邊界與翼型中線的法向距離為20倍翼型弦長，下游邊界距離翼型前緣40倍翼型弦長，翼型弦長取為1，后緣采用文獻［12］的處理方法(由于翼型后緣非零厚度，將后緣延長至x=1.0089處，并將原坐標同時除以新弦長1.008 9).圖3是網(wǎng)格在翼型附近的分布圖.計算工況為 Ma=0.8，α=1.25°.

圖3 NACA0012翼型網(wǎng)格分布圖

圖4給出了翼型表面壓力系數(shù)的分布情況.可以看出，新方法能夠很好地捕捉到翼型表面的強激波和弱激波，與文獻［12］吻合良好，適合模擬含有激波的翼型流場.進一步觀察發(fā)現(xiàn)，由于加權(quán)平均因子的選擇比較合理，兩種方法所引入的數(shù)值粘性均有效地抑制了激波前后的過沖和過膨脹.除了激波和后緣附近區(qū)域，兩計算結(jié)果與文獻的相對誤差小于2%.后緣處的誤差主要是網(wǎng)格密度不足所致.值得注意的是，新方法捕捉上表面強激波需要跨越3～4個空間網(wǎng)格點，而原方法需要5～6個;另外，新方法能夠成功地捕捉到下表面弱激波，而原方法的弱激波被抹平得很厲害，表現(xiàn)為一較弱的壓力陡升.究其原因，主要有以下兩個方面:①原方法采用交錯的時空網(wǎng)格布置，計算點分布在網(wǎng)格單元中心上，在每個時間步上僅有一半單元中心點參與計算，而采用節(jié)點/單元交互的新方法，計算點分布在網(wǎng)格節(jié)點和單元中心上，每一個時間步內(nèi)，全部單元中心點都參與計算，相當于增加了計算點的網(wǎng)格密度，因而，分辨諸如激波等間斷的能力強于原方法.②原方法的邊界條件賦值于與邊界毗鄰的虛單元中心處，邊界附近內(nèi)單元的守恒元跨越邊界，全局守恒區(qū)域大于計算域.而新方法邊界條件直接設置在邊界節(jié)點上，邊界附近內(nèi)單元的守恒元邊界剛好與計算域邊界重合，全局守恒區(qū)域即為計算域本身.因此，從守恒性角度出發(fā)，后者更好于前者.為更精確地反映CE/SE方法捕捉激波的能力，表1列出了新方法在翼型附近激波后馬赫數(shù)的計算值，并與由朗金-雨貢紐(R-H)關(guān)系得到的理論值進行對比.可以看出，新方法計算結(jié)果的相對誤差不超過5%.

表1 激波后馬赫數(shù)對比

2.2 RAE2822有攻角跨音流動

圖4 算例1的表面壓力系數(shù)對比曲線

為了進一步考察新方法的普適性，對典型的超臨界翼型RAE2822的激波流動特征進行數(shù)值模擬.超臨界翼型是一種為提高臨界馬赫數(shù)而采取的特殊翼型，能夠使機翼在接近音速時阻力劇增的現(xiàn)象推遲發(fā)生，因此，精確捕捉激波位置對超臨界翼型的研究具有十分重要的意義.本文采用C型網(wǎng)格，網(wǎng)格數(shù)為265×45(翼型周向布置147個節(jié)點)，翼型弦長取為1，計算域的上下邊界與翼型中線的法向距離為20倍翼型弦長，下游邊界距離翼型前緣20倍翼型弦長.圖5是網(wǎng)格在翼型附近的分布圖.計算工況為 Ma=0.75，α=3.0°.

圖5 RAE2822翼型網(wǎng)格分布圖

圖6和圖7分別是翼型表面壓力分布情況和新方法計算得到的馬赫數(shù)分布云圖.從中可以看出，本文的計算方法很好地模擬出了氣流在翼型上表面前緣膨脹加速，而后相對平緩流動，直至加速出現(xiàn)激波、壓力驟升的整個過程.兩種方法計算結(jié)果與文獻［12］吻合良好，在翼型前緣上表面，原CE/SE方法的結(jié)果有些膨脹過度，在后緣區(qū)域，由于本文的網(wǎng)格密度遠遠小于參考文獻的網(wǎng)格密度(320×64，O型網(wǎng)格)，兩種方法的結(jié)果均與文獻存在差距.而在其他區(qū)域，其相對誤差均小于1%.加權(quán)平均引入的數(shù)值耗散對激波附近解的過沖和過膨脹現(xiàn)象有明顯的抑制作用.新方法捕捉的激波跨越3～4個空間網(wǎng)格點，而原方法則需5～6個，進一步證實了新方法在激波分辨率方面的改進是有效的.從結(jié)果也可以看到，在不采用CNIS格式的情況下，新CE/SE方法的激波分辨率仍高于原方法.表2列出了激波位置計算值的對比情況，可以看出，兩種結(jié)果都十分接近文獻［12］，相對誤差不超過1%.

圖6 算例2的表面壓力系數(shù)對比曲線

圖7 新CE/SE方法計算得到的馬赫數(shù)分布云圖

表2 激波位置

3 結(jié)論

1)通過對不同跨音翼型流場的激波捕捉算例進行驗證，說明新CE/SE方法對翼型流場捕捉激波具有適用性.

2)對原CE/SE方法的改進是有效的，可以明顯提高捕捉激波的分辨率.

3)CE/SE方法能夠有效地抑制激波前后的數(shù)值震蕩，避免發(fā)生過沖和過膨脹現(xiàn)象.

References)

［1］Van Leer B.Flux-vector splitting for the Euler equations［J］.Lecture Notes in Physics ，1982，170:507－512

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

［3］Chang S C.The method of space-time conservation element and solution element—a new approach for solving the navier-stokes and euler equations［J］.Journal of Computational Physics，1995，119(2):295－324

［4］Chang S C.Courant number insensitive CE/SE schemes［R］.AIAA-2002-3890，2002

［5］劉敏，吳克啟.基于非結(jié)構(gòu)網(wǎng)格CE/SE算法圓柱繞流氣動聲場模擬［J］.工程熱物理學報，2009(2):227－229

Liu Min，Wu Keqi.The aero-acoustics simulation of flow around a cylinder using unstructured CE/SE scheme［J］.Journal of Engineering Thermophysics，2009(2):227－229(in Chinese)

［6］凌文輝，趙書苗，劉建文，等.時空守恒CE/SE數(shù)值方法在凹槽燃燒流動中的應用研究［J］.推進技術(shù)，2010(1):34－41

Ling Wenhui，Zhao Shumiao，Liu Jianwen，et al.Numerical investigations of the cavity combustion flow using the CE/SE method［J］.Journal of Propulsion Technology，2010(1):34－41(in Chinese)

［7］Ji Z，Zhou Y F.A comparison study of three CESE schemes in MHD simulation［J］.ChinesePhysicsLetters，2010，27(8):085201

［8］孫孔倩，樊未軍，楊茂林.CE/SE法模擬爆震波點火及傳播過程［J］.航空動力學報，2008，23(2):244－250

Sun Kongqian，F(xiàn)an Weijun，Yang Maolin.Numerical simulation of detonation ignition and propagation progress by CE/SE method［J］.Journal of Aerospace Power，2008，23(2):244－250(in Chinese)

［9］崔樹鑫，韓玉琪，高歌.二維翼型繞流的CE/SE方法［J］.北京航空航天大學學報，2011，37(1):21－24

Cui Shuxin，Han Yuqi，Gao Ge.Space-time conservation element and solution element method(CE/SE)applied to flows around 2D airfoil［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37(1):21－24(in Chinese)

［10］Zhang Z C，Yu S，Chang S C.A space-time conservation element and solution element method for solving the two-and three-dimensional unsteady Euler equations using quadrilateral and hexahedral meshes［J］.Journal of Computational Physics，2002，175(1):168－199

［11］Yen J C，Wagner D A.Computational aeroacoustics using a simplified courant number insensitive CE/SE method［R］.AIAA-2005-2820，2005

［12］Yoshihara H，Sacher P.Test cases for inviscid flow field methods［R］.Advisory Group for Aerospace Research and Development(AGARD)-AR-211，1985