基于微液层模型的单汽泡生长数值模拟研究

基于微液层模型的单汽泡生长数值模拟研究
王烨;蔡杰进
【摘要】Nucleate boiling heat transfer plays an important role in the heat and mass transfer. The mechanism and the heat and mass transfer of nucleate boiling are still the focus of the study. With the improvement of experimental methods, the microlayer model was confirmed and received extensive attention. The diffusion interface method was used to track the phase interface. Through the research on the heat and mass transfer in microlayer, the correlation of the microlayer thickness with the heat flux and gasification rate was established. The partial differential equation of the microlayer model was introduced at the bottom of the bubble. The effect of the change of the microlayer on the bubble evolution and the temperature distribution of the heating wall was studied in the condition of nucleate boiling. The results show that the bubble evolution process and the heating wall temperature distribution are in good agreement with those of the experiment, which indicates the accuracy of the model. Furthermore, the variation law of the microlayer, the dry zone and the bottom radius of the bubble and the heating wall temperature distribution were analyzed.%核态沸腾换热在传热传质方面有着重要的作用,其发生机理和传热传质过程仍是研究的重点.随着实验手段的提高,微液层模型得到了广泛的关注.通过对微液层中传热传质的分析,建立了微液层厚度与热流密度和气化率之间的关系.利用界面扩散法对汽液相界面进行追踪,并在汽泡与加热壁面之间构建微液层模型,研究在核态沸腾条件下,微液层的变化对汽泡生长和加热壁面温度分布的影响.结果
表明,数值模拟得到的汽泡生长过程和加热壁面温度分布与实验结果吻合得很好,初步验证了模型的正确性.并通过数值模拟,进一步分析了汽泡生长过程中微液层、干性区域和汽泡底部半径的变化规律以及壁面温度的分布情况.
【期刊名称】《原子能科学技术》
【年(卷),期】2018(052)004
【总页数】7页(P600-606)
【关键词】微液层模型;扩散界面法;核态沸腾;汽泡生长;数值模拟
【作者】王烨;蔡杰进
【作者单位】华南理工大学电力学院,广东广州 510640;华南理工大学电力学院,广东广州 510640
【正文语种】中文
【中图分类】TK284
由于沸腾换热可利用汽化潜热，小温差下有着很强的换热能力，在工业生产中有着广泛的应用。

在沸水堆和压水堆的蒸汽发生器二次侧中，核态沸腾换热为最主要、安全、高效的换热方式。

当处于核态沸腾换热阶段时，能保证良好的换热效果从而保证设备的安全稳定运行。

虽然目前对于工程应用，已建立了核态沸腾的相关经验公式和一般判别条件[1]，但由于核态沸腾汽泡产生速度快，发生尺度小，在很长的一段时间中，缺乏对核态沸腾过程进行机理性的研究。

随着实验手段的发展，高速相机和红外热像仪极大地提高了图像和温度测量的时间分辨率和空间分辨率。

实验结果表明微液层模型与实验结果吻合，使得微液层模型
成为国内外学者研究核态沸腾的主流模型[2-4]。

由于微液层理论较为复杂，如何
结合理论和实验，验证模型的正确性，数值模拟便是一种很好的方法。

相比于VOF法和Level-Set法，扩散界面法在界面处的处理使其更容易与其他添
加的偏微分方程和源项兼容，具有更好的收敛性和质量守恒特征。

扩散界面法由van de Waals[5]于1895年提出，Yue等[6-7]对扩散界面法在两相流中的应用进行了介绍，并利用界面扩散法对两相流的一些基本问题进行了模拟。

本文利用扩散界面法对相界面进行捕捉，并引入微液层模型，对核态沸腾条件下单汽泡的生长过程及加热壁面的温度分布进行数值模拟，并与实验结果进行比较。

1 数学物理模型
核态沸腾条件下的汽泡生长过程涉及多个控制方程。

首先最基本的是流动和传热控制方程；其次是相界面捕捉方程，最后是微液层内传热和传质过程的控制方程。

