相变传热广泛存在于自然界及工业生产过程，其中液滴蒸发是相变传热研究的重要分支。近年来，随着液滴蒸发在农业、医疗、冷却、印刷、镀膜等方面的广泛应用^［1−2］，液滴蒸发引起了越来越多学者的关注。国内外学者通过实验观察、理论推导和数值模拟对液滴蒸发进行了一系列的研究，得到了很多有意义的研究结果。

在过去的几十年里，格子玻尔兹曼方法（LBM）已经发展成为流动和传热领域一种非常有效的数值模拟方法^［3］。LBM由于具有边界处理简单、并行计算等优点，被广泛应用于模拟单相流、多相流、相变^［4−5］等流动和传热现象。目前存在多种多相流LBM模型，包括由Shan等^［6］提出的伪势模型，由Gunstensen等^［7］提出的颜色模型，由Swift等^［8］提出的自由能模型，以及由He等^［9］提出的相场模型。近年来，越来越多的学者致力于采用格子玻尔兹曼方法模拟蒸发、沸腾、冷凝等气液相变现象。Sun等^［10］采用一种三维混合格子玻尔兹曼模型模拟了水平加热壁面上单个气泡和多个气泡的生长等相变现象。Hazi等^［11］基于熵的能量方程和非理想气体的热力学关系导出了目标温度方程，提出了一种改进的伪势模型，给出了相应的热格子玻尔兹曼方程，并采用该模型模拟了缓慢流动的流体在水平平板上的沸腾等相变现象。Hatani等^［12］采用Lee提出的基于Cahn−Hilliard扩散界面理论的格子玻尔兹曼多相模型模拟了气体的冷凝等相变现象。Gong等^［13］基于单弛豫时间的伪势模型提出了一种新的模型用于模拟恒温壁面上和恒热流加热壁面上气泡的产生、长大和脱离壁面的相变现象，而后又模拟了在不同浸润性壁面上气泡的产生等相变现象。Li等^［14］采用多弛豫时间的伪势模型和有限差分法的耦合模型模拟了核态沸腾以及液滴的蒸发等相变现象。

在前人研究基础上，本文采用单组分伪势模型与有限差分耦合的混合热格子玻尔兹曼模型(TLBM)，采用真实流体P−R状态方程，对液滴的蒸发过程进行了研究。首先，通过模拟方腔内液滴蒸发过程，验证了所采用方法及程序的有效性。模拟了液滴撞击并附着在高温壁面后的蒸发相变过程，分析了壁面温度、液滴邦德数和雷诺数对液滴撞击壁面后蒸发相变过程的影响。

1 混合热格子玻尔兹曼模型

采用Gong等^［13］提出的伪势模型模拟流体的流动。该模型中流体粒子分布函数的演化方程为

$ \begin{split}&\quad\quad{f_i}(x + {e_i}{\delta _t},t + {\delta _t}) - {f_i}(x,t) = \\ &\quad\quad\;\;\;\;- \frac{1}{\tau }\left[ {{f_i}(x,t) - f_i^{\rm eq}(x,t)} \right] + \Delta {f_i}(x,t)\end{split} $

(1)

式中： ${f_i}(x,t)$ 为密度分布函数， $i$ 表示方向， $x$ 表示空间位置， $t$ 为格子时间； $\tau $ 为弛豫时间； ${\delta _t}$ 为时间步长； ${e_i}$ 为离散速度； $t + {\delta _t}$ 表示迭代中的下一个时刻； $x + e{}_i{\delta _t}$ 表示 $x$ 节点的相邻位置； $f_i^{\rm{eq}}(x,t)$ 为密度的平衡态分布函数； $\Delta f_i(x,t)$ 为体积力项。

对于单组分系统，流体黏度 $\nu $ 与弛豫时间 $\tau $ 的关系为

$ \qquad \nu = c_{\rm s}^2\left( {\tau - 0.5} \right){\delta _t} $

(2)

式中， ${c_{\rm s}}$ 为格子声速，且 $c_{\rm s}^2 = {{{c^2}} / 3}$ ，此处 $c = {{{\delta _x}} / {{\delta _t}}}$ ， ${\delta _x}$ 为空间步长，其值为1.0。

式（1）中 $f_i^{\rm eq}(x,t)$ 表达方程为

$ \quad\quad f_i^{\rm eq}(x,t) = {\omega _i}\rho \left(1 + \frac{{{e_i} \cdot u}}{{c_{\rm s}^2}} + \frac{{{{({e_i} \cdot u)}^2}}}{{2c_{\rm s}^4}} - \frac{{{u^2}}}{{2c_{\rm s}^2}}\right) $

