表面张力驱动对流对液滴散热影响的机理研究

The Influence of Surface Tension Driven Convection on Droplet Heat Dissipation

摘要: 为了揭示表面张力驱动对流对液滴散热影响的作用机理，本文采取数值模拟方法对热毛细对流(表面张力驱动对流)下的液滴散热进行研究。建立了以硅油作为工质的二维液滴模型，考虑在常重力和微重力两种情况下，对于液滴内部流动与传热特性的影响。结果表明，温差越大液滴内部流速越快，温度梯度变化越剧烈，其热毛细对流的换热效果越明显；热毛细对流使得液滴背部流速较大区域扩大，并且其内部的低速流动区域减少甚至消失，热毛细效应增强液滴换热效果得到明显提升；重力水平的增大，液滴的整体散热性能也得到增强。

Abstract: In order to reveal the mechanism of the effect of surface tension driven convection on droplet heat dissipation, this paper adopts numerical simulation methods to study droplet heat dissipation under thermocapillary convection (surface tension driven convection). A two-dimensional droplet model using silicone oil as the working fluid was established, considering the effects of constant gravity and microgravity on the internal flow and heat transfer characteristics of the droplet. The results indicate that the larger the temperature difference, the faster the internal flow velocity of the droplet, and the temperature gradient changes. The more intense the heat, the more obvious the heat transfer effect of capillary convection; The thermocapillary convection expands the area with high velocity on the back of the droplet, and reduces or even disappears the low-speed flow area inside. The thermocapillary effect enhances the heat transfer efficiency of the droplet and is significantly improved; The increase in gravity level enhances the overall heat dissipation performance of droplets.

1. 引言

近年来，关于热Marangoni效应的研究日渐成为热点，该效应驱动的流体流动也被称为热毛细对流。这一流动广泛存在于晶体生长、镀膜、焊接及空间微重力环境，其引起的对流换热及流动不稳定性引起了研究者的广泛关注。段文豪等[1]利用激光光热效应对热毛细对流进行主动控制，通过数值模拟的方法研究光热效应作用下矩形液池内热毛细对流热质传输特性，揭示激光光热效应对稳态热毛细对流的作用机理。陈上通等[2]研究了微重力下液体在椭圆截面管道内毛细驱动流动行为。徐泽来等[3]基于扩散界面方法，对追踪可混溶相的对流扩散方程进行改进，发展了可以求解带有Marangoni效应、界面运动、可混溶相的对流和扩散互相耦合的多相流动数值方法并进行了初步的研究。Xu [4]发现液滴表面应该存在双向流动，避免在其表面出现流动失衡，因此必然存在一个流动理论静止点，也可称为流动滞止点，在该点两侧的流动方向相反。金哲岩[5]发现在加热固体面，可能由于Marangoni效应的存在，液滴散热蒸发时内部存在涡流，表明内部是具有层流的某些特征。徐巧变等[6]在光电化学转化过程中气泡界面处建立Marangoni效应模型，阐明气泡对光电极表面Marangoni效应的作用机理，探Marangoni效应对气泡生离过程的影响。高鹏等[7]研究了液滴或气泡在不均匀表面张力的驱动下，发生热毛细迁移现象的问题。此外，Maxwell [8]指出，在空气环境中，研究了球形液滴蒸发速率与球半径R的关系，实验证明两者存在成正比函数变化。Picknet等[9]则分析了在同一散热通量条件下，液滴接触角和其半径大小的关系，设定两种蒸发散热模式：即固定接触角模式和固定接触线模式，前者接触半径减小而接触角不变，后者则接触角逐渐减小而接触半径不变，经过数据耦合，表明实际蒸发过程是其两者的混合模式。Lopes等[10]实验探究了接触基底热物性对液滴散热蒸发影响。Saada等[11]数值研究了受热基底热物性对液滴蒸发时，其表面温度、液相浓度以及流动速度分布的影响，这些因素均会影响Marangoni效应的强弱，改变热毛细对流的传热强度，从而影响液滴整体散热效果Deegan等[12]对含有固体微粒的液滴进行蒸干试验后，发现形成了咖啡环状沉积，认为该种现象的产生是液滴内部Marangoni对流效应的结果。Hu等[13]发现接触角的大小会影响液滴内部热毛细对流的方向。Girard等[14]研究了液滴加热半径和润湿半径的关系，发现当前者大于后者时，液滴的三相接触线温度明显高于其球形顶部温度；反之，则接触线呈现较高温度。Ghasemi等[15]发现在液滴不论是向内的受热内对流，还是向外的能量散失，即散热蒸发时的能量交换过程，热毛细对流在其中起着决定性作用，特别是液滴基底边缘固液气三相接触线附近，这种效应表现得尤为明显。了解到已有研究其中大部分是针对固体基底液滴薄液膜区间的蒸发传热特性的研究、对液滴温差加热时，液滴传热失稳的现象探究、液滴撞击受热界面时的换热性能以及液滴对流换热能力和其自身热流密度分布的情况分析等。