1.1 流动传热方程
核态沸腾中，与无相变的流体流动不同的是，在发生相变过程时，由于两相间的密度差异较大，由密度差所引起的速度、压力变化和由汽化潜热所引起的传热变化可通过在原有的动量方程、连续性方程和能量方程中添加适当的源项进行处理。

改写后的连续性方程、动量方程和能量方程如下：
(1)
(2)
(3)
式中：u为速度矢量；t为时间为汽化率；ξ为界面函数；ρ为密度；I为单位向量；
F为体积力；T为温度；cp为比定压热容；k为导热系数；L为汽化潜热。

1.2 扩散界面法
两相流数值模拟中，相界面追踪的方法尤为重要，是保证相界面演变准确、计算收敛的重要手段。

目前常见的界面追踪法有VOF、Level-Set、MPS和扩散界面法等[8-11]。

扩散界面法以相界面处理光滑，收敛性、守恒性好，易于添加和修改界面控制源项等优势，被逐步应用于两相流数值模拟研究[12-14]。

相变量的演化由Cahn-Hilliard方程控制，可分解成两个二阶偏微分方程：
(4)
(5)
(6)
式中：φ为相位量，用于标识不同区域；∂f/∂φ为外部自由能；εpf为界面厚度；λ为混合能量密度；σ为表面张力系数；χ为移动调节参数；γ为迁移率。

1.3 微液层模型
微液层模型认为，在核态沸腾产生汽泡的初始阶段，汽泡的底部存在一层非常薄的微液层(图1)。

由于受到壁面的加热，微液层内发生相变使汽泡不断生长，微液层也会随之变薄直至变为干性区域，且微液层内的传热传质对汽泡的生长和加热壁面的温度分布具有重要的影响。

Jung等[15-16]的实验获得了汽泡生长过程的形态和汽泡底部壁面的温度分布，该实验结果与微液层模型吻合较好，也证实了微液层模型的正确性。

图中Tsat为饱和温度，mml为微液层向汽泡的相变质量转换量，qml为微液层向汽泡内的热流密度，Twall为加热壁面温度，qbk为液层向汽泡的热流密度，mbk为液层向汽泡的相变质量转换量，Tbk为液层温度。

图1 微液层示意图Fig.1 Scheme of microlayer
由于微液层的厚度只有几μm，与汽泡生长的尺度相比十分微小，不宜通过界面追踪的方法直接进行模拟。

所以通过建立汽泡底部微液层的传热传质模型，并利用偏微分方程耦合到整个计算过程中是一种简便且有效的方法。

Utaka等[17]的实验给出了微液层模型的初始厚度，即：
(7)
式中：为不同位置微液层的初始厚度，m；r为汽泡与壁面的接触半径，m。

汽泡内的温度认为保持饱和温度Tsat，在微液层发生的热量传递可表示为：
(8)
式中：qml为通过微液层的热流密度；Tw为加热壁面温度；Ri为由加热壁面到汽泡内部的等效热阻。

由于微液层内的温度也处于饱和温度，在确定微液层的热流密度后，便可确定微液层内液相向汽相的质量转换率
(9)
由式(9)可推出微液层随时间变化的偏微分方程：
(10)
由式(8)～(10)可建立微液层内的传热传质模型，并将相应的方程通过数值方法离散耦合进入模型中便可计算含微液层模型的单汽泡核态沸腾生长过程。

但需要注意的是，这些方程只适用于微液层内，所以需要加入相应的限制和适用条件。

1) 微液层只存在于汽泡与加热壁面之间，在汽泡外部不存在微液层。

(11)
2) 当汽泡内微液层的变化降低至0时，便不再发生传热传质过程，且微液层也始
终保持为0(即干性区域)。

2 数值模拟计算
图2 计算域(a)和网格(b)示意图Fig.2 Scheme of calculation domain (a) and mesh (b)
计算模型采用2D轴对称模型(图2)。

计算域包括两相流区域和固体传热区域。

图
2中底部区域为固体传热区域，上部区域为两相流区域，左侧点划线位置为对称轴。

为保证计算结果不受壁面和出口作用影响，计算域需保证足够高度及宽度。

