¹上海理工大学能源与动力工程学院，上海

²上海市动力工程多相流动与传热重点实验室，上海

收稿日期：2022年2月15日；录用日期：2022年3月7日；发布日期：2022年3月15日

摘要

采用最近提出的格子Boltzmann方法，在具有规则排列的锥形结构表面模拟了核态沸腾。结果表明，在锥形结构的表面，障碍物角点处的液体先发生相变，相变只会发生在障碍物形成的凹槽内及角点处。与矩形结构表面相比，在锥形结构表面形成的气泡会完全从壁面脱离。另一方面，对比了不同表面结构间距和结构尺寸对气泡生长和传热性能的影响。当结构间距较大时，气泡从壁面脱离较早，壁面的核态沸腾换热性能较好，也具有较高的临界热流密度。当结构尺寸较大时，气泡从壁面脱离较晚，壁面换热性能变差，临界热流密度减小。

关键词

锥形结构，格子Boltzmann方法，核态沸腾

Numerical Simulation of Nucleate Boiling on Conical Surfaces

Yingxue Li^1,2, Haoyuan Wang^1,2*, Qin Lou^1,2

¹School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai

²Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, Shanghai

Received: Feb. 15^th, 2022; accepted: Mar. 7^th, 2022; published: Mar. 15^th, 2022

ABSTRACT

The lattice Boltzmann model is used to simulate nucleate boiling on a conical surface with regular arrangement. The results show that on the conical structure surface, the phase change process occurs firstly at the corners and in the grooves of the obstacle, and does not take place on the obstacle surface. In contrast to the rectangular surface, bubbles which appear on the conical surface completely depart from the wall. On the other hand, the effects of different surface structure space and structure size on bubbles growth and heat transfer performance are compared. For the large structure space, bubbles depart from the wall earlier, the nucleate boiling heat transfer performance is better, and the critical heat flux is higher. For the large structure size, bubbles depart from the wall later, the heat transfer performance becomes worse, and the critical heat flux decreases.

Keywords:Conical Structure, Lattice Boltzmann Method, Nucleate Boiling

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

核态沸腾(Nucleate boiling heat transfer)是一种公认的高效换热方式。近些年，许多研究者采用各种方法提高沸腾换热效率，例如改变流体性质 [1]，改变表面润湿度 [2]，或者对表面进行改造 [3] 等。其中，表面改造方法较为简单，成本较低，因此得到越来越多的关注。

通常来说，表面改造是将原本平整光滑的表面，通过一系列加工方式，在表面形成空腔或者障碍物，这些空腔或者障碍物可以作为沸腾中气泡的成核位点 [4]。Kandlikar [5] 的实验结果表明，和平直表面相比，空腔中的传热系数提高了8倍，临界热流密度(Critical heat flux, CHF)提高了2.5倍。不平整的表面促进了表面液体和蒸汽的流动，更有利于加热表面被液体润湿，因此改造后表面有更高的换热性能。尽管表面改造对核态沸腾换热有明显增强，但是此方面研究尚有许多不足。由于表面结构复杂，难以定量研究表面特征对核态沸腾的影响。因此，许多研究者对表面进行规则化处理。Dong等人 [6] 在具有规则排列的圆形纳米空腔表面进行了沸腾实验，结果表明，减小空腔直径能使成核位点增多，但提高传热效果并不显著。类似地，Ma和Cheng [7] 对比了规则排列的微腔表面和微柱表面，发现微柱表面更容易产生毛细润湿现象，即加热表面更容易被液体润湿，因此具有更高的CHF。Shen等人 [8] 研究了规则排列的矩形柱表面的核态沸腾换热，发现柱高度对换热性能的影响大于柱宽度的影响。通过对已有文献的回顾，可以看出在研究规则排列表面的沸腾时，大多数文献考虑了矩形和圆柱形。然而，其他形状障碍物对核态沸腾换热也有明显影响，例如锥形障碍物。Liter和Kaviany [9] 在沸腾表面添加了锥形结构的障碍物，通过实验和理论分析，发现和平直表面相比，锥形结构障碍物表面液体流动阻力更小，即气泡容易从表面脱离并且壁面更容易被液体润湿，所以锥形结构表面有更好的换热性能和更高的CHF。然而，由于实验的局限性，他们并没有揭示气泡在锥形结构表面的成核，生长和脱离过程，一些锥形结构的参数例如结构间距，结构尺寸对核态沸腾的影响尚未得知。因此，本文对锥形结构表面的沸腾研究做进一步完善。

