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

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

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

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng) 621000）

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

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

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

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

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

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

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

1.1 CFD數(shù)值模擬

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

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

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

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

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

層流時(shí)[6]，有：

湍流時(shí)[6]，有：

式中：為普朗特?cái)?shù)。

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

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

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

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

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

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

D-B方程[8]為：

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

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

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

Gnielinski方程為：

其中：

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

2 數(shù)值模擬

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

圖2 平板溫度分布

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

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

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

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

3 高雷諾數(shù)算例

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

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

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

4 超音速等直徑管換熱

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

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

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

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

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

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

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

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

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

5 收縮噴管段實(shí)際模型換熱分析

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

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

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

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

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

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

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

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

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

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

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

6 總結(jié)

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

此外，本文對(duì)低溫風(fēng)洞收縮噴管段實(shí)際模型進(jìn)行了CFD仿真，與經(jīng)驗(yàn)公式的算例結(jié)果進(jìn)行了對(duì)比分析，驗(yàn)證了采用恢復(fù)溫度修正傳熱換熱系數(shù)定義公式進(jìn)行CFD數(shù)值模擬計(jì)算換熱系數(shù)的科學(xué)性，以及兩個(gè)公式在高雷諾數(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]賴(lài)歡，陳振華，高榮，等. 大型高速低溫風(fēng)洞冷量回收的方法研究[J]. 西安交通大學(xué)學(xué)報(bào)，2016，50（6）：136-142.

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

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

[5]嚴(yán)慧芳，曾軍，婁德倉(cāng)，等. 高速?lài)姽軆?nèi)部對(duì)流換熱計(jì)算方法研究[J]. 航空科學(xué)技術(shù)，2015（10）：68-72.

[6]毛旭. 對(duì)流換熱邊界下高超音速平板邊界層相似性解研究[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]許鵬博，孫曉玲，章錦威，等. 高超聲速氣動(dòng)加熱/結(jié)構(gòu)傳熱耦合計(jì)算研究[C]. 全國(guó)計(jì)算流體力學(xué)會(huì)議，2014.

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

[12]杜廣生. 工程流體力學(xué)[M]. 北京：中國(guó)電力出版社，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]文華，王玲，殷國(guó)富，等. 插裝閥集成塊內(nèi)部流道流場(chǎng)分析研究[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－），男，四川仁壽人，碩士，工程師，主要研究方向?yàn)轱L(fēng)洞設(shè)備結(jié)構(gòu)設(shè)計(jì)等，E-mail：gaokerry@126.com；聶旭濤（1979－），男，安徽宿松人，博士，副研究員，主要研究方向?yàn)轱L(fēng)洞設(shè)備結(jié)構(gòu)設(shè)計(jì)等。