超臨界二氧化碳水平管內層流混合對流換熱數(shù)值模擬

楊傳勇 徐進良 王曉東 張 偉

(1新能源電力系統(tǒng)國家重點實驗室 北京 102206)(2能源的安全與清潔利用北京市重點實驗室 北京 102206)(3華北電力大學低品位能源多相流與傳熱北京市重點實驗室 北京 102206)

1 引言

超臨界流體由于其特殊的熱物理性質被廣泛應用于超臨界萃取、火箭和航天器的冷卻、空調和熱泵系統(tǒng)等。由于CO2無毒、不可燃，化學性質非常穩(wěn)定，且臭氧消耗潛能ODP=0、全球變暖潛能GWP=1，對環(huán)境具有友好性，因此開始受到人們更多重視。超臨界流體在微通道中的應用較多，如汽車空調采用跨臨界CO2制冷系統(tǒng)時，其氣體冷卻器要做成直徑小于1 mm的微通道管［1］，于是對超臨界流體在微細管道內基于其物性變化、浮升力等諸多因素對對流換熱的影響做出研究是非常有必要的。

Liao 等［2－3］對超臨界 CO2在水平和豎直微細圓管內換熱做了實驗和數(shù)值模擬，分析了管徑0.5—2.16 mm的管道傳熱系數(shù)和Nu數(shù)隨著主流溫度的變化，得到流體雷諾數(shù)Re高達105時管內的浮升力作用仍不可以忽略。張麗娜［4］、曹侃［5］等也對超臨界CO2在管內換熱進行了數(shù)值研究。因實驗條件的限制對微細管內的二次流研究較少，于是對超臨界CO2管內換熱時浮升力引發(fā)的二次流及其強化換熱的機理仍需要開展更深入的研究。本文通過數(shù)值計算對管內因浮升力引發(fā)的二次流進行定量分析，發(fā)現(xiàn)浮升力對管內換熱具有強化作用，分析了重力水平、壁面溫度、入口Rein等對流動和換熱的影響，可以為設計新型高效緊湊的換熱器提供理論依據(jù)。

2 物理模型和數(shù)學模型

2.1 物理模型

選用的物理模型及坐標系統(tǒng)如圖1所示，管道直徑d=0.5 mm，長L=1 000 mm，壁面采用恒壁溫冷卻。在水平管中，重力引發(fā)的浮升力會引起二次流，所以管道壁面參數(shù)沿周向角度θ的分布是不同的，其中管道的上母線和下母線分別對應周向角度90°和－90°。文中所有工況均采用Fluent軟件計算，模型采用變物性，不同溫度和壓力下超臨界CO2的熱物性參數(shù)通過NIST網(wǎng)站獲得，軟件中CO2的熱物性參數(shù)按照piecewise－liner方法輸入。壓力和速度耦合采用SIMPLEC算法，動量和能量方程采用二階迎風格式的差分方法，采用質量入口和自由出流邊界條件，壁面溫度設為T=Tw，且采用無滑移壁面邊界條件。為了消除網(wǎng)格密度對計算結果的影響，通過對逐漸加密的網(wǎng)格進行計算得到近似的網(wǎng)格無關解。

圖1 物理模型及坐標系統(tǒng)Fig.1 Physical model and coordinate system

2.2 數(shù)學模型

水平管中流體因浮升力的影響會產(chǎn)生分層現(xiàn)象，于是本文采用三維模型計算更接近實際。連續(xù)方程、動量方程、能量方程可參見文獻［1］，動量方程考慮重力加速度g’的變化，其中g=9.81 m/s2，文中所有參數(shù)均取國際單位。

主流的截面平均溫度Tb為:

通過管道壁面的熱流量qw為:

壁面的局部傳熱系數(shù)h為:

壁面的局部努塞爾數(shù)為:

通過引入相對二次流動能將二次流大小定量的表示出來。K定義為:

局部范寧摩擦系數(shù)Cf為:

其中:τw為壁面剪切力，由牛頓粘性定律公式來求:由于Cf的值比較小難以區(qū)分，因此選用Cf×Reb來代替Cf，其中截面平均Reb數(shù)定義為:

式中:u、v、w 分別對應 x、y、z三個方向的速度;d為管道直徑;λ為導熱系數(shù);ρ為流體密度;μ為動力粘度;下標b為取截面平均;A為管道橫截面積;w為壁面;g’為重力加速度;Tb為截面平均溫度;qw為熱流量，W/m2;h為傳熱系數(shù)，W/(m2·K)。

3 數(shù)值模擬結果及分析

3.1 模型驗證

壓力為8 MPa時，超臨界CO2的cp、λ、μ和ρ隨T的變化如圖2所示。8 MPa時超臨界CO2在307.6 K時cp達到峰值，這個峰值點就是超臨界CO2在8 MPa時對應的擬臨界點，此點的溫度即為擬臨界溫度。在擬臨界點附近超臨界CO2的物性變化最為顯著，特別是cp和ρ的值變化很大，而擬臨界點附近以外的區(qū)域物性隨溫度的變化很小。

圖2 8.0 MPa時CO2的熱物性參數(shù)隨溫度變化Fig.2 Thermophysical properties of CO2at 8.0 MPa

為驗證模型的正確性，首先對文獻［1］中的工況進行計算，其中超臨界CO2壓力為8 MPa，質量流量m=1.6 ×10－5kg/s，壁面溫度 Tw=298.15 K，入口溫度Tin=393.15 K，將計算得到的不同位置處壁面局部Nu數(shù)隨周向角θ的變化與文獻［1］中的數(shù)據(jù)對比，如圖3所示。通過對比可以發(fā)現(xiàn)周向角度θ從－90°到0°時模擬與文獻結果基本相同，當θ＞0°時Nu的值有些偏差，但總的來說模擬與文獻結果相差不大，且兩者的變化趨勢是相同的，于是本模擬所建模型是可靠的，可以得到超臨界CO2在管內流動的換熱特征。

圖3 模擬結果與文獻對比Fig.3 Comparison between simulation and literature data

3.2 局部傳熱系數(shù)

圖4為改變壁面溫度時上下母線的傳熱系數(shù)隨x/d的變化規(guī)律，其中重力加速度取g’=1 g，Tin=393.15 K，m=1.6×10－5kg/s。圖中上母線的傳熱系數(shù)明顯高于下母線，上母線傳熱系數(shù)先增加到峰值，跨過峰值之后傳熱系數(shù)迅速降低。由于冷卻時近壁面流體溫度低、密度大，中心流體溫度高、密度小，密度差產(chǎn)生浮升力使得近壁面處流體向下運動，而中間溫度較高的流體向上運動，于是在管道橫截面上產(chǎn)生二次環(huán)流，增強了壁面處的換熱。二次流使傳熱能力強的流體聚集在上母線處，而傳熱能力較差的流體分布在下母線處，于是上母線傳熱系數(shù)高于下母線。當流體處在擬臨界溫度附近區(qū)域時，比熱容會達到很大值，密度變化也非常大，此時引起的浮升力最強烈，使得在該區(qū)域傳熱系數(shù)達到峰值。壁面溫度越低上母線傳熱系數(shù)峰值越大，這主要是由于冷壁面溫度越低主流和壁面之間的溫差越大，引起的二次流越強，使傳熱系數(shù)越大?？邕^擬臨界溫度區(qū)域后壁面溫度越高上母線的傳熱系數(shù)越大。

圖4 壁面溫度對上下母線h隨x/d變化的影響Fig.4 Effect of wall temperature on variation of top and bottom generatrix heat transfer coefficient with x/d

圖5所示的是改變重力加速度g’時上下母線的傳熱系數(shù)隨Tb的變化規(guī)律，其中 Tw=298.15 K，Tin=393.15 K，m=1.6 ×10－5kg/s。由圖 4、圖 5 可見，在管道的入口處因入口效應傳熱系數(shù)均是先降低后升高，因為入口處溫度梯度較大，邊界層很薄，使傳熱系數(shù)較大。圖中重力加速度為零時上下母線傳熱系數(shù)相同，重力加速度不為零時，上母線傳熱系數(shù)高于下母線，隨著重力加速度的增加上母線的傳熱系數(shù)逐漸增加。重力加速度越大流體分層越明顯，管內上、下部分流體的溫差越大，二次流越強，使傳熱系數(shù)增加。傳熱系數(shù)在擬臨界溫度區(qū)域達到峰值，且重力加速度為零時傳熱系數(shù)在擬臨界溫度附近區(qū)域也存在一個峰值，這均是由超臨界CO2熱物性的巨大變化所引起的。

