重力熱管流動與傳熱特性的數值研究

戰洪仁， 張海松,2， 李春曉， 吳 眾

(1. 沈陽化工大學 能源與動力工程學院， 沈陽 110142；2. 中鋼集團鞍山熱能研究院有限公司， 遼寧鞍山 114000)

?

重力熱管流動與傳熱特性的數值研究

戰洪仁1， 張海松1,2， 李春曉1， 吳 眾1

(1. 沈陽化工大學 能源與動力工程學院， 沈陽 110142；2. 中鋼集團鞍山熱能研究院有限公司， 遼寧鞍山 114000)

采用計算流體力學(CFD)軟件，在流體體積(VOF)模型中引入Boussinesq近似模型，以水為工質對不同加熱功率下的重力熱管進行了數值模擬，將實驗、理論與數值結果三者結合對其內部流動與傳熱過程進行了分析.結果表明：CFD模擬得到的壁溫與實驗壁溫吻合較好，可以用CFD可視化定性分析加熱功率對重力熱管傳熱特性的影響；在10～80 W范圍內，隨著加熱功率的增大，蒸發段對流傳熱系數先增大后減小，冷凝段對流傳熱系數一直增大；當加熱功率超過一定值時，蒸發段傳熱性能惡化，CFD可視化結果顯示液池干涸，達到傳熱極限.

重力熱管； 蒸發； 冷凝； VOF模型； CFD可視化； 數值模擬

兩相閉式熱虹吸管(TPCT)又稱重力熱管，作為一項綠色節能技術，廣泛應用于能源領域尤其是工業余熱回收利用中[1].該技術包括兩相流和相變過程，傳統的實驗方法和理論分析難以對其內部熱質傳遞、蒸發冷凝及兩相流過程的具體細節進行分析研究.Shiraishi等[2]對TPCT的蒸發段和冷凝段分別進行了實驗研究，并針對其傳熱機理提出了比較簡明的傳熱模型，得到的分析結果與實驗結果吻合較好.Paramat thanuwat等[3]以自制納米流體為工質，分析了不同充液率、長徑比和加熱功率對重力熱管傳熱性能的影響.由于受客觀條件的制約，大部分理論分析和實驗探究沒有對其內部具體過程進行研究.為進一步探究封閉式TPCT的內部工作機理，越來越多的研究者開始嘗試設計可視化實驗裝置，來觀察其內部狀況.衛紅等[4]在重力熱管絕熱段設置了觀測口，對冷凝段是否存在攜帶作用進行了研究.劉國維等[5]設計了透明膜石英玻璃閉式熱管和非可視不銹鋼熱管，通過可視化實驗對不銹鋼熱管傳熱性能的影響因素進行了分析.韓振興等[6]應用電容層析成像技術對重力熱管冷凝段的流動換熱過程進行了可視化研究.然而，這些可視化實驗均基于一定的假設，且受操作環境制約，與熱管真實情況有一定的差距.近年來，有學者開始使用計算流體力學(CFD)軟件對熱管系統進行數值研究[7]，但還鮮有利用成熟的CFD模型來研究其內部傳熱傳質特性的報道[8].Schepper等[9]用流體體積(VOF)模型計算了水平熱管內兩相并流時的流動形態，發現貝克圖中的所有水平流態都可以使用CFD軟件進行模擬和計算.de Schepper等[10]在VOF模型中添加適當源項，對蒸汽裂化器的流動沸騰現象進行了建模研究，但該方法沒有應用到熱管系統.Alizadehdakhel等[11]結合重力熱管傳熱實驗，基于CFD軟件建立了熱管的二維模型，模擬得到的溫度分布與實驗測得的壁溫分布吻合較好.

筆者采用VOF模型探究了加熱功率對TPCT傳熱特性的影響規律，得到的壁面溫度與Liu等[12]的實驗壁溫分布吻合較好.為了使結果更加精確，在VOF模型中引入Boussinesq近似模型，通過Fluent軟件對TPCT內部具體傳熱傳質細節進行了探究，將實驗、理論分析以及數值研究三者結合，對重力熱管流動與傳熱特性進行了研究.

1 物理模型

根據典型的重力熱管結構和工作狀態，建立了重力熱管物理模型，如圖1所示.采用Shiraishi等[2]建立的經典傳熱模型，將TPCT分為三部分，即蒸發段、絕熱段和冷凝段.熱管模型大小與Liu等[12]的相同，長度分別為100 mm、100 mm和150 mm，直徑為8 mm.根據TPCT的實際情況進行如下假設[13]：(1) 熱管處于穩定工作狀態時，管內流體為穩態；(2) 忽略不凝結氣體和實驗槽道環境的影響；(3) 冷凝段冷凝只發生在近壁面處，且液膜沿壁面向下流動.

