, 杨星团1, 屠基元1,2, 姜胜耀11. 清华大学 核能与新能源技术研究院, 先进核能技术协同创新中心, 先进反应堆工程与安全教育部重点实验室, 北京 100084, 中国;
2. 皇家墨尔本理工大学 工程学院, 墨尔本 VIC 3083, 澳大利亚
收稿日期：2022-12-06
基金项目：国家自然科学基金项目（51576211）
作者简介：刘倩(1996-), 女, 博士研究生
通讯作者：桂南, 副教授, E-mail: guinan@mail.tsinghua.edu.cn
摘要：在依靠自然循环驱动的小型模块化反应堆主回路以及反应堆非能动安全系统中，冷凝换热是热交换过程的重要一环，因此对冷凝换热过程的深入研究和分析对提升换热效率、保障堆芯安全至关重要。该文基于格子Boltzmann方法，采用伪势模型，模拟研究了二维通道内静止饱和蒸汽在凝结过程中的流动和换热特性。结果表明：蒸汽冷凝会自发驱动蒸汽流动，蒸汽质量流速与通道宽度和凝结过程中的热流流量有关，保持壁面温度和通道宽度恒定状态下，液膜发展阶段壁面热流较大时，蒸汽质量流速增长较快；通道较窄时，入口处蒸汽平均质量流速初期增速较快但迅速达到稳态，通道宽度为150时的稳态入口平均质量流速约是通道宽度为500时的80.0%。对接触角的分析表明：亲水壁面上的液膜厚度受壁面亲水程度影响较小，壁面接触角为51°时出口处液膜厚度与接触角为72°时的相等。普通疏水壁面上珠状凝结难以维持，被液膜覆盖后相较于亲水壁面传热速率较慢，液膜滑移出计算域之前，壁面接触角为127°时壁面平均热流密度最大值约是接触角为51°时的75.8%，并随液膜滑移逐渐降低，但液膜受重力去除后再形成的过程能在一定范围内强化传热速率。
关键词：冷凝传热格子Boltzmann方法蒸汽流速接触角
Numerical simulation of saturated steam condensation heat exchange in a vertical channel
LIU Qian1, GUI Nan1
, YANG Xingtuan1, TU Jiyuan1,2, JIANG Shengyao11. Key Laboratory of Advanced Reactor Engineering and Safety of Ministry of Education, Collaborative Innovation Center of Advanced Nuclear Energy Technology, Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China;
2. School of Engineering, RMIT University, Melbourne VIC 3083, Australia
Abstract: [Objective] Condensing heat exchange is a crucial process in the primary circuit of small modular reactors and passive safety systems that rely on natural circulation as the driving force. With the higher requirements for heat exchange efficiency and reactor safety, in-depth research and an understanding of the condensation heat exchange process are needed. Therefore, numerical simulations of the condensing heat exchange process have attracted increasingly more interest. However, due to the complex phase change, the condensing heat exchange process is difficult to model using analytical equations. Traditional numerical simulation methods use the empirical equations summarized in experiments, and their universalities are controversial. In contrast, the lattice Boltzmann method is a mesoscopic-level numerical simulation method that tracks particle clusters and uses probability density functions to describe their distribution, resulting in a simple and clear structure that appropriately ignores the details of molecular motion. Moreover, it allows direct iterative solving of the probability density distribution function without relying on empirical equations. In previous studies, the feasibility of using the lattice Boltzmann pseudopotential model in condensation process simulation was verified. Subsequently, numerous researchers have used this model to analyze the condensation mechanism. [Methods] This study is based on the lattice Boltzmann method and uses a dual distribution function to simulate the condensation process of stationary saturated vapor within a vertical channel. To analyze the fluid flow characteristics, a pseudopotential model is used to simulate the density field variations during the vapor condensation. Additionally, a temperature distribution function is employed to simulate the temperature field changes during the vapor condensation, allowing for an examination of heat transfer efficiency. Throughout the simulation, we analyze the effects of channel width and the hydrophilicity and hydrophobicity of wall conditions on the condensate flow and heat transfer rate. [Results] The results showed that: 1) When saturated vapor encountered a hydrophilic wall, it first condensed to form a thin liquid film covering the entire wall surface and then formed a steady liquid film from the top of the vertical channel, gradually expanding downward. Due to the pressure difference caused by the vapor condensation, the saturated vapor flowed down into the channel from the inlet at the top of the channel. 2) Under hydrophilic wall conditions, decreasing the channel width from 500 to 150 decreased the steady-state average mass flow rate at the inlet by approximately 20% and decreased the steady-state average heat flux density on the wall by approximately 6.5%. 3) The simulation results under different hydrophobic and hydrophilic characteristics were consistent with the theoretical analysis, indicating that the stronger the wall hydrophobicity was, the later the starting time of droplet nucleation and the lower the starting point of the vertical liquid film. On the ordinary hydrophobic wall surface, the droplet condensation was difficult to sustain, and after this surface was covered by a liquid film, the heat transfer rate was slower compared to the hydrophilic wall surface. Before the liquid film slipped out of the computational domain, the maximum average wall heat flux at an angle of 127° was approximately 75.8% of that at an angle of 51°. [Conclusions] The lattice Boltzmann pseudopotential model can simulate the condensation process of stationary saturated vapor within a vertical channel. During the simulation process, the effects of channel width and the hydrophilicity and hydrophobicity of wall conditions on condensate flow and heat flux density are important. In general, wider channel widths lead to higher wall heat flux density, and a higher inlet mass flow rate of steam is achieved at a steady state for saturated vapor initially in a stationary state within the channel. However, ordinary hydrophobic wall surfaces cannot sustain droplet condensation and do not demonstrate enhanced heat transfer. These findings have certain reference values for designing and optimizing heat transfer systems involving condensation in vertical channels.
Key words: condensation heat transferlattice Boltzmann methodsteam velocitycontact angle
国家和社会发展离不开能源电力的支持，在低碳减排的发展背景下，清洁、高效的核能在生产和生活中发挥着日益重要的作用，提高换热效率、保障堆芯安全是核电技术发展历程中持续追求的目标。目前的第3代、第4代反应堆设计中广泛运用了基于自然循环技术的非能动系统，这些系统不借助机械动力，仅依靠冷热端密度差和高度差实现自然驱动，是事故工况下保障堆芯热量有效导出的关键[1-2]。除安全保障系统外，一些小型模块化反应堆(small modular reactor，SMR)不需要主泵，仅依靠自然循环即可实现全功率范围内堆芯冷却，简化反应堆结构的同时提高了反应堆的固有安全性[3]。竖直通道内冷凝换热是反应堆内自然循环回路换热方式中主要的方式之一，以非能动余热排出系统为例，李常伟等[1]的实验表明：在直管式冷凝换热器的换热环节中，冷凝换热的热阻约占总热阻的60%，且随着蒸汽压力的升高，冷凝换热的热阻呈现增大趋势，因此对竖直通道内冷凝换热过程的深入分析，对提升换热效率和保障堆芯安全至关重要[4]。
由于相变过程十分复杂，传统的数值计算方法主要依赖实验总结的经验关系式实现宏观物理量的求解，计算精度和适用性受限于经验关系式的准确性和适用范围[5-6]。格子Boltzmann方法作为一种介观层面的数值模拟方法，从粒子出发，构造简单的同时适当忽略分子运动细节，不依赖连续介质假设和经验关系式而直接进行求解，在处理多相、多尺度问题上计算效率和准确度都有明显的优势[7]。目前已经提出了许多应用于多相流的格子Boltzmann模型，大致分为以下几类[7]：颜色模型、伪势模型、自由能模型和相空间格子Boltzmann模型等，其中伪势模型因其概念清晰、易于执行的优势在多相方法中备受关注。Liu等[8]利用格子Boltzmann伪势模型模拟了稳态条件下层流区壁面膜状凝结，模拟结果与经典的解析解吻合程度较好，验证了利用格子Boltzmann伪势模型模拟冷凝过程的可行性，自此之后许多****利用格子Boltzmann模型对冷凝过程进行模拟分析。Mukherjee等[9]利用改进的伪势模型模拟了由微纳结构实现的超疏水表面，分析了具有不同微纳结构的水平表面上单液滴的生长模式、传热效率和成核时间等特征。Liu等[10]分析了过冷壁面的表面形态为锥形、三角形等微观结构时对于液滴冷凝效率的影响。Cai等[11]采用分层网格的设计将液滴的停留时间和半径分别减少了18%和17%。以上研究主要关注于滴状凝结过程的机理研究及传热优化。
本文将采用格子Boltzmann的数值模拟方法针对竖直通道内静止饱和蒸汽受冷凝驱动的流动特性和传热效率进行研究，分析通道直径和亲疏水特性对换热过程的影响。
1 格子Boltzmann方法模型模拟采用二维和9个离散速度方向(简称D2Q9)的双分布函数分别模拟流场演化的密度场和热量传递的温度场，详细介绍如下。
1.1 流场演化和伪势模型利用粒子分布函数来表征整个空间内的密度场，在模拟过程中粒子的运动被拆分为碰撞和迁移2个部分，其演化方程为
| $\begin{gathered}f_i\left(\boldsymbol{x}+\boldsymbol{e}_i \delta_t, t+\delta_t\right)=f_i(\boldsymbol{x}, t)-\frac{1}{\tau}\left(f_i(\boldsymbol{x}, t)-\right. \\\left.f_i^{\mathrm{eq}}(\boldsymbol{x}, t)\right) .\end{gathered}$ | (1) |
 |
