非均勻加熱管內(nèi)超臨界CO2傳熱特性研究

吳家榮，李紅智，楊 玉，韓萬龍

（西安熱工研究院有限公司，陜西 西安 710054）

超臨界CO2布雷頓循環(huán)(SCO2-BC)因布局簡單靈活，循環(huán)效率高，設(shè)備緊湊，工質(zhì)無毒、廉價(jià)等優(yōu)勢在核電、煤電、余熱回收和艦船推進(jìn)等多個(gè)領(lǐng)域得到了越來越多的研究[1-4]。作為SCO2-BC 發(fā)電系統(tǒng)的關(guān)鍵環(huán)節(jié)，超臨界CO2在鍋爐氣冷壁內(nèi)的傳熱流動規(guī)律關(guān)系著整個(gè)系統(tǒng)的安全和效率。不同于亞臨界流體，超臨界CO2的物性隨著壓力和溫度的變化呈現(xiàn)強(qiáng)烈的非線性變化，尤其在擬臨界點(diǎn)附近，流體的比定壓熱容和導(dǎo)熱系數(shù)會發(fā)生明顯突變，如圖1所示。物性變化及其引起的自然對流使流動狀態(tài)由單一的強(qiáng)制對流發(fā)展為復(fù)雜的混合對流，伴隨而來的是傳熱惡化和強(qiáng)化的發(fā)生，這給鍋爐氣冷壁的熱力設(shè)計(jì)造成了困難。

目前，比定壓熱容和導(dǎo)熱系數(shù)的增大被認(rèn)為是傳熱強(qiáng)化的主要原因；而關(guān)于傳熱惡化，一種觀點(diǎn)是Ackerman 等人[5]根據(jù)超臨界水傳熱提出的擬膜態(tài)沸騰理論，另一種則是Jackson 等人[6-8]提出的傳熱惡化與浮力效應(yīng)和加速效應(yīng)對湍流的產(chǎn)生和傳熱有關(guān)。

圖1 CO2 物性Fig.1 Physical properties of CO2

針對上述觀點(diǎn)，各國學(xué)者開展了大量超臨界流體的傳熱研究。但是，這些研究[9-12]絕大多數(shù)都是在均勻熱流條件下開展的，這與氣冷壁在爐膛內(nèi)單側(cè)受熱的實(shí)際情況不符。因此，也有學(xué)者開始重點(diǎn)研究超臨界流體在非均勻熱流條件下的流動傳熱。Bai 等人[13]比較了超臨界水在均勻和非均勻加熱圓管內(nèi)不同母線處的壁溫、溫差、對流傳熱系數(shù)隨主流焓的變化，模擬結(jié)果表明，非均勻加熱時(shí)，圓管周向壁溫的分布是不均勻的，傳熱強(qiáng)化也僅發(fā)生在局部區(qū)域，某些母線處甚至出現(xiàn)流體溫度高于壁溫的情況。Fan 等人[14]利用數(shù)值模擬對高質(zhì)量流速下超臨界CO2在豎直非均勻加熱圓管中的傳熱特性進(jìn)行了研究，結(jié)果表明，熱流密度很大的頂部母線處浮力和加速效應(yīng)對傳熱的影響遠(yuǎn)大于熱流密度較小的底部母線處。強(qiáng)化傳熱的主要原因是擬臨界區(qū)比定壓熱容的增大；而傳熱惡化除受浮力和加速效應(yīng)的影響外，還源于管內(nèi)黏性底層的增厚。Gu 等人[15]比較了超臨界水在內(nèi)螺紋管內(nèi)均勻和非均勻受熱的情況，數(shù)值模擬結(jié)果顯示，非均勻受熱時(shí)，因管內(nèi)旋流的影響，流體混合更好，使得管內(nèi)最高內(nèi)壁溫比均勻受熱時(shí)的最高內(nèi)壁溫低了10 K，且最高內(nèi)壁溫并不總出現(xiàn)在受熱面的中點(diǎn)。上述對超臨界流體在非均勻熱流條件下的研究都集中在內(nèi)徑＞ 10 mm、質(zhì)量流量較大的情況，對小管徑情況下超臨界CO2傳熱特性的研究還存在不足。

