降膜蒸发冷凝器是一种高效的蒸发换热设备,具有温差推动力小、工作寿命长、结构紧凑及压降小等优点,被广泛应用于石油、化工和海水淡化等领域。按液体分布壁面的形式,降膜蒸发冷凝器主要分为板式、竖直管式及水平管降膜蒸发器等。近年来,随着能源的短缺和环境的恶化,竖管降膜蒸发器因其换热效率高、占地面积小备受人们的关注。

针对气-液界面上复杂的相互作用导致的液膜无规则变形、波动等复杂运动现象,国外学者对降膜蒸发过程中波的特性进行了大量研究。Sekine与Kumada等[1]通过实验采用容积法测量了液膜的波形和速度特性。Mudawar 等[2]通过实验评估大波在降膜液膜质量和动量传递中的作用,结果表明在降膜的薄膜层中,大部分流量通过大波传输。Kil 等[3]使用积分方法分析了在垂直圆管上具有波浪的降膜波长,研究显示波的特性与参数有关。Miraya[4]最初提出用数值模拟的方法研究倾斜板上降液膜交界处波的流动和传热特性,提出孤立波由大振幅滚波和小振幅毛细波组成。Shmerler等[5]通过实验研究了自由下落湍流液膜的显式加热,给出了平均全展开传热系数与雷诺数和普朗特数的关系。国内学者就降膜管型做了大量的研究,朱晓静、陈自刚等[6,7]对水平圆管的瞬态降膜过程与影响降膜的变量进行了数值模拟和试验验证。赵志祥、谭起滨、罗林聪和张猛等[8,9,10,11,12]对异型管进行降膜分析,指出异型管降膜在降膜分布和传热传质效率上存在很大优势。刘升等[13]模拟了较高流速下强迫竖直圆管的降膜,分析了液膜前段的托举现象。

实际的气-液两相逆流传热传质过程中的液膜形态变化、速度变化、温度变化和流场中的压力分布相当复杂,前人对降膜冷却的研究主要集中在横管实验研究和模拟上。文中采用FLUENT软件对竖管降膜的结构变量进行数值模拟,为竖直管式降膜蒸发冷凝器结构和性能的优化提供理论支持。

以锥形液体分布器做液膜分布器,模拟不同液膜分布器结构参数对降膜膜厚、速度与对流换热强度的影响。竖管降膜结构见图1。以图1中锥形分布器下部竖直圆管为计算模型,圆管直径45 mm。计算模型作如下假设,①圆管壁厚对两相传热、传质影响很小,故忽略壁厚。②圆管为轴对称图形。

降膜蒸发冷凝器工作时,锥形分布器的插件与圆管形成一定的环形缝隙,降膜冷却液从间隙处流过覆盖在圆管壁面上, 使得液膜均匀分布。 锥形分布器中心为空的圆柱体,也是空气的出口。锥形分布器下端为凹型,可防止液体在下端集聚。

图1所示的竖管降膜结构为旋转轴对称图形,建立二维平面模型并采用四边形网格对其进行划分,所得模型网格划分图见图2。

为清晰捕捉气-液与液-固交界面的流动,对管壁附近作加密处理,比较了不同网格数下的液膜厚度值,见图3。

经网格无关性验证,网格数为54 649、103 681和192 743时的液膜厚度值相差不大,综合考虑模拟的精度和计算机的运行速度,最终采用网格数103 681进行模拟计算。

连续性方程(质量守恒):

∂ρ∂t+∇(ρV)=0$\frac{\partial \rho }{\partial t}+\nabla (\rho V)=0$ (1)

动量守恒方程:

∂∂t(ρv⃗ )+∇(ρv⃗ υ⃗ )=−∇p+μ(∇v⃗ +∇v⃗ T)+ρg⃗ +F⃗          (2)$\begin{array}{l}\frac{\partial }{\partial t}(\rho \overrightarrow{v})+\nabla (\rho \overrightarrow{v}\overrightarrow{\upsilon })=\\ -\nabla p+\mu (\nabla \overrightarrow{v}+\nabla \overrightarrow{v}{}^{\text{Τ}})+\rho \overrightarrow{g}+\overrightarrow{F}(2)\end{array}$

能量守恒方程:

∂∂t(ρE)+∇[v⃗ (ρE+p)]=∇(keff∇T)+Sh$\frac{\partial }{\partial t}(\rho E)+\nabla [\overrightarrow{v}(\rho E+p)]=\nabla (k{}_{\text{e}\text{f}\text{f}}\nabla {\rm T})+S{}_{\text{h}}$ (3)

对于大流量、小缝隙的环流降膜,采用湍流计算更贴合实际情况,本模拟中根据流质的湍动特性采用RNG k-ε模型。

湍动能k方程:

