不混溶液体表面上蒸发液滴的动力学特性

李春曦，庄立宇，施智贤，叶学民

（华北电力大学动力工程系，河北保定071003）

引 言

近年来，一种液滴在另一种不混溶液体表面上的铺展受到了广泛关注，并在各个领域得到了应用，如石油回收、支气管药物的雾化输送、灌注液体型超滑表面(slippery liquid-infused porous surfaces，SLIPS)的相变换热等[1-5]。了解蒸发液滴在不混溶液体表面上的动力学特性对上述应用具有重要意义。

铺展系数通常用于分析液滴的润湿行为和几何结构。Harkins[3]提出了一种液体在另一种不混溶液体基底上的铺展系数，定义如下:

式中，σLi为纯液体i的表面张力；σLij为液体i和j 之间的界面张力，1、2 分别表示液体基底和蒸发液滴。当液滴放置于另一种液体上时，是否会自发铺展形成单分子薄膜还是会去润湿形成稳定的液滴取决于铺展系数S 的正负。正铺展系数(S ＞0)表明有足够的力驱动液滴自发铺展并润湿液体基底[4]，而负铺展系数(S ＜0)则会导致一种液体在另一种液体表面形成稳定形状的液滴。图1(a)显示了部分润湿(partial wetting)的情形下，上层液滴的接触角不为零且稳定漂浮在下层液体上，称为液体透镜；图1(b)为假部分润湿(pseudo-partial wetting)状态下，上层液滴形成了一个接触角为有限值的球形帽(有学者将这种状态也称为液体透镜[5])，其边缘在底层液体表面铺展形成了厚度为分子水平的薄膜；图1(c)显示的是在完全润湿(total wetting)的系统中，上层液滴表面平坦且完全覆盖下层液体，类似于一个扁平的“煎饼”。Sebilleau[5]指出还存在一种准完全润湿状态(pseudo-total wetting)，其图形与图1(a)近似，不同之处在于液体2不仅形成稳定形状的液滴漂浮于液体1之上，而且其表面完全被液体1覆盖。

蒸发对于不混溶液体表面上液滴的运动特性有重要影响，这方面研究以实验居多。Nosoko 等[6]研究了施加温度梯度的去离子水表面上挥发性正戊烷透镜的蒸发特性，发现透镜纯度和温差对蒸发时间和透镜几何结构有显著影响。Rahman等[7]研究了黏性液滴在液-液界面上铺展早期的动力学，发现由于液滴底部出现变形，导致铺展半径与时间的关系不再符合指数律。Shabani等[8]观察到油气界面漂浮水滴的两种稳定结构：非凝聚液滴和帽/珠液滴。Sun等[9]、Liu等[10]对不同挥发性物质的液滴在水表面的铺展行为进行了实验研究，他们均发现液滴的形状变化经历了铺展主导阶段(半径增大)和蒸发(半径收缩)主导阶段。Dussaud等[11]研究了挥发性油透镜在密度较大液体基底上自发铺展的动力学特性，指出挥发性的、不混溶液滴的铺展半径R与时间t 的关系为R - t1/2，其铺展指数1/2 小于非挥发性液滴的3/4。Buffone 等[12]对蒸发或非蒸发液体基底上蒸发液体透镜的形成、稳定性和温度分布进行了实验研究，揭示了透镜与液体基底间的相互作用。Shimizu 等[13]研究了压力对正戊烷和R113 液滴在水中蒸发的影响，认为蒸发液滴由气相和未蒸发的液相组成，指出正戊烷和R113气泡在上升过程中会出现形状振荡和轨迹摆动。

