基于OpenFOAM的含不凝气蒸汽气泡冷凝模拟

陈家帅 曾义凯 都 宇

（1.西南交通大学机械工程学院 成都 610031；2.中国核动力研究设计院中核核反应堆热工水力技术重点实验室 成都 610213）

0 引言

直接接触冷凝具有较高的传热效率，在很多工业过程中都能看到他的身影，比如天然气锅炉烟气余热回收以及沸水堆核电站抑压水池。近年来，人们对于直接接触冷凝进行大量的研究[1-3]。在天然气余热回收领域，由于气泡的冷凝行为直接影响冷凝热水器和冷凝锅炉等设备对烟气中余热的回收效率，因此在此类烟气换热器的研发设计中，直接接触冷凝一直被视作核心问题而被研究[4-6]。含不凝结气体的蒸汽气泡表面发生传热和传质现象的同时，在气泡内部还存在蒸汽-不凝结气体组成的扩散层，因此其运动行为和传热特性十分复杂，对直接接触冷凝的研究是一个难度很大的挑战。为了能够更深入的了解含有不凝结气体的蒸汽气泡直接接触冷凝现象，有必要对其冷凝行为进行深入的研究。

至今，已经有许多学者对于气泡冷凝行为进行了大量的实验研究分析。在屈晓航[7]的研究中，通过对实验数据的拟合得到了冷凝Nu 数的关联式。在Kim 和Park[8]的研究中，通过实验分析得到低压条件下过冷沸腾中纯蒸汽气泡表面的换热系数关联式。随着计算机技术的发展和数值模拟方法的改进，越来越多的研究者采用CFD 来研究流动和换热类的问题[9,10]，两相流的模拟得到了进一步的发展。很多研究者采用VOF 方法对蒸汽气泡的冷凝过程进行数值模拟，并从气泡大小、气泡形状、上升轨迹和气泡速度等方面分析气泡的行为[11,12]。Welch 和Wilson[13]、屈晓航[14]和Samkhaniani[15,16]等人的研究中检验了采用VOF 模型分析过冷液中气泡冷凝的可行性。

本文主要研究了如何利用开源软件包OpenFOAM 中的程序模拟空气-水蒸气混合气泡的直接接触冷凝过程。本文对OpenFOAM[17]中的多相流多组分求解器icoReactingMultiphaseInterfoam进行了改进以对含有不凝结气体的蒸汽气泡进行数值模拟。将数值模拟结果与实验结果进行了比对，验证了该求解器的适用性。通过该求解器，研究了不同初始条件下空气-水蒸气混合气泡的冷凝特性，并与纯蒸汽气泡的冷凝特征对比，分析不凝结气体对冷凝气泡行为的影响。

1 数值模型与验证

1.1 控制方程

在本文中，由Hirt 和Nichols 提出Volume of Fluid（VOF）界面追踪方法[18]将被用来对含有干空气的蒸汽气泡在过冷水中冷凝这个两相三组分系统中的气-液相界面进行追踪。VOF 方法的控制方程如下：

式中，αL= 1表示的是液相区域，αL= 0表示的是气相区域，0 ＜ αL＜ 1表示气相和液相交界的区域。气泡的边界位于气相和液相交界的区域。在实际工程应用中，通常将αL= 0.5的等值面作为气相和液相的交界面。

随着流体的流动，气液交界面的位置随之发生改变，既流场中气相体积分数分布应该同步发生改变，需要对相体积分数输运方程进行求解。对于存在相变的多相流系统，其对应的体积分数方程为：

式中，m&表示单位体积冷凝率，kg/(m3·s)。数值扩散的存在导致求解上述方程得到的相体积分数场中两相的交接区域过大，相界面不够尖锐。前人已提出多种可用于VOF 模型中以保证界面尖锐的方法。在OpenFOAM 中采用Weller[19]提出的在体积分数输运方程中添加人工对流项的方法，来对两相界面附近的相分数进行压缩，削弱数值耗散的影响，保证求得相界面足够锐利。这个额外的对流项为 ∇· （α1（1 - α1）Uc ），其只在两相的交接区域起作用，在纯液相区域和纯气相区域该项的值为零。加入人工对流项后体积分数输运方程变为：

Cα表示压缩系数，当压缩系数的取值在1 到4 之间时，对于大多数工况都能得出较好的结果，但是对于某些工况压缩系数的值需要大于4，本文中的压缩系数Cα取值为1。