在实际的沸腾实验中，由于实验条件的限制，气泡的形成、生长和脱离很难被观察到，难以准确地分析浮力、壁面等因素对气泡和壁面传热的影响。因此，许多研究者使用数值模拟的方法来研究核态沸腾。近些年，格子玻尔兹曼方法(lattice Boltzmann method, LBM)被广泛用于沸腾研究 [10] [11] [12]。在LBM中，每个计算格点的状态(气态或液态)由所给的状态方程决定，因此不需要额外追踪相界面。本文选择的相变LBM模型，主要优点包括 [11]：1) 可以模拟整个沸腾过程，包括气泡的成核，2) 边界条件易于实现，3) 模拟中无需假设沸腾等待时间，4) 计算量较小。

在本文中，采用最近提出相变LB模型，在具有锥形结构的加热表面模拟了核态沸腾。具体讨论了包括锥形结构表面气泡的成核、生长及脱离过程，锥形结构间距对气泡生长和CHF的影响以及锥形结构尺寸对气泡生长和CHF的影响。

2. 数值方法与模型验证

2.1. 数值方法

本文采用了Gong和Cheng提出的相变LB模型 [11]，在此模型中，密度分布函数的演化方程为：

$f_i\left(x+e_i{\delta }_t,t+{\delta }_t\right)-f_i\left(x,t\right)=-\frac{1}{\tau }\left[f_i\left(x,t\right)-f_i^{eq}\left(x,t\right)\right]+\Delta f_i\left(x,t\right)$ (1)

其中 ${\delta }_t$ 是格子时间， $f_i\left(x,t\right)$ 是密度分布函数， $e_i$ 是离散速度，i是离散速度方向， $\Delta f_i$ 是外力项， $\tau$ 是松弛时间， $f_i^{eq}\left(x,t\right)$ 是平衡态分布函数，其表达式为：

$f_i^{eq}={\omega }_i\rho \left[1+\frac{e_i\cdot u}{c_s^2}+\frac{{\left(e_i\cdot u\right)}^2}{2c_s^4}-\frac{u\cdot u}{2c_s^2}\right]$ (2)

在公式(2)中，ρ为流体密度，u为平衡态速度，c_s为格子声速，其值为 $c_s=\sqrt{1/3}$，${\omega }_i$ 是权系数，其取值与选取的离散模型有关，本文选取D2Q9 (二维九速)模型，其对应权系数为：当i = 0时， ${\omega }_i=4/9$ ；当i = 1~4时， ${\omega }_i=1/9$ ；当i = 5~8时， ${\omega }_i=1/36$。e_i为离散速度，D2Q9的离散速度模型为：

