1.   引　言

沸腾换热是具有相变的气液两相流换热, 也是最有效的散热方式之一, 被广泛应用于日常生活和工业设备中, 如热交换器^[1]、核电站^[2]、火电厂^[3]等. 然而, 其换热性能受到临界热流密度(critical heat flux, CHF)的限制, 当热流密度超过CHF值, 加热器表面温度会急剧增大, 导致沸腾设备烧毁^[4]. 在沸腾传热研究中, 众多学者致力于提高CHF, 从而提高沸腾换热性能.

一般来说, 现有的强化沸腾换热方法可分为被动方法和主动方法两大类^[5]. 被动方法主要通过改变流体性质^[6]、微结构表面^[7]、表面润湿性^[8]等来实现. 以往的沸腾研究通常在平滑表面开展, 近些年, 随着微加工技术的发展, 众多学者针对微结构表面强化沸腾换热展开了大量的研究. Kong等^[9]将方形微针鳍结构与光滑结构结合, 实验研究了基于微针鳍结构的双结构表面池沸腾换热性能. 结果表明, 双结构表面相比于光滑表面和微针鳍表面, 气泡的合并和脱离频率增大, 其CHF明显增加. Kim等^[10]研究了不同直径、高度和间距的圆柱状微结构在提高沸腾换热和CHF方面的有效性, 相比裸露表面, 圆柱微结构表面的传热系数(heat transfer coefficient, HTC)最大提高$200{\text{%}}$, CHF改善$350{\text{%}}$. 除了规则结构, Elkhol等^[11]研究了使用撇除工艺从铜表面制造的凸起钩子形金属特征的沸腾性能, 与光滑铜表面相比, 由于凸起钩子结构增大了换热面积和成核点密度, 并协助液体通过结构底部凹槽, 其HTC最大增强了$96{\text{%}}$, CHF增加了$67{\text{%}}$. 不同于被动方法, 主动强化沸腾换热方法通常利用外部输入功率, 如机械混合、旋转或者施加磁场、电场等手段增加沸腾换热性能. 自从Clubb^[12] 1916年提出电流体动力学(electrohydrodynamics, EHD)强化流体换热以来, 通过外加电场强化沸腾换热的研究受到了广泛关注. Dong等^[13]研究了直流电场对单个气泡的生长、变形和脱离等气泡动力学行为的影响. 结果表明, 受电场作用的气泡形态变为细长球体, 并且在一定的临界电压下, 气泡在脱离壁面时往往会发生断裂. Gao等^[14]用针型电极研究了电场对小型加热器上单个气泡的成核及生长的影响. 他们发现随着电场强度的增加, 气泡脱离直径和脱离频率减小, 而气泡生长时间和等待时间增加. Hristov等^[15]利用棒状电极和网状电极施加电场, 研究了电场对常压下R123的池沸腾换热性能的影响, 发现直流电场增强了核沸腾的HTC和CHF.

为了进一步发挥被动和主动方法的优势, 众多研究人员提出了结构化表面和电场协同强化沸腾换热的理念, 并进行了深入研究. Garivalis等^[16]研究了微重力条件下电场与方柱微结构表面协同强化沸腾传热的问题, 发现相比于地球重力, 微重力下电场对方柱微结构表面池沸腾的CHF提高更加有效. 然而, 微结构表面结合电场并不只产生积极作用, 在某些尺寸结构下电场可能抑制沸腾换热. Quan等^[17]研究了电场作用下光滑表面和肋高为0.2, 0.5, 2.0 mm的直肋表面的池沸腾换热性能, 发现在低过热度区, 电场对光滑表面和三种肋表面的沸腾传热均表现促进作用; 中过热度区, 热边界层随着壁温的升高而变薄, 电场力阻止气泡离开加热表面, 出现传热恶化; 而在高过热度区, 电场增强了光滑表面和肋高为0.2, 0.5 mm的肋表面沸腾换热, 对于2 mm的肋高度, 电场抑制了蒸气的离开, 产生场阱效应, 导致传热恶化. Liu等^[18]在研究电场和微针鳍结构表面的池沸腾耦合效应过程中也观察到了抑制现象. 这说明电场与结构表面的结合对沸腾传热的影响较为复杂.

随着计算机技术的发展, 计算流体力学(computational fluid dynamics, CFD) 在研究相变传热方面得到了广泛运用. 通过数据可视化技术, 数值模拟能够清晰地观测到沸腾过程的气泡动力学行为、温度场、流场以及电场力等, 为沸腾传热的机理研究进一步提供信息. 众多模拟方法中, 相变格子玻尔兹曼(lattice Boltzmann, LB)方法由于无需追踪相界面以及模型易于实现等优点被广泛用于模拟多相流过程. Liu等^[19,20]运用LB方法对相场分布、流场分布和电势分布进行求解, 研究了等温条件下电场中的液滴形变问题, 结果与前人工作符合较好. 近几年, LB方法还被应用于研究电场对沸腾传热的影响. Feng等^[21]提出了耦合电场模型的相变LB模型, 研究电场对平滑表面的沸腾传热影响. 他们发现, 施加电场使气泡发生形变, 将过渡沸腾甚至膜态沸腾转变为核态沸腾, 沸腾在更高过热度下达到CHF, 从而提高CHF. 除了模拟平滑表面上电场对沸腾换热的影响, 最近, Li等^[22] 运用此模型研究了电场对柱状结构表面沸腾传热的影响, 他们也观察到了电场与柱状结构表面结合表现出的场阱效应. 指出在柱状结构顶部的端点位置, 电场强度成倍增加, 使柱根部气泡无法顺利脱离, 最终导致传热恶化. 针对场阱效应, 他们将复合润湿性表面运用到柱顶部分, 以利用电场的作用促进顶部气泡离开.