∂∂t(ρk)+∂∂xi(ρkui)=∂∂xj(αkμeff∂ε∂xj)+Gk+Gb−ρε−YM+Sk         (4)$\begin{array}{l}\frac{\partial }{\partial t}(\rho k)+\frac{\partial }{\partial x{}_i}(\rho ku{}_i)=\frac{\partial }{\partial x{}_j}\left(\alpha {}_{\text{k}}\mu {}_{\text{e}\text{f}\text{f}}\frac{\partial \epsilon }{\partial x{}_j}\right)+\\ G{}_{\text{k}}+G{}_{\text{b}}-\rho \epsilon -Y{}_{\text{Μ}}+S{}_{\text{k}}(4)\end{array}$

扩散ε方程:

∂∂t(ρε)+∂∂xi(ρεui)=∂∂xj[(μ+μtσε)∂ε∂xj]+C1ε(Gk+C3εGb)−C2ερε2k+Sε         (5)$\begin{array}{l}\frac{\partial }{\partial t}(\rho \epsilon )+\frac{\partial }{\partial x{}_i}(\rho \epsilon u{}_i)=\frac{\partial }{\partial x{}_j}\left[\left(\mu +\frac{\mu {}_{\text{t}}}{\sigma {}_{\text{ε}}}\right)\frac{\partial \epsilon }{\partial x{}_j}\right]+\\ C{}_{1\text{ε}}(G{}_{\text{k}}+C{}_{3\text{ε}}G{}_{\text{b}})-C{}_{2\text{ε}}\rho \frac{\epsilon {}^2}{k}+S{}_{\text{ε}}(5)\end{array}$

式(1)～(5)中,ρ为流体密度,kg/m3;V为流体体积,m3;p为压力,Pa;μ为动力黏度,μeff为有效黏度,kg/(m·s);g⃗ $\overrightarrow{g}$为重力加速度,m/s2;F⃗ $\overrightarrow{F}$为体积力源项,kg/(m2·s2);v为流体速度沿重力方向的速度分量,u为垂直重力方向的速度分量,m/s;E为能量,J;T为温度,K;keff为有效热传导系数;μt为时间t下的流体黏性系数;Sh为流体源项;k为紊流脉动动能;Gk为由层流速度梯度而产生的湍流动能;Gb为由浮力产生的湍流动能;ε为紊流脉动动能的耗散率,%;YM为可压缩湍流中过渡扩散产生的波动;Sk与Sε为用户自定义参量,C1ε、C2ε和C3ε为常量,αk与σε为k方程与ε方程的湍流普朗特数。

液相入口为速度入口,气相入口为压力入口,液相和气相出口均设为压力出口,参考压力均为大气压。液体与壁面的接触角设为0°,液相入口的体积分数为1,管内充满空气。气相与液相的进口温度均设置为300 K,壁面作为交换热源的源头,温度设置为350 K。圆管中心轴线作为对称边界条件。

多相流模型选用数值模拟方法中应用较为广泛的VOF算法。定义空气为主相,冷却水为第二相。湍流模型选用RNG模型,采用有限体积法(FVM)控制方程的离散,压力-速度耦合算法适用于非稳态计算的PISO算法,压力离散方式选用Presto算法,动量方程与气-液面追踪方法均选用精确度很高的二阶离散格式。降膜流动的流体物性参数见表1。

在竖管内冷凝水与管壁的蒸发冷凝换热过程中,布膜缝隙处的宽度对液膜膜态、膜厚分布及传热和传质效果的影响不明朗。以文献[14]中的降膜宽度为依据,采用控制变量法,保持其它变量不变,以不同的环缝隙宽度(液体分布器下端与管壁的距离)W为模拟研究变量。喷淋量为0.74 kg/(m·s),环缝隙宽度W为1.2 mm、1.5 mm、2.0 mm时,竖管管壁降膜的液膜膜厚δ与x方向分速度vx沿管长L的分布情况分别见图4和图5。

图4中δ曲线的波峰对应图5中vx曲线的波谷,波谷对应着波峰。图4中δ曲线的波峰越高,图5中相应的vx越小。图4中上段液膜平稳,附着良好,下段为不规则的液膜凸起,凸起由孤立波、小张力波和层流底层组成[15]。由此也印证了文献[16]中的膜厚分布并非简单的简谐振动的结论。随着环缝隙宽度的增大,液膜在重力、表面张力和黏附力的作用下,膜态分布先稳定后趋于波动,基层膜厚逐渐减薄,液膜凸起增高。图5c中vx在W=2.0 mm时出现方向向左的最大值,是因为降膜量增加,加速度增大,阻力也随之增大。

数值模拟得到的努塞尔数Nu与环缝隙宽度关系曲线见图6。从图6可知努塞尔数随环缝隙宽度的增大而增加,液膜对壁面的换热能力加强。