近年来，也有很多学者对不混溶液体表面上液滴的铺展过程展开了理论研究，但现有的理论模型较少涉及蒸发。对液体透镜的其他研究主要集中在铺展/去润湿动力学特性[14-16]、线张力[17]、润湿转变[18-19]和可变形基底上的润湿行为等[20]。Burton等[21]根据Young-Laplace 方程开发了一个分析透镜几何形态的代码，通过求解界面接触线上的Young-Laplace 微分方程，计算出三个界面的轴对称轮廓。Craster 等[1]运用润滑理论建立了透镜在各种不同参数组合下(包括两种液体的密度比、黏度比和表面张力比)的铺展，计算了拉伸应力在透镜中的作用，同时也为后续的研究奠定了基础，但其所建模型未涉及蒸发。Greco 等[22]研究了悬浮在两种流体界面上的液滴在平行于界面的温度梯度作用下的热毛细驱动运动，指出界面液滴的迁移方向与温度梯度相反。George 等[23]研究了在不满足Neumann 三角的系统中较重液滴在较轻液体上的漂浮，并将这种现象归因于线张力的作用，它可防止液滴完全浸没于液体基底。Phan等[24]证明了较重的液体透镜可以漂浮在较轻的液体基底上。Guan 等[25]发现在SLIPS 表面上，小固着液滴(直径约为2 mm)的蒸发遵循理想的恒定接触角模式。Brochard-Wyart 等[26]研究了黏性液滴在非黏性液体上的铺展情况，并确定液滴铺展遵循R - t1/4。Gelderblom 等[27]研究了液滴在液体基底上蒸发的情况，对两种液体的速度场进行了参数化研究，揭示了两种液体之间的接触角及黏度比对速度场的影响。另外，在微观尺度，采用分子动力学模拟的方法对此类问题进行的研究成果也有报道[28]。

图1 三流体系统可能的润湿情况Fig.1 Three possible wetting conditions for a three fluid system

综上所述，目前关于液滴在不混溶液体表面铺展的研究大多偏重于实验方向，但其结果受实验条件所限，存在许多不一致的地方，也不便于深入了解相关机理，而现有的理论模型较少涉及蒸发。为了更好地理解液滴在液体表面的蒸发机理，以促进相关技术的开发和利用，有必要采用理论建模和数值计算的方法对此类现象进行全面的研究。本文采用润滑理论建立蒸发液滴铺展的理论模型，通过数值计算，以接触线的运动等为着眼点，对不混溶液体表面上蒸发液滴的动力学特性进行探究，以增进对其物理机制的了解。

1 理论模型

1.1 物理模型

图2 液滴放置于均匀加热的液体基底上的示意图Fig.2 Sketch of a drop resting over uniform heated liquid substrate

液体i(i=1,2)铺展和蒸发过程满足连续方程、动量方程和能量方程：

式中，T*w为固体壁面温度。

在液-液界面处，满足速度场、温度场连续性，切向和法向应力平衡，运动边界条件以及热流量连续性：

润湿性液滴会铺展成接近分子厚度的薄膜，而非润湿性液滴则会形成稳定的透镜形态。本文主要研究非润湿挥发性的纯液滴在不混溶液体表面的铺展蒸发过程。当界面温度不均匀时，液膜的表面张力发生变化，产生Marangoni效应并导致液膜产生流动和明显变形。上层液滴位于h1≤z ≤h2和-∞≤x ≤+∞约束的区域，为了避免三相接触线处的应力集中，仿照Craster 等[1]的做法，假定上层液滴可在下层液体表面铺展成一层极薄的前驱膜，即图1(b)的情形，此时式(1)所描述的三种界面张力模型可简化为两种界面张力模型，即液1 与液2 以及液2与空气的界面张力。该前驱层厚度h∞由分离压维持，分离压Π表示为[1]：

式中，h ≡h2- h1；h∞为前驱膜厚度；n ＞m ＞1，指数(n,m)取(4,3)。式(26)中右侧第一项和第二项分别代表分子间斥力和引力对分离压的不同贡献。

1.2 初始条件和边界条件

采用PDECOL程序模拟上述演化方程。由于演化方程中的变量及物性参数均为无量纲量，因此除非特别说明，下文给出的表达式及变量分布曲线图均为无量纲形式。选取-5 ≤x ≤5 的一维区域进行计算，将该区域划分为均匀网格，并进行网格无关性验证，在满足精度和减少计算时间的条件下设定网格数为2000个。

计算初始条件为[1]：

图3 初始时刻界面高度分布Fig.3 Initial conditions of the interface height distribution