| 图 1 D2Q9模型离散速度方向示意 |
| 图选项 |
格子密度ρ和局部速度u分别为概率密度分布函数的零阶矩和一阶矩：
| $\rho=\sum\limits_i f_i, \rho \boldsymbol{u}=\sum\limits_i f_i \boldsymbol{e}_i.$ | (2) |
| $\begin{aligned}\boldsymbol{F}_{\mathrm{int}} & =-\beta \psi(\boldsymbol{x}) \sum\limits_i G\left(\boldsymbol{x}+\boldsymbol{e}_i\right) \psi\left(\boldsymbol{x}+\boldsymbol{e}_i\right) \boldsymbol{e}_i- \\& \frac{1-\beta}{2} \sum\limits_i G\left(\boldsymbol{x}+\boldsymbol{e}_i\right) \psi^2\left(\boldsymbol{x}+\boldsymbol{e}_i\right) \boldsymbol{e}_i.\end{aligned}$ | (3) |
| $\psi(\boldsymbol{x})=\sqrt{\frac{2\left(\boldsymbol{p}_{\mathrm{EOS}}-\rho \boldsymbol{c}_{\mathrm{s}}^2\right)}{6 g_0}} .$ | (4) |
| $\begin{gathered}p=\frac{\rho R T}{1-b \rho}-\frac{a \rho^2 \alpha(T)}{1+2 b \rho-b^2 \rho^2}, \\a=0.45724 R^2 T_{\mathrm{c}}^2 / P_{\mathrm{c}}, \\b=0.0778 R T_{\mathrm{c}} / P_{\mathrm{c}} . \\\alpha(T)=\left[1+\left(0.37464+1.54226 \omega_{\mathrm{p}}-\right.\right. \\\left.\left.0.26992 \omega_{\mathrm{p}}^2\right)\left(1-\sqrt{T / T_c}\right)\right]^2 .\end{gathered}$ | (5) |
重力的计算公式为
| $\boldsymbol{F}_{\mathrm{g}}(\boldsymbol{x})=\left(\rho(\boldsymbol{x})-\rho_{\mathrm{v}}\right) \boldsymbol{g}$ | (6) |
上述伪势力、重力和3.1节中涉及的壁面黏附力等作用力均需通过合适的方式引入上述方程中以恢复正确的流体力学控制方程。本文中采用Shan等[14]提出的速度偏移法引入作用力项，通过速度偏移反映作用力项的影响，偏移后的速度ueq计算公式如下：
| $\boldsymbol{u}^{\mathrm{eq}}=\boldsymbol{u}+\frac{\tau \boldsymbol{F}}{\rho} \delta_t$ | (7) |
1.2 温度场模拟模型传热模型利用温度分布函数来描述整个温度场，离散速度方向与密度场相同，为D2Q9模型，温度分布函数的演化方程利用单松弛时间方法表示如下：
| $\begin{gathered}h_i\left(\boldsymbol{x}+\boldsymbol{e}_i \delta_t, t+\delta_t\right)=h_i(\boldsymbol{x}, t)- \\\frac{1}{\tau_{\mathrm{T}}}\left(h_i(\boldsymbol{x}, t)-h_i^{\mathrm{eq}}(\boldsymbol{x}, t)\right)+\delta_t \omega_i \phi .\end{gathered}$ | (8) |
| $\phi=T\left[1-\frac{1}{\rho C_v}\left(\frac{\partial p}{\partial T}\right)_\rho\right] \nabla \cdot \boldsymbol{U} .$ | (9) |
| $\left(\frac{\partial p}{\partial T}\right)_\rho \approx \frac{\rho R_{\mathrm{P}}}{1-b_{\mathrm{P}} \rho} .$ | (10) |
在整个传热模型中，每个格子处的温度值可利用该点处的温度分布函数相加。
| $T=\sum\limits_i h_i.$ | (11) |
2.1 汽液共存密度曲线问题在汽液共存密度曲线问题的模拟中，计算域设置为200(Lx)×200(Ly)(Lx和Ly分别表示x和y方向表征长度的格子单位，如无特殊说明，下文中各物理参数均为无量纲化的格子单位)，在计算域中心设置半径为50的圆，圆外密度略低于临界密度，圆内密度略高于临界密度，充分演化后，得到稳定后的汽液共存密度，即为饱和密度。将饱和密度与根据Maxwell等面积理论求得的解析解[13]进行对比，如图 2所示为汽液共存曲线对比图，纵坐标为T与临界温度TC的比值，即约化温度Tr，横坐标为汽液共存的两相密度。Tr>0.77时，程序计算结果与理论分析结果基本吻合，说明程序模拟很好地保证了热力学一致性。当温度继续降低时，气相密度较低，汽液密度比增大，很容易发生数值发散，因此下文的模拟均在Tr>0.77温度区间内。
 |
