水平圆管外降膜流动与传热过程研究

Study of Falling Film Flow and Heat Transfer Process outside Horizontal Circular Tubes

摘要: 对水平圆管管外降膜流动过程进行数值研究，分析了换热管径为19 mm时不同喷淋密度和壁面温差对管外液膜厚度、液膜表面速度、温度分布以及传热系数的影响。研究结果表明，液膜厚度及管外局部传热系数受轴向位置、周向角度、壁面温差和喷淋密度等因素影响较大。液膜干涸区域易出现在周向角度为130˚150˚附近；液膜厚度随着受热而减小，壁面温差的增加会加剧Marangoni效应，使局部液膜厚度有所增加。管外局部传热系数沿周向角度的分布呈递减趋势，在周轴向角度为130˚150˚时最低，而在管底部时略有回升；平均传热系数随喷淋密度和壁面温差的增加而增大;当喷淋密度为0.022 kg(ms)⁻¹时，壁面温差从10℃增加到20℃，平均换热系数增加了37%。

Abstract: A numerical study of the falling film flow process outside the horizontal circular tube was carried out to analyze the effects of different spray densities and wall temperature differences on the thickness of the liquid film outside the tube, the surface velocity of the liquid film, the temperature distribution, and the heat transfer coefficient for a heat exchanger tube diameter of 19 mm. The results show that the liquid film thickness and local heat transfer coefficient outside the tube are greatly affected by axial position, circumferential angle, wall temperature difference and spray density. The liquid film drying region is easy to appear near the circumferential angle of 130° 150°; the liquid film thickness decreases with the heat, and the increase of the wall temperature difference will exacerbate the Marangoni effect, which will increase the local liquid film thickness. The local heat transfer coefficient outside the pipe shows a decreasing trend along the distribution of circumferential angle, which is the lowest when the circumferential angle is 130°150°, and slightly recovered at the bottom of the pipe; the average heat transfer coefficient increases with the increase of the spray density and the wall temperature difference; when the spray density is 0.022 kg(ms)^-1, the wall temperature difference increased from 10°C to 20°C, and the average heat transfer coefficient increased by 37%.

1. 引言

液体降膜蒸发技术因其传热系数高，在海水淡化、食品加工、电子元器件冷却、温湿度控制等领域被广泛应用。在降膜蒸发式冷却器中，换热管型、制冷工质的流动形态、液膜厚度和液膜波动程度等对热质交换特性有着重要影响，引起了许多学者对降膜蒸发机理的重视与研究 [1] [2] [3] 。

液膜厚度影响着液膜的热阻和气液界面的蒸发过程。Nusselt [4] 假设液膜为层流流动时表面是无波动的，并给出平均液膜厚度和雷诺数的准则关系式。Gstoehl等 [5] 通过实验得出了低雷诺数流动状态下管外液膜厚度随周向角度的变化规律。何茂刚等 [6] 采用激光诱导荧光法测量了常温常压下水平光滑管和强化管外液膜厚度，得出液膜在上半周较厚并且变化较大的结论。陈自刚 [7] 在不考虑传热情况下研究了换热管间为柱状流时管外液膜分布情况及液膜厚度变化规律，并得出液膜厚度沿轴向近似呈现波峰–平稳–波峰的周期性规律分布的结论。Li等 [8] 针对三种强化管束在10到110的雷诺数范围内进行对流蒸发换热实验，发现随着管外液膜雷诺数的增加，换热系数达到一个峰值后逐渐降低；但是随着管内热水温度增加，这个趋势会发生改变。牟兴森等 [9] 研究了不同布管方式所对应的三种管间距对管外降膜传热系数的影响，发现转角正方形布置下的传热系数最高，而转角三角形布置的传热系数最低，不同管间距传热系数产生差异的主要区域为管周向0˚~135˚区域。全晓宇等 [10] 研究了Marangoni热效应对降膜流动的影响，得出较大的加热温差和较低液相流速会使液膜收缩更明显，受热液膜边缘处湍动更强，而液膜中部所受影响较小。郑毅等 [11] 利用红外示踪技术，研究了低喷淋密度下超亲水表面滴状流时液膜温度的分布特点，发现在撞击区内温度变化剧烈，非撞击区内温度变化较小，且在撞击区形成一个马鞍形高温环。孙玉祥等 [12] 通过试验测量不同喷淋水温和喷淋密度对液膜表面蒸发率的影响，发现在入口水温20℃~60℃范围内，随着入口温度提高，液膜表面蒸发率显著提高，但在不同的周向角度上没有太大变化，而不同的喷淋密度对于不同周向角度上的液膜表面蒸发率有较大影响。

