李超， 王宗一， 廖敏， 趙欽新， 邵懷爽

(西安交通大學(xué) 熱流科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

超臨界二氧化碳sCO2布雷頓循環(huán)在核電領(lǐng)域具有諸多優(yōu)勢，循環(huán)效率能夠比蒸汽工質(zhì)高12%左右，安全性好，并且設(shè)備更加緊湊，具有良好的應(yīng)用前景[1]。然而超臨界二氧化碳物性隨溫度劇烈變化，其流動(dòng)換熱特性與常物性流體迥然不同，充分認(rèn)識(shí)并掌握sCO2流動(dòng)傳熱特性對(duì)于其在第4代核電站的實(shí)際應(yīng)用具有重要意義。目前針對(duì)sCO2的換熱特性已有大量的實(shí)驗(yàn)和模擬研究，一些學(xué)者通過實(shí)驗(yàn)研究了加熱條件下sCO2的傳熱惡化現(xiàn)象[2-3]。Wang等[4]通過實(shí)驗(yàn)研究了加熱條件下sCO2在垂直管中的對(duì)流換熱特性，并與水平管進(jìn)行比較，分析了壓力、熱流密度、質(zhì)量流速、入口溫度以及管徑對(duì)傳熱的影響，研究結(jié)果表明浮升力會(huì)強(qiáng)化垂直向下流動(dòng)的傳熱，惡化垂直向上流動(dòng)的傳熱。Kim等[5]通過實(shí)驗(yàn)分析了sCO2在向上流動(dòng)的管道中的傳熱特性，熱流密度大于某一值時(shí)管道會(huì)出現(xiàn)傳熱惡化現(xiàn)象，且環(huán)形通道能夠減弱惡化程度。Liu等[6]通過實(shí)驗(yàn)研究了浮升力對(duì)sCO2在豎直加熱管道內(nèi)流動(dòng)換熱的影響，實(shí)驗(yàn)結(jié)果表面浮升力在小于某一值時(shí)傳熱惡化消失，且傳熱惡化隨熱流密度的增大而向進(jìn)口處移動(dòng)。對(duì)于水平管，Liao等[7]通過實(shí)驗(yàn)研究得出在微細(xì)管中水平與垂直向下流動(dòng)的傳熱得到強(qiáng)化。Zhang等[8]同樣通過實(shí)驗(yàn)對(duì)不同流動(dòng)方向的4 mm加熱管進(jìn)行了研究，通過管道壁溫變化實(shí)驗(yàn)數(shù)據(jù)說明了水平方向與垂直向上流動(dòng)出現(xiàn)了傳熱惡化現(xiàn)象，但水平方向的惡化程度弱于垂直向上流動(dòng)。

另外，隨著計(jì)算流體力學(xué)的快速發(fā)展，數(shù)值模擬方法作為一種更加經(jīng)濟(jì)的研究方法被廣泛用于超臨界流體傳熱特性的研究中[9]。Cai等[10]通過數(shù)值模擬研究，表明sCO2在浮升力作用下產(chǎn)生的二次流可以增強(qiáng)管道頂部流體的傳熱并削弱管底部的傳熱。閆晨帥等[11]通過數(shù)值模擬的方法，解釋了加熱工況下頂母線和底母線傳熱出現(xiàn)明顯差異的機(jī)理，研究了水平管中“類液-類氣分層”現(xiàn)象對(duì)傳熱惡化的影響。Zhao等[12]對(duì)大直徑水平加熱管(22.14 mm)內(nèi)超臨界二氧化碳傳熱惡化現(xiàn)象進(jìn)行了數(shù)值模擬研究，并修正了更加適用于大直徑管的浮升力準(zhǔn)則數(shù)。此外，根據(jù)場協(xié)同原理[13]，對(duì)流換熱的增強(qiáng)不僅取決于物理性質(zhì)和流型，還取決于流體的速度矢量與溫度梯度之間的夾角。一些研究者利用場協(xié)同原理分析了sCO2管內(nèi)流動(dòng)換熱[14-15]。Xiang等[16]對(duì)水平圓管中sCO2的冷卻傳熱性能進(jìn)行了數(shù)值模擬研究，并利用場協(xié)同原理分析了截面上傳熱性能的差異。Yang等[17]分析了sCO2在帶有肋片的豎直管道內(nèi)的傳熱，證明了場協(xié)同原理可以成功地預(yù)測傳熱惡化或增強(qiáng)的程度。

綜合以上研究可知，sCO2管內(nèi)流動(dòng)換熱規(guī)律復(fù)雜，現(xiàn)有的實(shí)驗(yàn)與數(shù)值模擬研究并未完全闡明水平管中傳熱惡化的物理機(jī)制。本文通過數(shù)值模擬的方法研究了加熱工況下管徑對(duì)sCO2在水平管中的流動(dòng)換熱特性影響，并結(jié)合場協(xié)同理論，通過分析管內(nèi)流場、浮升力參數(shù)以及努塞爾數(shù)的變化，揭示了管徑對(duì)于浮升力效應(yīng)的影響，并量化了新渦流對(duì)于頂部的傳熱惡化程度。

1 水平換熱管物理模型與數(shù)值方法

1.1 幾何模型

本文的管道幾何模型如圖1所示，包括入口段和加熱段2部分，為了減小入口效應(yīng)對(duì)流動(dòng)換熱的影響，設(shè)有200 mm長的絕熱段作為入口段；加熱段為均勻熱流密度條件，長度為500 mm；為分析管徑對(duì)超臨界二氧化碳流動(dòng)換熱的影響，分別選用管徑為6、8和10 mm的管道。

圖1 水平換熱管物理模型Fig.1 Physical model of horizontal heat exchange tube

1.2 控制方程

穩(wěn)態(tài)湍流對(duì)流換熱的控制方程組在直角坐標(biāo)系中為[18]：

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

能量方程：

(3)

式中μeff是有效粘性系數(shù)，其計(jì)算公式為：

(4)

式中μt為湍流動(dòng)力粘度，其計(jì)算公式為：

(5)

(6)

(7)

式中：ρ為密度，kg/m3；u為速度，m/s；T為溫度，K；p為壓力，Pa；μ為動(dòng)力粘度，Pa·s；y為到最近壁面的距離；a1為常數(shù)，值為0.31；k為湍動(dòng)能，m2/s2。

本文基于Fluent商業(yè)CFD軟件，選擇SSTk-ω湍流模型作為數(shù)值計(jì)算中的湍流模型[19]。該模型通過使用混合函數(shù)，在近壁面處采用k-ω模型，其他區(qū)域采用k-ε模型。與標(biāo)準(zhǔn)k-ω模型相比，SSTk-ω湍流模型中使用了不同的湍流常數(shù)項(xiàng)，引入了橫向耗散導(dǎo)數(shù)項(xiàng)，并且在湍流粘度的定義中考慮了湍流剪應(yīng)力的輸運(yùn)過程，使得SSTk-ω模型適用性更高。Wang等[20]通過對(duì)比SSTk-ω模型與其他湍流模型得出SSTk-ω湍流模型能夠?qū)鳠嵯禂?shù)、壁面溫度和主流溫度進(jìn)行更好的預(yù)測，除此之外，也有許多研究采用了該湍流模型進(jìn)行研究并得到了準(zhǔn)確的結(jié)果[16,21-22]。對(duì)于SSTk-ω湍流模型，湍動(dòng)能k的運(yùn)輸方程為：