为了简化边界条件的定义，在OpenFOAM 中动量方程的总压（P）转换为了动压（Prgh）。求解器中求解的动量方程如下：

等式左边的最后一项表示气相和液相之间的表面张力。表面张力通过连续表面应力模型[20]计算。σ表示表面张力系数，k 表示气液交界面的曲率。

由于相变伴随能量传递，需要求解能量方程。在OpenFOAM 的多相流求解器中，为了削弱相界面两侧气相和液相物性参数的不平衡性，将能量方程表示为方程（7）所示的形式，以增强数值算法的鲁棒性。

式中，HLG表示单位质量蒸汽冷凝所释放的相变潜热，J/kg。

气泡中水蒸气在干空气中的扩散，由如下方程描述：

等式右边第二项是组分扩散方程的源项，既蒸汽的冷凝量。由于混合气体中只有水蒸气发生冷凝，水蒸气的冷凝量即是气相质量的变化量。除了不凝气-蒸汽混合气体密度Gρ 是由理想气体公式计算，其余的物性参数根据蒸汽和干空气的质量分数加权求得。

为了闭合控制方程，需要采用合适的相变模型计算相界面出的质量通量。本文中采用的相变模型是Lee 模型，其对应的数学方程为：

式中，rc是一个经验系数称为传质强度因子，通常是通过和实验数据对比得到。本文通过和已有的实验数据对比将其取值为500。

1.2 模型验证

为了验证本文所建立的数值模型的合理性，本文对屈晓航等人的实验工况进行数值模拟，并将数值模拟结果与实验结果比对。对含有不凝结气体的蒸汽气泡冷凝的模拟流程如图1 所示，具体步骤如下：

图1 直接接触气泡冷凝模拟流程Fig.1 Simulation flow of direct contact bubble condensation

第一步：在三维笛卡尔坐标系下，采用blockMesh 生成几何模型和网格。如图1 中所示，创建的计算域的尺寸为30×30×60mm3气泡具有向上的初始速度。计算域的上边界设置为压力出口边界，下边界和侧面则设置为壁面边界。

第二步：采用OpenFOAM 中自带的小工具setFields 对计算域中各个物理量进行初始化。对于气相区域的蒸汽质量分数、温度和速度等物理量的初始条件根据实验中测得的数据进行设置。所要模拟的两个工况的初始条件如表1 中所示，气泡1 对应的过冷度为14K，气泡2 对应的过冷度为10K。

表1 含不凝气蒸汽气泡的初始条件Table 1 Initial conditions of bubbles air-vapor mixture bubble

第三步：为了保证结果的合理性，需要保证气相区域有足够数量的网格，并且为了避免数值扩散对相界面的过度涂抹往往需要保证界面区域的网格尺寸小于气泡直径的1/16。本文为以较小的数值消耗计算出较为尖锐的气液两相界面，在界面区域采用了自适应网格技术。

第四步：在三维笛卡尔坐标系下，采用本文改进后的icoReactingMultiphaseInterfoam 求解器对单个含有不凝结气体的蒸汽气泡的冷凝过程进行模拟计算。

图2 定性的展示了在两组不同的初始工况下，数值模拟结果和实验结果的对比，包括冷凝过程中气泡形状，大小以及穿透距离等。气泡在上升的过程中随着气泡中蒸汽的冷凝，其体积减小、外形也随之变化。如图2 中所示，对于直径稍小的气泡，其外形由圆球形变为半球形接着再转变为椭圆形，对于初始直径较大的气泡，其外形由初始时的球形，随后下表面内凹，接着转变为半球形，最后转变为椭球形。这是由于，随着气泡直径的增大，表面张力对气泡变形的影响力减弱。而最后都变为椭球形则是由于受力的原因。由于数值模拟的初始条件和实验的真实条件有区别，故而数值模拟的结果和实验结果在气泡外形、气泡穿透距离等方面有细微的差异。但是鉴于实验中初始条件的不确定性，可以认为本文的数值模型计算得到的结果和实验结果是一致的。

图2 数值模拟和实验的结果对比Fig.2 Comparison of numerical simulation and experimental results

图3 定量展示了气泡的体积随时间变化的规律。数值模型计算出的体积变化曲线和实验数据基本吻合。对于含不凝结气体的蒸汽气泡，随着蒸汽在气泡表面处凝结，气泡中不凝气体含量增高，组分扩散层中的浓度梯度越发明显，传热和传质阻力也越发的大，最终会达到稳定的状态，气泡的体积基本维持稳定。由图2 和图3 展示的结果可看出，本文建立的数值模型中的质量源项，贴合实际相变过程。

