刘华^1,2, 龚路远¹, 沈胜强¹, 陈学¹
1. 大连理工大学 海洋能源利用与节能教育部重点实验室， 辽宁 大连 116024
2. 燕山大学 车辆与能源学院， 河北 秦皇岛 066004
收稿日期: 2014-11-20
基金项目: 国家自然科学基金重点项目(51336001)； 国家科技支撑计划项目(2014BAB09B00); 河北省高等学校科学技术研究项目(QN2016158).
作者简介: 刘 华(1984-),男,山东泰安人,大连理工大学博士研究生，燕山大学讲师;
沈胜强(1961-),男,河北保定人,大连理工大学教授。

摘要: 为研究水平管外液体和气体相互作用下的两相流动特性，选择蒸发器中广泛应用的转角正方形排列管束为物理模型，采用流体体积函数(VOF)方法追踪气液界面，提出液体在管间以液滴形式存在的模型假设，结合管间空隙率数据来初始化水平管外液膜厚度和液滴直径，模拟蒸汽在管外液膜和管间液滴作用下的流动过程，分析气液两相的压力场和速度场.结果表明：小喷淋密度下，进出口压降计算值和实验吻合良好；在计算域内，下部区域压力值高于上部区域，且最小压力分布在液滴附近的右下侧区域；压力分布的不均会造成液滴在下落过程中的变形.
关键词：VOF模型流动特性空隙率液膜液滴
Characteristic Research on Steam Flowing Across Tubes Under the Interaction of Liquid Film and Liquid Droplet
LIU Hua^1,2, GONG Lu-yuan¹, SHEN Sheng-qiang¹, CHEN Xue¹
1. Key Laboratory of Ocean Energy Utilization and Energy Conservation of Ministry of Education, Dalian University of Technology, Dalian 116024, China
2. College of Vehicles and Energy, Yanshan University, Qinhuangdao 066004, China
Corresponding author: SHEN Sheng-qiang, E-mail: zzbshen@dlut.edu.cn
Abstract: To investigate the two-phase flow characteristics under the interaction of liquid and gas over horizontal tubes, the rotated-square arranged tube bank was chosen and the VOF method was employed to propose an intertube droplet model. The thickness of liquid film and the droplet diameter was initialized according to the void fraction. The steam flowing process under the interaction of liquid film and droplet was simulated. The analyses of velocity field and pressure field were also illustrated. The results indicate that the calculated value of pressure drop coincides well with the experimental data at the small spray density. In the computational zone, the pressure in the bottom area is higher than that in upper area, and the minimum pressure distribution is in the lower right side near the droplet. The deformation of the liquid droplet in falling process is caused by the uneven distribution of pressure.
Key Words: VOF modelflow characteristicvoid fractionliquid filmliquid droplet
水平管降膜蒸发器在小温差和小喷淋密度下具有高传热系数的优点使其在制冷、化工和海水淡化等领域得以广泛应用^[1-2].水平管降膜蒸发器管束区内饱和液体以给定的流速喷淋到管子顶部并沿管子外表面向下流动，在降膜流动过程中吸收管内蒸汽凝结释放的热量而产生二次蒸汽^[3-4].因为两相流体黏性、界面表面力及气液界面的剪切力的存在，管外二次蒸汽的流动和液体降膜流动之间的相互作用必然对两相流动特性有所影响.描述管束间两相流动的一个重要参数是压降^[5].Ilyushchenko等^[6]给出了特定管束结构下蒸汽冲刷降膜管束时的压降系数关联式，但其与现在工业中广泛应用的管径、管束结构、相对管间距都有差异，甚至气相工质都不同；国内也有学者就其压降特性开展了相应实验研究^[7].但因实验研究耗资巨大，可查相关研究文献较少.
近代发展起来的计算流体力学(CFD)方法能直观体现两相流场，有助于两相流动动力学机理的研究.流体体积函数(VOF)方法^[8]特别适合计算有明显相界面的流场，适用于管外降膜流动的气液两相分层流动过程.本文基于ANSYS14.0计算平台，采用VOF模型并通过有限体积法数值模拟水平管管外降膜流动液体及蒸汽水平方向流动时的流动特性，进一步揭示降膜流动管束间二次蒸汽和液体流动的动力学特性机理.
1 模型与方法1.1 物理模型管外径25.4 mm的转角正方形排列管束当相对管间距为1.3D时，垂直和水平方方向的两排管心距离为46.7 mm.如图 1所示，选取二维模型进行计算，以便减小计算量.蒸汽绕圆管流动的雷诺数小于2 000，可认为是层流状态.本文做如下假设：流体不可压缩；流体物性均为相应饱和温度下的物性；流体工质液相为50~70 ℃的饱和水，气相为50~70 ℃的饱和水蒸气.
图 1(Fig. 1)

