高雷諾數(shù)低溫風(fēng)洞內(nèi)傳熱特性數(shù)值研究

高鑫宇，聶旭濤，蔡清青，張偉

高雷諾數(shù)低溫風(fēng)洞內(nèi)傳熱特性數(shù)值研究

高鑫宇，聶旭濤，蔡清青，張偉

（中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

低溫風(fēng)洞的超音速及其高雷諾數(shù)工況下的換熱特性是分析低溫風(fēng)洞組件溫度場分布的關(guān)鍵因素。對超音速及高雷諾數(shù)工況換熱特性進(jìn)行數(shù)值模擬，通過采用恢復(fù)溫度修正傳熱系數(shù)的方法，并且與Dittus-Boleter（D-B）公式和Gnielinski公式的計算結(jié)果進(jìn)行對比分析。結(jié)果顯示，采用恢復(fù)溫度修正后的D-B公式和Gnielinski公式在超音速及其高雷諾數(shù)工況下仍然具有較高準(zhǔn)確性。最后將噴管段實際模型工況數(shù)值模擬結(jié)果與理論預(yù)測值進(jìn)行對比分析，進(jìn)一步驗證了在溫度分布不均勻時采用D-B公式和Gnielinski公式求得的平均換熱系數(shù)在工程計算中的可行性，以及使用恢復(fù)溫度修正傳熱系數(shù)計算公式的有效性。

低溫風(fēng)洞；高雷諾數(shù)；超音速；換熱特性；恢復(fù)溫度

隨著空氣動力學(xué)的飛速發(fā)展，對作為試驗設(shè)備的風(fēng)洞有了更高的要求。目前常規(guī)風(fēng)洞已經(jīng)無法完成雷諾數(shù)高于50×106的試驗要求，而低溫風(fēng)洞[1]通過降低氣流溫度的方法達(dá)到提高風(fēng)洞雷諾數(shù)模擬的目的，已經(jīng)成為目前實現(xiàn)高雷諾數(shù)模擬的最佳方案。由于液氮在低溫工質(zhì)中廉價和安全，汽化后可迅速吸收大量熱量從而噴射液氮，成為目前低溫風(fēng)洞中實現(xiàn)風(fēng)洞空氣回路快速降溫的主要途徑。目前，世界上兩座大型的跨聲速低溫風(fēng)洞——美國國家跨聲速風(fēng)洞（NTF）和歐洲跨聲速風(fēng)洞（ETW）均采用液氮工質(zhì)制冷實現(xiàn)低溫運(yùn)行[2-3]。低溫風(fēng)洞在運(yùn)行過程中，低溫流體在風(fēng)洞中快速流動，通過對流換熱與風(fēng)洞中的結(jié)構(gòu)建立熱平衡，從而實現(xiàn)低溫實驗條件，氣流溫度甚至最低可達(dá)到90 K[4]，在這樣低的氣流溫度下，整個風(fēng)洞的雷諾數(shù)很高，可達(dá)到108量級，使得氣流與風(fēng)洞的對流換熱劇烈，由于熱脹冷縮的作用，對風(fēng)洞各部件的安全性和密封性提出了更高的要求。為保證風(fēng)洞的安全運(yùn)行，需要對風(fēng)洞傳熱進(jìn)行全面的計算分析，尤其是風(fēng)洞處于高雷諾數(shù)和超音速時的換熱特性。

國內(nèi)外對高速流動換熱特性進(jìn)行了研究。文獻(xiàn)[5]在對高速噴管內(nèi)部的對流換熱計算方法進(jìn)行研究時，采用絕熱壁面溫度（恢復(fù)溫度）求解對流換熱系數(shù)。文獻(xiàn)[6]在進(jìn)行高速平板邊界層的對流換熱研究時，引入了恢復(fù)溫度求解對流換熱系數(shù)。文獻(xiàn)[7]在研究收斂擴(kuò)張噴管湍流邊界層的熱量測量方法時，引入了恢復(fù)溫度的概念。

本文采用Fluent軟件對高雷諾數(shù)及其超音速工況下的管內(nèi)流動對流換熱進(jìn)行研究，引入恢復(fù)溫度的概念進(jìn)行對流換熱求解，并與管內(nèi)對流換熱公式進(jìn)行了對比驗證。此外，對風(fēng)洞收縮噴管段模型換熱系數(shù)進(jìn)行模擬求解，與管內(nèi)流動強(qiáng)制對流換熱經(jīng)驗公式進(jìn)行對比分析。驗證了管內(nèi)流動在引入恢復(fù)溫度修正后的可行性和準(zhǔn)確性，并且驗證了D-B公式在雷諾數(shù)為107的量級也是適用的。

