趙璇，孫智，呂潤(rùn)民，孫建紅

(南京航空航天大學(xué) 飛行器環(huán)境控制與生命保障工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，南京 210016)

0 引 言

隨著科技的不斷進(jìn)步，對(duì)飛行器的結(jié)構(gòu)及性能的需求大幅提高，但是在超音速巡航狀態(tài)下，機(jī)翼蒙皮表面氣流因受到強(qiáng)烈壓縮與摩擦，極高的動(dòng)能將轉(zhuǎn)化為熱能，對(duì)機(jī)翼結(jié)構(gòu)及內(nèi)部燃油產(chǎn)生嚴(yán)重影響，進(jìn)而可能影響其飛行性能、隱身性及燃油系統(tǒng)的可靠性。且機(jī)翼蒙皮表面的溫度場(chǎng)較流場(chǎng)更為復(fù)雜，外流場(chǎng)的熱流與蒙皮結(jié)構(gòu)的溫度相互影響，因此對(duì)于氣動(dòng)加熱問(wèn)題的研究需要對(duì)流-固-熱多物理場(chǎng)進(jìn)行耦合計(jì)算。

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

本文在分析外流場(chǎng)與外壁面的耦合換熱基礎(chǔ)上，對(duì)內(nèi)壁面設(shè)置不同的邊界條件，得到外壁面的熱流或溫度分布，研究?jī)?nèi)壁面邊界條件對(duì)氣動(dòng)熱耦合數(shù)值計(jì)算的影響。預(yù)期為研究機(jī)翼熱防護(hù)及燃油箱的熱負(fù)荷提供一定的參考基礎(chǔ)。

1 數(shù)值方法

1.1 流場(chǎng)求解方程

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

(1)

式中：F、G、H分別為三個(gè)坐標(biāo)方向上的通量矢量；下標(biāo)c為對(duì)流通量；下標(biāo)v為粘性通量。

具體形式為

(2)

(3)

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

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

1.2 溫度場(chǎng)求解方程

當(dāng)不考慮固體形變時(shí)，結(jié)構(gòu)內(nèi)部主要為熱量傳遞。結(jié)構(gòu)內(nèi)部溫度場(chǎng)的三維非穩(wěn)態(tài)導(dǎo)熱微分方程在笛卡爾坐標(biāo)系下的形式為

(4)

式中：ρ為固體材料的密度；Cp為材料的比熱容；λ為材料的導(dǎo)熱系數(shù)；qv為固體內(nèi)熱源密度。

1.3 耦合傳熱求解方法

耦合傳熱發(fā)生在流場(chǎng)與結(jié)構(gòu)場(chǎng)的交界面上，傳統(tǒng)方法是將結(jié)構(gòu)體表面設(shè)為絕熱或等溫壁，而流場(chǎng)與結(jié)構(gòu)場(chǎng)之間是相互作用與影響的關(guān)系，流場(chǎng)對(duì)結(jié)構(gòu)體進(jìn)行加熱的同時(shí)，結(jié)構(gòu)內(nèi)溫度場(chǎng)時(shí)刻發(fā)生變化，進(jìn)而影響流場(chǎng)向結(jié)構(gòu)場(chǎng)內(nèi)部的傳熱量，因此，需要將流場(chǎng)與結(jié)構(gòu)場(chǎng)之間進(jìn)行實(shí)時(shí)的信息傳遞，才能更加接近于真實(shí)的物理情況。

邊界上進(jìn)行的熱量傳遞式為

(5)

式中：qf為流場(chǎng)通過(guò)界面向結(jié)構(gòu)場(chǎng)傳熱的熱流；k為流場(chǎng)中流體的導(dǎo)熱率。

根據(jù)能量守恒，在忽略表面輻射影響的前提下，流場(chǎng)向結(jié)構(gòu)場(chǎng)傳熱的熱流應(yīng)與結(jié)構(gòu)場(chǎng)內(nèi)部導(dǎo)熱熱流相等，由此建立了二者的耦合關(guān)系。即給定t=0 s時(shí)遠(yuǎn)場(chǎng)的初始參數(shù)及壁面的溫度分布，通過(guò)求解流場(chǎng)控制方程得到t時(shí)刻后壁面的熱流分布，作為條件計(jì)算結(jié)構(gòu)內(nèi)部新的溫度場(chǎng)，將得到的壁面新溫度分布作為流場(chǎng)下一時(shí)間步計(jì)算的初始條件，由此進(jìn)行迭代計(jì)算，直至收斂。

