大平板上液膜冲击附板行为的数值模拟研究

胡 真，胡 珀，杨小杰，杜卡帅

(1.上海交通大学 核科学与工程学院，上海 200240；2.上海核工程研究设计院有限公司，上海 200233)

在第3代自主先进压水堆CAP1400的设计中，非能动余热排出(PRHR)系统在事故状态下的长期冷却过程中将发挥至关重要的作用。该系统中的安全壳内置换料水箱(IRWST)在此过程中将蒸发产生的水蒸气通过冷凝回流系统形成冷凝液并最终返回IRWST，可保证IRWST的水量可利用率，从而支持一回路的长期冷却。在此过程中，安全壳垂直壁面上的附板和返回槽会造成一定量的冷凝液损失，而该损失将直接影响冷凝回流的份额，因此预测冷凝回流的损失量对于PRHR长期冷却分析具有重要意义。

降液膜在安全壳外部壁面上的研究目前已有显著成效[1-7]，但关于降液膜在安全壳内壁面上冲击附板时的流动行为却鲜有研究，而对于流体冲击板面行为的研究也多为研究单个液滴的溅射和铺展行为[8-9]，因此当前国内外尚无可行的计算方法来有效估算该损失。商用计算软件FLUENT中的EWF模型(欧拉壁面模型)具有计算效率高、能很好模拟液膜在壁面的流动等优势[4,10]，在模拟全尺寸的模型时计算成本相对较小，故本文将以EWF模型为工具，结合相关实验数据，通过构建合理的数值模型，研究大平板上液膜冲击附板时的溅射行为，旨在通过研究液膜冲击附板的过程为冷凝回流损失率做出相关预测。

1 数学模型

1.1 模型建立及边界条件

图1 WAFT实验台架[2] Fig.1 WAFT experiment bench[2]

液膜实验装置(WAFT)如图1所示。实验环境为常温常压，实验所用平板涂层为Carbozinc 11 HS有机富锌涂料。实验过程中液膜分配箱将水均匀分配到平板上端，水以漫溢的方式沿平板平铺而下，冲击附板，发生溅射后的部分水流经返回槽收集后进行称重以评估实际的损失。

实验所用实验板及其焊接方式如图2所示，板长5 m、宽1.2 m，板端末设有高0.1 m的返回槽，距返回槽2 m处置有附板。

图2 实验板设置及附板焊接示意图Fig.2 Plate setting and welding diagram of attached plate

根据实验要求选用不同厚度和宽度的附板，附板通过角焊的方式固定在板面上，焊脚高度为l，附板高度为H，焊脚高比β定义为：

β=l/H

(1)

在实验过程中，通过控制变量法改变液膜分配箱流出水的流量、温度、附板尺寸以及焊脚高比来分析各因素对液膜溅射的影响。

在数值研究过程中，选择实验板及其上方0.1 m高处气-液流动的三维空间区域作为计算域，模拟计算中各尺寸采用与实验1∶1的比例进行建模计算，模型切面示于图3。

图3 三维模型切面示意图Fig.3 Section diagram of 3-D model

具体模型尺寸列于表1。由于实验板上方未加玻璃罩遮挡，故将上方区域边界设为压力出口。在计算区域内，液膜由上方入口处流入后沿壁面铺展开，待其发展充分后冲击位于板面下部的附板，飞溅损失的液体由上方压力出口射出，其余液体通过高0.1 m的返回槽流出。具体边界条件设置列于表2。

表1 三维模型的几何尺寸Table 1 Geometry dimension of 3-D model

表2 边界条件Table 2 Boundary condition

为方便计算，本文在建模过程中做了如下假设：1) 计算过程中不考虑液膜的波动特性；2) 不考虑液膜与空气之间的质量能量交换，即不考虑液膜在冲击附板和流动过程中的蒸发损失；3) 不考虑实验中返回槽上的滤网损失以及钢板失配影响。

1.2 控制方程

由于本文只研究大平板上降液膜的流动和溅射行为，因此传热情况和能量方程不在考虑之中。假设流体不可压缩，其稳态下连续性方程为：

∂(ρui)/∂xi=Smax

(2)

式中：ρ为流体密度；ui为流体在xi方向的速度；Smax为源项。

动量方程可描述为：

(3)

式中：i,j为卡笛尔坐标系中的不同坐标方向；u为速度，m/s；p为压力，Pa；μ为动力黏度，Pa·s；δ为Kronecker函数；Smon-i为源项；ρu′iu′j为基于Boussinesq假设的雷诺应力。

本文采用Realizablek-ε湍流模型描述液膜的湍流流动，该湍流模型的模拟输运方程k和ε方程分别如式(4)、(5)所示。

Gk+Gb-ρε+Sk

(4)

