不同熱邊界對超音速機翼氣動熱數值分析的影響

趙璇，孫智，呂潤民，孫建紅

(南京航空航天大學 飛行器環境控制與生命保障工業和信息化部重點實驗室，南京 210016)

0 引 言

隨著科技的不斷進步，對飛行器的結構及性能的需求大幅提高，但是在超音速巡航狀態下，機翼蒙皮表面氣流因受到強烈壓縮與摩擦，極高的動能將轉化為熱能，對機翼結構及內部燃油產生嚴重影響，進而可能影響其飛行性能、隱身性及燃油系統的可靠性。且機翼蒙皮表面的溫度場較流場更為復雜，外流場的熱流與蒙皮結構的溫度相互影響，因此對于氣動加熱問題的研究需要對流-固-熱多物理場進行耦合計算。

國外對飛行器氣動加熱方面的研究開展較早。例如，1958年J.R.Davidson等[1]對馬赫數為2的多翼機翼的氣動加熱問題進行試驗研究，分析了氣動加熱下熱應力及剛度的關系；1987年A.R.Wieting[2]應用8 ft(1 ft=0.304 8 m)高溫風洞(8-foot High Temperature Tunnel，簡稱8’ HTT)進行激波干擾對圓管前緣高超音速氣動加熱影響的試驗；1989年P.Dechaumphai等[3]對高超音速條件下流場-氣動熱-固體結構的相互作用進行有限元分析；S.A.Berry等[4]通過耦合算法對X-34模型進行了氣動熱模擬研究，結果顯示在層流與湍流狀態下與實驗值的吻合度較好。國內，朱增起[5]采用直接數值模擬方法計算了高超音速鈍頭體圓錐壁面在不同邊界條件下駐點及非駐點熱流，指出對流換熱公式較傅里葉熱傳導公式對于熱流計算更加合理，但未涉及壁面耦合效應；Chen Fang等[6]提出了一種用于流體-熱-結構建模與分析的時間自適應多物理耦合方法，且解決了高超聲速氣動熱力學與有限元分析軟件的結合，為分析高超聲速機翼結構的熱-結構-振動響應提供了有效的預測；李佳偉等[7]提出了基于有限體積的氣動熱-結構傳熱一體化計算方法，提高了計算效率，但未考慮內壁面邊界條件對于氣動熱耦合計算的影響；王洋等[8]建立了一種快速求解氣動熱參數的降階模型，熱邊界傳遞基于換熱系數的線性近似方法可以減少因壁面溫度變化所需的流場迭代，但忽略了壁面溫度對換熱系數的影響。

本文在分析外流場與外壁面的耦合換熱基礎上，對內壁面設置不同的邊界條件，得到外壁面的熱流或溫度分布，研究內壁面邊界條件對氣動熱耦合數值計算的影響。預期為研究機翼熱防護及燃油箱的熱負荷提供一定的參考基礎。

1 數值方法

1.1 流場求解方程

采用三維粘性可壓縮流動的Navier-Stokes(N-S)方程，在笛卡爾坐標系下的形式為

(1)

式中：F、G、H分別為三個坐標方向上的通量矢量；下標c為對流通量；下標v為粘性通量。

具體形式為

(2)

(3)

式中：u，v，w分別為x，y，z三個方向上的速度分量。

本文采用Spalart-Allmaras一方程湍流模型對控制方程進行封閉，該湍流模型對高速流動下近壁面的激波捕捉有較好的能力。并基于有限體積法求解方程，對流通量采用Roe-FDS格式進行離散，粘性通量采用二階迎風格式進行離散。

1.2 溫度場求解方程

當不考慮固體形變時，結構內部主要為熱量傳遞。結構內部溫度場的三維非穩態導熱微分方程在笛卡爾坐標系下的形式為

(4)

式中：ρ為固體材料的密度；Cp為材料的比熱容；λ為材料的導熱系數；qv為固體內熱源密度。

