高超聲速熱流計算湍流模型性能評估

張翔，閻超

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

高超聲速飛行器氣動熱的精確預(yù)測是計算流體力學(xué)(CFD)最具挑戰(zhàn)性的難題之一［1］.熱流是由黏性起主導(dǎo)作用的物理現(xiàn)象，它的計算精度與物理模型、數(shù)值格式、計算網(wǎng)格、收斂過程、熱流后處理等密切相關(guān)，這些多重因素的交錯影響導(dǎo)致了熱流計算的復(fù)雜性［2］.壁面熱流依賴溫度梯度在壁面上的精確計算，而高超聲速邊界層內(nèi)的溫度分布具有強(qiáng)非線性，要想獲得準(zhǔn)確的溫度分布需要在邊界層內(nèi)布置大量的計算網(wǎng)格［3］，以往的文獻(xiàn)均將網(wǎng)格雷諾數(shù)［2］(Recell)作為衡量網(wǎng)格因素的重要參數(shù)，并研究了保證熱流計算精度所需要的網(wǎng)格雷諾數(shù)條件［3-6］.對于計算格式和收斂判別方面，姜振華和閻超［2］研究了目前廣泛使用的兩類上風(fēng)格式:Roe的 FDS(Flux Difference Splitting)和AUSM(Advection Upstream Splitting Method)系列搭配不同限制器以及高階格式對熱流計算精度的影響;潘沙等［3］研究了熱流計算收斂性的判別標(biāo)準(zhǔn);李君哲和閻超［4］研究了幾類上風(fēng)格式對熱流計算精度的影響.

而高超聲速流動的雷諾數(shù)很大，在一般的飛行器上均會發(fā)生轉(zhuǎn)捩［7］，當(dāng)流動從層流轉(zhuǎn)變?yōu)橥牧鲿r，壁面摩阻和熱流會迅速增大，壁面熱流可以增大一個量級［8］.可以肯定的是，盡管 DNS(Direct Numerical Simulation)能夠給出簡單流動問題的計算結(jié)果，但對于工程常見的復(fù)雜外形，DNS高昂的計算成本和時間成本讓人望而生畏.本文結(jié)合已有文獻(xiàn)對計算網(wǎng)格、數(shù)值格式以及限制器選取的研究結(jié)果，分析了4種工程常用的湍流模型:BL，SA，k-ω，SST 對熱流計算精度的影響，并驗(yàn)證了計算精度較高的兩方程模型(k-ω與SST)搭配耗散性小的AUSMPW格式與混合限制器時對網(wǎng)格的依賴性.

另外，對于一般的吸氣式高超聲速飛行器，由于飛行高度較低，控制舵處流態(tài)大多已發(fā)展為湍流，此時舵前緣的氣動加熱異常嚴(yán)峻［9］.文獻(xiàn)［10］中詳細(xì)介紹了平板鈍舵模型的研究現(xiàn)狀，分析了不同因素對鈍舵熱流的影響.本文依據(jù)對湍流模型的評價結(jié)果，研究了不同后掠角對鈍舵熱流的影響，得到了鈍舵前緣最大熱流隨后掠角的變化趨勢.

1 數(shù)值計算方法

1.1 控制方程

本文熱流通過求解Reynolds平均N-S方程，守恒形式為［11］

各參數(shù)的意義詳見文獻(xiàn)［11］.為了使式(1)封閉，需要對式(1)中的雷諾應(yīng)力τij做出各種假設(shè).從對模式處理的出發(fā)點(diǎn)不同，一般可將湍流模式分為雷諾應(yīng)力模型和渦黏性模型兩類.受計算條件的約束，雷諾應(yīng)力模型計算量巨大，使其應(yīng)用范圍受到限制，在工程湍流問題中廣泛應(yīng)用的是渦黏性模型［12］.因此本文選取4種渦黏性湍流模型進(jìn)行計算.

壁面熱流qw由下式計算:

式中，Ma∞和Re∞為來流馬赫數(shù)和雷諾數(shù);μ為氣體動力黏性系數(shù);普朗特數(shù)Pr=0.72;兩者比熱之比γ=1.4;n為壁面法向.

1.2 數(shù)值離散方法

采用有限體積法求解雷諾平均可壓縮N-S方程和模型輸運(yùn)方程;無黏通量采用Roe的FDS格式和AUSMPW格式，并通過限制器［13］來抑制振蕩;黏性通量采用二階中心差分格式;時間推進(jìn)采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)隱式方法.

2 湍流模型對熱流計算的影響

對文獻(xiàn)［10］中的雙橢球風(fēng)洞實(shí)驗(yàn)進(jìn)行了數(shù)值模擬.該模型在低焓高超聲速風(fēng)洞中測得了詳實(shí)的氣動力及熱流數(shù)據(jù)，且其流動不受真實(shí)氣體效應(yīng)的影響.由于實(shí)驗(yàn)雷諾數(shù)較高，因此可將其作為驗(yàn)證高超聲速湍流模型性能的標(biāo)準(zhǔn)算例.

