水平波纹管外降膜蒸发的液膜流动与换热特性分析

王迎慧, 邢聪骢, 刘建停

(江苏大学 能源与动力工程学院， 江苏 镇江 212013)

水平管外降膜蒸发是一种高效的相变换热方法，多见于制冷热泵、石油化工、海水淡化等工业领域.相比满液式蒸发器，降膜式蒸发器的工质用量少，流动阻力小，换热系数高，换热面更为紧凑[1-3].研究发现，异形管管外的降膜蒸发可进一步强化换热效果.文献[4]以水为工质的试验和数值模拟研究发现，不同离心率椭圆管外的平均换热系数均比圆管提高20%～22%.文献[5]的数值模拟结果显示，水在水滴形管、卵形管外的液膜厚度分布不同于圆管，液膜分布的变化可以起到换热强化的作用.文献[6-7]基于流体体积(vo-lume of fluid，VOF)模型，构建水在椭圆管、圆管外降膜蒸发的三维模型，发现二者在横截面圆周方向上的液膜流速、液膜厚度存在差异.文献[8]试验研究发现，圆管外局部换热系数随周向角的增大而减小.文献[9]试验研究表明，水在圆管外降膜蒸发的强化换热源于管外液膜内较多的气泡.在制冷工质方面，文献[10-12]以R134a(C2H2F4)、R123(CF3CHCl2)和R290(CH3CH2CH3)等制冷剂为工质的试验研究发现，圆管外平均换热系数随降膜雷诺数、热流密度的增大而增大.文献[13]通过数值模拟R410A(CF3CFH2)在圆管外的降膜蒸发过程，发现圆管外液膜厚度随周向角呈先减小后增大的趋势，在周向角为110°处达到最小值，局部换热系数则随周向角逐渐减小.

现有研究多数关注水在光管外的降膜蒸发行为，关于制冷工质在管外降膜蒸发的报道则较少.波纹管作为工程中常见的强化换热管之一，其沿轴向的截面呈周期性变化，具有同时强化管内外流体的对流传热作用.与椭圆管、水滴形管、卵形管等异形管不同，波纹管不仅涉及换热管周向上流动及换热特性的变化对液态工质成膜相关特性的影响，且沿轴向流动及换热特性的变化也会影响液膜流动与换热.为此，笔者针对低温制冷剂R134a在换热管(水平波纹管和光管)外发生降膜蒸发的差异，重点关注液膜沿轴向、周向的流动与换热情况，研究液膜在管外的流动、液膜厚度及分布、液膜内气泡成核及分布等，探讨波纹管外换热性能的提升效果.

1.1 物理模型与参数设定

波纹管一般是经过成型工艺制成，波纹管模型如图1所示.水平管外的降膜蒸发一般是指液态工质自布液孔垂直降落于换热管外，在管外壁面上铺展形成较薄的液膜，吸收管内流体的热量，从而发生相变的过程.

图1 波纹管模型

鉴于波纹管几何结构的周期性和对称性特征，为节省计算资源，取波纹管沿轴向3个周期的长度为计算区域.波纹管计算区域的三维模型和波纹管计算区域尺寸分别如图2、3所示，图中D为波纹管直径(取值19.00 mm)，L为计算域总长(取值60.00 mm)，vin为入口速度，λp为波距(取值20.00 mm)，A为波幅(取值1.00 mm)，H为布液孔至波纹管中心垂直距离(取值20.25 mm)，θ为周向角.

图2、3中，液态R134a从波纹管中心上方20.25 mm处的布液孔(φ1.50 mm)以vin=1 m·s-1的流速垂直下落至换热管外壁.设定管内为冷冻水，管壁热流密度qw=2×104W·m-2，以R134a流出布液孔的时间设定为t=0 ms.为便于说明，讨论轴向坐标x分别为20.0(波谷处)、25.0和30.0 mm(波峰处)时截面液膜在周向上的相关特性.