圖5 重力水平g’對上下母線h隨Tb變化的影響Fig.5 Effect of gravity acceleration on variation of top and bottom generatrix heat transfer coefficient with Tb

3.3 管內二次流變化

圖6給出的是改變入口Rein數(shù)時K隨x/d的變化規(guī)律，其中Tw=298.15 K，g’=1 g，Tin=393.15 K。隨著入口Rein數(shù)的增加K是逐漸減小的，因為入口Rein數(shù)越大對應的質量流量越大，浮升力引發(fā)的二次流較大，但是軸向速度增加的更多，慣性力更強，使得K越小，K在靠近入口達到峰值，跨過峰值后K有一個極值點，因為此點處在擬臨界點附近，物性變化劇烈，二次流較強。K最終趨于零，主要是由于流體已經(jīng)接近充分發(fā)展，流體之間的溫差很小使得密度差很小，于是浮升力引發(fā)的二次流幾乎為零。

圖6 入口Rein數(shù)對K隨x/d變化的影響Fig.6 Effect of Reinnumber on variation of K with x/d

圖7所示的是改變壁面溫度時K隨x/d的變化關系，其中 g’=1 g，Tin=393.15 K，m=1.6 ×10－5kg/s?？梢娫诳拷艿廊肟谔幊霈F(xiàn)K峰值，在未達到峰值前隨著冷壁面溫度的降低K是逐漸升高的，壁面溫度越低K更快的趨于零，因為在沒有達到峰值前，壁面溫度越低，壁面與流體之間的溫差越大，最終引發(fā)的二次流越大，K越大。當?shù)竭_峰值后，冷壁面溫度越低管內二次流越小，K越小。

圖7 冷壁面溫度對K隨x/d變化的影響Fig.7 Effect of wall temperature on variation of K with x/d

圖8中所示的是改變重力加速度g’時K隨x/d的變化規(guī)律，其中 Tw=298.15 K，Tin=393.15 K，m=1.6×10－5kg/s。K先增加到峰值然后減小到趨于零，且在靠近管道入口處達到峰值。隨著重力加速度的增加K也增大，因為重力加速度越大浮升力越大，產(chǎn)生的二次流越強。重力加速度為零時K在靠近入口處就降到了零，因零重力時不會產(chǎn)生浮升力，也就不會有二次流。零重力時二次流趨于零，不會增強傳熱，表現(xiàn)為圖5中零重力時上母線傳熱系數(shù)遠低于有重力時。圖6—圖8中可見，K均是在管道入口處較大，與文獻［6］結果相同，因為入口處溫度梯度大，物性變化劇烈，浮升力大，二次流強，于是K較大。

圖8 重力水平g’對K隨x/d變化的影響Fig.8 Effect of gravity acceleration on variation of K with x/d

3.4 局部范寧摩擦系數(shù)

圖9給出的是不同入口Rein數(shù)時上下母線的Cf×Reb隨x/d的變化規(guī)律，其中Tw=298.15 K，g’=1 g，Tin=393.15 K?？拷肟谔幧夏妇€的 Cf×Reb先升高到峰值后迅速降低，最后趨于恒定值，且上母線的Cf×Reb高于下母線，下母線的Cf×Reb在靠近入口處就從最大值降到很低，然后產(chǎn)生范圍不大波動，當流體接近充分發(fā)展后上下母線的Cf×Reb維持16不變，這與文獻［1］的結果是一致的。因為浮升力的作用使得上母線處流體的溫差大于下母線，于是靠近上母線處流體的物性變化要強于下母線，當流體接近充分發(fā)展時，溫差幾乎不變，流體的物性參數(shù)不再變化，接近為常物性參數(shù)流動，于是上下母線的Cf×Reb維持16不變。隨著入口Rein數(shù)的增加Cf×Reb的峰值變化不大，但是出現(xiàn)峰值的位置要靠后，主要是由于物性劇烈變化的位置隨Rein數(shù)的增加而靠后。

圖9 入口Rein數(shù)對Cf×Reb隨x/d變化的影響Fig.9 Effect of Reinnumber on variation of Cf×Rebwith x/d