计算中PDECOL程序的边界条件设置如下：

在加热条件下，以演化方程式(23)和式(24)为基础，计算存在蒸发情形时的上层液滴在下层液体上的演化过程。结合相关资料[9,11]可得有量纲基本参数的取值范围(表1)，并由此推算出模拟中所用无量纲参数的取值(表2)。模拟中涉及的参数取值如下：Ca1= Ca2= 0.1，Bo = 10-1，0.2 ≤σ ≤5，0.25 ≤μ ≤10，0.1 ≤ρ ≤10，γ = 5，λ= 1，E = 0.001，Σi= 0.001，h1∞= 1，h∞= 0.01。

表1 典型参数的取值或变化范围Table 1 Order of magnitude estimates for the relevant physical parameters

表2 无量纲参数定义及其数量级Table 2 Representation and order of magnitude estimates of the relevant dimensionless groups

2 结果与讨论

2.1 演化过程

图4 显示了不混溶液体表面上液滴蒸发的时间-空间演化过程，图中液滴下半部分曲线对应表示液-液界面h1，上半部分曲线对应液-气界面h2，h2和h1之间的部分即为由液体2 形成的液滴（其边缘存在一层厚度极薄的预置液膜），后面的图也是与此相同的表示方法。液滴在蒸发过程中接触线后移，液滴厚度和体积逐渐减小。由式(19)和式(20)可知在加热条件下，液滴界面温度与液滴厚度呈负相关。蒸发过程中液滴下表面与液滴上表面温度差引起表面张力梯度差，产生的Marangoni 效应使液体分子从低表面张力处流向高表面张力处，即从热区域流向冷区域。由于蒸发对液滴变形的影响，从液-气界面蒸发去除的液体分子必须得到补充以满足界面力平衡，导致液滴内部产生从下向上的质量迁移，并达到瞬时准平衡态[9]。随着蒸发的进行，更多的液体分子进入周围的空气中，因此液滴上表面的液膜高度随时间降低，下表面的液膜高随时间先降低后升高，液滴体积不断减小，最终完全蒸发。蒸发、重力、毛细力和Marangoni 效应的耦合作用是导致蒸发液滴出现此种现象的主要因素。

图4 液滴快速蒸发到最终形状的时间-空间演化Fig.4 Spatial-temporal evolution of the fast evaporation of a drop to a final shape

从图5可以看出，演化过程分为两个阶段：①t =0～0.1，“铺展主导”阶段，在该阶段内，液滴在毛细力和重力作用下缓冲到液体基底中，初始条件的影响被消除；②t = 0.1～100，“蒸发主导”阶段，液滴体积逐渐减小直至完全蒸发。液滴在液体基底上铺展，达到最大铺展半径，之后由于蒸发铺展半径逐渐减小，蒸发液滴在不混溶液体表面上的铺展是非钉扎的。在“铺展主导”状态下，接触线的前进速率大于接触线的后退速率，液滴的直径随时间而增大，而液滴体积由于蒸发而减少。这种趋势导致接触角减小，接触线运动阻力增大，从而导致“铺展主导”状态向“蒸发主导”状态过渡。

图5 h1最小值、h2最大值以及液滴半径xmax随时间的演化Fig.5 Temporal evolution of minimal h1，maximal h2，and the drop half-width xmax

2.2 蒸发数的影响

蒸发数E 对铺展过程的影响如图6 所示。从图中可以看出，随着蒸发数E 的增加，液-气界面高度逐渐降低，液-液界面高度逐渐升高。蒸发数E表示铺展时间尺度与蒸发时间尺度的比值，而蒸发时间尺度为液滴自初始静止时刻开始蒸发至完全干涸所需时间。E 越大，液滴蒸发得越快，质量损失也就越快，相同时刻液滴体积越小，液滴寿命也就越短。驱动液滴铺展的主要因素包括温度梯度导致的Marangoni 效应以及蒸发诱导的外向流动。由于蒸发对液滴变形的影响，从液-气界面蒸发去除的液滴分子必须得到补充以满足界面力平衡，液滴内部产生从下向上的质量迁移，导致液-液界面不断升高。

