1.   引　言

当液滴处于非均匀加热固体表面上时, 温度梯度引起表面张力梯度驱动液滴从高温区域移向低温区域, 该现象称为热毛细迁移^[1]. 深入了解液滴热毛细迁移特性对于提高喷墨打印质量、控制润滑剂移动和提高冷凝效率等领域的应用至关重要^[2-5].

近年来, 通过改变固体表面性质, 如润湿性^[6]和粗糙度^[7]对液滴热毛细迁移进行调控成为研究热点. 目前可通过两种方法来改变固体表面的润湿性质: 一是在固体表面添加不同润湿性或润湿性梯度的涂层以驱动液滴自发移动^[8-11]; 二是在固体表面制造不同的结构图案实现加速或抑制液滴迁移^[11-16]. Dai等^[16]开展了固体表面改性操控液体润滑剂热迁移的实验研究, 其在超疏油表面制作了楔形超亲油槽道, 超亲油槽道的操控能力得到验证, 迁移距离相较于无涂层情况增加了2—3倍, 且减小楔形扩张角可强化液滴热迁移, 最后通过理论推导揭示了楔形槽道操控液滴的机制为热毛细力和结构力的综合作用. 除实验外, 研究者也尝试采用数值计算方法复现上述实验现象以揭示表面改性操控液滴热毛细迁移的机理. Fu等^[17]使用格子玻尔兹曼方法(LBM)复现了Dai等^[16]的实验现象, 系统地研究了液滴黏度、表面张力、初始半径和楔形轨道扩张角、润湿性对液滴热迁移的影响. 研究结果表明, 较高的表面张力和较低的黏度使液滴迁移速度更快, 初始驱动力随扩张角增大而增加, 使液滴初始移动速度加快, 而迁移阻力也随之增大, 导致最终迁移距离则随扩张角增大而减小. 所得计算结果与实验结果定性一致, 但计算所得迁移速度远高于实验结果, 其原因可能是受计算方法所限, 液滴与周围气体的密度比和黏度比无法按照实验中的真实情况设置. 因此, 有必要采取其他方法对液滴在润湿性差异表面上的热毛细迁移现象进行进一步深入研究, 以探明其内在机理.

基于润滑理论推导液膜厚度演化方程, 从而研究特定参数对液滴外形演化特征影响的方法在研究液滴热毛细迁移问题时被广泛采用. 该方法具有建模简便、计算量小等优点^[18], 同时无需预先设定接触线或接触角状态(接触线恒定或接触角恒定)就可计算出液滴铺展过程可能经历的任一阶段, 如定接触角铺展/收缩阶段、定接触线阶段或黏滑阶段^[19,20]. 该方法假设液滴厚度$\stackrel{~}{H}$相对于铺展方向特征尺度$\stackrel{~}{L}$为小量, 在此前提下将连续性方程、动量方程和能量方程进行数量级分析、略去高阶小量后推导出液膜厚度的演化方程, 借助数值求解可得到液滴接触角和接触半径及接触线移动速度. 其中接触线处的奇异性可通过两种方法解决. 一是在接触线处引入前驱膜, 如Gomba和Homsy^[21]所开展的液滴热毛细迁移研究, Beltrame等^[22]所开展的液滴在非均匀润湿性表面上的运动研究. 该方法在考虑壁面润湿性时引入分离压, 用润湿性对比强度系数定性表明基底亲疏水性. 二是在固液界面处采用滑移边界条件, 并在迭代计算中引入与动态接触角和平衡接触角有关的关系式来描述接触线的移动速度(下文简称滑移模型). 如Ehrhard和Davis^[23]研究了液滴在非均匀加热平面上的铺展和迁移特征; Smith^[24]建立了二维液滴在水平表面上的热毛细迁移模型, 并考虑接触角滞后的影响, 分析了液滴内部流场、自由界面形状和接触线运动; Karapetsas等^[25]分析了非等温倾斜固体基板上的液滴迁移动力学特性, 计算接触线移动速度时引入不同界面张力对温度的敏感性, 证明了温度引起的平衡接触角变化会引起复杂的动力学过程. 本课题组在此前研究中^[26,27]采用该方法分析了二维液滴在受热固体壁面上的润湿特性、铺展特征及传热特性, 从平衡接触角角度揭示了温度影响壁面润湿性及铺展过程的内部机理.

上述基于滑移模型的理论研究均为二维液滴情形, 但当液滴在润湿性受限轨道上进行铺展时, 在受限轨道两个方向上的铺展动力学特征势必有所不同. 因此, 二维液滴铺展模型不能准确反映其实际过程, 此时需考虑三维液滴在此平面上的运动过程, 但目前尚未有这方面研究的报道. 为此, 本文在课题组前期工作基础上, 针对液滴在润湿性受限轨道上的铺展过程, 将滑移模型拓展到三维情形. 分析液滴在非轴对称铺展时的迁移速度变化和热毛细迁移特性, 探究轨道宽度和轨道润湿性对热毛细迁移过程的影响.