圖1 閉式重力熱管工作原理示意圖Fig.1 Working principle of the closed gravity heat pipe

2 數學模型

2.1 VOF模型

VOF的基本原理[14]是通過計算網格單元中的流體與網格的體積，構造出一個體積比函數F，從而確定出自由面流體的變化.由于單元控制容積中所有相的體積分數之和為1，故對汽液兩相有：

(1)

式中：φ為體積分數；下標v、l分別代表汽相和液相.

連續性方程：

(2)

式中：Sk為相變質量源項，kg/(m3·s)；ρ為密度，kg/m3；v為相的實際速率，m/s；t為時間，s；下標k代表第k相.

動量方程：

(3)

式中:p為壓強，N/m2；μ為動力黏度系數，Pa·s；g為重力加速度，m/s2；Fs為單位流體所受表面張力大小，N·m3；v為實際速度矢量，m/s.

能量方程：

(4)

式中:E為控制體比能，J/kg；Sh為相變能量源項，W/m3；Keff為有效導熱系數，W/(m2·K)；T為溫度.

混合相的物性參數取決于相的體積分數，則密度、動力黏度系數分別由下式確定：

(5)

(6)

2.2 Boussinesq模型

為了便于處理溫差引起的浮升力作用，對動量方程的重力項采用Boussinesq模型[15]進行近似處理.動量方程體積力變成：

(7)

式中:β為體脹系數，K-1；θ為溫差，K.

2.3 控制方程各源項的描述

2.3.1 質量源項和能量源項

對于汽液兩相流所涉及的控制方程，由軟件Fluent 14.5通過有限單元法求解，從而得到相應的溫度場和速度場.為了實現TPCT內部的相變過程，需要自定義函數(UDF).根據Schepper等[10]的研究結果，質量源項和能量源項如表1所示.表1中，β1、β2和β為質量轉移時間松弛因子，反映了蒸發與冷凝速率，其中β1、β2的符號與相間轉移方向有關，根據文獻[10]的研究，β1、β2的數值取0.1；下標s表示飽和水或飽和蒸汽溫度，K；ΔH為相變過程中的焓變，反映了單位物質的能量變化，J.

表1 自定義函數Tab.1 User-defined functions

2.3.2 動量源項

表面張力模型是由Brackbill等[16]提出的連續表面力模型，將表面張力在汽液界面上的面積分轉變為體積分，并添加到動量方程的源項中：

(8)

式中：C為曲率，m-1；σ1,v為表面張力系數，N/m.

2.4 邊界條件

TPCT底部和頂部的邊界條件為絕熱邊界條件；蒸發段為恒熱流密度邊界條件；絕熱段為絕熱的非滑移邊界條件；冷凝段為自然對流邊界條件，則有：

(9)

(10)

qw=hΔt

(11)

(12)

(13)

(14)

式中：λ為導熱系數，W/(m2·K)；h為對流傳熱系數，W/(m2·K)；q為熱流密度，J/(m2·s)；下標const表示定值；n為換熱表面的外法線；v和w分別表示壁面切向和法向速度，m/s.

3 數值計算結果及分析

3.1 CFD熱管壁溫模擬值與實驗值的比較

模擬熱管運行壓力為7.4 kPa、蒸發段加熱功率為80 W，當質量、能量和速度的殘差小于10-4時，認為數值計算過程已經基本收斂.當熱管壁溫Tw趨于穩定時，將CFD計算得到的熱管壁溫模擬值與實驗值進行對比，如圖2所示.

圖2 模擬與實驗壁溫的比較Fig.2 Comparison of wall temperature between simulated and experimental results

由圖2可知，蒸發段壁溫平均誤差為9%，絕熱段壁溫平均誤差為2.5%，冷凝段壁溫平均誤差為3%.由上面比較可知，CFD模擬值與實驗值偏差較大，這是由于實驗測得的壁溫并非熱管實際壁溫，而是金屬絲表面的溫度.同時，冷凝段壁溫與實驗值吻合較好，驗證了上述分析.

3.2 蒸發段與冷凝段的對流傳熱系數

根據對流傳熱速率計算公式，可以將蒸發段和冷凝段的對流傳熱系數分別定義為：

(15)

(16)

3.3 加熱功率對傳熱性能的影響

圖3給出了不同充液率R下，TPCT蒸發段和冷凝段的對流傳熱系數隨加熱功率的變化.由圖3可知，在不同充液率條件下，隨著加熱功率的增大，蒸發段對流傳熱系數先增大后減小，冷凝段對流傳熱系數一直增大.