綜上所述，本文開展了超臨界CO2在非均勻熱流密度下圓管內(nèi)傳熱的數(shù)值模擬研究。所選取的參數(shù)為：質(zhì)量流速100 kg/(m2·s)、熱流密度28～50 kW/m2、壓力7.6～8.6 MPa。通過在加熱段一側(cè)施加熱流以模擬氣冷壁受熱情況，比較了均勻和非均勻加熱的差別，分析了非均勻加熱情況下熱流密度和入口壓力對對流傳熱系數(shù)、浮力效應(yīng)和加速效應(yīng)的影響，總結(jié)并檢驗(yàn)了4 種常用傳熱關(guān)聯(lián)式對本文數(shù)據(jù)的適用性，最后提出并檢驗(yàn)了新的傳熱關(guān)聯(lián)式。

1 物理模型和控制方程

圖2為加熱圓管的幾何模型。圖中豎直圓管外徑為8 mm，內(nèi)徑為5 mm，加熱段長660 mm，進(jìn)出口段長110 mm。在圓管外壁的一側(cè)施加熱流，由于管壁周向和軸向?qū)岬拇嬖冢麄€(gè)圓管的熱流密度是非均勻的。

圖2 加熱圓管幾何模型Fig.2 Geometry model of the heated circular pipe

笛卡爾坐標(biāo)系下可壓縮流體的連續(xù)性方程為

式中：ρ為密度，kg/m3；ui為速度分量，m/s。

動量方程為

式中：gi為重力加速度分量，m/s2；為黏性系數(shù)，kg/(m·s)；uj、uk為速度分量，m/s；u′為速度脈動量，m/s。

能量方程為

式中：λeff為有效導(dǎo)熱系數(shù)，W/(m·K)；h為焓，J/kg；cp為比熱容，J/(kg·K)。

湍動能方程為

式中：μt為湍流黏性系數(shù)，kg/(m·s)；k為湍動能，J/kg；σk、Gk、Yk為模型系數(shù)。

耗散率方程為

式中：ω為耗散率；σω、Gω、Yω、Dω為模型系數(shù)。

式(1)—(5)中，上標(biāo)“—”表示時(shí)間平均項(xiàng)，上標(biāo)“～”表示Favre 平均項(xiàng)，各項(xiàng)具體的表達(dá)式可參考文獻(xiàn)[16]。

2 邊界條件和網(wǎng)格無關(guān)性檢驗(yàn)

在ICEM CFD軟件中劃分六面體結(jié)構(gòu)網(wǎng)格并加密流固耦合面的網(wǎng)格。流體域第1 層網(wǎng)格高度設(shè)為0.000 2，以保證第1 層網(wǎng)格質(zhì)心到內(nèi)壁面的無量綱距離y+小于1。數(shù)值模擬借助ANSYS Fluent 軟件進(jìn)行，固體域材料為316L 不銹鋼，流體域材料選取ANSYS Fluent 軟件內(nèi)嵌套的NIST 真實(shí)氣體物性庫中的實(shí)際CO2，開啟Z方向的重力加速度。入口邊界條件選取可壓縮流體常用的質(zhì)量入口，出口邊界條件選取壓力出口。湍流模型選取SSTk-ω兩方程模型，采用壓力速度耦合的SIMPLIC 算法，動量、能量、湍動能和耗散率方程的離散均采用二階迎風(fēng)格式，松弛因子保持默認(rèn)。連續(xù)性方程、動量方程、能量方程、湍動能和耗散率方程的收斂殘差設(shè)置為1.0×10-6，當(dāng)殘差曲線和進(jìn)出口平均溫度不發(fā)生明顯變化，且進(jìn)出口質(zhì)量流量的誤差小于1%時(shí)，可認(rèn)為計(jì)算收斂。

選取質(zhì)量流速G=100 kg/(m2·s)、熱流密度qwi=50 kW/m2、入口壓力p=7.6 MPa、入口溫度Tin=295.15 K 的工況，分別采用4 套網(wǎng)格計(jì)算所給幾何模型的加熱段壓降，結(jié)果見表1。

