基于雷諾時均模型的高溫空氣流動傳熱數值研究

焦守華，王金雨，曾 未，孫 暢，李文強，楊 婷，柴 翔

(1.上海交通大學 核科學與工程學院，上海 200240；2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610213；3.上海電力大學 自動化工程學院，上海 200090)

持續穩定的動力供應將極大提高飛行器的作戰能力和續航能力。使用化石燃料作為動力源，會對飛行器的飛行速度、高度和范圍產生一定的限制。核能作為一種新型能源，與傳統的化石燃料相比具有能量密度高、體積占比小、清潔無污染等一系列優點。將核動力應用到這些飛行器，將極大提升它們的作戰能力和續航能力[1]。20世紀，美蘇在核動力發動機方面的技術發展路線上都不約而同地指向了核渦輪噴氣發動機。核渦輪發動機反應堆其堆芯通常采用高溫空氣作為冷卻劑，高溫空氣的流動傳熱特性直接關系到渦輪發動機的運行性能。因此數值研究反應堆內高溫空氣流動換熱特性，可以為設計高溫空氣流動傳熱實驗和核渦輪發動機反應堆提供一些參考和指導。

在氣體流動傳熱方面，國內外已經開展了很多研究， Locke J M和Landrum D B[2]對再生冷卻通道中的超臨界氫氣的流動傳熱特性進行了相關研究，并評估了相關的對流傳熱關系式對氫氣的適用性。Ehsan M M等[3]對處在加熱或者冷卻條件下的超臨界二氧化碳的流動傳熱和壓降特性以及相關的對流換熱關系式進行了全面綜述。Yangyu Guo和Moran Wan[4]系統地研究了微尺度氣體流動和傳熱的第二定律，闡明了在微觀尺度熱和流體流動中熵產生的物理性質和數學表達式。盧芳等[5]采用FLUENT軟件對板式換熱器通道內的氣體流動特性進行了數值計算。上述學者多是對非高溫空氣的流動傳熱特性研究，涉及高溫空氣的研究比較少。

本文使用開源CFD軟件OpenFOAM,對反應堆冷卻劑通道內高溫空氣的流動傳熱進行數值模擬，探究了直徑為4.24 mm、長度為760 mm的圓形通道內均勻加熱和非均勻加熱情況下高溫空氣的流動傳熱特性。研究結果對于設計高溫空氣流動傳熱實驗具有一定指導意義。

1 數值方法簡介

1.1 控制方程

本文研究的所有流動和傳熱過程都受三大基本定律支配，即質量守恒定律、動量守恒定律和能量守恒定律。在三維空間連續介質內，牛頓黏性流體遵循的三大控制方程為[6]：

質量方程為：

(1)

動量方程為：

(2)

(3)

(4)

能量方程為：

(5)

式中：t——時間；

x,y,z——空間X、Y、Z方向的位置；

u,v,w——流體流速在空間X、Y、Z方向上的分量；

p——作用在流體微元上面的壓力；

Fi——流體微元所受體積力沿空間X、Y、Z方向上的分量；

τij——表面黏性應力沿空間X、Y、Z方向上的分量；

Cp——流體的比定壓熱容；

T——流體的溫度；

λ——流體的導熱率；

ST——黏性耗散項。

1.2 空氣物性參數

本文選用理想氣體模型來研究高溫空氣在圓形通道內的流動傳熱。在XNJ-140E核渦輪發動機反應堆中[7]，冷卻劑的溫度在570～1 250 K范圍內，壓強在0.1～0.5 MPa范圍內。在相同的溫度下，壓強為0.1 MPa和壓強為0.5 MPa對應的空氣動力黏度、比定壓熱容和導熱系數的相對誤差均在0.1%范圍內，因此壓強對空氣的上述物性影響可以忽略不計，只需考慮溫度對上述物性的影響。通過NIST數據庫[8]得到不同溫度下空氣物性參數，然后進行多項式擬合，得到的空氣物性隨溫度變化關系式如下所示：

(6)

λ=-4.676×10-15×T4+

2.257×10-11×T3-4.617×10-8×T2+

