超臨界壓力下二氧化碳在水平管內(nèi)傳熱異常特性研究

宿詩(shī)雨， 姜文全， 李 琳， 石杰峰， 李 洋， 楊 帆

(1． 遼寧石油化工大學(xué) 石油天然氣工程學(xué)院， 遼寧撫順 113001；2． 遼寧石油化工大學(xué) 機(jī)械工程學(xué)院， 遼寧撫順 113001)

CO2因具有無(wú)毒、穩(wěn)定性強(qiáng)、對(duì)環(huán)境影響小等特點(diǎn)，近年來(lái)在能源發(fā)電、核反應(yīng)堆、航空航天、熱泵與制冷等領(lǐng)域有較好的應(yīng)用前景[1-2]。與傳統(tǒng)蒸汽朗肯循環(huán)相比，超臨界壓力下的CO2循環(huán)系統(tǒng)具有熱功轉(zhuǎn)換效率高、靈活性強(qiáng)和經(jīng)濟(jì)性好的優(yōu)勢(shì)，近年來(lái)得到了廣泛的商業(yè)應(yīng)用，也為換熱裝置的設(shè)計(jì)及系統(tǒng)安全運(yùn)行提供了理論基礎(chǔ)[3]。超臨界壓力下的CO2在壓力p=8 MPa時(shí)的熱物性如圖1所示。由圖1可知，其熱物性在溫度T=307.8 K附近發(fā)生劇烈變化，尤其是比熱容，在此溫度點(diǎn)達(dá)到極大值，將此點(diǎn)稱為擬臨界溫度點(diǎn)Tpc。

圖1 p=8 MPa時(shí)的CO2熱物性

超臨界壓力下CO2的傳熱可分為正常傳熱、傳熱強(qiáng)化和傳熱惡化3種模式。一般來(lái)說(shuō)，將傳熱系數(shù)出現(xiàn)峰值的情況視為傳熱強(qiáng)化，將壁溫飛升現(xiàn)象視為傳熱惡化[4-6]。在已有文獻(xiàn)中，超臨界壓力下CO2大多數(shù)是在豎直管中開(kāi)展實(shí)驗(yàn)與模擬研究的。在豎直流動(dòng)中，流場(chǎng)對(duì)稱，管壁的傳熱能力相同，但受浮升力效應(yīng)影響，向上流動(dòng)和向下流動(dòng)不同。一般來(lái)說(shuō)，在加熱條件下，向下流動(dòng)的傳熱能力比向上流動(dòng)強(qiáng)[7-9]。

超臨界壓力下CO2在水平管中的流動(dòng)相比在豎直管中的流動(dòng)更加復(fù)雜。在水平流動(dòng)中，因重力產(chǎn)生的浮升力發(fā)生在流動(dòng)的垂直方向上，導(dǎo)致水平管上下壁傳熱特性出現(xiàn)差異。眾多學(xué)者對(duì)水平管壁傳熱特性不同的解釋也不盡相同：Adebiyi等[10]通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)水平管頂部與底部傳熱能力不同，并發(fā)現(xiàn)底部傳熱效果較好；Liao等[11]通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)當(dāng)雷諾數(shù)(Re)高達(dá)105時(shí)，浮升力在水平流動(dòng)中仍起重要作用；相夢(mèng)如等[12]通過(guò)場(chǎng)協(xié)同原理解釋了上下壁面換熱不均勻的原因；閆晨帥等[13]著重研究水平管上壁的傳熱特性，發(fā)現(xiàn)在質(zhì)量流量較低時(shí)，周向上壁較低的密度及湍動(dòng)能與軸向較小的速度共同導(dǎo)致上壁出現(xiàn)高溫現(xiàn)象；Kim等[14]通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，現(xiàn)有的浮升力判別因子臨界值不適用于低熱流密度工況，所以需要新的浮升力因子對(duì)水平流動(dòng)作進(jìn)一步評(píng)價(jià)。綜上可得，浮升力與物性是造成水平管管壁傳熱特性不同的主要因素，且還需對(duì)水平管內(nèi)浮升力的判定標(biāo)準(zhǔn)開(kāi)展深入的研究。

筆者對(duì)超臨界壓力下CO2在4 mm水平管內(nèi)傳熱過(guò)程進(jìn)行數(shù)值模擬，分析了熱流密度、壓力和浮升力對(duì)傳熱的影響，通過(guò)對(duì)比確定了更適合本文工況下的浮升力標(biāo)準(zhǔn)，并根據(jù)截面的物性分布獲得了水平管壁不同位置的傳熱特性。

