超臨界壓力CO2在水平圓管內流動傳熱數值分析*

閆晨帥 徐進良

(華北電力大學低品位能源多相流與傳熱北京市重點實驗室, 北京 102206)

采用SST k-ω低雷諾數湍流模型對加熱條件下超臨界壓力CO2在內徑di = 22.14 mm, 加熱長度Lh =2440 mm水平圓管內三維穩態流動與傳熱特性進行了數值計算.通過超臨界CO2在水平圓管內的流動傳熱實驗數據驗證了數值模型的可靠性和準確性.首先, 研究了超臨界壓力CO2在水平圓管內的流動傳熱特點,基于超臨界CO2在類臨界溫度Tpc處發生類液?類氣“相變”的假設, 揭示了水平圓管頂母線和底母線區域不同的流動傳熱行為.然后, 分析了熱流密度qw和質量流速G對水平圓管內超臨界壓力CO2流動換熱的影響,通過獲取流體域內的物性分布、速度分布和湍流分布等詳細信息, 重點解釋了不同熱流密度qw和質量流速G下頂母線內壁溫度Tw,i分布產生差異的傳熱機理, 分析結果確定了類氣膜厚度δ、類氣膜性質、軸向速度u和湍動能k是影響頂母線壁溫分布差異的主要因素.研究結果可以為超臨界壓力CO2換熱裝置的優化設計和安全運行提供理論指導.

1 引 言

作為一種新型動力循環技術, 超臨界CO2(supercritical carbon dioxide, S?CO2)布雷頓循環在核能、太陽能、化石能源和地熱能等領域具有廣闊的應用前景[1?4].S?CO2流動傳熱是 S?CO2布雷頓循環的關鍵技術之一, 充分地了解并掌握S?CO2傳熱機理對于換熱裝置設計和安全穩定運行是至關重要的.自20世紀50年代以來, 各國學者對S?CO2傳熱進行了大量的實驗和數值分析, 但主要圍繞S?CO2在垂直圓管內向上和向下流動傳熱進行研究[5?11].然而, 相對于垂直管內的流動傳熱, 加熱條件下S?CO2在水平管內垂直方向會受到重力和浮升力的雙重作用, 發生流動分層現象,這將引起圓管頂部和底部不同的傳熱行為, 使得S?CO2在水平管內的流動傳熱特性更為復雜.

Adebiyi和 Hall[12]對水平圓管內 S?CO2流動傳熱進行了實驗研究, 通過觀察管壁溫度變化趨勢發現在圓管底部的傳熱能力較強, 而頂部的傳熱能力有所弱化, 這與Pidaparti等[13]和Chu與Laurien[14]的研究結果是一致的.Wang等[15]實驗研究了S?CO2在水平圓管內的對流傳熱特性, 通過采用浮升力準則數和流動加速準則數Kv對實驗數據進行分析認為, 在高熱流密度和低質量流速下,傳熱惡化與浮升力效應影響有關, 而流動加速效應可以忽略不計.Tanimizu和 Sadr[16]采用 BuC,BuJ和BuP三個浮升力準則參數預測了水平圓管內浮升力效應對S?CO2流動換熱的影響, 雖然研究結果顯示浮升力準則數的大小和變化趨勢與傳熱變化關聯性較小, 但Tanimizu和Sadr仍認為浮升力效應存在并且對流動傳熱具有重要影響.Kim等[17]實驗研究也發現浮升力準則數(BuC和BuJ)的變化趨勢與傳熱變化的關聯性并不確切,雖然Kim等同樣認為浮升力效應存在, 但是由于Kim等引入了流動加速準則數q+并發現其與傳熱變化關聯性較好, 認為流動加速與傳熱變化密切相關, 浮升力效應可能對于流動傳熱影響極小.