图2 波纹管计算区域三维模型

图3 波纹管计算区域尺寸(单位：
mm)

图2中，定义布液孔为速度入口，入口温度为Tin=274.95 K.换热管材料为铜，管壁为恒热流边界条件.管外工作压力p0=0.314 62 MPa，对应于R134a饱和状态的温度Ts=275.15 K.换热管下方底面设为压力出口，pout=p0=0.314 62 MPa，Tout=275.25 K，取气液固之间接触角为15°.考虑计算区域的不规则性以及工质在管外流动成膜的特点，计算区域分块划分网格，并对换热管外邻近区域网格做加密处理.经网格无关性检验，网格总数约为7×106个.波纹管外计算区域的网格划分见图4.

图4 波纹管外计算区域的网格划分

1.2 数学模型

工质在换热管外的降膜蒸发涉及气液两相流动传热，可基于多相流模型中的VOF模型构建数学模型.由于VOF模型本身没有相变模块，并不能反映工质因相变而发生的传热和传质过程.为此，需要自行开发、编译用户自定义函数(user-defined functions，UDF)，以确定气液两相间传递的质量源项和能量源项.考虑波纹管结构会加剧管外流体的湍动强度，结合液膜流动的雷诺数(Re=4Γ/μl=875，其中Γ为喷淋密度，μl为动力黏度)，管外液膜流动可按湍流处理.为反映近壁区流体的剪切输运效应，此处选取SST(shear-stress transport)k-ω湍流模型[6].换热管外降膜蒸发的控制方程组如下：

连续性方程为

(1)

(2)

动量方程为

(3)

k方程为

(4)

ω方程为

(5)

能量方程为

(6)

式中：αg和αl分别为气、液相体积分数，每个计算单元内，αg+αl=1；
ρg、ρl和ρ分别为气、液相密度和气液相平均密度,kg·m-3；
v为液膜平均速度，m·s-1；
v′为液膜脉动速度，m·s-1；
S为液相向气相转移的质量源项，kg·(m3·s)-1；
μ为动力黏度，Pa·s；
p为液膜平均压力，Pa；
Fvol为液体张力产生的体积力，N·m-3；
k为湍动能；
ω为特定耗散率；
Gk和Gω分别为湍动能和特定耗散率的产生项；
Γk和Γω分别为湍动能和特定耗散率的扩散系数；
Yk和Yω分别为湍动能和特定耗散率因湍流产生的损耗；
Dω为正交扩散项；
E为广义内能，J·kg-1；
λ为平均导热系数，W·(m·K)-1；
T为平均温度，K；
Q为能量源项，W·m-3.该控制方程组中的μ、ρ、T和E等分别为按气、液两相体积分数进行加权平均处理后的参数值.其中式(1)和(2)中的质量源项S、式(6)中的能量源项Q由UDF实时赋值.

动量方程中，因表面张力产生的体积力Fvol为

(7)

式中：σ为气液两相界面的张力系数，N·m-1；
κl为气液界面的液相曲率，m-1，满足

(8)

离散通用方程(式(1)-(6))时，压力速度耦合采用具有分裂算子的压力隐式(pressure-implicit with splitting of operators，PISO)算法，对流项和扩散项分别采用二阶迎风格式和中心差分格式.

2.1 数值计算结果验证

为验证数值计算的准确性和可靠性，图5为水在光管外发生降膜蒸发时局部换热系数h、R134a在光管外发生降膜蒸发时液膜厚度δ随周向角θ变化的关系曲线[8].由图5可知，数值计算与试验结果之间的平均误差小于10.00%.

图5 局部换热系数、液膜厚度随周向角变化的关系曲线

2.2 波纹管外的液膜流动