(8)

比耗散率ω的運(yùn)輸方程為：

Gω-Yω+Dω+Sω

(9)

式中：Gk和Gω為湍動(dòng)能和比耗散率的產(chǎn)生項(xiàng)；Гk和Гω為湍動(dòng)能和比耗散率的有效擴(kuò)散系數(shù)；Yk和Yω是湍動(dòng)能和比耗散率的耗散項(xiàng)；Dω為交叉擴(kuò)散項(xiàng)，Sk和Sω為湍動(dòng)能k和比耗散率ω的自定義源項(xiàng)。

1.3 網(wǎng)格劃分及計(jì)算方法

本文通過ANSYS ICEM軟件劃分模型網(wǎng)格，截面通過O型切分生成結(jié)構(gòu)化六面體網(wǎng)格，由于超臨界二氧化碳的物理性質(zhì)隨溫度發(fā)生劇烈變化，因此對(duì)近壁面處的網(wǎng)格進(jìn)行局部加密，由于使用了SSTk-ω湍流模型，因此需保持第1層網(wǎng)格y+<1，網(wǎng)格劃分如圖2所示。為確定網(wǎng)格數(shù)量對(duì)于計(jì)算結(jié)果的影響，需進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證。根據(jù)網(wǎng)格疏密程度不同，生成了5套網(wǎng)格，網(wǎng)格數(shù)量分別為659 610、1 074 920、1 790 370、3 015 360、4 149 610。不同網(wǎng)格下管道總壓降與流體出口溫度的計(jì)算結(jié)果如圖3所示，管道壓降隨網(wǎng)格數(shù)的增加而減小，出口溫度隨網(wǎng)格數(shù)的增加而增加。從圖中數(shù)據(jù)走勢可以看出，當(dāng)網(wǎng)格數(shù)量增大到一定程度后，壓降與出口溫度隨網(wǎng)格數(shù)量的進(jìn)一步增加基本趨于穩(wěn)定。綜合考慮計(jì)算資源和計(jì)算精度，本文選用第3套網(wǎng)格生成方案執(zhí)行相應(yīng)的計(jì)算任務(wù)。