以上研究表明, 结构表面会影响电场强度的分布以及气泡形态, 进而影响沸腾换热. 显然, 不同的结构下电场对沸腾换热的影响是不同的, 目前关于电场与结构表面协同强化沸腾换热的研究大多使用柱状或者肋状结构. 近年来, 一些学者发现锥翅结构对沸腾换热有一定的强化效果, Lou等^[23]在研究微腔和微翅表面的饱和池沸腾传热工作中发现, 锥翅表面根部相比于平滑表面更容易成核, 成核数量增加, 且CHF大大提高. Ezzatneshan等^[24] 在研究微腔表面的成核过程中, 在锥形腔体两侧增加两个锥翅结构, 发现锥翅结构阻碍了流体的流动, 气泡脱离时间是无锥翅结构情况的1.4倍.

综上所述, 结构化表面和电场作用都可以在一定程度上提高沸腾换热性能, 然而结构化表面耦合电场作用的强化沸腾换热机理的研究尚不充分. 鉴于此, 本文以沸腾换热工业应用中较常见的锥翅结构为例, 采用相变LB模型对电场作用下锥翅表面不同过热度时气泡动力学行为以及沸腾换热性能进行研究, 旨在探索电场与锥翅表面协同强化沸腾换热的机理, 并提高锥翅表面的沸腾换热性能.

2.   数值方法

2.1   LB方法

本文采用文献[25]中的LB模型, 流体的密度分布函数为

其中$i=0, 1, 2, \cdots, Q-1, Q$为离散方向; $f_i({{\boldsymbol{x}}}, t)$是$t$时刻${{\boldsymbol{x}}}$位置的分布函数; ${{\boldsymbol{c}}}_i$是$i$方向的离散速度; $\tau$ 是松弛时间, 它与运动黏度$\nu$有关($\nu=c_{\mathrm{s}}^2(\tau- 0.5){\text{δ}} t$), $\; c_{\mathrm{s}}=c/\sqrt{3}$, $c={\text{δ}} x/ {\text{δ}} t=1$, 其中${\text{δ}} x$是网格步长, ${\text{δ}}t$是时间步长. 在(1)式中, $f_i^{{\mathrm{eq}}}({{\boldsymbol{x}}}, \; t)$是由平衡速度$\; {{\boldsymbol{u}}}$和密度$\; \rho$计算的平衡分布函数:

其中$w_i$是权重系数. 本文采用D2Q9模型, 其离散速度为

其对应的权重系数为$w_0 = 4/9;\; w_i = 1/9$ (i = 1—4); $\; w_i=1/36$ (i = 5—8). 在(1)式中, $\; {{{F}}}_i$为力项

其中$\Delta {{\boldsymbol{u}}}={{\boldsymbol{F}}}{\text{δ}} t/\rho$为${{\boldsymbol{F}}}$控制的速度变化量, ${{\boldsymbol{F}}}$定义为

其中${{\boldsymbol{F}}}_{\rm{int}}$为粒子间相互作用力, ${{\boldsymbol{F}}}_{\rm{ext}}$为外力, 流体所受电场力${{\boldsymbol{F}}}_{\rm{e}}$将在下节具体给出.

其中$\psi({{\boldsymbol{x}}})$为有效质量, 与局部密度和压力有关: $\psi(\rho)=\sqrt{{2(p-\rho c_{\mathrm{s}}^2)}/{G c_{\mathrm{s}}^2}}$. $G$是相互作用强度, 局部压力$p$由状态方程(EOS)得到. 本研究采用Redlich-Kwong (R-K)状态方程:

其中$R$为气体常数, $T$为流体温度, $a = 0.42748\times R^{2} T_{\rm{c}}^{2.5}/ p_{\rm{c}}$, $b = 0.08664 RT_{\rm{c}}/p_{\rm{c}}$. $T_{\rm{c}}$和$p_{\rm{c}}$分别代表临界温度和临界压力. 外力${{\boldsymbol{F}}}_{\rm{ext}}$为

其中${{\boldsymbol{g}}}$为重力加速度, $\rho_{\rm{ave}}$为每一步计算时整个流体区域的平均密度.

能量方程用于求解计算域中的温度分布, 其对应的LB演化方程为^[25]

其中$g_i({{\boldsymbol{x}}}, t)$是温度分布函数. $\tau_{\rm{T}}$为温度的无量纲松弛时间, 与热扩散系数$\alpha=c_{\mathrm{s}}^2(\tau_{\rm{T}}-0.5){\text{δ}} t$有关. 温度的平衡分布函数为

其中${{\boldsymbol{U}}}$是实际的流体速度. (9)式中的源项$\phi$表示为

$c_{\rm{v}}$为流体的比热容. 宏观量表示为

特别指出, 模拟中有关气液界面处的物性${{\boldsymbol{\chi}}}$(如热扩散系数、介电常数、动力黏度、比热容)可表示为

此外, 方程 (11)中的速度梯度算子既可以采用各向同性中心差分算子处理^[26], 也可以采用非平衡态部分局部计算^[27], 由于本工作采用的密度平衡态分布函数中的速度是平衡速度而不是流体的宏观速度, 因此速度梯度算子采用各向同性中心差分^[26]算子处理.

2.2   电场模型

根据电流体动力学理论, 电介质在电场中所受的电场力${{\boldsymbol{F}_{\rm{e}}}}$为^[28]

其中$\rho_{\rm{e}}$是自由电荷密度, $\rho$是流体密度, ${{\boldsymbol{E}}}$是电场强度, $\varepsilon_{0}$为真空介电常数, $\varepsilon$为流体介电常数. 等式(17)右边第一项和第三项分别为库仑力和电致伸缩力, 理想流体中的自由电荷为零, 且流体不可压缩, 因此可以忽略库仑力与电致伸缩力. 第二项为介电泳力, 表示电场中介电常数分布不均匀所产生的力, 其对气泡的的动力学行为有较大影响. 因此流体所受电场力${{\boldsymbol{F}_{\mathrm{e}}}}$为

对于电场强度${{\boldsymbol{E}}}$可由下式求解:

电场强度${{\boldsymbol{E}}}$表示为电势$V$的梯度, 即${{\boldsymbol{E}}}=-\nabla V$, 则(19)式可以改写为

其对应的LB方程如下^[29]:

其中$\tau_{\rm{s}}=3\varepsilon\varepsilon_{0}+0.5$, $\eta_i$和$\eta_i^{{\mathrm{eq}}}$分别表示电势的分布函数和平衡分布函数. 电势的平衡分布函数$\eta_i^{{\mathrm{eq}}}=w_iV$, 其中$V$由(22)式给出:

这里需要指出, (20)式可以看成Poisson方程的特殊形式, 除了上文提到的方法外, 还可以用Chai等^[30,31]发展的适用于复杂管道以及高维问题的方法求解, 有兴趣的读者可以参考文献[30—32].

3.   物理模型和计算参数

图 1为本研究的物理模型, 整个计算域设置为$L_x \times L_y=400\times600$, 红色部分为固体加热台, 位于计算域底部中心区域, 尺寸$L_{\rm{H}} \times H=200 \times 20$, 加热表面亲水, 对应的接触角$\theta=52.5^\circ$, 关于不同润湿性对沸腾过程的影响可以参考文献[8, 33]. 固体以外区域均为流体区域. 在研究锥翅表面时, 固体表面设置5个锥翅结构(居中放置), 相邻锥翅顶点横向距离$L=40$, 为便于处理边界, 将锥翅结构设置为等腰直角三角形, 其尺寸$l \times h=12\times6$. 在计算域顶部和底部连接电压为$V$的导线, 上界面接地以保证上界面电势为0, 加热台上表面电势为$V_0=V$, 如此便会在加热表面与上界面之间形成电场强度为$E_0=(V_0-0)/(L_y-H)$的电场区域(蓝色阴影). 在左右边界, 密度和温度分布函数均使用周期边界条件; 上边界使用对流边界条件^[34], 保证上升气泡离开计算域; 对于下边界, 密度分布函数均采用半反弹格式^[35], 加热台两侧边以及两侧流体区域底部采用绝热边界条件, 加热台底边设置等温加热方式, 另外, 在加热台上表面对温度分布函数采用流固耦合边界^[36]; 对于电场区域, 左右边界设置周期边界条件, 上下边界为非平衡态外推格式^[21].

图 1  物理问题示意图

Fig. 1.  Schematic diagram of the physical problem.

下载: 全尺寸图片 幻灯片

初始时刻, 流体区域充满饱和温度$T_{\rm{sat}}=0.9 T_{\rm{c}}$密度为$\rho_{\rm{l}}$的静态饱和液体. 加热台内部温度为$T_{\rm{sat}}$, 下边界加热温度$T_{\rm{b}} > T_{\rm{sat}}$. 本文LB沸腾模型中, 由于演化方程 (1)中包含非理想气体分子间相互作用力, 该作用力刻画了压力、温度和密度的热力学关系, 当加热器温度高于饱和温度时, 表面自发成核. R-KEOS中的参数设置为$T_{\rm{c}}=0.1961$, $a= 2/49$, $b=2/21$, $R=1.0$, 饱和液体和饱和气体的密度分别为$\rho_{\rm{l}}=5.426,\; \rho_{\rm{v}}=0.8113$, 比热容为$c_{\rm{pl}}= c_{\rm{vl}} = 4.0,\; c_{\rm{pv}} = c_{\rm{vv}} = 2.0$, 动力黏度$\mu_{\rm{l}} = 0.3257,\; \mu_{\rm{v}}= 0.0487$, 热扩散系数$\alpha_{\rm{l}}=0.05,\; \alpha_{\rm{v}}=0.06$, 液体和蒸气的导热比为$\lambda_{\rm{l}}/\lambda_{\rm{v}}$=11.15. 文献[37]在研究加热台物性对沸腾传热的影响时, 取$\; (\rho{c_{\mathrm{p}}})_{\rm{s}}/(\rho{c_{\mathrm{p}}})_{\rm{l}}$为0.4, 1.0和1.5, 发现三者CHF差别不大. 因此, 本文取$\; (\rho{c_{\mathrm{p}}})_{\rm{s}}/(\rho{c_{\mathrm{p}}})_{\rm{l}}$=1.0, 即$\; (\rho{c_{\mathrm{p}}})_{\rm{s}}$=21.704, $\alpha_{\rm{s}}=1.5$, 则固体/液体和固体/气体的导热系数比为$\lambda_{\rm{s}}/\lambda_{\rm{l}}$=30, $\lambda_{\rm{s}}/\lambda_{\rm{v}}$=334.4. 真空介电常数取$\varepsilon_0=2.2360$, 蒸气与液体的介电常数比为$\varepsilon_{\rm{v}}/\varepsilon_{\rm{l}}=0.4472$.

考虑到锥翅结构边界的复杂性, 本文润湿性边界条件选用文献[38]的曲面边界条件实现. 本文涉及物理量用无量纲化表示

其中$h_{\rm{fg}}$为流体的潜热. 无量纲长度、无量纲时间, 无量纲热流分别定义为$L^*=L/l_ 0$, $t^*=t/t_ 0$, $q^*=q/q_ 0$, 相关参数的格子单位与物理单位转换见表 1. 另外绘制沸腾曲线需要时间和空间平均热流

表 1  格子单位与物理单位转换

Table 1.  The unit conversion from lattice unit to physical unit.