图3 气泡体积变化对比Fig.3 Comparison of bubble volume changes

图4 为图2 中冷凝气泡的横截面处温度分布和蒸汽质量分数分布。在图4 中红实线外的区域，红实线和红虚线之间的区域，以及红虚线以内的区域分别表示过冷水区，蒸汽-干空气的组分扩散层以及气泡内部蒸汽聚集的区域。如图所示，含不凝结气体的蒸汽气泡中温度连续分布，而纯蒸汽气泡在冷凝过程中的温度分布在相界面处有明显的间断。这是由于在蒸汽-干空气混合的气泡中，从气泡中心到气泡边界，水蒸气的质量分数不断的下降，水蒸气的分压也随之下降，进而导致对应的饱和温度也逐渐降低，且低于气泡中心的饱和温度。

图4 含不凝结气体蒸汽气泡的温度分布和质量分数分布Fig.4 Temperature distribution and mass fraction distribution of air-vapor bubble

2 几何模型

为了深刻理解含不凝气蒸汽气泡冷凝过程中各种因素的影响，本文在直角坐标系下，对图5 中所示的二维气泡进行了数值模拟。在计算域底部释放含不凝气蒸汽气泡，研究各种参数对其冷凝行为的影响。情形一，如图5（a）所示，研究不凝结气体含量、气泡初始直径、流场速度以及过冷度等因素对气泡冷凝行为的影响。情形二，如图5（b）和图5（c）所示，将两个气泡在计算域底部并排或者竖排放置，研究相邻气泡对冷凝行为的影响。

图5 干空气-蒸汽混合气泡冷凝模拟示意图Fig.5 Schematic diagrams of air-vapor mixture bubble condensation simulations

为了求得不受网格密度影响的数值解，在五组不同的网格密度下，对图5（a）中的工况进行了数值计算。图6 展示了在不同的网格密度下，计算所得的气泡等效直径变化曲线以及气泡在同一时刻的外形。通过对五组网格计算结果的比较，可以看出250×500 所对应的网格密度计算出的结果是较为合理的，因此其被用于本文的后续计算。

图6 网格无关性验证Fig.6 Obtaining a grid independence solution

3 结果讨论

3.1 单个含不凝气蒸汽气泡运动特征分析

为了研究流场速度以及不凝气含量对于迁移轨迹的影响，在过冷度为15K 的条件下，对初始直径为8mm 的具有不同的初始蒸汽质量的气泡在不同的流场速度下模拟。如图7 中所示，图中对五种工况进行了计算，（a）、（b）、（c）、（d）、（e）对应的初始蒸汽质量和流场速度分别为0.9、0.7、0.7、0.6、0.6 和0m/s、0m/s、0.1m/s、0.1m/s、0.2m/s。对每个工况每间隔0.01s 提取一次气泡的边界坐标。对于含不凝结气体的蒸汽气泡，在相同的过冷度下，不凝结气体的含量主要影响气泡体积基本稳定后气泡的变形，对于气泡迁移轨迹影响甚微。流场速度对于气泡的迁移轨迹影响较大，但是增大流场速度不能够加快气泡中的蒸汽凝结。

图7 气泡迁移轨迹Fig.7 Bubble migration trajectory

3.2 单个含不凝气蒸汽气泡冷凝行为分析

对于含不凝气蒸汽气泡直接接触冷凝，需要注意其完成冷凝需要的时间，本文将气泡体积基本不变视为完成冷凝，此时气泡对应的等效直径称为终端等效直径。冷凝时间受到气泡初始尺寸，气泡中初始蒸汽质量分数以及过冷度等因素的影响。由图8（a）可知，随着直径的增大，气泡的冷凝时间增大，但是气泡的冷凝时间与气泡初始直径之间不呈现明显的线性关系，这点与纯蒸汽气泡不同。随着过冷度的增大，冷凝时间随之增大，并且随着气泡直径的增大，在一定的过冷度范围内，气泡的冷凝时间分布范围减小。如图8（b）所示，终端等效直径随过冷度的增大而减小，随气泡初始直径的增大而增大，并且与气泡的初始直径呈线性关系。此外，在一定的过冷度范围内，随着气泡直径的增大，其对应的终端等效直径的范围也随之扩大。已有文献指出，纯蒸汽气泡的冷凝时间与气泡的初始直径在一定程度上呈现线性关系，而对于含不凝气蒸汽气泡，不同初始直径或者不同初始蒸汽质量分数的气泡，其在冷凝过程中形成的不凝气-蒸汽组分扩散层也随之出现差异，进而对蒸汽冷凝阻碍程度不同，导致气泡的冷凝时间与气泡尺寸和不凝结气体的含量呈非线性关系。