| 图 2 汽液共存曲线对比 |
| 图选项 |
2.2 Nusselt层流降膜凝结问题参考Liu等[8]的研究，为获得网格无关性的模拟结果，计算域设置为150×500，在计算域内充满蒸汽温度Ts=0.96Tc的饱和蒸汽，左边界为恒温壁面边界温度Tw= 0.91Tc的竖直平板，上边界和右边界使用Guo等[17]提出的非平衡外推法实现等压边界，下边界采用Neumann边界，左边界为无滑移壁面。蒸汽遇冷在过冷壁面凝结，受重力作用, 液相向下流动形成稳定的液膜。将饱和汽液密度的均值ρave=(ρl+ρv)/2作为两相分界线，ρl为液相密度，得到图 3a所示的液膜厚度曲线，出口处两相分界线在水平方向上的位置x约为28.9。分析解为Nusselt层流模型的结果，相关参数的计算参考Liu等[8]的研究；修正分析解为采用Ja修正液膜过冷度的分析解[9]。结果表明：数值模拟解与修正分析解基本吻合[16]。图 3b中展示了水平方向上出口约化温度Tr的分布曲线和约化密度ρr=ρ/ρave的分布曲线，x < 28.9的液膜内温度沿水平方向近似呈线性分布，x>28.9的气体域内温度基本保持不变，且在两相分界线处，液膜存在一定的过冷度。由于存在壁面亲水效应，近壁面处液体密度略微增大[8]。
 |