虽然以上研究已经取得了显著的成果，然而针对表面张力驱动对流对液滴散热影响的机理性研究，相对比较缺乏。为此本文在温差加热状态下，通过对表面张力驱动对流，在不同基底温度、常重力条件以及区别量级微重力条件下，通过研究Marangoni对流效应中度场，甚至浓度场的分布特性，将有助于提高焊接工艺、晶体生长和涂层技术等领域的效率。

2. 物理数学模型

图1为物理模型，液滴底部半径为r，液滴基底温度为T_w (即半球底部温度)。液滴周围环境温度设定为25℃时的环境温度，即T_a = 298.15 K。液滴表面为自然对流换热，其对外界的辐射散热忽略不计。液滴的接触角α为90˚，即液滴为半球形状，半径为10 mm。为了简化计算模型，我们假定：1) 液滴内部流体流速较低，流动为层流；2) 流体视为单相不可压缩牛顿流体；3) 忽略液滴蒸发效应，其形状不随时间发生变化，且T_w > T_a；4) 加热基板温度为液滴下边界温度。在数值模拟中，忽略硅油材料的黏性耗散导致的发热(称为黏性加热)，并假设流体为牛顿流体，以便于在流场求解过程中，可以通过质量、动量和总能量的平衡来表征，因此这种平衡所用的三维控制方程如下：

Figure 1. Physical model

图1. 物理模型

连续性方程：

$\nabla \cdot u=0$ (1)

纳维–斯托克斯方程：

$\rho \frac{\partial u}{\partial t}+\rho u\cdot \nabla \left(u\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla u+\nabla u^{\text{T}}\right)\right]$ (2)

总能量方程:

$\rho C_p\frac{\partial T}{\partial t}+\rho u{C}_{p}u\nabla T=\nabla \cdot \left(k\nabla T\right)$ (3)

其中，u是流体速度，p是流体压力，ρ是流体密度，μ是流体黏度。

2.1. Marangoni效应

Marangoni效应的强弱尺度，最主要的影响因素是表面张力梯度的变化尺度，其中对于表面张力的变化主要有热效应和成分(浓度)效应：

$\text{d}\sigma =\frac{\partial \sigma }{\partial T}\text{d}T+\frac{\partial \sigma }{\partial c}\text{d}c$ (4)

对于本文的研究不涉及浓度效应，所以主要研究热效应的影响：

$\gamma =\left|\frac{\text{d}\sigma }{\text{d}T}\right|$ (5)

通过用表面张力随温度变化的绝对值，表示表面张力对温度的依赖性。热诱导的马兰戈尼效应将剪切应力的法向分量与温度的切向导数相关联。对于在层流中没有压力贡献的不可压缩流体，可以写为：

$\left[\mu \left(\nabla u\right)+{\left(\nabla u\right)}^{\text{T}}\right]\cdot n=\gamma {\nabla }_{t}T$ (6)

这是作用在通过纳维–斯托克斯方程建模的流体自由表面(通常为气–液界面)上的边界条件。

Marangoni数是一个无量纲数，给出热毛细效应与黏性力的比值。

$Ma=\frac{\left|\gamma \right|L\Delta T}{\mu k}$ (7)

其中L是系统的长度尺度(单位为m)，是整个系统的最大温差(单位为K)，是动力黏度(单位为kg/s/m)，是热扩散系数(单位为m²/s)。

2.2. 边界条件

上边界条件：Marangoni效应在界面(液体/空气)上诱导出的力的公式：

$\eta \frac{\partial u}{\partial y}=\gamma \frac{\partial T}{\partial x}$ (8)

这里是表面张力的温度导数(N/(mK))。该公式指出，液体表面上的表面张力与温度梯度是成正比的。

下边界条件：$u=v=0$ ，$T=T_w$ 液滴初始温度为$T_a$ (9)

2.3. 硅油材料的物性参数

多种具有不同特性的硅油材料的物性参数如表1所示。

Table 1. Physical properties of silicone oil

表1. 硅油材料的物性参数

