摘要:采用Level Set-VOF方法建立单液氮滴撞击壁面的数值模型, 探索壁面润湿性(30°—150°)、撞击速度(0.1和1.6 m/s)及壁面温度(300—500 K)对液滴撞壁演化过程中相变行为的影响, 并理论推导了气膜生长数学模型. 结果表明: 增强壁面润湿性、提高撞击速度有利于液滴沿径向铺展, 从而增大了换热面积并降低热阻, 使换热性能得到显著提升; 提高壁面温度增大了换热温差, 热流密度随之上升; 三相接触线处热阻较小导致边缘处热流密度高于中心处, 不同润湿壁面上热流分布的差异性因初始速度的增大而缩小, 呈现明显的速度效应; 在膜沸腾区, 传热过程主要集中在撞击初期, 气膜是主要换热热阻; 基于质量守恒和能量守恒建立气膜生长数值模型, 模型预测结果与本文模拟结果和其他研究结果非常吻合.
关键词: Leidenfrost/
液氮/
液滴撞击/
Level Set模拟

English Abstract

全文HTML

--> --> -->

喷雾冷却技术作为一种高效的热控制技术引起众多****的广泛关注, 其冷却过程涉及单液滴撞击壁面、液滴与液滴、液滴与液膜之间的相互作用, 属于复杂的多相流问题^[1,2].
由于液氮具有冷却温度低、成本低和无残留等优点成为低温控制领域(低温风洞、低温试验条件的开发、大功率电子设备的热管理、食品速冻以及低温手术等)最理想的工质之一^[3-6]. 然而液氮沸点较低, 雾化液滴与壁面接触后瞬间蒸发并在液滴底部和壁面之间产生一层连续的气膜, 发生膜沸腾, 即出现Leidenfrost现象^[7,8], 热阻较大的气膜往往会导致换热效果剧烈下降. 研究低温液滴撞击壁面的铺展动态及相变行为过程有助于深入理解低温喷雾冷却技术的换热机理, 尤其是发生Leidenfrost效应时的热传输机理.
液滴撞壁是跨尺度、多物理场耦合控制的强瞬变过程, 其影响因素除液滴自身的热物理性质外, 还包括热沉面结构或理化性质等方面^[9]. 许多****利用理论分析、实验研究及数值模拟手段开展了大量细致的研究工作. Karl和Frohn^[10]研究了在Leidenfrost温度以上液滴撞击热壁的力学行为, 根据实验结果建立的关系式可用于改善两相流的数值模拟, 并且从理论方面利用动量损失的相关性, 推导出液滴最大变形的理论近似值, 其测量结果与其他文献报道的结果非常吻合. Scheller和Bousfield^[11]通过实验研究撞击条件(初始直径、初始速度)和壁面条件(壁面种类)的改变对撞击行为的影响, 得出液滴撞击壁面扩展过程不发生回缩时无量纲润湿长度最大值的经验公式.
在液滴物性方面, 刘海龙等^[12]采用碳纳米管、石墨烯、纳米石墨粉制备三种稳定的纳米流体, 利用显微高速数码摄像技术捕捉上述三种纳米流体撞击固体壁面的动态过程, 结果表明, 三种纳米颗粒均能促使基液表现出显著的剪切变稀特性; 流体的剪切黏度在液滴撞击壁面的扩展过程中起到了重要的作用; 液滴的无量纲高度和铺展因子的变化与纳米流体剪切黏度的变化呈负相关; 在液滴的最初铺展阶段, 惯性力起主导作用, 随着液滴撞击速度的增大液滴的扩展范围和速度亦增大. 沈胜强等^[13]采用高速摄像仪观察了水和乙醇液滴撞击不锈钢壁面的蒸发过程, 分析了液滴撞击加热壁面后的蒸发特性参数, 实验得到水滴的平均热流密度为0.014—0.110 W·mm^–2.
在热沉面结构或理化性质等方面, 已有研究表明热沉面润湿性对液滴的蒸发和二次雾化影响较大, ?ikalo等^[14]研究发现已有的动态接触角经验模型(Hoffman-Voinov-Tanner模型)不能很好地预测动态接触角的变化, 特别是在毛细数较大的情况下. 梁刚涛等^[15]实验研究低速液滴撞击润湿球面的现象, 发现液滴撞击润湿球面上发生反弹、局部反弹和铺展等现象, 考察了黏度对撞击过程的影响, 定量讨论了液滴铺展特征参数随撞击速度、球体直径和黏度的变化规律. 另外, 壁面温度是影响相变过程的直接因素. Tran等^[16]利用高速摄像机拍摄水滴(0.5 ≤ We ≤ 500)撞击高温光滑壁面(250 ℃ ≤ T ≤ 560 ℃)过程呈现3个不同的机制: 接触沸腾、温和膜沸腾、雾化膜沸腾. 温和膜沸腾、雾化膜沸腾机制下的液滴无量纲最大扩展因子与We^2/5成正比, 随着We增加T_DL也随之增大; 从温和膜沸腾到雾化膜沸腾的过渡温度与液膜内的沸腾气泡有关, 并且随着We的增加过渡温度降低.
综上所述, 国内外****对液滴撞击壁面的实验和数值研究, 多集中于液滴撞壁动力学特性, 如三相接触线行为、动态接触角变化、自由表面流动等, 而往往忽略液滴撞击壁面的相变行为, 尤其是缺乏低温液滴撞壁后的热传输机制研究.
本文利用复合Level Set-VOF相界面追踪的方法, 建立单个液氮液滴撞击不同润湿壁面的计算模型, 并进行了模型的可靠性验证. 进一步从微观角度探究壁面热流密度分布和气膜厚度随壁面润湿性、初始速度和壁面温度等因素的变化规律, 进而揭示Leidenfrost效应下低温液滴撞击壁面时的运动形态和蒸发换热耦合机制.

2

-->

2.1.复合Level Set-VOF相界面追踪方法