圖2 水平管網(wǎng)格劃分Fig.2 Meshing diagram of horizontal tube

考慮到超臨界二氧化碳物性變化劇烈，選用FLUENT中的NIST 實(shí)際氣體模型以反映流體物性的劇烈變化給傳熱帶來的影響。模擬入口設(shè)置為質(zhì)量流速入口，出口為壓力出口，入口段設(shè)置為絕熱條件，加熱段為均勻熱流密度加熱。數(shù)值模擬計(jì)算采用壓力基、穩(wěn)態(tài)求解器進(jìn)行求解，求解方法采用基于壓力-速度耦合的COUPLED算法。壓力標(biāo)準(zhǔn)差分、連續(xù)性方程、動(dòng)量方程、能量方程、湍動(dòng)能以及比耗散率項(xiàng)均采用二階迎風(fēng)格式進(jìn)行離散。

圖3 水平管道模型網(wǎng)格獨(dú)立性驗(yàn)證Fig.3 Mesh independence verification of horizontal tube model

2 數(shù)值模擬計(jì)算結(jié)果及分析

2.1 數(shù)據(jù)處理方法

超臨界二氧化碳流體在管內(nèi)流動(dòng)的局部對(duì)流換熱系數(shù)h為：

(10)

(11)

格拉曉夫數(shù)Gr和主流區(qū)雷諾數(shù)Reb分別表示為：

(12)

(13)

努塞爾數(shù)表示為：

(14)

式中：q為壁面熱流密度，W/m2;T為溫度，K；A為管道換熱面積，m2；ρ為密度，kg/m3;u為速度，m/s；d為管道直徑，m；μ為動(dòng)力粘度，Pa·s；λ為導(dǎo)熱系數(shù)，W/(m·K)；下標(biāo)b表示主流區(qū)參數(shù)，下標(biāo)w表示壁面，下標(biāo)in表示流體入口，下標(biāo)out表示流體出口；Tw為加熱壁面平均溫度。

2.2 可靠性驗(yàn)證

為驗(yàn)證數(shù)值模型的準(zhǔn)確性，將數(shù)值模擬結(jié)果與Dang和Hihara的實(shí)驗(yàn)結(jié)果[23-24]進(jìn)行對(duì)比實(shí)驗(yàn)工況進(jìn)口溫度為350 K, 質(zhì)量流速為200 kg/(m2·s)，壁面熱流密度為12 kW/m2，壓力為8 MPa，實(shí)驗(yàn)測量結(jié)果與數(shù)值模擬計(jì)算結(jié)果如圖4所示?？梢钥闯觯M結(jié)果和實(shí)驗(yàn)結(jié)果的吻合較好，具有相同的變化趨勢，最大誤差發(fā)生在跨臨界區(qū)，誤差為18.84%，最小誤差為0.23%，平均誤差為7.49%，說明本文采用的數(shù)值方法準(zhǔn)確可行。

圖4 數(shù)值方法驗(yàn)證Fig.4 Numerical method validation

2.3 管徑對(duì)換熱效果的影響