物性参数	密度$\rho$	表面张力系数$\sigma$	表面张力温度系数${\sigma }_T$	动力粘度$\mu$
数值	950 kg/m^3	0.0159 N/m	−6.23 × 10^−^5 N/(m∙K)	0.019 Pa/s
物性参数	热膨胀系数$\beta$	热导率$\lambda$	恒压热容$C_{pl}$
数值	0.00107 K^−^1	0.14 W/(m∙K)	1.63 × 10^3 J/(Kg∙K)
物性参数	比热率	普朗特数$P_r$
数值	1	10.3

3. 计算方法

本文数值模拟计算基于Navier-Stokes方程开展，在该模拟过程中，使用Boussinesq近似。在COMSOL中创建仿真模型，在其内部软件界面中可以提供“模型开发器”并且可以支持灵活地耦合多个物理场，对于解决一些实际工程问题和准确模拟真实场景下的物理现象。所以可以通过完成设置硅油工质和基本传热流动所需参数，以及流动状态、传热过程与多物理场的设定，在此基础上对模型设置合理的网格划分，最后进行稳态研究。得到相关的液滴表面的速度变化图、压力等值线图、表面温度差异图和整个液滴等温线图等模拟结果示意图。在基底温度为300 K、环境温度为298.15 K的温差加热，且不存在重力因素的影响的条件下。在材料库中，查找到模拟所需的硅油材料，添加到列表中并定义参数。

4. 计算结果

4.1. 温差对液滴换热特性的影响

图2是基于模型的硅油流体的二维液滴的表面速度大小和速度场的分布，液滴初始温度298.15 K，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。由流线上箭头的方向可知，液滴表面及内部是双向流动，明显地内部存在对称的流动涡流。在常重力条件温差为1.85 K时，其液滴内部存在两组流向相反且对称的涡胞。液滴上边沿明显存在理论滞止点。其两侧流体流动方向相反。液滴上边沿涡流流动方向相反两侧边界处的流动速度较大，而滞止点附近、液滴底部以及对称涡流中心流动速度较小。随着与基底温差的不断增大，以上流体速度差异的变化不断减小，而且涡流流速增加且有逐渐消失的趋势。基底温度为305 K时为两个涡流，且流线较为密集；温差为11.85 K时也为两个对称涡流，但其内部流线明显减少。温差加热下，可以观察到液滴的内部流动存在明显的区域差别，可大致分为对流区、滞止区、边缘区等。其中表面速度变化明显的区域逐渐增加，温差的不断增大，其液滴内部的浮力对流强度和液滴表面的热毛细对流强度都有明显增加。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 2. Flow field and velocity distribution under different temperature differences under constant gravity

图2. 常重力下不同温差下流场和速度分布

图3给出了基于模型的硅油流体的二维液滴的表面温度场和温度场的分布，液滴初始温度为298.15 K，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。液滴上边沿的滞止点附近区域温度最低，液滴基底部分、以及触角上部带状区域的温度较高。在对称的涡流中心处温度等值线变化明显，温度梯度差异大。温差为1.85 K时，滞止点两侧及垂直方向上，温度呈等梯度变化；而温差达到6.85 K时，等值线弯曲变化程度增加显著。随着加热温差的不断增大，这种变化趋势更加明显，其液滴内部温度等值线变化差异大；液滴内涡流交界线上，温度等值线密集程度存在显著差异，温差越大，等值线越稀疏，表明该处温度变化越剧烈。说明此时，液滴内部流体流动状态较为活跃，表面浮力对流和热毛细对流对液滴传热流动效果的提升显著。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 3. Temperature and contour distribution under different temperature differences under constant gravity

图3. 常重力下不同温差下温度及等值线分布

图4给出了基于模型的硅油流体的二维液滴的上边沿边界的温度梯度变化曲线，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。液滴上边沿温度梯度变化呈对称分布，以温差为6.85 K时为例，其温度梯度变化区域可分为三部分，区域Ⅰ：从基底左右各上移8˚的对称区，温度梯度急剧增加。区域Ⅱ：8˚~80˚之间的对称区，温度梯度只是略微增大。区域Ⅲ：80˚~90˚之间的对称区，温度梯度存在先下降后急剧增加的特殊现象，但最终的温度梯度变化数值仍低于区域Ⅱ。随着基底温度的不断增加，可明显观察发现，区域Ⅰ的面积和区域Ⅲ的面积不断减少，区域Ⅱ所占面积不断增加。整个液滴的温度趋向于稳定，表现出均一性。其中三条曲线随各自温度场变化趋势基本一致，然而由于基底温度的不同，其温度梯度的数值差异巨大。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 4. Changes in temperature gradient at the upper boundary under different temperature differences under constant gravity

图4. 常重力下不同温差上边界温度梯度变化

