水平管外降膜流动液膜厚度数值模拟

邱庆刚，吕 多

（大连理工大学 能源与动力学院，辽宁 大连 116024）

0 引 言

水平管外液体的降膜流动具有传热温差小、传热效率高、低温传热性能优良等优点，因此被广泛应用于海水淡化、制冷系统、石油化工等行业中，在节能、环保方面有重大的意义［1－2］.国内外许多学者对降膜流动过程进行了卓有成效的研究.Ribatski等［3］的研究表明在水平管外液体降膜流动充分发展的层流区域，液膜与管壁的换热是以导热过程为主的对流换热，液膜厚度与表面传热系数成反比关系，因此，研究液膜厚度对于了解液膜流动机理以及提高换热系数有重大意义.Nusselt［4］通过理论分析给出了液膜厚度的经典表达式ρl（ρl－ρg）gsinβ，Rogers［5］应用动量方程积分方法得出了液膜厚度与喷淋密度、黏度、管周向角度关系的表达式；Gstoehl等［6］应用光学测量法测得了沿管壁周向80°范围内的液膜厚度分布，并发现水平管外周向下半区域同理论预测值不符，液膜厚度偏小；许莉等［7］测量了无换热情况下管外液体成膜的平均厚度及膜的几率分布与流率、管径的关系；何茂刚等［8］通过实验测量得出了流量及布液高度对液膜厚度的影响规律；Harikrishnan等［9］通过二维数值计算方法计算了R134a水平管外降膜吸收过程的传热传质特性；Li等［10］通过实验方法测试了雷诺数在21.6～108.1的水平管束降膜传热特性；Hou等［11］利用实验测量了水平管壁周向不同角度的液膜厚度；Kostoglou等［12］利用统计学方法研究了湍流状态下液体降膜流动的表面形态及液膜重组问题.

上述的研究工作虽取得了各种成果，但仍有许多工作要做，比如水平管外降膜流动过程中结构及流动参数对液膜厚度与分布影响规律的相关文献就较少.本文在前人研究基础上，首先建立二维水平管外降膜流动模型，根据Rogers对水平管外降膜流动的分析，把液膜分为发展区与充分发展区，然后利用VOF 方法模拟不同雷诺数状态下的水平管外降膜流动过程，通过精确读取管壁周向180°的液膜厚度，分析管间距、雷诺数等结构及流动参数对液膜厚度与分布的影响规律.

1 数值模拟方法

1.1 模型与基本参数

本文所研究的水平管外液体降膜流动模型参见图1.图2为水平管降膜蒸发器内部管束排列示意图，管束采用等边三角形排列，管束上方设有布液器.本文主要研究雷诺数及管间距对液膜的影响，故选取垂直于管排方向一横截面作为计算区域，为消除布液方式对液膜的影响，取第2根管读取周向液膜厚度；为消除管排方式对液膜的影响，管间距s为变化值，布液器孔径为3mm，布液高度为3mm，管径R为19mm，数值计算区域为图中阴影部分.

图1 水平管外液体降膜流动示意图Fig.1 Diagram of liquid falling film flow out a horizontal tube

图2 管排方式及计算区域Fig.2 Arrangement of tube bundles and computational domain

1.2 流体物性与基本假设

液体降膜流动雷诺数定义为Re＝4Γ／μ，其中Γ为液体的喷淋密度，为管径方向单位长度上流体的质量流量，μ为液体的动力黏度.本文所计算的Re范围均小于1 000，因此可认为本文所涉及的液膜流动均在稳定层流的范围内.

假设蒸发器内气相全部为空气，液相为饱和水，基本物性参数见表1.

表1 流体热物性参数Tab.1 Thermo physical properties parameters of fluid

1.3 网格及边界条件

本文采用精度较高的四边形网格对计算区域进行网格划分，同时在气液交界面及液相区进行了网格加密.计算域的网格划分如图3所示，计算中比较了网格数分别为2 127、5 083、12 149的计算模型，结果差异较小，最终选取网格数为2 127的计算模型，在达到精度要求的前提下节省了计算时间；同时比较了时间步长分别为0.05、0.10、0.20ms的计算结果，经比较0.10ms的时间步长能较为准确且高效率地得出计算结果，残差收敛于10－5.

图3 网格划分及局部网格放大图Fig.3 Sketch of grid and local magnified grid

图中速度入口为不同喷淋密度下的流体入口流速，压力条件为1.013×105Pa，其余为对称边界.由于模拟条件无换热过程，边界条件中无温度及热流密度边界.

1.4 控制方程及数值计算方法

微元控制体质量守恒方程为

