竖直管降膜蒸发器流动与传热过程的数值模拟

郝　丽，肖晓明，张吕鸿，姜　斌，2*，孙永利

(1.天津大学化工学院，天津 300072； 2. 精馏技术国家工程研究中心，天津 300072)

降膜蒸发器具有传热效率高，传热温差小，停留时间短等特点，特别适用于热敏性物料的蒸发，是目前被广泛使用的高效蒸发设备[1-4]。原料液由加热管的顶部加入，在自身重力作用下沿管内壁呈膜状下流，并被蒸发浓缩，气液混合物由加热管底部进入分离室，经气液分离后，将完成液由分离器的底部排出。为使溶液能在壁上均匀成膜，在每根加热管的顶部均需设置液体布膜器。由于影响其蒸发传热性能的主要因素是液膜的均匀分布状况和液体在管内结垢的情况，实际生产中对降膜蒸发器的研究主要围绕这两个方面进行。许多学者已对降膜蒸发器的传热性能和流体动力学进行了大量的实验研究和理论探索[5-8]。对降膜蒸发传热的理论研究最早是Nusselt[9]提出来的，在1916年对光滑层流下降液膜进行理论模拟，提出了液膜流动及其传热机理，得到液膜传热系数及液膜厚度。Moon-Hyun Chund等[10]采用改进的湍流降膜模型预测降膜过程的热传导系数，模型对比了有无界面剪切力时液膜被加热和冷却的情况，研究预测了较大的雷诺数范围内(Re=3 000～65 000)与流动相关的膜厚度和轴向传热系数。Feddaoui等[11]报道了用于研究降膜蒸发的液体沿竖直管蒸发冷却的数学模型，结果表明较好的液膜冷却系统具有较高的进液温度，较高的气体雷诺数或较低的液相流量。

在蒸发过程的相同热负荷作用下，给液不足的换热管可能会结垢、烧焦，甚至出现干斑或烧毁现象，而液膜过厚的换热管因传热量不足而不能充分传热，导致传热情况恶化。由此可见，对加热管液膜分布的研究至关重要。综合降膜蒸发器的模拟研究现状，对于降膜蒸发的温度-速度-传质耦合的膜运动过程，采用三维数值模拟的方法研究降膜蒸发过程的传质传热，并分析液膜厚度变化的研究还鲜有报道。因此，本研究在前人理论研究及实验的基础上，用FLUENT软件对三维降膜蒸发器加热管的两相流流场进行了稳态模拟，使用UDF编写混合物能量方程源项、两相质量传递方程源项，建立了降膜流动及传热模型，采用Mixture方法建立了气-液两相CFD模型，通过模拟得到管内流体的速度分布及降膜蒸发过程的液体分布，进而得出加热管内的液膜厚度。

1　竖直管降膜蒸发器CFD模拟

1.1　物理模型

由于每根加热管内的传热、传质情况类似，且原料液经过降膜蒸发器的液体分布装置均匀分布，所以单管的数值模拟即可基本代表蒸发器内所有加热管内的传热、传质情形。本研究的单管传热-流动的物理模型，受液槽上设计一个向下的圆管进口，保证液体进入导流管前能够在受液槽上有一定的持液量。加热管的计算区域如图1所示，圆管进口和螺旋槽结构如图2，其中Z轴为轴向，X、Y轴为径向，表1是螺旋型导流管和加热管的结构和几何参数。

表1　布膜器和加热管结构参数Table 1　Structural parameters of the film distributor and heating tube

图1　加热管Fig.1　Heating tube

图2　圆管进口和螺旋槽Fig.2　Inlet of cyclic tube and spiral flute

1.2　数学模型

降膜蒸发器加热管模拟涉及流体流动、传热和两相流传质，选择两相流模型为Mixture模型。流体流动状态由雷诺数决定，流体的雷诺数由进口质量流量计算得出，本研究中流体以层流形式流动，在此选用层流模型描述液膜的流动及传热过程。

混合模型[12]的方程主要包括混合物的连续方程，混合物的动量方程以及每相分别的能量方程。混合模型允许相之间互相贯穿，所以对一个控制容积的体积分数αq和αp可以是0和1之间的任意值，取决于相q和相p所占有的空间。混合模型使用了滑移速度的概念，允许相间以不同的速度运动。