2.   理论模型

2.1   问题描述

如图1(a)所示, 一液滴放置于沿x方向具有均匀温度梯度的壁面上, 且壁面润湿性在y方向存在差异, 以x轴为对称轴, $-\dfrac{1}{2}\widetilde{E}\leqslant y\leqslant \dfrac{1}{2}\widetilde{E}$为亲水区域(图1(b)中绿色部分), 其余为疏水区域(图1(b)中灰色部分), 形成润湿性受限轨道. 液滴为不可压缩牛顿流体, 在温度梯度和润湿性受限的作用下从高温区向低温区移动. 液滴最大厚度为$\widetilde{H}$, 流动方向尺度为 $\widetilde{L}$, 壁面温度为${\widetilde{T}}_{\rm{w}}$, 因$\varepsilon =\widetilde{H}/\widetilde{L}$很小(图1中的$\widetilde{H}, \widetilde{L}$仅为清晰起见, 不代表真实比例), 因此可应用润滑理论(上标˜表示有量纲, 下同). 假设液滴运动区域为$\widetilde{B}\times \widetilde{l}$, 轨道宽度为$\widetilde{E}$. 初始时液滴置于高温区, 定义无量纲相对轨道宽度$E=\widetilde{E}/\widetilde{d}$, $\widetilde{d}$为液滴初始直径. 如无特殊说明, 本文所给出的各种方程、边界条件、初始条件以及数值计算结果均为无量纲形式.

2.2   演化方程

液滴运动的控制方程包括连续性方程、动量方程和能量方程, 其无量纲形式为^[27]

$$\frac{\partial u}{\partial x}+\frac{\partial \upsilon }{\partial y}+\frac{\partial w}{\partial z}=0\text{, }$$

(1)

$$\frac{\partial p}{\partial x}=\frac{{\partial }^{2}u}{\partial {z}^{2}}\text{, }$$

(2)

$$\frac{\partial p}{\partial y}=\frac{{\partial }^{2}\upsilon }{\partial {z}^{2}}\text{, }$$

(3)

$$\frac{\partial p}{\partial z}=-\varepsilon Bo\text{, }$$

(4)

$$\frac{{\partial }^{2}T}{\partial {z}^{2}}=0\text{, }$$

(5)

式中, Bo为邦德数, 无量纲参数的定义式参见文献[27].

在气-液界面z = h处, 满足切向应力和法向应力平衡、热平衡及运动学边界条件,

$$\frac{\partial u}{\partial z}=\frac{\partial {\sigma }_{\rm{lg}}}{\partial x}\text{, }$$

(6)

$$\frac{\partial \upsilon }{\partial z}=\frac{\partial {\sigma }_{\rm{l}\rm{g}}}{\partial y}\text{, }$$

(7)

$$p=-C\left(\frac{{\partial }^{2}h}{\partial {x}^{2}}+\frac{{\partial }^{2}h}{\partial {y}^{2}}\right)\text{, }$$

(8)

$$\frac{\partial T}{\partial z}=-BiT\text{, }$$

(9)

$$\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}+\upsilon \frac{\partial h}{\partial y}=w\text{, }$$

(10)

式中, σ_lg表示气-液界面张力; Bi为毕渥数, 表示液滴内部导热热阻与液滴外气体侧对流换热热阻之比, $Bi={\widetilde{h}}_{\rm{g}}\widetilde{H}/\widetilde{\lambda }$, $\widetilde{\lambda }$表示液滴内部导热系数, ${\widetilde{h}}_{\rm{g}}$为气液界面对流换热系数.

在固体壁面z = 0处, 满足有滑移、无渗透边界条件及热平衡条件:

$$w=0\text{, }$$

(11)

$$\begin{split}& u=\beta \frac{\partial u}{\partial z}\text{, }~~~ \upsilon =\beta \frac{\partial \upsilon }{\partial z}\text{, } \end{split}$$

(12)

$$T={T}_{\rm{w}}=\varGamma x\text{, }$$

(13)

式中, C为毛细数; β和Г分别为滑移长度和温度梯度.

结合(9)式和(13)式对(5)式积分可得, 无量纲形式的气液界面温度为

$$T= \dfrac{\varGamma x} {1+Bih}.$$

(14)

无量纲形式的界面张力与温度关系为