1 物理模型與數(shù)值模擬方法

1.1 物理模型

三維物理模型如圖2所示。水平管直徑(d)為4 mm，管長(zhǎng)(L)為1 000 mm，計(jì)算時(shí)忽略壁厚。水平管周向角度等于0°時(shí)定義為上母線(top)，周向角度為180°時(shí)定義為下母線(bottom)。

圖2 水平管物理模型

1.2 網(wǎng)格劃分與數(shù)值方法

計(jì)算區(qū)域內(nèi)利用ICEM軟件劃分網(wǎng)格，圓形截面采用O型剖分保證網(wǎng)格質(zhì)量，近壁面處網(wǎng)格進(jìn)行加密處理，且第1層網(wǎng)格無(wú)量綱高度y+<1，如圖3所示。采用FLUENT軟件進(jìn)行數(shù)值模擬，其中定壓力和定溫度下的物性參數(shù)通過(guò)REFPROP軟件計(jì)算得到，并以線性插值方法輸入。

圖3 網(wǎng)格劃分截面

邊界條件如下：入口質(zhì)量流量G為140～260 kg/(m2·s)；出口壓力p為8 MPa；加熱壁面恒熱流密度(即壁面熱流密度)qw為20～40 kW/m2；入口溫度為 302 K。采用Standardk-ε增強(qiáng)壁面函數(shù)湍流模型，壓力與速度耦合采用SIMPLEC算法，當(dāng)連續(xù)性方程的殘差小于10-5時(shí)，認(rèn)為計(jì)算收斂。采用的數(shù)值模型中包括連續(xù)性方程、動(dòng)量方程和能量方程，在三維直角坐標(biāo)系下分別用以下方程式表示。

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

能量方程：

(3)

式中：u為流體速度；E為比內(nèi)能；μt為湍流黏度；g為重力加速度；φ為黏性引起的能量耗散；下標(biāo)i、j為方向。

1.3 網(wǎng)格無(wú)關(guān)性與模型驗(yàn)證

p=8 MPa，G=260 kg/(m2·s)，qw=35 kW/m2時(shí)不同網(wǎng)格節(jié)點(diǎn)數(shù)的上母線軸向壁溫(Tw)見(jiàn)圖4。由圖4可知，當(dāng)網(wǎng)格節(jié)點(diǎn)數(shù)大于329萬(wàn)時(shí)，壁溫曲線幾乎無(wú)偏差，對(duì)后續(xù)計(jì)算結(jié)果影響不大，所以最終選定網(wǎng)格節(jié)點(diǎn)數(shù)為329萬(wàn)。

圖4 不同網(wǎng)格節(jié)點(diǎn)數(shù)的上母線軸向壁溫

為了驗(yàn)證數(shù)值計(jì)算的合理性和準(zhǔn)確性，選用Adebiyi等[10]的實(shí)驗(yàn)工況(p=7.59 MPa,G=384.4 kg/(m2·s)，qw=15.1 kW/m2)進(jìn)行模型驗(yàn)證。母線溫度模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比如圖5所示。由圖5可知，模擬結(jié)果與實(shí)驗(yàn)結(jié)果的變化趨勢(shì)基本一致，因此所采用的數(shù)值模型可以較好地預(yù)測(cè)水平管中超臨界壓力下CO2的傳熱特性。

圖5 母線溫度模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比

2 計(jì)算結(jié)果與討論

2.1 熱流密度對(duì)傳熱的影響

為了獲得不同運(yùn)行參數(shù)下水平管內(nèi)超臨界CO2的傳熱特性，本節(jié)分別考察了2種質(zhì)量流量下不同熱流密度對(duì)管壁傳熱的影響。

對(duì)流傳熱系數(shù)hw定義為：

(4)

(5)

式中：Tb為主流溫度；qw為壁面熱流密度；A為水平管橫截面面積；i為當(dāng)?shù)亓黧w焓值。

不同熱流密度下壁溫與傳熱系數(shù)隨主流焓值的變化曲線如圖6所示。