采用复合Level Set-VOF方法追踪不可压缩两相流相界面, 可有效解决VOF和Level Set两种气液两相流动相界面追踪方法存在的计算收敛性、稳定性及准确性不高等问题^[17]. 复合Level Set-VOF相界面追踪方法计算过程包括以下几个步骤: Level Set函数和VOF函数初始化、流动控制方程求解、Level Set函数和VOF相函数($\phi $函数和F函数)的对流输运方程求解、相界面重构以及函数再次初始化等. 而气液相界面重构和函数再次初始化使得函数满足质量守恒特性, 并且能够更加精确地计算表面张力. 图1为复合Level Set-VOF气液两相流动计算流程图, 虚线部分为气液相界面追踪和壁面润湿模型的计算流程.
图 1 复合Level Set-VOF相界面追踪方法计算流程图
Figure1. Coupled Level Set-VOF phase interface tracking method calculation flow chart.

复合Level Set-VOF法通过$\phi $函数和F函数共同构造相界面, 采用分段线性界面重构思想, 在法向方向上移动界面使得单元液态区域面积比率和F函数值相匹配, 单元中心到相界面的垂直距离由割线法迭代求解得到.
界面法向量${{n}}$计算方法为

${{n}} = {{\nabla \phi }}/{{\left| {\nabla \phi } \right|}}.$

界面曲率$\kappa \left( \phi \right)$的表达式为

$\kappa \left( \phi \right) = \nabla \cdot ({{\nabla \phi }}/{{\left| {\nabla \phi } \right|}}).$

复合Level Set-VOF方法的相界面的追踪主要通过求解相函数对流输运方程来实现.

$\frac{{{\rm{D}}F}}{{{\rm{D}}t}} = \frac{{\partial F}}{{\partial t}} + ({{u}} \cdot \nabla )F = 0,$

$\frac{{{\rm{D}}\phi }}{{{\rm{D}}t}} = \frac{{\partial \phi }}{{\partial t}} + ({{u}} \cdot \nabla )\phi = 0,$

式中${{u}}$为速度. $\phi $函数经过对流输运方程求解后, 将不再保持距离函数的性质, 因此必须对$\phi $函数重新初始化, 再次重新初始化主要包括$\phi $函数符号以及$\phi $函数值的确定.
$\phi $函数值是单元中心到相界面的最小距离d, 而$\phi $函数的符号由F函数值决定:

${S^\phi } = {\rm{sign}}(0.5 - F).$

当F < 0.5时, 单元中心在气态区域, $\phi $函数为正; 当F > 0.5时, 单元中心在液态区域, $\phi $函数为负; 当F = 0.5时, 表示相界面, $\phi = 0$; $\rm sign(\cdot)$是符号函数.
2 -->

2.2.控制方程

1)连续性方程和动量方程

$\nabla \cdot {{u}} = 0,$

$\begin{split} &\rho \left( \phi \right)\left( {\partial {{u}}/\partial t + {{u}} \cdot \nabla {{u}}} \right) \\=\, & - \nabla p - {{{F}}_{\rm{s}}} + \nabla \cdot {{\tau }} + \rho \left( \phi \right){{g}},\end{split}$

式中, $\rho $为密度, p为压强, ${{\tau }}$为黏性应力张量, ${{g}}$为重力加速度, ${{{F}}_{\rm{s}}}$为源项.
采用CSF (continuum surface force)模型对表面张力进行求解, 在CSF模型中附加的表面张力通过在动量方程中添加源项的方式实现:

${{{F}}_{\rm{s}}} = \sigma \kappa \left( \phi \right)\nabla H\left( \phi \right),$

式中, $\sigma $为表面张力系数, $H\left( \phi \right)$为Heaviside函数.
Heaviside函数表达式如下:

$H( \phi ) = \begin{cases} {1,} \qquad\qquad {\phi \geqslant + h},\\{0,} \qquad\qquad {\phi \leqslant - h},\\0.5 + \dfrac{\phi}{3h} + \dfrac1{2\text{π}}\sin \left( \dfrac{2{\text{π}}\phi}{3h} \right),\\\qquad\qquad \quad {| \phi | \leqslant h},\end{cases}$

式中, h为网格大小. 黏性应力张量${{\tau }} = 2\mu \left( \phi \right){{S}}$, 应变率张量${{S}}$表达式为

${{S}} = \frac{1}{2}\left[ {\nabla {{u}} + {{\left( {\nabla {{u}}} \right)}^{\rm{T}}}} \right].$

不同区域$\rho $和$\mu $可由Heaviside函数计算:

$\rho \left( \phi \right) = {\rho _{\rm{l}}}\left[ {1 - H\left( \phi \right)} \right] + {\rho _{\rm{g}}}H\left( \phi \right),$

$\mu \left( \phi \right) = {\mu _{\rm{l}}}\left[ {1 - H\left( \phi \right)} \right] + {\mu _{\rm{g}}}H\left( \phi \right).$

2)能量方程

$\frac{{\partial \rho {c_{\rm{p}}}T}}{{\partial t}} + \nabla (\rho {c_{\rm{p}}}T{{u}}) = \nabla \cdot \lambda \nabla T,$

式中, ${c_{\rm{p}}}$为计算单元内定压热容, $\lambda $为计算单元内热导率, T为计算单元内温度.
3)壁面黏附
模型中还考虑了壁面黏附作用(wall adhesion), 如图2所示, 通过给壁面接触角赋值的形式来定义壁面的润湿性, 当液滴在固体表面上铺展时, 定义壁面静态接触角为θ, 用于调整壁面附近单元表面的法向量. ${{{n}}_{\rm{w}}}$和${{{\tau }}_{\rm{w}}}$分别为单位法向量和壁面切向量, 壁面黏附作用被定义在表面向量中:
图 2 壁面黏附
Figure2. Wall adhesion diagram.

${{{n}}_{\rm{s}}}={{{n}}_{\rm{w}}}\cos \theta +{{{\tau }}_{\rm{w}}}\sin \theta .$

2 -->

2.3.蒸发模型