图6 t = 100时E对液滴形状的影响Fig.6 The effect of E on the drop shape at t = 100

2.3 表面张力比的影响

图7为t = 100时不同张力比σ对液滴形状的影响。如图所示，随着张力比的减小，液滴上表面比下表面变形更大[29]，接触角θ1逐渐减小，θ2逐渐增加。因此，在较小的σ值下，液滴大部分伸出液体基底，而对于较大的σ，液滴则“下沉”到液体基底某一深度。液滴在一定程度上穿透了液体基底，深度取决于表面张力比。

图7 t = 100时表面张力比对液滴形状的影响Fig.7 The effect of the surface tension ratio σ on the drop shape at t = 100

2.4 密度比的影响

密度比ρ 的影响如图8 所示。当液滴密度小于基底液体密度时，液体基底整体向上拱起呈凸形，三相接触点位于液体基底最高点。而对于ρ ＞1，液体基底整体向下凹陷，三相接触点位于液体基底的局部最低点，液滴在重力作用下会下沉到液体基底中，密度比越大，下沉深度越大。上层液滴对液体基底的挤压，导致液体基底厚度从液滴边缘开始向外逐渐增加。对于负铺展系数的液滴，考虑到液滴的尺寸足够小，无论其密度如何，液滴都能够漂浮[30]。当ρ = 1 情况下，接触线两边基底高度保持不变。显然，重力是导致液滴呈现图中所示现象的主要原因。

图8 t = 100时密度比对液滴形状的影响Fig.8 The effect of the density ratio ρ on the drop shape at t = 100

2.5 黏度比的影响

黏度比μ 对液滴动力学的影响如图9 所示。假设液滴上方的气体具有足够小的黏度和密度，不会对系统产生影响。图9(a)表明，液滴中心处的厚度变化过程均要经历先下降后上升然后再下降三个阶段，这是因为重力和毛细力、热毛细力和蒸发分别在不同阶段起主要作用。液滴在蒸发过程中逐渐变薄并产生不规则振荡运动[6]。早期阶段液滴在低黏度比下的流动性更强，导致液滴铺展更加迅速，对于较大的μ值，黏度比的增加会导致铺展和收缩速率的降低，同一时刻液滴最大厚度也相对较大。从图9(b)可以看出，铺展半径随着时间先增加后减小，系统从“铺展主导”阶段逐渐向“蒸发主导”阶段过渡。接触线以恒定速度运动的区域有两个：一个处于前进阶段，另一个处于后退阶段[9]。当液滴的黏度增加时，由于黏性耗散和摩擦而产生的阻尼效应会变大[31]，导致铺展和收缩速率的降低。同时由于蒸发的存在，阻碍了液滴的铺展，从而缩短了铺展过程。接触角的时间演化取决于液滴的铺展和蒸发以及液滴内部的质量传输。对于低黏度液体，初始运动主要由表面张力和初始液滴形状影响，由于黏度会导致表面张力的变化，因此液滴的初始接触角是不同的。图9(c)表明在铺展的早期阶段，随着液滴铺展，接触角减小，当减小至某一极限时，液滴前缘变得不稳定。结果表明，两种液体间的黏度比对最大液滴厚度、铺展半径和接触角均有一定影响。模拟所得接触角随时间的变化趋势与Sun 等[9]将体积10 μl 甲苯透镜放置在21℃去离子水表面后的观察结果相同，均为铺展阶段接触角减小较快，在蒸发过程中大部分时间保持不变，蒸发过程遵循恒定的接触角模式。

图9 黏度比对最大液滴厚度(h2 - h1)max、铺展半径xmax和接触角θ1和θ2的影响Fig.9 The effect of the viscosity ratio on the maximal drop thickness(h2 - h1)max，spreading radius xmax and contact angles θ1 and θ2

2.6 界面温度分布

通过润滑理论忽略量级较小的项，消去能量方程中的非稳态项，并计算液滴温度随高度的相对变化率。由式(18)和式(19)可知液-液界面(z=h1)和液-气界面(z=h2)温度分布和蒸发量为：

