李 皓

（中國直升機設(shè)計研究所，江西 景德鎮(zhèn) 333000）

0 引言

隨著計算機性能的飛速發(fā)展以及計算流體力學(xué)（Computational Fuild Dynamics，CFD）技術(shù)的不斷進(jìn)步，數(shù)值模擬方法獲得了蓬勃發(fā)展。鑒于其經(jīng)濟(jì)快捷的特性，在設(shè)計初期使用數(shù)值模擬來指導(dǎo)設(shè)計過程能大大降低設(shè)計成本，減少設(shè)計周期。在飛機防除冰系統(tǒng)的設(shè)計過程中，首先需要通過數(shù)值模擬來確定防護(hù)范圍、防冰所需熱流等，進(jìn)而才能開始整個防除冰系統(tǒng)的初步設(shè)計[1]。

由于飛機飛行速度快，對流換熱是熱量傳遞過程中的主要熱流項，對結(jié)冰過程、加熱除冰過程有至關(guān)重要的影響。能否獲得準(zhǔn)確的對流換熱系數(shù)會直接關(guān)系到防除冰系統(tǒng)設(shè)計的準(zhǔn)確性，因此進(jìn)行對流換熱系數(shù)計算準(zhǔn)確性的研究是非常必要的[2]。

該文的數(shù)值計算工作基于FLUENT 軟件。作為如今CFD 領(lǐng)域最全面的軟件之一，F(xiàn)LUENT 提供了多種求解方法和湍流模型，包括Spalart-Allmaras（S-A）湍流模型和Transition SST 湍流模型等，具有較高的求解精度。該文通過計算不同底層網(wǎng)格厚度、不同疏密網(wǎng)格、不同參考溫度、不同換熱溫差下的對流換熱系數(shù)，分析參數(shù)變化對計算結(jié)果的影響，確定了更精確的計算方法。最后將采用S-A 湍流模型和Transition SST 湍流模型的計算結(jié)果與國外文獻(xiàn)中的試驗數(shù)據(jù)進(jìn)行對比，驗證了該數(shù)值計算結(jié)果的準(zhǔn)確性，可為數(shù)值仿真工作人員提供一定參考。

1 對流換熱系數(shù)的計算

通過牛頓冷卻公式可以得到對流換熱系數(shù)h，如公式（1）所示。

式中：q為壁面總熱流；Tw為壁面溫度；Tref為參考溫度。

根據(jù)公式（1）可知，選取不同的參考溫度會獲得不同的對流換熱系數(shù)。為盡可能獲得準(zhǔn)確值，應(yīng)該選取更加貼近實際的計算方法。因此該文分別選取來流靜溫T∞、來流總溫Tt和絕熱壁溫Tad作為參考溫度，來計算對流換熱系數(shù)，如公式（2）～公式（4）所示。

式中：r為恢復(fù)系數(shù)，其值為1.4；Ma為馬赫數(shù)。

2 網(wǎng)格無關(guān)性驗證

采用FLUENT軟件計算對流換熱系數(shù)實質(zhì)上是計算空氣流場，并求解Navier-Stokes 方程?？諝饬鲌龅臄?shù)值計算受到計算域的空間離散條件的影響，在貼近壁面附近需要足夠細(xì)密的網(wǎng)格來解析黏性底層。為了獲得網(wǎng)格無關(guān)解，需要進(jìn)行網(wǎng)格無關(guān)性驗證。

首先，該文選擇半圓型的試驗件作為仿真模型，試驗件半圓直徑為44.5 mm，并在后面連接一個錐狀的后體以保持連續(xù)的氣動外形，易于計算收斂。其次，利用Gridgen 軟件來劃分結(jié)構(gòu)網(wǎng)格，如圖1 所示。計算域尺寸約為7 m，只有范圍足夠大，才能保證遠(yuǎn)場邊界條件成立，不會影響到壁面附近的空氣流場。并選取0.001 mm、0.005 mm、0.01 mm、0.05 mm 和0.1 mm 共5 種不同的底層網(wǎng)格厚度進(jìn)行計算，不同網(wǎng)格的無量綱壁面距離y+在0.1～45.0。

圖1 Gridgen 劃分C 字型結(jié)構(gòu)網(wǎng)格

再次，選用Spalart-Allmaras 湍流模型計算空氣流場?？諝獠捎美硐霘怏w，比熱cp=1 006.43 J/（kg·K），空氣導(dǎo)熱系數(shù)λ=0.023 56 W/（m·K），動力黏性系數(shù)μ=1.789 4×10-5kg/（m·s）。并設(shè)置壓力遠(yuǎn)場和等溫壁面，環(huán)境壓力P=34 500 Pa，馬赫數(shù)Ma=0.25，壁面溫度Tw=273.15 K，環(huán)境溫度Tinf=253.15 K。最后，通過壁面總熱流來分析網(wǎng)格無關(guān)性，通過計算得到的壁面總熱流如圖2 所示。