液氮沸点较低, 当液氮滴撞击常温壁面时, 液滴撞壁过程存在相变, 需考虑液滴撞壁后的中相变行为. 本文采用应用最广泛的Lee^[18]模型, 其具体表达式如下.
蒸发过程,

${\dot m _{\rm{g}}} = \gamma {\alpha _{\rm{l}}}{\rho _{\rm{l}}}\frac{{T - {T_{{\rm{sat}}}}}}{{{T_{{\rm{sat}}}}}},~~~\;T > {T_{{\rm{sat}}}},$

冷凝过程

${\dot m _{\rm{l}}} = \gamma {\alpha _{\rm{g}}}{\rho _{\rm{g}}}\frac{{{T_{{\rm{sat}}}} - T}}{{{T_{{\rm{sat}}}}}},~~~\;T < {T_{{\rm{sat}}}},$

式中, ${\dot m _{\rm{g}}}$和${\dot m _{\rm{l}}}$分别为气相和液相的相变速率(kg·m^–3·s^–1); ${\alpha _{\rm{l}}}$和${\alpha _{\rm{g}}}$分别为液相和气相体积分数; ${\rho _{\rm{l}}}$和${\rho _{\rm{g}}}$分别为液相和气相密度; ${T_{{\rm{sat}}}}$为饱和温度; $\gamma $是控制相变的因子, 该因子的设定主要依据为使包含相界面的计算网格的温度维持在饱和温度, 本文$\gamma $设定为0.1 s^–1, de Schepper等^[19]和Alizadehdakhel等^[20]采用该模型对热管内的相变换热过程进行了模拟, $\gamma $同样设定为0.1 s^–1.
蒸发的气体质量作为质量源项离散后加载到控制方程中. 计算涉及的物性参数如密度、热导率、比热容等均取两相数学的平均值. 当蒸发发生在交界面时, Hertz Knudsen方程可以表达为

$F = \sqrt {\frac{M}{{2{\text{π}}R{T_{{\rm{sat}}}}}}} (p - {p_{{\rm{sat}}}}),$

式中, M为气体分子质量, R为理想气体常数. Clapeyron-Clausius建立了压力和饱和温度的状态关系:

$\frac{{{\rm{d}}p}}{{{\rm{d}}T}} = \frac{L}{{T\left({1}/{{{\rho _{\rm{v}}}}} - {1}/{{{\rho _{\rm{l}}}}}\right)}},$

式中L为汽化潜热. 当温度和压力接近饱和状态时, (18)式可表达为

$(p - {p_{{\rm{sat}}}}) = \frac{L}{{T\left(\dfrac{1}{{{\rho _{\rm{v}}}}} - \dfrac{1}{{{\rho _{\rm{l}}}}}\right)}}(T - {T_{{\rm{sat}}}}).$

则Hertz Knudsen方程可表达为

$F = \sqrt {\frac{M}{{2{\text{π}}R{T_{{\rm{sat}}}}}}} L\left( {\frac{{{\rho _{\rm{v}}}{\rho _{\rm{l}}}}}{{{\rho _{\rm{l}}} - {\rho _{\rm{v}}}}}} \right)\frac{{T - {T_{{\rm{sat}}}}}}{{{T_{{\rm{sat}}}}}}.$

为了正确添加源项, 将(20)式乘以单位体积内的表面积数${A_i} = {{6{\alpha _{\rm{g}}}{\alpha _{\rm{l}}}}}/{{{D_{{\rm{sm}}}}}}$, ${D_{{\rm{sm}}}}$为Sauter平均直径, 则蒸发源项可表达为

${\dot m _{\rm{g}}} = F{A_i} = \frac{{6{\alpha _{\rm{g}}}{\alpha _{\rm{l}}}}}{{{D_{{\rm{sm}}}}}}\sqrt {\frac{M}{{2{\text{π}}R{T_{{\rm{sat}}}}}}} L\left( {\frac{{{\rho _{\rm{v}}}{\rho _{\rm{l}}}}}{{{\rho _{\rm{l}}} - {\rho _{\rm{v}}}}}} \right)\frac{{T - {T_{{\rm{sat}}}}}}{{{T_{{\rm{sat}}}}}}.$

2 -->

2.4.模型设置及网格无关性验证

图3为单液氮液滴撞击壁面和发生Leidenfrost蒸发后的形貌演变图, 其计算域为3 mm × 6 mm的二维四边形均匀网格. 定义初始直径为D₀, 初始速度为U₀, 方向垂直壁面向下, 环境压力为0.1 MPa, 接触直径为D_t, 最大铺展直径为D_max, ${\delta _{\rm{t}}}$为中心气膜厚度, 考虑重力g的影响, 环境气体为氮气, 顶部及两侧为压力出口边界, 出口边界压力梯度为零, 底部为无滑移壁面.
图 3 初始时刻液滴撞击壁面几何模型及边界条件
Figure3. Model of droplet impact on wall at initial time.

实际喷雾液滴的几何形状较为复杂, 本文仅研究完美球形液滴的撞壁过程, 且初始时刻液滴底部与壁面相切; 由于液滴与环境及壁面之间的温差在222—422 K, 此时辐射换热较弱, 因此不考虑低温液氮液滴与环境、壁面之间的热辐射效应; 模拟液滴撞击速度较低, 环境压力为常压, 可将气液两相均视为不可压缩流体, 且忽略环境气体与液滴间的剪切力.
利用软件Fluent 17.0求解数值模型, 控制方程采用有限体积法进行离散, 压力与速度耦合采用Coupled方法, 压力求解采用PRESTO方法, 采用QUICK算法对Level Set方程和对流项进行求解, 采用Geometric-Reconstruction Scheme求解VOF方程, 动量和能量方程采用二阶迎风格式, 时间步长Δt = 10^–6 s, Δt内迭代次数为20. 计算相关参数见表1.