9.25×10-5×T+0.002 17

(7)

(8)

2 計算過程

2.1 幾何模型

本文以XNJ-140E核渦輪發動機反應堆[7]的單個冷卻劑通道作為計算對象。該對象是直徑D為4.24 mm，長度L為760 mm的圓柱管。為了在加熱段進口得到充分發展的流動，在入口處增加了一長為200 mm的輔助進口段；為了降低出口邊界條件對加熱段結果的影響，在加熱段出口增加了一長為40 mm的延伸出口段，因此整個幾何長度為1 000 mm。本文研究的對象幾何結構如圖1所示，x方向為流體的流動方向，y方向為徑向。為了減少計算量，將網格創建為二維軸對稱。

圖1 幾何模型Fig.1 Geometric model

2.2 邊界條件

計算模型的邊界條件設置如下：

(1)進口邊界：給定質量流量和溫度;

(2)出口邊界：給定壓力;

(3)非加熱壁面邊界：無滑移和絕熱;

(4)加熱壁面邊界：給定熱流密度。

根據邊界條件的不同，本文的計算工況如表1所示，表1中的熱流密度為管道的平均熱流密度。其中，工況1、工況2、工況3用來研究均勻加熱下熱流密度的影響，工況2、工況4、工況5用來研究均勻加熱下進口溫度的影響，工況2、工況6、工況7用來研究均勻加熱下進口雷諾數的影響。工況8-13是非均勻加熱工況，工況8、工況9、工況10是余弦功率分布，歸一化功率最大值為1.57，工況11、工況12、工況13的功率沿軸向先不變后線性降低，這和美國研究人員建成的Tory試驗堆的軸向功率分布是一致的，本文將這種功率分布稱之為Tory型功率分布。

表1 計算工況Table 1 Case for simulation

2.3 網格敏感性分析

為了在有限的計算資源下，更高效地得到合理精確的計算結果，本文首先使用三套網格進行預計算。預計算設置相同的邊界條件，僅在網格數量上有所不同。并取加熱段出口的速度及溫度分布情況和加熱壁面溫度分布情況作為網格敏感性分析的判據。具體的網格劃分方案如表2所示。本文計算工況的Y+平均值為0.773，符合湍流模型適用范圍。

表2 不同的網格劃分方案Table 2 Different meshing schemes

不同網格數目下，工況2計算結果如圖2至圖4所示。由圖可知，三種網格數量計算的溫度和速度分布相差很小，相對誤差在1%范圍內，因此選用網格3的網格數目作為最終計算的網格標準。

圖2 加熱段出口溫度分布圖Fig.2 Temperature distribution of outlet

圖3 加熱段出口速度分布圖Fig.3 Velocity distribution of outlet

圖4 加熱壁面溫度分布圖Fig.4 Temperature distribution in heated wall

2.4 湍流模型選取

在各湍流數值模型中，k-ε模型能較好地模擬充分發展的湍流流動，k-ω模型則能較好地模擬各壓力梯度下近壁面處的流動情況。Menter提出k-ωSST模型，不僅在近壁面處保留原始的k-ω模型，也對充分發展的湍流流動狀況運用k-ε模型[9]。k方程和ω方程分別為：

(9)

(10)

(11)

(12)

(13)

式中：t——時間；

ρ——流體密度；

Gk——湍動能k的產生項；

Gω——比耗散率ω的產生項；

U——流體速度；

μ——動力黏度；

x——三維坐標系；

μt——流體的比定壓熱容,μt=ρkT;

T——湍流時間尺度；

S——流體平均應力張量的模；

y——距壁面的距離；

ν——流體的運動黏度。

σk、σω、β和β*可通過混合函數求解得到；γ為混合系數；CDkw與流體的橫向擴散項相關。

針對圓管內均勻加熱的氣體流動傳熱， Taylor等[10]總結氫氣、氦氣、氮氣和空氣在均勻加熱條件下的實驗數據，提出的均勻加熱條件下的氣體換熱關系式，如公式(14)所示，關系式計算值與實驗值偏差在25%以內。本文選用k-ωSST湍流模型對工況2進行預計算，計算得到沿流動方向上不同位置處的努塞爾數，并將模擬結果與Taylor關系式計算結果進行對比，結果如圖5所示。