符号	格子单位	物理单位	转换因子
$\rho_{\rm{l}}$	5.426	570.02 ${\rm{kg/m^3}}$	106.16 ${\rm{kg/m^3}}$
$\rho_{\rm{v}}$	0.8113	86.13 ${\rm{kg/m^3}}$	106.16 ${\rm{kg/m^3}}$
$l_0$	16	$4.72\times 10^{-6}\;{\rm{m}}$	$2.95\times 10^{-7}\;{\rm{m}}$
$u_0$	0.0358	38.56 ${\rm{m/s}}$	1077.09 ${\rm{m/s}}$
$t_0$	447.8	$1.224\times 10^{-7}\;{\rm{s}}$	$2.734\times 10^{-10}\;{\rm{s}}$
$\nu$	0.06	$0.19\times 10^{-4}\;{\rm{m^2/s}}$	$3.18\times 10^{-4}\;{\rm{m^2/s}}$
$T_{\rm{c}}$	0.1961	647.2 ${\rm{K}}$	3300.36 ${\rm{K}}$
$p_{\rm{c}}$	0.1784	$0.221\times 10^{8}\;{\rm{Pa}}$	$1.24\times 10^{8}\;{\rm{Pa}}$
$c_{\rm{vl}}$	4.0	1405.9 ${\rm{J/(kg\cdot K)}}$	351.48 ${\rm{J/(kg\cdot K)}}$
$h_{\rm{fg}}$	0.624	$0.726\times 10^{6}\;{\rm{J/kg}}$	$1.16\times 10^{6}\;{\rm{J/kg}}$
$\lambda_{\rm{s}}$	32.556	390.67 ${\rm{W/(m\cdot K)}}$	12.0 ${\rm{W/(m\cdot K)}}$
$q_0$	0.01269	$1.69\times 10^{9}\;{\rm{J/(m^2\cdot s)}}$	$1.33\times 10^{11}\;{\rm{J/(m^2\cdot s)}}$
$\varepsilon_0\varepsilon_{\rm{l}}$	2.236	$1.98\times 10^{-11}\;{\rm{F/m}}$	$8.85\times 10^{-12}\;{\rm{F/m}}$
$\varepsilon_0\varepsilon_{\rm{v}}$	1	$8.85\times 10^{-12}\;{\rm{F/m}}$	$8.85\times 10^{-12}\;{\rm{F/m}}$
$V$	1	1096.96 ${\rm{V}}$	1096.96 ${\rm{V}}$

下载: 导出CSV 

| 显示表格

其中$t_a$到$t_b$时间段是沸腾开始后的很长一段时间, $q$是空间平均热流

计算模拟中使用的参数设置为: $h_{\rm{fg}}=0.624$, $\sigma= 0.095$, ${{\boldsymbol{g}}}=(0, -0.00008)$.

4.   模拟结果

4.1   无电场作用时平滑表面与锥翅表面的池沸腾性能和机理

在研究均匀电场作用下的锥翅结构表面池沸腾换热之前, 有必要对无电场作用的平滑表面(即图1无锥翅结构情况下的模型)和锥翅结构表面的池沸腾换热进行研究. 数值模拟中, 锥翅结构尺寸为$h \times l=6 \times 12$.

4.1.1   沸腾曲线

为了直观表示平滑表面和锥翅表面的沸腾换热性能, 绘制了两种表面的沸腾曲线, 如图2所示, 其中壁面过热度$\Delta T=(T_{\rm{b}}-T_{\rm{sat}})/T_{\rm{c}}$, 其中ONB(onset of nucleate boiling)表示核态沸腾起始点. 观察发现: 当$0.06\leqslant \Delta T\leqslant 0.14$, 沸腾处于核态沸腾阶段, 且此阶段锥翅表面的热流密度大于平滑表面; 当$\Delta T > 0.14$, 随着加热温度的增加, 热流密度呈现先减小再缓慢上升的趋势, 沸腾状态分别对应过渡沸腾阶段和膜态沸腾阶段, 这两个阶段锥翅表面的热流密度小于平滑表面, 且过渡沸腾阶段相差较大, 膜态沸腾阶段差距微小; 此外, 锥翅表面的CHF大于平滑表面. 这说明对于不同的沸腾阶段, 平滑表面和锥翅表面的换热性能有所差异. 为了探究不同沸腾阶段两种表面的沸腾换热机理, 以下将从核态沸腾、过渡沸腾和膜态沸腾3个阶段进行分析.

图 2  平滑表面和锥翅表面的沸腾曲线

Fig. 2.  Boiling curves for smooth and conical surfaces.

下载: 全尺寸图片 幻灯片

4.1.2   核态沸腾

本节模拟平滑表面和锥翅表面的核态沸腾过程. 图3表示加热温度$T_{\rm{b}}=0.96 T_{\rm{c}}$, 两种表面在相同时刻的气泡成核状态. 观察图3(a)发现, 在$t^*=33.50$时刻, 在平滑表面两侧生成两个小尺寸气泡, 气泡间相互作用力弱, 沸腾处于初始核态沸腾阶段. 随着成核过程的发展, 在$t^*=37.96$时刻观察到气泡尺寸增大, 并在$t^*=42.43$时刻在体积力(浮力和重力)的作用下脱离壁面. 对于锥翅表面(图3(b)), $t^*=33.50$时刻在翅与翅之间共生成4个气泡, 且气泡尺寸大于平滑表面. $t^*=37.96$时刻, 锥翅表面气泡尺寸增大, 并在$t^*=42.43$时刻出现两侧气泡脱离壁面, 且中间两气泡在表面持续增大的现象. 对比两种表面的气泡形态, 可以发现锥翅表面增加了成核位点, 且气泡尺寸增大.

图 3  $T_{\rm{b}}=0.96 T_{\rm{c}}$时, 平滑表面和锥翅表面的沸腾过程

Fig. 3.  Snapshots of the boiling processes on smooth and conical surfaces under $T_{\rm{b}}=0.96 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