1.3 耦合傳熱求解方法

耦合傳熱發生在流場與結構場的交界面上，傳統方法是將結構體表面設為絕熱或等溫壁，而流場與結構場之間是相互作用與影響的關系，流場對結構體進行加熱的同時，結構內溫度場時刻發生變化，進而影響流場向結構場內部的傳熱量，因此，需要將流場與結構場之間進行實時的信息傳遞，才能更加接近于真實的物理情況。

邊界上進行的熱量傳遞式為

(5)

式中：qf為流場通過界面向結構場傳熱的熱流；k為流場中流體的導熱率。

根據能量守恒，在忽略表面輻射影響的前提下，流場向結構場傳熱的熱流應與結構場內部導熱熱流相等，由此建立了二者的耦合關系。即給定t=0 s時遠場的初始參數及壁面的溫度分布，通過求解流場控制方程得到t時刻后壁面的熱流分布，作為條件計算結構內部新的溫度場，將得到的壁面新溫度分布作為流場下一時間步計算的初始條件，由此進行迭代計算，直至收斂。

耦合傳熱根據時間步長的選取又分為緊耦合與松耦合。緊耦合采用流場的特征時間作為耦合計算的共同時間，松耦合是先計算出流場的穩定解并作為耦合的邊界條件，再采用結構的時間步長作為耦合的特征時間進行耦合計算。但緊耦合的計算時間相對較長，效率較松耦合低，而松耦合在滿足計算結果正確的前提下可以提高計算速度。因此，本文采用松耦合的計算方法。

2 圓管氣動熱計算

本文選取文獻[2]中的試驗條件對圓管進行氣動熱計算。因為重點關注圓管前部分的流場及氣動加熱情況，所以將試驗模型簡化為試驗模型的1/4。二維結構網格劃分為383×251，如圖1所示，圓管外直徑及壁厚分別為：D=0.076 2 m，d=0.012 7 m，基于圓管直徑的第一層網格無量綱高度為1.3×10-6。

圖1 結構網格劃分Fig.1 Structural meshing

邊界條件設定為壓力遠場、無滑移等溫壁(Tw=294.44 K)及對稱邊界條件。計算穩定后外流場的密度云圖如圖2所示，與試驗符合較好。

圖2 流場密度與試驗對比Fig.2 Comparison of flow field density with experimental

圓管駐點前中心線上溫度分布如圖3所示，可以進一步得出激波位置(x≈-0.055 m)、波后溫度(T=2 263 K)與文獻[3]中的結果較為接近。壓力與表面熱流密度分布如圖4～圖5所示，其中，p0為駐點處壓力，q0為駐點處熱流密度，θ為圓管外表面與對稱面的夾角。

圖3 圓管駐點前中心線上溫度分布Fig.3 Temperature distribution on the centerline before the stagnation point of the tube

圖4 t=0 s時表面無量綱壓力分布Fig.4 Surface pressure distribution at t=0 s

圖5 t=0 s時表面無量綱熱流密度分布Fig.5 Heat flux distribution at t=0 s

從圖4～圖5可以看出：計算結果與試驗結果擬合較好，其中駐點熱流密度的理論值為4.8×105W/m2，試驗值為6.7×105W/m2，計算結果為5.2×105W/m2，介于理論值與試驗值之間，考慮到試驗氣體為甲烷燃燒后產物，與理想空氣的熱力學性質并非完全相同，因此計算與試驗的結果會存在一定的差異。

計算2 s后圓管內部溫度云圖如圖6所示，文獻[3]中駐點溫度為388.88 K，本文計算結果為393.3 K，誤差為1.5%，結果吻合較好。

計算2 s后圓管外表面的熱流密度分布如圖7所示。

圖6 t=2 s時圓管內溫度云圖Fig.6 Temperature of the tube at t=2 s

圖7 t=2 s時表面無量綱熱流密度分布Fig.7 Heat flux distribution at t=2 s