(3)

式中： ${\omega _i}$ 为权系数； $u$ 为流体粒子的速度； $\rho $ 为流体密度。

$\rho $ 和 $\rho u $ 的计算式为

$ \quad\quad\rho = \sum\limits_{i = 0}^8 {{f_i}},\;\;\;\;\rho u = \sum\limits_{i = 0}^8 {{e_i}{f_i}} $

(4)

此处 $u$ 并非真实宏观速度 $U$ 。

$ \quad\quad\rho U = \sum\limits_{i = 0}^8 {{e_i}{f_i} + \frac{{{\delta _t}}}{2}} F $

(5)

式中，F为流体受到的合力。

式（1）中 $\Delta {f_i}(x,t)$ 是通过精确差分方法^［15］进行求解，即

$ \quad\quad\Delta {f_i}(x,t) = f_i^{\rm eq}[\rho (x,t),u + \Delta u] - f_i^{\rm eq}[\rho (x,t),u]\!\!\!\! $

(6)

式中， $\Delta u = F {{{\delta _t}} / \rho }$ 为速度变量，流体受到的合力 $F$ 由流体粒子间的相互作用力 ${F_{\rm{int} }}$ 和流体粒子受到的重力 ${F_{\rm g}}$ 组成。

${F_{\rm{int} }}$ 采用Gong等^［13］提出的形式，即

$ \begin{split} \quad\quad {F_{\rm{int} }} =& - \beta G\psi \left( x \right)\sum\limits_{x'} {\psi (x')(x' - x) - }\\ &{\frac{{1 - \beta }}{2}} G\sum\limits_{x'} {{\psi ^2}(x')(x' - x)} \end{split} $

(7)

式中： $\beta $ 为用于提高计算稳定性的权重因子； $G$ 为相互作用力系数； $x'$ 为 $x$ 的相邻格子点坐标； $\psi \left( x \right)$ 为和密度有关的有效质量。

$ \qquad\psi (x) = \sqrt {\frac{2(p - \rho c_{\rm s}^2)}{{c_0}g}} $

(8)

式中： $g$ 为粒子间作用强度系数； ${c_0}$ 为格子结构确定系数； $p$ 为流体压力。

本文采用的真实流体P−R状态方程^［13］为