圖5 工況2對應的Nu隨x/D變化Fig.5 Nu vs x/D for case2

(14)

式中：Tb——沿流動方向某一位置處的流體主流溫度；

Ts——與Tb相對應的壁面溫度。

由圖5可知，模擬計算值與由Taylor關系式計算的值在x<20D吻合很好，誤差極低，在x>20D趨勢一樣，誤差略微增大，但誤差總體還是在10%以內。因此本文選用k-ωSST湍流模型進行均勻和非均勻加熱下的數值計算。

3 計算結果及分析

3.1 均勻加熱情況下結果分析

圖6 均勻加熱情況下模擬值與Taylor關系式值對比Fig.6 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under uniform heating

選取工況1、工況2和工況3分析熱流密度對對流換熱的影響，得到的結果如圖7所示。由圖7可知，沿流動方向，努塞爾數一直在減小，即對流換熱能力越來越弱。此外，隨著熱流密度的增加，努塞爾數在減小，但幅度并不大。

圖7 不同熱流密度下努塞爾數沿軸向分布Fig.7 Profiles of Nu along the axial direction at different heat flux

選取工況2、工況4和工況5分析進口溫度對對流換熱的影響，計算結果如圖8和圖9所示。由圖8可知，沿流動方向，努塞爾數一直在減小，對流換熱能力越來越弱。在進口雷諾數一定的情況下，隨著進口溫度增加，努塞爾數會略微降低，溫度升高40 K，平均努塞爾數僅降低2%。

圖8 不同進口溫度下努塞爾數沿軸向分布Fig.8 Profiles of Nu along the axial direction at different inlet temperature

選取工況2、工況6和工況7來探究進口雷諾數對對流換熱的影響。計算結果如圖9所示。由圖9可知，沿流動方向，努塞爾數一直在減小，進口速度越大，對應的努塞爾數越高，進口雷諾數每增加7 500，平均努塞爾數會增加30%左右。

圖9 不同進口雷諾數下努塞爾數沿軸向分布Fig.9 Profiles of Nu along the axial direction at different Re of inlet

3.2 非均勻加熱情況下結果分析

選取工況1、工況2、工況3和工況8、工況9、工況10來研究均勻加熱和余弦加熱對對流換熱的影響，計算結果如圖10和圖11所示。由圖10可知，在余弦功率加熱下，沿流動方向，由于進口效應，對流換熱系數會先下降，之后由于熱流密度增加，對流換熱系數會上升。由圖11可知，與均勻加熱一樣，沿流動方向，努塞爾數也是一直在下降。在x<60D，均勻加熱下的努塞爾數要低于余弦加熱的，之后則相反。冷卻劑在入口段后形成了充分發展的流動。均勻加熱下進口段處加熱功率較高，近壁面處冷卻劑的密度和流速受溫度的變化而劇烈變化；在余弦加熱下，進口段處的加熱功率較低，湍流邊界層隨著加熱功率的增加而逐漸演變。

圖10 對流換熱系數沿軸向分布Fig.10 Profiles of convective heat transfer coefficient along the axial direction

圖11 努塞爾數沿軸向分布Fig.11 Profiles of Nu along the axial direction

將所有余弦加熱工況的模擬值與Taylor關系式計算值相比較，結果如圖12所示。由圖12可知，模擬計算值與Taylor關系式計算值偏差在10

圖12 余弦加熱情況下模擬值與Taylor關系式值對比Fig.12 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under cos-type heating