计算域高度及宽度均取15 mm，加热板厚度取0.1 mm。

2.1 计算条件与网格
为验证模型计算的准确性，取文献[15]中的实验条件作为计算边界条件，并与实验结果进行对比。

计算边界条件为：顶部压力出口边界，竖直壁面为无滑移边界，底部壁面为润湿壁面，表现接触角θa取自文献[3]，如图3所示，其中ΔTi为过热度。

加热壁面初始温度为393.15 K，给定热流密度为53 kW/m2。

图3 接触角随生长时间的变化Fig.3 Change of contact angle with evolution time
在两相流计算域中，既存在相界面这样狭小的区域，也有不发生相变的单一相区域，因此网格需在汽泡会生长的区域进行细化。

细化的区域为r从0至3 mm、z从0到6 mm的区域，以分裂最长边的方法细化5次。

图4示出最大网格尺寸分别为0.25、0.3、0.5 mm的网格划分策略下，汽泡等效直径随生长过程的变化。

由图
4可看出，当网格较粗时，汽泡等效直径更大，当网格为0.3 mm和0.25 mm时，网格影响已不明显，综合考虑计算效率，以下数值模拟均在网格最大尺寸为0.3
mm下进行。

图4 不同网格下等效半径随汽泡生长的变化Fig.4 Equivalent radius with bubble evolution under different mesh generations
2.2 结果与验证
本文数值模拟结果与实验结果[15]的对比列于表1。

可看出，汽泡的生长周期和脱离半径的模拟结果与实验结果相差较小，底部加热壁面温度场与实验结果相似，并与微液层模型预测相一致，即随着汽泡中心微液层的蒸发，微液层区域变为干性区域，加热壁面底部会呈现环形的温度分布。

表1 数值模拟与文献[15]实验结果比较Table 1 Results comparison between numerical simulation and literature [15] experiment项目文献[15]实验数值模拟汽泡生长周期，ms1518汽泡脱离直径，mm474675ms时壁面温度，K汽泡脱离形态
3 讨论与分析
利用已验证的模型，对单汽泡生长过程进行数值模拟，得到了汽泡在不同生长阶段下的形态，并分析了汽泡生长过程中汽泡底部半径、微液层厚度及加热壁面温度的变化。

图5为壁面汽泡在生长过程中不同时刻的形态变化。

由图5可见，在汽泡生长的初始阶段，汽泡主要会沿着壁面生长。

随着汽泡体积的增大，汽泡开始逐步上浮，汽泡与壁面的接触面积开始缩小，汽泡形状由半球形逐渐变化为倒水滴形态。

最后当汽泡脱离壁面时，受浮力作用不断上升。

图5 不同生长阶段汽泡形态Fig.5 Bubble shape at different evolution stages 图6为7.5 ms时汽泡底部微液层在各位置处的厚度。

内圈圆形为干性区域即微液层厚度为0的界线，外圈圆形为汽泡底部半径界线。

由图6可看出，汽泡在生长
过程中，中心的微液层首先蒸发变薄至0，成为干性区域。

随着汽泡底部半径的生长，微液层也会出现，且厚度沿半径增加。

R1表示汽泡底部半径长度，R2表示干性区域半径长度，ΔR表示微液层环形长度。

图7为不同区域长度随着汽泡生长的变化过程。

可看出，在初始阶段，汽泡中心的微液层快速蒸发，转变为干性区域，干性区域的范围快速增大。

在汽泡的发展阶段，汽泡底部半径先缓慢增大后缓慢降低，微液层的变化与汽泡底部半径的变化相似，而干性区域保持缓慢增大。

在汽泡收缩阶段，汽泡底部半径和微液层范围快速减少，而干性区域保持增加。

微液层范围降低为0，意味着汽泡底部全部成为干性区域，最后干性区域随着汽泡半径的减小而减小，直到降低为0，即汽泡完全脱离壁面。

图6 汽泡底部各位置微液层厚度Fig.6 Microlayer thickness at bottom of bubble
图7 汽泡底部半径、微液层及干性区域随时间的变化Fig.7 Radius of bubble, microlayer and dry region vs evolution time
图8为不同时刻微液层厚度与汽泡底部半径的关系。