$\qquad\left\{ \begin{align} & {{p}_{{}}}=\frac{\rho RT}{1-b\rho }-\frac{a{{\rho }^{2}}\varepsilon (T)}{1+2b\rho -{{b}^{2}}{{\rho }^{2}}} \\ & \varepsilon (T)=\Big{[}1+(0.374\ 64+1.542\ 26\omega - \\ & \qquad0.269\ 92{{\omega }^{2}})(1-\sqrt{T/{{T}_{\rm{c}}}})\Big{]} \\ & a=0.457\ 24{{R}^{2}}{{T}_{\rm{c}}}^{2}/{{p}_{\rm{c}}} \\ & b=0.077\ 8R{{T}_{\rm{c}}}/{{p}_{\rm{c}}} \end{align} \right.$

(9)

式中： $R$ 为气体常数； $a$ 为引力参数； $b$ 为斥力参数； $T$ 为温度； $\omega $ 为偏心因子； $ {T_{\rm c}} $ 、 ${p_{\rm c}}$ 分别为临界温度和临界压力。

流体粒子受到的重力表示为

$ \qquad{F_{\rm g}} = (\rho - {\rho _{\rm ave}}){g_{\rm s}} $

(10)

式中： ${\rho _{\rm ave}}$ 为计算域内流体的平均密度； ${g_{\rm s}}$ 为重力加速度。

另外，根据Ba等^［16］和Liu等^［17］所采用的方法，计算壁面一层的密度，并给定壁面浸润性。

除了流体流动的模型，本文参考Márkus等^［18］使用的能量方程模拟相变过程，即

$ \;\;{\partial _t}T = \frac{1}{{\rho {c_v}}}\nabla \cdot (\lambda \nabla T) - U \cdot \nabla T - \frac{T}{{\rho {c_v}}}\left( {\frac{{\partial p}}{{\partial T}}} \right)\nabla \cdot U $

(11)

式中： ${\partial _t}T$ 为在当前时刻对温度求偏导； $\nabla T$ 为温度梯度； $\lambda $ 为导热系数； ${c_v}$ 为定压比热容。

将式(11)中等号右边定义为 $K(T)$ 。对于该方程的求解，空间离散化采用各向同性中心差分法^［19］，对于一个变量 $\phi $ ，其空间梯度和拉普拉斯算子可以采用二阶各向同性差分法求解，即

$ \qquad {\partial _i}\phi (x) \approx \frac{1}{{c_{\rm s}^2{\delta _t}}}\sum\limits_i {{\omega _i}} \phi \left( {x + {e_i}{\delta _t}} \right){e_i} $

(12)

$ \qquad {\nabla ^2}\phi (x) \approx \frac{2}{{c_{\rm s}^2{\delta _t}^2}}\sum\limits_i {{\omega _i}\left[ {\phi \left( {x + {e_i}{\delta _t}} \right) - \phi (x)} \right]} $

(13)

时间离散化采用四阶龙格库塔法^［14］，即

$ \qquad \left\{ \begin{array}{l} {T^{t + \delta_t}} = {T^t} + \dfrac{{{\delta _t}}}{6}({h_1} + 2{h_2} + 2{h_3} + {h_4}) \\ {h_1} = K\left( {{T^t}} \right) \\ {h_2} = K\left( {{T^t} + \dfrac{{{\delta _t}}}{2}{h_1}} \right) \\ {h_3} = K\left( {{T^t} + \dfrac{{{\delta _t}}}{2}{h_2}} \right) \\ {h_4} = K\left( {{T^t} + {\delta _t}{h_3}} \right) \end{array}\right. $

(14)

式中， ${T^t}$ 、 $T^{t + \delta_t}$ 分别为当前时刻、下一时刻的温度。

2 模拟结果与分析 2.1 液滴蒸发的 ${D^2}$ 定律

为了验证本文所采用方法及程序的有效性，首先对方形区域内液滴蒸发过程进行模拟，考察液滴直径随时间的变化是否符合 ${D^2}$ 定律^［20］。

本算例中重力为0，导热系数和定压比热容为常数。计算域为方形区域，网格数151 × 151，初始化计算域内液滴直径 ${D_0}$ 为40个格子单位，液滴温度设置为饱和温度 ${T_{\rm s}} = 0.86{T_{\rm c}}$ ，周围汽相温度和边界温度均设置为高温 ${T_{\rm v}} = 0.90{T_{\rm c}}$ 。

图1给出了 ${\left( {{D / {{D_0}}}} \right)^2}$ 随时间的变化，其中D为蒸发过程中液滴直径， ${D_0}$ 为液滴初始直径。从图中可以看到， ${\left( {{D / {{D_0}}}} \right)^2}$ 随时间几乎呈线性变化，这与Nishiwaki等^［21］的实验结果相吻合，也印证了本文所采用模型的正确性。

2.2 液滴撞击高温壁面并蒸发

计算域为方形区域，网格数为151 × 151，左右边界设置为周期性边界条件，上边界设置为压力出口，下边界设置为高温壁面。

为了更好地描述撞击高温壁面液滴的物性，引入雷诺数Re和邦德数Bo，其中Re和液滴黏性有关，Re越大，液滴黏性越小。

$\qquad Re = \frac{{{u_0}{D_0}}}{\nu }$

式中： $\nu $ 为液滴运动黏度； ${u_0}$ 为液滴初始速度。

邦德数Bo和液滴所受重力有关，液滴所受到的重力越大，Bo越大。

$\qquad Bo = \frac{{\left| {{g_{\rm s}}} \right|\rho {D_0}^2}}{\sigma }$

式中， $\sigma $ 为液滴的表面张力。

为了更加清晰地描述液滴的蒸发过程，引入液滴蒸发速率 ${V_{\rm s}}$ 来反映液滴蒸发的快慢。

$\qquad {V_{\rm s}} = \frac{{{V_0}}}{{{t_{\rm z}}}} $

式中： ${V_0}$ 为液滴蒸发过程中的体积变化量； ${t_{\rm z}}$ 为液滴蒸发完全所需的时间。

液滴在蒸发过程中，若液滴初始体积相同，蒸发完全时间越短，蒸发速率越大。

该算例中导热系数 $\lambda $ 设置为 $\lambda = \rho {c_v}\chi $ ，其中 $\chi $ 为热扩散率。计算中设置底面的浸润性保持不变，即液滴撞击壁面后在壁面上蒸发的过程中浸润角始终保持 ${113^ \circ }$ ，但液滴与壁面的接触直径不断减小，初始化液滴直径大约为40个格子单位， $Re = 10.8$ ， $Bo = 15.15$ ，液滴温度为 $0.86{T_{\rm c}}$ ，下壁面温度为 $0.896{T_{\rm c}}$ 。液滴撞击高温壁面后蒸发过程如图2所示。初始液滴位于方形区域的中心位置，在重力作用下液滴下落并撞击高温壁面。从图2中可以看出，液滴撞击高温壁面后会先向周围铺展，然后再回缩，稳定后在高温壁面上蒸发，体积不断减小。

2.2.1 壁面温度对液滴撞击高温壁面后蒸发过程的影响

在 $Re = 10.8$ 、 $Bo = 15.15$ 、温度为 $0.86{T_{\rm c}}$ 、初始直径相同的液滴撞击不同温度的高温壁面的条件下，分析壁面温度对液滴撞击高温壁面后蒸发过程的影响。高温壁面与液滴之间温差 $\Delta T$ 分别取 $0.036{T_{\rm c}}$ 、 $0.038{T_{\rm c}}$ 、 $0.040{T_{\rm c}}$ 。图3为液滴撞击不同温度高温壁面后液滴体积随时间的变化。可以看到， $\Delta T$ 越大，液滴蒸发速率越大。这是由于下壁面温度越高，液滴与高温壁面的换热越剧烈，因此液滴撞击壁面后的蒸发速率越大。

2.2.2 重力对液滴撞击高温壁面后蒸发过程的影响

考虑 $Re = 10.8$ 、温度为 $0.86{T_{\rm c}}$ 、初始直径相同，Bo分别为15.15、20.19、25.24的液滴撞击温度为 $0.896{T_{\rm c}}$ 的高温壁面，分析了不同Bo的液滴对液滴撞击高温壁面后蒸发过程的影响。为了便于分析，定义液滴撞击壁面后与壁面的接触直径为铺展直径 ${D_{\rm s}}$ ，其示意图如图4所示。图5为不同邦德数的液滴撞击同一温度壁面过程中 ${D_{\rm s}}$ 和液滴体积随时间的变化。由图5（a）中可以看到，邦德数越大，液滴撞击到壁面后的最大铺展直径越大。同时，邦德数越大，液滴撞击到壁面所需的时间越短。这是因为液滴所受重力越大，邦德数越大，进而液滴下落速度就越快。从图5（b）中可以看到，随着邦德数的增加，液滴撞击壁面后的蒸发速率也增加。这是由于邦德数的增加导致液滴撞击壁面后最大铺展直径增加，即液滴与高温壁面接触得更多，换热面积增大，换热增强，所以增大了液滴在壁面上的蒸发速率。

2.2.3 黏性对液滴撞击高温壁面后蒸发过程的影响

考虑 $Bo = 15.15$ 、温度为 $0.86{T_{\rm c}}$ 、初始直径相同，Re分别为10.08、12.60、16.80的液滴撞击温度为 $0.896{T_{\rm c}}$ 的高温壁面，雷诺数越大，液滴的黏性就越小。分析了不同雷诺数的液滴对液滴撞击高温壁面后蒸发过程的影响。图6给出了不同雷诺数的液滴撞击同一高温壁面后 ${D_{\rm s}}$ 、液滴体积随时间的变化。从图6（a）中可以看到，由于重力相同，所以液滴撞击到壁面所需时间几乎相同。另外，液滴雷诺数越大，液滴黏性越小，液滴撞击到壁面上的最大铺展直径越大。从图6（b）中可以看到，随着雷诺数的增大，液滴撞击到壁面后蒸发速率也增大。这同样是由于雷诺数越大，液滴的黏性越小，液滴撞击到壁面后与高温壁面的接触直径越大，即换热面积增大，换热增强，液滴蒸发得越快。

3 结　论

本文采用混合热格子玻尔兹曼模型模拟了液滴在方形区域内的蒸发过程以及液滴撞击高温壁面并蒸发的过程，得出以下结论：

(1)在方形区域内的蒸发过程中，液滴直径随时间呈线性变化，这与实验结果相吻合。

(2)在液滴撞击高温壁面并蒸发的过程中，壁面温度对液滴撞击壁面后的蒸发速率有显著的影响。壁面温度越高，液滴在撞击壁面后与壁面的换热就越剧烈，液滴在壁面上的蒸发速率就越大。

(3)在液滴撞击高温壁面并蒸发的过程中，重力越大或液滴黏性越小，液滴撞击壁面后的最大铺展直径越大，即液滴与壁面的换热面积就越大，换热就越强，液滴在壁面上的蒸发速率就越大。