為研究加熱工況下管道直徑大小對(duì)sCO2管內(nèi)流動(dòng)換熱特性的影響，分別對(duì)3種不同管徑的管道進(jìn)行數(shù)值模擬分析。管道直徑取為6、8和10 mm，入口質(zhì)量流量為0.006 kg/s，加熱段的壁面熱流密度為50 kW/m2，流體入口壓力為7.5 MPa，入口溫度為310 K。表1為3種管徑對(duì)應(yīng)的模擬結(jié)果，在相同的入口質(zhì)量流率下，管徑越大，管內(nèi)流速越小，而對(duì)流傳熱系數(shù)受到管內(nèi)流速的影響，因此3種管徑的平均傳熱系數(shù)相差較大。為了研究管徑對(duì)傳熱效果的影響，本文通過改變重力加速度的大小來分析不同管徑中浮升力在加熱情況下對(duì)整體傳熱性能的影響。浮升力是由于截面內(nèi)密度分布不均勻，在重力的影響下，密度較大的流體向下運(yùn)動(dòng)，而密度較小的流體則向上運(yùn)動(dòng)，在截面內(nèi)形成二次環(huán)流，在重力加速度為0時(shí)截面內(nèi)則不會(huì)出現(xiàn)二次環(huán)流現(xiàn)象。從表中數(shù)據(jù)可以看出，重力加速度為0.5g時(shí)，管徑為6、8和10 mm的管內(nèi)平均傳熱系數(shù)比0g時(shí)分別增加了3.45%、12.86%和23.52%；重力加速度為g時(shí)，管徑為6、8和10 mm的管內(nèi)平均傳熱系數(shù)比0g時(shí)分別增加了4.80%、15.73%和25.84%。模擬結(jié)果表明，在本文選定實(shí)驗(yàn)工況范圍內(nèi)，浮升力效應(yīng)會(huì)強(qiáng)化管道整體的換熱效果；并且在等質(zhì)量流率條件下，管徑增大會(huì)讓浮升力對(duì)于整個(gè)管道的流動(dòng)換熱特性影響增大；同時(shí)也可以看出重力加速度增大后，浮升力對(duì)于整個(gè)管道的流動(dòng)換熱特性影響也隨之增大。

表1 3種管徑的傳熱效果數(shù)值模擬結(jié)果

2.4 管徑對(duì)流場的影響

圖5為不同管徑對(duì)應(yīng)的密度分布云圖和比熱容分布云圖。受重力影響，sCO2在水平管道內(nèi)流動(dòng)時(shí)呈現(xiàn)出上下分層的流動(dòng)形態(tài)，頂部流體密度較小，底部流體密度較大，其流體物性分別與氣態(tài)和液態(tài)sCO2相近，也可稱為“類氣”和“類液”sCO2[11]。由比熱容云圖分布可以看出，管道頂部“類氣”流態(tài)的sCO2不僅密度較小，其比熱容同樣較低，隨管徑增大管道流通截面積增加，相同質(zhì)量流量下sCO2工質(zhì)的平均流速降低，管道頂部積聚的低密度“類氣”流體相應(yīng)增多，超臨界流體上下分層愈加明顯，其聚集在管道頂部會(huì)使得壁面局部傳熱惡化現(xiàn)象更加嚴(yán)重。

圖5 不同管徑密度和比熱容分布Fig.5 Density and specific heat capacity distribution of different pipe diameters

2.5 浮升力對(duì)傳熱性能的影響

(15)

圖6 浮升力參數(shù)沿流向變化Fig.6 Variation of buoyancy parameters along the flow direction

為研究不同管道直徑下浮升力效應(yīng)的影響程度，進(jìn)一步分析了壁面局部傳熱系數(shù)隨流動(dòng)方向的變化，如圖7所示。有重力時(shí)的平均傳熱系數(shù)和底母線上局部傳熱系數(shù)均大于無重力時(shí)的傳熱系數(shù)，但頂母線的傳熱系數(shù)小于無重力時(shí)的傳熱系數(shù)，說明在加熱情況下傳熱強(qiáng)化與傳熱惡化同時(shí)存在，在重力作用下產(chǎn)生的浮升力，強(qiáng)化了整個(gè)管道和底母線的傳熱效果，但同時(shí)也惡化了管道壁面頂母線附近的傳熱；通過對(duì)比不同管徑的傳熱情況可發(fā)現(xiàn)，隨著管徑的增大，頂母線和底母線的傳熱系數(shù)相差越來越大，表明管徑增加引起的浮升力增強(qiáng)，會(huì)同時(shí)加劇傳熱強(qiáng)化和傳熱惡化這2種現(xiàn)象。

圖7 不同管徑表面?zhèn)鳠嵯禂?shù)沿流向變化Fig.7 Variation of heat transfer coefficient of different pipe diameters along the flow direction

為分析管徑對(duì)管道頂部傳熱惡化的影響，本文取Z=375 mm截面處的流線及密度分布云圖進(jìn)行分析，如圖8所示。隨著浮升力的增大，截面內(nèi)流體的密度分布更加不均勻，并且在管徑為8和10 mm時(shí)，管道頂部出現(xiàn)了新的渦流，而這部分低密度流體導(dǎo)熱能力較差且不參與主流區(qū)的流動(dòng)，聚集在管道頂部使得傳熱惡化進(jìn)一步加劇。