从图8可看出，在汽泡生长的初始阶段，汽泡底部半径会不断生长，即随着时间的推移，汽泡底部半径变大。

对于微液层，在汽泡生长的初始阶段，汽泡中心的微液层会快速汽化使微液层变薄直至消失，其余位置的微液层也会由于相变的作用变薄，但在汽泡边界的微液层会随着汽泡半径的生长不断产生新的微液层。

在汽泡生长的中期，受到浮力、表面张力和接触角的共同作用，汽泡会开始浮升，汽泡底部半径开始收缩，且随着汽泡底部半径的收缩，边界处的微液层虽未完全汽化，但也会消失，不再产生相变。

在汽泡生长的末期，汽泡继续浮升，汽泡底部半径快速收缩，并由于微液层已完全相变为汽相，在最后阶段，汽泡底部完全为干性区域。

图8 微液层厚度与汽泡底部半径随时间的变化Fig.8 Microlayer thickness and
bubble radius at different evolution time
图9 加热壁面温度分布随时间的变化Fig.9 Temperature distribution on heating wall at different evolution time
图9为汽泡底部加热壁面的温度分布随时间的变化。

在汽泡开始生长的阶段，汽
泡底部半径和微液层范围较小，温度只在很小区域有较大的变化。

随着汽泡的生长，汽泡底部半径变大，微液层范围也相应变大，温度降低的区域也增大。

在微液层内，由于汽泡中心的微液层最薄，在很短的时间便完全消失，不再产生相变，也不会有传热传质的过程发生，而在中心向外的一定区域微液层存在的时间更多，所以会有更多的热量通过相变换热被带走，而在汽泡底部半径外的区域受相变换热的影响很小，所以一直维持较高的温度。

由于在不同的区域壁面的换热能力不同，所以使得壁面的温度分布呈现汽泡中心高，之后降低，汽泡边缘处再升高的分布。

壁面的最低温度将会随着汽泡底部半径的生长而变化，且由于壁面热传导的影响，温度分布将会逐渐趋于平缓。

4 结论
本文通过建立汽泡核态沸腾的微液层模型，利用扩散界面法对相界面进行捕捉，对单汽泡生长过程进行了数值模拟研究。

数值模拟结果与实验结果在汽泡生长周期和脱离半径上吻合较好，脱离形状与壁面温度分布规律相同。

通过分析汽泡微液层在整个生长周期中的变化，得到以下结论。

1) 在汽泡生长的初始阶段，汽泡底部中心的微液层快速蒸发，使汽泡底部中心区
域快速转变为干性区域，且伴随大量的传热传质现象的发生，不仅使汽泡快速长大，同时也使加热壁面温度快速下降。

2) 在汽泡的发展阶段，汽泡中心的干性区域持续扩大，微液层范围也随着汽泡底
部半径的增大而增大。

在汽泡底部中心区域，由于微液层蒸发转为干性区域，所以中心温度较两侧温度更高，底部温度分布呈现圆环形。

当汽泡长大至一定程度后，
受浮力、表面张力和接触角等因素的共同作用，汽泡向上浮升，底部开始收缩。

3) 在汽泡的收缩阶段，微液层会受到汽泡底部半径的缩小和汽泡中心部分干性区域的增大的影响而快速缩小，直到消失。

此时，汽泡底部完全是干性区域。

汽泡底部半径会继续缩小直至汽泡脱离壁面。

壁面温度也会在热传导的作用下趋于一致。

本文所建立的微液层数值模拟模型能对单汽泡生长过程进行较好的模拟，初步验证了模型的正确性，对单汽泡核态沸腾的进一步研究也正在展开。