图5给出了基于模型的硅油流体的二维液滴的下边界的传导热通量的变化曲线，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。下边界热通量的变化以底部中心为呈对称分布，液滴触角外侧附近与环境存在换热，因而无法确定其数值变化。以温差为11.85 K为例，从图中明显地发现存在三个区域，区域Ⅰ：下边界弧长0.008~0.012 m之间，传导热通量增加速率迅速提高，随着接近底部中心，几乎接近垂直状态。区域Ⅱ：下边界弧长0.002~0.008 m与0.012~0.018 m的两部分，同时在增加，不过明显低于前一部分，区域Ⅲ：0~0.002 m与0.018~0.02 m的部分，传导热通量依旧在增加，但速度较缓慢。而温差为1.85 K时的变化曲线可以近似的看作两个图5(b)连接在一起，因为其液滴内部存在两组对称的涡流，其他变化特征类似。在液滴基底中心的传导热通量最大，在不断靠近中心的边界上，热通量也随之增加。随着基底温度的增加，以上三个区域的分布特点更加明显，而且其中心点的最大热通量也不断增大。

4.2. 相同温差下不同对流对换热特性的影响

常重力条件下，模拟一个温差作用下的两种情况：(1) 单纯的常重力驱动浮力对流，(2) 常重力驱动浮力对流 + 热毛细对流，通过比较两者的传热对流特点，分析热毛细对流对液滴换热的影响。给到液滴下边界温度为305 K，通过进行模拟研究计算，得到如下结果。其中环境温度298.15 K，并且液滴初始温度等于环境温度。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 5. Lower boundary heat flux at different temperatures under constant gravity

图5. 常重力下不同温度时下边界热通量

图6给出了基于模型的硅油流体的二维液滴的表面速度大小和速度分布，温差为6.85 K。重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。液滴内部均存在对称且流向相反的涡流场，表明不同对流形式下，液滴内部始终存在对称涡流，使得液滴自身散热达到最大程度。纯浮力对流下，液滴内部流体高速区较多，尤其在液滴中间处存在明显的较高速区，且液滴表面的较高速度区域较小。浮力对流+热毛细对流时，液滴上表面边沿的流体流动高速区明显增加，其内部较高速区域显著减少，液滴内部流动速度普遍均一。并且在比较了液滴区域的流速情况后，发现存在热毛细对流时，液滴内部和边沿的流速，无论是整体的情况，还是速度最值均有明显的增加。

图7给出了基于模型的硅油流体的二维液滴的表面温度大小及分布，温差为6.85 K。重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。在常重力条件下，液滴内部只存在纯浮力对流时，高温区在液滴垂直中线上及底边界温度集中，温度等值线分布稀疏，说明温度梯度变化显著，而且低温点则向上边沿两侧分离；液滴内部对流形式是浮力对流 + 热毛细对流时，高温区只分布在基底上方，低温点处于上边沿中心，与此同时，在中线上的温度等值线稀疏。

图8给出了基于模型的硅油流体的二维液滴的下边界底部热通量变化曲线，温差为6.85 K。重力加速度为−9.8 m/s²，液滴半径为0.01 m的情况。在常重力条件下，不同对流形式下的下边界热通量的变化存在明显区别，处于纯浮力对流时，下边界底部热通量呈现由两接触角向中心递减的变化特征，且在中心处达到最低值。在这个过程中，呈现出向底部中心靠近，其热通量随胡长的变化率先增大后减小，最后又增大的特点。而在浮力对流和热毛细对流共同作用下，下边界热通量则从接触角向中心处递增的变化特点。两者在中心处的热通量变化存在显著差异，这表明在基底温度为305 K时，关于液滴内的流通传热特性由热毛细对流主导，所以在模拟过程中热通量的变化表现出与纯浮力对流时的趋势不同。

(a) 纯浮力对流 (b) 浮力对流 + 热毛细对流

Figure 6. Distribution of different flow velocities at constant gravity with ∆T = 6.85 K

图6. 常重力下∆T = 6.85 K时不同流动速度分布

(a) 纯浮力对流 (b) 浮力对流 + 热毛细对流

Figure 7. Different flow temperature sizes and distributions at ∆T = 6.85 K under constant gravity

图7. 常重力下∆T = 6.85 K时不同流动温度大小及分布

(a) 纯浮力对流 (b) 浮力对流 + 热毛细对流

Figure 8. Boundary heat flux under different flows at ∆T = 6.85 K under constant gravity

图8. 常重力下∆T = 6.85 K时不同流动下边界热通量

4.3. 不同温差，两种对流的换热特性