根據(jù)圖2 可知，基于不同底層網(wǎng)格厚度計算出來的壁面總熱流存在較大的差異，可見底層網(wǎng)格厚度能顯著影響計算結(jié)果。而隨著底層網(wǎng)格厚度減少，計算結(jié)果也逐漸趨于穩(wěn)定，厚度為0.001 mm 和0.005 mm 網(wǎng)格的計算結(jié)果幾乎完全重合。由于0.005 mm 厚度的網(wǎng)格的y+≈1，因此可以認(rèn)為0.005 mm 厚度的網(wǎng)格已經(jīng)足夠細(xì)密，計算結(jié)果即表現(xiàn)出網(wǎng)格無關(guān)性。

圖2 不同底層網(wǎng)格厚度下的壁面總熱流

除了底層網(wǎng)格厚度對計算結(jié)果有較大影響之外，網(wǎng)格的疏密也會影響計算的準(zhǔn)確性。因此為了獲得網(wǎng)格無關(guān)解，還需要驗證網(wǎng)格疏密對計算結(jié)果的影響。該文選取3套不同疏密的網(wǎng)格，底層網(wǎng)格厚度均為0.005 mm，網(wǎng)格數(shù)量分別為70 850 個、58 000 個和30 240 個。計算得到的壁面總熱流結(jié)果如圖3 所示。

根據(jù)圖3 可知，基于3 套不同疏密的網(wǎng)格計算得到的壁面總熱流在總體趨勢和范圍上差距不大，但細(xì)節(jié)上存在較大的不同。粗糙網(wǎng)格的網(wǎng)格數(shù)量少，劃分不夠精細(xì)，熱流曲線呈折線狀，過渡不夠平滑。精細(xì)網(wǎng)格的網(wǎng)格數(shù)量較多，反而出現(xiàn)了鋸齒形的波動，可能是由網(wǎng)格劃分的質(zhì)量引起的?？梢娭械葦?shù)量網(wǎng)格的計算結(jié)果的連續(xù)性最好。

圖3 不同疏密網(wǎng)格下的壁面總熱流

3 參考溫度對對流換熱系數(shù)的影響

該文選擇來流靜溫作為參考溫度來計算對流換熱系數(shù)h，計算條件模型和網(wǎng)格同第2 節(jié)，底層網(wǎng)格厚度0.005 mm，網(wǎng)格數(shù)量適中。采用Transition SST 湍流模型，空氣采用理想氣體，并設(shè)置壓力遠(yuǎn)場和等溫壁面，環(huán)境壓力P=34 500 Pa，馬赫數(shù)Ma=0.25，攻角為0°，環(huán)境溫度T∞=253.15 K。只改變壁面溫度，分別取Tw為263.15 K、273.15 K、283.15 K、293.15 K、303.15 K 和313.15 K。計算結(jié)果如圖4 所示。

根據(jù)圖4 可知，基于不同溫差（Tw-Tinf）計算得到的對流換熱系數(shù)相差極大，即隨著壁面溫度的升高，駐點附近的對流換熱系數(shù)逐漸增大，而遠(yuǎn)離駐點的下游處對流換熱系數(shù)逐漸減小。并且隨著溫差的增大，計算結(jié)果的相對誤差逐漸減小。產(chǎn)生該結(jié)果的原因主要在于未考慮氣動加熱的影響，由于駐點附近的空氣速度顯著降低，動能轉(zhuǎn)化為熱能，導(dǎo)致溫度有所升高，因此實際的換熱溫差小于計算值。

顯然，簡單地將來流靜溫作為參考溫度來計算對流換熱系數(shù)會存在較大誤差，為了消除氣動加熱的影響，考慮將來流總溫作為參考溫度，在相同條件下計算得到的對流換熱系數(shù)曲線如圖5 所示。

根據(jù)圖5 可知，將來流總溫作為參考溫度進(jìn)行計算后，消除了駐點附近的計算誤差，基于不同壁面溫度計算得到的駐點附近的對流換熱系數(shù)完全相同。因為在駐點附近空氣速度接近于0 m/s，可以認(rèn)為實際溫度就是來流總溫。而在遠(yuǎn)離駐點的下游位置，將來流總溫作為參考溫度進(jìn)行計算則不夠準(zhǔn)確。根據(jù)圖5 也可知，遠(yuǎn)離駐點位置在換熱溫差較小時存在較大相對誤差?？偨Y(jié)圖4 和圖5 的方法可知，采用來流靜溫作為壁面溫度會導(dǎo)致駐點附近的溫度與實際不符，而采用來流總溫作為壁面溫度則會導(dǎo)致遠(yuǎn)離駐點區(qū)域的溫度與實際不符。也即這2 種計算方法都存在一定不足，均會引入較大的相對誤差，因此采用等溫壁面進(jìn)行計算是不合適的。

圖4 不同換熱溫差下對流換熱系數(shù)計算結(jié)果（Transition SST 模型，壓力遠(yuǎn)場和等溫壁面邊界條件）

