空氣流動換熱數值計算模型適用性分析

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

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

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

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

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

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

1.1 物理及物性模型

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

圖1 核動力發動機Fig.1 Schematic of a nuclear engine

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

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

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

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

1.2 驗證工況范圍

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

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

1.3 計算方法

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

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

2 模型適用性分析

2.1 湍流模型適用性分析

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

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

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

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

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

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

2.2 定性溫度選擇

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

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

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

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

2.3 粘性加熱效應影響

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

圖7 粘性加熱效應影響Fig.7 The effect of viscous heating

3 結論

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

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

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

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