1.   引　言

沸腾是一种典型的气液相变现象, 被广泛应用于微纳机电强化换热、高热流密度电子设备冷却、蒸气发电、医学消毒和药物制取、化工杂质去除、食品加工等领域^[1-3]. 如何增强沸腾传热性能以及提高临界热流密度(critical heat flux, CHF)是沸腾研究的核心目的^[4-6]. 从1961年Clubb^[7]提出施加外部电场可以增强流体的传热后, 研究者们展开了大量的研究工作以了解电流体动力学(electrohydrodynamic, EHD)增强沸腾传热的机理. 20世纪80年代以来, EHD强化沸腾传热的实验研究受到了广泛关注.

在EHD强化沸腾实验中发现电场的放电源和加热表面特性的不同对沸腾换热具有很大影响. Madadnia和Koosha^[8]利用电极丝施加电场, 在直径为1 mm的加热丝上研究了电场作用下孤立气泡的离核直径、成核速率(频率)和成核点密度等特性. 发现当典型热流密度为35.4 kW/m²并保持恒定时, 电极电压为6.0, 6.6和7.5 kV时气泡的尺寸和成核频率减小, 而电极电压在8 kV及以上时未观察到成核现象. 在5 mm × 5 mm矩形加热器表面, Gao等^[9]研究了沸腾过程中EHD对R113制冷剂成核及其后续生长动态的影响. 他们观察到随着外加电场强度的增大, 气泡离开直径(体积)和离开频率减小, 而气泡生长时间和等待时间延长. 他们还观察到对于单个R113气泡的情况, 电场的存在可增强沸腾传热. 实验中除了研究EHD对沸腾过程的影响, 还有学者通过在电场中放置单个气泡, 在等温条件下研究EHD对气泡动力学行为的影响. Dong等^[10]调查了直流电场对单个惰性气泡生长、变形和脱离等行为的影响. 发现在电场作用下, 气泡被电应力的水平分量强烈压缩, 并被垂直分量拉长. Zu和Yan^[11]也在EHD对气泡形变影响的实验研究中观察到相同的现象.

除了电场对气泡成核过程的影响, 一些学者还研究了电场对整个沸腾曲线以及CHF的影响. Quan等^[12]利用网状电极在直径约为18 mm的圆形加热表面对光滑表面和肋表面沸腾过程的电流体动力学进行了实验研究, 并探讨了不同壁面过热度情况下EHD对沸腾换热性能的影响机制. 他们发现电对流是低过热区增强换热的主要机制, 且电场抑制了气泡成核. 在中过热区, 电对流效应减弱, 电场力阻止气泡从受热表面生长和脱离, 导致传热恶化和壁温升高. 在高过热区, 电场力将气柱解体成更小的气泡, 并增大气泡脱离的频率, 有助于提高沸腾传热性能. Liu等^[13]研究了10 mm × 10 mm的矩形表面上均匀电场对不同尺寸微针鳍结构表面池沸腾换热的影响, 发现蒸气层的破裂有利于CHF的增强, 而场陷阱效应不利于CHF的增强. 因此, 电场对CHF的影响由以上两个因素共同决定. 在1 cm × 1 cm的加热表面, Garivalis和Manfredini^[14]研究了微重力条件下电场对微结构表面(方柱)和平面表面池沸腾的影响. 结果表明, 电场和微结构的结合使微重力条件下的CHF比普通表面提高了114%—144%.

以上在一定的电场强度和加热器尺寸下研究了沸腾换热性能, 极大地提高了对该问题的认识. 近年来, 随着计算机技术和数值方法^[15-19]的快速发展, 采用模拟方法研究EHD及其增强沸腾换热受到大量学者的广泛关注. 在6 mm × 6 mm硅橡胶绝缘疏水表面, Songoro^[20]采用有限元法对电场作用下单个水滴的运动进行数值模拟. 结果表明, 在水平和垂直施加交流电场的情况下, 液滴振动的预测共振频率与实验共振频率吻合. 采用耦合水平集(level set)和流体体积(volume-of-fluid)的CLSVOF(coupled level set and volume-of-fluid)方法^[21], Pandey等^[22]研究了在6.825 mm加热器表面上电场和过热度对气泡成核方式与传热特性的影响, 发现在低过热度下, 随着电场强度的增大, 气泡释放速率增大, 同时气泡间距减小, 形成了离散气泡. 当过热度较大时, 较高的电场强度导致形成大量的蒸气柱. 然而, 以上方法无法模拟气泡成核, 也无法对整个沸腾曲线进行刻画. 最近, 格子Boltzmann (LB)方法^[23]由于可以从微观尺度描述分子之间的相互作用、不需要追踪相界面、自动成核、可以方便描述沸腾的所有阶段等优点被广泛应用于气-液相变的问题研究中^[4,5,24-26]. 采用格子Boltzmann方法并耦合电场作用力模型, 张浏斌等^[27]模拟了均匀电场作用下在5个格子单位的小型加热器表面上, 单气泡成核、生长、脱离过程, 重点研究了重力加速度和电场强度对气泡动力学的影响, 发现了均匀电场有助于低重力加速度下池沸腾气泡的脱离, 从而提高换热效率. Yao等^[28]研究了14.4 mm加热器表面上电场强度对气泡行为的影响, 发现了电场可以将气泡变形为扁长形和扁圆形, 分别导致气泡加速和减速离开热表面. 最近, Li等^[29]用格子Boltzmann方法研究了电场作用下在柱状结构表面上的沸腾传热性能, 发现外加电场对柱状结构表面的池沸腾有正、负两方面的影响, 并将润湿性修饰区域应用于柱状结构表面的柱顶, 进一步增强了电场作用下柱状结构表面的沸腾换热.