参数	数值
液氮密度/kg·m–3	806.88
液氮表面张力系数/N·m–1	8.22 × 10–3
液氮黏度/Pa·s	160.08 × 10–6
液氮滴初始速度/m·s–1	0.1—1.6
液氮滴初始直径/mm	0.5
液氮滴初始温度/K	78
壁面温度/K	300—500
液氮潜热/kJ·kg–1	199.3
氮气导热系数/W·cm–1·K–1	0.002475
氮气黏度/Pa·s	17.54 × 10–6
氮气密度/kg·m–3	1.16

表1相关参数
Table1.Related parameters.

网格无关性验证分别采用网格数200 × 400, 250 × 500, 300 × 600及350 × 700的四边形结构化网格, 具体工况为: 撞击初速度为0.1 m/s, 液滴直径为0.5 mm, 壁面温度为300 K, 壁面静态接触角为90°. 对比分析四种网格的铺展因子β随时间的变化关系, 结果如图4所示. 其中液滴铺展因子定义为液滴接触直径D_t与液滴初始直径D₀的比值:
图 4 网格无关性验证
Figure4. Verification of grid independence.

$\beta = {{{D_{\rm{t}}}}}/{{{D_0}}}.$

从图4可以看出, 随着网格数密度增大, 铺展因子β逐渐提高, 当网格数增加至300 × 600后, 铺展因子β增加趋势平缓, 网格数对铺展因子影响较小. 综合考虑计算精度及计算效率, 本文模拟采用300 × 600网格数.
2 -->

2.5.模型验证

为了对模型进行验证, 本文选取文献[21]中的两种润湿壁面上(静态接触角θ = 57°, θ = 25°)的液滴撞击实验数据进行模拟计算(工况: 壁面温度为500 ℃, 液滴初始直径为2 mm, 液滴初始温度为20 ℃, 液滴初始速度为2.05 m/s), 并将计算结果与实验结果进行对比(图5(a)). 对比文献[21]实验的结果可以看出, 随着时间推移, 铺展因子逐渐增大, 且在相同撞击速度和初始直径的情况下, 接触角较小的壁面上铺展因子要高于接触角较大的壁面, 两种工况下模拟结果与实验结果高度吻合; 另外, 为进一步验证模型的准确性, 选取文献[22]中一组数据进行液滴撞壁后的形貌演变对比(工况: 壁面温度384 ℃, 液滴初始直径为2 mm, 液滴初始温度为20 ℃, We = 2), 对比结果见图5(b), 可以看出液滴撞击光滑壁面后先后经历了铺展、回缩和反弹等过程, 模拟液滴撞壁后的运动行为与实验结果较吻合.
图 5 实验结果与模拟结果对比　(a)铺展因子随时间变化; (b)液滴形貌对比
Figure5. Comparison of experimental and simulated results: (a) Spreading factor changes with time; (b) comparison of droplet morphology evolution.

Chandra和Aziz^[23]设计了一种低温液滴发生装置, 利用该装置实验观测到直径约为1.9 mm的液氮撞击铜面后发生了Leidenfrost现象, 液滴撞击壁面的速度由下落高度控制, 实验发现液氮滴在距壁面27 mm处自由下落撞击铜面后在其底部形成气膜, 经8.2 s后液氮滴才蒸发完全, 而本文模拟的最大撞击时间尺度t_max = D₀/U₀ = 5 ms, 即液滴完全蒸发时间尺度远大于液滴撞击时间尺度; 另外, 实际液氮喷雾冷却过程中, 喷嘴距热沉面距离较短, 液氮滴初始速度较大, 液氮滴在撞击壁面之前不会因为环境温度过高而完全蒸发. 通过上述分析可看出本文的数值模拟方法及建立具有一定初速度和初始直径的液氮滴撞击水平壁面的物理模型是可行的.

2

-->

3.1.壁面润湿性和初始速度的影响

液滴撞击壁面的铺展行为是润湿现象的一种. 壁面的润湿性是影响表面沸腾和凝结性能的重要参数. 一般壁面润湿性可通过静态接触角来表征: 亲水壁面(θ < 90°)、疏水壁面(θ > 90°)和超疏水壁面(θ > 150°).
图6为0.5 mm的液滴以0.1和1.6 m/s分别撞击亲水壁面(θ = 30°)、普通壁面(θ = 90°)和疏水壁面(θ = 150°)的铺展行为随无量纲时间(t^* = tU₀/D₀)的变化趋势. 分析图6(a)可知, 当初速度为0.1 m/s, θ = 30°时, 撞击初期液滴底部与壁面直接接触, 接触区域发生剧烈相变, 从而在液滴底部形成空穴(t^* = 0.02), 与此同时液滴动能向表面能和黏性耗散转化, 三相接触线处的合力指向外侧, 三相接触线不断外延, 固液接触面积增大, 促使更多液相向气相转移, 空穴体积逐渐增大; 在t^* = 0.1时, 液膜的压迫使空穴内部蒸气沿径向逃逸, 液滴底部初步形成连续气膜. 亲水壁面具有较大的表面能有利于液滴撞壁后液膜边缘沿径向向外铺展, 而液体的表面张力则阻碍这种趋势的发生, 导致液膜颈部出现极薄液膜区, 最终伴随着换热过程的进行, 在颈部发生断裂(t^* = 0.16).
图 6 初始速度为(a) 0.1 m/s和(b) 1.6 m/s的液滴撞击不同润湿壁面的形貌演变(D₀ = 0.5 mm)
Figure6. Morphology evolution of droplets impinging on different wetted walls with initial velocity (D₀ = 0.5 mm): (a) 0.1 m/s; (b) 1.6 m/s.

在普通和疏水壁面上, 撞击前期液滴在壁面上的铺展状态以及连续气膜的生长与亲水壁面较为相似, 但由于接触角增大, 壁面对液体吸附力减弱, 液滴在壁面上的铺展速度相对减慢, 铺展面积相对减小, 液滴边缘难以发生破碎; 尤其是当液滴撞击疏水壁面时, 液体在其表面难以润湿, 黏性耗损减少, 使得液膜表面能较大, 从而表现出在t^* = 6.4时液膜开始回缩, 液体向液膜中心聚集.
液滴在回缩过程中表面能的变化可表示为