式中：k为湍流动能；Gk为由主流速度梯度产生的湍流动能；Gb为由浮升力产生的湍流动能；σk为k方程的普朗特数，σk=1.0；ε为湍流耗散率；Sk为用户自定义的源项。

(5)

EWF模型可用来预测壁面上薄液膜的产生和流动情况，此外，该模型还可用来模拟液滴飞溅、剥离、液膜分流等现象。该模型适用于液膜厚度较薄、膜内液体流动与壁面平行的情况[12]。

液膜的质量守恒方程为：

(6)

液膜的连续性方程为：

(7)

式(7)右边第1项表示动压气体压力在壁面方向上的重力分量以及液体表面张力对液膜的影响，第2～4项分别表示在液膜流动方向上的重力作用、气膜交界面处的黏性剪切力作用和膜内的黏性作用。

本文所用的液膜分离模型[13-14]为O’Rouke模型。FLUENT中对液膜分离的判定准则为：

(8)

在底板和附板的边界设置中均需开启EWF模型，其中附板和焊缝的DPM模型选项设置为reflect，而底板的该选项设置为trap，确保液膜在冲击附板前贴板流动。

本文模型选用基于压力的稳态求解器，由于底面设有附板，为更好地实现计算收敛，选用标准近壁面模型进行计算。

1.3 网格无关性分析

考虑到计算的准确性，采用结构化网格划分的方式对计算区域进行网格划分。由于应用了标准近壁面模型，贴近实验板的第1层网格的y+值保持在20左右，其他壁面附近的第1层网格的y+值保持在20以上。本文选取24万到96万不同数目的网格模型，对每10 min的损失水量进行检测，网格无关性验证结果如图4所示。为综合计算效率和计算准确性，本文最终选用48万网格的模型。

图4 网格无关性验证Fig.4 Grid independence validation

2 结果与分析

2.1 冷凝回流损失比较

图5 冷凝回流损失结果比较Fig.5 Comparison of condensation reflux loss results

为验证数值模型在预测冷凝回流损失结果上的正确性，本文选取焊脚高比为1.0时不同工况下的模拟计算值与实验值进行比较。图5为不同流量的降液膜冲击不同厚度和长度附板的计算值与实验值对比。从图5可明显看出，小流量下液膜冲击附板后损失率很小，且计算值与实验值极为接近；中流量(500 m3/h)和大流量(1 000 m3/h)情况下，液膜损失率显著增大，且计算值和实验值间的误差也有明显提高，个别工况的计算值与实验值间的误差超过5%。总体而言，通过本文所建数值模型得到的冷凝回流损失量与实验误差能控制在8%以内，即本文所建数值模型能较好地预测冷凝回流损失量。

2.2 降液膜冲击附板时的溅射行为分析

降液膜冲击附板的整个行为过程如图6所示，降液膜在平板上经过充分发展，冲击附板后的行为过程可分为4个阶段，即液膜冲击附板阶段、第1次液膜再聚合阶段、液膜分离阶段、第2次液膜再聚合阶段。

液膜厚度和液膜流动的速度场如图7所示，图中前3阶段液膜行为的速度场速度范围均为0～1 m/s，第2次液膜再聚合阶段的速度场速度范围为0～2 m/s，且速度方向均为沿板面方向。

在液膜冲击附板阶段，液膜冲击附板时，其流动方向发生改变，并沿焊缝方向继续流动。在液膜冲击附板后第1次发生分离时，液膜主要分成3个方向流动：根据无滑移条件假设，黏性底层因与壁面保持相对静止状态，因此仍紧贴壁面流动，湍流核心区的液体受表面张力和液体内部作用力影响，大部分液体在沿焊缝方向流动一定距离后继续回到附板表面贴壁面流动，少部分液体克服表面张力影响脱离液膜飞溅射出。

图6 液膜冲击附板的行为过程分析Fig.6 Analysis of film impinging on attached plate

图7 液膜厚度和速度场分析Fig.7 Analysis of falling film thickness and velocity field

在液膜第1次聚合阶段，当附板宽度足够短、焊缝高度足够高时，液膜在冲击附板后可能无法在附板上完成第1次再聚合而直接以射流形式流入返回槽中，这从附板上液膜厚度小于附板前的液膜厚度分布可看出。

在降液膜冲击附板的过程中，液膜冲击附板后部分克服表面张力作用飞溅射出的液体是造成液膜损失，即冷凝回流损失的主要原因。

(9)

在液膜的第2次聚合阶段，由于液膜在脱离附板的过程中速度存在差异，因此脱离壁面射出的液体回到主板进行聚合所需时间不同，导致附板后的平板上液膜厚度分布呈现不均匀的现象。而从附板后平板上液膜厚度大于附板上液膜厚度可发现，参与第2次聚合的液滴中有部分液滴在冲击附板后直接越过附板飞回到平板上。

