不同管徑下超臨界CO2流動(dòng)傳熱的數(shù)值分析

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

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

超臨界二氧化碳sCO2布雷頓循環(huán)在核電領(lǐng)域具有諸多優(yōu)勢(shì)，循環(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é)者通過(guò)實(shí)驗(yàn)研究了加熱條件下sCO2的傳熱惡化現(xiàn)象[2-3]。Wang等[4]通過(guò)實(shí)驗(yàn)研究了加熱條件下sCO2在垂直管中的對(duì)流換熱特性，并與水平管進(jìn)行比較，分析了壓力、熱流密度、質(zhì)量流速、入口溫度以及管徑對(duì)傳熱的影響，研究結(jié)果表明浮升力會(huì)強(qiáng)化垂直向下流動(dòng)的傳熱，惡化垂直向上流動(dòng)的傳熱。Kim等[5]通過(guò)實(shí)驗(yàn)分析了sCO2在向上流動(dòng)的管道中的傳熱特性，熱流密度大于某一值時(shí)管道會(huì)出現(xiàn)傳熱惡化現(xiàn)象，且環(huán)形通道能夠減弱惡化程度。Liu等[6]通過(guò)實(shí)驗(yàn)研究了浮升力對(duì)sCO2在豎直加熱管道內(nèi)流動(dòng)換熱的影響，實(shí)驗(yàn)結(jié)果表面浮升力在小于某一值時(shí)傳熱惡化消失，且傳熱惡化隨熱流密度的增大而向進(jìn)口處移動(dòng)。對(duì)于水平管，Liao等[7]通過(guò)實(shí)驗(yàn)研究得出在微細(xì)管中水平與垂直向下流動(dòng)的傳熱得到強(qiáng)化。Zhang等[8]同樣通過(guò)實(shí)驗(yàn)對(duì)不同流動(dòng)方向的4 mm加熱管進(jìn)行了研究，通過(guò)管道壁溫變化實(shí)驗(yàn)數(shù)據(jù)說(shuō)明了水平方向與垂直向上流動(dòng)出現(xiàn)了傳熱惡化現(xiàn)象，但水平方向的惡化程度弱于垂直向上流動(dòng)。

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

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

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

1.1 幾何模型

本文的管道幾何模型如圖1所示，包括入口段和加熱段2部分，為了減小入口效應(yīng)對(duì)流動(dòng)換熱的影響，設(shè)有200 mm長(zhǎng)的絕熱段作為入口段；加熱段為均勻熱流密度條件，長(zhǎng)度為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]。該模型通過(guò)使用混合函數(shù)，在近壁面處采用k-ω模型，其他區(qū)域采用k-ε模型。與標(biāo)準(zhǔn)k-ω模型相比，SSTk-ω湍流模型中使用了不同的湍流常數(shù)項(xiàng)，引入了橫向耗散導(dǎo)數(shù)項(xiàng)，并且在湍流粘度的定義中考慮了湍流剪應(yīng)力的輸運(yùn)過(guò)程，使得SSTk-ω模型適用性更高。Wang等[20]通過(guò)對(duì)比SSTk-ω模型與其他湍流模型得出SSTk-ω湍流模型能夠?qū)鳠嵯禂?shù)、壁面溫度和主流溫度進(jìn)行更好的預(yù)測(cè)，除此之外，也有許多研究采用了該湍流模型進(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ì)算方法

本文通過(guò)ANSYS ICEM軟件劃分模型網(wǎng)格，截面通過(guò)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)格無(wú)關(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ù)走勢(shì)可以看出，當(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í)際氣體模型以反映流體物性的劇烈變化給傳熱帶來(lái)的影響。模擬入口設(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)測(cè)量結(jié)果與數(shù)值模擬計(jì)算結(jié)果如圖4所示?？梢钥闯?，模擬結(jié)果和實(shí)驗(yàn)結(jié)果的吻合較好，具有相同的變化趨勢(shì)，最大誤差發(fā)生在跨臨界區(qū)，誤差為18.84%，最小誤差為0.23%，平均誤差為7.49%，說(shuō)明本文采用的數(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ì)傳熱效果的影響，本文通過(guò)改變重力加速度的大小來(lái)分析不同管徑中浮升力在加熱情況下對(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ì)流場(chǎng)的影響

圖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ù)均大于無(wú)重力時(shí)的傳熱系數(shù)，但頂母線的傳熱系數(shù)小于無(wú)重力時(shí)的傳熱系數(shù)，說(shuō)明在加熱情況下傳熱強(qiáng)化與傳熱惡化同時(shí)存在，在重力作用下產(chǎn)生的浮升力，強(qiáng)化了整個(gè)管道和底母線的傳熱效果，但同時(shí)也惡化了管道壁面頂母線附近的傳熱；通過(guò)對(duì)比不同管徑的傳熱情況可發(fā)現(xiàn)，隨著管徑的增大，頂母線和底母線的傳熱系數(shù)相差越來(lái)越大，表明管徑增加引起的浮升力增強(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與無(wú)重力工況下Nuzero的比值。在加熱時(shí)，頂母線的傳熱由于浮升力的出現(xiàn)而發(fā)生惡化，而無(wú)重力工況時(shí)沒(méi)有浮升力，本文用努塞爾數(shù)的比值來(lái)評(píng)估頂母線附近流體傳熱的惡化程度，該比值越小，頂母線的傳熱能力越差，說(shuō)明傳熱惡化越劇烈。由頂母線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，說(shuō)明了新渦流的出現(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 場(chǎng)協(xié)同分析

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

(16)

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

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

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

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

3 結(jié)論

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

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

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