正方形小通道内的冷凝换热数值模拟

李盼盼　 陈振乾

(东南大学能源与环境学院, 南京 210096)



正方形小通道内的冷凝换热数值模拟

李盼盼 陈振乾

(东南大学能源与环境学院, 南京 210096)

采用VOF模型对制冷剂R134a在水平放置的边长为1 mm、长度50 mm正方形小通道内的凝结换热过程进行数值模拟计算,研究了入口质量流速、热流密度和冷量施加方式对小通道内冷凝换热的影响．结果表明:冷凝液最先出现在通道拐角处,随后在表面张力的作用下逐渐布满整个通道截面;沿程换热系数在入口处较大,然后沿流向急剧减小,最后趋于稳定;在恒定热流密度条件下,随着质量流速的增加,气液界面的剪切力逐渐增加,冷凝液膜随之变薄,从而导致换热系数变大;当入口质量流速恒定时,热流密度的改变对冷凝换热系数影响不明显;冷量均衡施加于通道边界比集中施加于某一边界液膜分布更均匀,更有利于冷凝换热．

正方形小通道;两相流;冷凝换热;数值模拟;VOF模型

小通道中的流动冷凝换热是微热管、紧凑型传热器以及热控制系统中的关键热物理过程[1],研究小通道中的冷凝换热过程,对于相关能量传输装置的优化设计,提高装置运行过程中的换热效率,具有重要的学术意义和市场价值．

冷凝过程的换热效率依赖于如通道的截面形状、结构尺寸、流体性质、质量流量等诸多因素．da Riva等[2-4]对制冷剂R134a在半径为1 mm的小通道内的冷凝过程进行了一系列的研究，结果表明:① 低质量流速时,重力主导冷凝过程,表面张力的作用可以忽略,而高质量流速时,剪切力作用凸显;② 在圆形截面通道内,通道上半部分的液膜厚度沿流向基本不变,在重力的作用下,通道下半部分的液膜厚度沿流向增加;③ 矩形通道的冷凝换热系数明显大于圆形通道．王勋等[5]完成了制冷剂R134a在微通道中的相变传热特性实验,认为质量流速、微通道的尺寸等是影响传热的主要因素．Shin等[6]实验研究了圆形截面和矩形截面小通道内的冷凝过程,研究结果表明热流密度对通道的换热系数和压降均没有明显影响;冷凝换热系数随着通道半径的减小和质量流速的增加而显著增大．杨英英等[7]采用可视化的方法,观测了工质R32在内径2 mm的水平光滑圆管内冷凝换热的流型,随着流量的增大气液界面剪切力增大,流型由环波状流转换成间歇流的干度推迟,环状流的流型区域增加．胡灿等[8]还收集了R22,R134a等6种工质的1 183个小通道内的冷凝换热试验数据点,对小通道内的冷凝换热模型进行分析,表明应根据工质和工况来选择合适的模型．

Wang等[9-10]提出了一个水平放置的小通道膜状冷凝的理论模型,并利用此模型对工质R134a,R22和R410A在半径为0.5～5.0 mm的正方形和三角形截面的小通道冷凝进行了研究.结果表明该理论模型可以较好地预测通道内的冷凝换热和液膜分布,短通道内重力对液膜在通道底部汇聚作用并不明显,通道足够长时通道底部液膜厚度明显高于通道顶部,垂直于流向截面上的气相分布近似为圆形;对于非圆形截面,通道拐角处最先出现冷凝液,此处的液膜厚度最厚,换热系数较小．陈永平等[11]对矩形微通道中的环状冷凝进行三维数值模拟,结果表明薄液膜区液膜将沿程逐渐增厚,到达极值后再逐渐变薄．在通道截面中薄液膜区的传热系数大于弯月面,最大局部传热系数及壁面最高温度皆位于薄液膜区和弯月面的连接处．在冷凝起始段,通道横截面平均传热系数沿程急剧减小至一极值;在此后很长的一段距离内,则基本保持不变;直至接近环状冷凝终点时又再次沿程减小．Ganapathy等[12]采用VOF方法数值模拟了高质量流速(245～615 kg/(m2·s))的制冷剂R134a在半径为100 μm小通道内且高热流密度(200～800 kW/m2)条件下的冷凝换热过程,模拟所得到的流型与实验所得流型吻合较好;在恒热流时质量流速对冷凝换热具有显著影响,而在恒定质量流速时,热流密度的改变并不明显影响冷凝换热．