$\Delta E = \int_{{A_1}}^{{A_2}} {\sigma \cos \theta {\rm{d}}A},$

式中, A₁和A₂表示液滴回缩过程中的润湿面积, 可以看出壁面润湿性直接影响液滴形貌演变过程中表面能变化的大小, 若润湿面积固定, 随着接触角增大(0°—180°), 液滴收缩需要克服的表面能减小^[24], 因此液滴碰撞疏水壁面后更容易发生回缩, 导致液体向中心聚集.
分析图6(b)可知, 当速度为1.6 m/s时, 液滴在三种润湿性壁面上快速铺展, 并且在边缘处发生二次小液滴飞溅(t^* = 0.16), 原因是液滴速度增大后We随之增大, 惯性力起主导作用并能够完全克服液体的表面张力和黏性应力, 因此在边缘处产生破碎的小液滴, 体积较小的二次液滴受到周围气体的阻滞和液膜内部径向流动的作用致使其斜向上运动; 在t^* = 1.1时, 三种润湿壁面连续气膜均初步形成, 伴随液滴在三种壁面上的快速铺展, 液体在壁面上的覆盖面积逐渐增大; 在t^* = 3.2时, 铺展成极薄的液膜, 由于表面张力和惯性力的拮抗作用导致颈部液滴不断剥落, 此时二次液滴粒径较大, 周围气体的阻滞作用相对减弱, 剥落的小液滴受液膜径向流动作用而不断外移; 在t^* = 6.4时, 液滴破碎完全形成链状小液滴群. 重要的是, 与低速液滴相比, 液滴在三种润湿壁面上形貌演变的差异性逐渐减弱, 在壁面的铺展面积逐渐接近, 且润湿面积均大于同等条件下的低速液滴, 表现出明显的速度效应.
图7为0.5 mm的液滴以0.1和1.6 m/s撞击不同润湿壁面的热流分布, 壁面温度均为300 K. 图8对应速度为0.1 m/s的液滴撞壁面后的温度分布, 静态接触角θ = 30°.
图 7 初始速度为(a) 0.1 m/s和(b) 1.6 m/s的液滴撞击不同润湿壁面的热流分布
Figure7. Heat flux density distribution of droplets impinging on different wetted walls with initial velocity: (a) 0.1 m/s; (b) 1.6 m/s.

图 8 不同时刻的温度分布(θ = 30°, U₀ = 0.1 m/s, T_w = 300 K)
Figure8. Temperature distribution at different times (θ = 30°, U₀ = 0.1 m/s, T_w = 300 K).

结合图6(a)液滴的形貌演化和图7(a)热流分布图可以看出, 当初始速度为0.1 m/s时, 在撞击初期(t* = 0.01), 液滴底部局部与壁面有短暂接触, 壁面中心接触区域传热方式以热传导为主, 因此三种润湿壁面的热流分布集中在壁面中心, 疏水、普通和亲水壁面中心的热流峰分别为2 × 10⁶, 4 × 10⁶和5 × 10⁶ W/m², 热流密度随着润湿性增强而增大; 随后(t^* = 0.02—0.1), 液滴撞壁后在壁面上沿径向逐渐铺展并伴随着剧烈相变, 液滴底部空穴内部压力逐渐增大直至将液膜托起, 发生Leidenfrost效应, 热流密度急剧下降, 亲水壁面上的最高热流峰值仅为初期的10%. 由图8(a2)可以看出, 在t* = 0.1时气膜的存在显著增大了热阻, 壁面中心处A点温度梯度仅为370 K/mm, 远小于边缘处B点的2469 K/mm, 这对于喷雾冷却是非常不利的. 由图6(a)的形貌图可以看出, 壁面上铺展的液膜中心和四周凸起, 壁面各处热阻不同, 反映在图7(a2)中为热流分布不均匀. 当t* = 0.3时, 整体热流峰值显著下降, 最高峰值不足35 W/m², 分析图8(a3)发现, 在撞击后期气膜充分生长, 厚度显著增大, 壁面中心处A点与边缘处B点的温度梯度较为接近约为222 K/mm, 可见在撞击后期换热被严重削弱. 另外, 随着接触角的增大, 液滴收缩需要克服的表面能变小, 因此液滴碰撞疏水壁面后更容易发生回缩, 导致液体向中心聚集, 中心处换热加强, 从而呈现出热流峰宽逐渐减小且热流向壁面中心转移的趋势.
分析图7(a)和图7(b)发现, 当撞击初速度增至1.6 m/s时, 壁面热流密度在不同时期(t^* = 0.01—0.3)均大幅提升; 在撞击初始阶段(t^* = 0.01), 液滴在三种润湿壁面上的形貌演变较为相似, 壁面热流曲线几乎重合, 并且热流集中在液滴与壁面接触区, 且热流峰面积明显大于同期的低速液滴; 当t^* = 0.1时, 液滴的动能向持续表面能和黏性耗损转化, 液膜边缘沿径向铺展, 铺展面积增大, 热流分布区域也随之增大, 此时液膜边缘处较薄, 边缘处换热得到加强, 从而在两侧边缘处出现两个热流峰, 但此时液相主体仍集中在壁面中心区, 中心处最高热流峰仍维持在5 × 10⁶ W/m², 而同期低速液滴的最高热流峰仅为其1/10; 当t^* = 0.3时, 液滴底部已形成连续气膜, 此时亲水壁面的热流密度峰值比初期下降了近4 × 10⁶ W/m², 但仍高于此时低速液滴的最高热流峰, 而且三相接触线处热流峰值(1250 W/m²)明显高于液膜中心(75 W/m²), 热流由中心向边缘转移; 原因主要有两点, 一是在三相接触线附近壁面温度与液膜之间的温差大于液膜中心的温差, 二是在壁面中心区域气膜较厚热阻较大.
观察图6和图7的液滴形貌演变和热流密度分布可以看出, 提高撞击速度后液滴在不同润湿壁面上的运动行为和相变行为差异性逐渐减弱, 壁面润湿作用呈现出较强的速度效应. Kim等^[25]实验发现, 随着壁面静态接触角由63°降低至2°, 液滴的蒸发速率由0.0208 mg·s^–1增大至0.0483 mg·s^–1, 可见增强壁面的润湿性有利于液滴在壁面的蒸发. 李大树等^[24,26]同样发现, 壁面润湿作用随初始速度的增加而减弱, 呈现出较强的速度效应, 并将其归因于液滴的初始速度增大导致壁面接触角所产生的毛细力作用降低所致, 并且还发现增大液滴初始速度, 热流密度随之增大, 这与本文的模拟结果是一致的.
图9给出了两种初速度的液滴撞击亲水壁面后在0.05 ms时的压力分布图; 速度为1.6 m/s的液滴撞壁后在壁面中心A (26000 Pa)和边缘区B (2000 Pa)的压力远大于速度为0.1 m/s的液滴在中心A (14000 Pa)和边缘处B (1000 Pa)的压力, 可见高速液滴沿径向的压力梯度较大可以快速沿径向铺展, 导致换热面积相对增大, 热流密度因而随着撞击速度的增大而增大.
图 9 不同速度的液滴撞壁后在0.05 ms时刻的压力分布
Figure9. Pressure distribution at 0.05 ms after droplets impact wall with different velocities.