混合模型的连续方程为：

(1)

(2)

混合模型的动量方程：

(3)

(4)

(5)

混合模型的能量方程：

(6)

(7)

这里keff是有效热传导率(k+kt，这里kt是紊流热传导率，根据使用的紊流模型定义)，方程(7)右边的第一项代表了由于传导造成的能量传递。SE包含了所有的体积热源。

FLUENT中的混合模型使用了代数滑移公式，代数滑移混合模型的基本假设是规定相对速度的代数关系，相之间的局部平衡应在短的空间长度标尺上达到。

混合物能量方程源项、两相质量传递方程源项使用UDF编写定义。其中源项定义式[13-14]为：T>TSAT时，即液相温度高于蒸发温度时，

(8)

(9)

SE=-RVhfg

(10)

其中hfg为界膜导热系数。

液相向气相的质量转移由UDF编写。

1.3　边界条件

1)进口条件：加热管进口定义为全液体为进口，液体速度为0.000 769 3 m/s，进料处温度设置为363 K；2)出口条件：加热管出口设定为质量出口边界条件；3)壁面条件：加热管壁面定义为无滑移，无渗透，0厚度，壁面温度保持为447 K，其余壁面定义为无滑移，无渗透及绝热。

1.4　加热管初始设置

液体设为第一相，气体设为第二相；模拟的气体由于在壁面附近为分散相，而在加热管管中心为连续相，汽液间曳力系数采用对称模型(symmetric)；操作压力定义为400 Pa，重力方向为Z轴的负方向，大小为9.8 m/s2，环境密度为0.03 kg/m3，与气相密度一致。

1.5　网格生成

在FLUENT的前处理软件Gambit中进行网格划分，本研究中需要重点考察加热管管内壁膜状流动的的流体，壁面网格的划分采用边界层网格。管内采用结构化网格，生成六面体网格，和螺旋槽联接处以及进口圆管采用混合网格画法，生成四面体网格，总网格数为50万。图3为网格划分示意图。图4和图5为局部放大的加热管横截面和纵截面网格示意图。

图3　圆管进口和螺旋槽网格分布Fig.3　Mesh of inlet of cyclic tube and spiral flute

图4　加热管横截面网格分布Fig.4　Mesh in cross section of heating tube

图5　加热管纵界面网格分布Fig.5　Mesh in longitudinal section of heating tube

2　结果与讨论

2.1　加热管内的流场分布

通过模拟得到加热管内液体流动的流线图和速度矢量图，分别如图6和图7所示。

图6　加热管内液体流线图Fig.6　Stream-line of heating tube

从图6中可以看出，液体在进入加热管之前，在受液槽(a)中有一定的持液量，如图6中深蓝色迹线所示，液体汇集到进口圆管(b)处，通过螺旋型导流管(c)形成一定的扰动，均匀分成3股流体，在重力的作用下流入加热管(d)中，液体在加热管中均匀流动，加热管中的流体迹线靠近壁面附近比较平缓的向下流动，中心处出现扰动较强烈的迹线。

图7　加热管内速度矢量图Fig.7　Velocity vectors in heating tube

由图7可知，流体在管中心的流速大于管壁附近的速度，且越靠近中心处流速越大，如图7中红色线条部分所示，中心流体发生了比较强烈的扰动，主要由于中心流体是液相汽化后生成的气体，随着气体的不断生成而出现扰动。加热管壁面由于无滑移边界条件，因此壁面速度接近于0，图7中蓝色显示，对应标尺的速度值接近于0。由于加热管较长，图7将整个加热管分成3部分显示，在入口段(a)流体速度在中心和壁面处的差异较小，且流体在中心速度较高的区域范围均匀的向下流动，主要由于刚进入加热管，流动平缓，而到达中部段(b)之后，出现了明显的波浪状的高速区，中心速度比入口段高出许多，湍动明显，说明汽化生成的气体主要出现在这一部分，到了出口段(c)速度沿主流方向慢慢下降，较平稳的流出加热管，其原因是沿着主流方向液膜不断蒸发，流量减小。

2.2　加热管内液膜的形成

Mixture模型中气液共存在每一个计算网格中，从液体体积分数定义的角度出发，计算区域内，液体体积分数大的地方视为液相区，液体体积分数小的地方视为气相区。沿液体流动的方向取6个横截面，做每个横截面的液体体积分率分布如图8所示。