| 图 3 Nusselt层流降膜凝结数值模拟结果 |
| 图选项 |
综合上述验证可知，模型程序适用于冷凝过程的模拟。
3 冷凝相变两相流模拟3.1 黏附力系数与接触角的对应关系不同材料和结构对不同流体的黏附力有差异，称为材料或结构的亲疏水特征。在格子Boltzmann方法中，一般采用给紧邻壁面的流体层施加作用力(即固壁黏附力)的方式实现亲疏水特性的调节，宏观表现为固、液、气三相点处接触角的变化。在本文模拟中用以描述这种表面黏附力Fads的公式为[18]
| $\boldsymbol{F}_{\mathrm{ads}}=-g_{\mathrm{s}} \psi(\boldsymbol{x}) \sum\limits_{i=0}^8 \omega_i s\left(\boldsymbol{x}+\boldsymbol{e}_i \delta t\right) \boldsymbol{e}_i.$ | (12) |
在无重力等其他外力项干扰时，模拟伪势力与壁面黏附力同时作用时，近壁处的单液滴稳定后可得黏附力系数与固壁静态接触角的对应关系，如图 4所示，其比例系数与黏附力的描述形式和伪势的强度系数的选择有关。
 |
| 图 4 黏附力系数与固壁接触角的对应关系 |
| 图选项 |
由于壁面处施加了反弹的无滑移边界，三相点周围格点的伪势在壁面切线方向上对称相等，因此在gs=0时，三相点周围格点的液相与气相必须对称分布才满足受力平衡条件，此时表现为接触角呈90°。当gs≠0时，黏附力的引入破坏了上述的受力平衡，三相点处周围汽液两相重新分布改变伪势力的大小和方向以重新达到受力平衡状态，宏观表现为三相点处气液界面的切线方向发生变化呈现出不同的接触角。
2.2节的液膜模拟中，亲水壁面完全被液相覆盖，由于gs < 0，在黏附力的作用下近壁处流体的密度有突增现象，如图 3b所示，采用引入黏附力的方式改变壁面亲疏水特性必然会导致近壁层流体密度改变。但图 3b表明，该影响波及的范围较小且相对于本文计算问题中密度的变化幅度，此变化量不显著。
3.2 竖直通道内饱和蒸汽凝结竖直通道内冷凝换热是反应堆内自然循环换热过程的重要一环，除经过特殊强化技术处理的超疏水表面外[5]，工程实际中常见的过冷壁面上主要通过膜状凝结导出热量。考虑与重力方向和液膜厚度方向均垂直的方向上没有壁面条件变化(如温度、壁面亲疏水性等)的工况，将模拟问题近似简化为二维竖直通道内的饱和蒸汽凝结以提高计算效率。如图 5所示为计算域示意图，参考2.2节模型验证时的计算域设置，本节模拟中计算域尺寸设置为150×600，初始时计算域内充满静止的饱和蒸汽(Ts=0.96Tc)，当y < 500时，在左右两侧壁面设置低温的无滑移边界，下边界为绝热的Neumann边界，上边界为定压的自由速度流入边界[17]。相较于2.2节计算域，在上边界附近增设了100个网格的绝热区域，是由于在实际模拟中(见3.2.1节所述)发现通道内的蒸汽凝结与2.2节的竖壁上凝结不同，凝结产生的压差会驱动饱和蒸汽从入口流入且流速不可忽略，增设绝热区域可降低饱和蒸汽在通道内凝结时的入口效应[1]，模拟的重点区域尺寸仍然为150×500。
 |