以上对于小通道内冷凝换热的研究大多基于恒温或者高热流密度边界条件下进行的,随着小通道冷凝在电子元器件散热等领域应用的发展,需要对较小热流密度的冷凝进行探索,同时也需要对冷量的施加方式进行研究．为此,本文基于VOF模型对制冷剂R134a在正方形小通道内的冷凝过程进行数值模拟,结果表明冷凝最先出现在通道拐角处,质量流速对冷凝过程影响显著,热流密度的影响可以忽略,冷量均衡施加于通道截面时,液膜分布均匀更有利于冷凝换热．

1　模型

1.1VOF模型

VOF模型目前已被广泛用于数值模拟两相流动过程相界面的追踪．在VOF模型中相的体积分数作为变量,在每个计算单元中各相的体积分数总和为1,采用Nvaier-Stokes方程组描述其流动问题,质量、动量和能量守恒方程为

(1)

(2)

(3)

(4)

(5)

式中,ρ为密度;t为时间；V为速度矢量;u, v, w分别为X, Y, Z轴的速度分量;P为压力;μ为动力黏度;Cp为定压比热;k为湍流动能;T为温度;fX, fY, fZ为X, Y, Z轴的体积力;QT为冷凝过程中单位体积产生的热量．

本文中,设定气相为第一相,液相为第二相,VOF方程可表示为

(6)

1.2相变模型

(7)

(8)

(9)

式中,db为液滴直径;β为调节系数;βc为传质时间弛豫系数[14];Ts为蒸汽饱和温度．

1.3表面张力CSF模型

小通道内的冷凝换热机理不同于大通道,主要是因为小通道内的重力、切应力和表面张力的相互作用不同．小通道内的表面张力采用连续表面力模型(CSF),附加的表面张力通过在动量方程中添加源项的方法来实现[15]．交界面的压力差取决于表面张力系数和表面曲率,由垂直方向的2个半径r1和r2决定,即

ΔP=P2-P1=σκ

(10)

(11)

曲率κ可表示为

(12)

表面上的力可表示为体积力fvol,即

(13)

本文中只有气相和液相,为此方程(13)可推导为

(14)

2　模型建立

本文对制冷剂R134a在边长d=1mm、通道长度L=50mm的正方形小通道内的冷凝过程进行数值模拟．计算过程中,气相和液相均采用K-ε湍流流动模型,通道入口采用质量流量进口边界,出口采用压力出口边界．为了模拟热流量大小和冷量施加方式对冷凝换热的影响,其余4个通道边界采用定热流边界条件,工质流动方向沿Z轴．采用VOF模型追踪气液界面,为了保证计算精度及尽可能减少计算量,选取HRIC格式,其他方程采用third-orderMUSCL格式．模拟过程中考虑重力、表面张力和剪切力的作用,并对模拟过程的网格无关性进行验证．由于冷凝成核过程物理机制较为复杂,初始假定环通道有一层极薄的液膜,并对初始液膜厚度的模拟结果的无关性进行了验证,假定初始液膜厚度为10μm．

3　结果分析与讨论

3.1模型验证

通道出口蒸汽干度χ是衡量通道内冷凝量的一项参数,将采用本文模型计算的蒸汽干度与相同边界条件下文献[16]的出口蒸汽干度进行比较,结果如表1中所示,二者的平均绝对误差为7.6%．而且矩形截面小通道出口蒸汽干度小于圆形截面通道出口蒸汽干度,即矩形截面通道内的冷凝量更大,更有利于冷凝换热．

表1　蒸汽出口干度验证模型准确性

3.2冷凝过程分析

制冷剂R134a在正方形通道中流动时,速度分布如图1所示,在通道的拐角处流速最低,截面中心处速度最高．通道边界各处的热流密度相同时,通道横截面上拐角处的温度最低,如图2所示．

图1　通道横截面上的速度分布示意图

由于通道横截面上温度最低出现在拐角处,使得拐角处最先出现冷凝液,如图3(a)的气液分布

图2　通道横截面的温度分布示意图

(a) 通道横截面上的气液界面分布(Z=0,10,20,30,40,50 mm)

(b) Y=0.5 mm截面上的气液界面分布

所示,图中红色部分为液相,蓝色部分为气相．由于通道入口处的冷凝量极小,很难观察拐角中间位置处的液膜．沿蒸汽流动方向,流动速度逐渐减小,横截面上的温度也逐渐降低,冷凝量逐渐增多,液膜逐渐增厚,但是相对于通道宽度,此时的液膜仍然极薄,如图3(b)所示．