图8　沿管长不同横截面上的液体体积分率分布(Z为所选截面距加热管入口端长度)Fig.8　Liquid volume fraction at different cross section (Z is the length from the inlet to the defined cross section)

图8a)中，液体在入口段由3个螺旋型导流管流入管内，在3个导流管与加热管管壁连接处液体体积分数最大(红色显示)，管中心和壁面其余处液体体积分数较小(蓝色显示)，代表液体主要聚集在这3个连接处，其余计算区域主要是气相。由于计算初始化中计算区域中并无气体，表明气相是从液相受热蒸发而来。图8b)和图8c)中，液体在管壁处的体积分数明显大于管中心处体积分数，代表液体开始大量附着在壁面上，管中心则主要为气相，且在管壁上，体积分数沿周向分布逐渐均匀，代表周向上的液体逐渐充满管壁，开始形成液膜。图8d)中，在管壁的液体体积分数沿周向变化不明显，可视为均匀分布，管中心处液体体积分数相对很小，代表在周向上液相充满了管壁，形成了稳定的液膜，而在管中心处则为气相区。图8e)和图8f)中，管内各处液体体积分数基本保持不变，表明液膜流动状况稳定。

本研究的模拟表现了流体在加热管内沿主流流动方向在管壁上逐渐发展形成液膜，并且受热汽化生成气体充满管中心这一过程。首先液体在导液管的作用下喷入管内打到壁面上，沿液体流动方向在周向上逐渐布满管壁，在某一高度(Z=550～900 mm)完全布满，开始形成稳定的液膜，其后液体继续沿流动方向流动，液膜状况保持不变，在整个过程中液体加热生成气体充满管中心。

2.3　液膜厚度的理论计算与模拟验证

采用数值分析法[15]和设计计算法[16]计算该降液过程液膜厚度，分别得出液膜厚度为0.727 mm和0.760 mm。2种方法第1种算法为数值计算，进行了部分简化，第2种算法以工程半经验式为基础。

2.3.1理论计算出的液膜厚度

设计计算法所用数据如表2所示。

表2　设计计算法所用参数Table 2　Parameters of design calculation method

表2中，Reff=2460Pr-0.616=2460×937.3-0.616=36.33，Reff是完全层流的临界雷诺数，Re是液膜的雷诺数，由于Re

(11)

2.3.2CFD模拟得出的液膜厚度

CFD模拟计算中，液膜在管内充分发展后，将发展段内2个不同横截面(Z=1.0和=1.8 m)上液体体积分数沿径向的分布作图9和图10。

图9　Z=1.0 m和Z=1.8 m横截面径向液体体积分数分布Fig.9　Liquid volume fraction at radial line in cross section (Z=1.0 m and Z=1.8 m)

由图9分析可知：1) CFD模拟得出的液膜厚度为1.0 mm，图9中两侧壁面附件液体体积分率接近1，经过一定厚度，液体体积分率剧烈下降，趋于平缓，这一厚度即为所形成的液膜厚度，与数值分析法及设计计算法得到的液膜厚度0.727 mm和0.760 mm相差不大，相对误差分别为27.3%及24%。CFD模拟采用三维计算方法，所取截面为Z=1.0 m和Z=1.8 m，来考察液膜厚度，得出厚度为1.0 mm，沿主流方向随着蒸发程度的不同，液膜厚度是有变化的，而两种理论计算方法作了一些简化和假设，所得到的液膜厚度属于平均值，不随主流方向变化，因此数值模拟结果与理论计算存在一定的误差。2) 在实际过程中，由于液体在流动过程中会受热蒸发，距离进口较近的截面上液膜厚度要大于距离较远的截面上的液膜厚度，Z=1.0 m处的液膜厚度应略大于Z=1.8 m的液膜厚度，但液膜厚度的减小并不明显，和主流平均速度减小的原因基本一致。3)Z=1.0 m处的截面中心的液体体积分数比Z=1.8 m的截面中心的液体分数大，在液体降膜蒸发的过程中，沿液体主流方向，管中心区域会被更多汽化的气体占据。

3　结论