| 图 5 计算域示意图 |
| 图选项 |
3.2.1 冷凝过程中蒸汽质量流速及壁面热流的变化规律本节以接触角θ=51°为例分析亲水壁面上蒸汽凝结过程中的流动特征和传热速率变化。
冷凝过程中流场变化趋势示意图如图 6所示，为便于观察，将密度场进行二值化处理，红色区域为液相，蓝色区域为气相，x、y分别代表水平和竖直方向的位置，t代表格子单位下当前的时间步。近壁处蒸汽能很快液化形成薄液层覆盖过冷壁面，蒸汽冷凝后形成压差驱动蒸汽向壁面迁移。与2.2节中的半无限大平板上蒸汽凝结过程不同，通道内的蒸汽只能从上边界流入，因此液膜率先在竖直方向位置y=500处开始构建，随着凝结过程的进行，饱和蒸汽的流速逐渐增大并向下发展，使下端蒸汽得到补充，凝结强度提升，且上方凝结液受重力作用开始向下流动，下方的液膜逐渐发展。液相受重力作用逐渐向下流动。
 |
| 图 6 冷凝过程中流场变化趋势示意图 |
| 图选项 |
在整个冷凝过程中随着液膜厚度和蒸汽流速的发展，壁面传热速率经历2个波次变化后逐渐达到稳态(见图 7)。初始状态(0~1 000步)下饱和高温蒸汽直接接触低温的过冷壁面，蒸汽凝结释放汽化潜热，瞬态热流很高，但蒸汽质量流速仍然处于较低水平，蒸汽未能及时补充，而后(1 000~2 500步)温度场逐步扩散为稳定的温度梯度，热流密度逐渐减小。忽略前期瞬态过程热流密度的小幅波动，约2 500步后，随着计算域内蒸汽质量流速逐渐增大并向下发展，壁面整体的凝结强度得到强化，传热速率不断增大，整个过冷壁面上的液膜厚度逐渐增加；在4 000步时壁面平均热流密度开始下降，最终14 000步之后壁面平均热流密度趋于稳定，形成稳定的液膜。如图 7所示，在形成稳定液膜之前，热流密度一直保持较高水平，凝结强度较高，驱动蒸汽自发向凝结区域流动，且热流密度越大蒸汽质量流速的增长越快，形成稳定液膜之后，蒸汽的入口平均质量流速也趋于稳定。
 |
| 图 7 蒸汽入口平均质量流速及壁面平均热流密度的变化趋势 |
| 图选项 |
3.2.2 通道宽度对蒸汽凝结的影响为研究通道宽度对凝结过程的影响，对计算域进行一定调整，将Lx分别设置为150、200、250、300、500和700，针对3.2.1节的工况进行模拟。
如图 8所示，竖直窄通道内的液膜厚度远小于2.2节中Nusselt层流的Ja修正分析解，且通道宽度越窄液膜越薄。液膜变薄的原因之一为计算域内蒸汽流动对液膜施加黏滞应力，凝结过程中蒸汽流动方向与重力方向相同使液膜拉薄。蒸汽质量流速如图 9所示，在初始阶段，各算例壁面过冷度和过冷面积相同，蒸汽的凝结强度一致，通道越窄入口的蒸汽流速越快。凝结过程中，窄通道的入口蒸汽质量流速先达到稳态，通道宽度为150时的稳态入口平均质量流速约是通道宽度为500时的80%，稳态壁面平均热流密度为93.5%，之后，随着通道宽度的增加，稳态入口蒸汽质量流速增量减小。稳态时通道宽度与壁面平均热流密度及主流蒸汽出口约化温度的关系见图 10，通道越宽主流蒸汽出口约化温度越高、传热速率也越高。
 |