以上研究表明, 电场对沸腾过程的气泡动力学行为和传热性能的影响较为复杂, 其与电极施加方式、加热表面设置、过热度等有密切联系. 然而, 大部分关于电场增强池沸腾换热的研究都是在特定的加热器长度基础上进行的, 很少有研究人员考虑电场对不同长度加热器池沸腾传热的影响. 相关研究表明: 加热器长度对气泡动力学、流场分布和沸腾传热性能有显著影响. Rainey和You^[30]和Lee等^[31]在实验研究中发现, 随着长度增大, 加热器侧面冷流体与加热表面上方热流体的对流阻力增加, 从而导致CHF逐渐降低. Zhang等^[32]采用二维格子Boltzmann方法模拟不同长度加热器的池沸腾过程也得出类似结论. 施加电场时, 加热器长度对气泡动力学行为以及沸腾传热性能的影响尚不清晰. 因此, 本文将研究不同长度加热器下均匀电场对气泡动力学行为、沸腾曲线、以及临界热流密度的影响. 由于实验工作中不同尺寸加热器样本与电极设置较为复杂, 且基于放电源与电场之间的复杂关系, 很难保证施加于不同加热器样本的电场保持一致. 而LB方法在研究气相相变问题时具有不需要追踪界面、成核自发发生、流固耦合易于处理等优点^[33]. 本文采用二维伪势LB和电场的耦合模型, 数值模拟均匀电场作用下不同长度加热器的池沸腾过程, 详细研究不同长度加热器下均匀电场对沸腾曲线, 气泡动力学行为和CHF的影响.

2.   数值方法

2.1   格子Boltzmann方法

在LB模型^[34]中, 流体的密度分布函数为

$$\begin{split}& f_i({\boldsymbol{x}}+\delta t{\boldsymbol{c}}_i,t+\delta t)-f_i({\boldsymbol{x}},t) \\ =\;& -\frac{1}{\tau}[f_i({\boldsymbol{x}},t)-f_i^{\rm{eq}}({\boldsymbol{x}},t)]+{\boldsymbol{F}}_i, \end{split}$$

(1)

其中$i = 0, 1, 2, \cdots, Q-1, Q$为离散速度数, $f_i({\boldsymbol{x}}, t)$是t时刻x位置的流体密度分布函数, i是分布函数的方向, ${\boldsymbol{c}}_i$是i方向上的离散速度, τ是松弛时间, 它与运动黏度ν的关系为$\nu = c_{\rm s}^2(\tau-0.5)\delta t$, $c_{\rm s} = c/\sqrt{3}$是模型常数, $c = \delta x/ \delta t = 1$, 其中$\delta x$是网格间距步长, $\delta t$是网格时间步长. 在(1)式中, $f_i^{\rm{eq}}({\boldsymbol{x}}, \; t)$是由平衡速度u和密度ρ决定的平衡态分布函数:

$$f_i^{\rm eq} = w_i\rho\left[1+\frac{{\boldsymbol{c}}_i\cdot{\boldsymbol{u}}}{c_{\rm s}^2}+\frac{({\boldsymbol{c}}_i\cdot{\boldsymbol{u}})^2}{2c_{\rm s}^4}-\frac{{\boldsymbol{u}}^2}{2c_{\rm s}^2}\right],$$

(2)

式中$w_i$是与离散速度模型相关的权重系数. 在(1)式中, 力项${\boldsymbol{F}}_i$表达式如下所示:

$${\boldsymbol{F}}_i = f_i^{\rm eq}(\rho({\boldsymbol{x}},t),{{\boldsymbol{u}}+{\boldsymbol{\Delta}} {\boldsymbol{u}}})-f_i^{\rm eq}(\rho({\boldsymbol{x}},t),{\boldsymbol{u}}),$$

(3)

其中$\Delta {\boldsymbol{u}} = {\boldsymbol{F}}\delta t/\rho$, F为相互作用力, 其包含粒子间相互作用力${\boldsymbol{F}}_\mathrm{int}$, 外力${\boldsymbol{F}}_\mathrm{ext}$以及电场力${\boldsymbol{F}}_\mathrm{e}$, 具体表达式如下:

$${\boldsymbol{F}} = {\boldsymbol{F}}_\mathrm{int}+{\boldsymbol{F}}_\mathrm{ext}+{\boldsymbol{F}}_\mathrm{e},$$

(4)

式中粒子间相互作用力${\boldsymbol{F}}_\mathrm{int}$可以写为

$${\boldsymbol{F}}_\mathrm{int}({\boldsymbol{x}}) = -G\psi({\boldsymbol{x}})\sum\nolimits_iw_i\psi({\boldsymbol{x}}+{\boldsymbol{c}}_i\delta_t){\boldsymbol{c}}_i,$$

(5)

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

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

(6)

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

$${\boldsymbol{F}}_{\rm ext}({\boldsymbol{x}}) = {\boldsymbol{g}}(\rho({\boldsymbol{x}})-\rho_\mathrm{ave}),$$