為了減少誤差，如果采用絕熱壁面計算出在空氣流場作用下壁面的實際溫度，并認(rèn)為在壁面附近的空氣溫度與壁面溫度相同，再將該溫度作為參考溫度來計算壁面的對流換熱系數(shù)，就可以消除氣動加熱引起的誤差，改善結(jié)果的準(zhǔn)確性。計算結(jié)果如圖6 所示。

根據(jù)圖6 可知，與圖4、圖5 的2 種方法相比，圖6 的計算結(jié)果有極大改善。駐點附近計算結(jié)果完全相同，遠(yuǎn)離駐點的下游區(qū)域的相對誤差也顯著減小。因此采用該種方法作為對流換熱系數(shù)的簡化算法，無須空氣流場的反復(fù)迭代，誤差又能控制在可接受范圍內(nèi)。其對流換熱系數(shù)h的計算如公式（5）所示。

圖5 不同換熱溫差下對流換熱系數(shù)計算結(jié)果（Transition SST 模型，壓力遠(yuǎn)場和等溫壁面邊界條件，考慮滯止溫升）

圖6 不同換熱溫差下對流換熱系數(shù)計算結(jié)果（Transition SST 模型，壓力遠(yuǎn)場和絕熱壁面邊界條件）

駐點無量綱距離弧長/弦長

式中：q為壁面總熱流；Tad為絕熱條件下的壁面溫度。

4 算例驗證

為了驗證該文對流換熱系數(shù)計算的準(zhǔn)確性，該文選取了美國國家航空航天局（NASA）在Lewis 冰風(fēng)洞中對NACA0012翼型所做的對流換熱試驗作為對比。該翼型弦長為0.533 m，取其中2 個算例作為驗證，見表1。計算同時采用了S-A 湍流模型和Transition SST 湍流模型，并采用遠(yuǎn)場邊界條件，調(diào)節(jié)馬赫數(shù)使自由流速度與試驗數(shù)據(jù)吻合，空氣采用理想氣體，密度為1.326 kg/m3，動力黏度為1.676 7×10-5kg/（m·s），比熱為1 006.43 J/（kg·K），空氣導(dǎo)熱系數(shù)為0.023 56 W/（m·K），環(huán)境溫度為266.45 K。為了減少相對誤差，盡可能增大換熱溫差，取壁面溫度比環(huán)境溫度高60 ℃。計算所用網(wǎng)格滿足第2 節(jié)網(wǎng)格無關(guān)性驗證中的條件y+＜1，第一層網(wǎng)格厚度為0.005 mm。

表1 對流換熱系數(shù)計算條件

圖7、圖8 分別為算例1、算例2 對流換熱系數(shù)計算結(jié)果與文獻(xiàn)試驗結(jié)果的比較。根據(jù)圖7、圖8 可知，不同的湍流模型計算結(jié)果存在較大差異。在駐點附近，S-A 湍流模型和Transition SST 湍流模型的計算結(jié)果基本相同，同時也與文獻(xiàn)試驗結(jié)果基本吻合，很好地描述了實際的換熱情況。但隨著位置往后移動，S-A 模型的轉(zhuǎn)捩很快出現(xiàn)，對流換熱系數(shù)不再下降，反而有微小回升，比文獻(xiàn)試驗結(jié)果偏高。而Transition SST 模型則在整個試驗測量范圍內(nèi)有良好的仿真效果。因此該文后面對對流換熱系數(shù)的計算全部采用Transition SST 湍流模型。

圖7 算例1 對流換熱系數(shù)計算結(jié)果

圖8 算例2 對流換熱系數(shù)計算結(jié)果

5 結(jié)論

該文采用FLUENT 軟件計算對流換熱系數(shù)，根據(jù)上述分析和比較可知，網(wǎng)格疏密、第一層網(wǎng)格的厚度、參考溫度、換熱溫差和湍流模型都會對對流換熱系數(shù)的計算結(jié)果產(chǎn)生影響。要想得到盡可能準(zhǔn)確的計算結(jié)果，應(yīng)遵循以下規(guī)律：1） 保證較好的網(wǎng)格質(zhì)量。過于密集的網(wǎng)格可能造成網(wǎng)格質(zhì)量的下降，過于稀疏的網(wǎng)格則會引入較大的計算誤差。2） 底層網(wǎng)格厚度對對流換熱系數(shù)計算準(zhǔn)確性影響很大，應(yīng)盡可能保證網(wǎng)格的無量綱壁面距離y+接近或小于1。3） 在仿真時盡可能增大換熱溫差，以縮小相對誤差。4） 選擇不同的參考溫度會影響對流換熱系數(shù)計算結(jié)果。應(yīng)采用絕熱壁面計算氣動加熱作用下的壁面溫度，用該溫度作為參考溫度來計算對流換熱系數(shù)，這樣能消除由于氣動加熱引起的誤差。5） 不同湍流模型對計算結(jié)果也有影響。在該文的算例下，Transition SST 模型的計算結(jié)果表現(xiàn)出了更好的符合性。