圖9為頂母線Nutop和周向平均Nuave與無重力工況下Nuzero的比值。在加熱時(shí)，頂母線的傳熱由于浮升力的出現(xiàn)而發(fā)生惡化，而無重力工況時(shí)沒有浮升力，本文用努塞爾數(shù)的比值來評(píng)估頂母線附近流體傳熱的惡化程度，該比值越小，頂母線的傳熱能力越差，說明傳熱惡化越劇烈。由頂母線Nu比值可以看出，管徑為8和10 mm的比值要遠(yuǎn)小于管徑為6 mm時(shí)的Nu比值，這與圖8中管道頂部新渦流的出現(xiàn)有關(guān)。管徑為6 mm時(shí)，Nu的平均值為0.55，管徑為8和10 mm時(shí)平均值分別為0.35和0.33，說明了新渦流的出現(xiàn)會(huì)使管道頂部傳熱惡化進(jìn)一步加劇，Nu的比值降低了約0.2。從周向平均Nu比值可以看出，管徑越大，浮升力對(duì)管內(nèi)傳熱影響越大，這與上文所得結(jié)論一致，在加熱時(shí)管內(nèi)傳熱惡化和傳熱強(qiáng)化現(xiàn)象同時(shí)存在，管道整體傳熱效果得到強(qiáng)化，而頂部則會(huì)發(fā)生劇烈的傳熱惡化現(xiàn)象。

圖8 不同管徑下Z=375 mm截面處的流線及密度分布云圖Fig.8 Streamline diagram and density distribution at Z=375 mm under different pipe diameters

2.6 場協(xié)同分析

場協(xié)同原理[13]從流體速度場和溫度場相互作用的角度研究了對(duì)流傳熱機(jī)理，認(rèn)為對(duì)流傳熱的強(qiáng)化或弱化不僅和物性、流動(dòng)狀態(tài)有關(guān)，還與速度矢量和溫度梯度場的協(xié)同性有關(guān)。局部場協(xié)同角θ為：

(16)

式中：V為流體的速度矢量；▽T為流體的溫度梯度。場協(xié)同角表示速度場和溫度場之間的協(xié)同程度，場協(xié)同角為0°時(shí)，場協(xié)同程度最好。圖10為不同管徑對(duì)應(yīng)的局部場協(xié)同角分布圖。

圖9 頂母線、周向平均與無重力工況的Nu比值Fig.9 The Nu ratio of top bus to zero gravity and circumferential average to zero gravity

圖10 局部場協(xié)同角分布Fig.10 Local field synergy angle distribution

本文沿流向分別取Z=125,250,375 mm 3個(gè)截面進(jìn)行分析，具有較小場協(xié)同角的流體主要聚集在管道的底部，說明管道底部的傳熱得到了強(qiáng)化，同時(shí)結(jié)合圖7、8進(jìn)行分析，超臨界二氧化碳在加熱時(shí)傳熱強(qiáng)化與惡化現(xiàn)象同時(shí)出現(xiàn)在管道底部和頂部；隨著管徑的增加，場協(xié)同角分布變得更不均勻，具有較小場協(xié)同角的流體更加靠近壁面；并且可以看出在管徑增大后，管道下游的局部場協(xié)同角增大，接近90°，即速度場和溫度梯度的夾角趨近于90°，這是由于管徑增大后，x方向與y方向的速度分量減小，即截面內(nèi)主渦流的強(qiáng)度減小，因此在管道底部流體的局部場協(xié)同角隨之增大。

3 結(jié)論

1) 浮升力效應(yīng)會(huì)強(qiáng)化管道整體的換熱效果，并且在等質(zhì)量流率條件下，管徑越大，浮升力對(duì)于整個(gè)管道的流動(dòng)換熱特性影響越大。

2) 隨著管徑的增大，管道頂部積聚的低密度“類氣”流體相應(yīng)增多，超臨界流體上下分層愈加明顯，頂母線和底母線的傳熱系數(shù)相差越來越大，表明管徑增加引起的浮升力增強(qiáng)，會(huì)同時(shí)加劇傳熱強(qiáng)化和傳熱惡化這2種現(xiàn)象。并且管徑較大時(shí)，“類氣”流體會(huì)在管道頂部形成新的渦流，本文通過努塞爾數(shù)的比值得出新渦流的出現(xiàn)會(huì)使管道頂部傳熱惡化進(jìn)一步加劇。

3) 通過場協(xié)同原理分析了管內(nèi)速度場與溫度場的協(xié)同程度，管徑增大后會(huì)使主渦流強(qiáng)度減小，并且場協(xié)同角分布更加不均勻，具有較小場協(xié)同較的流體更加靠近壁面。