重力热管流动与传热特性的数值研究

战洪仁， 张海松,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