为了进一步研究锥翅表面促进成核的机理, 图4展示了与图3中$t^*=33.50$时刻对应的光滑表面和锥翅表面的温度分布. 可以发现温度在固体内部扩散较为均匀, 而在流-固界面附近有较大的温度梯度, 这是因为固体的导热率远大于流体. 正是因为固体内部温度高, 凸起的锥翅结构内部(图中椭圆标记处)温度远高于其两侧流体, 增大了加热面与流体的换热面积, 从而使锥翅之间饱和液体更易成核. 此外, 成核过程为相变吸热过程, 使沸腾过程的潜热换热增强, 出现气泡所在区域温度较低的现象. 对比发现, 由于光滑表面上方流体主要吸收来自下方平板的热量, 对流换热面积较小, 最终出现成核位点少, 气泡尺寸较小的情况.

图 4  $t^*=33.50$时刻, 加热温度为$T_{\rm{b}}=0.96 T_{\rm{c}}$的温度分布

Fig. 4.  Temperature distribution during the boiling process with $T_{\rm{b}} = 0.96 T_{\rm{c}}$ and $t^*=33.50$.

下载: 全尺寸图片 幻灯片

接下来模拟加热温度为$T_{\rm{b}}=1.04 T_{\rm{c}}$的光滑表面和锥翅表面的核态沸腾过程. 考虑到高过热度下沸腾较为剧烈, 以及加热表面结构的复杂性, 图5除了展示出两种表面的气泡形态演变, 还增加了当前时刻的流场分布. 观察发现, 由于过热度较高, 在两种表面均生成若干大尺寸气泡, 气泡间相互作用力增强, 但并未出现气泡大面积合并现象, 此时沸腾处于充分核态沸腾阶段. 此外, 对比两种加热表面, 平滑表面气泡呈现出随机分布的状态, 并在加热表面随流场向中间区域滑动; 而锥翅表面气泡成核点位于锥翅之间的空隙, 且由于结构表面的存在, 气泡不易在加热表面滑动, 相比于平滑表面, 流场分布较为杂乱. 锥翅表面沸腾过程中流场的无序性增大了表面气泡的扰动, 出现局部气泡合并的现象, 如图5(b)所示, $t^*=33.50$时刻由于锥翅结构两侧气泡的合并, 出现锥翅结构被气泡覆盖的现象(椭圆标记位置). 值得注意的是, 在$t^*=29.03, 31.26$时刻锥翅表面有6个气泡生成, 而平滑表面气泡数量在所示3个时刻均为5个, 说明在加热温度$T_{\rm{b}}=1.04 T_{\rm{c}}$条件下, 锥翅结构仍然具有增加成核点的作用. 为了进一步研究$T_{\rm{b}}=1.04 T_{\rm{c}}$加热温度下平滑表面和锥翅表面的池沸腾换热性能, 图6 给出了此温度下两种表面空间平均热流密度随时间的变化, 为保证沸腾过程处于充分发展阶段, 选取$t^*=20.00$时刻之后的数据. 可以发现锥翅表面的空间平均热流密度在时间尺度上整体高于平滑表面, 证明锥翅表面在较高过热度($T_{\rm{b}}=1.04 T_{\rm{c}}$)下相比于平滑表面仍然能够促进沸腾换热.

图 5  $T_{\rm{b}}=1.04 T_{\rm{c}}$时, 平滑表面 (a)和锥翅表面 (b)的沸腾过程及流场分布

Fig. 5.  Snapshots of the boiling process and the flow field distributions on smooth (a) and conical (b) surfaces under $T_{\rm{b}}=1.04 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

图 6  加热温度$T_{\rm{b}}=1.04 T_{\rm{c}}$下, 平滑表面和锥翅表面空间平均热流密度随时间的变化

Fig. 6.  The time histories of space-averaged heat flux on smooth and conical surfaces under $T_{\rm{b}}=1.04 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

4.1.3   过渡沸腾

本节研究光滑表面和锥翅表面的过渡沸腾过程, 加热温度取$T_{\rm{b}}=1.06 T_{\rm{c}}$. 图7同样给出了此工况下两种表面气泡形态演变以及当前时刻流场分布, 可以发现光滑表面和锥翅表面均存在气泡大面积覆盖加热表面的现象, 气泡相互作用力强. 对于平滑表面, 气泡在$t^*=82.63$时刻占据大部分加热表面, 并在$t^*=84.86$时刻发生合并. 由于流场的作用, 平滑表面合并生成的大尺寸气泡两侧的三相点向加热面中心移动, 并在$t^*=87.09$时刻脱离加热表面, 此时由于气泡的滑动和脱离过程增强了加热表面的气液交换, 促进了加热面两侧新气泡的生成. 对于锥翅表面(如图7(b)), 由于其具有增加成核点的作用, $t^*=82.63$时刻加热表面生成较多气泡, 且气泡形态复杂, 在$t^*=84.86$时刻合并. 与平滑表面不同, 图中标注的三相点从$84.86$到$87.09$时间段仅向中心移动约10个格子单位距离, 并且$t^*=89.33$时刻整个加热表面被蒸气膜覆盖, 这说明锥翅结构表面增大了加热面气泡的滑动阻力. 关于锥翅结构阻碍流体流动的现象, Ezzatneshan等^[24]在研究锥形腔与锥形翅结合表面对气泡动力学行为影响时也有相似的结论.

图 7  $T_{\rm{b}}=1.06 T_{\rm{c}}$, 平滑表面 (a)和锥翅表面 (b)的沸腾过程和对应的流场分布

Fig. 7.  Snapshots of the boiling process and the flow field distributions on smooth (a) and conical (b) surfaces under $T_{\rm{b}}=1.06 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

以上现象表明: 在过渡沸腾阶段, 锥翅结构一定程度上阻碍了加热表面上方流体向中间区域流动, 进而影响气泡在加热表面的滑动和脱离. 气泡长时间停留在加热表面阻碍了加热面上方的热量交换, 最终演变为蒸气膜, 导致传热恶化. 值得一提的是, 关于过渡沸腾阶段锥翅表面气泡更易成膜的现象在沸腾曲线(图2)中也有所体现. 观察沸腾曲线发现, 随着过热度增大, 锥翅表面过渡沸腾阶段热流密度下降梯度远大于平滑表面, 并在$\Delta T = 0.18$先于平滑表面达到膜态沸腾.