如圖3(a)所示，蒸發段對流傳熱系數隨加熱功率的增大而增大，當超過一定值時傳熱系數不再增大.這主要是由于隨著加熱功率的增大，TPCT壁面溫度升高，沸騰汽化核心增多，單位時間汽泡生成率變大，而加熱面上汽泡的生長和脫離對加熱面附近的液體產生了強烈的擾動，使蒸發段的對流換熱得到強化.當加熱量超過一定值時，對流傳熱系數減小，壁面溫度逐漸升高而達到了熱管傳熱極限，傳熱性能惡化.如圖3(b)所示，冷凝段對流傳熱系數隨著加熱功率的增大而增大.蒸發段蒸汽生成率的增大使得冷凝段液膜厚度不斷增大，理論上會導致冷凝對流傳熱系數的減小，但是蒸汽速度的增加會在汽液交界面處形成界面波，界面波不僅增大了傳熱表面積，而且加劇了對冷凝段液膜的擾動，故也可能強化了換熱.

(a) 蒸發段

(b) 冷凝段圖3 加熱功率對傳熱系數的影響Fig.3 Effect of heating power on the heat transfer coefficient

根據Gross等[17]的研究結果，當雷諾數Re1>5時，界面波所帶來的增強作用足以抵消液膜厚度增大所帶來的消弱作用，總體上傳熱系數增大.冷凝段出口處的液膜雷諾數Re1由下式確定：

(17)

式中：v1為液膜的速度，m/s；hfg為汽化潛熱，J/kg.

經計算可知，當加熱功率由10 W增加到80 W時，冷凝段出口處的液膜雷諾數在20.8～285內，而計算得到液膜厚度約為0.8×10-3～2.5×10-3mm，故隨著加熱功率的增大，冷凝對流傳熱系數一直呈增大趨勢.CFD計算得到的冷凝段換熱規律與Gross等的理論分析結果一致.

3.4 CFD可視化分析

對充液率為30%、加熱功率為50 W的TPCT模擬結果進行CFD可視化分析.圖4為熱管內部過程的冷凝段液膜體積分數云圖和速度矢量圖.由圖4(a)可知，蒸汽在頂端沿軸向向兩邊運動，在近壁處冷凝成液體，之后逐漸形成一定厚度的液膜.因為沿壁面方向蒸汽冷凝成水放出潛熱，所以壁面溫度逐漸降低.由圖4(b)可知，液膜在重力作用下沿壁面向下流動，從矢量圖可以看到，液膜在壁面兩側的速度方向向下.流動CFD可視化結果與熱管穩定時的工作基本特征是一致的[18].

(a)體積分數云圖(b)速度矢量圖

圖4 冷凝段液膜體積分數云圖和速度矢量圖

Fig.4 Volume fraction contour and velocity vector chart of liquid film at condensation section

一段時間后，體積分數云圖和溫度分布圖如圖5所示.在TPCT蒸發段壁面有微小汽泡形成，靠近壁面及液體內部存在明顯的溫度梯度.根據分子動力學理論可知[19]，液體中分子能量分布不均勻，能量較大的活化分子隨機聚集形成了暫時的局部微小低密度區，這些低密度區被認為是具有一定半徑的微小汽泡.經過一段時間后可以清晰地看出，溫度梯度明顯增大，汽化核心數增多且產生的汽泡互不干擾.此時，主流液體溫度沒有達到飽和溫度，壁面溫度大于飽和溫度，二者之差約為2 K.

(a)體積分數云圖(b)溫度分布云圖

圖5 蒸發段水的體積分數云圖和溫度分布云圖

Fig.5 Volume fraction contour and temperature field distribution of water at evaporation section

再經過一段時間，管內液體溫度逐漸達到飽和溫度，甚至超過飽和溫度，此時，液體過熱度達到5 K左右.這時熱量首先由加熱面傳給液體，然后再由過熱液體傳給微小汽泡，另一部分熱量對流傳給液體本體，此時有更多的微小汽泡脫離壁面進入液體，正是這種擾動強化了對流換熱.

除此之外，可以看出汽泡明顯變大，并觀察到汽泡的脫離與合并過程.由成核理論和汽泡動力學基本理論可知[19]，在活化凹坑上首先形成汽化核心，在各種力和熱的作用下，汽化核心會繼續長大.當表面張力平衡不了汽泡內外壓差時，汽泡繼續長大，然后脫離壁面進入液體，形成下一個汽泡，進行周期性替換.在上升過程初始階段，汽泡先相互接觸，隨著熱流密度的增大，逐漸形成規則的汽泡.蒸發段液池水的體積分數隨時間的變化如圖6所示.由圖6可知，CFD流動可視化觀察到的汽泡運動與汽泡動力學基本理論一致，最后汽泡到達液面進入汽相中.