綜上所述, 現有研究主要通過實驗方法分析加熱條件下S?CO2在水平圓管內的流動傳熱特性,而且各學者對S?CO2在水平圓管內流動傳熱發生變化原因的意見并不一致.此外, 由于實驗方法僅能測量表觀數據, 不能獲得流體域內部詳細的參數分布, 因此采用數值模擬的方法, 對S?CO2在水平圓管內流動傳熱機理進行進一步研究是十分必要的.鑒于此, 本文主要包括以下研究內容.首先, 通過S?CO2在水平圓管內流動傳熱實驗數據驗證數值模型的可靠性和準確性, 并建立拓展計算工況,分析水平圓管內S?CO2在恒定熱流密度加熱條件下的流動傳熱特點.然后, 分析熱流密度qw和質量流速G對S?CO2在水平圓管內流動傳熱的影響,通過獲得流體域詳細的溫度分布、物性分布和流場分布等信息, 揭示不同熱流密度qw和質量流速G條件下S?CO2在水平圓管內流動傳熱機理及管壁溫度Tw,i產生差異的原因.研究結果可以為S?CO2布雷頓循環換熱裝置設計和安全穩定運行提供理論參考.

2 物理模型與數值方法

2.1 幾何模型

水平圓管的計算幾何模型見圖1.計算域包括入口絕熱段和加熱段兩部分, 入口絕熱段長度Lad=1220 mm, 以確保加熱段入口處的工質已處于充分發展狀態; 加熱段長度Lh= 2440 mm, 是研究S?CO2流動傳熱特性的主要部分, 圓管內徑di=22.14 mm.計算過程中忽略圓管壁厚, 圓管周向角度定義為 φ, φ = 0°的位置定義為頂母線, φ =180°的位置定義為底母線.

2.2 控制方程

圖1 幾何模型示意圖Fig.1.Schematic diagram of the geometric model.

本文數值計算模型包括連續性方程、動量方程和能量方程, 并采用以下表示形式對上述方程進行描述.

連續性方程:

動量方程:

(2)式中,

能量方程:

(1)—(3)式中, u, ρ, μ和 h 分別表示工質的速度、密度、黏度和焓值.

選擇合適的湍流模型對S?CO2流動傳熱數值計算結果的準確性至關重要.對于超臨界流體流動換熱計算, Wang等[18]的計算結果顯示, 采用SST k?ω低雷諾數湍流模型能夠得到比其它湍流模型更準確的計算結果.因此, 本文采用SST k?ω低雷諾數湍流模型進行數值計算.SST k?ω湍流模型簡述如下.

k輸運方程:

ω輸運方程:

(4)和(5)式中的函數項和常數項詳見文獻[19].

2.3 網格劃分及計算方法

計算域網格由ANSYS ICEM軟件生成, 采用結構化六面體網格, 圓管橫截面采用O型網格劃分.由于S?CO2物性變化劇烈, 在圓管近壁區域流動換熱情況極其復雜, 因此在流體域近壁面處對網格進行局部加密, 并始終保持第一層網格無量綱高度y+＜ 1.為檢驗計算網格獨立性, 圖2比較了不同網格總節點數時, 水平圓管頂母線內壁溫度Tw,i沿圓管軸向z/d的分布情況.由圖2可以看出, 當網格總節點數大于269.7萬時, 計算結果隨節點總數的增加已無明顯偏差.因此, 計算網格總節點數最終確定為269.7萬個.

圖2 網格無關性驗證Fig.2.Verification for grid independence.

采用ANSYS Fluent 15.0計算水平圓管內S?CO2三維穩態流動傳熱的情況.入口設置為質量流速入口, 出口設置為壓力出口, 入口絕熱段管壁設置為絕熱邊界, 加熱段管壁設置為均勻熱流密度且為無滑移剪切力邊界條件.控制方程采用有限體積法進行離散, 采用二階迎風差分方法以提高計算精度, 壓力?速度耦合方程應用SIMPLEC算法求解.采用NIST實際氣體模型以準確反映超臨界流體物性變化對流動傳熱的影響.

3 計算結果及分析

3.1 S-CO2熱物性

