盧瑞博， 寧可為， 邱志方， 余霖， 趙富龍， 譚思超

(1.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150001； 2.黑龍江省核動力裝置性能與設(shè)備重點實驗室，黑龍江 哈爾濱 150001； 3.中國核動力研究設(shè)計院,四川 成都 610000)

臨近空間距離地面20～100 km，近年來因其在天氣預(yù)報，電磁通信，高分辨成像等方面具有重要應(yīng)用前景而引起廣泛關(guān)注[1]。常規(guī)動力飛機和衛(wèi)星無法在臨近空間穩(wěn)定運行，而采用核動力發(fā)動機的飛行器由于具有功率密度高、續(xù)航時間長、推力大、覆蓋范圍廣等優(yōu)點，可作為一種臨近空間飛行器。該飛行器以反應(yīng)堆作為熱源，高溫高速空氣為冷卻劑兼推進劑[2]。反應(yīng)堆處于高熱流密度，不均勻功率分布的條件，與可壓縮空氣間的流動換熱特性直接影響核動力發(fā)動機的安全和運行性能；現(xiàn)有實驗條件較難模擬堆內(nèi)換熱工況，因此選擇和構(gòu)建合理的空氣流動換熱數(shù)值計算模型，對預(yù)測復(fù)雜條件下空氣換熱特性具有重要意義。學(xué)者對圓管內(nèi)空氣流動換熱開展了試驗和理論研究，在圓管平均換熱系數(shù)方面，Humble 等[3]研究了內(nèi)徑10.21 mm，有效加熱長度609.6 mm的加熱管的流動換熱特性，分析了普朗特數(shù)對流動換熱關(guān)系式的影響，實驗中雷諾數(shù)為1萬～25萬。Lowdermilk等[4]基于Humble的研究，將雷諾數(shù)擴展到50萬，內(nèi)壁面溫度到1 139 K，得到的換熱關(guān)系式較前者偏大。Desmon 等[5]采用了耐高溫的鉑金試驗段，擴展內(nèi)壁面溫度到1 696 K，探究了定性溫度對換熱關(guān)系式的影響，研究表明，定性溫度的選擇可在壁面溫度修正和平均膜溫修正之間。Deissler[6]通過理論分析得到了普朗特數(shù)為0.73條件下的傳熱關(guān)系式，結(jié)果表明，合理選擇定性溫度，物性對換熱系數(shù)的影響可以忽略。Humble[7]開展的實驗中，內(nèi)壁面溫度達到1 694 K，擴展入口溫度到833 K，在高熱流密度和入口溫度條件下，由于空氣溫度和速度梯度大、物性差異明顯，傳統(tǒng)的經(jīng)驗關(guān)系式適用性較差，需要引入溫度修正系數(shù)。文獻[3-7]考慮了空氣的可壓縮性對換熱的影響，加熱管中的空氣密度沿流動方向變化較大。由于等截面加熱管道的影響，管內(nèi)流動均為亞音速，管道出口空氣流速接近音速。

針對空氣流動換熱開展的研究中，傳熱關(guān)系式大多符合D-B公式的形式，關(guān)系式系數(shù)有所差異；已有研究采用不同的定性溫度對空氣物性進行處理，但無法就處理方法達成一致；實驗中采用了測量試驗管道進出口儲氣罐中空氣滯止參數(shù)的方法測量溫度和壓力，與管道進出口空氣參數(shù)相比有一定差異，需要得到管內(nèi)空氣參數(shù)的精細(xì)化模型；同時，實驗中較難模擬反應(yīng)堆通道加熱條件，且高加熱功率工況下的試驗成本較大。相比于試驗和理論研究方法，在誤差允許的范圍內(nèi)，通過合理設(shè)置邊界條件，數(shù)值模擬方法能夠得到試驗管道全場數(shù)據(jù)，預(yù)測所關(guān)心工況的流動換熱特性，為進一步的細(xì)致化試驗提供指導(dǎo)，具有諸多優(yōu)勢。因此，在空氣流動換熱數(shù)值計算中，確定合適的湍流模型、優(yōu)化全場數(shù)據(jù)的分析處理方法、開展高還原度的仿真試驗、提高多工況數(shù)值計算計算效率，以達到較精確模擬空氣流動換熱過程的目的，對于開展復(fù)雜換熱設(shè)備的初步設(shè)計具有重要意義。