圖6 蒸發段液池水的體積分數隨時間的變化Fig.6 Change of water volume fraction in liquid pool with time at evaporation section

從以上分析可以看出，蒸發段汽泡行為和冷凝段冷凝過程與經典的熱力學成核理論相符合，蒸發段CFD可視化所呈現的流動過程也與陳崗等[20]的實驗結果一致.陳崗等在鍍膜玻璃管下粘貼了金屬管，自制TPCT并對其加熱，觀察到隨熱流密度的增大，蒸發段大量汽泡脫離壁面并發生了長大聚合過程.對TPCT繼續加熱，發現當加熱功率超過60 W時，TPCT的傳熱系數減小.此時，繼續增大加熱功率，觀察到TPCT蒸發段中汽泡聚合連成一片，此時壁面溫度升高，傳熱性能惡化.

分析原因可知，此時汽泡的生長速度大于汽泡躍離加熱面的速度，致使汽泡聚集覆蓋在加熱面上，形成一層蒸汽膜，增大了熱阻，導致管壁面溫度升高達到傳熱極限，進而使蒸發段傳熱性能惡化，這也解釋了圖3(a)中蒸發段存在最大加熱功率的原因.數值結果所呈現的換熱規律與Liu等[12]的實驗結果一致，熱管整個過程不同時刻的溫度分布云圖如圖7所示.

圖7 溫度分布云圖隨時間的變化Fig.7 Change of temperature field distribution with time

4 結 論

(1) 在VOF模型中引入Boussinesq近似模型，將已有的傳熱傳質關系式通過自定義函數轉化為相應控制方程源項，實現了TPCT內部流動與傳熱過程的模擬，較好地預測了TPCT蒸發段汽泡成核和運動的過程.

(2) 在加熱功率為12～80 W的范圍內，TPCT蒸發段對流傳熱系數隨加熱功率的增大先增大后減小，冷凝段對流傳熱系數隨加熱功率的增大一直呈增大趨勢.

(3) 當充液率為定值時，蒸發段傳熱系數隨加熱功率增大到一定值時不再增大，液池內部發生干涸現象，此時達到了TPCT的傳熱極限，最佳加熱功率在60 W左右.

(4) CFD可視化較好地解釋了加熱功率引起傳熱系數變化的原因，可以用經典熱力學成核理論來定性分析蒸發段的汽泡行為和冷凝成核過程，數值可視化研究結果與實驗探究、理論分析結果一致.

[1] BARZI Y M, ASSADI M. Evaluation of a thermosyphon heat pipe operation and application in a waste heat recovery system[J]. Experimental Heat Transfer, 2014, 28(5): 493-510.

[2] SHIRAISHI M, KIKUCHI K, YAMANISHI T. Investigation of heat transfer characteristics of a two-phase closed thermosyphon[J]. Journal of Heat Recovery Systems, 1981, 1(4): 287-297.

[3] PARAMATTHANUWAT T, BOOTHAISONG S, RITTIDECH S,et al. Heat transfer characteristics of a two-phase closed thermosyphon using deionized water mixed with silver nano[J]. Heat Mass Transfer, 2010, 46(3): 281-285.

[4] 衛紅，馬同澤，陳煥倬. 兩相閉式熱虹吸管內凝結換熱的研究[J]. 工程熱物理學報，1990，11(2)：182-187.

WEI Hong, MA Tongze, CHEN Huanzhuo. Condensation heat transfer in two-phased closed thermosyphon[J].Journal of Engineering Thermophysics, 1990, 11(2): 182-187.

[5] 劉國維，單巖昆，黃鴻鼎．兩相閉式熱虹吸管強化傳熱研究[J]. 化工學報，1991，11(1)：66-71.

LIU Guowei, SHAN Yankun, HUANG Hongding. Heat transfer enhancement in a two-phased closed thermosyphon[J].CIESC Journal, 1991,11(1): 66-71.

[6] 韓振興，王冬驍，王飛.重力熱管冷凝段運行特征的可視化實驗研究[J]. 化工學報，2014，65(8)：2934-2938.

HAN Zhenxing,WANG Dongxiao,WANG Fei. Visual experimental study on operation characteristics of condensation segment of gravity-assisted heat pipe[J]. CIESC Journal, 2014, 65(8): 2934-2938.

[7] 閆偉偉，葛仕福，李揚. 槽式太陽能 DSG 系統集熱管內強化傳熱的數值模擬[J]. 動力工程學報，2013，33(8)：550-555.