(7)

其中g为重力加速度, $\rho_\mathrm{ave}$为每一步计算时整个计算域的平均密度. 方程(4)中电场力${\boldsymbol{F}}_\mathrm{e}$将在下节具体给出.

温度分布函数$g_i({\boldsymbol{x}}, t)$其对应的LB演化方程为^[34]

$$\begin{split} &g_i({\boldsymbol{x}}+\delta t{\boldsymbol{c}}_i,t+\delta t)-g_i({\boldsymbol{x}},t) \\ =\;& -\frac{1}{\tau_\mathrm{T}}[g_i({\boldsymbol{x}},t)-g_i^{\rm eq}({\boldsymbol{x}},t)]+\delta t w_i\phi, \end{split}$$

(8)

其中$\tau_\mathrm{T}$为温度的无量纲松弛时间, 其由热扩散系数α决定, $\alpha = c_{\rm s}^2(\tau_\mathrm{T}-0.5)\delta t$, 方程(8)中的$g_i^{\rm eq}$为温度平衡分布函数, 表达式如下所示:

$$g_i^{\rm eq} = w_iT\left[1+\frac{{\boldsymbol{c}}_i\cdot{\boldsymbol{U}}}{c_{\rm s}^2}+\frac{({\boldsymbol{c}}_i\cdot{\boldsymbol{U}})^2}{2c_{\rm s}^4}-\frac{{\boldsymbol{U}}^2}{2c_{\rm s}^2}\right],$$

(9)

其中U是实际的流体速度. (8)式中的源项ϕ表示为

$$\phi = T\left[1-\frac{1}{\rho c_v}\left(\frac{\partial{p}}{\partial{T}}\right)\right]\nabla\cdot{\boldsymbol{U}},$$

(10)

其中$c_v$是流体的比热容. 流体的宏观量可由密度分布函数和温度分布函数的矩得到:

$$\rho = \sum\nolimits_i f_i,$$

(11)

$$\rho{\boldsymbol{u}} = \sum\nolimits_i {\boldsymbol{c}}_if_i,$$

(12)

$$\rho{\boldsymbol{U}} = \sum\nolimits_i {\boldsymbol{c}}_if_i+\frac{\delta t{\boldsymbol{F}}}{2},$$

(13)

$$T = \sum\nolimits_i g_i.$$

(14)

特别指出, 模拟中气液界面处的物性χ (如热扩散系数、介电常数、动力黏度、比热容)由(15)式根据气液份额给出：

$${\boldsymbol{\chi}} = {\boldsymbol{\chi}}_\mathrm{l} \frac{\rho-\rho_\mathrm{v}}{\rho_\mathrm{l}-\rho_\mathrm{v}}-{\boldsymbol{\chi}}_\mathrm{v} \frac{\rho_\mathrm{l}-\rho}{\rho_\mathrm{l}-\rho_\mathrm{v}}.$$

(15)

二维格子Boltzmann方法的D2Q9模型^[35,36]对模拟沸腾传热过程具有良好的适用性^[5,29]. 因此, 本文采用D2Q9模型, 其离散速度为