計算條件為:Ma∞=7.8，T∞=56 K，壁溫Tw=288 K，來流雷諾數(shù) Re∞=1.67 ×107m－1，攻角為0°.計算網(wǎng)格共267萬，壁面法向第1層網(wǎng)格雷諾數(shù)Recell為10.采用耗散較小的AUSMPW格式和混合限制器計算.

圖1給出了雙橢球?qū)ΨQ面上下表面處的壓力分布與實(shí)驗(yàn)值的比較，其中壓力分布通過來流靜壓Pinf無量綱化.可以看出4種湍流模型均可給出正確的壓力分布，并且與實(shí)驗(yàn)結(jié)果吻合良好.圖2給出了雙橢球?qū)ΨQ面上下表面的熱流沿軸向的分布.由圖可以看到，無論上下表面，BL模型與實(shí)驗(yàn)值相差較大，這是由于BL模型為零方程模型，它完全由當(dāng)時、當(dāng)?shù)氐钠骄鲃訁?shù)的函數(shù)來決定，對存在較大逆壓梯度以及分離再附等流動，該模式不能給出合理的結(jié)果.SA模型雖然能較好地模擬下表面的熱流分布，但對上表面兩橢球相貫處的激波-邊界層干擾模擬不足.兩種兩方程湍流模型:Wilcox k-ω模型和Menter SST模型對上下表面的熱流分布模擬較好，但k-ω模型對流動分離處的熱流計算較大.SST模型計算結(jié)果與實(shí)驗(yàn)吻合較好，這是由于SST模型是k-ε模型和k-ω模型的混合模型，同時保持了k-ω模型近壁面特性和k-ε在尾跡區(qū)的特性，其借鑒了J-K模型的思想，對雷諾剪切應(yīng)力進(jìn)行了經(jīng)驗(yàn)式的約束，計算效果較好［14］.

圖1 對稱面上下表面壓力分布Fig.1 Pressure distribution on upper and lower symmetrical surface

圖2 對稱面上下表面熱流分布Fig.2 Heat flow distribution on upper and lower symmetrical surface

3 網(wǎng)格雷諾數(shù)對熱流計算的影響

由以上分析可知，兩方程模型對熱流計算具有較高的精度.為考察湍流模型下熱流計算精度與網(wǎng)格雷諾數(shù)的關(guān)系，依舊選取雙橢球模型，取網(wǎng)格雷諾數(shù) Recell依次為 30，10，5 進(jìn)行計算.Recell為30的計算網(wǎng)格法向共80個點(diǎn)，為保證法向的過度小于1.1，Recell為10和5的計算網(wǎng)格在法向上依次增加到100和120.并且，這3種網(wǎng)格雷諾數(shù)下y+(平板自由來流結(jié)果預(yù)估)均小于1.依舊選取耗散較小的AUSMPW格式和混合限制器，湍流模型選取前文計算結(jié)果較好的Wilcox k-ω模型和Menter SST模型.

圖3給出了3種網(wǎng)格雷諾數(shù)Recell下k-ω模型對稱面上下表面熱流分布.可以看出對于上、下表面，3種不同Recell下熱流計算結(jié)果幾乎相同.這是由于:首先，無論哪種Recell，y+均小于1;其次，文獻(xiàn)［2］認(rèn)為，用耗散性較大的格式和限制器會導(dǎo)致網(wǎng)格依賴性較小、較穩(wěn)定的、略大于實(shí)驗(yàn)值的熱流結(jié)果，盡管計算格式和限制器選取了耗散性較小的格式，但是由于湍流模型較大的耗散性，導(dǎo)致熱流計算結(jié)果較穩(wěn)定，特別在y+＜1時，對網(wǎng)格雷諾數(shù)的敏感性較弱.而且可以看到，CFD計算結(jié)果總是略高于實(shí)驗(yàn)值.

圖3 k-ω模型上下表面熱流分布Fig.3 Heat flow distribution on upper and lower surface of k-ω model

圖4 SST模型上下表面熱流分布Fig.4 Heat flow distribution on upper and lower surface of SST model