式(29)表明，蒸发通过影响液滴界面的温度分布进而影响表面张力以及液滴铺展。无论液滴是浸没还是漂浮在液体基底上，在垂直于基底的温度梯度下温度场需关于z 轴对称。空气中的热通量和动量通量较小，对液滴和液体基底的温度场和速度场无显著影响，液滴内部温度波动较小。Buffone[12]进行了乙醇表面FC-72 液滴的蒸发铺展实验，并用红外热像仪测得了FC-72 液滴表面温度的分布图。为了与Buffone[12]的实验结果进行对比，根据其文献中提供的物性参数计算了密度比ρ=2.1，热导率比λ=0.34，黏度比μ=0.43，表面张力比σ=0.46 时，液滴蒸发演化过程中的液-气界面温度分布和液-液界面温度分布以及温度梯度分布，如图10所示。图10(a)为不同时刻液滴界面温度分布，由于温度与高度呈反比，故液-液界面温度高于液-气界面，且在液滴中心处达到峰值，液-气界面则与之相反。图10(b)、(c)表明液-液界面的温度梯度大于液-气界面。由于边缘处蒸发更快，蒸发冷却效应使得液滴边缘较中心区域温度更低。计算所得液-气界面温度的数量级与Buffone[12]实验结果无量纲化之后液滴表面温度数量级一致，且液滴表面温度随时间的变化趋势也与实验结果一致，即表面温度均随时间而升高。

图10 液-液和液-气界面温度与温度梯度的演化过程Fig.10 Evolutions of liquid-liquid and liquid-gas interfaces temperatures and temperature gradients

2.7 与实验结果的比较

Dussaud 等[11]在23℃时观测了表面张力比分别为1和0.4，密度比分别为1和0.9的辛烷和甲苯液滴在去离子水表面的蒸发行为，发现液滴前缘以t1/2的幂律向前推进，即xmax- t1/2。为了与实验进行对比，采用相同的表面张力比和密度比计算所得的不混溶液体表面上蒸发液滴铺展早期阶段铺展半径xmax随时间变化关系如图11 所示(箭头指向σ 增大的方向，图中的实验结果采用本文的无量纲定义式进行了无量纲化)。图11 表明，实验与数值模拟结果具有相同的变化趋势，但存在一定的误差，原因可能在于模拟液滴是二维的，无法完全与实验进行比较。此外，由于实验未涉及加热，而加热会使液滴蒸发加快进而导致铺展速率减小，故在铺展的早期阶段，铺展速率模拟值比实验值要小，随着时间的进行，蒸发对液滴铺展的影响越来越大，在铺展的后期阶段铺展半径模拟值小于实验值。

图11 不同表面张力比σ时铺展半径的时间演化对数图Fig.11 Log-log plot of the temporal evolution of spreading radius for different surface tension ratio σ

3 结 论

采用润滑理论推导了液-气和液-液界面的演化方程。 通过数值模拟研究了上层液滴外形以及铺展半径的时间演化规律，以及蒸发、重力、毛细力和Marangoni效应之间的耦合作用对接触角、铺展半径、温度场和界面几何形状的影响。模拟结果表明演化过程存在两种状态，一种是“铺展主导”的状态，另一种是“蒸发主导”的状态。随着时间的推移，接触角逐渐减小，接触线运动阻力增大，“铺展主导”状态逐渐向“蒸发主导”状态过渡。蒸发量与液滴高度呈反比，蒸发数越大，液滴蒸发得越快。液滴高度随时间逐渐降低，表面温度则越来越高。液滴在一定程度上会伸出和穿透液体基底，深度取决于表面张力比和密度比。低黏度比下的流动性增强会导致更快的铺展，黏度比的增加会导致铺展和收缩速率的降低。液滴前缘首先沿液体基底表面前进，在达到最大接触直径后由于蒸发而逐渐后退。在铺展的早期阶段，铺展半径xmax随时间变化关系与实验中得到的xmax- t1/2规律具有相同趋势。