式中：ρ为流体密度；t为时间；x、y、z分别表示坐标方向；u、v、w分别为速度矢量u在x、y、z方向的分量.由于ρ为常数，模型简化为二维，故该式可转化为＝0.

本文中流体为牛顿流体，N－S方程可简化为

式中：p为压强；μ为流体的动力黏度；ρgx、ρgy为微元体上的体积力，即液膜重力在各方向的分力.

本文采用VOF方法追踪两相流自由界面的流动情况，VOF方法的基本思想是引入某一相的体积分数αi，对于某一计算单元，αi＝0说明第i相在该单元内是空的；αi＝1说明第i相充满该单元；0＜αi＜1说明该单元包含了多相流界面.αi的控制方程为

本文中设定空气为第一相，水为第二相.计算单元中密度及黏度值为体积分数加权平均值，由下式决定：

2 结果及分析

通过数值计算得出不同参数下的气液两相分布图与速度云图，并读取管壁不同位置的液膜厚度值.图4为管间距s＝6.4mm，Re＝574时的气液两相分布图与速度云图.从气液两相分布图中可以看出此时液膜流动状态较稳定，液膜铺展均匀连续；同时从速度云图中也可看出管壁下方135°区域液膜流速相对较大.

图4 s＝6.4mm，Re＝574时的气液两相分布图与速度云图Fig.4 Vapor－liquid phase and velocity distribution chart while s＝6.4mm，Re ＝574

图5 液膜厚度随角度变化图Fig.5 Film thickness versus angle

本文在Re为574、744 时计算了s分别为6.4、9.5、19.4 mm 条件下管外液体降膜流动情况，读取了管壁周向180°范围内的液膜厚度值（d），并同文献［4］的理论值及文献［6］的实验值进行对比，如图5所示.从图中可以看出，当Re为574时，模拟值变化趋势同实验值比较吻合，液膜平均厚度随管间距s的增加而减小，其原因可能是随着布液高度的增大，液体流速变大，液膜形成初期滴状流变得更加不稳定，产生的液滴喷溅对膜流量存在较大影响，这也验证了液膜厚度与液体流速成反比关系；同一条件下沿管壁周向液膜分布趋势同实验值较吻合，同文献［4］的理论值在90°以上范围内存在一定偏差，液膜厚度总体上在沿管壁周向180°范围内先减小后增大，趋势基本符合文献［4］的抛物线理论值，但是在管壁下端，液膜厚度增长的趋势较小，其原因可能是管壁上方90°区域属于液膜层流充分发展区，液膜成膜较好，但波动较大；而在管壁下方90°区域内，体积力占主导作用，液膜在黏滞力、表面张力与体积力共同作用下流动趋于稳定，其中黏滞力使流体能均匀地附着在管壁上，表面张力维持液膜的形态，体积力保持液膜的流动，这可能是由文献［4］分析中忽略液体动量的变化这一条件所造成的.同时，20°～45°液膜厚度波动相对较大，在45°～90°液膜厚度变化较为平稳，这一点同文献［5］把管外液体降膜流动分为发展区与充分发展区相符合，可近似地认为在该种情况下45°以上区域为液膜流动的发展区，液膜厚度波动较大；之后区域为充分发展区，液膜厚度平稳减小.

当Re为744时，液膜厚度总体上增大，同时波动也增大，随着雷诺数及管间距的增大，气液两相区边界液膜厚度波动较大.在此雷诺数下，当s＝19.4mm 时，液膜厚度沿管壁方向呈现较大的波动趋势，模拟值普遍比实验值以及理论推导值大，其原因可能主要有两点：一是本文中采用的是二维简化模型，计算时忽略了沿管长方向液膜间的相互影响；二是计算中边界条件较为理想，并未估计扰动值，导致液体成膜情况较好，与实验值存在一定差异.

Re为100～1 000，管壁周向角45°、90°、135°处液膜厚度随Re变化如图6所示.

图6 液膜厚度随雷诺数变化Fig.6 Film thickness versus Re

从图中可以看出，管外壁周向某一角度处液膜厚度随雷诺数的增大而增大，增长的速率近似为线性，管壁周向角为45°时，液膜厚度随管间距的增大而减小，液膜增长曲线较平滑，其原因可能是液膜流动处于充分发展区，流动稳定；周向角为90°时，液膜厚度增长曲线开始出现波动，管间距较大时波动较为明显，其原因可能是此时液膜流速增大，流动不稳定性增大；周向角为135°时，液膜厚度波动加剧，局部出现了“干区”，可见液膜厚度波动最为剧烈，此时，从速度云图中可以看出液膜流速较大.重力与表面张力共同作用，导致液膜厚度波动较大，因此，水平管降膜流动液膜最薄处应该出现在管壁周向90°～135°.