采用计算流体力学的方法对降膜蒸发器进行了三维数值模拟研究，考察了降膜蒸发器加热管内的两相流动流场分布，研究了加热管内液体沿壁面膜状下降并受热汽化的现象，可以得出以下结论：

1) 降膜蒸发器的受液槽上设计一个向下的圆管进口，圆管上配有螺旋槽能够保证液体进入导流管前在受液槽上有一定的持液量，研究了单管的流场和两相流分布情况。

2) 通过计算流体力学模拟可以准确表现流体在加热管内沿流动方向在管壁上逐渐发展形成液膜，并且受热汽化生成气体充满管中心的过程。沿液体流动方向在周向上逐渐布满管壁，在某一高度开始形成稳定的液膜，其后液体继续流动，液膜厚度保持不变。

3) 计算流体力学模拟得出的液膜厚度约为1.0 mm，和由数值分析法及设计计算法的理论计算得出的液膜厚度大小基本一致，误差在可接受的范围内，并且CFD模拟采用三维模拟，更准确的反应了加热管径向的气液两相变化。

参考文献：

[1]Brotherton F. Alcohol recovery in falling film evaporators[J]. Applied Thermal Engineering, 2002, 22: 855-860

[2]Prost J S, Gonzaleza M T, Urbicain M J. Determination and correlation of heat transfer coefficients in a falling film evaporator[J]. Journal of Food Engineering, 2006, 73: 320-326

[3]Bu X, Ma W, Huang Y. Numerical study of heat and mass transfer of ammonia-water in falling film evaporator[J]. Heat and Mass Transfer, 2012, 48:725-734

[4]王永福，周荣琪，段占庭. 垂直管降膜蒸发传热研究进展[J]. 山东化工, 2001, 30: 27-30

Wang Yongfu, Zhou Rongqi, Duan Zhanting. Progress of study on heat transfer of falling film evaporation on the vertical tube[J]. Shandong Chemical Industry, 2001, 30: 27-30(in Chinese)

[5]Fujita I, Hihara E. Heat and mass transfer coefficient of falling-film absorption process[J]. International Journal of Heat and Mass Transfer, 2005, 48: 2 779-2 786

[6]Medrano M, Bourouis M, Coronas A. Absorption of water vapor in the falling film of water-lithium bromide inside a vertical tube at air-cooling thermal conditions[J]. International Journal of Thermal Sciences, 2002, 41: 891-898

[7]Jabrallah S B, Belghith A, Corriou J P. Convective heat and mass transfer with evaporation of a falling film in a cavity[J]. International Journal of Thermal Sciences, 2006, 45: 16-28

[8]Medrano M, Bourouis M, Perez-Blanco H,etal. A simple model for falling film absorption on vertical tubes in the presence of non-absorbable[J]. International Journal of Refrigeration, 2003, 26: 108-116

[9]Nusselt W. Die oberflachenkodensation des wasserdamqfes[J]. Zeitschrift Verein Deutscheringenieure, 1916, 60: 541-546, 569-575

[10]Chun M H, Park S J. Effects of turbulence model and interfacial shear on heat transfer in turbulent falling liquid films[J]. International Communication in Heat and Mass Transfer, 1995, 22(1): 1-22

[11]Feddaoui M, Belahmidi E M, Mir A,etal. Numerical study of the evaporative cooling of liquid film in laminar mixed convection tube flows[J]. International Journal of Thermal Sciences, 2001, 40: 1 011-1 020

[12]Fluent Inc. FLUENT User’s Guide[M]. 2005

[13]Lee W H. A pressure iteration scheme for two-phase flow modeling[J]. Multiphase Transport: Fundamentals, Reactor Safety, Applications, 1980: 407-432

[14]Brain S T. Thermal-Fluid analysis of lithium vaporizer for a high power magneto plasma dynamic thruster[D]. Worcester Polytechnic Institute, 2006

[15]柴诚敬，张国亮. 化工流体流动与传热[M]. 北京：化学工业出版社，2004

Chai Chengjing, Zhang Guolang. Fluid flow and heat transfer[M]. Beijing: Chemical Industry Press, 2004

[16]刘先开. 最新热交换器技术手册[M]. 北京：中国知识出版社，2009

Liu Xianka. The new heat exchanger technology handbook[M]. Beijing: China Knowledge Press, 2004