1 對流換熱系數(shù)的求解方法

對于管內(nèi)流動對流換熱系數(shù)主要通過CFD數(shù)值模擬和經(jīng)驗公式兩個方面進(jìn)行求解。

1.1 CFD數(shù)值模擬

換熱系數(shù)的求解公式[8]如下：

式中：q為對流換熱熱流密度，W/m2；T、T分別為物體表面和流體的溫度，K；為對流換熱系數(shù)，W/(m2·K)。

式（1）適用于低速流體，當(dāng)流體速度很高時則要考慮粘性流體摩擦引起的溫升，即[9-10]：

式中：T為壁面恢復(fù)溫度，K；為恢復(fù)系數(shù)，其物理意義為絕熱壁面由摩擦而引起的溫升與高速氣體駐點(diǎn)的絕熱溫升之比；為流體比熱比，無量綱；∞為流體主流溫度，K；∞為流體主流流動馬赫數(shù)，Ma。

當(dāng)∞為無窮小時T＝∞，高速流體的公式與低速流體的公式完全吻合。

層流時[6]，有：

湍流時[6]，有：

式中：為普朗特數(shù)。

用CFD數(shù)值模擬求解換熱系數(shù)，在給定恒定固體壁面溫度時，可通過CFD求解出流體與壁面的熱流密度，并通過式（3）求解出T，進(jìn)而通過式（2）求解出換熱系數(shù)。

本文使用Fluent軟件進(jìn)行CFD數(shù)值模擬，選用湍流模型、標(biāo)準(zhǔn)壁面函數(shù)，求解方式采用壓力基SIMPLE進(jìn)行求解，流體物性采用NIST物性數(shù)據(jù)庫。

1.2 采用經(jīng)驗公式求解換熱系數(shù)

根據(jù)相關(guān)研究[3]影響對流換熱系數(shù)的主要因素可定性地用函數(shù)形式表示為：

式中：為流體流動速度；、、、、c為幾何結(jié)構(gòu)換熱表面的特征長度、密度、粘度、導(dǎo)熱系數(shù)和定壓比熱。

在低溫風(fēng)洞中流體高速流動，主要涉及的換熱類型是管內(nèi)流動強(qiáng)制對流換熱，而常用管內(nèi)流動強(qiáng)制對流換熱經(jīng)驗公式有Dittus-Boelter（D-B）公式和Gnielinski方程。

D-B方程[8]為：

式中：當(dāng)流體被加熱時取0.4，當(dāng)流體被冷卻時取0.3；＞10000。

該經(jīng)驗公式一般要求≥10（其中：為流體管道直徑，m；為流體流動管長，m）。

當(dāng)/＜10時，需要考率入口效應(yīng)，工程上常采用的入口效應(yīng)修正系數(shù)為[9]：

Gnielinski方程為：

其中：

從圖1可以看出D-B公式和Gnielinski公式在104～108范圍內(nèi)整體吻合得很好，在106.5～108范圍，D-B公式的值較大，但目前兩個公式的實驗驗證范圍還沒有達(dá)到107～108量級。但當(dāng)風(fēng)洞處于高壓及高速流動甚至風(fēng)速在實驗段達(dá)到超音速時，管內(nèi)流動的雷諾數(shù)可以達(dá)到107～108這個量級。本文通過Fluent進(jìn)行數(shù)值模擬分析并通過與常用管內(nèi)流動經(jīng)驗公式進(jìn)行對比分析，給出一個可以預(yù)測風(fēng)洞處于高雷諾數(shù)時強(qiáng)制對流換熱可應(yīng)用的經(jīng)驗公式。

2 數(shù)值模擬

為了驗證CFD數(shù)值模擬的正確性，本文通過外掠平板換熱系數(shù)結(jié)果進(jìn)行驗證。本文在進(jìn)行數(shù)值模擬過程中取板的特征長度為1 m，氮?dú)鉁囟葹?00 K，壓強(qiáng)為1 bar，平板溫度取恒定壁溫220 K，材料為不銹鋼。邊界條件為：進(jìn)口采用速度進(jìn)口，出口采用流量出口，平板采用恒定壁溫，無滑移，其余壁面采用絕熱，無滑移。圖2為流體外掠平板的溫度分布，可以看出平板表面的熱邊界層厚度逐漸增厚，這符合熱力學(xué)一般現(xiàn)象。