$${{\boldsymbol{c}}_i = }\left\{ \begin{aligned} &(0,0)c, &{ i = 0,}\qquad\\ &\left(\cos{\left[ {(i - 1)\frac{\pi }{2}} \right]}, \sin{\left[ {(i - 1)\frac{\pi }{2}} \right]}\right)c, &{i = 1—4,}\;\\ &\sqrt{2}\left(\cos{\left[ {(2i - 1)\frac{\pi }{4}} \right]}, \sin{\left[ {(2i - 1)\frac{\pi }{4}} \right]}\right)c, &{i = 5—8.}\; \end{aligned} \right.$$

(16)

其对应的权重系数为$w_0 = 4/9$; $w_i = 1/9$($i = 1— 4$); $w_i = 1/36$($i = 5—8$).

2.2   电场模型

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

$${{\boldsymbol{F}}_{\rm{e}}} = \rho_{\rm{e}}{\boldsymbol{E}}-\frac{1}{2}{\boldsymbol{E}}^{2}\nabla(\varepsilon\varepsilon_{0})+\frac{1}{2}\varepsilon_{0}\nabla\left[{\boldsymbol{E}}^{2}\left(\frac{\partial{\varepsilon}}{\partial{\rho}}\right)_t\rho\right],$$

(17)

其中$\rho_{\rm e}$是自由电荷密度, ρ是流体密度, E是电场强度, $\varepsilon_{0}$为真空介电常数, ε为流体介电常数. 等式(17)右边第一项为库仑力, 理想流体中的自由电荷为零, 因此可以忽略库仑力. 第二项为介电泳力, 表示电场中介电常数分布不均匀所产生的力. 在两相流中, 介电泳力主要由气相与液相的介电常数差引起, 这决定了气泡在电场中的运动方向. 第三项为电致伸缩力, 对于不可压缩流体, 此项可以忽略. 因此流体所受电场力${{\boldsymbol{F}}_{\rm{e}}}$为

$${{\boldsymbol{F}}_{\rm{e}}} = -\frac{1}{2}{\boldsymbol{E}}^{2}\nabla(\varepsilon\varepsilon_{0}),$$

(18)

对于电场强度E可由下式求解:

$$\nabla\cdot(\varepsilon\varepsilon_{0}{\boldsymbol{E}}) = 0.$$

(19)

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

$$\nabla\cdot(\varepsilon\varepsilon_{0}\nabla V) = 0.$$

(20)

为了求解(20)式, 本文采用文献[26, 38]中的LB方程:

$$\eta_i({\boldsymbol{x}}+\delta t{\boldsymbol{c}}_i,t+\delta t)-\eta_i({\boldsymbol{x}},t) = -\frac{1}{\tau_{\rm s}}[\eta_i({\boldsymbol{x}},t)-\eta_i^{\rm eq}({\boldsymbol{x}},t)],$$

(21)

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

$$V = \sum\limits_i \eta_i({\boldsymbol{x}},t).$$

(22)

3.   物理模型和计算参数

如图1所示, 整个计算区域设置为$L_x \times L_y = 400\times1000$格子单位, 红色部分为加热台, 位于底部中心, 高度$H = 30$格子单位, 长度为$L_H$, 加热器表面具有不同润湿性, 由润湿性边界条件^[39]描述, 计算区域内红色以外部分为流体区域. 另外, 在加热台下表面及计算域上边界施加大小为$V_0$的电压, 并在计算域上界面接地线, 以保证上界面电势为0, 加热台上表面的电势为$V_0$. 由于加热器上表面与上边界平行, 类似于两平行带电平板, 会在蓝色阴影区域形成电场强度为$E = V_0/(L_y-H)$的均匀电场. 在计算区域的左右边界, 对密度分布函数和温度分布函数采用周期边界格式; 在计算区域的上边界, 对密度和温度分布函数采用适用于自由边界的对流边界格式, 以消除边界处可能给内部流场带来的影响; 在下边界对密度分布函数采用半步长反弹格式, 对温度分布函数在加热台上边界采用流固耦合边界条件^[40], 考虑固体导热率以及流体物性差异, 为确保流固边界温度和热通量传递的连续性, 加热台左右两边及流体侧下边界采用绝热边界条件, 在加热台内部下边界设置等温加热方式. 对于理想电介质模型, 上下边界采用非平衡外推边界, 左右边界为周期边界条件^[26]. 初始时刻, 流体温度为饱和温度$T_\mathrm{sat}$, 密度为$T_\mathrm{sat}$对应的饱和密度$\rho_\mathrm{l}$, 加热器下边界温度为$T_\mathrm{b}$. 在数值模拟中, R-K状态方程中的参数设置为$a = 2/49$, $b = 2/21$, $R = 1.0$, 计算得$T_\mathrm{c} = 0.1961$, $p_\mathrm{c} = 0.1784$. 取饱和温度$T_\mathrm{sat} = 0.9 T_\mathrm{c}$, 则饱和液体和饱和气体的密度分别$\rho_\mathrm{l} = 5.426$, $\rho_\mathrm{v} = 0.8113$. 比热容为$c_{p\mathrm{l}} = c_{v\mathrm{l}} = 4.0$, $c_{p\mathrm{v}} = c_{\rm{v}}{{v}} = 2.0$, 动力黏度$\mu_\mathrm{l} = 0.3257$, $\mu_\mathrm{v} = 0.0487$, 热扩散系数$\alpha_\mathrm{l} = 0.05$, $\alpha_\mathrm{v} = 0.06$, 液体和蒸气的导热比为$\lambda_\mathrm{l}/\lambda_\mathrm{v} = 11.15$, 流体的潜热$h_\mathrm{fg} = 0.624$, 表面张力系数$\sigma = 0.095$, 重力加速度${\boldsymbol{g}} = (0, -0.00008)$. 文献[24]在研究加热器物性对CHF影响时, 对流固耦合边界条件加热器物性取$(\rho{c_p})_\mathrm{s}/(\rho{c_p})_\mathrm{l}$ = 0.4,1.0, 1.5, 发现加热器物性对CHF值影响不明显. 因此, 本文取$(\rho{c_p})_\mathrm{s}/(\rho{c_p})_\mathrm{l}$ = 1.0, 即$(\rho{c_p})_\mathrm{s} = 21.704$, $\alpha_\mathrm{s} = 1.5$, 则固体/液体和固体/气体的导热系数比为$\lambda_\mathrm{s}/\lambda_\mathrm{l} = 30$, $\lambda_\mathrm{s}/\lambda_\mathrm{v} = 334.4$. 真空介电常数取$\varepsilon_\mathrm{0} = 2.2360$, 蒸气与液体的介电常数比为$\varepsilon_\mathrm{v}/\varepsilon_\mathrm{l} = 0.4472$. 需要特别指出的是, 以上单位均为格子单位, 其与实际物理单位的转换关系如表1所示.

为了便于与前人工作进行对比, 如无特别说明均采用无量纲单位, 所采用的特征长度、特征速度、特征速度和特征热流密度如下:

$$\begin{split} &l_0 = \sqrt{\frac{\sigma}{g(\rho_\mathrm{l}-\rho_\mathrm{v})}},\; u_0 = \sqrt{gl_0},\;\\ &t_0 = l_0/u_0,\; q_0 = \frac{\nu \rho_\mathrm{l} h_{\mathrm{fg}}}{l_0}. \end{split}$$

(23)

经过无量纲化得到的无量纲长度、无量纲时间、无量纲热流分别定义为$L^* = L/l_ 0$, $t^* = t/t_ 0$, $q^* = q/q_ 0$. 另外绘制沸腾曲线需要的时间和空间平均热流如下:

$$Q = \left(\int_{t_\mathrm{a}}^{t_\mathrm{b}}{q}{\rm d}t\right)/(t_\mathrm{b}-t_\mathrm{a}),$$

(24)

其中$t_\mathrm{a}$到$t_\mathrm{b}$时间段是沸腾开始后的很长一段时间, q是空间平均热流:

$$q = \left[\int_{L_H}\left(-\lambda\frac{\partial{T({\boldsymbol{x}})}}{\partial{y}}\right)_{y = 0}{\rm d}x\right]/L_H.$$

(25)

4.   模拟结果

4.1   无电场作用时气泡行为的变化对不同长度加热器表面沸腾传热性能的影响

在研究均匀电场对不同加热器长度的沸腾传热影响之前, 首先对不同加热器长度下的沸腾传热进行研究. 本节模拟了不同长度加热器的沸腾过程, 考虑到整个计算域尺寸以及物理模型设置, 格子单位长度取$L_H$ = 50, 100, 120, 150, 200, 300, 对应无量纲长度$L_H^*$ = 3.125, 6.25, 7.5, 9.375, 12.5, 18.75. 为了研究润湿性的影响, 本文分别考虑了亲水表面和疏水表面, 气泡与加热表面静态接触角分别设置为$\theta = 45^\circ$, $\theta = 120^\circ$.

图2(a)和图2(b)分别展示了亲水表面和疏水表面不同长度加热器下的沸腾曲线. 从图2(a)可以观察到对于亲水表面, 当$L_H^* = 3.125$时, 随着壁面过热度增大, 空间时间平均热流密度越来越大, 整个沸腾曲线并没有出现过渡沸腾阶段, 即不存在临界热流密度. 为了揭示这一现象出现的原因, 图3给出了$L_H^* = 3.125$时不同加热温度下同一时刻的气泡形态. 从图中可以发现，当壁面过热度较小, 即$T_\mathrm{b} = 0.98 T_\mathrm{c}$时, 加热表面每次仅生成一个气泡. 当壁面温度升高到$T_\mathrm{b} = 1.00 T_\mathrm{c}$时, 尽管壁面过热度的增大导致生成气泡的尺寸增大, 然而由于加热器尺度较小, 加热表面同样只能生成一个气泡, 即加热表面生成的气体始终以孤立的形式存在, 并没有发生由于气泡的合并导致的沸腾换热恶化的现象. 而当壁面温度升高到$T_\mathrm{b} = 1.02 T_\mathrm{c}$, 加热表面形成稳定的蒸气柱, 蒸气柱顶端持续形成气泡并继续上升, 此沸腾状态为膜态沸腾. 由于加热器尺寸较小, 在加热表面形成的小尺寸蒸气膜在浮力和重力的共同作用下整体处于上升状态. 由以上结果可知, 由于加热器尺寸较小, 没有为出现多个气泡提供充足的空间, 因此没有出现因为多气泡合并和脱离速率降低导致的过渡沸腾阶段, 此时随着壁面过热度的增加, 沸腾由核态沸腾阶段直接转变为膜态沸腾阶段, 空间时间平均热流密度在整个沸腾阶段逐渐升高.

从图2(a)还可以发现, 当$L_H^*\geqslant 6.25$时, 沸腾曲线表现为经典形状, 随着壁面过热度的增加沸腾过程存在过渡沸腾阶段和CHF, 且随着加热器长度增大, 核态沸腾阶段和膜态沸腾阶段的时间空间平均热流密度逐渐减小. Rainey等^[30]和Lee等^[31]在实验研究加热器尺寸对沸腾换热的影响时也观察到了类似的现象, 并将引起该现象的原因归结于较冷流体的再润湿阻力随着加热器长度的增大而增大. 然而, 他们并没有给出引起再润湿阻力增大的原因. 为了揭示这一现象, 图4给出了$t^* = 55.83$, $T_\mathrm{b} = 0.98 T_\mathrm{c}$条件下, 加热器长度分别为$L_H^* = 6.25$和$L_H^* =18.75$时的气泡形态和对应的流场分布, 其中加热器两端到白色标记的距离表示平行于加热表面的流线长度. 从图中可以发现当$L_H^* = 6.25$时加热器表面上生成的气泡数量明显少于$L_H^* = 18.75$时的气泡数量, 且当$L_H^* = 18.75$时, 靠近并且平行于加热表面的流线长度明显大于$L_H^* = 6.25$. 平行于加热表面的流线越长, 加热器两侧较冷流体与加热表面较热流体之间的对流换热越困难, 即增大了较冷流体的再湿阻力, 该现象证实了文献[30]中的结论. 进一步观察图2(a)还可以发现, 当$6.25\leqslant L_H^*\leqslant 9.375$时, CHF随着加热器长度的增大逐渐减小; 而当$L_H^* > 9.375$时, CHF随加热器长度增大变化不大, 这与Zhang等^[32]模拟研究的加热器长度对CHF影响的结果一致. Zhang等^[32]指出由于流阻只发生在加热器边缘附近而不是加热器中心附近, 因此当加热器尺寸足够大, 流阻效应不再随加热器长度的增大而增大, 此时, 加热器尺寸的增大对CHF影响不大显著.

当表面疏水时不同加热器长度对沸腾曲线的影响规律(如图2(b)所示)与亲水表面时类似, 在加热器长度较小时, 即$L_H^* = 6.25$时沸腾曲线同样未出现过渡沸腾阶段. 随着加热器长度的增大, 核态沸腾阶段和膜态沸腾阶段热流密度均减小, CHF同样由逐渐降低到维持相对稳定. 然而当表面疏水时, CHF由逐渐降低到维持相对稳定的阈值增加到为$L_H^* = 12.5$.

4.2   均匀电场作用下气泡行为的变化对亲水表面沸腾传热的影响

本节模拟了在均匀电场作用下亲水表面不同长度加热器的池沸腾传热过程. 加热器长度和亲水表面润湿性与4.1节相同, 即$L_H^*$ = 3.125, 6.25, 7.5, 9.375, 12.5, 18.75, $\theta = 45^\circ$. 均匀电场强度E取0.05442, 0.10884, 0.16327, 这里无电场作用时的结果($E = 0$)用于对比分析电场作用对沸腾过程的影响.

图5为不同长度加热器在不同电场强度下的沸腾曲线, 其中ONB (onset of nucleate boiling)表示核态沸腾起始点. 从图5(a)和图5(b)可以发现，当$L_H^*$ = 3.125, 6.25时, 均匀电场对核态沸腾阶段和膜态沸腾阶段均表现为抑制作用, 随着电场强度的增大, 核态沸腾和膜态沸腾的热流密度有所下降. 为了揭示导致这一现象的原因, 图6为$L_H^* = 6.25$, 电场强度$E = 0$, $E = 0.16327$作用下, 核态沸腾(图6(a))和膜态沸腾(图6(b))气泡形态对比. 从图6(a)可以发现, 核态沸腾时电场强度$E = 0.16327$作用下的气泡尺寸比$E = 0$的气泡尺寸小. 沸腾过程主要通过饱和液体相变吸热过程进行热量传递, 气泡尺寸减小意味着由饱和液体相变为气体的总质量减小, 即潜热换热减小, 从而导致核态沸腾阶段热流密度降低. 从图6(b)可以发现，在膜态沸腾阶段, 蒸气在加热表面铺展成膜. 无电场作用时($E = 0$)气体由两边向加热器中心聚集并在浮力的作用下向上运动, 在加热表面中心位置形成向上的蒸气柱, 铺展在加热表面蒸气膜的气液界面比较平滑. 当存在电场作用时($E = 0.16327$), 同样在加热表面形成蒸气膜, 但在电场力的作用下, 铺展在加热表面蒸气膜的气液界面发生形变, 蒸气膜更加均匀地铺展在加热表面, 增大了蒸气膜的整体厚度, 导致气膜施加在整个加热表面的热阻增大, 膜态沸腾的热流密度降低. 可以得出结论: 当加热器尺寸较小时, 核态沸腾阶段均匀电场减小气泡尺寸, 降低了沸腾换热过程中的潜热换热, 从而降低热流密度; 膜态沸腾阶段均匀电场使加热表面蒸气膜更加均匀地铺展, 增大了热阻从而降低膜态沸腾阶段热流密度.

观察图5(c)—(f)可以发现, $L_H^*\geqslant 7.5$时, 均匀电场对整个沸腾过程的影响与壁面过热度相关. 当过热度较低时, 均匀电场抑制沸腾, 当过热度较高特别是沸腾处于核态沸腾向膜态沸腾过渡的阶段, 均匀电场对沸腾表现为促进作用, 且CHF均有所提高. 为了分析不同长度加热器下均匀电场提高热流密度的机理, 图7列出了加热器长度分别为$L_H^*$ = 7.5, 9.375, 12.5, 18.75, $E = 0.10884$与$E = 0$条件下同一时刻加热表面气泡形态对比图, 图中加热温度为电场强度$E = 0.10884$时CHF点所对应温度. 对比图5可知, 当电场强度为零, $L_H^*$ = 7.5, 9.375, 12.5三种情况下沸腾状态为过渡沸腾, $L_H^*$ = 18.75时沸腾状态为膜态沸腾. 而当存在电场力时($E = 0.10884$), 4种加热器长度情况下加热表面有多个气泡产生, 且气泡相互作用力较弱, 说明电场力延迟了临界热流密度. 为此, 可以确定当$L_H^*$ = 7.5, 9.375, 12.5, 18.75时, 均匀电场通过阻止气泡在加热表面合并, 将过渡沸腾状态甚至是膜态沸腾状态转变为核态沸腾, 使沸腾在更高的过热度达到CHF, 从而达到提高CHF的效果. 这也进一步解释了当$L_H^*$ = 3.125, 6.25时电场的存在无法提高热流密度和CHF, 反而会抑制沸腾. 因为当加热器尺寸较小时, 加热表面无法生成多个气泡, 因此不存在气泡在加热表面合并的现象, 此情况下均匀电场对沸腾的影响主要表现为减小气泡尺寸如图6(a)所示.

为了进一步研究在均匀电场作用下, 加热器长度与CHF的关系, 图8展示了不同加热器长度下均匀电场对CHF提高的百分比. 首先, 在相同加热器长度下, 随着电场强度的提高, CHF提高的百分比越来越大. 此现象与图5和图7中规律一致, 即电场强度越大, 加热表面气泡越不易合并, 会使沸腾在更高过热度下达到CHF, 从而使CHF越来越大. 另外, 当$L_H^*\leqslant 9.375$时, 在同一电场强度下, 随着加热器长度的增大, CHF提高的百分比逐渐增大; 当$9.375 < L_H^*\leqslant 18.75$, 在同一电场强度下, 随着加热器长度的增大, CHF提高的百分比逐渐减小. 为了解释这两种现象, 下面分别选取符合这两种条件的工况进行研究.

图9展示了在电场强度$E = 0.10884$, 加热器长度分别为$L_H^* = 7.5$和$L_H^* = 9.375$时, CHF对应加热温度($T_\mathrm{b} = 0.995 T_\mathrm{c}$)下空间平均热流密度随时间的变化, 为保证沸腾过程处于充分发展阶段, 图中展示了$t^* = 20$时刻以后的数据. 可以发现当$L_H^* = 7.5$时, $q^*$基本在0.10和0.20之间波动, 而$L_H^* = 9.375$的$q^*$整体处于0.15和0.20之间, 其沸腾传热性能提高. 为了进一步分析引起该现象的原因, 图10(a), (b)分别展示了在$t^* = 69.32$、$t^* = 71.46$时刻两种加热器长度下的气泡形态演变、流场分布和气泡所受电场力的方向以及位置(图中白色箭头). 从图中结果可以发现, 当$L_H^* = 7.5$时(如图10(a)所示)$t^* = 69.32$时刻加热表面有3个气泡生成, 这3个气泡在$t^* = 71.46$时刻合并为一个大气泡. 通过分析以上两个时刻的流场图和气泡的受力可以发现, 虽然均匀电场具有横向压缩气泡阻止气泡在加热表面合并的作用, 但由于加热台长度较小导致气泡间距很小, 并且根据流场的方向得知气泡倾向于向加热表面中心部分聚集, 最终仍然会出现气泡合并现象. 气泡在加热表面合并会导致其施加在加热表面的换热热阻增大, 使热流密度陡然降低. 当$L_H^* = 9.375$时(如图10(b)所示), $t^* = 69.32$时刻加热表面同样生成3个气泡, 由于此时加热台长度增大, 生成的3个气泡之间的间距较大, 在电场力作用下这3个气泡在$t^* = 71.46$时刻并未发生合并, 并且在加热器两端生成又生成两个新的气泡, 此时沸腾换热热阻较小, 并且由于新气泡的产生, 沸腾换热潜热增加, 从而使$q^*$维持在较高水平. 由此可以得出结论: 当$6.25 < L_H^*\leqslant 9.375$, 更长的加热器为加热表面气泡的生成提供了充分的空间, 均匀电场作用下的气泡间距增大、气泡数量增加. 加热表面气泡不易合并, 降低换热热阻, 同时增加了气化核心, 提高了潜热换热.

图11展示了在电场强度$E = 0.16327$, 加热器长度分别为$L_H^* = 12.5$, $L_H^* = 18.75$时, CHF对应加热温度($T_\mathrm{b} = 0.995 T_\mathrm{c}$)下空间平均热流密度随时间的变化. 同样为了保证沸腾到达了充分发展阶段, 图中仅展示了$t^*\geqslant 20$的数据. 从图中结果可以发现在$t^* = 77$之前, 两种长度下的$q^*$相差不大, 而当$t^*\geqslant 77$时, $L_H^* = 18.75$的$q^*$明显小于$L_H^* = 12.5$时的数值. 为进一步分析$L_H^* = 18.75$的$q^*$较小的原因, 图12展示了在相同工况下, 两种长度加热器在不同时刻的气泡形态、气泡所受电场力和流场分布. 从图12(a)可以发现, 当$L_H^* = 12.5$, 加热表面生成的气泡向加热器中心移动, 且电场力主要作用于气液界面并指向气体内部, 此时气泡更容易离开加热表面. 当$L_H^* = 18.75$(如图12(b)所示)时, 由于加热器较长, 平行于加热表面的流线较长(图12(b)的白色矩形区域为加热器右侧流线与加热表面平行区域), 气泡运动趋势与图12(a)一致. 同时由于加热表面附近流体密度变化梯度大, 在加热表面存在较多方向向下的电场力, 使得部分蒸气被紧贴于加热表面, 并在气泡根部形成逆时针方向流动的涡, 此现象不利于加热表面气泡向中间移动, 而且贴于加热表面的蒸气增大了固体与流体之间的热阻, 最终出现$q^*$降低的现象. 因此可以得出结论: 由于$L_H^* = 18.75$的加热器存在更长的平行于加热表面的流线, 此区域蒸气在电场力作用下容易被紧贴于加热表面, 增大了固体与流体之间的换热热阻, 并在气泡根部形成不利于气泡向中间移动的涡, 减缓了加热表面热流体与两侧较冷流体的热质交换, 从而导致均匀电场对CHF的提高百分比低于$L_H^* = 12.5$的加热器.

根据前面的分析和讨论, 对于亲水表面, 均匀电场对不同长度加热器池沸腾的影响不尽相同. 当$L_H^*\leqslant 6.25$时, 由于加热表面产生气泡较少, 均匀电场对池沸腾的影响表现为减小气泡尺寸从而抑制沸腾. 当$7.5\leqslant L_H^*\leqslant 18.75$时, 均匀电场总是能够提高CHF, 且CHF随着电场强度的增大而增大. 另一方面, 相同均匀电场强度对不同长度加热器池沸腾的CHF提高的百分比也有较大的差异. $L_H^*\leqslant 9.375$时, 长度的增大为电场力作用下的加热表面气泡更大的间距提供了空间条件, 气泡不宜合并, 均匀电场对CHF提高的百分比逐渐增大. $L_H^* > 9.375$时, 过大的加热器长度具有更长的水平流线区域, 均匀电场作用下此区域蒸气容易贴附于加热表面, 沸腾换热热阻增大, 导致CHF提高的百分比降低.

4.3   均匀电场作用下气泡行为的变化对疏水表面沸腾传热的影响

本节研究疏水表面均匀电场作用下, 加热器长度对池沸腾传热性能的影响. 为了与亲水表面的结果做对比, 这里加热器长度和均匀电场强度与上节亲水表面情况的数据相同, 而疏水表面润湿性取$\theta = 120^\circ$.

图13为疏水表面不同长度加热器在不同电场强度下的沸腾曲线. 观察发现当$L_H^* = 3.125$时均匀电场对整个沸腾阶段均有抑制作用, 当$L_H^*\geqslant 6.25$时, 均匀电场对较高过热度下的沸腾有促进作用, 并且其此时CHF随电场强度的增大而有所均有提高, 这与亲水表面得到的结果类似. 不同的是当$L_H^* = 6.25$时亲水表面的CHF被抑制(图5(b)), 而疏水表面的CHF被促进. 为了分析引起这一现象的原因, 图14展示了$L_H^* = 6.25$时, 加热温度$T_\mathrm{b} = 0.99 T_\mathrm{c}$的疏水表面气泡动力学行为. 从图中结果可以发现当电场强度为零时, 由于疏水表面相变速率较大^[33], $t^* = 80.39$时刻形成了一个占据大部分加热表面的气泡, 蒸气施加在加热表面的热阻较大, 导致热流密度减小. 当$E = 0.16327$, 在相同时刻加热表面的气泡在电场力作用下发生横向挤压形变, 并在两侧生成两个小气泡, 成核点的增加进一步使沸腾换热性能增强, 因此CHF被促进.

图15给出了不同电场强度下疏水表面加热器长度与CHF的关系. 可以发现在同一加热器长度下, CHF随着电场强度的增大而增大. 相同电场强度下, 当$L_H^* < 12.5$, CHF提高百分比随着长度的增大逐渐增大, 而当$L_H^*\geqslant 12.5$, 随着长度的增大逐渐减小. 对比图8发现, 亲水表面和疏水表面关于加热器长度与CHF的关系具有相同的趋势, 不同的是CHF提高百分比的最大值所对应的加热器长度不同, 如$E = 0.10884$时, 亲水表面的阈值为$L_H^* = 9.375$, 疏水表面的阈值为$L_H^* = 12.5$. 另外, 均匀电场对疏水表面CHF的提高程度远高于亲水表面, 比如当$E = 0.16327$时, 所有列举的加热器长度, 亲水表面最大提高25.30%, 疏水表面最大提高73.07% . 这是因为一方面疏水表面的CHF本身就小于亲水表面; 另一方面, 疏水表面的气泡基底比亲水表面大, 气泡容易在加热表面合并, 而均匀电场提高CHF的机理为通过施加于气液界面的电场力阻止气泡合并, 使沸腾在更高过热度下达到CHF. 这两个原因导致均匀电场对疏水表面的CHF提高程度大于亲水表面.

5.   结　论

采用格子Boltzmann直接模拟数值方法研究了均匀电场作用下不同加热器表面润湿性以及加热器长度的池沸腾传热过程. 重点研究了气泡行为变化、电场强度、加热器长度以及临界热流密度之间的关系, 并通过气泡动力学行为、气泡所受电场力和流场分布对其内在机理进行综合分析, 得出以下结论.

1) 无电场作用时, 在沸腾过程中生成的气泡之间的相互作用随着加热器长度的增大先增大再不变, 对应的临界热流密度也随着加热器长度的增大先降低再趋于稳定, 且亲水表面临界热流密度保持稳定时对应的加热器长度值小于疏水表面的情况.

2) 当存在均匀电场时, 不论是亲水表面还是疏水表面, 在核态沸腾阶段初期, 电场的存在会小幅度降低热流密度. 然而在核态沸腾阶段的中后期, 电场的存在会增大热流密度, 且疏水表面电场作用增大的热流密度比例大于亲水表面的情况.

3) 随着加热器长度的增大, 均匀电场对亲水表面和疏水表面CHF提高百分比均先增大后减小. 电场强度E = 0.05442, 0.10884, 0.16327时, 亲水表面CHF分别在$L_H^*$ = 12.5, 9.375, 9.375达到增大的最大百分比, 而疏水表面在$L_H^*$ = 9.375, 12.5, 12.5达到增大的最大百分比.

本文的研究在微观尺度揭示了电场力和加热器特性对池沸腾换热性能的影响机理, 所得到的结论可为微电子器械加热器的设计提供理论基础.