圖3為S?CO2在不同運行壓力下物性參數(定壓比熱cp、密度ρ、黏度μ和導熱系數λ)隨溫度的變化曲線.由圖3可以看出, 不同壓力下S?CO2物性參數隨溫度的變化趨勢是相似的, 對于給定超臨界壓力, cp總是存在一個極大值, cp極大值對應溫度定義為類臨界溫度Tpc, 隨著溫度的升高,密度、導熱系數和黏度會大大降低, 特別是在類臨界溫度Tpc附近.Holman等[20]對垂直加熱管環內超臨界R12湍流傳熱進行了可視化實驗研究, 在近壁區域發現了類似于蒸汽的跡象, Holman等認為產生這一現象的原因可能是由于近壁區劇烈的密度梯度導致的, 即使流體在超臨界狀態下不存在亞臨界狀態時的氣?液界面, 通過“類沸騰”命名這一行為也是合理的.Banuti[21]認為, 在Tpc附近超臨界流體熱物性的劇烈變化相似于典型的亞臨界蒸發行為, 并基于“類沸騰”的概念對超臨界流體在Tpc附近的變化行為進行了詳細的理論說明.因此, 為了便于分析, 本文將溫度低于Tpc的工質狀態定義為類液相, 溫度高于Tpc的工質狀態定義為類氣相, Tpc視為類氣?類液的“相變”溫度.

需要說明的是, 與亞臨界條件下相變過程以及相界面分界氣?液相的現象不同, S?CO2在Tpc附近“類相變”過程是連續的, 而且超臨界狀態下CO2的表面張力σ為零[22].對于傳統的亞臨界氣?液兩相流傳熱數值計算問題, 引入的多相流計算模型需考慮相界面的表面張力效應[23?24].顯然, 亞臨界氣?液兩相流傳熱計算模型并不適用于S?CO2.因此, S?CO2流動傳熱數值模擬屬于單相流體流動傳熱的計算范疇, 引入兩相流計算模型并跟蹤“相界面”運動對傳熱的影響是難以實現的.

水平圓管內主流工質溫度Tb定義:

對流傳熱系數h定義:

主流焓值ib定義:

(6)—(8)式中, ρ 為當地流體密度, kg·m—3; u 為當地流體速度, m·s—1; T為當地流體溫度, K; A為水平圓管流體域橫截面面積, m2; qw為壁面熱流密度, kW·m—2; i為流體當地焓值, kJ·kg—1.

3.2 計算模型驗證

圖3 不同壓力下, S?CO2熱物性曲線 (a) 定壓比熱cp; (b) 密度ρ; (c) 導熱系數λ; (d) 黏度μFig.3.Thermophysical properties of S?CO2 under different pressures: (a) Specific heat at constant pressure cp; (b) density ρ;(c) thermal conductivity λ; (d) viscosity μ.

圖4 數值預測壁溫與Adebiyi和Hall[12]實驗數據比較Fig.4.Comparison of wall temperature predicted by simulation against experimental data by Adebiyi and Hall[12].

為了驗證本文所采用計算方法和湍流模型的可靠性與準確性, 本文基于Adebiyi和Hall[12]的實驗工況運行參數對水平圓管內S?CO2流動傳熱進行了數值模擬.水平管內徑為22.14 mm, 均勻加熱.圖4對模擬計算值與實驗測量結果進行了比較, 受計算模型理想假設、幾何模型簡化以及實驗測量誤差等因素的影響, 計算值與實驗值之間存在一定的偏差, 但是兩者沿軸向變化趨勢總體上保持一致.圖4(a)和圖4(b)中計算最大偏差值分別是1.0%和0.3%, 說明計算結果與實驗數據吻合良好.

3.3 傳熱特性

本節以計算工況P = 8.0 MPa, G = 600 kg·m—2·s—1, qw= 40 kW·m—2為例, 主要分析了加熱條件下水平圓管內S?CO2典型的流動傳熱特點.圖5給出了S?CO2在水平圓管內流動換熱過程中頂母線與底母線軸向內壁溫度Tw,i和對流換熱系數h、周向類氣膜厚度δ分布的典型結果.其中,圖5(b)中類氣膜厚度δ周向分布曲線對應圖5(a)中主流焓值ib= 244.8 kJ·kg—1處的圓管橫截面.需要說明的是, 為了便于分析, 以圓管橫截面內工質當地溫度T = Tpc所處位置表征穩態時類液相與類氣相的“界面”, 氣?液“界面”至圓管內壁的徑向距離定義為類氣膜厚度δ.