耦合傳熱根據(jù)時(shí)間步長(zhǎng)的選取又分為緊耦合與松耦合。緊耦合采用流場(chǎng)的特征時(shí)間作為耦合計(jì)算的共同時(shí)間，松耦合是先計(jì)算出流場(chǎng)的穩(wěn)定解并作為耦合的邊界條件，再采用結(jié)構(gòu)的時(shí)間步長(zhǎng)作為耦合的特征時(shí)間進(jìn)行耦合計(jì)算。但緊耦合的計(jì)算時(shí)間相對(duì)較長(zhǎng)，效率較松耦合低，而松耦合在滿足計(jì)算結(jié)果正確的前提下可以提高計(jì)算速度。因此，本文采用松耦合的計(jì)算方法。

2 圓管氣動(dòng)熱計(jì)算

本文選取文獻(xiàn)[2]中的試驗(yàn)條件對(duì)圓管進(jìn)行氣動(dòng)熱計(jì)算。因?yàn)橹攸c(diǎn)關(guān)注圓管前部分的流場(chǎng)及氣動(dòng)加熱情況，所以將試驗(yàn)?zāi)Ｐ秃?jiǎn)化為試驗(yàn)?zāi)Ｐ偷?/4。二維結(jié)構(gòu)網(wǎng)格劃分為383×251，如圖1所示，圓管外直徑及壁厚分別為：D=0.076 2 m，d=0.012 7 m，基于圓管直徑的第一層網(wǎng)格無(wú)量綱高度為1.3×10-6。

圖1 結(jié)構(gòu)網(wǎng)格劃分Fig.1 Structural meshing

邊界條件設(shè)定為壓力遠(yuǎn)場(chǎng)、無(wú)滑移等溫壁(Tw=294.44 K)及對(duì)稱邊界條件。計(jì)算穩(wěn)定后外流場(chǎng)的密度云圖如圖2所示，與試驗(yàn)符合較好。

圖2 流場(chǎng)密度與試驗(yàn)對(duì)比Fig.2 Comparison of flow field density with experimental

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

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

圖4 t=0 s時(shí)表面無(wú)量綱壓力分布Fig.4 Surface pressure distribution at t=0 s

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

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

計(jì)算2 s后圓管內(nèi)部溫度云圖如圖6所示，文獻(xiàn)[3]中駐點(diǎn)溫度為388.88 K，本文計(jì)算結(jié)果為393.3 K，誤差為1.5%，結(jié)果吻合較好。

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

圖6 t=2 s時(shí)圓管內(nèi)溫度云圖Fig.6 Temperature of the tube at t=2 s

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

從圖7可以看出：氣動(dòng)加熱在圓管外壁產(chǎn)生的熱量隨著時(shí)間逐漸向圓管內(nèi)部傳導(dǎo)，使圓管壁面溫度升高，從而使駐點(diǎn)熱流密度有5%的降低。

為了進(jìn)一步分析結(jié)構(gòu)場(chǎng)傳熱，得到了0～6 s內(nèi)不同時(shí)刻圓管駐點(diǎn)處熱流密度和外壁的溫度分布，如圖8～圖9所示。

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

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

從圖8～圖9可以看出：隨著時(shí)間推進(jìn)，因氣動(dòng)熱效應(yīng)向壁面不斷傳熱，駐點(diǎn)處熱流密度值逐漸降低，而壁面受到加熱溫度不斷升高，且駐點(diǎn)溫度升高較快。

綜合得到圓管外流場(chǎng)和結(jié)構(gòu)場(chǎng)的計(jì)算結(jié)果，可以得到該方法對(duì)于超音速繞流及氣動(dòng)熱計(jì)算具有有效性，可運(yùn)用此方法對(duì)超音速飛行下機(jī)翼的氣動(dòng)熱效應(yīng)進(jìn)行數(shù)值計(jì)算。

3 不同熱邊界對(duì)機(jī)翼氣動(dòng)熱影響

為驗(yàn)證該網(wǎng)格對(duì)機(jī)翼外流場(chǎng)計(jì)算的正確性，本文針對(duì)文獻(xiàn)9中的實(shí)驗(yàn)條件對(duì)稱雙弧翼型進(jìn)行數(shù)值計(jì)算。來(lái)流馬赫數(shù)為2.13，攻角為-10°，計(jì)算得到的機(jī)翼表面壓力系數(shù)分布如圖11所示，與實(shí)驗(yàn)結(jié)果擬合良好，證明了網(wǎng)格和方法對(duì)于超音速機(jī)翼繞流計(jì)算的可靠性。

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

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