2.3 液膜流量对液膜损失的影响

图8 液膜流量对液膜损失的影响Fig.8 Influence of film flow rate on film loss

图8为液膜流量对液膜损失的影响。可看出，在模拟范围内，小流量下液膜损失率较小，随着流量增大，2、3号附板的损失率先增大后略有减小，而1号附板的损失率在大流量下出现“突增”现象。附板尺寸列于表3。

表3 附板尺寸Table 3 Size of different attached plates

通过增加不同流量的数值模拟点，发现随着流量的逐渐增大，1号附板出现损失率变化的突增点。分析其原因为：由于1号附板的宽度最小、厚度最大，从而液膜在到达附板时具备的速度最大，且沿附板流动时获得的垂直于附板方向的速度最大。

根据O’Rourke等[14]和Wu等[15]的研究，液膜的分离过程和液滴从液膜中剥离主要受We影响，而上述过程是影响液膜在冲击附板后损失的主要原因。由此，定义附板与返回槽间的高度差和返回槽高度的比值为液膜溅射空间数He：

He=(L-H)/L

(10)

式中：L为返回槽高度(100 mm)；H为附板高度。

从而通过比较液膜损失率ηloss和We以及He之间的关系，得到关系式(11)：

ηloss∝ΔWex·Hey0

(11)

式中：ΔWe=|We-Wecr|，Wecr为液膜冲击附板时导致液膜损失率出现“突增”现象时的We，在本文模拟范围内取1.5；x和y均为常数。由于小流量下液膜损失率较小，且在模拟过程中未表现出明显规律，预测价值较小，故本文主要针对较大流量(Q=500 kg/h和1 000 kg/h)分析液膜损失率和以上二者的关系。

在本文模拟范围内，分别取y=-1.2，x=0.82。ηloss与Wex·Hey的关系如图9所示，误差线为10%的相对误差，可看出关系式与数据点之间拟合关系良好。

2.4 附板焊缝高度对液膜损失的影响

图10为附板焊脚高比对液膜损失的影响。可看出，在实验和模拟范围内，随着焊脚高比的不断减小，模拟计算获得的液膜损失率与实验值间的误差越来越大，本文所用计算模型的预测准确性越来越低。

图9 不同流量下损失率和Wex·Hey之间的关系Fig.9 Relationship between film loss rate and Wex·Hey under different flow rates

图10 附板焊脚高比对液膜损失的影响Fig.10 Influence of attached plate’s welding height on film loss

图11a、b为液膜冲击焊脚高比β为1.0和0.3的附板时垂直于板面方向的速度分布场，图11c、d为液膜的速度矢量场。从图11a可看到，EWF模型模拟β=1.0的附板时对速度场的模拟与实际情况较接近，能较好地模拟液膜冲击时飞溅射出的现象。图11b显示，液膜冲击附板时液体在拐角处获得较高焊脚附板工况更大的z方向速度，而由于EWF模型基于液膜内速度与壁面平行的基本假设，因此在经过拐角处时，可明显看到液膜z方向速度以极高的速度梯度减小。

从图11c、d中圆圈标记处能更加清晰地看到EWF模型对不同于壁面方向速度的修正。对于焊脚高比较大的工况，由于液膜垂直于运动方向的速度不大，EWF模型的这种修正对液膜的流动影响较小；而对于焊脚高比较小的工况，EWF模型对液膜垂直于运动方向的速度进行修正，使液膜运动方向强行贴合壁面，导致液膜的流动不符合惯性和实际作用，最后导致模型的计算结果与实验结果偏差较大。

图11 液膜冲击不同焊脚高度的附板时速度流场Fig.11 Velocity field of falling film impinging on attached plate with different welding heights

3 总结与讨论

本文对大平板上降液膜冲击附板的流动行为进行分析，建立了相关几何模型，利用商用CFD程序FLUENT建立了降液膜冲击附板行为的相关数值模拟方案并进行了数值模拟研究。经过对计算结果的分析和与相关实验结果的比较，发现大平板上降液膜冲击附板的流动行为过程主要可分为4个阶段：液膜冲击附板阶段、液膜在附板上第1次再聚合阶段、液膜脱离附板分离阶段、液膜在平板上第2次再聚合阶段。通过对4个阶段研究，主要得到以下结论。

1) 液膜在小流量下冲击附板时普遍损失率较小，而在中大流量下冲击附板有较大损失，且损失率受附板尺寸影响；损失率和韦伯数、受附板高度影响的液膜空间数之间的关系为ηloss∝Wex·Hey(0

2) FLUENT内EWF模型对于液膜冲击高焊脚的附板能实现较好的预测；而对于低焊脚的附板，EWF模型对液膜冲击附板后的速度有较大修正，导致EWF模型对液膜冲击低焊脚的附板后的损失无法很好预测。