| 图 8 通道宽度对稳态液膜厚度的影响 |
| 图选项 |
 |
| 图 9 通道宽度对蒸汽质量流速发展趋势的影响 |
| 图选项 |
 |
| 图 10 稳态时通道宽度对主流蒸汽出口约化温度及壁面平均热流密度的影响 |
| 图选项 |
综合以上结果可知：初始凝结时虽然窄通道的入口蒸汽流速较高，但受限于通道宽度，蒸汽质量流量依然较低，计算域内的低压会抑制蒸汽相变释放相变潜热，热传导在换热过程中所占的比例逐渐升高，主流蒸汽在通道内冷却、密度下降，膜状凝结的过冷度降低，液膜薄、传热速率低。
3.2.3 壁面亲疏水特性对蒸汽凝结的影响根据3.1节的分析，壁面亲疏水特性可通过改变式(10)中的gs进行调节。以图 5所示的计算条件为基础，设置平整、均匀、光滑的无滑移壁面，不添加壁面结构或温度波动等非对称条件，这种仅依赖固体表面能改变壁面亲疏水特性的方式，一般在普通亲水/疏水的调节范围。设置gs分别为-0.4、-0.2、0.2、0.4和0.6，对应接触角θ约为51°、72°、105°、117°和127°(见图 4)。
不同接触角条件下液膜形态的变化趋势如图 11所示，在亲水和弱疏水壁面(接触角θ=105°)上，液膜发展与3.2.1节所述基本一致，靠近上边界的过冷壁面上率先形成稳态的液膜，后逐渐向下扩展，最终达到稳态。θ=117°和θ=127°的壁面为疏水壁，相变起始时刻较晚，且凝结初期不会形成薄液层覆盖壁面，而是在靠近上边界的过冷壁上形成珠状凝结，液滴合并后三相线在过冷壁面向外扩展逐渐形成液膜，由于疏水壁面表面黏附力较小，液膜受重力作用向下滑落，疏水程度越高，液膜滑移速度越快。
 |
| 图 11a 和 11b 为亲水壁面，图 11c 、11 d和11e为疏水壁面。 图 11 不同亲疏水表面冷凝液膜变化趋势 |
| 图选项 |
模拟过程中蒸汽凝结形成的两相分界面为曲面，由于表面张力的存在，蒸汽需要一定的过冷度才能维持液滴，过冷度与液滴半径有关[19]。根据施明恒等[20]分析，温度为Ts的蒸汽中液滴核化存在临界半径rc，只有半径大于临界半径的液滴才能进一步长大，因此在二维模拟时，蒸汽中液滴核化的必要条件为：以核化点为圆心，在半径为rc的圆形范围内，蒸汽温度都必须低于液滴核化所需的过冷度。在平直壁面附近[16]，固壁黏附力的引入改变了液滴的力平衡状态，降低了核化所需的表面能，即降低了核化所需的过冷度。其中冷凝过程的减小系数f(θ)为
| $f(\theta)=\frac{1}{4}\left(2-3 \cos \theta+\cos ^3 \theta\right) .$ | (13) |
壁面亲疏水性对液膜特征的影响符合理论分析解的结果[19-20]。亲疏水性对核化起始时刻的影响如图 11e所示，接触角较大时，过冷度较高，在壁面核化临界半径范围内蒸汽温度下降的导热过程时间更长，液滴核化的起始时刻较晚；针对液膜起始点的分析发现，图 11c、11d和11e中，疏水壁面被液膜覆盖时液膜起始点低于过冷壁面的起始点，且随着疏水性的增加，液膜形成的起始点逐渐降低，原因是入口高温饱和蒸汽来流的直接冲击使上方过冷壁面起始点附近温度较高，难以达到液滴核化的过冷度要求；亲水壁面上的液膜发展如图 11a和11b所示，亲水壁面上的液膜发展受壁面亲水程度影响较小，壁面接触角θ为51°和72°时，出口处液膜厚度相等。原因是壁面被液膜覆盖，相变主要发生在汽液交界面上，热量通过液膜传递到壁面，壁面特性对凝结的影响较小。
凝结过程中，不同亲疏水表面平均热流密度的变化如图 12所示。当接触角θ为51°、72°和105°时，壁面上平均热流密度变化趋势与3.2.1节中的描述一致，壁面传热速率经历2个波次变化后逐渐达到稳态。接触角对传热速率的影响主要体现在换热初期，接触角越小，核化的过冷度要求越低，可在更早的时刻有更多区域满足核化条件，传热速率迅速增大；同时根据3.1节中的分析可知，亲水条件下，近壁面处液体的密度轻微增大，也会产生强化传热的效果，所以接触角越小，换热初期热流密度增速更快、极值更高。对于弱疏水壁面，由于部分壁面上未发生相变，通过蒸汽导热的传热速率较小，因此各阶段的平均热流密度相对更低。
 |
