張翔，閻超

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

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

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

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

1 數(shù)值計算方法

1.1 控制方程

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

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

壁面熱流qw由下式計算:

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

1.2 數(shù)值離散方法

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

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

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

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

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

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

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

圖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種網格雷諾數(shù)Recell下SST模型的計算結果與實驗數(shù)據的對比.同樣，3種不同網格雷諾數(shù)下的熱流計算結果相近且均較實驗值偏大.另外，相比k-ω模型，SST模型對雙橢球相貫處，流動分離后再附時的熱流峰值預測能力依舊優(yōu)秀.由圖3(a)與圖4(a)對比看出，k-ω模型對分離流動再附后的熱流峰值預測偏高，而SST模型的熱流峰值計算結果與實驗值吻合非常好，這兩種湍流模型計算的熱流分布在y+＜1時均不隨網格雷諾數(shù)的變化而改變.圖5給出了在網格雷諾數(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é)對氣動熱計算中湍流模型性能和網格雷諾數(shù)影響的研究，將其結論應用于下文中平板鈍舵熱流的計算，以揭示舵前緣熱流峰值隨后掠角的變化趨勢.鈍舵前緣半徑為5 mm，高40 mm，后掠角從0°每隔10°～60°，共 7 種模型.計算條件為來流馬赫數(shù)Ma∞=5，自由來流參數(shù)取15 km高度處的美國標準大氣參數(shù)，壁溫Tw=300 K，攻角為0°，計算網格約210萬，保證Recell小于5，同時網格過度小于1.1.采用 Roe的FDS和AUSMPW格式搭配混合限制器，湍流模型選取Menter SST模型.

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

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

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

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

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

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

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

5 結論

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

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

3)鈍舵駐點線熱流峰值隨后掠角呈現(xià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］.北京航空航天大學學報，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ù)值模擬中的網格相關性及收斂［J］.航空學報，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計算的格式效應研究［J］.北京航空航天大學學報，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ù)值模擬的幾個關鍵因素［C］//第三屆高超聲速科技學術會議會議文集.北京:中國力學學會，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ù)值模擬方法及應用［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］ 李素循.激波與邊界層主導的復雜流動［M］.北京:科學出版社，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］ 閻超.計算流體力學方法與應用［M］.北京:北京航空航天大學出版社，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］.北京航空航天大學學報，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).