通道边界上施加的冷量在向通道内部传递蒸汽过程中必须先通过冷凝液膜,此时液膜厚度成为热量传递过程中的主要阻力．由通道横截面上气液界面分布可知,在通道拐角处液膜最厚,相应的传热阻力最大,因此局部换热系数最小,而在2个拐角相连处液膜最薄,产生局部换热系数最大值．

图4是液膜沿通道增长示意图．如图所示,最初在拐角处形成的液膜气液界面近似为一条斜线,

图4　通道内液膜增长示意图

然后在表面张力的作用下,拐角处被不断拉伸,气液界面由直线拉伸为近似弧形,此时由于截面上的冷凝量极小,拐角连接处的液膜几乎不可见．沿通道流向,冷凝量增加,通道拐角处和拐角连接处的液膜厚度都不断增加．

如果没有表面张力的作用,则冷凝量沿流向增加,通道内液膜增厚,在拐角处的气液界面会一直保持近似直线的形状[17]．在表面张力的作用下,拐角连接处的液膜被拉至拐角处,其他位置处的液膜相对变薄,这有利于热量传递．矩形截面通道相对于圆形截面更有利于冷凝换热的主要原因在于表面张力的存在．

图4中冷凝液膜在整个通道截面内的分布较为均匀,并未看到冷凝液在通道底部汇聚的现象,主要原因是冷凝通道较短,冷凝量较少,重力对液膜的作用不明显;当通道足够长时会出现在通道底部冷凝液膜厚度明显较其他边界厚的现象．Wang等[9-10]模拟600 mm长的矩形通道冷凝过程结果也说明了这一现象,在50 mm位置截面上不能观察到重力使底部液膜增厚的现象,而在600 mm位置截面上则可以明显观察到通道底部液膜厚度大于顶部液膜．

3.3质量流速和热流密度对冷凝换热的影响

为了研究质量流速和热流密度对冷凝换热强度的影响,分别计算了质量流速G=50,300,500 kg/(m2·s)和热流密度q=-4,-10,-100 kW/m2时通道内平均换热系数和沿程换热系数,如图5和图6所示．

图5和图6表明,质量流速对通道内平均换热系数h和沿程换热系数hZ的影响都比较显著,而热流密度对两者的影响皆不明显．在图6中,由于入口效应的存在,在通道初始阶段换热系数较大,随后沿流向换热系数趋于稳定．

在恒定热流密度的工况下,入口质量流速对通

图5　质量流速、热流密度对通道内平均换热系数的影响

(a) 热流密度对通道沿程换热系数的影响(G=50 kg/(m2·s))

道内冷凝换热影响显著的原因在于：随着质量流速的增加,存留于通道内的液膜厚度变薄,如图7(a)所示,质量流速为500 kg/(m2·s)时的液膜厚度比50 kg/(m2·s)时的液膜厚度明显要薄,这有利于热量传递．在较低的入口质量流速时,通道内液相在表面张力的作用下被拉伸至拐角处,使得通道内的液膜相对减薄,从而使热量传递的阻力减小,通道内的凝结换热被强化．当入口质量流速较高时,表面张力的作用被弱化,气液界面的剪切力主导了通道内的冷凝过程．增大的气液界面剪切力使得伴随蒸汽带出通道的液相增加,通道内的残留冷凝液减少,通道边界的液膜厚度减薄,增大了通道的换热系数．

在恒定入口质量流速的工况下,热流密度对通道内冷凝换热的影响不明显．如图7(b)所示,热流密度为-4 kW/m2时的液膜厚度和-10 kW/m2时的液膜厚度差距不明显,从而换热强度近似．

(a) 质量流速对通道横截面液膜厚度的影响(Z=50 mm)

(b) 热流密度对通道横截面液膜厚度的影响(Z=50 mm)

3.4冷量施加方式对冷凝换热的影响

为了研究冷量施加方式对冷凝换热的影响,模拟了G=50 kg/(m2·s)时通道施加相同大小冷量但施加方式不同对冷凝换热强度的影响．冷量相同且均衡施加时,热流密度q=-4 kW/m2,而当冷量仅施加在其中一个通道边界时热流密度q=-16 kW/m2．

图8为不同冷量施加方式对通道内沿程换热系数的影响．当冷量均衡施加在底、顶、前、后4个表面时的沿程换热系数明显比冷量只施加于其中

图8　不同冷量施加方式对通道沿程换热系数的影响