表1 網(wǎng)格數(shù)與加熱段壓降Tab.1 The grids number and the pressure drop of heating section

由表1可知，加熱段壓降隨網(wǎng)格數(shù)量增大而增大，但當(dāng)網(wǎng)格數(shù)增大到848 194 時(shí)，壓降的變化很小。考慮到計(jì)算成本和求解精度，本文選取網(wǎng)格數(shù)848 194 進(jìn)行數(shù)值模擬。

3 模型驗(yàn)證

為驗(yàn)證數(shù)值方法的可靠性，需進(jìn)行模型驗(yàn)證。因?yàn)楝F(xiàn)存文獻(xiàn)中鮮有非均勻熱流條件下超臨界CO2的傳熱試驗(yàn)，所以選取Kim 等人[12]在均勻熱流條件下進(jìn)行的試驗(yàn)中的1 組數(shù)據(jù)(G=238 kg/(m2·s)，qwi=52 kW/m2，p=7.621 MPa)進(jìn)行模型驗(yàn)證，得到內(nèi)壁溫和主流溫度沿管長的分布如圖3所示。結(jié)果表明，試驗(yàn)和模擬所得內(nèi)壁溫和主流溫度的誤差均在3%以內(nèi)，認(rèn)為模型可靠。

圖3 內(nèi)壁溫和主流溫度沿管長的分布Fig.3 The distribution of inner wall temperature and bulk temperature along the tube length

4 結(jié)果與討論

4.1 均勻與非均勻加熱

與均勻加熱的情況不同，非均勻加熱時(shí)，管壁加熱面和未加熱面的局部熱流密度不同，壁溫和對流傳熱系數(shù)也不同。圖4給出了p=7.6 MPa、G=100 kg/(m2·s)、qwi=50 kW/m2工況下非均勻加熱和均勻加熱的對流傳熱系數(shù)和內(nèi)壁溫隨流體主流溫度的變化曲線。

為便于比較，對于均勻加熱，熱流密度即實(shí)際施加于加熱面的熱流密度；對于非均勻加熱，熱流密度為整個(gè)圓周面的平均熱流密度，實(shí)際每處的熱流密度并不相同。其中，加熱面對流傳熱系數(shù)定義如下：

未加熱面對流傳熱系數(shù)定義如下：

總對流傳熱系數(shù)定義如下：

式中，hh、huh、h為加熱面對流傳熱系數(shù)、未加熱面對流傳熱系數(shù)、總對流傳熱系數(shù)，kW/(m2·K)；Twi、Tb為內(nèi)壁面溫度和主流溫度，K。

由圖4a)可知，2 種加熱情況對應(yīng)的對流傳熱系數(shù)變化趨勢相同：擬臨界點(diǎn)前，都隨比定壓熱容的增大而增大；擬臨界點(diǎn)后，隨比定壓熱容的減小而減小。非均勻加熱時(shí)，未加熱面的對流傳熱系數(shù)最大，加熱面的對流傳熱系數(shù)最??；均勻加熱時(shí)，加熱面的對流傳熱系數(shù)則介于兩者之間。

圖4b)的內(nèi)壁溫變化反映著圖4a)中對流傳熱系數(shù)的變化：非均勻加熱時(shí)，加熱面的局部熱流密度最大，對流傳熱系數(shù)最小，內(nèi)壁溫最高。因此，較之均勻加熱，非均勻加熱受熱面壁溫更容易超溫。

圖4 對流傳熱系數(shù)和內(nèi)壁溫變化曲線Fig.4 Change curves of the convective heat transfer coefficient and inner wall temperature

4.2 非均勻加熱熱流密度的影響

圖5為G=100 kg/(m2·s)、p=7.6 MPa 工況下，熱流密度對非均勻傳熱的影響。