| 图 12 不同亲疏水表面平均热流密度变化趋势图 |
| 图选项 |
疏水壁面上热流密度的变化趋势与亲水壁面的差别较大，壁面传热速率随液膜的生长和滑移脱落波动变化。初始状态下饱和高温蒸汽直接接触低温壁面，瞬态热流较高，由于没有相变释放热量，高热流密度不能维持，逐步构建蒸汽导热的温度梯度，热流密度不断下降。达到核化条件后，热阻降低，热流密度开始增大，凝结液累积到一定尺寸后受重力作用沿壁面向下滑移脱落，壁面暴露在蒸汽中通过蒸汽导热，平均热流密度逐渐降低，接触角越大，核化过冷度越高、液膜滑移速度越快，因此在各阶段的热流密度极值更小、下降更快，但能够更快进入液膜再覆盖阶段。以接触角θ=127°为例，当步数为10 000步时，壁面平均热流密度达到极大值，约是接触角为51°时壁面平均热流密度第一个极大值的75.8%，当步数为69 000时，液膜滑移出计算域，平均热流密度达到最低，随后经过10 000步的等待时间，壁面处液滴核化，液膜再覆盖过程传热速率迅速提升。
4 结论为更好地认识竖直通道内凝结换热过程，本文采用伪势模型对二维通道内初始静止的饱和蒸汽进行模拟。分析了通道宽度及壁面亲疏水性对流动和传热特性的影响，主要结论如下：
1) 亲水壁面上饱和蒸汽遇冷迅速形成薄液层覆盖壁面，靠近上边界的过冷壁面上率先形成稳态的液膜，后逐渐向下扩展。
2) 壁面上蒸汽冷凝会自发驱动蒸汽向下流动，保持亲水壁面上壁温和通道宽度恒定的条件下，在液膜发展过程中，当壁面热流较大时，蒸汽质量流速增长较快。通道较窄时，会抑制蒸汽相变释放潜热，蒸汽质量流速初期增速较快但迅速达到稳态，通道宽度为150时的稳态入口平均质量流速约是通道宽度为500时的80.0%，稳态壁面平均热流密度约为93.5%。
3) 将壁面上饱和蒸汽凝结过程分为蒸汽导热和液滴核化2个阶段进行分析，在不同亲疏水特性下，模拟结果与理论分析解的结果一致，即壁面疏水性越强，液滴核化的起始时刻越晚、液膜在竖直方向上的起始点越低。
4) 对不同亲疏水壁面上传热速率的分析表明：亲水壁面上传热速率受壁面亲水程度影响较小，普通疏水壁面上珠状凝结难以维持，被液膜覆盖后与亲水壁面相比，传热速率较慢，液膜滑移出计算域之前，壁面接触角为127°时的壁面平均热流密度最大值约是接触角为51°时的75.8%，并随液膜滑移逐渐降低，但液膜受重力去除后再形成的过程能在一定范围内强化传热速率。
参考文献
| [1] | 李常伟. 非能动余热排出换热器换热特性研究及换热模型验证[D]. 哈尔滨: 哈尔滨工程大学, 2013. LI C W. Study the heat transfer characteristics and verify the heat exchanger model of passive residual heat removal heat exchanger[D]. Harbin: Harbin Engineering University, 2013. (in Chinese) |
| [2] | WANG M J, ZHAO H, ZHANG Y P, et al. Research on the designed emergency passive residual heat removal system during the station blackout scenario for CPR1000[J]. Annals of Nuclear Energy, 2012, 45: 86-93. DOI:10.1016/j.anucene.2012.03.004 |
| [3] | MARCEL C P, FURCI H F, DELMASTRO D F, et al. Phenomenology involved in self-pressurized, natural circulation, low thermo-dynamic quality, nuclear reactors: The thermal-hydraulics of the CAREM-25 reactor[J]. Nuclear Engineering and Design, 2013, 254: 218-227. DOI:10.1016/j.nucengdes.2012.09.005 |
| [4] | BHOWMIK P K, SCHLEGEL J P, REVANKAR S. State-of-the-art and review of condensation heat transfer for small modular reactor passive safety: Experimental studies[J]. International Journal of Heat and Mass Transfer, 2022, 192: 122936. DOI:10.1016/j.ijheatmasstransfer.2022.122936 |
| [5] | 黄潇立, 陈泽亮, 桂南, 等. 石墨烯强化沸腾传热研究进展及应用综述[J]. 清华大学学报(自然科学版), 2022, 62(10): 1681-1690. HUANG X L, CHEN Z L, GUI N, et al. Review of graphene enhanced boiling heat transfer[J]. Journal of Tsinghua University (Science and Technology), 2022, 62(10): 1681-1690. (in Chinese) |
| [6] | 马喜振, 贾海军, 刘洋. 竖直管道内强迫循环下非凝性气体对蒸汽冷凝的影响[J]. 清华大学学报(自然科学版), 2017, 57(5): 530-536. MA X Z, JIA H J, LIU Y. Effect of non-condensable gases on steam condensation in a vertical pipe with forced convection[J]. Journal of Tsinghua University (Science and Technology), 2017, 57(5): 530-536. (in Chinese) |
| [7] | LI Q, LUO K H, KANG Q J, et al. Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[J]. Progress in Energy and Combustion Science, 2016, 52: 62-105. DOI:10.1016/j.pecs.2015.10.001 |
| [8] | LIU X L, CHENG P. Lattice Boltzmann simulation of steady laminar film condensation on a vertical hydrophilic subcooled flat plate[J]. International Journal of Heat and Mass Transfer, 2013, 62: 507-514. DOI:10.1016/j.ijheatmasstransfer.2013.03.002 |
| [9] | MUKHERJEE A, BASU D N, MONDAL P K, et al. Characterization of condensation on nanostructured surfaces and associated thermal hydraulics using a thermal lattice Boltzmann method[J]. Physical Review E, 2022, 105(4): 045308. DOI:10.1103/PhysRevE.105.045308 |
| [10] | LIU L J, LIANG G J, ZHAO H Q, et al. Effect of surface micromorphology and hydrophobicity on condensation efficiency of droplets using the lattice Boltzmann method[J]. Thermal Science, 2022, 26(4B): 3505-3515. |
| [11] | CAI J J, CHEN J T, DENG W, et al. Towards efficient and sustaining condensation via hierarchical meshed surfaces: A 3D LBM study[J]. International Communications in Heat and Mass Transfer, 2022, 132: 105919. DOI:10.1016/j.icheatmasstransfer.2022.105919 |
| [12] | GONG S, CHENG P. Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows[J]. Computers & Fluids, 2012, 53: 93-104. |
| [13] | YUAN P, SCHAEFER L. Equations of state in a lattice Boltzmann model[J]. Physics of Fluids, 2006, 18(4): 042101. DOI:10.1063/1.2187070 |
| [14] | SHAN X W, CHEN H D. Lattice Boltzmann model for simulating flows with multiple phases and components[J]. Physical Review E, 1993, 47(3): 1815-1819. DOI:10.1103/PhysRevE.47.1815 |
| [15] | GONG S, CHENG P. A lattice Boltzmann method for simulation of liquid-vapor phase-change heat transfer[J]. International Journal of Heat and Mass Transfer, 2012, 55(17-18): 4923-4927. DOI:10.1016/j.ijheatmasstransfer.2012.04.037 |
| [16] | VASYLIV Y, LEE D, TOWER T, et al. Modeling condensation on structured surfaces using lattice Boltzmann method[J]. International Journal of Heat and Mass Transfer, 2019, 136: 196-212. DOI:10.1016/j.ijheatmasstransfer.2019.02.090 |
| [17] | GUO Z L, ZHENG C G, SHI B C. Non-equilibrium extra- polation method for velocity and pressure boundary conditions in the lattice Boltzmann method[J]. Chinese Physics, 2002, 11(4): 366-374. DOI:10.1088/1009-1963/11/4/310 |
| [18] | LI Q, LUO K H, KANG Q J, et al. Contact angles in the pseudopotential lattice Boltzmann modeling of wetting[J]. Physical Review E, 2014, 90(5): 053301. DOI:10.1103/PhysRevE.90.053301 |
| [19] | 杨世铭, 陶文铨. 传热学: 3版[M]. 北京: 高等教育出版社, 1998. YANG S M, TAO W Q. Heat transfer: 3rd ed[M]. Beijing: Higher Education Press, 1998. (in Chinese) |
| [20] | 施明恒, 甘永平, 马重芳. 沸腾和凝结[M]. 北京: 高等教育出版社, 1995. SHI M H, GAN Y P, MA C F. Boiling and Condensation Heat Transfer[M]. Beijing: Higher Education Press, 1995. (in Chinese) |