4.1.4   膜态沸腾

为了研究膜态沸腾阶段锥翅表面和光滑表面的池沸腾换热机理, 取加热温度$T_{\rm{b}}=1.12 T_{\rm{c}}$. 图8 展示了不同时刻两种表面的气泡形态, 加热表面均被稳定的蒸气膜覆盖. 对比图8(a)和图8(b), 明显看出锥翅表面蒸气膜厚度大于平滑表面蒸气膜厚度, 沸腾过程中较厚的蒸气膜会增大施加在加热表面的沸腾换热热阻. 根据4.1.3节得出的结论, 锥翅结构一定程度上会阻碍加热表面流体的流动, 这进一步增大了对流换热热阻. 为直观表示膜态沸腾阶段锥翅表面和光滑表面的沸腾换热性能差异, 图9 展示了$t^*=96.03$时刻两种表面的局部热流密度, 其中$x=100$和$x=300$分别对应加热器左侧和右侧水平坐标. 可以看出锥翅表面的局部热流密度略低于平滑表面, 且由两侧到中间区域差距逐渐增大, 这是流体从两侧向中间区域流动过程中锥翅结构对流体流动的阻碍效应逐渐叠加所导致.

图 8  $T_{\rm{b}}=1.12 T_{\rm{c}}$, 平滑表面 (a)和锥翅表面 (b)的沸腾过程

Fig. 8.  Snapshots of the boiling processes on smooth (a) and conical (b) surfaces under $T_{\rm{b}}=1.12 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

图 9  加热温度为$T_{\rm{b}}=1.12 T_{\rm{c}}$, $t^*=96.03$时刻平滑表面和锥翅表面的局部热流密度

Fig. 9.  Local heat flux on smooth and conical surfaces with $T_{\rm{b}}=1.12 T_{\rm{c}}$ and $t^*=96.03$.

下载: 全尺寸图片 幻灯片

根据前面的分析和讨论, 相比于平滑表面, 不同沸腾阶段锥翅表面对池沸腾换热性能影响的机理不同. 在初步核态沸腾阶段和充分核态沸腾阶段, 锥翅表面通过增加对流换热面积和成核点促进沸腾换热, CHF得到提高. 过渡沸腾阶段, 锥翅结构增加了加热表面上方流体流动阻力, 阻碍气泡在加热表面滑动和脱离, 热流密度陡然降低. 膜态沸腾阶段, 锥翅结构导致较厚的蒸气膜以及较大的对流换热热阻, 抑制沸腾传热.

4.2   均匀电场对锥翅表面池沸腾换热性能的影响及机理

在沸腾传热工业中, 为防止高热流密度下传热恶化, 壁温极具升高的现象出现, 往往不涉及过渡沸腾阶段和膜态沸腾阶段. 为了获得更好的换热性能, 人们期望将过渡沸腾或者膜态沸腾转变为核态沸腾. 根据4.1节对所有沸腾阶段的模拟研究, 发现锥翅表面在核态沸腾阶段换热性能良好, 而过渡沸腾阶段和膜态沸腾阶段换热性能均比平滑表面弱. 为了进一步增强锥翅表面沸腾换热性能, 本节将通过在锥翅表面施加均匀电场, 研究均匀电场与锥翅结构协同强化沸腾换热. 模拟所用到的电场强度$E_0=0.0517, \; 0.0862, \; 0.1207$, 表面润湿性和锥翅结构尺寸与4.1节相同.

4.2.1   均匀电场对核态沸腾的影响

图10表示加热温度$T_{\rm{b}}=0.96 T_{\rm{c}}$时不同电场强度下的气泡形态演变, 首图的时刻代表各工况下的初始成核时间. 从图10可以发现, 随着电场强度的增大, 初始成核时间增加, 比如无电场作用下, 初始成核时间为$t^*=33.50$, 而当$E_0=0.1207$时, 初始成核时间延长至$t^*=147.39$. 说明均匀电场延迟了锥翅表面初始核态沸腾阶段成核, 且电场强度越大, 延迟作用越强. 其次, 观察气泡大小可以发现, 电场强度越大, 气泡尺寸越小. 另外, 随着电场强度的增加, 成核数量并没有改变, 成核位置(锥翅间隙)变化不大. 以上现象表明, 均匀电场会延迟锥翅表面成核, 减小气泡尺寸, 而对成核位置和成核数量影响不大.

图 10  $T_{\rm{b}}=0.96 T_{\rm{c}}$时, 不同电场强度下锥翅表面的沸腾过程　(a) E₀ = 0; (b) E₀ = 0.0517; (c) E₀ = 0.0862; (d) E₀ = 0.1207

Fig. 10.  Snapshots of boiling processes on the conical surface at $T_{\rm{b}}=0.96 T_{\rm{c}}$ under different electric field intensities: (a) E₀ = 0; (b) E₀ = 0.0517; (c) E₀ = 0.0862; (d) E₀ = 0.1207

下载: 全尺寸图片 幻灯片

为研究均匀电场对锥翅表面初始核态沸腾阶段换热性能的影响, 图11给出了不同电场强度下空间平均热流密度随时间的变化. 图中热流密度第一次急剧上升随后降低的现象对应初始成核气泡的增大和脱离过程, 此时沸腾换热潜热增大, 热流密度大幅度提升. 可以发现电场强度越大, 热流密度首次急剧上升所处的时间越晚, 这进一步证明了均匀电场会延后锥翅表面成核. 另外, 随着电场强度增大, 空间平均热流密度在时间尺度上总体降低, 沸腾换热性能降低.

图 11  加热温度$T_{\rm{b}}=0.96 T_{\rm{c}}$, 不同电场强度作用下锥翅表面空间平均热流密度随时间的变化

Fig. 11.  The time histories of space-averaged heat flux on the conical surface with $T_{\rm{b}}=0.96 T_{\rm{c}}$ under different electric field intensities.