图 1 物理模型Fig.1 Physical model

1.2 数学模型质量连续方程：

(1)

动量方程：

(2)

对于存在自由表面的降膜流动，尤其是在水平蒸汽流动作用下，影响其运动形式的重要因素包括气液两相间的摩擦力和表面张力.式(2)中，动量源项F可表示成气液界面摩擦力源项F_LG与表面张力源项F_VOL之和.即

(3)

气液两相界面摩擦力可表示为

(4)

式中，f_LG表示气液界面摩擦系数：

(5)

式中，f_SC表示蒸汽流过与计算域相同直径圆管时的摩擦系数^[9].
表面张力表示为

(6)

式中，界面曲率κ由自由表面处单位曲面法向量的散度表示如下^[10]：

(7)

(8)

体积分率连续方程：VOF模型通过追踪计算单元内的体积分率 α_i的分布而得到相界面位置.式(9)和(10) 用于计算相界面的分布.

(9)

(10)

控制方程中的物理性质由系统各相共同确定.
1.3 计算设置及模型求解1) 边界条件.如图 2所示，对模型的边界进行定义，相关的数学模型如下.
图 2(Fig. 2)

图 2 边界条件及网格图Fig.2 Boundary conditions and grids

无滑移边界：u=0.速度进口：u=u_Lin.
压力出口：p=p₀,.
对称边界：,.
2) 网格独立性.数值计算采用的是四边形结构化网格，网格尺寸如下：①管径外壁液膜分布区：最底层(壁面处)网格为0.01 mm并以1.1倍的比例因子向外扩展10层；②蒸汽流通区的网格尺寸比较了4种规格：0.05,0.1,0.2,0.3 mm.为验证网格无关性，分别选用4种网格进行计算，以进出口的压降数值为参考依据进行比较得到网格无关性结论：0.05，0.1，0.2 mm 三种网格尺寸下的计算数据(压差值)最大相对误差小于2.1%.本文综合考虑网格无关性和减小计算量等因素，择优选取网格数为0.1 mm的网格模型.
3) 气液两相组分设置.数值计算开始前，在圆管外壁液膜区及在上下两管垂直管间距上初始化设置添加液膜、液滴来体现不同喷淋密度的影响.图 3a为设置完成后的液膜和液滴的组分图，图中红色表示液体，蓝色表示气体；圆管外壁的一层红色区域代表液膜，圆管下垂线和上垂线上的红色圆球即为液滴.图 3b为实验中拍摄的不同喷淋密度下的降膜流型图.
图 3(Fig. 3)

图 3 液膜液滴的初始设置Fig.3 Initial install of liquid film and droplet (a)—组分设置；(b)—实验图片.

液膜厚度设置:学者们从理论分析、实验测量甚至数值模拟方面对液膜厚度的大小、液膜沿管周方向的分布及其影响因素等做了大量详实的工作.现成可用的液膜厚度计算公式有Nusselt理论分析公式、Roger积分公式及Hou等依据实验数据的改进公式^[11]. 为便于计算，作者采用液膜沿圆周方向厚度均匀分布的假设.
液滴大小设置:在不同喷淋密度工况下拍摄了大量的管间降膜流动照片，有足够的图片展示了与所选物理模型相同的两管之间的降膜流动，典型的流型见图 3b.
众所周知，依据图像分析原理可以将同一张图片中的不同影像所占像素区分开.作者采用Labview(NI,Austin,USA)中的Vision Assistant工具将液体和气体的像素点进行处理后加以区分，进一步处理可获得液体和气体分别所拥有的像素点的数量.选取上下两水平管间面积作处理区域，进而求出该降膜流动区域的空隙率ε(气体区域像素点占整个区域像素点的比值).
基于图 4所示空隙率的数据，将各个工况下( Γ=0.02,0.04,0.06,0.08 kg/(m·s) )的空隙率值求期望值可得该喷淋密度下的平均值ε_Γ.由ε_Γ可计算上下两管间的液体所占比例，进而能求出该模型下液滴的等效直径：
图 4(Fig. 4)