圖4給出了3種網(wǎng)格雷諾數(shù)Recell下SST模型的計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對比.同樣，3種不同網(wǎng)格雷諾數(shù)下的熱流計算結(jié)果相近且均較實(shí)驗(yàn)值偏大.另外，相比k-ω模型，SST模型對雙橢球相貫處，流動分離后再附時的熱流峰值預(yù)測能力依舊優(yōu)秀.由圖3(a)與圖4(a)對比看出，k-ω模型對分離流動再附后的熱流峰值預(yù)測偏高，而SST模型的熱流峰值計算結(jié)果與實(shí)驗(yàn)值吻合非常好，這兩種湍流模型計算的熱流分布在y+＜1時均不隨網(wǎng)格雷諾數(shù)的變化而改變.圖5給出了在網(wǎng)格雷諾數(shù)10的情況下，SST模型計算的對稱面等馬赫線圖及雙橢球壁面極限流線圖.整體來看，流場等馬赫線清晰無抖動，可以清晰看到激波-邊界層干擾引起的流動分離，以及兩球相貫處的二次激波與頭部弓形激波的相交.

圖5 雙橢球等馬赫線圖和壁面極限流線圖Fig.5 Lines counter of equal Mach number and limit stream lines on wall of double ellipsoid model

4 后掠角對鈍舵熱流峰值的影響

通過第2節(jié)和第3節(jié)對氣動熱計算中湍流模型性能和網(wǎng)格雷諾數(shù)影響的研究，將其結(jié)論應(yīng)用于下文中平板鈍舵熱流的計算，以揭示舵前緣熱流峰值隨后掠角的變化趨勢.鈍舵前緣半徑為5 mm，高40 mm，后掠角從0°每隔10°～60°，共 7 種模型.計算條件為來流馬赫數(shù)Ma∞=5，自由來流參數(shù)取15 km高度處的美國標(biāo)準(zhǔn)大氣參數(shù)，壁溫Tw=300 K，攻角為0°，計算網(wǎng)格約210萬，保證Recell小于5，同時網(wǎng)格過度小于1.1.采用 Roe的FDS和AUSMPW格式搭配混合限制器，湍流模型選取Menter SST模型.

圖6給出了30°后掠角鈍舵壁面熱流等值線.可以看出舵前緣受熱嚴(yán)重，尤其是流動再附處.這表明鈍舵前緣駐點(diǎn)線的熱流是衡量飛行器熱環(huán)境的重要標(biāo)桿.圖7給出了在各后掠角情況下，前緣駐點(diǎn)線熱流隨z軸的變化趨勢.計算格式選取AUSMPW格式.可以看到駐點(diǎn)線熱流均在106量級，其中，分離區(qū)處駐點(diǎn)的熱流達(dá)到最大值，20°后掠角的熱流峰值約是流動非分離處的3倍.并且，熱流峰值的位置隨后掠角的增加而單調(diào)下降.這是由于后掠角的增大導(dǎo)致流動滯止點(diǎn)下移，進(jìn)而使駐點(diǎn)線最大熱流下移.

圖6 壁面熱流等值線Fig.6 Heat flow equal lines on the wall

圖7 不同后掠角鈍舵前緣熱流分布Fig.7 Heat flow distribution on leading edge of different angles

圖8給出了駐點(diǎn)線上熱流峰值隨后掠角的變化規(guī)律.計算分別采用Roe格式和AUSMPW格式配合SST模型，以及層流Roe格式.可以看出，湍流模型下，駐點(diǎn)線的熱流峰值隨后掠角的變化不是單調(diào)的，先隨后掠角的增大而增大，到后掠角大概為22°時達(dá)到最大值，然后熱流又隨后掠角的增大而逐漸下降.

圖8 不同后掠角下鈍舵熱流峰值Fig.8 Peak heat of blunt fin in different sweep angles

按文獻(xiàn)［2］的觀點(diǎn)，Roe格式的數(shù)值耗散要大于AUSM系列格式的，較大的耗散會導(dǎo)致計算的熱流值偏高，并且這一偏高的規(guī)律不隨計算網(wǎng)格的變化而改變.圖8中湍流模型計算結(jié)果表明，采用Roe格式計算的熱流在所有后掠角下均大于AUSM格式，在小后掠角情況下，這一表現(xiàn)更為顯著.

由于湍動能加熱的原因，湍流流態(tài)下熱流要遠(yuǎn)大于層流流態(tài)下的熱流.圖8也給出了湍流結(jié)果與層流結(jié)果的比較.可以看出，層流流態(tài)下熱流峰值隨后掠角的變化較平緩，隨后掠角的增大，熱流峰值呈現(xiàn)緩慢的增加趨勢，后掠角達(dá)到22°時，熱流達(dá)到最大值，隨后又緩慢減小.并且，湍流熱流峰值約為層流熱流峰值的2倍.

5 結(jié)論

1)對于雙橢球類存在激波-邊界層干擾流動的模型，4種湍流模型均能準(zhǔn)確模擬壁面的壓力分布，但BL模型無法正確模擬壁面的熱流分布，SA模型雖能較正確刻畫下表面的熱流分布，但是對上表面的激波-邊界層干擾處流動分離刻畫能力表現(xiàn)不足.兩方程模型表現(xiàn)優(yōu)秀，k-ω和SST模型的計算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，但是k-ω模型高估了流動再附處的熱流.SST模型總體表現(xiàn)優(yōu)異.