由圖5(a)可以看出, 在加熱條件下, 水平圓管加熱段入口區域頂母線壁溫Tw,i劇烈升高, 而底母線壁溫變化則比較緩和.隨著加熱的不斷進行,頂母線壁溫變化也逐漸趨于穩定, 但是頂母線內壁溫度Tw,i明顯高于底母線, 并形成較大的壁溫差值.此外, 不難發現, 頂母線與底母線換熱系數h的相對大小則與壁溫分布呈相反的變化趨勢.由圖5(b)可以看出, 水平圓管內類氣膜厚度δ沿周向是非均勻分布的, 頂母線處類氣膜厚度δ最大,底母線處類氣膜厚度δ最小甚至為零.

圖5 水平圓管內S?CO2傳熱特性典型曲線 (a) 頂母線與底母線內壁溫Tw,i和傳熱系數h軸向分布; (b) 類氣膜厚度δ周向分布Fig.5.Typical curves of heat transfer characteristics for S?CO2 in horizontal tube: (a) Axial distribution of inner wall temperature and heat transfer coefficient at top generatrix and bottom generatrix; (b) circumferential distribution of vapor?like film thickness.

圖6 特征橫截面內矢量速度分布Fig.6.Vector velocity distribution in the characteristic cross?sections.

圖7 特征橫截面內溫度分布Fig.7.Temperature distribution in the characteristic cross?sections.

圖6和圖7分別顯示了水平圓管特征橫截面內S?CO2矢量速度與溫度分布情況, 特征截面A,B和C的選取詳見圖5(a).其中, ib,A= 241.1 kJ·kg—1,ib,B= 244.8 kJ·kg—1和 ib,C= 251.8 kJ·kg—1.由圖6可以看出, 隨著加熱的不斷進行, 圓管近壁區工質密度減小甚至發生類相變, 近壁區域密度低、傳熱能力差的工質受浮升力的影響沿管壁向頂母線區域流動, 圓管核心區域密度高、傳熱能力強的工質受重力影響向底母線區域流動, 橫截面內由此形成二次流.此外, 對比圖6(a)—(c)可以看出, 在傳熱階段A → B, 二次流速度整體有所增加, 表明近壁區類氣工質將迅速匯聚至頂母線區域形成高溫度、低密度、傳熱能力弱的類氣膜, 嚴重抑制了頂母線區域管壁熱量向主流工質的傳遞; 在傳熱階段B → C, 二次流速度逐漸減小, 主要原因是在加熱段下游區域, 近壁區工質密度與主流核心工質密度差異逐漸減小.二次流的形成使得圓管橫截面內發生流體溫度和類氣?類液分層現象, 而且橫截面內工質的“熱力學干度”沿著工質流動方向呈逐漸增加的趨勢, 詳見圖7.不難發現, 以上現象與亞臨界條件下水平圓管內氣?液兩相流動出現的現象是相似的[25].

圖8顯示了水平圓管特征截面內軸向速度u沿徑向的分布情況.由圖8可以看出, S?CO2在水平圓管內頂母線至底母線的軸向速度u相對于圓管中心呈不對稱分布, 頂母線區域軸向速度u低于底母線, 因此頂母線處工質沿軸向轉移來自于管壁熱量的能力比底母線是更差的.此外, 一個有趣的現象是, 隨著主流焓值ib的增加, 在發生壁溫飛升的傳熱階段A → B, 頂母線附近工質軸向速度u總體逐漸減小, 而在傳熱階段B → C, 軸向速度u逐漸增大; 然而底母線附近工質軸向速度u則隨主流焓值ib整體總是逐漸增加的.因此, 頂母線較厚的類氣膜厚度δ和較低的軸向速度u是頂母線壁溫Tw,i高于底母線的重要原因.

圖8 水平圓管內S?CO2軸向速度沿徑向分布Fig.8.Radial distribution of axial velocity for S?CO2 in horizontal tube.

3.4 熱流密度對傳熱的影響

在保持壓力P = 9 MPa和質量流速G =500 kg·m—2·s—1不變條件下, 分別采用 qw= 30,50和70 kW·m—2評價熱流密度對水平圓管內S?CO2流動換熱的影響.需要說明的是, 3.3節中的計算結果顯示頂母線壁溫變化特征比底母線更為顯著, 因此本文主要通過分析頂母線附近的物性分布和流場分布的特點, 探究熱流密度和質量流速對傳熱流動的影響及原因.此外, 為了定量地說明圓管橫截面內二次流對傳熱的影響, 本文引入二次流動能參數Ke以表征二次流強度的大小[26?28], 其表達式為