图6为波纹管和光管外的液膜沿轴向和周向的流动速度变化云图(取αl=0.9等值面上的速度)，其中v为流速.由图6a可知：在轴向上，液柱从波峰截面向两侧铺展，相邻2个单元液膜在波谷截面处汇聚，形成液膜较厚的液膜聚集区；
在周向上，液膜在重力和表面张力的主导作用下沿周向铺展，最终当t=170.0 ms时在换热管底部形成新液柱.由图6可知：在沿轴向铺展的过程中，由于波纹管外液膜受重力沿轮廓线切向分力的作用，波纹管外液膜可以在较短时间内铺展至波谷截面(对应t=23.0 ms时的I处)，相同的轴向长度，光管外液膜铺展则需要25.0 ms；
在周向上，由于不同截面处液膜沿周向的铺展速度不同，导致换热管外液膜形成不同的前端轮廓.t=36.0 ms时，液膜在波纹管铺展至周向角θ=60°附近，各处截面上的铺展速度相对均衡，液膜前端呈“一”字形轮廓；
光管外对应液柱中心截面的铺展速度较快，远离液柱的截面处速度则较慢，且距离越远，周向速度越低，t=36.0 ms时液膜前端呈现的“V”字形轮廓越清晰.

图6 换热管外不同时刻的液膜流动

由图6还可知，在t=56.0 ms时，波纹管外的液膜基本铺满换热管迎面区(θ=0°～90°)，光管外液膜则整体越过半圆(θ≥90°)，部分液膜已铺展至θ=120°附近.此外，液膜在换热管外脱离，形成的新液柱也存在差异.对应t=170.0 ms时刻，液膜已完成周向铺展，光管底面脱离区形成多个分散的小液柱，而波纹管外的脱离区则形成几股较大的新液柱，这对管束外的降膜蒸发十分有利.

2.3 波纹管外的液膜厚度

由于实际情况下换热管外液膜沿轴向、周向流动存在差异，二者的液膜厚度分布也会有所不同.取液膜在管外完全铺展的t=170.0 ms时刻，关注图3中x=20.0，25.0，30.0 mm截面的液膜铺展情况，取αl=0.9等值面上的液膜厚度.考虑液膜流动中液膜表面的扰动，此处取3次液膜厚度计算结果的平均值.可知换热管外在x=25.0，30.0 mm截面的液膜较薄，厚度分别为0.096～0.166 mm和 0.099～0.144 mm.而液膜聚集区x=20.0 mm截面的液膜较厚，厚度达到0.565～0.899 mm，且在液膜脱离时伴有不连续的液体溅落.由于换热管外液膜聚集区所占面积较小，因此主要讨论x=25.0，30.0 mm截面的液膜厚度及分布.图7分别为换热管在x=25.0，30.0 mm截面的液膜厚度随周向角变化的关系曲线.

图7 换热管在x=25.0，30.0 mm截面的液膜厚度随周向角变化的关系曲线

由图7可知：
在x=25.0，30.0 mm截面上，液膜厚度均随周向角增大呈先减小后增大的趋势；
θ=100°时的波纹管液膜最薄，厚度分别为0.096、0.099 mm；
θ=110°时的光管液膜最薄，厚度分别为0.106、0.110 mm.经计算可得：波纹管外的液膜平均厚度为0.120 mm，比光管薄了7.00%；
波纹管迎面区的液膜平均厚度为0.119 mm，比光管薄了13.14%；
背面区的液膜平均厚度为0.121 mm，比光管薄了0.24%.

2.4 波纹管外液膜内的气泡