降膜宽度的增加,液膜凸起个数即回流区个数增加,加之基层膜厚的减薄是液膜对流换热强化的主要原因。但对流换热强度在宽度W=1.5 mm相比W=1.2 mm时增加了4.3%,在宽度W=2.0 mm比W=1.5 mm时增加4.8%,可见增大降膜宽度一定程度上可增大液膜的冷却能力,但效果不明显。

逆向流动的空气作用是把液膜表面热对流传递得到的热量与蒸发的水蒸气分子带走,从而起到冷却作用。文中模拟了在喷淋量为0.74 kg/(m·s)和环缝隙宽度为1.2 mm时不同空气进口宽度B下的液膜分布情况,见图7。

图7表明,随着B的增大,液膜波动程度整体逐渐增大,液膜凸起增高,凸起的形态逐渐呈现出不规则性,出现液滴飞溅的现象。B=6 mm时液膜的形态为低振幅对称波,B为12 mm与18 mm时液膜形态转化为高振幅侧向波,这与文献[17]中实验研究的膜形态结构相一致。

从图7a～图7c和图4c中可看到,随着B的增加,整体上的液膜波动程度或范围逐渐增大。原因是B增大,空气通量增加,蒸汽压差逐渐增大, 液膜蒸发加快,下段管壁的液膜持液量减少,再加之空气阻力的增大,使得液膜波动增大并向上端移动。在B≥12 mm时,基层膜厚出现最小值,但无干斑,说明液膜的附着性良好。

数值模拟得到努塞尔数Nu与空气进口宽度关系曲线见图8。

从图8可知,液膜与壁面的对流传热强度随空气进口宽度的增大而增大,空气进口宽度增大,即空气通量增大,蒸发的水汽所携带的热量能很快被带走,液膜的表面更新速率加快,尤其在降膜的下端,有较薄的基层液膜、大的波动和波动范围,且波动凸起中存在着回流区,故而管壁对液膜的传热热阻较小,传热得到加强。B=18 mm时的Nu比无空气通过时的Nu增加14.8%。

管长对液膜的均匀性分布起着至关重要的作用,模拟在环缝隙宽度W=1.2 mm和空气进口宽度B=18 mm条件下液膜厚度δ与液膜在x方向分速度vx沿管长的分布,分别见图9和图10。

由图9和图10可知,L=0.3 m时液膜的润湿性良好,δmin=0.47 mm;0.4 m≤L≤0.6 m时液膜的分布情况基本一致,均是在L=0.3 m后出现大的波动。由膜厚的波动范围分别为0.18～2.6 mm、0.19～2.8 mm和0.10～5.1 mm可知,这种波动程度随着管长逐渐增大。在L=0.7 m时,液膜凸起的驼峰状波最高、最多,波动程度最剧烈,基层液膜厚度最薄,在降膜距离的L=0.52 m、0.6 m、0.64 m和0.69 m处,此时的vx达到最大值,δ几乎为0,管壁形成干斑,无法完整地润湿壁面,壁面干烧严重,设备的使用寿命降低。

分析其原因为,①管长增加,液膜的覆盖面积增大,冷却液在惯性力的作用下难以在大面积的管壁上附着。②流体在管壁下端的流动阻力逐渐加大,并产生较小的vx值,使得液膜的局部堆积严重,液膜不能及时输送,液膜出现断裂。③降膜过程中随着蒸发不断进行,管长越长,管壁下端的持液量越少。④液膜流动的距离越长,温度就越高,液膜的毛细吸附能力减弱[15],液膜极易脱离壁面。故而降膜管越长,液膜的稳定性越差。由此可知在L≤0.6 m时,液膜的附着性是良好的。

数值模拟得到的Nu随管长的变化关系曲线见图11。

从图11可知, Nu随着降膜管长的增加而增大,即对流换热强度在逐渐增大,一方面管长越长接触越充分,液膜吸热越多,换热越充分。另一方面降膜下段波动程度逐渐增大,液膜更新速率加快,从而液膜对管壁的吸热增强。L=0.7 m时的换热强度比L=0.3 m时增长了1.72倍,相比降膜宽度与空气进口宽度,管长因素在降膜换热中起着更重要的影响作用。

通过数值模拟的控制变量法对降膜宽度、空气进口宽度和管长进行了管壁上液膜的膜态分布和传热特性研究。数值模拟研究得出的结果表明,在一定的喷淋量下,环缝隙宽度越大,液膜的波动程度先稳定后波动,Nu逐渐增大。在环缝隙宽度W=1.2 mm条件下,随着空气入口宽度的增大,液膜的波动性增大,B=18 mm时液膜的传热特性最好。在环缝隙宽度W=1.2 mm、空气进口宽度(模拟区域)B=18 mm时,随着管长的增大,液膜的稳定性逐渐降低,但液膜对管壁的对流换热强度逐渐增大。从Nu可知,在设计降膜蒸发冷凝器时,应优先考虑管长,空气进口宽度次之,降膜宽度最后,管长能够使管壁对液膜的传热传质机理充分实现,实现冷却效率的最大化。