$e_i=\{\begin{array}{l}0i=0\\ \left(\cos \left[\left(i-1\right)\pi /2\right],\sin \left[\left(i-1\right)\pi /2\right]\right)ci=1\text{-}4\\ \sqrt{2}\left(\cos \left[\left(i-5\right)\pi /2+\pi /4\right],\sin \left[\left(i-5\right)\pi /2+\pi /4\right]\right)ci=5\text{-}8\end{array}$ (3)

其中c为格子速度取c = 1，在公式(1)中，外力项 $\Delta f_i$ 表达式为：

$\Delta f_i\left(x,t\right)=f_i^{eq}\left(\rho \left(x,t\right),u+\Delta u\right)-f_i^{eq}\left(\rho \left(x,t\right),u\right)$ (4)

其中 $\Delta u=F{\delta }_t/\rho$，$F$ 为合力，其由三部分组成：

$F=F_{int}\left(x\right)+F_g\left(x\right)$ (5)

其中 $F_g$ 是重力或者浮力， $F_{int}$ 为液-液相互作用力。 $F_{int}$ 表达式为：

$F_{int}\left(x\right)=-G\psi \left(x\right){\displaystyle \underset{i}{\sum }{\omega }_i\psi }\left(x+e_i{\delta }_t\right)e_i$ (6)

其中G为相互作用力强度， $\psi \left(x\right)$ 为“有效密度”， $\psi \left(x\right)=\sqrt{2\left(p-\rho c_s^2\right)/c_{0}G}$，p为压力。在本文中，选取Redlich-Kwong (R-K)状态方程因为其有良好的精确性和简便性，其表达式为：

$p=\frac{\rho RT}{1-b\rho }-\frac{a{\rho }^2}{\sqrt{T}\left(1+b\rho \right)}$ (7)

重力(浮力) $F_g$ 的表达式为：

$F_g\left(x\right)=\left(\rho \left(x\right)-{\rho }_{ave}\right)g$ (8)

其中， $g$ 是重力加速度， ${\rho }_{ave}$ 是整个计算区域每个计算步的流体平均密度。

另一方面，温度分布函数的演化方程可以由以下公式给出：

$g_i\left(x+e_i{\delta }_t,t+{\delta }_t\right)-g_i\left(x,t\right)=-\frac{1}{{\tau }_T}\left[g_i\left(x,t\right)-g_i^{eq}\left(x,t\right)\right]+{\delta }_t{\omega }_i\varphi$ (9)

其中， $g_i\left(x,t\right)$ 是温度分布函数， ${\tau }_T$ 是温度松弛时间， $g_i^{eq}\left(x,t\right)$ 是温度平衡态分布函数，其表达式为：

$g_i^{eq}={\omega }_{i}T\left[1+\frac{e_i\cdot U}{c_s^2}+\frac{{\left(e_i\cdot U\right)}^2}{2c_s^4}-\frac{U\cdot U}{2c_s^2}\right]$ (10)

其中T为温度，U为流体速度。在公式(10)中，源项 $\varphi$ 与相变有关，其表达式为：

$\varphi =T\left[1-\frac{1}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\right]\nabla \cdot U$ (11)

其中c_v为比热。宏观密度，速度和温度由以下公式给出：

$\rho ={\displaystyle \underset{i}{\sum }{f}_i},\rho u={\displaystyle \underset{i}{\sum }{e}_i}{f}_i,\rho U={\displaystyle \underset{i}{\sum }{e}_i}{f}_i+0.5{\delta }_{t}F,T={\displaystyle \underset{i}{\sum }{g}_i}$ (12)

2.2. 模型验证

根据前人的工作 [13]，在水平光滑表面，气泡从表面脱离时，其脱离时的直径与重力大小有关，并满足 $D\propto {\left|g\right|}^{-0.5}$。接下来，将利用上述相变LB模型，模拟水平光滑表面单个气泡的生长脱离过程，对该经验关系进行验证。

初始时，整个计算区域充满饱和液体，其密度和温度分别为ρ_l和T_sat。计算区域底部为无滑移壁面，在底部中心有一个长度为5格子的热源，其温度恒定为T_w，底部其他区域为无滑移绝热边界。上表面是等温边界并且采用对流边界使流体可以自由流出 [14]。左右两侧为周期边界。在本文中，无滑移边界和等温(绝热)边界均采用半反弹处理格式 [15] [16]。计算区域大小为150 × 450，其他参数设置为：T_sat = 0.9T_c，a = 2/49，b = 2/21，R = 1.0，τ = 1.0，τ_T = 1.0，c_v = 6.0，ρ_l = 5.426，T_w = 0.98T_c。

气泡脱离直径(气泡从壁面脱离时的当量直径)和重力关系如图1所示。图中各点代表了利用上述LB模型在不同重力下得到的脱离直径，红色曲线为拟合曲线。脱离直径与重力满足 $D\propto {\left|g\right|}^{-0.4884}$，与实验关系吻合良好 [13]。

图1. 气泡脱离直径和重力关系

3. 初始条件及边界条件

如图2所示，计算区域大小为L_x × L_y = 490 × 1000，其中L_x为横向长度，L_y为纵向长度。初始时，计算区域内充满静止的饱和液体，温度为T_sat，密度为ρ_l。底部有一个表面被加工为锥形结构的加热器，表面锥形结构等间距排列，加热器底部采用恒温加热，温度为T_w，两端为绝热，加热器上表面及两侧均为无滑移边界。加热器总长为L，绝热边高为H，表面的锥形结构数量为N，宽为D，高为0.5 D，锥形结构的间距为P。计算区域底部无加热器部分为绝热边界，顶部为自由出口采用对流边界条件 [14]，左右两边为周期边界。本文中，无滑移边界和等温(绝热)边界均采用半反弹处理格式 [15] [16]。若无其他声明，本文剩余部分参数设置均为：T_sat = 0.9T_c，ρ_l = 5.426，a = 2/49，b = 2/21，R = 1.0，τ = 1.0，τ_T = 1.0，c_v = 6.0，L = 390，H = 10。

图2. 模型示意图

4. 结果与讨论

4.1. 锥形结构表面气泡的生长脱离

图3给出了气泡在锥形表面的形成，生长及脱离过程，其他参数为：DT_w = 0.095T_c，D = 30，P = 54，其中DT_w = T_w − T_sat。如图所示，在t = 6000时刻，气泡成核出现在锥形结构的两侧角点处，因为固体的热扩散率大于液体的热扩散率，这导致角点处侧壁的温度大于凹槽内表面液体的温度。因此，角点处的液体先发生了相变。t = 13,000时刻，可以看到最两侧的气泡先发生了脱离而且有向中间流动的趋势，因为加热器两侧的液体向加热器中心流动，从而对气泡产生向中心拖拽力 [17]。另一方面，此时凹槽内未脱离的气泡被浮力拉伸，凹槽内没有发生气泡合并。在t = 17,000时刻，中心的气泡从壁面脱离并且发生合并，两侧的气泡直接脱离没有发生合并。这与常见的矩形结构表面不同 [18]，矩形结构的凹槽中气泡发生合并，合并后的气泡不会完全从壁面脱离，会堆积在凹槽中。和矩形结构相比，可以看出锥形结构表面的气泡更容易脱离，因此更有利与液体对加热表面的润湿，锥形结构表面具有更佳的传热效果，这和实验得到的结果相同 [19] [20]。

图3. 锥形结构表面气泡生长脱离过程

为了进一步分析锥形结构表面气泡生长过程，图4给出了气泡在角点处生长时的受力情况。如图所示，水平加热壁面和倾斜加热壁面分别对气泡产生蒸发动量力，F₁和F₂，方向均为平行于壁面且向气泡外部。当气泡生长时，水平壁面温度高于倾斜壁面温度，即F₁ > F₂，所以，F₂在水平方向上的分力小于F₁，导致气泡在水平方向上受力不平衡，从而气泡在壁面发生移动。在矩形结构表面，F₂在水平方向上的分力为0，气泡在F₁的作用下更容易向凹槽内生长，因此气泡在矩形结构表面的凹槽内更容易发生合并 [18]。

图4. 气泡蒸发动量力示意图

4.2. 锥形结构间距的影响

图5给出了锥形结构间距P = 110时的气泡生长脱离过程，其他参数为：T_w = 0.095 T_c，D = 30。在t = 4100时刻，障碍物角点和凹槽内液体都发生了相变，并且角点处的气泡较大，这是因为角点处温度较高。需要注意的是，由于障碍物间距较大，凹槽内的壁面能加热较多的液体，因此凹槽内也产生了相变。凹槽内的蒸汽将两边角点处的气泡连通，形成一个大气泡，如t = 6000时刻所示。随着气泡的生长，t = 9000时，气泡从壁面脱离。与此同时，在脱离的气泡两侧产生了新的气泡。与P = 54的情况相比，气泡脱离时刻发生了提前，因为表面障碍物数量减少，促进了表面液体的流动。可以看出，气泡总是在凹槽内(包括角点)内产生，锥形障碍物表面不会发生气泡成核，因此，凹槽内的换热效率大于障碍物表面。

图5. 锥形障碍物间距110时的气泡生长脱离过程

为了进一步研究锥形结构表面的沸腾现象，下面我们得到了锥形结构表面的沸腾曲线。其中为时间–空间平均热流密度 ${\bar{q}}_t$ 的表达式为：

${\bar{q}}_t=\frac{{\int }_{t_1}^{t_2}{\bar{q}}_w(t)dt}{t_2-t_1}$ (13)

其中t₁和t₂表示沸腾开始后的很长一段时间， ${\bar{q}}_w$ 为空间平均热流密度，其表达式为：

${\bar{q}}_w=\frac{{\displaystyle {\int }_L{\left(-\lambda \frac{\partial T\left(x\right)}{\partial y}\right)}_{s-l}\text{d}x}}{L}$ (14)

其中角标s-l表示固液界面，l为固体的导热率。

图6给出了不同锥形结构间距的沸腾曲线，图中曲线的峰值为临界热流密度CHF，CHF之前为核态沸腾阶段。如图所示，在核态沸腾阶段，壁面的传热性能随着间距的增加而增加。从前文分析可以得知，障碍物间距增加，表面凹槽的总长度增加，有效换热区域增加，因此障碍物间距较长的表面有较好的换热性能。图中每条曲线的峰值即对应的CHF，可以看出，CHF随着间距的增加而增加，这是因为障碍物间距较大时，气泡更容易从表面脱离，液体更容易润湿表面，因此障碍物间距较大的表面有着较高的CHF。

图6. 不同间距下的沸腾曲线

4.3. 锥形结构尺寸的影响

图7给出了锥形结构尺寸为D = 20时的气泡生长过程，其他参数为：T_w = 0.095 T_c，P = 62。t = 4100时，障碍物角点和凹槽内液体都发生了相变，随着沸腾过程的进行，角点处的气泡与凹槽内的气泡发生相互作用，合并成一个大气泡，如t = 6000时刻所示。最终，在t = 8000时刻气泡发生了脱离。与D = 30的情况相比，气泡脱离发生了提前，因为障碍物尺寸减小，液体更容易润湿表面，加快了气泡的脱离。

图8给出了不同锥形结构尺寸下的沸腾曲线，图中曲线的峰值为临界热流密度CHF，CHF之前为核态沸腾阶段。如图所示，在核态沸腾阶段，壁面的传热性能随着尺寸的增加而降低。从前文分析可以得知，障碍物尺寸增加，导致凹槽间距减小，相应的表面凹槽的总长度减小，有效换热区域减小，因此障碍物尺寸较大的表面换热性能较差。图中每条曲线的峰值即对应结构尺寸的CHF，可以看出，CHF随着结构尺寸的增加而减小，这是因为障碍物尺寸较小时，凹槽间距变大，气泡更容易从表面脱离，液体更容易润湿表面，因此障碍物尺寸较小的表面有着较高的CHF。

图7. 锥形障碍物尺寸20时的气泡生长脱离过程

图8. 不同尺寸下的沸腾曲线

5. 结果与讨论

本文采用相变LBM模型，在锥形结构表面模拟了气泡的成核、生长和脱离，并且得到了锥形结构表面的沸腾曲线，具体结论如下：

1) 在锥形结构表面，气泡成核发生在表面障碍物角点处和凹槽内。与矩形结构表面相比，锥形结构表面气泡脱离更容易，液体更容易润湿表面。

2) 随着锥形结构间距的增加，气泡脱离时间提前，核态沸腾阶段的传热效率增加，临界热流密度也增加。

3) 随着锥形结构尺寸的增加，气泡脱离时间变长，核态沸腾阶段的传热效率变差，临界热流密度减小。

基金项目

本研究由国家自然科学基金(51976128)及上海市自然科学委员会项目(19ZR1435700)资助。

文章引用

李迎雪,王浩原,娄 钦. 锥形结构表面核态沸腾数值模拟研究
Numerical Simulation of Nucleate Boiling on Conical Surfaces[J]. 建模与仿真, 2022, 11(02): 373-381. https://doi.org/10.12677/MOS.2022.112034

参考文献