目前的研究中，管外对流传热系数主要由试验获得。而在管外降膜厚度分布规律的数值分析中，往往忽略传热和相变。但在实际应用中，降膜过程中的传热和相变过程不可忽视，且与液膜厚度互相影响。本文将建立管外降膜流动的三维模型进行数值模拟，分析不同加热条件和喷淋密度对水平管外降膜流动过程中液膜厚度沿管周向和轴向分布特点，以及对传热系数的影响规律。

2. 数值模拟

2.1. 物理模型

本文研究的水平管液体降膜流动换热模型如图1所示。定义自换热管顶部为起点，沿顺时针转过的角度为周向角度b (˚)，液膜厚度记为δ (mm)，管外液体喷淋密度记为G (kg·m⁻¹·s⁻¹)，为管径方向单位长度上流体的质量流量，降膜流动雷诺数为Re = 4 G/μ，其中μ为液体的动力黏度。本文所计算的雷诺数范围是300~1000，所涉及的液膜流动属于层流的范围 [13] 。

Figure 1. Horizontal tube liquid falling film flow heat transfer diagram

图1. 水平管液体降膜流动换热示意图

2.2. 数学模型

2.2.1. 体积分数方程

本文采用VOF (Volume of Fluid)方法来计算管外气液两相的流动情况，通过跟踪计算每个控制体积内第i相的体积分数 ${\alpha }_i$ 得到相与相的界面分布，控制方程如下：

$\frac{\partial {\alpha }_i}{\partial \tau }+u\cdot \nabla {\alpha }_i=0$ (1)

${\displaystyle \underset{i=l,a}{\sum }{\alpha }_i}=1$ (2)

式中， $u$ 为速度矢量， $\text{m/s}$ ；在每个计算单元中各相总和为1， $0\le {\alpha }_i\le 1$ 。

2.2.2. 质量守恒方程

$\frac{\partial \rho }{\partial \tau }+\nabla \cdot \left(\rho u\right)=0$ (3)

式中， $\rho$ 为密度， ${\text{kg/m}}^3$ 。

2.2.3. 动量方程

$\frac{\partial \left(\rho u\right)}{\partial \tau }+\nabla \cdot \left(\rho uu\right)=-\nabla p+\nabla \cdot \left[\mu \cdot \left(\nabla u+{\left(\nabla u\right)}^{\text{T}}\right)\right]+\rho g+F$ (4)

式中， $P$ 为压力， $\text{Pa}$ ； $\mu$ 为动力黏度， $\text{Pa}\cdot \text{s}$ ； $F$ 为表面张力源项。

2.2.4. 能量方程

$\frac{\partial \left(\rho E\right)}{\partial \tau }+\nabla \cdot \left[u\left(\rho E+p\right)\right]=\nabla \cdot \left(k_{eff}\nabla T\right)+S_e$ (5)

式中， $S_e$ 为能量方程源项， ${\text{W/m}}^{\text{3}}$ ； $k_{eff}$ 为有效热传导率， $\text{W/}\left(\text{m}\cdot \text{K}\right)$ 。

2.2.5. 相变模型

相变模型采用Lee [14] 模型来计算，该模型提出了一种体积传质的处理方法。