总体看来, 壁面润湿性越强换热效果越好, 撞击初速度增大带来的速度效应不仅提高了壁面的热流密度, 也缩小了不同润湿壁面上的相变行为的差异性; 另外, 两种速度的液滴撞击壁后的热流密度随时间推移均剧烈下降, 分析认为Leidenfrost效应是诱发上述现象的主要原因. Leidenfrost效应涉及的气膜和液膜都是动态的, 气膜通过液体的蒸发供给, 并且由于液体施加在其上方的压力而流动, 其特征是在液滴压在蒸气层之上, 形成液滴悬浮现象, 并且由于蒸气反向压力和液滴惯性力二者的拮抗作用固定了蒸气层的厚度, 可见气膜的厚度依赖于反向压力和惯性力的相对大小. 初始速度的改变直接影响惯性力的大小, 图10给出了两种初始速度下无量纲气膜厚度δ_t/D₀随无量纲时间t^*的变化曲线, 无量纲气膜厚度定义为壁面中心气膜厚度δ_t与液滴初始直径D₀的比值.
图 10 初始速度为(a) 0.1 m/s和(b) 1.6 m/s液滴撞击不同润湿壁面的无量纲气膜厚度随无量纲时间的变化
Figure10. Dimensional film thickness with dimensionless time curve of droplets impinging on different wetted walls with initial velocity: (a) 0.1 m/s; (b) 1.6 m/s.

结合图6(a)和图10(a)可以看出, 在气膜生长初期, 液滴惯性力相对较大, 气膜厚度较小, 随着铺展过程的进行, 液滴垂直向下的惯性力逐渐削弱, 同时部分液相逐渐向气相转化, 进入底部空穴, 气化反向压力对惯性力的比值快速增大, 反映在图10(a)中为t^* < 0.06阶段气膜厚度快速增长; 当t^* = 0.06时, 接触角对气膜厚度的影响逐渐显现, 随着接触角的减小气膜厚度逐渐增大, 亲水壁面上稳定的气膜厚度约为疏水壁面的1.3倍, 这与增强润湿性促进换热效果是符合的, 观察此时液滴的形貌可以发现亲水壁面上液膜铺展的面积明显大于普通壁面和疏水壁面, 使得亲水壁面上液相蒸发速率远高于普通壁面和疏水壁面; 当0.1 < t^* < 0.16时, 连续气膜已形成, 这意味着气化反向压力起主导作用, 液膜被托起, 在疏水和普通壁面上空穴内部分气相沿横向逃逸, 气膜厚度增长趋势减缓; 当t^* > 0.16时, 液滴的动能完全转化为液体的表面能和黏性耗损, 此时表面能达到最大并驱使液体向液膜中心聚集, 换热面积随之减小, 因此在普通和疏水壁面中心气膜厚度出现轻微下降; 但是此时亲水壁面上气膜厚度仍持续增加, 原因是亲水壁面上液相向气相迁移的质量较大, 并且在液膜边缘处生成破碎的小液滴, 液相主体惯性力下降显著, 有利于气化反向压力将液相主体托起.
分析图6(b)和图10(b)发现, 当初始速度增大至1.6 m/s时, 液滴动能大幅增加, 液滴在三种壁面上快速铺展, 固液接触面积相对增大, 气液相变剧烈, 因此反映在图10(b)中为气膜厚度短期内(t^* < 0.75)快速增大; 当t^* > 0.75时, 气膜厚度增长趋势变得缓慢, 随时间推移至t^* = 4.5时附近三者气膜厚度相当(无量纲气膜厚度约为0.23), 并在随后的过程中气膜厚度出现轻微波动, 这可能是因为在撞击后期, 换热被严重恶化, 液滴在壁面破碎成链状的液滴群, 导致气膜内部各处压力不均匀. 另外, 对比图10(a)和图10(b)发现高速液滴撞壁后气膜厚度整体上显著小于低速度的液滴, 原因是撞击速度增大后液体的惯性力远大于蒸气反向压力^[27]. Chandra和Aziz^[23]开展了液氮滴撞击壁面的实验研究, 基于动量守恒和能量守恒建立了气膜厚度与壁面/液体之间的温差$\Delta T = {T_{\rm{w}}} - {T_{\rm{l}}}$, 铺展直径D_t以及撞击初速度U₀的数学模型:

${\delta _{\rm{t}}} \sim {\left[ {{{3{{\left( {\frac{{{D_{\rm{t}}}}}{2}} \right)}^2}{\lambda _{\rm{v}}}{\eta _{\rm{v}}}\Delta T} / ({U_0^2{\rho _{\rm{v}}}{\rho _{\rm{l}}}L}}}) \right]^{{1 / 4}}},$