$$\begin{split}&{\sigma }_{\rm{sg}}={\delta }_{\rm{sg}}\left(1-{{\varOmega }}_{\rm{sg}}T\right) \text{, } \\ &{\sigma }_{\rm{ls}}={\delta }_{\rm{ls}}\left(1-{\varOmega }_{\rm{ls}}T\right) \text{, } \\ &{\sigma }_{\rm{lg}}={\delta }_{\rm{lg}}\left(1-{\varOmega }_{\rm{lg}}T\right)\text{, }\end{split}$$

(15)

式中, 角标sg, ls, lg分别代表气-固界面、液-固界面、气-液界面; 其中${\delta }_{i}={\sigma }_{i, T0}/{\sigma }_{{\rm{l}\rm{g}}, T0}$(i = sg, ls, lg), σ_{lg, T0}为气-液界面在T₀下的界面张力, T₀为原点处的温度, ${\varOmega }_{\rm{lg}}$为气-液界面温度敏感系数.

对上述方程组和边界条件推导, 可得液膜厚度演化方程为(详见文献[27])

$$\begin{split}{h}_{\rm{t}}=&-\nabla \cdot \left[\left(\frac{{h}^{3}}{3}+\beta {h}^{2}\right)\nabla \left({C\nabla }^{2}h-\varepsilon Boh\right)\right.\\ &-\left.{\varOmega }_{\rm{l}\rm{g}}T\left(\frac{{h}^{2}}{2}+\beta h\right)\right]\text{, }\end{split}$$

(16)

其中$\nabla =\dfrac{\partial }{\partial x}{{\boldsymbol{i}}}+\dfrac{\partial }{\partial y}{{\boldsymbol{j}}}$, 式中右侧三项分别代表毛细力、重力、热Marangoni应力在演化过程中的作用.

2.3   接触线动力学

在理想光滑表面上, 液滴达到平衡状态时, 在三相接触线处满足Young方程^[28]:

$${\widetilde{\sigma }}_{\rm{sg}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right)={\widetilde{\sigma }}_{\rm{lg}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right)\cos {\widetilde{\theta }}_{\rm{e}}+\widetilde{\sigma }_{\rm{ls}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right)\text{, }$$

(17)

式中, ${\widetilde{\theta }}_{\rm{e}}$为有量纲的平衡接触角. (17)式相应的无量纲形式为

$${\sigma }_{\rm{s}\rm{g}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right)={\sigma }_{\rm{l}\rm{g}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right)\cos{\theta }_{\rm{e}}+\sigma_{\rm{l}\rm{s}}\left({x}_{\rm{c}},{y}_{\rm{c}}\right).$$

(18)

液滴铺展/迁移过程中, 接触线处的液膜厚度始终为零, 即h = 0; 沿接触线处做液滴外形的切线, 可得液滴的动态接触角θ. 三维液滴情形下, x方向和y方向接触线的移动速度u, υ与动态接触角θ和平衡接触角θ_e间的关系如下:

$$u=\frac{{\rm{d}}x}{{\rm{d}}t}=\left\{ \begin{array}{ll}{B\left(\theta -{\theta }_{\rm{e}}\right)}^{m},&\theta > {\theta }_{\rm{e}},\\ {-B\left({\theta }_{\rm{e}}-\theta \right)}^{m}, & \theta < {\theta }_{\rm{e}},\end{array}\right.$$

(19)

$$\upsilon =\frac{{\rm{d}}y}{{\rm{d}}t}=\left\{ \begin{array}{ll}{B\left(\theta -{\theta }_{\rm{e}}\right)}^{m},&\theta > {\theta }_{\rm{e}},\\ {-B\left({\theta }_{\rm{e}}-\theta \right)}^{m},&\theta < {\theta }_{\rm{e}},\end{array}\right.$$

(20)

式中, B和m为常数, 通常取B = 0.001, m = 3; u(υ) > 0代表向外铺展, u(υ) < 0代表向内收缩.

值得注意的是, 在流动控制方程简化过程中采用了润滑理论, 液滴厚度方向上的无量纲参数是基于有量纲参数与$\varepsilon \widetilde{L}=\widetilde{H}$相比得到的. 而液滴x方向、y方向上的无量纲参数是基于有量纲参数与液滴半径$\widetilde{L}$相比得到的, 且假定$\varepsilon =\widetilde{h}/\widetilde{L}\rm{很}\rm{小},$这使得无量纲情形下, 液滴的纵横比被放大, 同时液滴的动态接触角θ也被相应放大^[27]. 因此液滴的平衡接触角θ_e也需要放大, 以便采用(19)式和(20)式确定液滴的铺展/定扎/收缩等状态. 平衡接触角的有量纲形式${\widetilde{\theta }}_{\rm{e}}$与无量纲形式θ_e的关系为: $\cos{\widetilde{\theta }}_{\rm{e}}\to \cos\left(\varepsilon {\theta }_{\rm{e}}\right)$, 又有$1- ({\varepsilon }^{2}{{\theta }_{\rm{e}}}^{2} )/2 \approx \cos (\varepsilon {\theta }_{\rm{e}})$, 结合(15)式、(17)式和(18)式可得

$$\begin{split} {\theta }_{\rm{e}} ^{2}=&-\frac{2}{{\varepsilon }^{2}}\left[\frac{A\left(1-{\varOmega }_{\rm{sg}}\varGamma {x}_{\rm{c}}\right)}{\left(1-{\varOmega }_{\rm{lg}}\varGamma {x}_{\rm{c}}\right)}\right.\\ &\left.+\frac{\left({\varOmega }_{\rm{ls}}-{\varOmega }_{\rm{sg}}\right){\delta }_{\rm{ls}}\varGamma {x}_{\rm{c}} + \left({\varOmega }_{\rm{lg}} - {\varOmega }_{\rm{sg}}\right)\varGamma {x}_{\rm{c}}}{\left(1-{\varOmega }_{\rm{lg}}\varGamma {x}_{\rm{c}}\right)}\right]\text{, }\end{split}$$

(21)

式中, A = δ_sg – δ_ls – 1. 本文的计算中假定$\varOmega_{\rm{ls}}$, $\varOmega_{\rm{sg}}$和$\varOmega_{\rm{lg}}$始终等于1, 平衡接触角则为定值:

$${\theta }_{\rm{e}}^{2}=-\frac{2}{{\varepsilon }^{2}}A.$$

(22)

通过A的不同取值来改变平衡接触角大小, 进而划分为亲疏水区域. 本文计算中, 亲水区域平衡接触角θ_eq分别取值为0.1410, 316和0.446 (单位为弧度), 对应的角度值分别为8.1°, 18.1°和25.6°; 疏水区域平衡接触角θ_es取值为1.9 rad, 对应的角度值为109°.

液滴铺展/迁移的计算过程主要步骤如下. 首先设定时间迭代步长Δt, 通过结合上一时刻液膜厚度外形, 求解液膜厚度演化方程(16), 得到当前时刻接触线内侧邻近点的液膜厚度, 进而求得气-液界面切线斜率tanθ及动态接触角θ. 因为液滴沿x, y两个方向铺展, 所以选取M₁—M₄四个点作为接触角监测点, 如图2所示, 记录下各时刻四监测点处的θ值. 最后通过(19)式—(20)式得到液滴在不同方向上的接触线移动速度和下一时刻接触线位置.

2.4   数值计算方法

初始时设置液滴为三维轴对称几何体, 与固体壁面接触的底面为半径等于1的圆形, 垂直于壁面的剖面具有抛物线形状, 最大液滴厚度为1, 即

$$\begin{split} & h\left(x,y,t=0\right) \\ =&\left(1-{x}^{2}-{y}^{2}\right)\left[F\left(1-\sqrt{{x}^{2}+{y}^{2}}\right)\right.\\ &\left.-F\left(-1-\sqrt{{x}^{2}+{y}^{2}}\right)\right]\text{, }\\[-15pt]\end{split}$$

(23)

式中, F(x) = 0.5[1 + tanh(20x)]为海氏阶跃函数.

基于Karapetsas等^[25]提出的坐标变换法, 将其扩展到二维, 将液滴底面不断变化的瞬态物理域$(x, y, t)$映射到固定计算域$(x', y', {t}')$, 以避免迭代过程中移动网格, 从而提高计算效率:

$${x}'=\left(x-{x}_{\rm{r}}\right)\frac{{x}_{\rm{a}}\left(t=0\right)-{x}_{\rm{r}}\left(t=0\right)}{{x}_{\rm{a}}-{x}_{\rm{r}}}\text{, }$$

(24)

$${y}'=\left(y-{y}_{\rm{c}\rm{l}}\right)\frac{{y}_{\rm{c}\rm{r}}\left(t=0\right)-{y}_{\rm{c}\rm{l}}\left(t=0\right)}{{y}_{\rm{c}\rm{r}}-{y}_{\rm{c}\rm{l}}}\text{, }$$

(25)

式中, x_r和x_a分别为液滴后退接触线(沿x方向左侧接触线)和前进接触线(沿x方向右侧接触线)位置, y_cl和y_cr分别为液滴沿y轴运动方向左侧和右侧接触线位置. 因计算中涉及坐标变换, 需将(16)式中的时间导数用下式替换:

$$\frac{\partial }{\partial t}=\frac{\partial }{\partial {t}'}-\frac{{\rm{d}}{x}'}{{\rm{d}}t}\frac{\partial }{\partial {x}'}-\frac{{\rm{d}}{y}'}{{\rm{d}}t}\frac{\partial }{\partial {y}'}.$$

(26)

采用有限元分析软件FreeFEM++模拟液滴热毛细迁移的动态演化过程. 将半径为1的圆形固定计算域$(x', y', {t}')$的周长分别均匀划分成120, 240和300段, 从而在圆形区域内自动生成4584, 10028和15578个网格, 并进行了网格无关性验证, 如表1所列. 以t = 40000时, 网格数为15578的计算结果为基准, 网格数为4584和10028对应的液滴前进接触角偏差分别为1.9%和0.2%, 接触线移动速度偏差分别为3.3%和0.1%; 但与网格数为10028相比, 网格数为15578时计算时间增加约40%, 故下文计算采用10028个网格.

模拟所需物性参数的典型数量级列于表2, 其无量纲参数取值范围见表3. 如无特殊说明, 各参数取值如下^[27]: C = 0.1, ε = 0.1, A = –0.001, β = 1 × 10^–5, Bo = 0.5, Bi = 0.01, K = 0.001, m = 3, Г = 0.02.

2.5   计算模型验证

为验证本文模型的正确性, 模拟了Dai等^[29]开展的3 μL石蜡油在非均匀加热无润湿性差别固体表面上的迁移过程, 并与其实验结果进行对比. 为此, 将Dai等^[29]的实验数据根据本文的无量纲参数定义式换算后得到各无量纲参数为C = 0.03, ε = 0.1, Bi = 0.01, A = –0.0005, Bo = 0.5, Г = 0.005. 其余计算所需无量纲参数取值为$\beta = 1 \times 10^{-5}$, K = 0.001, m = 3. 结果如图3所示, 图3中上半部分为实验结果, 下半部分为模拟结果. 图3中时间进行了归一化处理, 即分别将实验和模拟中涉及的各迁移时间除以迁移至最终位置的对应时间, $\widehat{{t}_{\rm{s}}}$和$\widehat{{t}_{\rm{e}}}$分别表示归一化后的模拟时间和实验时间. 由图3可见, 液滴的轮廓从开始的圆形逐渐演化为最终时刻的近似椭圆形, 迁移过程中形成了一层薄液膜拖在液滴主体之后, 模拟所得液滴外形轮廓达到与实验结果一致时, 内部的液滴高度变化也基本与实验结果一致, 由此可确定本文所用计算方法的可靠性.

模拟所得液滴迁移距离随时间的演化结果与Dai等^[29]的实验结果对比如图4所示, 图中时间和迁移距离按照图3中处理迁移时间的方法进行了归一化处理. 由图4可知, 模拟结果与实验结果总体吻合良好, 但计算所得液滴迁移距离变化较为均匀, 而实验中前期液滴迁移距离增长较快, 后期较慢. 同一时刻, 迁移距离计算结果相对于实验结果的最大偏差为18%. 引起偏差的原因可能为: 采用润滑理论时假定液膜厚度远小于它的迁移方向尺度, 即$\varepsilon\ll 1$, 这是一种近似计算方法, 目的是略去运动控制方程中以惯性项为代表的非线性项及高阶小量. 有研究表明, 液固接触角越小, 润滑理论的计算结果越准确^[18], 实验中液滴的初始接触角为11°—18°, 相当于液滴纵横比ε = 0.10—0.15, 略高于本文取值0.1. 而且本文在计算时假定$\varOmega_{\rm{ls}}$, $\varOmega_{\rm{sg}}$和$\varOmega_{\rm{lg}}$始终等于1, 也就是壁面润湿性不受温度影响(平衡接触角为定值), 而实际上壁面润湿性在温度不同时可能会发生变化, 因此使得计算结果存在一定偏差.

为了进一步验证本文计算方法的正确性, 模拟了Karapetsas等^[25]所研究的非等温基板上的液滴迁移过程, 参数取值与文献[25]所用一致, 即C = 0.1, ε = 0.1, Bi = 0.01, A = –0.001, Bo = 0.5, Г = 0.005, m = 3. 液滴轮廓演化过程对比结果如图5(a)所示. 从图5(a)可以看出, 模拟复现的结果与文献计算结果十分吻合, 在各个时刻液滴外形轮廓基本一致. 图5(b)为液滴左侧接触线(x_cl)和右侧接触线(x_cr)移动距离随时间变化, 模拟得到的x_cl变化趋势与文献结果基本一致, x_cr的最大偏差为22%. 产生偏差的原因主要是, 文献[25]的计算采用的二维液滴, 而本文模型为三维液滴, 因此, 结果会存在一定偏差.

3.   结果和讨论

3.1   液滴演化历程

当无量纲相对轨道宽度(下面简称轨道宽度)E为0.8, 亲水区域平衡接触角θ_eq为0.446, 液滴的演化历程如图6所示. 下文所有描述均以液滴迁移方向(x轴正方向)称为前, 反之称为后, y轴正方向称为左, 负方向称为右. 液滴刚开始铺展时(t =100), 呈轴对称分布, 如图6(a)所示; 在终了时刻t = 40000, 如图6(b)所示, 液滴主体向前迁移一段距离, 液滴高度下降, 后面尾部有一明显的“拖尾”现象. 图6(c)表明, x方向上液滴外形演化过程呈现两个重要特征: 一是由热毛细力形成的Marangoni效应促使液滴主体部分向低温区移动, 同时液滴在移动过程中后缘液膜厚度随时间逐渐降低, 与液滴主体部分间形成一层薄液膜(图6(c)中区域A), 这与Mukhopadhyay等^[30]的实验结果相符; 二是液滴中心厚度(最大液膜厚度) 先迅速下降(t < 5000), 之后下降速度明显减缓(5000 ≤ t ≤ 20000), 最后趋于不变(20000 < t ≤ 40000). 值得注意的是, 液滴在y方向上的铺展受到了明显抑制, 如图6(d)所示, 初始时刻液滴左右两侧位于疏水区域, 受壁面润湿性限制, 液滴左右两侧接触线迅速向内收缩到轨道边缘的位置(y = ±0.8), 继而一直保持定扎状态.

为更清楚地分析液滴接触线的移动特征, 图7给出了液滴接触线位置随时间的变化. 为了深入分析液滴迁移过程中各主导因素的作用机理, 将液膜厚度演化方程(16)式分解成重力项(F_G)、毛细力项(F_cap)和热Marangoni应力项(F_thermal), 即

$${F}_{\rm{G}}=-\nabla \cdot \left[C\left( \frac{{h}^{3}}{3}+\beta {h}^{2}\right)\nabla {\nabla }^{2}h\right],$$

(27)

$${F}_{\rm{c}\rm{a}\rm{p}}=\nabla \cdot \left[\varepsilon Bo\left( \frac{{h}^{3}}{3}+\beta {h}^{2}\right)\nabla h\right],$$

(28)

$${F}_{\rm{t}\rm{h}\rm{e}\rm{r}\rm{m}\rm{a}\rm{l}}=\nabla \cdot \left[{\varOmega }_{\rm{l}\rm{g}}T\left(\frac{{h}^{2}}{2}+\beta h\right)\right].$$

(29)

在不同时刻分别计算上述三项在后退接触线和前进接触线处的值, 绘制成图7(c)和图7(d). 图7(a)表明, 液滴迁移时的前进接触线x_a始终向低温区移动, 移动距离与时间满足x_a~ t^0.7关系, 相应的图7(d)中前进接触线处F_G和F_cap与F_thermal均为正值, 表明三者共同驱动液滴铺展, 三者之和渐趋减小; 从图7(a)还可以看出, 液滴后退接触线x_r变化呈现三个阶段: 首先在t ≤ 2500时x_r向后(高温区)移动, 相应的图7(c)中(虚线左侧)后退接触线处F_G和F_cap小于0, 驱动后退接触线向后移动, F_thermal大于0, 驱动后退接触线向前移动, 但F_thermal + F_G + F_cap小于0, 说明初始阶段重力和毛细力在此区域占据主导地位, 驱动液滴向后迁移; 其次在2500 < t ≤ 6000时图7(a)中x_r保持不变, 相应的图7(c)中F_thermal + F_G + F_cap近似为0, 三者达到平衡, 表明重力、毛细力和热Marangoni应力在此阶段达到了平衡, 致使后退接触线处于定扎状态; t > 6000后, 图7(a)中x_r转为向前 (低温区)移动, 图7(c)中后退接触线处三项作用力均较弱, 结合前文的图6(b)和图6(c)可知, 此时的液滴主体部分已经向前移动, 与液滴后缘之间拖拉出一层薄液膜, 且随时间延续, 薄液膜覆盖面积不断增大, 说明x_r的增长是被液滴主体拖动的结果. 此外, 图7(a)中的液滴前后铺展范围x_domain (= x_a – x_r) 随时间始终增大, 在所有时间范围内满足~t ^0.36关系. 由图7(b)可以看出, 受润湿性限制, 液滴左右两侧的铺展受到抑制, 即液滴左右两侧的接触线在开始时均迅速向内收缩(阶段Ⅰ), 移动至轨道边缘(y = ±0.8), 之后保持定扎(阶段Ⅱ). 液滴两侧与壁面接触范围y_domain在阶段Ⅰ迅速减小, 之后则一直保持不变.

图8为动态接触角随时间的变化. 由图8(a)可知, 前进接触角θ_a开始时随时间迅速下降至θ_a=0.9 (0 ≤ t ≤ 6000), 之后缓慢降低, 接近于恒定. 对比θ_a与平衡接触角θ_eq可以发现, θ_a始终大于θ_eq, 所以x_a始终向前移动. 在0 ≤ t ≤ 7500内, 后退接触角θ_r随时间逐渐减小, 在t > 7500后则基本不变. 这与液滴后缘薄液膜的形成有关, 尤其是当后退接触线开始向低温区移动, 此薄液膜就一直“拖”在液滴主体部分之后基本不变(如图6(c)所示), 导致后退接触角θ_r基本不再变化. 对比θ_r与θ_eq可以发现, 在0 ≤ t < 4400内, θ_r > θ_eq, 当t ≥ 4400后, θ_r降至θ_eq以下, 于是后退接触线先向后移动, 之后随液滴主体部分向前移动, 这与上文所说重力和毛细力、热毛细力先后主导该区域铺展的特征相吻合. 因液滴关于x轴对称, 所以左右两侧接触角的变化完全一致, 故此处仅给出液滴左侧接触角θ_yl的时间演化特征, 如图8(b)所示. 可以看出, θ_yl在开始时短暂升高, 之后逐渐下降. 受润湿性限制, θ_yl始终大于θ_eq.

3.2   轨道宽度E的影响

Chowdhury等^[31]发现, 轨道宽度对液滴在润湿性限制表面上的自发迁移过程有很大影响, 且影响规律并非单调关系, 即相对于更宽的轨道或更窄的轨道, 液滴迁移速度在中等宽度的轨道上达到最快. 为探究润湿性受限轨道宽度对液滴热毛细迁移特性的影响, 分别对E = 0.4, 0.6, 0.8, 1.0和2.0等情形下液滴的迁移过程进行模拟(液滴的初始半径均为1), 其中E = 1.0时轨道宽度与液滴直径相等, E = 2.0时轨道宽度为液滴直径的2倍, 亲水区域平衡接触角 θ_eq为0.446.

图9为不同轨道宽度时接触线位置、接触线移动速度以及接触角随时间的变化曲线. 图9(a)表明, 在不同E情形下, 后退接触线x_r随时间的变化大体相同, 均经历了先向后铺展、后定扎、最后随液滴整体向前移动的过程, 轨道宽度仅影响液滴向后铺展阶段的延续时间和铺展距离. 图9(b)表明, 液滴前进接触线x_a随时间的增长关系近似满足指数关系x_a ~ t^b, 指数b与轨道宽度E负相关, 轨道宽度E从0.4增长到2时, 指数b分别为0.72, 0.71, 0.70, 0.69, 0.68, 即轨道越窄, 液滴热毛细迁移速度越快.

图9(c)表明, 多种轨道宽度情形下液滴后退接触线移动速度随时间的变化规律总体相同, 均经历速度绝对值减小→停滞→速度增加→趋于稳定的过程. 在迁移早期, 后退接触线移动速度绝对值|u_r|均从2.4×10^–4迅速减小至0, 且随E增大, |u_r|在更短时间内减小至0. 图9(d)表明, 多种轨道宽度情形下液滴前进接触线的移动速度主要经历快速下降和缓慢下降两个阶段, 且随E增加, u_a减小, 并更早结束u_a快速下降阶段进入缓慢下降阶段.

结合图9(e)和图9(f), 可以看出轨道宽度对接触线、接触线移动速度的影响机理. 以E = 0.4的情形为例, 在最初的0—0.28 × 10⁴时间段, 液滴左右两侧边缘从疏水区域收缩到轨道边缘, 导致y_cr下降(图9(e)), 液滴y方向截面变窄, 液滴高度增加. 由于液滴为不可压缩流体, 液滴左右两侧多余的液体除了令液滴增高, 也被排挤到x方向, 加剧液滴前后向外铺展, 即x_r减小(图9(a)), x_a增大(图9(b)), 直到t = 0.28 × 10⁴时刻收缩结束. 此时图9(c)中|u_r|减小的曲线出现了第一个拐点, 图9(d)中u_a出现小幅上升. 此后, 0.28 × 10⁴—0.5 × 10⁴时间段, 毛细力驱动液滴克服热毛细力继续向后小幅铺展, 但铺展速度|u_r|继续减小, 直到t = 0.5 × 10⁴时刻u_r = 0, 与此相对应, 图9(a)中x_r持续减小. 之后0.5 × 10⁴—0.7 × 10⁴时间段, 液滴后退接触线经历短暂的定扎阶段, u_r = 0, x_r = –1.67. 由于对于液滴的前进接触线, 毛细力和热毛细力对于液滴向前迁移铺展均起促进作用, 因此x_a的增长在0—0.28 × 10⁴时间段(受到左右两侧收缩影响)快于之后的时间段. 从t = 0.7 × 10⁴时刻开始, 在热毛细力作用下, 液滴后退接触线随同液滴主体向前移动, u_r增长至0.5 × 10^–4左右, 经历小幅振荡后, 保持0.5 × 10^–4的移动速度基本不变, 而前进接触线移动速度u_a在t = 0.7 × 10⁴时刻也进入缓慢下降阶段. 由于排挤到x方向的液体随轨道宽度的增加而减小, 液滴高度随之降低, 前进接触角θ_a也随之减小, 因而图9(d)中u_a随E的增加而减小. 同时, 液滴左右两侧收缩到轨道边缘所需时间随E增加而缩短, 图9(c)中|u_r|减小曲线的拐点位置则随E增加向左移动, |u_r|减小到0的时刻逐渐提前, x_r随同液滴主体向前移动的时刻也相应提前, u_r开始增长的时刻也提前. E = 1和E = 2时, 液滴左右两侧不经历收缩阶段, |u_r|减小曲线也就不存在拐点. E = 2时, y_cr最大移动距离为1.6, 未超出轨道边缘, 也就是说液滴始终在亲水区域铺展, 相当于无润湿性限制情形. 轨道宽度E对θ_r稳定值的影响较小, E = 0.4时, θ_r = 0.08; E = 2时, θ_r = 0.1. 因此E = 0.4时u_r略大(图9(c)).

3.3   轨道润湿性的影响

以上分析均基于亲水区域接触角θ_eq = 0.446情形, 下面在E = 0.6时, 将θ_eq分别设置为0.446, 0.316和0.141, 进一步探究轨道润湿性对液滴迁移特性的影响, 其中接触角越小代表润湿性越强.

图10为轨道润湿性对液滴接触线移动速度的影响. 图10(a)和图10(b)表明, |u_r|和u_a总的变化趋势与3.2节中的描述相同, 其中 |u_r|和u_a的初始值随润湿性增强而增大. 这是因为当轨道润湿性增强, 平衡接触角减小, 液滴初始接触角与平衡接触角之间的差值增加, 由(20)式可知接触线初始移动速度加快. 图10(a)表明, 不同润湿性时u_r最终稳定值随润湿性增强而略有减小, 原因是u_r进入稳定状态时, θ_r < θ_eq, 由(19)式, θ_eq越大, u_r越大. 图10(a)说明轨道润湿性变化对u_r初始值影响更显著. 图10(b)表明, u_a随润湿性的增强而加快, 原因同样可用(19)式解释, 当θ_a > θ_eq, θ_eq越大, u_a越小. 这与Fu等^[17]所研究的液滴在楔形轨道上的热毛细迁移结果一致. 由图10(b)还可以看出, 在u_a缓慢下降阶段, t = 4 × 10⁴时, 与θ_eq = 0.446和0.141对应的u_a分别为0.78 × 10^–4, 1.32 × 10^–4, 即润湿性增强3倍后, u_a提高约70%. 而从图9(d)中可知, 相较于E = 1, E = 0.4时u_a提高约90%. 由此可知相比于改变轨道润湿性, 改变轨道宽度对u_a的提升效果更显著; 而且在实际应用中, 制备亲水和超亲水表面的难度往往比改变轨道宽度要大得多^[32], 因此改变轨道宽度可能是调控液滴热毛细迁移更好的方法.

4.   结　论

液滴在润湿性受限轨道的热毛细迁移过程呈现三个特征: 液滴主体由高温区向低温区移动, 前进接触线移动距离随时间成指数关系; 液滴后缘在移动过程中液膜厚度随时间逐渐降低, 与液滴主体部分间形成一层薄液膜; 液滴中心厚度先迅速下降, 后趋于不变, 轨道的润湿性限制导致液滴在垂直于轨道方向(y方向)上的铺展受到抑制, 接触线移动至轨道边缘后保持定扎.

液滴前进接触线移动速度与轨道宽度成负相关. 在不同轨道宽度下, 前进接触线移动速度开始时均随时间迅速减小, 之后缓慢下降趋于平稳; 后退接触线均经历向后铺展、然后定扎、最后随液滴整体向前移动这一特征. 液滴直径大于等于轨道宽度时, 液滴y方向接触线均随时间收缩到轨道边缘后保持定扎. 垂直于轨道方向(y方向)的壁面润湿性限制导致的排挤作用, 在初始的短暂时刻对液滴在轨道(x方向)上的热毛细迁移具有加速作用, 且这种作用随轨道宽度的减小而增强.

增强轨道润湿性使得后退接触线移动速度的初始值增大, 但对其稳定值影响不大; 前进接触线移动速度和轨道润湿性呈正相关. 相较于改变轨道宽度, 改变轨道润湿性对前进接触线移动速度的影响较小, 因此改变润湿性轨道宽度可能是调控液滴热毛细迁移更好的方法.