由圖6(a)可知，上母線的壁溫在入口處突增，出現(xiàn)尖銳峰值，即發(fā)生局部傳熱惡化，隨后壁溫下降，在擬臨界焓值(ipc)附近時(shí)趨于平緩；下母線的壁溫隨著主流焓值ib的增大而單調(diào)上升。隨著熱流密度的增大，上母線的壁溫峰值逐漸上升，與下母線的壁溫差也隨之增大。當(dāng)qw=40 kW/m2時(shí)，上下母線壁溫差達(dá)到186 K。在圖6(a)中，Tb為302～378 K，Tw為304～520 K。流體受重力影響產(chǎn)生浮升力，浮升力對(duì)換熱的影響隨著熱流密度的增大而增強(qiáng)，導(dǎo)致上下母線出現(xiàn)較大溫差。

由圖6(b)可知，傳熱系數(shù)與溫度的變化趨勢(shì)相反，上母線傳熱系數(shù)(hw,t)始終低于下母線的傳熱系數(shù)(hw,b)，上母線傳熱系數(shù)在壁溫峰值處達(dá)到谷值，隨著主流焓值的增大，換熱逐漸恢復(fù)，并在擬臨界焓值附近達(dá)到最大值。熱流密度越大，傳熱惡化程度越高，傳熱系數(shù)越小。隨著主流焓值ib的增大，下母線傳熱系數(shù)呈先增大后減小的趨勢(shì)，在主流焓值低于擬臨界焓值處出現(xiàn)傳熱系數(shù)峰值，即發(fā)生傳熱強(qiáng)化，其出現(xiàn)峰值的位置比上母線更早，隨著熱流密度增大，傳熱系數(shù)減小，峰值降低。

由圖6(c)和圖6(d)可知，在G=260 kg/(m2·s)條件下，上母線壁溫同樣在入口處達(dá)到峰值，對(duì)應(yīng)的傳熱系數(shù)位于最低點(diǎn)，隨后換熱恢復(fù)，在ib≈330 kJ/kg處達(dá)到峰值。在圖6(c)中，Tb為302～320 K，Tw為304～377 K。隨著熱流密度的增大，換熱恢復(fù)程度降低，傳熱系數(shù)減小。下母線的壁溫和傳熱系數(shù)與G=140 kg/(m2·s)時(shí)的變化規(guī)律相同，不再贅述。與G=140 kg/(m2·s)相比，G=260 kg/(m2·s)時(shí)的上下母線壁溫差減小，浮升力效應(yīng)減弱，上母線傳熱惡化現(xiàn)象有所緩解。

2.2 浮升力對(duì)傳熱的影響

上下母線的溫度分布不均勻可充分說(shuō)明浮升力對(duì)換熱的影響，針對(duì)水平管內(nèi)浮升力效應(yīng)，目前仍沒(méi)有統(tǒng)一的判別式和閾值。呂海財(cái)?shù)萚15]通過(guò)對(duì)比發(fā)現(xiàn)Kim等[16]提出的浮升力(BK)判別式能較好地反映超臨界壓力下CO2在小管徑通道內(nèi)的傳熱過(guò)程；蒲星宇等[17]認(rèn)為判斷浮升力產(chǎn)生的原因應(yīng)分為qw/G較小和qw/G較大這2種情況，計(jì)算浮升力時(shí)應(yīng)選取上母線溫度，并得出BP[18]判別式更適合壁溫差較大工況的結(jié)論。選取BK和BP來(lái)評(píng)估管內(nèi)浮升力效應(yīng)，計(jì)算公式如下：

(6)

(7)

(8)

(9)

(10)

(11)

當(dāng)BP>1時(shí)，不可忽略浮升力。BP公式構(gòu)成中包含d4，管徑對(duì)浮升力產(chǎn)生較大影響。在已有文獻(xiàn)中，普遍得出了大管徑模型中浮升力較強(qiáng)的結(jié)論[19]。對(duì)于小管徑的幾何模型，浮升力對(duì)傳熱也有影響[20]。G=140 kg/(m2·s)和G=260 kg/(m2·s)時(shí)，2種判別式下上母線浮升力隨熱流密度變化的曲線如圖7所示，其中ΔT為上下母線壁溫差。由圖7可知，2種判別式下浮升力均隨熱流密度的增大而增大，但曲線趨勢(shì)不同。BK的變化趨勢(shì)為浮升力峰值出現(xiàn)在擬臨界焓值處，而B(niǎo)P的變化趨勢(shì)為浮升力在發(fā)生傳熱惡化的位置處出現(xiàn)尖銳峰值，然后下降，直到加熱結(jié)束。另外，浮升力在擬臨界焓值點(diǎn)前下降較快，而在擬臨界焓值點(diǎn)后下降較平緩。