本文以工況2和工況9為例分析了在靠近進出口處和中間部分的流體徑向溫度分布情況，結果如下所示。由圖13可知，在x=10D處，均勻加熱下，流體徑向溫度分布很不均勻，壁面溫度是主流溫度的1.52倍，數值模擬值與Taylor關系式計算值偏差0.53%；余弦加熱下，該處由于熱流密度較低，流體徑向溫度分布相對比較均勻，壁面溫度是主流溫度的1.125倍，數值模擬值與Taylor關系式計算值偏差7.54%。均勻加熱和余弦功率加熱兩種加熱模式在x=90D處的熱流密度都較高，分別是277 155 W/m2和433 133 W/m2。由圖14可知，由于較高的熱流密度，壁面溫度分別是主流溫度的1.33和1.55倍。努塞爾數的數值模擬值與Taylor關系式計算值偏差分別是4.57%和2.86%。由圖15可知，在x=170D處，均勻加熱下，流體徑向溫度分布很不均勻，壁面溫度是主流溫度的1.24倍，數值模擬值與Taylor關系式計算值偏差5.12%；余弦加熱下，熱流密度較低，僅為均勻加熱下的1/4，流體徑向溫度分布相對比較均勻，壁面溫度是主流溫度的1.07倍，數值模擬值與Taylor關系式計算值偏差27.27%。由于Taylor關系式是針對均勻加熱的實驗數據總結而出，根據主流溫度確定的物性校正出來的，實驗中由于熱流密度較高，流體徑向溫度分布不均勻情況比較明顯，因此該關系式較適用于高熱流密度的加熱區域。但是在余弦加熱情況下，在進出口處，加熱功率較低，流體徑向溫度分布相對比較均勻，因此數值模擬值與Taylor關系式計算值偏差比較大。

圖13 x=10 D處溫度分布圖Fig.13 Temperature distribution at x=10 D

圖14 x=90 D處溫度分布圖Fig.14 Temperature distribution at x=90 D

圖15 x=170 D處溫度分布圖Fig.15 Temperature distribution at x=170 D

將所有Tory型加熱工況的模擬值與Taylor關系式值相比較，結果如圖16和圖17所示。Tory型加熱情況下，功率沿軸向先不變后線性降低，這和美國研究人員建成的Tory試驗堆的軸向功率分布是一致的。

圖16 Tory型加熱情況下模擬值與Taylor關系式值對比Fig.16 Comparison of Nu obtained by simulation against Nu calculated by Taylor correlation under Tory-type heating

圖17 Tory型加熱情況下模擬值與Taylor關系式值偏差Fig.17 Error of Nu obtained by simulation against Nu calculated by Taylor correlation under Tory-type heating

由圖16和圖17可知，在189區域，功率線性降低，屬于非均勻加熱情況，模擬計算值與Taylor關系式計算值的偏差沿流動方向逐漸增大，在出口處達到了最大值，最大值為15%。

3 結論

本文基于開源CFD軟件OpenFOAM，對直徑為4.24 mm、長度為760 mm的圓形通道內均勻加熱和非均勻加熱情況下高溫空氣的流動傳熱進行了研究，可以得到如下結論。

(1)均勻加熱情況下，k-ωSST湍流模型數值模擬結果與Taylor關系式計算結果吻合較好，趨勢一致，k-ωSST湍流模型適用于均勻加熱下的高溫空氣流動傳熱模擬。

(2)均勻加熱情況下，沿流動方向，努塞爾數一直在降低，對流換熱能力越來越弱。進口雷諾數對努塞爾數影響比較大，進口雷諾數每增加7 500，努塞爾數會增加30%。熱流密度和進口溫度對其影響較小。

(3)余弦加熱情況下，k-ωSST湍流模型數值模擬結果與Taylor關系式計算結果在進出口偏差比較大。

(4)由于Taylor關系式是針對均勻加熱的實驗數據總結而出，根據主流溫度確定的物性校正出來的，實驗中由于熱流密度較高，流體徑向溫度分布不均勻情況比較明顯，因此該關系式較適用于高熱流密度的加熱區域。但是在余弦功率加熱情況下，在進出口處，加熱功率較低，流體徑向溫度分布相對比較均勻，因此數值模擬值與Taylor關系式計算值偏差比較大。

(5)Tory型加熱情況下，沿流動方向，k-ωSST湍流模型模擬結果與Taylor關系式計算結果在前半段偏差比較小，在后半段偏差逐漸增大，并在出口處達到了最大值。