图8分别为换热管外在不同时刻的气液两相体积分数分布云图，以及Ⅰ-Ⅴ处周向上相应截面的局部放大图.由图8可知：在t=33.6 ms时，波纹管外x=21.5 mm截面上θ=15°处出现首个成核点(图8a中Ⅰ处)，光管尚无成核点出现，经计算，光管外气泡成核点始于t=36.2 ms，较波纹管推迟2.6 ms；
在t=56.0 ms时，换热管外的液膜均铺展至半圆(θ=90°)，经统计，波纹管和光管外出现的成核点数分别为90、76个.由图8还可知：t=56.0 ms时，在x=25.0 mm截面上，波纹管和光管外均出现成核点，分别位于周向角为56°和71°处(图8中Ⅱ、Ⅳ处)；
在x=30.0 mm截面上，波纹管和光管外也出现了成核点，分别位于周向角为71°和82°处(图8中Ⅲ、Ⅴ处).在液膜完全铺展的t=170.0 ms时刻，波纹管和光管外的成核点数分别为286、248个，波纹管成核点数比光管多15.32%.其中，波纹管、光管迎面区成核点数分别为114、79个，波纹管比光管多44.30%；
波纹管、光管背面区成核点数大致相当，分别为172、169个.分析计算结果可知，在x=25.0，30.0 mm截面上，波纹管外的气泡成核点多见于θ=55°～135°，光管则多见于θ=80°～135°.可见，波纹管外液膜内出现气泡成核点较早，且成核点多于光管，周向分布区域较大.

图8 换热管外的气液两相体积分数分布云图

进一步计算可知，对应于t=33.6 ms时波纹管外首个成核点处(位于x=21.5 mm截面的θ=15°处)的液膜厚度最小，为0.126 mm.分析认为，管外液膜内的成核点与液膜厚度有关，液膜较薄区域成核点数较多，较厚区域则相对较少.

图9为x=30.0 mm截面θ=135°处(该点附近成核点较多)气泡成核、生长和脱离过程.由图9可知：t=86.0, 94.0 ms时，气泡成核点出现，并长大；
t=110.0 ms时，气泡脱离壁面；
t=118.0 ms时，溢出液膜.

图9 x=30.0 mm截面的θ=135°处气泡成核、生长和脱离过程

2.5 降膜蒸发的换热系数

为考察波纹管外降膜蒸发的换热强化效果，可以先计算出t=170.0 ms(此时液膜完全铺展)时换热管外的平均换热系数与不同截面处的局部换热系数.图10为换热管在x=25.0，30.0 mm截面上局部换热系数随周向角变化的关系曲线.

图10 换热管在x=25.0，30.0 mm截面上换热系数随周向角变化的关系曲线

由图10可知，这两处截面的局部换热系数均随周向角的增大而减小，换热管迎面区(θ<90°)的局部换热系数降幅显著，背面区(θ≥90°)降幅有所减缓.统计换热管外计算区域内的平均换热系数可知：波纹管外平均换热系数为2.60 kW·(m2·K)-1，比光管大10.17%；
在x=30.0 mm截面上，波纹管外平均换热系数为3.15 kW·(m2·K)-1，比光管大11.70%；
在x=25.0 mm截面上，波纹管外平均换热系数为2.64 kW·(m2·K)-1，比光管大10.46%.分析结果表明：距离波峰越近，波纹管外的强化换热效果越好；
波纹管外迎面区和背面区平均换热系数分别为3.92、1.11 kW·(m2·K)-1，迎面区换热系数为背面区的3.53倍.可见，迎面区为管外换热的主要贡献区.

1) 波纹管外液膜的流动状况与光管存在差异.轴向上波纹管液膜铺展速度较光管快，周向上波纹管外液膜铺展速度较光管慢.波纹管各处截面上的铺展速度相对均衡，而光管离波峰截面越近，铺展速度越快.波纹管外液膜脱离时能形成几股较大的新液柱.

2) 波纹管与光管外液膜厚度均随周向角的增大呈先减小后增大的趋势.波纹管外液膜平均厚度为0.120 mm，比光管薄7.00%.

3) 波纹管外液膜内的气泡成核点及分布与液膜厚度相关，且波纹管外的起始成核点早于光管.液膜较薄的区域成核点数较多，较厚的区域则相对较少.波纹管外的气泡成核点数多于光管，且分布较广(θ=55°～135°).

4) 波纹管外的平均换热系数比光管大了10.17%，其迎面区为管外换热的主要贡献区.