式中, ${\lambda _{\rm{v}}}$为气相导热系数, ${\eta _{\rm{v}}}$为气相黏度, L为液相汽化潜热.
由(24)式可以看出气膜厚度与液滴初速度成反比, 这很好地验证了本文中液滴速度增大气膜厚度减小的结果.
液滴初始速度增大有利于降低气膜厚度, 并且促使液滴顶部温度较低的液滴能够沿径向快速扩展, 最终导致热流密度大幅增加, 这对换热行为是非常有利的. 对于低温喷雾冷却而言, 提高热沉面的润湿性和液滴撞击速度对喷雾冷却有促进作用; 在膜沸腾区换热过程主要集中在撞击初期, 随着时间推移气膜的厚度逐渐增大, 换热性能剧烈恶化, 因此保持新鲜液滴以一定撞击速度持续冲击热沉面是影响换热效果的关键^[28].
2 -->

3.2.壁面温度的影响

壁面温度是影响液滴撞壁相变行为的关键参数, 尤其是伴随Leidenfrost效应的低温液滴. 研究发现, 液滴撞击不同温度壁面的运动特性差异较大, 可能出现接触沸腾、接触反弹、反弹时液滴部分破碎、液滴下部蒸气吹扫导致破裂、从液滴上表面喷射、液滴完全分解以及完全分解时伴随的破碎片段的快速移动等现象^[29]. 图11和图12分别给出了直径为0.5 mm的液滴撞击温度分别为300, 400, 500 K的壁面后的热流密度分布和局部温度分布图, 静态接触角为30°.
图 11 (a) 0.05 ms和(b) 1.6 ms时不同温度壁面上的热流密度分布
Figure11. Heat flux distribution on different temperature walls at (a) 0.05 ms and (b) 1.6 ms.

图 12 (a) 0.05 ms和(b) 1.6 ms不同温度壁面上的液滴温度分布
Figure12. Droplets temperature distribution on different temperature wall at (a) 0.05 ms and (b) 1.6 ms.

综合分析图11(a)和图12(a)的热流分布、形貌演变以及壁面温度分布可以看出, 在0.05 ms时液滴底部未形成连续气膜, 液滴底部局部液相仍与壁面接触(2.7和3.3 mm), 换热以热传导为主从而表现出较高的热流峰, 而在2.85, 3.00和3.15 mm处出现空穴, 温度梯度沿径向分布不均匀, 导致壁面热流沿径向由内向外呈现出现多个热流峰, 而且提高壁面温度显然增大了壁面和液体之间的温差, 热流密度随之显著提高; 由图11(b)可以看出, 在1.6 ms时, 液滴在三种温度的壁面上均破碎完全形成链状小液滴群, 大量液相蒸发成气相, 并且形成较厚的气膜, 此时换热性能急剧下降, 而且不同温度壁面上的热流分布曲线差异性也逐渐缩小; 另外观察图12(b)发现, 在破碎小液滴区仍维持一定的温度梯度, 而在断裂区的温度梯度则显著下降, 导致壁面出现多个热流峰. 吴苏晨等^[30]发现提高壁面和液滴的温差, 在撞击开始阶段平均热流密度出现一个较大的峰值, 此后热流密度变化缓慢, 分析认为在随后的阶段气膜托起液滴, 液滴无法直接润湿壁面, 导致热流密度下降, 换热效果变差, 但未观察到液滴破碎, 并且未详细分析气膜厚度随时间的变化规律.
图13给出了1.6 m/s液滴撞击温度分别为300, 400和500 K的润湿壁面上气膜厚度随无量纲时间的变化, 可以看出液滴撞击不同温度壁面后气膜厚度随着时间推移呈现相同的变化趋势, 在撞击初始阶段(t* < 0.75)是气膜厚度快速增长区, 而后气膜增长趋于平缓; 观察t* = 0.75时液滴形态发现, 在500 K壁面上, 液体和壁面的温差最大, 空穴内部压力随着温度的升高和液相蒸发供给而快速增大, 致使此时底部空穴内部压力大于液滴惯性力, 足以将液滴托举, 而此时300和400 K壁面上连续气膜未完全形成, 因而在500 K的壁面上气膜较厚; 当t* > 1.1时, 观察500 K壁面上液滴的形貌可以看出, 由于表面张力的作用液膜边缘处产生了破碎小液滴, 此时气膜内部沿横向逃逸的蒸气对小液滴的阻滞作用和液膜内部流体的径向运动, 致使其斜向上运动, 最终液膜边缘不断向外剥落小液滴, 在壁面形成链状小液滴群, 换热性能急剧恶化, 可见Leidenfrost效应的存在导致液滴与壁面之间的换热行为主要发生在撞击初期.
图 13 不同壁面温度上无量纲气膜厚度随无量纲时间的变化
Figure13. Curves of dimensionless film thickness with dimensionless time on different wall temperatures.

2 -->

3.3.气膜生长行为

通过上述液滴形貌演变、壁面热流密度分布和气膜厚度随时间变化规律的分析可以看出, 气膜增长是影响低温液滴换热的关键, 主要受液体与壁面传热控制, 类似于过热液体中气泡的生长, 而换热性能与气膜厚度成反比, 且传热过程主要集中在撞击初期. 因此下文着重分析撞击初期气膜生长行为. 为便于理论分析, 假设液滴撞壁后底部为扁平状的球体(图3), 液滴和气膜之间存在连续气膜, 其厚度δ_t远小于液滴初始直径D₀和接触直径D_t, 接触直径D_t与液滴初始直径D₀存在如下关系^[31]:

${D_{\rm{t}}} \sim {\left( {\frac{{{D_0}}}{2}} \right)^2}\kappa,\quad \frac{{{D_0}}}{2} < {\kappa ^{ - 1}},$

${\kappa ^{ - 1}} = \sqrt {\sigma/({\rho _{\rm{l}}}g)},$