(9)式中, v和w分別表示圓管橫截面內沿x與y方向的速度分量.

圖9展示了不同熱流密度qw條件下頂母線壁溫Tw,i及二次流強度Ke隨主流焓值ib的分布曲線.由圖9(a)可以看出, 隨著熱流密度qw的增加,頂母線壁溫Tw,i明顯升高, 特別是加熱段入口區域的壁溫增加幅度逐漸劇烈, 但是在加熱段下游區域, 壁溫均體現出比較緩和的變化趨勢.圖9(b)顯示熱流密度qw越大, 二次流強度Ke總體水平越高, 說明橫截面內將有更多的低密度工質沿管壁流動至頂母線附近.而且, 不同熱流密度條件下加熱段入口區域Ke均先顯著升高, 然后逐漸減小并呈較小幅度的變化水平.綜合圖9(a)和圖9(b)可以看出, 圓管橫截面二次流強度Ke大小對管壁溫度變化具有重要影響.

圖9 熱流密度qw對壁溫和二次流強度的影響Fig.9.Influences of heat flux on wall temperature and second flow intensity.

圖10 不同熱流密度時, 水平圓管特征截面內熱物性、軸向速度和湍動能沿徑向分布曲線 (a) 定壓比熱cp; (b) 導熱系數λ;(c) 軸向速度u; (d) 湍動能kFig.10.Radial distribution curves of thermophyical properties in the characteristic cross?sections under different heat flux: (a) Spe?cific heat at constant pressure cp; (b) thermal conductivity λ; (c) axial velocity u; (d) turbulent kinetic energy k.

選取圖9(a)中熱流密度 qw= 30 kW·m—2和qw= 70 kW·m—2工況下特征截面 A 和 A' (ib,A=ib,A'= 236.8 kJ·kg—1), B 和 B' (ib,B= ib,B'=242.3 kJ·kg—1)分析 S?CO2在不同熱流密度加熱條件下水平圓管內流動換熱差異的原因.圖10顯示了不同熱流密度qw下, 特征截面內熱物性(包括定壓比熱cp和導熱系數λ)、軸向速度u和湍動能k沿徑向分布情況.由圖10(a)和圖10(b)可以看出, 熱流密度(qw= 30 kW·m—2)較低時, 頂母線區域仍為密度較大的類液工質, 并未形成類氣膜, 因此近壁區定壓比熱cp和導熱系數λ處于較高水平.然而, 當 qw= 70 kW·m—2時, 頂母線區域聚集了大量低密度類氣工質并形成類氣膜, 并且隨著加熱的進行, 類氣膜逐漸增厚, 大物性變化區域逐漸遠離管壁, 近壁區域表征類氣膜性質的定壓比熱cp和導熱系數λ均較小, 弱化了近壁區類氣膜的吸熱能力和導熱能力.圖10(c)顯示在等焓特征截面處, 相對于低熱流密度, 高熱流密度時頂母線區域近壁工質在局部熱加速效應作用下軸向速度u略大, 這有利于S?CO2沿軸向及時輸運轉移來自管壁的熱量.由圖10(d)可以看出, 高熱流密度條件下頂母線區域類氣膜內部及主流核心工質的湍動能k整體上總是較小的, 這對于管壁熱量在徑向向主流核心工質傳遞是具有抑制作用的.

綜上可以看出, 在高熱流密度條件下, 除了施加到管壁的熱量增加導致管壁溫度較高這一影響因素外, 頂母線區域類氣膜的形成和增厚以及類氣膜內部及主流核心工質較小的湍動能k將進一步削弱圓管壁面熱量向主流工質的擴散, 這是較高熱流密度條件下加熱段入口區域壁溫增加幅度更劇烈且管壁溫度總體更高的重要原因.

3.5 質量流速對傳熱的影響

圖11 質量流速G對壁溫和二次流強度的影響Fig.11.Effects of mass flux on wall temperature and second flow intensity.