3.1 不同翼型流場(chǎng)及氣動(dòng)熱的差異

在計(jì)算工況相同的情況下，三種翼型表面壓力系數(shù)分布如圖12所示，因?yàn)槿N翼型均為對(duì)稱翼型，所以上下翼面的壓力系數(shù)分布表現(xiàn)出一致性。

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

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

為進(jìn)一步分析不同翼型表面氣動(dòng)熱的差異，內(nèi)壁面選取無(wú)滑移絕熱壁進(jìn)行分析。計(jì)算得到三種翼型表面溫度分布，如圖13所示。

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

翼型表面溫度分布相似于壓力系數(shù)分布，雙弧翼型溫度分布近似線性，而NACA翼型前緣存在小部分高溫區(qū)，后沿弦向溫度分布較均勻，且NACA翼型的駐點(diǎn)溫度高于雙弧翼型的駐點(diǎn)溫度，這是由于NACA翼型在前緣附近產(chǎn)生的正激波強(qiáng)度較雙弧翼型產(chǎn)生的斜激波強(qiáng)度大，即來(lái)流受到激波壓縮后動(dòng)能消耗更多，進(jìn)而轉(zhuǎn)化的熱能也更多，當(dāng)壁面為絕熱時(shí)，翼型前緣的溫度也就更高。而真實(shí)情況的機(jī)翼壁面無(wú)法做到絕熱，因此絕熱條件下的機(jī)翼壁面溫度分布僅可為熱防護(hù)材料的選取提供一定的理論依據(jù)。真實(shí)情況下，外流場(chǎng)產(chǎn)生的氣動(dòng)熱對(duì)壁面進(jìn)行加熱，改變了壁面溫度，而新的壁溫又將影響外流場(chǎng)向壁面?zhèn)鲗?dǎo)的熱流密度，同時(shí)內(nèi)部環(huán)境因素與內(nèi)壁面的換熱也將使壁面溫度發(fā)生變化，因此，需要考慮內(nèi)壁面邊界條件對(duì)氣動(dòng)熱效應(yīng)的影響。

3.2 內(nèi)壁面不同熱邊界對(duì)翼型氣動(dòng)熱影響

為了研究?jī)?nèi)壁面不同邊界條件對(duì)機(jī)翼氣動(dòng)熱的影響，針對(duì)NACA0012翼型選取三種不同邊界條件，分別為Case1(內(nèi)壁面絕熱)、Case2(內(nèi)壁面等溫)及Case3(內(nèi)壁面對(duì)流換熱)，Case2與Case3均基于機(jī)翼內(nèi)部滿油作為條件進(jìn)行設(shè)置，其中燃油物性參數(shù)參考文獻(xiàn)10?？紤]到高空的低溫效應(yīng)，假設(shè)壁面及燃油的初始溫度均為283.15 K，且將實(shí)際內(nèi)部為滿油工況作為Case4進(jìn)行對(duì)比。三種邊界條件及真實(shí)情況下的翼型表面溫度和熱流密度分布如圖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

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

3.3 三維機(jī)翼燃油箱滿油狀態(tài)下的氣動(dòng)熱計(jì)算

考慮內(nèi)壁面為對(duì)流換熱的基礎(chǔ)上，對(duì)三維機(jī)翼油箱在滿油狀態(tài)下進(jìn)行氣動(dòng)熱數(shù)值計(jì)算。機(jī)翼蒙皮厚度為5 mm，蒙皮材料為鋁，導(dǎo)熱系數(shù)202.4 W/(m·k)。基于弦長(zhǎng)的無(wú)量綱雷諾數(shù)為1.09×107，三維結(jié)構(gòu)網(wǎng)格徑向、周向和展向的網(wǎng)格數(shù)分別為51、406和101。機(jī)翼表面壓力分布云圖如圖16所示，機(jī)翼前緣的壓力最高，沿徑向逐漸降低。

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

為了進(jìn)一步分析三維機(jī)翼表面的氣動(dòng)熱環(huán)境，沿展向截取三個(gè)截面，如圖17所示：中心截面A0、距離中心截面左右兩側(cè)各1 m的截面A1和A2。

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

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

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

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

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

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

4 結(jié) 論

(1) 圓管與翼型的超音速繞流數(shù)值模擬驗(yàn)證結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好，本文采用的耦合方法可以對(duì)翼型繞流進(jìn)行有效的氣動(dòng)熱耦合計(jì)算。

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

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