图8 过冷度对含不凝气蒸汽气泡冷凝的影响Fig.8 Influence of subcooling degree on air-vapor bubble condensation

图9（a）中展示了过冷度为15K，初始直径为6mm 的不同蒸汽质量分数的气泡在冷凝过程中的体积变化规律。图9（b）则展示了不凝气体含量以及气泡初始直径对于气泡终端等效直径的影响。随着气泡中不凝气体含量的增加，气泡的体积衰减速率减小，终端等效直径增大。并且，随着气泡初始直径的增大，在一定的过冷度范围内，气泡的终端等效直径的分布范围增大。不凝气体含量大的蒸汽气泡随着初始直径的增大，对应的终端等效直径增长的更快。

图9 不凝气含量对冷凝的影响Fig.9 Influence of non-condensable gas content on condensation

3.3 多气泡分析

图10 展现了两个并列和竖排的含不凝气蒸汽气泡同时冷凝时，相互之间的影响。图中列出了四中工况的模拟结果。情形一，模拟竖排的气泡对的冷凝过程，气泡中心之间的间距为1.1D0，如图5（c）所示。情形二和情形三分别是情形一中两个气泡的冷凝过程的单独模拟，以和情形作对照，探究相邻气泡对冷凝过程的影响。情形四是模拟并列的气泡对的冷凝过程，气泡中心的间距仍为1.1D0，如图5（b）所示。情形二中气泡在0.05s 对应的等效直径为5.223mm，情形三中气泡在0.05s 对应的等效直径为5.175mm。这是由于情形三中气泡所处的深度较情形二深，其气泡内部的压力较大，进而气泡内蒸汽对应的饱和温度越高，冷凝加快。另外，综合对比情形一、情形二和情形三，位于下侧的气泡的蒸汽冷凝率要低于其单气泡时的冷凝率。当气泡对竖直放置时，且气泡间隔较小时，下部气泡位于上部气泡的热边界层中，气泡中的蒸汽温度与气泡外侧流体温度的差值减小，蒸汽冷凝减缓。由情形四可知，当两个含不凝气蒸汽气泡并排放置时，在冷凝过程中受到来自过冷液的力的作用将向两侧旋转分离。

图10 多气泡冷凝过程Fig.10 Multi-bubble condensation process

4 结论

（1）含不凝气蒸汽气泡在冷凝过程中的变形过程为：球形，月牙形，半球形，椭圆形，当气泡的体积不再变化后，气泡的形状会发生震荡。对于初始直径较小的含不凝气蒸汽气泡，其在冷凝过程中不会出现下表面内凹。气泡边界处存在不凝气-蒸汽组分扩散层，存在明显的温度梯度。

（2）气泡中不凝气含量决定了气泡完成冷凝后的体积，进而直接影响气泡的变形，对气泡的移动轨迹影响相对较小。而流场速度主要影响气泡的迁移轨迹。

（3）含不凝气蒸汽气泡表面存在蒸汽冷凝的时长随过冷度的增大而增大，比如对于初始蒸汽质量分数为0.7，初始直径为6mm 的气泡，当过冷度为5K 时，冷凝时长为50ms，而当过冷度为20K时，冷凝时长为60ms。此外，含不凝结气体的蒸汽气泡的冷凝时长与气泡初始直径以及过冷度均呈现非线性关系，而终端等效直径则是与气泡的初始直径呈线性正相关，并随过冷度的增大而减小。初始不凝气含量较多的气泡对应的蒸汽冷凝率较小，其对应的终端等效直径随气泡初始直径的增大而增长较快。

（4）当两个含不凝气蒸汽气泡竖直排列，且间隔较小时，上面的气泡加热的液体位于下面的气泡的上部，进而对下面气泡的冷凝会起到一个抑制，故而对于利用直接接触气泡冷凝来进行换热的换热器，要保证生成的气泡之间有足够的间隔。