下载: 全尺寸图片 幻灯片

接下来研究均匀电场作用下充分核态沸腾阶段的气泡动力学行为变化. 图12表示$T_{\rm{b}}=1.04 T_{\rm{c}}$, 电场强度为$E_0=0, \; 0.0862$的沸腾过程中的气泡动力学行为, 以及当前时刻电场力的分布(图中白色箭头表示电场力的位置和方向, 且箭头越长, 当前位置所受电场力越大). 当电场强度为零时, 如图12(a), 锥翅间隙的气泡有相互合并的现象, 使部分锥翅结构被基底较大的气泡包裹, 造成大面积干斑. 干斑的扩散和蔓延会减小表面润湿面积, 从而导致CHF的产生^[39,40], 此现象在图5(b) 也有出现. 当$E_0=0.0862$时(图12(b)), 电场力主要作用在气液界面处, 气泡的形态被横向挤压变长. 受挤压的气泡之间距离增大, 气泡间相互作用力减弱, 从而阻止气泡合并.

图 12  $T_{\rm{b}}=1.04 T_{\rm{c}}$下, 电场强度$E_0=0$ (a) 和 $E_0=$$0.0862$ (b)时锥翅表面的沸腾过程及局部电场力分布

Fig. 12.  Snapshots of boiling processes and the distribution of localized electric field forces on the conical surface at $T_{\rm{b}}=1.04 T_{\rm{c}}$ with $E_0=0$ (a) and $E_0=0.0862$ (b).

下载: 全尺寸图片 幻灯片

值得注意的是, 图 12(b)中, 相比于锥翅表面上方其他位置, 锥翅结构尖端附近气泡所受电场力较大(此位置箭头更长), 说明锥翅结构本身对气泡所受电场力也有较大影响. 为了说明这一现象, 接下来仅对电场区域的电场强度进行演化求解, 以分析锥翅结构对其影响. 数值模拟中, 电势与图 12(b)相同, 整个过程不涉及相场和温度场的演化, 计算域温度保持$T_{\rm{sat}}=0.9 T_{\mathrm{c}}$, 流体区域为密度$\rho_{\rm{l}}= 5.426$ 的饱和液体. 图 13表示电势演化一段时间后, 锥翅结构周围电场强度的模$|{{\boldsymbol{E}}}|$的分布. 可以发现在锥翅结构顶端附近存在较大的电场强度, 其最大值约为一般性区域的$6$倍. 由此可以明确, 锥翅结构尖端附近存在较大的电场强度, 从而使尖端附近气泡所受电场力较大, 这种电场和锥翅表面耦合产生的尖端效应对气泡动力学行为有较大影响, 如图 12(b)的$t^*=40.20$时刻所示, 两侧气泡在锥翅结构顶端位置出现向中间弯折的现象, 延缓气泡向中间区域运动, 进一步避免了锥翅间隙气泡之间的合并, 从而增强了沸腾换热性能. 图 14为加热温度$T_{\rm{b}}=1.04 T_{\rm{c}}$, 施加电场强度$E_0=0.0862$与无电场作用时的空间平均热流密度随时间的变化曲线. 可以看出, 在$t^*$=20—60时段, 均匀电场作用下的空间平均热流密度更大, 说明施加电场促进了锥翅表面的沸腾传热.

图 13  锥翅结构周围电场强度的模$|\boldsymbol{E}|$的分布(电势差$V=$$50$, $t^*=44.66$)

Fig. 13.  Distribution of electric field strength $|\boldsymbol{E}|$ around the conical structure (potential difference $V=50$, $t^*=$$44.66$).

下载: 全尺寸图片 幻灯片

图 14  加热温度$T_{\rm{b}}=1.04 T_{\rm{c}}$时, 电场强度分别为$E_0=0$和$E_0=0.0862$的空间平均热流密度随时间的变化

Fig. 14.  The time histories of space-averaged heat flux with $T_{\rm{b}}=1.04 T_{\rm{c}}$ under different electric intensities of $E_0=0$ and $E_0=0.0862$.

下载: 全尺寸图片 幻灯片

由以上得出结论, 尽管在初步核态沸腾阶段, 施加均匀电场延迟了锥翅结构表面成核, 减小了气泡尺寸, 从而抑制沸腾换热. 而在充分核态沸腾阶段, 电场和锥翅结构协同增强沸腾换热. 一方面电场减小了气泡的横向尺寸, 一定程度上避免了加热表面干斑的扩散蔓延. 另一方面, 锥翅结构本身使得锥翅尖端附近的电场强度成倍增大, 有利于加热器表面再润湿, 阻止气泡在锥翅表面合并.

4.2.2   均匀电场对过渡沸腾和膜态沸腾的影响

本节模拟均匀电场作用下的锥翅表面过渡沸腾和膜态沸腾过程, 加热温度分别取$T_{\rm{b}}= 1.06 T_{\rm{c}}$和$T_{\rm{b}}=1.12 T_{\rm{c}}$.