圖2 平板溫度分布

這里定性溫度按210 K取，物性參數(shù)采用NIST REFPROP數(shù)據(jù)庫。

從圖3可看出，F(xiàn)luent仿真結(jié)果和解析解吻合得很好，并且誤差較小，因此可以采用CFD方法在給定恒定壁面溫度邊界條件，通過求解固體壁面和流體間換熱的熱流密度，從而對低溫流體求解換熱系數(shù)。

圖3 平板處于層流狀態(tài)時Nu模擬值和理論值對比

對于雷諾數(shù)＜500000，平板處于層流狀態(tài)，理論值根據(jù)如下經(jīng)驗公式計算：

3 高雷諾數(shù)算例

選取管道直徑為1 m，管長為20 m，取馬赫數(shù)為0.1、0.2、0.4、0.6的四種工況進(jìn)行驗證，流體溫度100 K，管壁溫度150 K，流體壓強(qiáng)為1 bar。邊界條件為進(jìn)口采用壓力進(jìn)口，出口采用壓力出口，壁面采用無滑移恒定壁面溫度。

經(jīng)過計算，管內(nèi)流動模擬值與理論值對比如圖4所示，其中CFD修正模擬值采用式（2）（3）即引入恢復(fù)溫度進(jìn)行修正，此時雷諾數(shù)范圍在107量級?？梢钥闯觯M值與經(jīng)驗公式計算值整體吻合得較好，尤其在低馬赫數(shù)時CFD未修正值和修正值差別不大、且均和D-B公式和Gnielinski公式吻合很好，這是由于低馬赫數(shù)時恢復(fù)溫度近似等于流體靜溫。而在高馬赫數(shù)，可以看出采用恢復(fù)溫度修正過后的CFD換熱系數(shù)模擬值與Gnielinski吻合很好，而未修正過的CFD模擬值與經(jīng)驗公式計算值差值較大，說明D-B公式和Gnielinski公式可以應(yīng)用于低溫風(fēng)洞中的高雷諾數(shù)范圍進(jìn)行工程計算，并且在采用CFD數(shù)值模擬進(jìn)行求解換熱系數(shù)時在高馬赫數(shù)時需要用恢復(fù)溫度對換熱系數(shù)定義式進(jìn)行修正。

圖4 換熱系數(shù)CFD模擬值與經(jīng)驗公式值對比

4 超音速等直徑管換熱

對超音速流速下的等直徑圓形管道進(jìn)行CFD數(shù)值模擬，用Fluent求解超音速換熱時采用密度基瞬態(tài)求解，直到算到穩(wěn)態(tài)，考慮氣體的可壓縮性，流體物性采用Refprop軟件擬合。考慮長徑比為10的等直徑管及長徑比為2的等直徑管，如表1所示。邊界條件為進(jìn)口采用壓力進(jìn)口，出口采用壓力出口，壁面采用無滑移恒定壁面溫度。

表1 等直徑管參數(shù)及其流體工況參數(shù)

4.1 算例1結(jié)果分析

由圖5可以看出，在0.3 m附近，即管長為直徑的3倍長度時，馬赫數(shù)有一個從超音速向亞音速的突變，出現(xiàn)管道堵塞現(xiàn)象[12]。即實際氣體在等直徑管中流動，等直徑管道相當(dāng)于漸縮管，考慮摩擦作用時由超音速氣流不能連續(xù)轉(zhuǎn)變?yōu)閬喴羲贇饬?，出口的極限狀態(tài)只能為聲速。

圖5 算例1的馬赫數(shù)分布

4.2 算例2結(jié)果分析

由圖6可以看出，管內(nèi)超聲速的流動現(xiàn)象基本上是符合理論基礎(chǔ)的，即等直徑管內(nèi)超聲速流動速度減小，但最小減小到跨音速，否則會出現(xiàn)堵塞現(xiàn)象[12]。整個管內(nèi)的流動均在超音速范圍內(nèi)，說明管長小于限制管長。

圖6 算例2的馬赫數(shù)分布

算例2中，?。? bar、＝100 K、＝1.21×108，通過式（2）（3）采用恢復(fù)溫度對處于超音速工況流體求解換熱系數(shù)的定義式進(jìn)行修正求得(CFD)＝986 W/(m2·K)，通過式（4）和（8）求得(D-B)＝928 W/(m2·K)，通過式（8）～（10）求得(Gnielinski)＝940 W/(m2·K)?？梢钥闯鯠-B經(jīng)驗公式和Gnielinski公式與CFD修正后的仿真結(jié)果均吻合較好，且兩個公式計算結(jié)果在雷諾數(shù)在108這個量級的計算結(jié)果相差不大，可以說明D-B公式和Gnielinski公式在超音速高雷諾數(shù)的適用性。