本文基于核動力發(fā)動機反應(yīng)堆冷卻劑通道，擬以單圓管內(nèi)變物性空氣為研究對象，以Fluent為計算平臺，確定合理的高溫高速可壓空氣流動換熱數(shù)值計算方法。

1 空氣流動換熱計算模型及方法

1.1 物理及物性模型

核動力發(fā)動機反應(yīng)堆中燃料形式為中空六棱柱(如圖1所示)，冷卻劑通道直徑在5.4～5.8 mm[8-9]。在本數(shù)值計算中，根據(jù)反應(yīng)堆通道物理尺寸，選擇了直徑為4、6和8 mm，長度為0.6 m的管道，加熱管長徑比均大于60，減少長徑比對換熱的影響。

圖1 核動力發(fā)動機Fig.1 Schematic of a nuclear engine

采用二維軸旋轉(zhuǎn)模型表征三維圓管模型開展數(shù)值計算，上述2種模型的計算結(jié)果表明：換熱系數(shù)的相對誤差不超過±0.5%[10]，該表征方式可行。二維軸旋轉(zhuǎn)模型網(wǎng)格繪制方案如圖2所示，網(wǎng)格形式為四邊形長0.6 m，通過改變邊界節(jié)點的數(shù)量，控制網(wǎng)格數(shù)，網(wǎng)格節(jié)點增長率小于1.1。

圖2 流動換熱模型局部網(wǎng)格示意Fig.2 Schematic of the local mesh in the model of flow and heat transfer

本文利用網(wǎng)格數(shù)為1.8萬、3.8萬、9.4萬和22.4萬的模型，開展無關(guān)性驗證。結(jié)果表明網(wǎng)格數(shù)在3.8萬以上，速度、溫度和對流換熱系數(shù)等參數(shù)的相對誤差在±5%以內(nèi)，滿足網(wǎng)格無關(guān)性的要求，并且最終本文選擇網(wǎng)格數(shù)為9.4萬。對比了網(wǎng)格長徑比對計算結(jié)果的影響，結(jié)果表明：在文中較大長徑比下，能夠滿足計算精度要求。此外，為了適應(yīng)壁面函數(shù)要求，對近壁面網(wǎng)格加密處理。

壓力和溫度影響空氣密度和物性參數(shù)，其中溫度對物性參數(shù)的影響更大。采用理想氣體模型，考慮可壓縮效應(yīng)的影響，通過氣體狀態(tài)方程得到其密度。通過調(diào)用REFPROP得到不同溫度和壓力下空氣的物性參數(shù)，數(shù)據(jù)結(jié)果表明，可將空氣物性參數(shù)表示為溫度的單值函數(shù)，并以分段多項式的形式加載到材料物性庫中。文中通過交互命令調(diào)用NIST物性庫，探究了物性對換熱特性的影響，結(jié)果表明，調(diào)用數(shù)據(jù)庫和擬合多項式得到的換熱系數(shù)偏差在±2%以內(nèi)，因此在數(shù)據(jù)后處理過程中，采用分段多形式計算空氣物性參數(shù)，能夠確保數(shù)據(jù)后處理過程中物性參數(shù)的正確性。

1.2 驗證工況范圍

高空中飛行器的運行速度為3 Ma，運行溫度為-56.5 ℃，此時氣流進入飛行器增壓段。絕熱滯止到速度為0時，空氣溫度可達606 K，即來流空氣總溫為606 K，實際進入增壓段后空氣具有一定的流速。若按照絕熱滯止到一定的馬赫數(shù)，計算得到加熱區(qū)域空氣的靜溫和流速如表1所示，堆芯前的空氣馬赫數(shù)在0.3以下時，其靜溫達到572 K以上。因此在計算過程中，確定進口溫度范圍在450～750 K，進口壓力在0.3～1.5 MPa，雷諾數(shù)范圍在104～4×105。加熱采用均勻熱流或均勻壁溫的方式，加熱功率范圍在2～5 kW，空氣出口溫度在550～1 200 K。

表1 不同馬赫數(shù)對應(yīng)的空氣靜溫和流速Table 1 The static temperature and velocity of air under different Mach number

1.3 計算方法