${\stackrel{\dot{}}{m}}_e=\{\begin{array}{l}-coeff{\alpha }_l{\rho }_l\frac{T_l-T_{sat}}{T_{sat}}T\ge T_{sat}\\ -coeff{\alpha }_a{\rho }_a\frac{T_{sat}-T_a}{T_{sat}}\text{T}\le T_{sat}\end{array}$ (6)

其中

$coeff=\frac{6}{d_b}\gamma \sqrt{\frac{M}{2\pi R{T}_{sat}}}L\left(\frac{{\alpha }_a{\rho }_a}{{\rho }_l-{\rho }_a}\right)$ (7)

式中， $coeff$ 为调节系数； $T_l$ 和 $T_a$ 分别为液相和气相的温度； $T_{sat}$ 为相变饱和温度； $d_b$ 为气泡直径； $\gamma$ 为传质时间弛豫系数。

2.3. 网格划分及算法设置

本文采用精度较高的六面体网格对计算区域进行网格划分，同时在气液交界面及液相区进行了网格加密。计算域的网格划分如图2所示，计算中比较了网格数为892490、938282、1027520、1229310的计算模型，结果差异较小，在达到精度要求的前提下最终选取网格数为892490的模型计算。基本计算方法如下：(1) 空间为三维直角坐标系，流场为非稳态模拟；(2) 压力–速度的耦合采用PISO算法，压力项离散方式选用PRESTO！算法，动量和能量方程离散均选用二阶迎风格式；(3) 模拟过程中考虑重力和表面张力的作用。

3. 模拟结果与分析

3.1. 数学计算方法验证

为验证模型的准确性，对文献 [5] 和文献 [6] 中的实验工况进行模拟计算，比较液柱截面的液膜厚度分布情况，比较结果如图3所示。从图中可以看出，液膜厚度的模拟值与文献中实验结果的变化趋势基本一致，与文献 [6] 的模拟值相差较大的原因是：文献中采用的模型是二维模型，而本文中采用的是三维模型，与实验值更为接近，误差更小。模拟结果表明，本文中所建立的数值模型可用于研究管外降膜流动过程。

Figure 2. Model grid division diagram

图2. 模型网格划分示意图

Figure 3. Comparison of simulated values of liquid film thickness with those in literature

图3. 液膜厚度模拟值与文献值对比

3.2. 管外液膜铺展过程

图4分别为流动时间0.012 s、0.03 s、0.048 s、0.192 s时的液膜铺展状况。由图中可以看到，液膜在降落初期，两侧流体铺展较快。随着两侧液膜向中间汇聚，管周下半部分开始出现“M”型分布，并从上层管的两液柱截面中间降落至下层管。

在T = 0.192 s时管外表面形成连续稳定的液膜。此时在速度分布云图中可以看到，液膜表面的速度并未均匀分布，但呈现出一定的规律性。在0~30˚区域液膜处于下落冲击阶段，液膜表面速度比降落初始流速要高；在40˚~70˚区间扩展过程中受到重力和壁面粘滞阻力作用，流速逐渐降低，最终两股液膜交汇后在重力作用下沿壁面向下流动。在接近90˚位置处由于液体所受重力沿流动方向的分力增大，此时液膜沿管周方向的流速逐渐增大；随着液膜进一步润湿表面，液体受到的表面张力和粘性效应增强，流速开始降低，并在底部积聚最终降落到下层管壁。而在下层管壁液膜波动较大，不能均匀润湿壁面，形成较多干涸区。

3.3. 常温下管外液膜厚度分布

管外降膜流动过程中液膜的厚度分布并不均匀，在液膜铺展过程中两液柱截面中由于液膜交汇而形成一个波峰。图5为喷淋密度为0.02 kg×(m×s)⁻¹时液膜在液柱截面和波峰截面的气液两相分布云图。从图

Figure 4. External liquid film spread and velocity distribution diagram

图4. 管外液膜铺展情况及速度分布图

Figure 5. Cloud image of gas-liquid two-phase distribution outside the pipe

图5. 管外气液两相分布云图

中可以看出波峰处液膜厚度明显大于液柱处，而且在圆管顶部和尾部区域液膜厚度较大，中间区域液膜较薄。

为研究液膜厚度沿换热管轴向的分布规律，本文分析了自圆管顶部周向角度为30˚、60˚、90˚、120˚和150˚时气液两相分布及液膜厚度，如图6和图7所示。从图中可以看出：(1) 在不同周向角度处，液膜厚度沿相邻两个布液孔中心对称分布；(2) 在b = 30˚处，由于尚未脱离液面冲击区域的影响，液膜厚度在波峰和液柱区域都较高；上半管周在b = 60˚~90˚区间液膜厚度最小；下半管周在b = 120˚附近的区域，液膜厚度较大，在b = 150˚附近液膜厚度降低。(3) 在距离液柱截面3~7 mm处下半管周液膜厚度迅速增大，而上半周的液膜厚度在此区间略有降低。(4) 在b = 120˚区域，液膜厚度在距离液柱4 mm和16 mm时有一个“小波峰”，由此可见下半周的液膜波动比上半周更为强烈。

Figure 6. The thickness of liquid film is distributed axially at different circumferential angles

图6. 不同周向角度下液膜厚度沿轴向分布

Figure 7. Gas-liquid two-phase distribution outside the pipe at different circumferential angles

图7. 不同周向角度下管外气液两相分布图

3.4. 壁面温差和喷淋密度对液膜厚度分布的影响

3.4.1. 壁面温差对液膜厚度分布的影响

为了分析加热条件对液膜厚度的影响，选取上排管波峰截面和液柱截面进行研究，获得了同一喷淋密度下，不同壁面加热温差下液膜厚度随周向角度变化的情况。

Figure 8. The change of liquid film thickness with circumferential Angle of wave crest section

图8. 波峰截面液膜厚度随周向角度的变化

由图8可知，冷却温差对液膜厚度影响较大，且随着位置的变化而变化。在波峰截面b = 0˚~120˚区间内，温差为10℃~15℃时的液膜厚度明显低于无温差时，主要因为管外表面液膜受热后，水的粘性降低，进而导致液膜流速增加，液膜变薄。温差为15℃时的液膜厚度比温差为10℃时大，这是因为受热液膜下降过程存在较强的Marangoni效应，温差越大，液膜从温度较低处向较高处收缩的现象越发严重；而波峰截面由于液膜较厚，温度更低，更易促进“波峰”两侧液膜向中心收缩，因而波峰处因壁面温差增加而增大。在b > 120˚时，从图中可以看到，波峰截面受热状态下的液膜厚度陡然上升并超过未加热时，这是因为液膜在加热状态时，周向角度130˚~150˚的液柱截面区间由于液膜太薄，容易干涸导致管壁过热，此时两侧液膜与中心液膜温差增大，Marangoni效应更加明显，液膜沿轴向收缩流向“波峰”处。

3.4.2. 喷淋密度对液膜厚度分布的影响

为确定不同喷淋密度对液膜厚度的影响，对加热温差为15℃时，不同降膜雷诺数下管外波峰截面液膜厚度的分布情况进行研究。图9的计算结果表明：喷淋密度在0.014~0.020 kg×(m×s)⁻¹范围内时，液膜厚度沿管周变化趋势一致，然而，当喷淋密度增加至0.022 kg×(m×s)⁻¹时，液膜厚度在b > 120˚之后的趋势发生了改变，并未继续增大；而是在b = 150˚左右略微下降，随后上升，并在较小范围内波动。这是由于随着液膜喷淋量增大，圆管外壁液膜速度较大，圆管底部积聚的液膜受到表面张力的影响减小，下半周液膜厚度反而降低。

3.4.3. 管外表面润湿情况分布

管外液膜破裂会造成液膜在传热管上无法均匀铺展，因此而导致的液膜过厚会使局部热阻增加，不利于液膜内的导热和气液界面蒸发的进行；而液膜过薄则会产生局部干涸，均严重影响传热。因此降膜流动过程中管壁的润湿面积也是评价换热的重要因素。图10给出了不同壁面温差下管壁表面被液膜覆盖的百分比。从图中可以看出，当壁面温差为15℃时，液膜表面覆盖率比壁面温差为10℃和20℃时更高，这是因为在一定范围内，随着液膜温度的升高，水的粘性下降，液膜变薄，铺展面积增大；但随着温差进一步增大，表面蒸发量变大，Marangoni效应更强，液膜更易破裂，干涸区域增大，液膜润湿面积减少。但是随着喷淋密度增加，温差对润湿面积的影响逐渐减弱。

Figure 9. Effect of different spray density on liquid film thickness

图9. 不同喷淋密度对液膜厚度的影响

Figure 10. Percentage of wetting area of liquid film surface under different wall temperature difference

图10. 不同壁面温差下液膜表面润湿面积百分比

3.4.4. 不同壁面温差下液膜内温度分布情况

管外液膜沿管周的不同角度，各个截面的液膜平均温升值如图11所示。由图可知，在管外液膜下降过程中，由于不断从管壁吸收热量，温度逐渐升高，且管壁与初始降膜温差值越大，液膜下降过程中温升则越大。在b = 120˚~150˚范围内，液膜平均温度略有降低，这是由于在此范围内，距离液柱4~8 mm处液膜较薄，液膜表面温度较高。此处液膜表面张力比液柱和波峰截面处更小，为了消除液膜内温度梯度的变化，表面温度较高的液膜在表面张力梯度的作用下向两侧收缩，使液膜较厚处的厚度增大，液膜较薄处厚度减小，导致液膜更容易破裂。由于液膜覆盖处膜厚较大，传热热阻也变大，此时该处液膜的平均温度反而有所降低。

Figure 11. Distribution of mean temperature rise of liquid film at different circumferential angles with heating temperature difference

图11.不同周向角度液膜平均温升随加热温差分布规律

3.5. 降膜传热分析

3.5.1. 管外局部传热系数

图12(a)为15℃加热温差条件下局部传热系数沿圆管周向的变化趋势。由图可见，当液膜流经圆管顶部时局部传热系数最大，随着角度增大，传热系数迅速减小，在120˚~160˚区域内降到最低，之后略有增高。这是由于圆管顶部受液体下落冲击的影响较大，而后热边界层沿着周向不断发展，热阻逐渐增加，管外液膜传热系数降低。在120˚~160˚区域，距离液柱截面0~6 mm时容易形成局部干涸，传热系数最低；而在8~10 mm处由于液膜厚度较大，液膜热阻也更大，传热系数也较低。在圆管底部，两侧液体向下流动发生相互碰撞、脱落，加剧了膜内扰动使传热性能突增。

图12(b)为管外局部传热系数沿圆管轴向的位置分布规律，从图中可以看出，在圆管上半周，局部传热系数从液柱截面处开始沿轴向先是缓慢降低，在距离液柱截面7 mm后迅速下降，在圆管中部呈现出一定的波动趋势，但波动范围不大，且整体趋势下降。在圆管下半周0~7 mm区域由于出现干涸区传热系数较低，在近波峰截面区域内，传热系数逐渐上升。圆管底部由于两侧液体的碰撞、汇聚扰动较大，传热系数反而较高，而且沿着轴向分布变化较小。

3.5.2. 管外平均传热系数

图13为不同加热温差和喷淋密度对壁面平均传热系数的影响。从图中可以看出，平均传热系数随加热温差和喷淋密度的提高明显增加。以喷淋密度0.022 kg×(m×s)⁻¹为例，加热温差从10℃提高到20℃，平均传热系数提高了37%。这主要是因为当壁面温度增加时，水的平均温度也升高，导热系数相应提高，同时粘度下降，导致表面张力下降，液膜内扰动增强，促进了热量的传递。

Figure 12. Local heat transfer coefficients are distributed along the circumference and axis of the tube

图12. 局部传热系数沿圆管周向和轴向分布

Figure 13. Influence of different wall temperature difference on the average heat transfer coefficient of the wall

图13. 不同壁面温差对壁面平均传热系数的影响

4. 结论

本文以水为工质，模拟研究了水平管间的降膜流动传热过程，获得了喷淋密度、壁面温差对管外液膜厚度的分布以及对传热系数的影响规律，在本文研究范围内得出以下结论：

1) 在轴向上，液膜厚度沿相邻两个布液孔中心线呈对称分布；在周向上，液膜厚度逐渐降低，当位于周向角60˚及150˚附近时液膜厚度最小，之后液膜在底部积聚，厚度升高。

2) 壁面温差对液膜分布有较大影响，随着壁面温差增大，在沿周向130˚~150˚的液柱截面处容易形成干涸区域，波峰截面则因Marangoni效应的影响厚度增加。

3) 管外液膜厚度及流动状态造成管外局部换热系数有较大差异。管外局部传热系数沿周向分布呈递减趋势，在周向角130˚~150˚时最低，之后到达底部时略有回升。

4) 在喷淋密度为0.014~0.022 kg×(m×s)⁻¹范围内，管外平均传热系数随喷淋密度和壁面温差的增加而增大。当喷淋密度为0.022 kg×(m×s)⁻¹时，壁面温差从10℃增加到20℃，平均换热系数增加了37%。

NOTES

^*通讯作者。