图15为电场强度$E_0=0.1207$作用下, 两种加热温度下的气泡形态演变和电场力分布, 为了与无电场作用时进行对比, 选取的时刻分别与图7和图8相同. 加热温度为$T_{\rm{b}}=1.06 T_{\rm{c}}$时, 如图15(a)所示, 电场作用下加热表面气泡位于锥翅结构的间隙, 没有气泡大面积合并的现象, 沸腾状态表现为核态沸腾, 而图2沸腾曲线所示无电场作用时锥翅表面在加热温度$T_{\rm{b}}=1.06 T_{\rm{c}}$的沸腾阶段处于过渡沸腾, 说明均匀电场能够将过渡沸腾转变为核态沸腾. 此外, 气泡在$t^*=82.63$和$t^*=89.33$时刻呈现出带状, 并且在锥翅尖端位置出现较大弯折, 根据4.2.1节分析这是尖端附近电场强度过大导致的. 尖端位置气泡受到的较大电场强度一定程度上限制了高过热度条件下气泡在加热表面的合并, 有利于液体对加热表面润湿. 图15(b)展示了均匀电场对锥翅结构膜态沸腾阶段气泡动力学行为的影响, 从图中可以发现蒸气膜在尖端附近断裂(图中椭圆标记处), 可见在膜态沸腾阶段均匀电场与锥翅结构协同表现出来的尖端效应依然存在. 与图15(a)不同, 由于过热度过高, 加热表面气化过于剧烈, 图15(b)中加热表面依然会出现气泡大面积覆盖的现象. 图16比较了加热温度分别为$T_{\rm{b}}=1.06 T_{\rm{c}}$, $T_{\rm{b}}=1.12 T_{\rm{c}}$的锥翅表面无电场作用和$E_0=0.1207$的空间平均热流密度随时间的变化. 可以发现施加均匀电场均能够大幅度增加热流密度, 这是均匀电场与锥翅结构协同作用的结果. 从定量上, 对于图16(a)施加电场后的核态沸腾阶段, 其空间和时间平均热流密度提高$75.94{\text{%}}$.

图 15  $E_0=0.1207$时不同壁面过热度得到的锥翅表面的沸腾过程及局部电场力分布　(a) T_b = 1.06 T_c; (b) T_b = 1.12 T_c

Fig. 15.  Snapshots of boiling processes and the distribution of localized electric field forces on the conical surface with $E_0=0.1207$ under different wall superheat degrees: (a) T_b = 1.06 T_c; (b) T_b = 1.12 T_c.

下载: 全尺寸图片 幻灯片

图 16  电场强度$E_0=0$和$E_0=0.1207$时(a)$T_{\rm{b}}=1.06 T_{\rm{c}}$ 和(b)$T_{\rm{b}}=1.12 T_{\rm{c}}$的空间热流密度随时间的变化

Fig. 16.  The time histories of space-averaged heat flux with $E_0=0$ and $E_0=0.1207$ under (a)$T_{\rm{b}}=1.06 T_{\rm{c}}$ and (b)$T_{\rm{b}}=1.12 T_{\rm{c}}$.

下载: 全尺寸图片 幻灯片

4.2.3   沸腾曲线

图17展示了锥翅表面不同电场强度下的沸腾曲线, 为了比较沸腾换热性能, 这里增加了无电场作用的情况. 根据前面的分析讨论, 在初步核态沸腾阶段, 施加均匀电场导致锥翅表面成核时间延后, 气泡尺寸减小, 沸腾被抑制; 对于充分核态沸腾阶段, 作用在气液界面的电场力以及均匀电场与锥翅结构协同作用表现出的尖端效应抑制了加热表面干斑的扩散和蔓延, 从而延后了CHF点的形成; 过渡沸腾和膜态沸腾阶段, 尖端效应更加明显, 气泡发生剧烈形变, 增强了表面润湿. 随着电场强度增大, 过渡沸腾甚至转化为核态沸腾, 沸腾在更高过热度达到CHF点, CHF逐渐提高, 前人工作^[21,41,42]也得出了电场强度与CHF属于递增关系的结论. 相比无电场作用的锥翅表面, 电场强度为$E_0=0.0517, \; 0.0862, \; 0.1207$时, 其CHF提高百分比分别为$2.68{\text{%}}$, $10.33{\text{%}}$, $14.72{\text{%}}$. 相比平滑表面, 如电场强度$E_0=0.1207$时, 电场结合锥翅表面的CHF提高了$28.76{\text{%}}$.

图 17  锥翅表面不同电场强度下的沸腾曲线

Fig. 17.  Boiling curves on the conical surface for different electric field intensities.

下载: 全尺寸图片 幻灯片

5.   结　论

采用格子Boltzmann直接模拟数值方法研究了平滑表面和锥翅表面的池沸腾过程, 通过分析两者不同沸腾阶段的气泡动力学行为、流场分布、温度场及热流密度, 得到无电场作用下锥翅表面沸腾传热机理. 接着基于无电场作用时的沸腾传热机理, 采用耦合电场模型的Boltzmann数值方法, 研究了均匀电场作用下锥翅表面的池沸腾换热. 旨在揭示电场和结构表面协同强化沸腾换热及其机理, 所得结论如下.

1)无电场作用下, 对于不同的沸腾阶段, 锥翅表面和平滑表面的沸腾换热性能和机理不同. 在核态沸腾阶段, 锥翅表面相比于平滑表面具有更多的成核点, 热流密度增大进而提高CHF; 过渡沸腾阶段, 由于锥翅结构的存在, 气泡不易在加热表面滑动脱离, 容易在加热表面形成蒸气膜, 导致热流密度陡然降低; 膜态沸腾阶段, 锥翅结构增大了加热面上方流体的流动阻力和蒸气膜厚度, 对流换热热阻增大, 热流密度略微降低.

2)均匀电场作用下, 初始核态沸腾阶段随着电场强度的增加, 锥翅结构表面初始成核时间逐渐变长, 气泡尺寸逐渐减小, 沸腾被轻微抑制. 充分核态沸腾阶段, 由于作用在气液界面的电场力, 以及均匀电场和锥翅结构协同表现出的尖端效应, 均匀电场对气泡动力学行为的影响表现为减小加热表面气泡横向尺寸, 阻止气泡在加热表面合并, 进而抑制干斑的扩散和蔓延, 增强沸腾换热.

3)过渡沸腾阶段和膜态沸腾阶段, 均匀电场结合锥翅结构表现出的尖端效应更加明显, 表现为抑制锥翅间隙气泡合并, 增大了加热表面的润湿面积, 以至于将沸腾状态转变为核态沸腾, 沸腾在更高过热度达到临界热流密度. 相比锥翅表面和平滑表面, 本文电场结合锥翅表面得到的CHF最大提高14.72%和28.76%.