從圖7可以看出：氣動加熱在圓管外壁產生的熱量隨著時間逐漸向圓管內部傳導，使圓管壁面溫度升高，從而使駐點熱流密度有5%的降低。

為了進一步分析結構場傳熱，得到了0～6 s內不同時刻圓管駐點處熱流密度和外壁的溫度分布，如圖8～圖9所示。

圖8 不同時刻圓管駐點處熱流密度分布Fig.8 Heat flux distribution at the stagnation point of the tube at different times

圖9 不同時刻外壁面溫度分布Fig.9 Outer wall temperature distribution at different times

從圖8～圖9可以看出：隨著時間推進，因氣動熱效應向壁面不斷傳熱，駐點處熱流密度值逐漸降低，而壁面受到加熱溫度不斷升高，且駐點溫度升高較快。

綜合得到圓管外流場和結構場的計算結果，可以得到該方法對于超音速繞流及氣動熱計算具有有效性，可運用此方法對超音速飛行下機翼的氣動熱效應進行數值計算。

3 不同熱邊界對機翼氣動熱影響

為驗證該網格對機翼外流場計算的正確性，本文針對文獻9中的實驗條件對稱雙弧翼型進行數值計算。來流馬赫數為2.13，攻角為-10°，計算得到的機翼表面壓力系數分布如圖11所示，與實驗結果擬合良好，證明了網格和方法對于超音速機翼繞流計算的可靠性。

圖10 三種翼型輪廓圖Fig.10 Three airfoil shapes

圖11 對稱雙弧翼型表面壓力系數分布Fig.11 Pressure coefficient distribution of the symmetrical arc airfoil

3.1 不同翼型流場及氣動熱的差異

在計算工況相同的情況下，三種翼型表面壓力系數分布如圖12所示，因為三種翼型均為對稱翼型，所以上下翼面的壓力系數分布表現出一致性。

圖12 三種翼型的壓力系數Fig.12 Pressure coefficient of three airfoils

由駐點處的壓力系數可以得出：雙弧翼型作為典型的超音速翼型，其尖前緣較其他兩種鈍前緣有效減小了波阻，且壓力系數近似線性分布。將厚度不同的兩個NACA翼型的壓力系數進行對比，可以看出：厚度越大，機翼前緣曲率越小，激波過渡較慢，而沿著弦長方向NACA0012的曲率大于NACA0006的曲率，隨之壓力系數降低較快。

為進一步分析不同翼型表面氣動熱的差異，內壁面選取無滑移絕熱壁進行分析。計算得到三種翼型表面溫度分布，如圖13所示。

圖13 翼型表面溫度分布Fig.13 Surface temperature distribution of airfoils

翼型表面溫度分布相似于壓力系數分布，雙弧翼型溫度分布近似線性，而NACA翼型前緣存在小部分高溫區，后沿弦向溫度分布較均勻，且NACA翼型的駐點溫度高于雙弧翼型的駐點溫度，這是由于NACA翼型在前緣附近產生的正激波強度較雙弧翼型產生的斜激波強度大，即來流受到激波壓縮后動能消耗更多，進而轉化的熱能也更多，當壁面為絕熱時，翼型前緣的溫度也就更高。而真實情況的機翼壁面無法做到絕熱，因此絕熱條件下的機翼壁面溫度分布僅可為熱防護材料的選取提供一定的理論依據。真實情況下，外流場產生的氣動熱對壁面進行加熱，改變了壁面溫度，而新的壁溫又將影響外流場向壁面傳導的熱流密度，同時內部環境因素與內壁面的換熱也將使壁面溫度發生變化，因此，需要考慮內壁面邊界條件對氣動熱效應的影響。

3.2 內壁面不同熱邊界對翼型氣動熱影響

