基于雷諾時(shí)均模型的高溫空氣流動(dòng)傳熱數(shù)值研究

焦守華，王金雨，曾 未，孫 暢，李文強(qiáng)，楊 婷，柴 翔

(1.上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240；2.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213；3.上海電力大學(xué) 自動(dòng)化工程學(xué)院，上海 200090)

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

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

本文使用開源CFD軟件OpenFOAM,對(duì)反應(yīng)堆冷卻劑通道內(nèi)高溫空氣的流動(dòng)傳熱進(jìn)行數(shù)值模擬，探究了直徑為4.24 mm、長(zhǎng)度為760 mm的圓形通道內(nèi)均勻加熱和非均勻加熱情況下高溫空氣的流動(dòng)傳熱特性。研究結(jié)果對(duì)于設(shè)計(jì)高溫空氣流動(dòng)傳熱實(shí)驗(yàn)具有一定指導(dǎo)意義。

1 數(shù)值方法簡(jiǎn)介

1.1 控制方程

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

質(zhì)量方程為：

(1)

動(dòng)量方程為：

(2)

(3)

(4)

能量方程為：

(5)

式中：t——時(shí)間；

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

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

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

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

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

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

T——流體的溫度；

λ——流體的導(dǎo)熱率；

ST——黏性耗散項(xiàng)。

1.2 空氣物性參數(shù)

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

(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 計(jì)算過程

2.1 幾何模型

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

圖1 幾何模型Fig.1 Geometric model

2.2 邊界條件

計(jì)算模型的邊界條件設(shè)置如下：

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

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

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

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

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

表1 計(jì)算工況Table 1 Case for simulation

2.3 網(wǎng)格敏感性分析

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

表2 不同的網(wǎng)格劃分方案Table 2 Different meshing schemes

不同網(wǎng)格數(shù)目下，工況2計(jì)算結(jié)果如圖2至圖4所示。由圖可知，三種網(wǎng)格數(shù)量計(jì)算的溫度和速度分布相差很小，相對(duì)誤差在1%范圍內(nèi)，因此選用網(wǎng)格3的網(wǎng)格數(shù)目作為最終計(jì)算的網(wǎng)格標(biāo)準(zhǔn)。

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

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

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

2.4 湍流模型選取

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

(9)

(10)

(11)

(12)

(13)

式中：t——時(shí)間；

ρ——流體密度；

Gk——湍動(dòng)能k的產(chǎn)生項(xiàng)；

Gω——比耗散率ω的產(chǎn)生項(xiàng)；

U——流體速度；

μ——?jiǎng)恿︷ざ龋?/p>

x——三維坐標(biāo)系；

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

T——湍流時(shí)間尺度；

S——流體平均應(yīng)力張量的模；

y——距壁面的距離；

ν——流體的運(yùn)動(dòng)黏度。

σk、σω、β和β*可通過混合函數(shù)求解得到；γ為混合系數(shù)；CDkw與流體的橫向擴(kuò)散項(xiàng)相關(guān)。

針對(duì)圓管內(nèi)均勻加熱的氣體流動(dòng)傳熱， Taylor等[10]總結(jié)氫氣、氦氣、氮?dú)夂涂諝庠诰鶆蚣訜釛l件下的實(shí)驗(yàn)數(shù)據(jù)，提出的均勻加熱條件下的氣體換熱關(guān)系式，如公式(14)所示，關(guān)系式計(jì)算值與實(shí)驗(yàn)值偏差在25%以內(nèi)。本文選用k-ωSST湍流模型對(duì)工況2進(jìn)行預(yù)計(jì)算，計(jì)算得到沿流動(dòng)方向上不同位置處的努塞爾數(shù)，并將模擬結(jié)果與Taylor關(guān)系式計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果如圖5所示。

圖5 工況2對(duì)應(yīng)的Nu隨x/D變化Fig.5 Nu vs x/D for case2

(14)

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

Ts——與Tb相對(duì)應(yīng)的壁面溫度。

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

3 計(jì)算結(jié)果及分析

3.1 均勻加熱情況下結(jié)果分析

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

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

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

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

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

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

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

3.2 非均勻加熱情況下結(jié)果分析

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

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

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

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

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

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

圖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關(guān)系式值相比較，結(jié)果如圖16和圖17所示。Tory型加熱情況下，功率沿軸向先不變后線性降低，這和美國(guó)研究人員建成的Tory試驗(yàn)堆的軸向功率分布是一致的。

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

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

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

3 結(jié)論

本文基于開源CFD軟件OpenFOAM，對(duì)直徑為4.24 mm、長(zhǎng)度為760 mm的圓形通道內(nèi)均勻加熱和非均勻加熱情況下高溫空氣的流動(dòng)傳熱進(jìn)行了研究，可以得到如下結(jié)論。

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

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

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

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

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