由圖5a)可以看出，3 種熱流密度條件下，總對流傳熱系數(shù)隨主流溫度的升高先增大后減小，在擬臨界點(diǎn)附近出現(xiàn)峰值。擬臨界點(diǎn)前，熱流密度越小，對流傳熱系數(shù)越大；擬臨界點(diǎn)后，3 種熱流密度下的對流傳熱系數(shù)相差甚微。這是因?yàn)椋簲M臨界點(diǎn)前流體因溫度升高比定壓熱容不斷增大，傳熱能力增強(qiáng)；而以上3 種工況對應(yīng)的Gr/Re2.7均大于文獻(xiàn)[17]所提到的臨界值10?5，傳熱受物性和浮力2 種因素的影響；由圖5b)可知，熱流密度越大，浮力越大，浮力效應(yīng)對傳熱的抑制越嚴(yán)重。因此：28 kW/m2熱流密度下的對流傳熱系數(shù)最大，39 kW/m2次之，50 kW/m2最小；擬臨界點(diǎn)附近，Gr/Re2.7驟減，浮力對傳熱的抑制驟減，對流傳熱系數(shù)達(dá)到峰值；擬臨界點(diǎn)后，流體比定壓熱容減小，傳熱系數(shù)減?。贿h(yuǎn)離擬臨界點(diǎn)的位置，流體的比定壓熱容隨主流溫度變化不大，3 種工況對應(yīng)的Gr/Re2.7也隨主流溫度變化不大，且相差甚微。

圖5 熱流密度對非均勻傳熱的影響Fig.5 Effects of heat flux on non-uniform heat transfer

圖5c)給出了3 種熱流密度條件下管壁加熱面和未加熱面對流傳熱系數(shù)的變化?？傮w趨勢和圖5a)相同，但擬臨界點(diǎn)前，未加熱面的對流傳熱系數(shù)大于加熱面的對流傳熱系數(shù)，擬臨界點(diǎn)后，兩者相差較小。這是因?yàn)?，擬臨界點(diǎn)前，加熱面內(nèi)壁溫和主流溫度相差較大而未加熱面內(nèi)壁溫和主流溫度相差較小，隨著管壁周向?qū)岷土黧w對流傳熱的不斷進(jìn)行，未加熱和加熱面的內(nèi)壁面溫差逐漸減小，因此，擬臨界點(diǎn)后對流傳熱系數(shù)的差值也逐漸減小。

圖5d)給出了3 種熱流密度下內(nèi)壁溫隨主流溫度的變化。熱流密度越大，內(nèi)壁溫越高。擬臨界點(diǎn)前，內(nèi)壁溫變化平緩；擬臨界點(diǎn)后，內(nèi)壁溫飛升；遠(yuǎn)離擬臨界點(diǎn)的位置，溫度變化又變得平緩。這與圖5a)中對流傳熱系數(shù)的變化一致。

4.3 非均勻加熱入口壓力的影響

選取p=7.6、8.0、8.6 MPa 共3 組工況（G=100 kg/(m2·s)，qwi=50 kW/m2），分析入口壓力對非均勻傳熱的影響，如圖6所示。

由圖6a)可以看出，不同于低熱流密度下壓力對對流傳熱系數(shù)的影響，在本文所給的高熱流密度下，入口壓力越大，對流傳熱系數(shù)越大。結(jié)合圖6b)可知，在當(dāng)前熱流密度下，Gr/Re2.7遠(yuǎn)大于文獻(xiàn)[17]提到的Gr/Re2.7臨界值10?5，浮力對傳熱的抑制作用十分明顯。擬臨界點(diǎn)前，流體因不斷吸熱，溫度升高，比定壓熱容逐漸增大，傳熱能力增強(qiáng)；擬臨界點(diǎn)處，流體比定壓熱容達(dá)到峰值，對流傳熱系數(shù)也達(dá)到峰值；擬臨界點(diǎn)后，流體比定壓熱容驟減，對流傳熱系數(shù)減?。贿h(yuǎn)離擬臨界點(diǎn)的位置，流體比定壓熱容變化緩慢，對流傳熱系數(shù)隨之變化緩慢。7.6 MPa 入口壓力下流體的物性隨主流溫度變化的劇烈程度大于8.0 MPa 和8.6 MPa 壓力下物性變化程度，管內(nèi)近壁區(qū)流體與主流區(qū)流體密度差更大，導(dǎo)致浮力更大，對湍流傳熱的抑制作用更強(qiáng)。所以，入口壓力較大的一組，對流傳熱系數(shù)更大。熱流密度較低時(shí)，物性變化引起的浮力效應(yīng)較弱，甚至可忽略其對湍流傳熱的抑制作用，此時(shí)，比定壓熱容的變化對湍流傳熱起主導(dǎo)作用。因此，在低熱流密度下，壓力越小，對流傳熱系數(shù)越大。