当管间距大于12mm，液膜雷诺数大于200时，管壁下部会出现“干区”，如图7所示，并且“干区”的范围随雷诺数及管间距的增大而增大.在水平管降膜蒸发器中，“干区”的出现会造成换热系数降低，局部过热引起设备的损坏，在工程中应当避免.可以通过在适当范围内减小液膜雷诺数及管间距来避免“干区”出现.

图7 水平管壁“干区”示意图Fig.7 The dryout region of horizontal tube

3 结 论

（1）在一定雷诺数下，水平管壁周向180°范围内液膜厚度的变化趋势为先减小后增大，同实验值及理论值吻合较好，证明了本文模型与算法的正确性和可靠性.

（2）随着雷诺数的增大，液膜厚度波动增大，管壁周向某一角度液膜厚度逐渐增大，增大速率近似线性.

（3）管壁周向某一角度的液膜厚度随着管间距的增大而减小.

（4）一定雷诺数与管间距条件下，管壁下部会出现“干区”，“干区”的范围随着雷诺数及管间距的增大而增大，工程实际中应在适当范围内尽量减小液膜雷诺数及管间距来避免“干区”的出现.

［1］ 王小飞，何茂刚，张 颖.水平管降膜蒸发器管外液体流动数值模拟［J］.工程热物理学报，2008，29（8）：1347－1350.WANG Xiao－fei，HE Mao－gang，ZHANG Ying.Numerical simulation of the liquid flowing outside the tube of the horizontal tube falling film evaporator ［J］.Journal of Engineering Thermophysics，2008，29 （8）：1347－1350.（in Chinese）

［2］ 杨 丽，王 文，白云飞，等.水平管降膜蒸发器传热优化研究［J］.工程热物理学报，2009，30（11）：1913－1916.YANG Li，WANG Wen，BAI Yun－fei，etal.Heat transfer optimization of the horizontal tube bundles in falling film evaporators ［J］.Journal of Engineering Thermophysics，2009，30（11）：1913－1916.（in Chinese）

［3］ Ribatski G，Jacobi A M.Falling－film evaporation on horizontal tubes — a critical review ［J］.International Journal of Refrigeration，2005，28（5）：635－653.

［4］ Nusselt W.Die oberflchenkondensation des wasserdampfes zeitschr ［J］.Verein Deutscher Ingenieure，1916，60（2）：541－546.

［5］ Rogers J T.Laminar falling film flow and heat transfer characteristics on horizontal tubes［J］.The Canadian Journal of Chemical Engineering，1981，59（2）：213－222.

［6］ Gstoehl D，Roques J F，Crisinel P，etal.Measurement of falling film thickness around a horizontal tube using a laser measurement technique［J］.Heat Transfer Engineering，2004，25（8）：28－34.

［7］ 许 莉 王世昌 王宇新 等.水平管外壁液膜流动状态及其对传热的影响 ［J］.化工学报，2002，53（6）：555－559.XU Li，WANG Shi－chang，WANG Yu－xin，etal.Flowing state of liquid films over horizontal tubes and its influences on heat－transfer characteristics［J］.Journal of Chemical Industry and Engineering，2002，53（6）：555－559.（in Chinese）

［8］ 何茂刚，王小飞，张 颖.制冷用水平管降膜蒸发器的研究进展及新技术［J］.化工学报，2008，59（2）：23－28.HE Mao－gang，WANG Xiao－fei，ZHANG Ying.Review of prior research and new technology for horizontal－tube falling－film evaporator used in refrigeration［J］.Journal of Chemical Industry and Engineering，2008，59（2）：23－28.（in Chinese）

［9］ Harikrishnan L，Shaligram T，Maiya M P.Numerical study of heat and mass transfer characteristics on a falling film horizontal tubular absorber for R－134a－DMAC ［J］.International Journal of Thermal Sciences，2011，50（2）：149－159.

［10］ LI Wei，WU Xiao－yu，ZHONG Luo，etal.Heat transfer characteristics of falling film evaporation on horizontal tube arrays［J］.International Journal of Heat and Mass Transfer，2011，54（9－10）：1986－1993.

［11］ HOU Hao，BI Qin－cheng，MA Hong，etal.Distribution characteristics of falling film thickness around a horizontal tube［J］.Desalination，2012，285：393－398.

［12］ Kostoglou M，Samaras K，Karapantsios T D.Reconstruction of film thickness time traces for wavy turbulent free falling films［J］.International Journal of Multiphase Flow，2010，36（3）：184－192.