5 收縮噴管段實際模型換熱分析

為進(jìn)一步驗證D-B公式和Gnielinski公式在超音速高雷諾數(shù)下的適用性，對低溫風(fēng)洞收縮噴管段實際模型針對靜溫變化比較大的工況進(jìn)行CFD模擬。如圖7所示的風(fēng)洞收縮噴管段結(jié)構(gòu)模型，分為收縮段、硬板噴管段、柔板噴管段三個部分。

圖7 風(fēng)洞收縮噴管段模型

氮?dú)鈱崿F(xiàn)跨音速流動，需要考慮工質(zhì)壓縮性，使用壓力邊界條件，進(jìn)口壓力為inlet,total＝3.5 bar，出口壓力為outlet,static＝1.263199 bar，出口馬赫數(shù)達(dá)到1.3。壁面邊界條件對于靜態(tài)流場分析，壁面采用恒壁溫邊界條件，壁面是氣動降溫過程，流體壁面溫度設(shè)置為比流體入口主流溫度高50 K。

計算結(jié)果如圖8所示，壓力沿流動方向逐漸減小，速度逐漸增大，在喉部位置氮?dú)饬魉龠_(dá)到音速，v＝194.05 m/s。噴管內(nèi)氮?dú)鈮毫Σ粩嘟档停魉俨粩嘣黾?，?nèi)腔體壁面換熱性能沿氣流方向不斷變化。

由圖9所示的換熱系數(shù)分布可以看出，縮放噴管換熱系數(shù)沿流動方向先增大后減小，頂部柔板的峰值在噴管喉部之前＝312 mm時最大值為1351.56 W/(m2·K)，側(cè)壁面的峰值也分布在喉部之前的位置，＝410 mm時最大值為1115.33 W/(m2·K)，有文獻(xiàn)[4]中研究發(fā)現(xiàn)CFD模擬值在喉部以前略小于實驗數(shù)據(jù)，在喉部以后與實驗數(shù)據(jù)吻合較好，且熱流密度的峰值在喉部以前。比較圖9（a）（b），可以看出本文模擬的噴管段換熱系數(shù)變化規(guī)律與文獻(xiàn)中的趨勢完全一致，且峰值的位置也一致，說明了模擬的正確性。

圖8 噴管流場仿真結(jié)果云圖

圖9 沿噴管型線換熱系數(shù)分布

使用經(jīng)驗公式對噴管柔板處換熱系數(shù)進(jìn)行工程計算，柔板段模型為縮放噴管，取平均當(dāng)量直徑為d＝4A/＝0.2928 m?？s放噴管氣流速度始終增大，取當(dāng)?shù)匾羲?i>v＝194.05 m/s作為平均速度，根據(jù)壓力分布云圖取柔板處壓力為1.988 bar，溫度取靜溫static＝94.184 K，計算物性根據(jù)工況條件通過NIST進(jìn)行查找，計算結(jié)果如表2所示。

表2 流體物性參數(shù)

在＝194.05 m/s條件下求取噴管柔板處流道的面平均換熱系數(shù)，將數(shù)值分析值與經(jīng)驗公式計算值進(jìn)行對比，分析兩者誤差，驗證數(shù)值分析的準(zhǔn)確性。對比結(jié)果如表3所示，考慮到速度、溫度、壓力均為變化數(shù)值，取平均值會引起誤差，則計算結(jié)果可靠。

表3 噴管柔板換熱系數(shù)對比結(jié)果

6 總結(jié)

通過對低溫流體超音速及其高雷諾數(shù)的換熱特性通過CFD和經(jīng)驗公式對比分析，得出了在超音速及其高雷諾數(shù)工況，使用恢復(fù)溫度修正后的Dittus-Boelter公式和Gnielinski公式在高雷諾數(shù)時仍然具有實用性，并且兩個公式在雷諾數(shù)在108這個量級誤差很小，在106～108量級，Dittus-Boelter公式的計算結(jié)果偏大。

此外，本文對低溫風(fēng)洞收縮噴管段實際模型進(jìn)行了CFD仿真，與經(jīng)驗公式的算例結(jié)果進(jìn)行了對比分析，驗證了采用恢復(fù)溫度修正傳熱換熱系數(shù)定義公式進(jìn)行CFD數(shù)值模擬計算換熱系數(shù)的科學(xué)性，以及兩個公式在高雷諾數(shù)范圍的工程適用性。