2)在y+＜1時，熱流計算對壁面法向網(wǎng)格的依賴性較小，網(wǎng)格雷諾數(shù)滿足10以內(nèi)即可得到與實(shí)驗(yàn)值吻合較好的熱流結(jié)果.

3)鈍舵駐點(diǎn)線熱流峰值隨后掠角呈現(xiàn)先增大后減小的趨勢.因此，在防熱設(shè)計中，可以考慮選擇較大的后掠角來減小鈍舵駐點(diǎn)線上的熱流峰值.

References)

［1］ Bertin J J，Cummings R M.Critical hypersonic aerothermodynamic phenomena［J］.Annual Review of Fluid Mechanics，2006，38:129-157.

［2］ 姜振華，閻超.高超聲速熱流高精度數(shù)值模擬方法［J］.北京航空航天大學(xué)學(xué)報，2011，37(12):1529-1533.Jiang Z H，Yan C.High order accurate methods in numerical simulation of hypersonic heat transfer［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37(12):1529-1533(in Chinese).

［3］ 潘沙，馮定華，丁國昊，等.氣動熱數(shù)值模擬中的網(wǎng)格相關(guān)性及收斂［J］.航空學(xué)報，2010，31(3):493-499.Pan S，F(xiàn)eng D H，Ding G H，et al.Grid dependency and convergence of hypersonic aerothermal simulation［J］.Acta Aeronautica et Astronautica Sinica，2010，31(3):493-499(in Chinese).

［4］ 李君哲，閻超.氣動熱CFD計算的格式效應(yīng)研究［J］.北京航空航天大學(xué)學(xué)報，2003，29(11):1022-2025.Li J Z，Yan C.Research on scheme effect of computational fluid dynamics in aerothermal［J］.Journal of Beijing University of Aeronautics and Astronautics，2003，29(11):1022-2025(in Chinese).

［5］ Klopfer G H，Yee H C.Viscous hypersonic shock/shock interaction on blunt line cowl lips，AIAA-1988-0233［R］.Reston:AIAA，1988.

［6］ Hoffmann K A，Siddiqui M S，Chiang S T.Difficulties associated with the heat flux computations of high speed flow by the Navier-Stokes equations，AIAA-1991-0457［R］.Reston:AIAA，1991.

［7］ Benard E，Cooper R K，Sidorenko A.Transitional and turbulent heat transfer of swept cylinder attachment line in hypersonic flow［J］.International Journal of Heat and Mass Transfer，2006，49(5-6):836-843.

［8］ 朱輝玉，王剛，孫泉華，等.高超聲速飛行器氣動熱數(shù)值模擬的幾個關(guān)鍵因素［C］//第三屆高超聲速科技學(xué)術(shù)會議會議文集.北京:中國力學(xué)學(xué)會，2010:0056-1-8.Zhu H Y，Wang G，Sun Q H，et al.Some critical factors in hypersonic aerothermal simulation［C］//Proceedings of Third Hypersonic Technology Academic Conference.Beijing:CSTAM，2010:0056-1-8(in Chinese).

［9］ 賀旭照，趙慧勇，樂嘉陵.吸氣式高超聲速飛行器氣動力氣動熱的數(shù)值模擬方法及應(yīng)用［J］.計算物理，2008，25(5):555-560.He X Z，Zhao H Y，Le J L.Aerodynamic force and heat of hypersonic laminar and turbulent flow［J］.Chinese Journal of Computational Physics，2008，25(5):555-560(in Chinese).

［10］ 李素循.激波與邊界層主導(dǎo)的復(fù)雜流動［M］.北京:科學(xué)出版社，2007:11-51，146-148.Li S X.The complex flow caused by shock/boundary interaction［M］.Beijing:Science Press，2007:11-51，146-148(in Chinese).

［11］ Wilcox D C.Turbulence modeling for CFD［M］.3rd ed.California:DCW Industries Inc，2006:243-249.

［12］ 閻超.計算流體力學(xué)方法與應(yīng)用［M］.北京:北京航空航天大學(xué)出版社，2006:221-223.Yan C.The methods and application of computational fluid dynamics［M］.Beijing:Beihang University Press，2006:221-223(in Chinese).

［13］ van Leer B.Towards the ultimate conservative difference scheme.V.A second-order sequel to Godunov’s method［J］.Journal of Computational Physics，1997，135(2):227-248.

［14］ 耿云飛，閻超，徐晶磊，等.高超聲速流動湍流模式評估［J］.北京航空航天大學(xué)學(xué)報，2011，37(8):907-911.Geng Y F，Yan C，Xu J L，et al.Turbulent modeling validation in hypersonic flow［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37(8):907-911(in Chinese).