圖6 入口壓力對非均勻傳熱的影響Fig.6 Effects of inlet pressure on non-uniform heat transfer

圖6c)對比了3 種入口壓力下，管壁加熱面和未加熱面對流傳熱系數(shù)的變化?？傏厔萃瑘D6a)。

圖6d)給出了3 種入口壓力下內(nèi)壁溫隨主流溫度的變化曲線。由圖6d)可知，8.6 MPa 入口壓力對應(yīng)的內(nèi)壁溫較低，但與其他2 個(gè)條件下的內(nèi)壁溫相差不大。

4.4 熱流密度和入口壓力對加速效應(yīng)的影響

除CO2物性變化和浮力影響外，加速效應(yīng)也是一個(gè)影響超臨界CO2在管內(nèi)傳熱流動的重要因素。由前文可知，因流動過程的壓降遠(yuǎn)小于入口壓力，可忽略壓降引起的加速效應(yīng)。但是，隨著管內(nèi)流體溫度的升高、密度的降低和熱膨脹系數(shù)的變化，軸向的密度差也會使流體產(chǎn)生加速效應(yīng)，進(jìn)而使流體速度梯度發(fā)生變化，切應(yīng)力隨之變化，影響湍流的生成、擴(kuò)散和傳熱能力，嚴(yán)重時(shí)易引發(fā)湍流的層流化。McEligot 等人[18]研究了湍流向?qū)恿鬓D(zhuǎn)變過程中的傳熱特性，提出了無量綱準(zhǔn)則數(shù)Kv的表達(dá)式（式(9)），以判別因密度差引起的加速效應(yīng)對管內(nèi)超臨界流體傳熱的影響。

當(dāng)Kv＜3×10?6時(shí)，流動狀態(tài)保持為湍流；當(dāng)Kv≥3×10?6時(shí)，湍流減弱，傳熱惡化。

圖7給出了熱流密度和入口壓力對加速效應(yīng)的影響。

圖7 熱流密度和入口壓力對加速效應(yīng)的影響Fig.7 Effects of heat flux and inlet pressure on acceleration effect

由圖7a)可知，Kv隨主流溫度的升高先減小后增大，擬臨界點(diǎn)為最小值。正如式(9)所示：隨著溫度的升高，流體的比定壓熱容增大，在擬臨界點(diǎn)處達(dá)到最大值，此時(shí)Kv最??；擬臨界點(diǎn)后，隨著流體比定壓熱容減小，Kv逐漸增大。而熱流密度越大，流體軸向密度場分布越不均勻，加速效應(yīng)就越強(qiáng)。

由圖7b)可知，Kv的變化與入口壓力對浮力效應(yīng)的影響類似：壓力越小，Kv越大，加速效應(yīng)就越強(qiáng)。原因是7.6 MPa 下CO2的物性變化較8.0 MPa和8.6 MPa 更為劇烈。比較熱流密度和入口壓力的影響可知，熱流密度對加速效應(yīng)的影響更大，這一點(diǎn)從Kv的表達(dá)式也可以看出。本文所有工況得到的Kv都遠(yuǎn)小于臨界值，因此，加速效應(yīng)的影響可忽略。

4.5 模擬結(jié)果與現(xiàn)有關(guān)聯(lián)式的比較