文中數(shù)值計算采用了質(zhì)量流量入口、壓力出口及無滑移壁面的邊界條件，壁面為均勻熱流密度加熱條件；采用二階迎風(fēng)格式，殘差收斂至10-5。為探究數(shù)值計算模型的適用性，綜合考慮適用于加熱空氣管流的湍流模型。Spalart-Allmaras模型適用于繞流，Reynolds Stress模型適用于復(fù)雜三維流場、強旋流，因此計算中不考慮上述2種模型；最終選用Standardk-ε(S-KE); RNGk-ε(RNG-KE); Realizablek-ε(R-KE)及Standardk-ω(S-KW); BSLk-ω(BSL-KW); SSTk-ω(SST-KW)等粘性模型對比分析[11]。

不同壁面函數(shù)模型對y+的要求有所區(qū)別，且文中雷諾數(shù)范圍較廣，壁面平均y+隨流速變化，文中繪制網(wǎng)格考慮了該因素的影響，最終輸出所有算例的結(jié)果表明，在增強壁面函數(shù)條件下，88.3%的工況平均y+在5以下，其余在個別高流速工況下y+值在15以下，可基本滿足壁面函數(shù)對不同工況y+的影響。文獻[7]流動換熱關(guān)系式的計算方法，基于能量守恒方法計算加熱管的平均換熱系數(shù)和通用傳熱關(guān)系式為：

(1)

Nu=hD/λa=f(Rea0.8Pra0.4)

(2)

式中:T為總溫；下標(biāo)o表示出口；下標(biāo)i表示入口；S為換熱面積；m為質(zhì)量流量；cp,b為根據(jù)流體平均溫度計算得到的定壓比熱容；cp為使用不同定性溫度得到的空氣的定壓比熱；λ為導(dǎo)熱系數(shù)；μ為粘度。

根據(jù)流動換熱通用傳熱關(guān)系式 (2)，采用不同的定性溫度如進出口平均溫度Tb和平均膜溫Tf，可得到不同的傳熱關(guān)系式，總結(jié)不同學(xué)者提出的空氣流動換熱關(guān)系式為：

Nu=c×Re0.8Pr0.4=c(md/Aμ)0.8(cpμ/λ)0.4

(3)

Tb=(Ti+To)/2，Tf=(Ts+Tb)/2

(4)

式中：c為努塞爾數(shù)表達式的系數(shù)，壁面與流體溫度修正及長徑比修正項可加入其中；m為空氣質(zhì)量流量；d為加熱管徑；A為加熱管流通面積。

在計算開展之前，對數(shù)值模擬算法和求解器適用性進行了對比。壓力基求解器已經(jīng)過優(yōu)化，能夠求解可壓流。在保證計算精度的前提下，計算可壓縮空氣問題，壓力基求解器擁有更小的收斂步數(shù)，和更好的收斂性，因此文中選用了該求解器。對比不同算法的求解精度和收斂步數(shù)結(jié)果表明，壓力基求解器中不同算法對加熱空氣管流問題求解精度誤差在±1.5% 以內(nèi)，因此可忽略算法對計算結(jié)果的影響。由于耦合算法能夠同時求解流體3大方程，在計算高速可壓流動時更有優(yōu)勢，因此文中較少網(wǎng)格數(shù)量的算例采用該算法收斂步數(shù)較少，求解過程所需要的時長較短。綜合考慮算法計算精度和效率，文中計算采用Coupled耦合求解算法。

2 模型適用性分析

2.1 湍流模型適用性分析

針對選用的不同直徑管道，參考實驗中空氣流量調(diào)節(jié)方法，在確定邊界條件合理性的基礎(chǔ)上，逐漸增大空氣質(zhì)量流量，以滿足不同雷諾數(shù)的要求；改變加熱管壁面熱流密度，以滿足不同加熱功率的要求；改變空氣進口溫度參數(shù)，以滿足加熱管不同入口溫度條件的要求。數(shù)值計算模型工況的選擇滿足單一變量準(zhǔn)則，不同湍流模型采用相同的工況。通過數(shù)值計算方法，得到了不同加熱功率、入口溫度和壓力條件下加熱管流的空氣參數(shù)，為獲取管道的平均傳熱關(guān)系式，依據(jù)式 (1)～(3) 的數(shù)據(jù)處理方法，得到了加熱管Nu、Pr和Re的關(guān)系趨勢如圖3所示。

圖3 不同模型傳熱關(guān)系式對比Fig.3 The comparation of heat transfer correlations with different models