為了研究內壁面不同邊界條件對機翼氣動熱的影響，針對NACA0012翼型選取三種不同邊界條件，分別為Case1(內壁面絕熱)、Case2(內壁面等溫)及Case3(內壁面對流換熱)，Case2與Case3均基于機翼內部滿油作為條件進行設置，其中燃油物性參數參考文獻10。考慮到高空的低溫效應，假設壁面及燃油的初始溫度均為283.15 K，且將實際內部為滿油工況作為Case4進行對比。三種邊界條件及真實情況下的翼型表面溫度和熱流密度分布如圖14～圖15所示。

圖14 不同邊界條件翼型表面溫度分布Fig.14 Surface temperature distribution of airfoils under different boundary conditions

圖15 不同邊界條件翼型表面熱流密度分布Fig.15 Heat flux distribution of airfoils under different boundary conditions

真實狀態下，靠近壁面的燃油受熱導致內部燃油出現溫差而流動，使得壁面溫度及熱流密度存在小幅度的波動；絕熱條件忽略傳熱導致翼型平均溫度高于真實情況約13 K，由于等溫條件保持較低的壁面溫度，導致與外流場換熱時的溫差較大，Case2的平均熱流密度為3 941.5 W/m2，Case3的平均熱流密度為1 966.3 W/m2，Case4的平均熱流密度為1 933.3 W/m2；等溫條件忽略了熱壁面效應，平均熱流密度約為真實狀態的2倍，邊界條件為對流換熱較真實狀態下平均熱流密度的誤差為1.7%。因此，內壁面邊界條件為對流換熱時與真實情況最為接近，且較真實情況下計算量更小。

3.3 三維機翼燃油箱滿油狀態下的氣動熱計算

考慮內壁面為對流換熱的基礎上，對三維機翼油箱在滿油狀態下進行氣動熱數值計算。機翼蒙皮厚度為5 mm，蒙皮材料為鋁，導熱系數202.4 W/(m·k)。基于弦長的無量綱雷諾數為1.09×107，三維結構網格徑向、周向和展向的網格數分別為51、406和101。機翼表面壓力分布云圖如圖16所示，機翼前緣的壓力最高，沿徑向逐漸降低。

圖16 機翼表面壓力云圖Fig.16 Surface pressure cloud map of the wing

為了進一步分析三維機翼表面的氣動熱環境，沿展向截取三個截面，如圖17所示：中心截面A0、距離中心截面左右兩側各1 m的截面A1和A2。

圖17 機翼不同截面位置Fig.17 Different cross-section positions of the wing

不同截面機翼表面的溫度與熱流密度分布如圖18～圖19所示。

圖18 機翼不同截面溫度分布Fig.18 Temperature distribution of different cross sections of the wing

圖19 機翼不同截面熱流密度分布Fig.19 Heat flux distribution of different cross sections of the wing

從圖18～圖19可以看出：機翼平行于翼跟的各截面前10%所受到的氣動熱較明顯。A0截面尾緣處的壓力云圖如圖20所示，尾緣的高壓導致溫度和熱流密度出現突增。

圖20 A0截面尾緣壓力云圖Fig.20 Trailing edge pressure cloud map on A0 section

4 結 論

(1) 圓管與翼型的超音速繞流數值模擬驗證結果與實驗數據吻合較好，本文采用的耦合方法可以對翼型繞流進行有效的氣動熱耦合計算。

(2) 雙弧翼型雖然可以有效降低激波阻力，但其前半翼型處于迎風面，能量易于積累，翼型表面溫度及熱流密度均高于NACA0006翼型和NACA0012翼型。

(3) 在所選取的工況下，翼型內壁面為絕熱時翼型表面溫度高于真實狀態下約13 K，邊界條件為等溫時翼型的平均熱流密度約為真實狀態下的2倍，而對流換熱的邊界條件下平均熱流密度為1 966.3 W/m2，與真實狀態下的1 933.3 W/m2誤差為1.7%，所以邊界條件為對流換熱與真實狀態最為接近。三維機翼平行于翼根的各截面前10%屬于氣動熱較嚴重區，且駐點熱流密度最高可達4 200 W/m2。