圖10給出的是不同壁面溫度時上下母線的Cf×Reb隨x/d的變化規(guī)律，其中 g’=1 g，Tin=393.15 K，m=1.6×10－5kg/s。上母線 Cf×Reb的峰值隨冷壁面溫度的降低而增大，因為壁面溫度越低時擬臨界區(qū)域物性變化越劇烈，浮升力越大，進而引起的壁面剪切力越大。冷壁面溫度越低，Cf×Reb達到峰值的位置越靠近管道入口，最終上下母線的Cf×Reb的值收斂于16。

圖10 冷壁面溫度對Cf×Reb隨x/d變化的影響Fig.10 Effect of wall temperature on variation of Cf×Rebwith x/d

圖11給出的是改變重力加速度g’時上下母線的Cf×Reb隨x/d的變化規(guī)律，其中Tw=298.15 K，Tin=393.15 K，m=1.6 ×10－5kg/s。零重力時上下母線的Cf×Reb是相同的。有重力加速度時靠近入口處上母線的Cf×Reb明顯高于下母線，下母線的Cf×Reb低于零重力時，這是因為在有重力時靠近下母線處的流體的物性變化相對較小，使得壁面剪切力較小，從而下母線的Cf×Reb也較小。隨著重力加速度的減小上母線Cf×Reb的峰值是減小的，因為重力加速度越小管內流體的溫差越小，物性變化不再劇烈，峰值減小。下母線的Cf×Reb在有重力時相差不大。圖9—圖11可見，由于入口效應Cf×Reb在入口處均是先降低后增加，這是因為入口處速度梯度比較大，使剪切力較大，Cf×Reb較高。

圖11 重力水平g’對Cf×Reb隨x/d變化的影響Fig.11 Effect of g’on variation of Cf×Rebwith x/d

4 結論

對超臨界CO2在長1 000 mm，管徑0.5 mm的水平圓管內層流混合對流換熱做了數(shù)值模擬，分析了入口Rein、冷壁面溫度及重力水平對管道內換熱的影響，有如下結論:

(1)因入口效應管道上母線的傳熱系數(shù)在靠近入口處先降低后升高，在擬臨界溫度附近區(qū)域達到峰值，上母線的傳熱系數(shù)要遠大于下母線。冷壁面溫度越低上母線傳熱系數(shù)峰值越大。重力加速度越大管內傳熱系數(shù)越大，零重力時管道上母線傳熱系數(shù)要遠低于有重力加速度時。

(2)K均在入口處達到最大，且流體在管內充分發(fā)展時均趨于零。入口Rein數(shù)越大K越小，冷壁面溫度越低時K的峰值越大，重力加速度越大K越大，重力加速度為零時K為零。

(3)Cf×Reb在入口處因入口效應先降低后增加，上母線的Cf×Reb先增加后減小且在靠近入口處遠大于下母線，下母線的Cf×Reb在入口處就降到很低而后產(chǎn)生波動，最后二者都收斂于16，這與常物性的流動換熱是一致的。入口Rein數(shù)變化時Cf×Reb的峰值變化不大，冷壁面溫度越低或重力加速度越大時，上母線Cf×Reb的峰值越大。

1 Cao X L，Rao Z H，Liao S M.Laminar convective heat transfer of supercritical CO2in horizontal miniature circular and triangular tubes［J］.Applied Thermal Engineering，2011，31:2374－2384.

2 Liao S M，Zhao T S.A numerical investigation of laminar convection of supercritical carbon dioxide in vertical mini/micro tubes［J］.Progress in Computational Fluid Dynamics，2002(2):144－152.

3 Liao S M，Zhao T S.An experimental investigation of convection heat transfer to supercritical carbon dioxide in miniature tubes［J］.2002，45:5025－5034.

4 張麗娜，劉敏珊，董其伍，等.超臨界二氧化碳微細管內冷卻換熱研究［J］.低溫工程，2010(1):38－42.

5 曹 侃，董其伍，劉敏珊，等.超臨界CO2冷卻換熱特性數(shù)值模擬［J］.低溫工程，2012(1):56－60.

6 徐 峰，郭烈錦，白博峰.管內超臨界壓力水的混合對流換熱［J］.工程熱物理學報，2005，26(1):76－80.