為了驗證模型的適用性，文中經(jīng)計算得到不同模型的傳熱關(guān)系式，并將其與較為成熟的空氣換熱公式[7](系數(shù)c為0.023，適用范圍：Re在39萬以下，入口溫度在297～833 K，壁面溫度在1 694 K以下) 和精度較高的Gnielinski公式 (適用范圍為：Re為2 300～100萬，Pr為0.6～10萬，為準(zhǔn)確度較高的換熱關(guān)系式，90%的數(shù)據(jù)與關(guān)系式偏差在±20%以內(nèi)，大部分在±10%以內(nèi)) 進行對比，以確定湍流模型的適用性，對比不同模型和經(jīng)驗關(guān)系式中的Nu與S-KE中Nu的偏差如圖4所示。通過數(shù)據(jù)后處理可知，不同模型的相同工況下，管內(nèi)氣體的流動參數(shù)相近，進出口溫差的相對偏差多數(shù)在±0.50%，少數(shù)在1.7%以內(nèi)。若僅關(guān)心管內(nèi)空氣流動參數(shù)，上述模型均能較好描述流體流動狀態(tài)參數(shù)。

圖4 不同模型Nu與S-KE模型的偏差Fig.4 The deviation of Nu from S-KE model

由圖3可知，上述6種模型下的Nu、Pr和Re的關(guān)系均滿足傳統(tǒng)傳熱關(guān)系式的一般形式，即式 (3)。不同模型能夠定性描述空氣加熱管流換熱效果，但不同模型對流體和固體間換熱的描述存在差異，固體壁面溫度求解不確定。因此，基于傳統(tǒng)試驗關(guān)系式，選擇適合的模型準(zhǔn)確模擬壁面溫度，是解決數(shù)值模擬流固耦合問題的關(guān)鍵。將不同模型的Nu與S-KE模型Nu對比如圖4所示，Nu主要受到固體壁面與流體的傳熱溫差的影響。幾種模型對比結(jié)果可知，RNG-KE模型對壁面溫度的預(yù)測值偏低，而R-KE對壁面溫度的預(yù)測結(jié)果更高，在文中工況范圍內(nèi)，不同模型對壁面溫度的計算偏差在±15.0%左右。半數(shù)湍流模型中的Nu相對偏高，S-KW、SST-KW和R-KE模型得到的Nu和2個經(jīng)驗關(guān)系式偏差均相對較小，初步滿足計算精度要求。

綜上，數(shù)值計算得到的傳熱關(guān)系式與Gnielinski公式和文獻[7]中公式的對比表明，在本文工況范圍內(nèi)，采用S-KW、SST-KW和R-KE均能較好描述空氣流動換熱特性，其中由于R-KE模型偏差在兩經(jīng)驗關(guān)系式之間，同時該模型可以較好解決粘性底層低雷諾數(shù)區(qū)域的計算[12]，其適用性更強?；诒竟?jié)定性分析結(jié)果，下文將采用R-KE模型計算空氣流動換熱工況，壁面處理采用增強壁面函數(shù) (y+值滿足要求)，對數(shù)值計算模型展開進一步分析。

2.2 定性溫度選擇

水的物性參數(shù)隨溫度變化較小，在加熱管足夠長、溫差范圍小的情況下，充分發(fā)展段壁面和流體的溫差幾乎不隨長度變化，可取充分發(fā)展段水的物性參數(shù)計算管道平均傳熱系數(shù)?？諝馕镄詤?shù)隨溫度變化較大，且在文中高加熱功率工況下，加熱管道進出口物性參數(shù)變化較大，考慮了加熱管為普通結(jié)構(gòu)(Type-A)、200 mm非加熱入口(Type-B)和圓柱形非加熱入口(Type-C)條件下空氣的局部換熱系數(shù)如圖5所示。入口形狀不能消除入口段效應(yīng)對空氣換熱強化的影響，在充分發(fā)展段，由于空氣加熱后膨脹，密度變小，換熱能力降低，局部對流換熱系數(shù)隨流向距離增加而減小，與文獻[13]的結(jié)果相一致。若選取充分發(fā)展段局部換熱系數(shù)作為平均換熱系數(shù)，會引入誤差。因此，對于加熱空氣管流，需選用合適的定性溫度計算平均換熱系數(shù)，以保證數(shù)值計算的后處理精度。本文依據(jù)進出口平均溫度、平均膜溫和壁面溫度，此3種定性溫度描述數(shù)值計算中空氣流動換熱關(guān)系式，其中，為了保證數(shù)據(jù)覆蓋較多工況范圍，每組關(guān)系式中的數(shù)據(jù)由120組工況構(gòu)成。