图 4 管间空隙率图Fig.4 Diagram of the intertube void fraction

(11)

式中，S表示上下相邻两管的最近距离.Γ=0.02 kg/(m·s)时，管间流态为滴状流，液滴直径即为

(12)

相比滴状流的设置，柱状流态下的区别就是增加液滴的数量.如图 3a，下垂线上主液滴的直径依然按Γ=0.02 kg/(m·s)时赋值：d₁=d_Γ=0.02；上垂线上的另外一个液滴直径d₂根据下式计算：

(13)

4) 数值格式设置.离散时，时间项采用隐格式；对流项采用二阶迎风格式；压力项采用BODY算法；压力-速度耦合方程选用PISO方法.各个输运方程计算精度为10^-3.时间步长为10^-5.气液界面追踪方法选用精度较高的Geo-Reconstruct格式.计算采用的是非稳态计算方法，因此计算结果和时间步长有一定的关系.在保证液滴不被吹散的前提下，选取5 000个时间步长统一比较.
1.4 模型验证液滴模型及初始化气液两相组分的设置，包含了一定的假设和等效简化.因此，数学模型、两相组分初始设置都关系到计算结果的准确性，通过实验结果来验证.
文献[1]获得了转角正方形排列管束前后的压降值.除以管列数(N)后可求得单列管压降值，选其作为参考值来验证计算结果的准确性.验证的数据工况为：t=70 ℃，Γ 分别为0.02和0.04 kg/(m·s)，见图 5和图 6.
图 5(Fig. 5)

图 5 喷淋密度为0.02 kg/(m·s)时的数据验证Fig.5 Verification in Γ=0.02 kg/(m·s)

图 6(Fig. 6)

图 6 喷淋密度为0.04 kg/(m·s)时的数据验证Fig.6 Verification in Γ=0.04 kg/(m·s)

由图 5可知，喷淋密度Γ=0.02 kg/(m(s)时压降计算值和实验值吻合非常好，最大误差在±3.5%以内，证实计算结果可信.计算采用的液滴模型和实际滴状流态很是相近，可以通过计算值来补充实验数据，为研究蒸汽横掠降膜管束时的流动特性丰富数据储备.
由图 6可知，喷淋密度Γ=0.04 kg/(m·s)，在蒸汽流速较小工况时，压降值与实验值基本一致.但从8 m/s开始，计算值和实验值出现规律性的偏移和分离，具体表现为计算值小于实验值.这是因为Γ =0.04 kg/(m·s)时管间降膜流动的存在形式多为液柱或尚未发展完全的液柱.当蒸汽流速较高时，液柱受到横向吹扫的作用更加明显，因此实验状态下测得的压降值要大.液滴相比液柱对横向气流的阻碍要小一些，蒸汽流速较大时，数据的偏移已经显现.因此，喷淋密度Γ=0.04 kg/(m·s)工况下，只能选取较低流速下的计算值补充实验数据.
通过验证可知：液滴模型的有效计算在Γ=0.02 kg/(m·s)下能够成立，另外还包括较低蒸汽流速下的Γ=0.04 kg/(m·s)工况.
2 结果与讨论2.1 两相组分及压力分布图 7为Γ=0.02和Γ=0.04 kg/(m·s)下的两相组分图和压力分布图.从两相组分图可知小液滴更容易受气流的影响，Γ=0.04 kg/(m·s)时，小液滴的偏移距离比主液滴的偏移距离要大，见图 7b，v=4 m/s时小液滴的偏移距离约是大液滴的2.37倍. 从压力分布图可知：1)液滴前后(水平方向)的压力变化很大，因为液滴位置和边界层分离的缘故，最小压力分布在液滴的右下侧；2)中垂线之前的下部区域要高于上部区域，这是主液滴置于管底下部所致.压力分布的不均会造成液滴在下落和偏移过程中的变形.
图 7(Fig. 7)

图 7 相组分及压力分布图Fig.7 Contours of phase and pressure (a)—Γ=0.02 kg/(m·s);(b)—Γ=0.04 kg/(m·s).