[1]Goodyer M J，Kilgore R A. High-Reynolds-Number Cryogenic Wind Tunnel[J]. Journal of the Japan Society for Aeronautical & Space Sciences，2013，42（480）：26-31.

[2]賴歡，陳振華，高榮，等. 大型高速低溫風(fēng)洞冷量回收的方法研究[J]. 西安交通大學(xué)學(xué)報，2016，50（6）：136-142.

[3]劉政崇. 風(fēng)洞結(jié)構(gòu)設(shè)計（精）[M]. 北京：宇航出版社，2005.

[4]麻越垠，陳萬華，王元興，等. 基于ABAQUS的低溫風(fēng)洞擴(kuò)散段結(jié)構(gòu)-熱耦合仿真研究[C].中國計算力學(xué)大會2014暨錢令希計算力學(xué)獎頒獎大會，2014.

[5]嚴(yán)慧芳，曾軍，婁德倉，等. 高速噴管內(nèi)部對流換熱計算方法研究[J]. 航空科學(xué)技術(shù)，2015（10）：68-72.

[6]毛旭. 對流換熱邊界下高超音速平板邊界層相似性解研究[D]. 天津：天津大學(xué)，2013.

[7]Back L H，Massier P F，Gier H L. Convective heat transfer in a convergent-divergent nozzle[J]. International Journal of Heat and Mass Transfer，1964，7（5）：549-568.

[9]楊世銘，陶文銓. 傳熱學(xué)（BZ）（TB）[M]. 4版. 北京：高等教育出版社，2010.

[10]許鵬博，孫曉玲，章錦威，等. 高超聲速氣動加熱/結(jié)構(gòu)傳熱耦合計算研究[C]. 全國計算流體力學(xué)會議，2014.

[11]李鵬飛，吳頌平. 類航天飛機(jī)前身結(jié)構(gòu)與高超聲速流場的耦合傳熱模擬分析[J]. 航空動力學(xué)報，2010，25（8）：1705-1710.

[12]杜廣生. 工程流體力學(xué)[M]. 北京：中國電力出版社，2005.

[13]Stoll J，Straub J. Film cooling and heat transfer in nozzles[J]. Journal of turbo machinery，1988，110（1）：57-65.

[14]吳加俊，戴恒震. 基于ABAQUS的導(dǎo)電滑環(huán)的單絲靜態(tài)接觸力學(xué)特性分析[J]. 機(jī)械，2018，7（18）：18-22.

[15]文華，王玲，殷國富，等. 插裝閥集成塊內(nèi)部流道流場分析研究[J]. 機(jī)械，2017，5（10）：10-13.

Numerical Study of Heat Transfer Characteristics in Cryogenic Wind Tunnel at High Reynolds Number

GAO Xinyu，NIE Xutao，CAI Qingqing，ZHANG Wei

(China Aerodynamics Research and Development Center, Mianyang 621000, China )

The heat transfer characteristics of cryogenic wind tunnel under supersonic and high Reynolds number conditions are the key factors to analyze the distribution of the temperature field in cryogenicwind tunnel. The heat transfer characteristics under supersonic and high Reynolds number conditions are simulated by using the recovery temperature correction heat transfer coefficient calculation method, and the results are compared with the calculations by Dittus-Boelter equation and Gnielinski equation. The results show that the application of D-B formula and Gnielinski equation under supersonic and high Reynolds number conditions after the recovery temperature correction still have a high validity. Finally, the numerical simulation results of the spray pipe section are compared with the predicated values and analyzed, which further verifies the applicability of the average heat transfer coefficient obtained by the above two equations and the validity of the recovery temperature correction heat transfer coefficient calculation formula with non-uniform temperature distribution.

cryogenic wind tunnel；high Reynolds number；supersonic；heat transfer characteristics；recovery temperature

TB69

A

10.3969/j.issn.1006-0316.2021.11.005

1006－0316 (2021) 11－0034－07

2021－01－22

高鑫宇（1985－），男，四川仁壽人，碩士，工程師，主要研究方向為風(fēng)洞設(shè)備結(jié)構(gòu)設(shè)計等，E-mail：gaokerry@126.com；聶旭濤（1979－），男，安徽宿松人，博士，副研究員，主要研究方向為風(fēng)洞設(shè)備結(jié)構(gòu)設(shè)計等。