微重力条件下，当底面分别为300 K、305 K、310 K温差加热过程时，与液滴初始温度存在差异。通过进行模拟研究计算，得到如下结果。其中环境温度298.15 K，并且液滴初始温度等于环境温度。

图9给出了基于模型的硅油流体的二维液滴的表面速度大小和速度场的分布，液滴初始温度298.15 K，温差分别为1.85 K、6.85 K、11.85 K，重力加速度−0.0098 m/s²，液滴半径为0.01 m的情况。在微重力条件下，液滴内部存在流向相反的两个对称涡流(其中基底温度300 K时，涡流由四个变为两个)。上述的流线形状和数量并无明显区别。温差加热下，可以观察到液滴的内部流动存在明显的区域差异，其中随着温差的增大，液滴内部流速较高区域逐渐减少，表面边沿出现的高速区增加，且整体流速和流速最值也有显著提高，说明此时热毛细对流强度逐渐增大，加强了液滴内部流动和散热效果。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 9. Flow field and velocity distribution under different temperature differences in microgravity

图9. 微重力下不同温差下流场和速度分布

图10给出了基于模型的硅油流体的二维液滴的温度场和等温线分布，液滴初始温度为298.15 K，温差分别为1.85 K、6.85 K、11.85 K，重力加速度−0.0098 m/s²，液滴半径为0.01 m的情况。在微重力条件下，观察液滴中心线上的等温线的变化形状和密集程度，分析可知，随着基底温度的不断增加，液滴顶部中心处的低温区域不断减小；且在液滴中心线上温度等值线的密集程度不断减小，这也表明此处的温度变化剧烈。随基底温度的提高，液滴整体表现出温度不断增加，高温区的面积逐渐扩大，趋于温度同一性。

图11给出了基于模型的硅油流体的二维液滴的上边沿边界的温度梯度变化曲线，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−0.0098 m/s²，液滴半径为0.01 m的情况。在微重力条件下，液滴上边沿温度梯度变化呈对称分布，随着基底温度的提高，区域Ⅰ温度梯度急剧增加，且所占的面积急剧减小；区域Ⅱ对称区内温度梯度率不断减小，最后趋于不变。区域Ⅲ的温度梯度存在先下降后急剧增加的特殊现象，但最终的温度梯度变化仍低于区域Ⅱ。微小的温度差，引起液滴上边沿温度梯度的改变是巨大的，当温度达到一定程度时，整个液滴温度保持一致，区域Ⅱ占整个温度梯度变化区域面积的接近100%。直至液滴顶部低温区消失。

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 10. Temperature and contour distribution under different temperature differences in microgravity

图10. 微重力下不同温差下温度及等值线分布

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 11. Variation of upper boundary temperature gradient at different temperature differences under microgravity

图11. 微重力下不同温差时上边界温度梯度变化

图12给出了基于模型的硅油流体的二维液滴的下边界的传导热通量的变化曲线，温差分别为1.85 K、6.85 K、11.85 K，重力加速度为−0.0098 m/s²，液滴半径为0.01 m的情况。下边界热通量的变化以底部中心为呈对称分布，液滴触角外侧附近与环境存在换热，因而无法确定其数值变化。随着基底温度的不断提高，区域Ⅰ的传导热通量增加速率迅速提高，区域Ⅱ的面积虽然也在增加，但相较于前一区域，热通量的增加明

(a) ∆T = 1.85 K (b) ∆T = 6.85 K

(c) ∆T = 11.85 K

Figure 12. Lower boundary heat flux at different temperatures under microgravity

图12. 微重力下不同温度时下边界热通量

显缓慢。区域Ⅲ传导热通量增加速度更缓慢，甚至出现在温差较大的时，液滴下边界中心附近传导热通量反而下降的现象。

5. 结论

本文数值研究了不同参数条件下，表面张力驱动对流对液滴内部流动与传热特性的影响。得出如下结论：

1) 常重力条件下，液滴内部存在流向相反的对称涡流，流动过程中具有明显的流速分区现象；随着温差的不断增大，这种流体速度差异的变化不断减小，而且涡流有逐渐消失的趋势。

2) 在浮力对流和热毛细对流共同作用下，下边界热量则从接触角向中心处递增的变化特点，两者在中心处的热通量变化存在显著差异，液滴内的流通传热特性由热毛细对流主导，所以热量的传递表现出与纯浮力对流时的趋势不同。

3) 常重力条件下，当液滴内部只存在纯浮力对流时，高温区在液滴垂直中线上及底边界温度集中，温度梯度较小；液滴内部对流形式是浮力对流 + 热毛细对流时，高温区只分布在基底上方，低温点处于上边沿中心。

NOTES

^*通讯作者。