對于亞臨界單相流體，經(jīng)典的Dittus-Boelter 關(guān)聯(lián)式[19]對傳熱有著較好的預(yù)測，而超臨界流體因其物性在擬臨界點(diǎn)劇烈變化，以及引起的浮力和加速效應(yīng)，使得Dittus-Boelter 關(guān)聯(lián)式所預(yù)測的傳熱特性和試驗(yàn)結(jié)果存在很大偏差。針對這一問題，許多學(xué)者在各自試驗(yàn)的基礎(chǔ)上，提出了許多超臨界CO2在不同條件下的傳熱關(guān)聯(lián)式。這些關(guān)聯(lián)式大多把主流區(qū)流體和近壁區(qū)流體物性的比值作為修正項(xiàng)引入Dittus-Boelter 關(guān)聯(lián)式。表2給出了4 種常見的傳熱關(guān)聯(lián)式。

表2 4 種常見的傳熱關(guān)聯(lián)式Tab.2 Four kinds of heat transfer correlations

圖8為G=100 kg/(m2·s)、qwi=50 kW/m2、p=7.6 MPa 條件下由表2中傳熱關(guān)聯(lián)式計(jì)算所得Nu和本文模擬值對比。由圖8可見，不同傳熱關(guān)聯(lián)式所得Nu相差很大，擬臨界點(diǎn)附近更為明顯。并且，4 種關(guān)聯(lián)式預(yù)測的Nu與本文模擬值相差很大，主要原因在于本文所給的計(jì)算模型為非均勻熱流密度下的豎直圓管內(nèi)超臨界CO2傳熱，而4 種關(guān)聯(lián)式是針對均勻熱流密度得到的，因此僅適用于各自的試驗(yàn)工況。

圖8 4 種常見傳熱關(guān)聯(lián)式計(jì)算Nu 與模擬值對比Fig.8 The Nu calculated by the above four common correlations and the simulated values

為較為準(zhǔn)確地預(yù)測依據(jù)本文計(jì)算模型所得Nu，在Dittus-Boelter 關(guān)聯(lián)式的基礎(chǔ)上，結(jié)合浮力效應(yīng)和物性變化，提出以下關(guān)聯(lián)式：

式(10)中，浮力效應(yīng)由近壁區(qū)流體和主流區(qū)流體的密度差產(chǎn)生，因此，用近壁區(qū)流體和主流區(qū)流體的密度之比作為浮力修正項(xiàng)。物性變化中，比定壓熱容的變化最明顯，采用平均比定壓熱容與主流區(qū)流體比定壓熱容之比作為物性修正項(xiàng)。由4.4 節(jié)可知，加速效應(yīng)的影響可忽略不計(jì)。利用多元線性回歸的方法求得以下傳熱關(guān)聯(lián)式：

任意選取3 組工況，比較新關(guān)聯(lián)式計(jì)算的Nu和本文模擬值，結(jié)果如圖9所示。由圖9可見，新傳熱關(guān)聯(lián)式與模擬值吻合較好。但應(yīng)注意，雖然超臨界CO2的傳熱規(guī)律具有一定的相似性，但CO2物性的劇烈非線性變化及其帶來的浮力效應(yīng)和加速效應(yīng)給研究帶來了很大的困難。目前，還不存在一個(gè)普遍適用且準(zhǔn)確的傳熱關(guān)聯(lián)式。因此，新關(guān)聯(lián)式也只適用于本文中的數(shù)據(jù)。

圖9 新關(guān)聯(lián)式計(jì)算Nu 與模擬值對比Fig.9 The Nu calculated by the new proposed correlation and the simulated values

5 結(jié) 論

1）相同加熱量條件下，與均勻加熱相比，非均勻加熱局部熱流密度大，對流傳熱系數(shù)小，內(nèi)壁溫更高。

2）在壓力和質(zhì)量流速一定的情況下，增大熱流密度，對流傳熱系數(shù)減小，浮力和加速效應(yīng)增強(qiáng)；在高熱流密度下，增大入口壓力，對流傳熱系數(shù)增大，浮力和加速效應(yīng)減弱。

3）對于非均勻加熱圓管，未加熱面的對流傳熱系數(shù)大于加熱面的對流傳熱系數(shù)，隨著周向?qū)岷蛯α鞯淖饔?，兩者之差逐漸減小。

4）新的傳熱關(guān)聯(lián)式能較好地預(yù)測文中超臨界CO2的傳熱規(guī)律。