[1] 李美琳，林萌，杨燕华，等. 应用于反应堆热工水力程序的核态沸腾传热关系式评价[J]. 核科学与工程，2015,35(1)：25-31.
LI Meilin, LIN Meng, YANG Yanhua, et al. Assessment of nucleate boiling correlations applied to thermal-hydraulic code of reactor[J]. Nuclear Science and Engineering, 2015, 35(1): 25-31(in Chinese).
[2] 高明，章立新，郑平，等. 沸腾汽泡微液层变化规律的实验研究[J]. 工程热物理学报，2013,34(5)：931-934.
GAO Ming. ZHANG Lixin, ZHENG Ping, et al. Experimental study of the microlayer beneath nucleation bubble[J]. Journal of Engineering Thermophysics, 2013, 34(5): 931-934(in Chinese).
[3] CHEN Z, UTAKA Y. On heat transfer and evaporation characteristics in the growth process of a bubble with microlayer structure during nucleate boiling[J]. Int J Heat Mass Transfer, 2015, 81: 750-759.
[4] MA H B, CHENG P, BORGMEYER B, et al. Fluid flow and heat transfer in the evaporating thin film region[J]. Microfluidics and Nanofluidics, 2008, 4(3): 237-243.
[5] van der WAALS J D. The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density[J]. Journal of Statistical
Physics, 1979, 20(2): 200-244.
[6] YUE P, ZHOU C, FENG J J, et al. Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[J]. Journal of Computational Physics, 2006, 219(1): 47-67.
[7] YUE P, FENG J J, LIU C, et al. A diffuse-interface method for simulating two-phase flows of complex fluids[J]. Journal of Fluid Mechanics, 1999, 515: 293-317.
[8] CHEN R, TIAN W, SU G H, et al. Numerical investigation on bubble dynamics during flow boiling using moving particle semi-implicit
method[J]. Nucl Eng Des, 2010, 240(11): 3 830-3 840.
[9] 赵云淦，牛风雷，单祖华. 气泡在液态铅铋合金内上升行为的数值模拟[J]. 原子能科学技术，2015，49(增刊)：278-282.
ZHAO Yungan, NIU Fenglei, SHAN Zuhua. Numerical simulation on bubble rising behavior in liquid lead-bismuth alloy[J]. Atomic Energy Science and Technology, 2015, 49(Suppl.): 278-282(in Chinese).
[10] WU J, DHIR V K, QIAN J. Numerical simulation of subcooled nucleate boiling by coupling level-set method with moving-mesh method[J]. Numerical Heat Transfer Fundamentals, 2007, 51(6): 535-563.
[11] 左娟莉，田文喜，秋穗正，等. 液态金属内单个气泡上升行为的MPS法数值模拟[J]. 原子能科学技术，2011，45(12)：1 449-1 455.
ZUO Juanli, TIAN Wenxi, QIU Suizheng, et al. Numerical simulation on single bubble rising behavior in liquid metal using moving particle semi-implicit method[J]. Atomic Energy Science and Technology, 2011, 45(12): 1 449-1 455(in Chinese).
[12] 王烨，蔡杰进. 基于扩散界面法的较大气泡上升过程数值模拟[J]. 原子能科学技术，2017，51(2)：269-274.
WANG Ye, CAI Jiejin. Numerical simulation of large bubble rising process based on diffuse interface method[J]. Atomic Energy Science and Technology, 2017, 51(2): 269-274(in Chinese).
[13] WANG Y, CAI J, LI Q. Numerical simulation of large bubble rising behavior using diffuse interface method[J]. International Journal of Hydrogen Energy, 2017(in press).
[14] WANG Y, CAI J. Numerical investigation on bubble evolution during nucleate boiling using diffuse interface method[J]. Int J Heat Mass Transfer, 2017, 112: 28-38.
[15] JUNG S, KIM H. An experimental method to simultaneously measure the dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface[J]. Int J Heat Mass Transfer, 2014, 73: 365-375.
[16] JUNG S, KIM H. An experimental study on heat transfer mechanisms in the microlayer using integrated total reflection[J]. Laser Interferometry and Infrared Thermometry Technique, 2015, 36(12): 1 002-1 012.
[17] UTAKA Y, KASHIWABARA Y, OZAKI M. Microlayer structure in nucleate boiling of water and ethanol at atmospheric pressure[J]. Int J Heat Mass Transfer, 2013, 57(1): 222-230.。