一个通道表面时的换热系数高，即冷量均衡施加于通道边界更有利于冷凝换热．产生这种现象的主要原因在于冷量均衡施加时通道横截面上的液膜近似均匀分布,而冷量仅施加于其中一个通道表面时,液相堆积在被施加冷量边界处,此处的液膜厚度明显较厚,不利于冷凝换热的进行,如图9(a)所示．

图9显示了冷量分别集中施加于通道顶部、底部和侧面3种方式下的液膜分布．冷量集中施加时,液膜分布不均匀,液相主要集聚在冷量施加侧的通道边界处,如图9(a)所示．冷量分别施加于通道顶部和侧面时的液膜分布形态不能显示出重力对液相的汇聚作用,即通道底部的液膜厚度并没有明显高于其他方向处,产生这一现象的主要原因在于通道长度L=50 mm相对较短,因此重力作用不明显．

(a) 对短通道的影响(Z=50 mm)

(b) 对长通道的影响(Z=500 mm)

图9(b)显示了不同冷量施加方式对长通道(Z=500 mm)内液膜分布的影响,结果表明在冷量分别施加在通道顶部、底部和侧面时,最终液膜分布形态相似,皆在重力的作用下汇聚在通道底部,即不同的冷量施加方式下液膜最终在通道横截面上的分布是相同的．

4　结论

1) 通道拐角处由于蒸汽流速较缓,温度较低,从而最先形成液膜,且此处液膜最厚,换热系数最小．

2) 恒定热流密度时,质量流速的增加导致气液界面剪切力增加,从而通道内的液膜厚度相对减薄,换热系数增大;恒定质量流速时,热流密度的改变对换热系数的影响不明显．

3) 冷量均衡施加于通道边界时,横截面上液膜厚度较均匀,相对于冷量集中施加于其中一个通道边界时,液膜厚度较薄,换热系数更大;短通道时,重力作用不明显,冷量分别集中施加于顶部、底部和侧面的液膜形态近似,在施加冷量的边界处液膜较厚,当通道足够长时,重力作用凸显,冷凝液均汇聚在通道底部．

References)

[1]陈永平, 肖春梅, 施明恒, 等.微通道冷凝研究的进展与展望[J]. 化工学报, 2007, 58(9): 2153-2160. DOI:10.3321/j.issn:0438-1157.2007.09.001.

Chen Yongping, Xiao Chunmei, Shi Mingheng, et al. Review of condensation in microchannels [J].JournalofChemicalIndustryandEngineering, 2007, 58(9): 2153-2160. DOI:10.3321/j.issn:0438-1157.2007.09.001. (in Chinese)

[2]da Riva E, del Col D. Effect of gravity during condensation of R134a in a circular minichannel [J].MicrogravityScienceandTechnology, 2011, 23(S1):87-97. DOI:10.1007/s12217-011-9275-4.

[3]da Riva E, del Col D. Numerical simulation of laminar liquid film condensation in a horizontal circular minichannel [J].JournalofHeatTransfer, 2012, 134(5): 807-824.

[4]del Col D, Bortolin S, Cavallini A, et al. Effect of cross sectional shape during condensation in a single square minichannel [J].InternationalJournalofHeatandMassTransfer, 2011, 4(17):3909-3920.DOI:10.1016/j.ijheatmasstransfer.2011.04.035.

[5]王勋, 王文, 耑锐. R134a制冷剂在微通道中的相变传热特性[J]. 低温工程, 2009, 2: 10-14.

Wang Xun, Wang Wen, Zhuan Rui. Experimental of condensation and heat transfer for R134a refrigerant in micro-channel heat exchanger [J].Cryogenics, 2009, 2: 10-14.(in Chinese)

[6]Shin J S, Kim M H. An experimental study of flow condensation heat transfer inside circular and rectangular mini-channels [J].HeatTransferEngineering, 2005, 26(3):36-44. DOI:10.1080/01457630590907185.

[7]杨英英, 李敏霞, 马一太. 水平光滑细管内R32冷凝换热的流行特性[J]. 化工学报, 2014, 65 (2): 445-452.

Yang Yingying, Li Minxia, Ma Yitai. Characteristics of flow pattern for condensation heat transfer of R32 in horizontal small tube [J].JournalofChemicalIndustryandEngineering, 2014, 65(2): 445-452. (in Chinese)

[8]胡灿, 李敏霞, 马一太, 等.小通道内的冷凝换热模型分析[J]. 机械工程学报, 2012, 48(24): 134-140. DOI:10.3901/JME.2012.24.134.