圖12 不同質量流速下, 水平圓管特征截面內熱物性、軸向速度和湍動能沿徑向分布曲線 (a) 定壓比熱cp; (b) 導熱系數λ;(c) 軸向速度u; (d) 湍動能kFig.12.Radial distribution curves of thermophyical properties in the characteristic cross?sections under different mass flux: (a) Spe?cific heat at constant pressure cp; (b) thermal conductivity λ; (c) axial velocity u; (d) turbulent kinetic energy k.

圖11展示了不同質量流速G條件下頂母線壁溫Tw,i及二次流強度Ke隨主流焓值ib的分布曲線.由圖11(a)可以看出, 隨著質量流速G的增加, 頂母線壁溫Tw,i顯著降低, 加熱段入口區域的壁溫增加幅度明顯減小, 說明較高的質量流速可以有效地降低壁溫并抑制加熱段入口處的壁溫飛升.在加熱段下游區域, 壁溫同樣表現出比較緩和的變化趨勢.圖11(b)顯示質量流速G越大, 二次流強度Ke總體水平越低, 而且指定質量流速G條件下二次流強度Ke變化趨勢與圖9(b)相似, 同樣反映出圓管橫截面二次流強度Ke大小與管壁溫度變化的關系.

選取圖11(a)中質量流速 G = 600 kg·m—2·s—1和 G = 1000 kg·m—2·s—1工 況 下 特 征 截 面 C 和C' (ib,C= ib,C'= 238.4 kJ·kg—1), D 和 D' (ib,D=ib,D'= 244.3 kJ·kg—1)分析 S?CO2在不同質量流速G條件下水平圓管內流動換熱差異的原因.圖12顯示了不同質量流速下, 特征截面內熱物性、軸向速度u和湍動能k沿徑向分布情況.由圖12(a)和圖12(b)可以看出, 相對于高質量流速(G =1000 kg·m—2·s—1), 低質量流速 (G = 600 kg·m—2·s—1)條件下, 頂母線區域更早地出現了類氣工質并形成熱阻較大的類氣膜, 而且類氣膜總是較厚的.類氣膜越厚, 大物性變化區域距離管壁越遠, 近壁區類氣膜性質越差.由于類氣膜不能及時有效地將管壁熱量傳遞至主流核心區, 將導致類氣溫度升高, 氣膜性質進一步變差.綜上可以看出, 頂母線區域類氣膜的出現和增厚是低質量流速條件下壁溫更高的主要原因之一.

圖12(c)顯示在等焓特征截面處, 相對于高質量流速, 低質量流速時頂母線區域工質軸向速度u總是較小的, 這不利于S?CO2沿軸向輸運轉移來自管壁的熱量.由圖12(d)可以看出, 相對于高質量流速, 低質量流速時頂母線區域類氣膜內部及主流核心工質的湍動能k整體上總是處于較低水平,這將抑制管壁熱量沿徑向向主流核心工質的傳遞.因此, 較低的軸向速度u和湍動能k水平是低質量流速條件下頂母線壁溫更高的另一個主要原因.

4 結 論

1) SST k?ω低雷諾數湍流模型能夠比較可靠和準確地預測S?CO2在水平圓管內的流動換熱情況.

2) 加熱條件下, S?CO2在水平圓管內的流動傳熱特點與亞臨界壓力下水平圓管內氣?液兩相流動傳熱現象是相似的.熱流密度qw越高或質量流速G越小, 頂母線管壁溫度Tw,i則越高, 而且加熱段入口區域管壁溫度飛升幅度越劇烈.

3) 高熱流密度qw或低質量流速G條件下, 圓管橫截面二次流強度Ke整體總是較大的, 反映出橫截面內將有更多的低密度類氣工質沿管壁流動至頂母線附近, 形成熱阻較大的類氣膜, 同時頂母線區域工質的湍動能k處于較低水平, 兩者綜合效應抑制了管壁熱量沿徑向向主流核心工質傳遞, 是管壁溫度較高的主要影響因素.

4) 低質量流速G條件下, 頂母線附近的軸向速度u整體總是較小的, 這將降低S?CO2沿軸向輸運轉移來自頂母線管壁熱量的能力, 是質量流速G較低時管壁溫度較高的另一個重要影響因素.