圖5 不同入口段條件下局部換熱系數(shù)Fig.5 Local heat transfer coefficient under different inlet section conditions

將仿真關(guān)系式仿真Nua與文獻 [7]中的關(guān)系式和Gnielinski公式進行對比如圖6所示。其中，基于Gnielinski公式數(shù)據(jù)，用±10%的誤差帶表示誤差范圍。經(jīng)過數(shù)據(jù)處理可知，3種定性溫度中平均壁面溫度數(shù)值最高、進出口平均溫度最低，平均膜溫同時考慮了流體和壁面溫度，因此數(shù)值在兩者之間。由于空氣熱導(dǎo)率、黏度和比熱容，在文中定性溫度范圍內(nèi)均為單調(diào)遞增?？諝馕镄詤?shù)影響傳熱關(guān)系式，根據(jù)式 (2)、 (3) 可知：在相同換熱量的條件下，以壁面溫度作為定性溫度計算得到的Nu和Re數(shù)值均偏小。定性溫度對Nu與Re的綜合影響效果表現(xiàn)為基于壁面溫度計算的傳熱關(guān)系式偏高，基于進出口平均溫度的傳熱關(guān)系式偏低，如圖6 (d) 所示。在誤差方面，對此不同定性溫度下的傳熱關(guān)系式與經(jīng)典關(guān)系式，在使用平均壁溫作為定性溫度計算空氣物性時，數(shù)值模擬結(jié)果仿真Nus與文獻[7]和Gnielinski 公式的偏差大于±10%；平均膜溫下仿真Nuf偏差略小，部分?jǐn)?shù)據(jù)誤差超過±10%；而進出口平均溫度下仿真Nub與Gnielinski公式誤差在±10%以內(nèi)。因此，在空氣流動換熱數(shù)值模擬中，可優(yōu)先選用進出口平均溫度作為定性溫度。

圖6 不同定性溫度下的傳熱關(guān)系式Fig.6 Conceptual schematic of a nuclear engine

2.3 粘性加熱效應(yīng)影響

流體粘度很大或流動速度很高的情況下需要考慮粘性加熱的影響[14]，文中雷諾數(shù)范圍較大，空氣流速范圍廣，因此在相同的邊界條件下，計算了考慮粘性加熱和不考慮條件下加熱管流的Nu如圖7所示。在Re小于20萬時，不考慮粘性加熱效應(yīng)條件下的Nu偏大，但較多數(shù)據(jù)與Gnielinski公式的相對誤差仍在10%的范圍內(nèi)；Re大于20萬時，由于空氣流速較快粘性加熱效應(yīng)明顯，兩者相對誤差超過10%。對于高速空氣加熱管流，屬于熱壁面問題，粘性加熱作用類似于內(nèi)熱源，表現(xiàn)為降低傳熱效果[14]。在文中開式堆工況范圍內(nèi)，不同雷諾數(shù)條件下，空氣流速高，粘性加熱效應(yīng)明顯，均須考慮該效應(yīng)對換熱的影響。

圖7 粘性加熱效應(yīng)影響Fig.7 The effect of viscous heating

3 結(jié)論

1) 在不考慮固體域的前提下，本文選用的6種湍流模型均能準(zhǔn)確預(yù)測空氣流動換熱參數(shù)，考慮流固耦合換熱特性時，R-KE模型對固體域溫度的求解精度更高，適用性更強。

2) 數(shù)值計算中推薦采用進出口平均溫度作為定性溫度，相對于平均膜溫和壁溫精度更高。

3)Re大于20萬時，粘性加熱效應(yīng)會使高雷諾數(shù)工況的Nu減小超過10%，與水相比影響程度明顯，空氣流動換熱數(shù)值計算，應(yīng)考慮粘性加熱效應(yīng)。

4) 本文可為復(fù)雜工況條件下，高溫高速空氣流動換熱特性研究中模型和數(shù)據(jù)處理方法選擇提供依據(jù)。