式中${\kappa ^{ - 1}}$为毛细特征长度.
ΔT在300—500 K时, 辐射换热约占对流换热的5%, 因此不考虑辐射换热的影响^[32], 壁面主要通过气膜的一维导热向液滴传递热量, 导热热通量与壁面和液体之间的温差ΔT、铺展面积${\text{π}} {( {{{{D_{\rm{t}}}}}/{2}})^2}$、蒸气导热系数λ_v成正比, 与气膜厚度δ_t成反比, 引入蒸发潜热L可以推导出单位时间内液体蒸发质量${{{\rm{d}}m}}/{{{\rm{d}}t}}$, 依据能量守恒:

$[L+c_{\rm p}(T_{\rm sat}-T_{\rm l})]\frac{{{\rm{d}}m}}{{{\rm{d}}t}} \sim {\lambda _{\rm{v}}}\frac{{\Delta T}}{{{\delta _{\rm{t}}}}}{\text{π}}{\left( {\frac{{{D_t}}}{2}} \right)^2}.$

在撞击初期气膜内部的压力主要由撞击惯性力决定, 且内部压力仅在液滴发生形变的撞击初期才有意义. 因此在撞击初期, 液膜压迫下方气膜沿径向逃逸, 逃逸的气膜通过液体蒸发产生的蒸气供给, 此时壁面和液相界面可近似看成无滑移条件, 单位时间内气膜质量增长率与液滴底部蒸发速率应该保持一致^[32]:

$\frac{{{\rm{d}}m}}{{{\rm{d}}t}} \sim {\rho _{\rm{v}}}\frac{{{\rm{d}}{\delta _{\rm{t}}}}}{{{\rm{d}}t}}{\text{π}}{\left( {\frac{{{D_t}}}{2}} \right)^2}.$

联立方程(27)和(28)可以得出

${\delta _{\rm{t}}} \sim {\left\{ {{{{\lambda _{\rm{v}}}\Delta Tt}/ {[L+c_{\rm p}(T_{\rm sat}-T_{\rm l})]{\rho _{\rm{v}}}}}} \right\}^{1/2}}.$

观察(29)式可以看出在撞击初期气膜厚度与温差和时间的1/2次方成正比. Breitenbach等^[33]基于能量守恒分析在撞击前期的气膜增长行为, 得出液滴撞击高温壁面后的气膜厚度与液滴物性、温差以及撞击时间的关系:

${\delta _{\rm{t}}} \sim \left[ {K{e_{\rm{w}}}\left( {{T_{\rm{w}}} - {T_{{\rm{sat}}}}} \right)/(L{\rho _{\rm{l}}}}) \right]{t^{{1 / 2}}},$

$K = \sqrt {{{(B - G)}^2} + \frac{{4G}}{{\text{π}}}} - B - G,$

$G = \frac{{\sqrt {\text{π}} {\lambda _{\rm{v}}}{\rho _{\rm{l}}}L}}{{2\left( {{T_{\rm{w}}} - {T_{{\rm{sat}}}}} \right)e_{\rm{w}}^2}},$

$B = \frac{{\sqrt 5 \Delta T{e_{\rm{l}}}}}{{\sqrt {\text{π}} \left( {{T_{\rm{w}}} - {T_{{\rm{sat}}}}} \right){e_{\rm{w}}}}},$

式中, K为无量纲系数, e_w为壁面蓄热系数, e_l为液相蓄热系数. Breitenbach等^[33]还发现其模型的预测结果与Tran等^[7]的实验结果非常吻合.
由(29)和(30)式可以看出Breitenbach等^[33]的理论模型与本文分析结果在表达形式上保持较高的一致性; 图14(a)拟合了在撞击初期300—500 K壁面上气膜厚度随时间的变化(0.01 < t* < 1.05), 这与图14(b) Breitenbach等^[33]的变化曲线非常相似, 可见在传热控制阶段气膜的厚度不仅与壁面/液体之间温差正相关, 还与时间的1/2次方呈正相关; 而在撞击后期, 主要换热过程基本结束, 换热性能恶化严重, 气膜内部压力的作用减小, 气膜的内部应力起主要作用, 并在最后阶段控制气膜的微小波动, 从而表现出图10和图13中撞击后期阶段的现象, 气膜厚度随时间变化曲线较为平缓.
图 14 不同壁面温度上气膜厚度随时间的变化　(a)拟合结果; (b) Breitenbach等^[33]的分析结果
Figure14. Curves of film thickness with time on different wall temperatures: (a) Fitting results; (b) Breitenbach et al.^[33] analysis results.

采用复合Level Set-VOF相界面追踪方法建立单液氮滴撞击不同润湿壁面的数值模型, 并结合实验数据验证了模型的可靠性, 为低温液氮喷雾冷却提供了数值研究基础, 得出如下结论:
1) 增强壁面润湿性(30°—150°)、提高液滴的撞击速度(0.1和1.6 m/s)以及提高壁面温度(300—500 K)均有利于换热性能的提升, 且不同润湿壁面上相变行为的差异性因撞击速度的增大而逐渐减弱, 表现出明显的速度效应, 换热过程主要集中在撞击初期, 随着时间推移, 换热性能急剧恶化;
2) 基于质量守恒和能量守恒建立了气膜生长的数学模型, 气膜厚度与壁面温度和撞击时间成正比, 模型预测结果很好地验证了数值分析结果, 并且与他人实验和理论分析结果非常吻合;
3) 保持新鲜液滴以一定的速度持续冲击热沉面是影响换热效果的关键, 在膜沸腾阶段, 通过强化壁面条件和提高撞击初速度来促进换热在撞击初始阶段仍可获得较为显著的效果, 这对于指导低温喷雾冷却有重要的参考意义;
实际低温液氮喷雾冷却的过程涉及多种复杂物理场, 耦合具有三维性, 而本文通过建立二维模型研究球形单液滴撞壁后的相变行为, 存在一定的误差与局限. 今后的相关研究工作将结合相关实验建立精度更高的三维模型开展相关撞击过程的数值研究.