2.2 速度流场图 8所示为两种工况下的速度矢量图.由图可知，液滴(甚至两个液滴)的存在使得圆管上下的速度分布并不均匀.数值计算中，喷淋密度的不同是由液滴数量和液滴尺寸来体现.在相同蒸汽入口速度下，提取Γ=0.02和Γ=0.04 kg/(m·s)两种喷淋密度下的速度计算值比较发现，最大值相差34.7%.另外，图 8右侧部分显示的是速度大于1/2倍速度最大值的“高流速”区域(相对高流速区).液滴的存在使得高流速区是一个沿管壁圆周方向的狭长区域.
图 8(Fig. 8)

图 8 速度矢量图Fig.8 Vectors of fluid velocity (a)—Γ=0.02 kg/(m·s);(b)—Γ=0.04 kg/(m·s).

相比实验测量，计算能大大节省所需时间、经济及人力成本.该计算模型的提出为进一步分析蒸汽和液体两相的流动特性提供了一种新的方法，为揭示降膜流动管束间二次蒸汽和液体流动的动力学特性机理奠定了一定的研究基础.
3 结论1) 为分析水平管外液体和气体的相互作用，提出了管外液滴模型，经过实验数据验证可知：小喷淋密度下，采用本文提出的液滴模型模拟计算结果与实验值吻合较好.
2) 在计算域内，下部区域压力值高于上部区域，且最小压力分布在液滴的右下侧；压力分布的不均会造成液滴在下落和偏移过程中的变形.
3) 在相同蒸汽入口速度下，Γ=0.02和Γ=0.04 kg/(m·s)两种喷淋密度，计算域内的速度最大值相差34.7%.
参考文献

[1]	Liu H, Shen S Q, Gong L Y, et al. Shell-side two-phase pressure drop and evaporation temperature drop on falling film evaporation in a rotated square bundle[J].Applied Thermal Engineering, 2014, 69(1/2): 214–220.(0)
[2]	Yang L, Shen S. Experimental study of falling film evaporation heat transfer outside horizontal tubes[J].Desalination, 2008, 220(1/2/3): 654–660.(0)
[3]	刘华, 沈胜强, 龚路远, 等. 水平管降膜蒸发器温度损失的计算与分析[J].西安交通大学学报, 2014, 48(4): 90–94. ( Liu Hua, Shen Sheng-qiang, Gong Lu-yuan, et al. Evaporation temperature loss evaluation in horizontal tube falling film evaporator[J].Journal of Xi’an Jiaotong Universtiy, 2014, 48(4): 90–94.)(0)
[4]	Shen S Q, Liu H, Gong L Y, et al. Thermal analysis of heat transfer performance in a horizontal tube bundle[J].Desalination and Water Treatment, 2015, 54(7): 1809–1818.(0)
[5]	Ribatski G, Jacobi A M. Falling-film evaporation on horizontal tubes—critical review[J].International Journal of Refrigeration, 2005, 28(5): 635–653.(0)
[6]	Ilyushchenko V V, Podbereznyi V L, Golub S I. Influence of spray density of coefficient of gas dynamic resistance for tube bundles of different configuration[J].Desalination, 1989, 74(1/2/3): 383–389.(0)
[7]	沈胜强, 刘华, 冯寅, 等. 真空条件下蒸汽横掠水平降膜管束的流动阻力[J].化工学报, 2013, 64(3): 886–890. ( Shen Sheng-qiang, Liu Hua, Feng Yin, et al. Steam flow resistance across horizontal tube bundle under vacuum condtion[J].CIESC Journal, 2013, 64(3): 886–890.)(0)
[8]	Hirt C, Nichols B. Volume of fluid (VOF) method for the dynamics of free boundaries[J].Journal of Computational Physics, 1981, 39: 201–225.(0)
[9]	谷芳, 刘春江, 余黎明, 等. 气-液两相降膜流动及传质过程的CFD研究[J].高校化学工程学报, 2005, 19(4): 438–444. ( Gu Fang, Liu Chun-jiang, Yu Li-ming, et al. The CFD simulation of mass-transfer process in falling film with countercurrent two-phase flow[J].Journal of Chemical Engineering of Chinese Universities, 2005, 19(4): 438–444.)(0)
[10]	Brackbill J U, Kohe D B, Zemach C. A continuum method for modeling surface tension[J].Journal of Computational Physics, 1992, 100(2): 335–354.(0)
[11]	Hou H, Bi Q, Ma H, et al. Distribution characteristics of falling film thickness around a horizontal tube[J].Desalination, 2012, 285(1/2/3): 393–398.(0)