Hu Can,Li Minxia,Ma Yitai, et al. Analysis of condensation heat transfer model in small channels [J].JournalofMechanicalEngineering, 2012, 48(24): 134-140. DOI:10.3901/JME.2012.24.134. (in Chinese)

[9]Wang H S, Rose J W. A theoretical model of film condensation in square section horizontal microchannels [J].ChemicalEngineeringResearchandDesign, 2004, 82(4):430-434.DOI:10.1205/026387604323050137.

[10]Wang H S, Rose J W. A theory of film condensation in horizontal noncircular section microchannels [J].HeatTransfer, 2005, 127: 1096-1105.

[11]陈永平, 吴嘉峰, 施明恒, 等.矩形微通道中环状冷凝的三维数值模拟[J]. 化工学报, 2008, 59(8): 1923-1929. DOI:10.3321/j.issn:0438-1157.2008.08.007.

Chen Yongping, Wu Jiafeng, Shi Mingheng, et al. Three dimensional simulation for steady annular condensation in rectangular microchannels [J].JournalofChemicalIndustryandEngineering, 2008, 59(8): 1923-1929. DOI:10.3321/j.issn:0438-1157.2008.08.007. (in Chinese)

[12]Ganapathy H, Shooshtari A, Choo K, et al. Volume of fluid-based numerical modeling of condensation heat transfer and fluid flow characteristics in microchannels [J].InternationalJournalofHeat&MassTransfer, 2013, 65(5): 62-72. DOI:10.1016/j.ijheatmasstransfer.2013.05.044．

[13]Lee W H.Apressureiterationschemefortwo-phaseflowmodeling[M]. Washington, USA: Hemisphere Publishing, 1980: 407-431.

[14]ANSYS Inc. ANSYS fluent theory guide 15.0 [EB/OL]. (2013-11)[2015-11]. http://www.ansys.com.

[15]Rose J W. Surface tension effects and enhancement of condensation heat transfer [J].ChemicalEngineeringResearchandDesign, 2004, 82(4): 419-429. DOI:10.1205/026387604323050128．

[16]Bortolin S, da Riva E, del Col D. Condensation in a square minichannel: Application of the VOF method [J].HeatTransferEngineering, 2014, 35(2):193-203. DOI:10.1080/01457632.2013.812493.

[17]Liu Z Y, Sunden B, Yuan J L. Numerical simulation of condensation in a rectangular minichannel using VOF model [C]//ProceedingsoftheASME2012InternationalMechanicalEngineeringCongress&Exposition. Houston, Texas, USA, 2012: 85422-01-85422-07.

Numerical simulation of condensation heat transfer inside single square minichannel

Li Panpan Chen Zhenqian

(School of Energy and Environment, Southeast University, Nanjing 210096, China)

A volume of fluid (VOF) model was adopted to simulate the condensation of R134a in a horizontal single square minichannel with side length as 1 mm and channel length as 50 mm. The effects of inlet mass velocity, heat flux, and mode of imposing cooling capacity on minichannel condensation were studied. The results show that condensation first appears at the corner of the channel, and then it is stretched on the effects of the surface tension until the whole channel boundary is covered. The local heat transfer coefficient decreases drastically on the head stream to a lowest value and then keeps almost unvaried. In the constant heat flux, the increased inlet vapor mass velocity has an impact on the enhancement of gas-liquid interface shear stress. As a result, the condensation film is thinner and the heat transfer coefficient gets higher value. While at constant inlet mass velocity, the effect of heat flux on heat transfer is unobvious. Compared with that exerted on a certain boundary intensively, cooling capacity uniformly exerted on channel boundary is more beneficial to condensation heat transfer with a uniform distribution of the liquid film.

square minichannel; two-phase flow; condensation heat transfer; numerical simulation; VOF (volume of fluid) model

10.3969/j.issn.1001-0505.2016.04.015

2015-11-20.作者简介: 李盼盼(1989—)，女，博士生；陈振乾(联系人)，男，博士，教授，博士生导师，zqchen@seu.edu.cn.

总装载人航天工程空间应用系统“天舟一号”货运飞船“两相系统实验平台关键技术研究”资助项目(TZYY08001).

10.3969/j.issn.1001-0505.2016.04.015.

TK124

A

1001-0505(2016)04-0763-07

引用本文: 李盼盼,陈振乾.正方形小通道内的冷凝换热数值模拟[J].东南大学学报(自然科学版),2016,46(4):763-769.