(a) G=140 kg/(m2·s)

比較圖7(a)與圖7(b)發(fā)現(xiàn)，G=140 kg/(m2·s)時(shí)的浮升力更強(qiáng)，峰值更明顯，所以壁溫差越大浮升力越強(qiáng)，且浮升力在峰值后的下降速度更快。而G=260 kg/(m2·s)時(shí)，3種熱流密度下的BP均小于1，所以此工況下的浮升力可忽略。

浮升力與壁溫差的曲線趨勢(shì)大致相同，壁溫差的峰值可大致代表浮升力的峰值。所以在本文2種工況下，BP更能反映浮升力的分布，壁溫差大小可反映浮升力的強(qiáng)弱。

2.3 壓力對(duì)浮升力的影響

在G=140 kg/(m2·s)、qw=30 kW/m2條件下，壓力分別為8.0 MPa、8.5 MPa和9.0 MPa時(shí)的浮升力及CO2物性參數(shù)如圖8所示。由圖8可以看出，浮升力隨著壓力的增大而下降，對(duì)傳熱的影響也逐漸減弱。這是由于當(dāng)壓力遠(yuǎn)離臨界壓力時(shí)，CO2的物性劇烈變化程度減小。從圖8可知，物性隨著壓力的變化而變化，比熱容在擬臨界溫度點(diǎn)附近突升幅度隨壓力的增大而減小，而壓力對(duì)傳熱的影響主要取決于比熱容。因此，增大壓力可抑制傳熱惡化，使浮升力效應(yīng)減弱。

圖8 不同壓力下浮升力隨主流焓值的變化曲線

2.4 水平管傳熱特性機(jī)理

選取G=140 kg/(m2·s)，qw=30 kW/m2的工況，取ib=280 kJ/kg(入口處)、287.3 kJ/kg(壁溫峰值處)、320 kJ/kg、340 kJ/kg(擬臨界焓值)、360 kJ/kg、400 kJ/kg和440 kJ/kg 7個(gè)典型截面，各截面的密度、溫度、湍動(dòng)能(k)及速度分布如圖9所示。

由圖9(a)和圖9(b)可知，上母線聚集著高溫且密度小的流體，產(chǎn)生了溫度梯度，而下母線處的流體因重力作用密度較大，溫度較低，其中在ib=287.3 kJ/kg處的壁溫差和密度差最大，對(duì)應(yīng)的截面密度分布如圖10所示。由圖10可知，因浮升力效應(yīng)各截面密度出現(xiàn)明顯的分層現(xiàn)象，隨著加熱的進(jìn)行，總密度逐漸減小，分層效果逐漸減弱，上下母線密度差也逐漸減小。

由圖9(c)和圖9(d)可知，湍動(dòng)能和速度不完全呈軸對(duì)稱，湍動(dòng)能隨著加熱的進(jìn)行逐漸增大，頂部湍動(dòng)能稍小于底部湍動(dòng)能，底部湍動(dòng)能強(qiáng)化了下母線的傳熱。ib=287.3 kJ/kg處的頂部密度低，所以對(duì)應(yīng)的頂部速度較小，在換熱恢復(fù)期間，流體速度逐漸增大，且速度梯度增加，上下母線的速度逐漸接近，在加熱即將結(jié)束時(shí)，頂部和底部的速度接近相同。

3 結(jié) 論

(1) 熱流密度的增大可導(dǎo)致傳熱惡化的發(fā)生。熱流密度越大，傳熱惡化現(xiàn)象越明顯，同時(shí)上下母線的壁溫差會(huì)隨熱流密度的增大而增大。

(2) 在G=140 kg/(m2·s)條件下，浮升力導(dǎo)致上下母線壁溫差較大；在G=260 kg/(m2·s)條件下，浮升力可忽略。BP判別式更能反映本文工況下浮升力的分布，同時(shí)壁溫差大小可大致反映浮升力的強(qiáng)弱。

(3) 壓力增大使CO2物性變化趨于平緩，浮升力隨著壓力的增大而下降，因此增大壓力可抑制傳熱惡化。

(4) 由典型截面物性分布可知，在發(fā)生傳熱惡化處(ib=287.3 kJ/kg)，密度分層明顯且溫度梯度較大，上母線的速度和湍動(dòng)能稍低于下母線，反映了上下母線不同的傳熱機(jī)理。