YAN Weiwei，GE Shifu，LI Yang. Numerical simulation on heat transfer enhancement in parabolic trough solar collector of DSF systems[J].Journal of Chinese Society of Power Engineering，2013，33(8)：550-555.

[8] KARTHIKEYAN V K, RANACHANDRAN K, PILLAI B C. Effect of nanofluids on thermal performance of closed loop pulsating heat pipe[J]. Experimental Thermal and Fluid Science, 2014, 54: 171-178.

[9] de SCHEPPER S C K, HEYNDERICKX G J, MARIN G B. CFD modeling of all gas-liquid and vapor-liquid flow regimes predicted by the Baker chart[J].Chemical Engineering Journal,2008, 138: 349-357.

[10] de SCHEPPER S C K, HEYNDERICKX G J, MARIN G B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker[J].Computers and Chemical Engineering, 2009, 33(1): 122-132.

[11] ALIZADEHDAKHEL A, RAHIMI M， ALSAIRAFI A A. CFD modeling of flow and heat transfer in a thermosyphon[J].International Communications in Heat and Mass Transfer, 2010, 37(3): 312-318.

[12] LIU Z H, LI Y Y, BAO R. Thermal performance of inclined grooved heat pipes using nanofluids[J]. International Journal of Thermal Sciences, 2010,49(9): 1680-1687.

[13] ASMAIE L, HAGHSHENASFARD M, MEHRABANI-ZEINABAD A，et al. Thermal performance analysis of nanofluids in a thermosyphon heat pipe using CFD modeling[J]. Heat Mass Transfer, 2013, 49(5): 667-678.

[14] FAGHRI A. Heat pipe science and technology[M]. Washington, USA: Taylor and Francis, 1995: 285-286.

[15] JOHANNSEN M. On the validity of the Boussinesq approximation for the elder problem[J].Computational Geosciences, 2003, 7(3):169-182.

[16] BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J]. Journal of Computational Physics, 1992, 100(2): 335-354.

[17] GROSS U. Reflux condensation heat transfer inside a closed thermosyphon[J]. International Journal of Heat of Mass Transfer, 1992, 35(2): 279-294.

[18] 莊駿，張紅. 熱管技術及其工程應用[M]. 北京：化學工業出版社，2000.

[19] 施明恒，甘永平，馬重芳. 沸騰和凝結[M]. 北京：高等教育出版社，1995: 24-54.

[20] 陳崗，辛明道，陳遠國. 兩相閉式熱虹吸管內的流動與傳熱[J]. 工程熱物理學報, 1987, 8(2): 149-152.

CHEN Gang, XIN Mingdao, CHEN Yuanguo. Flow and heat transfer in two-phase closed thermosyphons[J]. Journal of Engineering Thermophysics, 1987, 8(2): 149-152.

Simulation Study on Flow and Heat-transfer Characteristics of Gravity Heat Pipes

ZHANHongren1,ZHANGHaisong1,2,LIChunxiao1,WUZhong1

(1.School of Energy and Power Engineering, Shenyang University of Chemical Technology,Shenyang 110142, China; 2. Sinosteel Anshan Research Institute of Thermo-energy Co., Ltd., Anshan 114000, Liaoning Province, China)

By using computational fluid dynamics (CFD) software and introducing Boussinesq approximate model to the volume of fluid (VOF) model, numerical simulations were conducted on the gravity heat pipe under different heating power with water as the working medium, so as to study the flow and heat-transfer characteristics of the gravity heat pipe based on experimental, theoretical and simulation data. Results show that the wall temperature obtained by CFD simulation agrees well with that of experiment, and the CFD visualization can be used to qualitatively analyze the effects of heating power on the heat-transfer characteristics of gravity heat pipes. In the range of 10-80 W, with the increasing of heating power, the heat-transfer coefficient at evaporation section firstly increases and then decreases, while that at condensation section keeps increasing. Once the heating power exceeds a certain value, the heat-transfer performance at evaporation section deteriorates, when a dried up liquid pool could be observed through CFD visualization, indicating a heat transfer limit is already reached.

gravity heat pipe; evaporation; condensation; VOF model; CFD visualization; numerical simulation

2016-04-25

2016-06-17

國家自然科學基金資助項目(61473056，61104157)

戰洪仁(1964-)，女，山東蓬萊人，教授，博士，主要從事強化傳熱與節能技術方面的研究與利用. 張海松(通信作者)，男，碩士，電話(Tel.): 15040788554；E-mail：haisongz@yeah.net．

1674-7607(2017)07-0540-06

TK172.4

A

470.10