液滴在壁面上受到重力、离心力、剪切力等而发生定向运动是一个十分常见的物理现象，如下雨时，房间的窗户上、荷叶的叶面上、汽车的车窗表面都有各种因素导致的不同的液滴滑落现象。实际生产生活中，人们接触到的表面往往不是均质和绝对光滑的，通常都有因表面粗糙度、表面化学特性或表面形貌特征变化所导致的局部缺陷；尽管这些缺陷的尺度可能十分微小，但却能够影响其表面特性，并阻止或改变液滴的运动过程。

不少学者关注到此问题并试图将此类表面特性利用到微流体控制以及航天飞机表面除冰等工程领域。Mannetje等^[1]通过对壁面加载电压的方式使得倾斜的壁面上产生局部润湿特性改变，研究了液滴在倾斜壁面上跨越润湿缺陷区液滴运动变化过程，获得了液滴被润湿缺陷捕获的临界条件；Priest等^[2]研究了因表面化学特性异质而导致的多个离散区域的润湿性的变化，并将其与连续润湿变化进行了比较，研究发现只有在离散润湿变化的情况下液滴才能发生黏滑运动；Sbragaglia等^[3]通过对玻璃涂抹多条化学物质，使得玻璃表面的润湿特性发生周期性的变化，获得了液滴在非均质倾斜表面的黏滑运动规律。

不难发现，实验研究过程制备的非均质表面通常需要采用化学试剂涂抹、加载电压或改变微观结构等物理化学手段^[4]，而这些实验材料制备过程的难度较大，数据采集不够精确，实验周期较长。因此，仅仅依靠实验手段来进行研究，并不能够准确对其进行理论解释。而数值模拟方法却能够很好地解决实验中所存在的局限性，并且润湿特性在数值模拟领域已经取得了较多成果^[5]。Bussmann^[6]和Fukai^[7]等采用流体体积(Volume of Fluid, VOF)方法解决了流动过程中三相接触线的移动问题，并且讨论了界面润湿特性对液滴铺展过程的影响；Spelt^[8]采用水平集方法(level-set)、Khatavkar等^[9]采用扩散界面场(diffuse-interface)模型也都找到了因三相接触线导致的尖锐界面解决方法，并能够模拟出不同润湿特性下液滴的铺展过程。然而，上述方法对三相接触角的解决方案大都是在接触点处预置一个角度再进行后续计算，但这样会导致在接触点处的流场并不真实，并且处理壁面上的离散润湿缺陷并不方便。格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)能够通过设置壁面力来使得液滴在上面自然形成一个与壁面力相对应的接触角，其是能够解决离散非均质壁面的方法。Sbragaglia等^[3]采用LBM与其实验结果进行比较，并对液滴的黏滑运动进行了理论分析，但该方法的缺点在于界面没有精确定义，因而接触角大小及接触点处的速度变化也难以精确捕捉，无法弥补实验研究的缺陷。

界面追踪法(Front Tracking Method, FTM)^[10]是一种较为经典的数值模拟方法，不同于其他界面追踪方法，FTM将界面力都储存于界面标记点中，界面自然形成，无需计算界面曲率来构造界面，这也为其在三相接触线处的尖锐界面的处理变得简单。Yamamoto^[11-12]、Muradoglu^[13]、Hu-ang^[14]等都采用FTM研究了液体在固体壁面上运动时接触角的变化规律，并得到了非常准确的角度及清晰的界面。但上述研究均只针对FTM来构建实现三相接触角的模型，并对单一接触角情况进行模型验证，而基于FTM对润湿性间断变化的非均质表面的研究报道却并未发现。

因此，本文借助FTM界面描述优势，建立了基于FTM的二维接触角模型。通过该模型，针对非均匀润湿的壁面液滴运动规律进行探索，通过改变流体参数及表面非均质程度，分析不同因素导致的液滴运动规律改变，从而掌握非均质表面对液滴运动的控制规律，为表面修饰及微流控等实际工程应用领域提供帮助。

1 计算模型与验证 1.1 相界面跟踪的Navier-Stokes方程及处理

控制方程应用到整个计算区域时，表面张力都集中在界面上，考虑引入一源项并将其看作体积力再与δ函数相乘，该δ函数只有在界面上时才不为零。此时，二维不可压缩流体的动量方程可写为

式中：s为界面元；δ(x－x_f)为二维的狄拉克函数，x_f代表界面位置；n为界面法向量；σ为界面张力系数；κ为界面平均曲率；u为速度；ρ和μ分别为非连续的密度和黏度；p为压力；g为重力。

为了求解方程式(1)，将其分解，把压力项单独分离出来，并引入一个临时速度u^*，于是得到只包含对流项A和扩散项D的方程：

式中:Δt为时间步长。

求解后得到不含压力的临时速度场u^*，再将压力项加入得到通用的速度离散方程：

式中：▽_h为取步长为h的离散形式哈密顿算子，用于解决此处不连续的梯度分布。n+1时刻的速度应满足离散形式的质量守恒，即

将式(4)与式(3)联立可得压力Poisson方程：

通过超松弛迭代法(successive over relaxation)对压力进行迭代求解，迭代系数β=1.2。表面张力的计算采用连续表面力(continuous surface force)方法，表达式为

取单位线元上的表面张力来研究，对于二维流动，有

得到单位界面元上的表面张力为

式中：t为界面上的切向量；l为界面上的节点。

1.2 密度和黏度重构

对于不互溶、不可压缩的流体，相界面两侧流体保持各自的特性，因而界面上存在物性阶跃。界面移动时，界面上的密度和黏度分布也随之变化，故引入Heaviside函数作为指示函数来表征这种变化。

式中：φ为给定点到界面的距离；2α为两流体间过渡区域厚度的绝对值。

因此，密度、黏度的分布表达式为

式中：ψ₁和ψ₂分别为2种流体的密度和黏度。

1.3 界面及接触点移动

界面通常与固定网格结合在一起，界面与网格之间通过面积权重函数实现参数转换。采用双线性插值法使得从界面到固定网格中的突变密度和表面张力光滑变化，最后使用一个简单的一阶显示时间积分处理，便可获得下个时刻的界面位置。

式中：x_f⁽ⁿ⁺¹⁾为n+1时刻界面的位置；x_fⁿ为n时刻界面的位置；u_fⁿ为n时刻界面的速度。为保证数值模拟的计算精度，将式(11)中关于时间的微分都采用二阶精度的一阶向前差分。

在FTM中，界面标记点位置由标记点速度随时间积分确定。因此，无需提前设置气液界面与壁面处夹角，由不平衡力导致的流动过程会使得界面形状自然形成。为使壁面处标记点获得速度从而进行移动，本文使用广义滑移边界(generalized Navier boundary condition)^[15]处理壁面处接触点的移动过程，表达式为

式中：u_cl为接触点滑移速度；λ为壁面滑移长度；μ_l为液相黏度；x^w和x_f分别为壁面处流体单元和界面坐标；f ^Y为因壁面润湿特性而产生的不平衡弹性力，即

其中：θ_e和θ_d分别为静态和动态接触角。

1.4 模型验证

先对本文所采用的方法进行验证。初始时刻设置一个半圆液滴置于固体水平面(θ=0)，在重力的作用下液滴会在固体平面上自由铺展，铺展程度用稳定时液滴高度与初始液滴高度比值H/H₀来衡量；引入Bo数作为无量纲参数，其表示重力和表面张力之比。此时，界面的形状由Bo数及稳定时的接触角θ_e决定。当Bo≪1时，重力作用几乎忽略不计，因此界面高度仅由式(14)决定：

当Bo≫1时，重力对界面高度产生影响，高度为

为方便研究，参照文献[13]的方法，采用静态接触角θ_e=135°的情况进行验证，并绘制了液滴形变H/H₀随Bo数变化曲线，如图 1(a)所示。图 1(b)为文献[13]的结果，通过比较发现结果基本符合公式，并与文献[13]的结果十分吻合(厄特沃什数Eo=Bo)。说明本文方法能够很好地描述重力场及表面张力作用下液滴在壁面上的运动过程。

1.5 初始条件

本文简化模型如图 2所示，计算域为虚线所含2×1的长方形区域，其中底面由均匀润湿区(长度L₁=0.5)和非均匀润湿区(长度L₂=1.5)共同组成，底面与水平面夹角γ=45°。平行于x轴方向的上下边界为广义滑移边界条件，垂直于x轴方向的左右边界为周期性边界条件。为了满足稳定性，指示函数的迭代精度为10^-4，压力迭代的精度为10^-3，时间步长取同时满足对流项和扩散项收敛的最小值。初始时刻，定义一个半径为R₀=0.1的半圆液滴置于均匀润湿部分中间，即L₁/2处。非均匀润湿部分在图 2中由黑白相间的方格表示，定义方格边长为W；白色区域表示亲水区，接触角大小与初始段一致，为θ₁，黑色区域表示疏水区，接触角为θ₂。

分别选取不同网格密度进行模拟，发现当网格密度大于256×512后，液滴位置随时间变化曲线基本重叠。因此，为保证计算结果与网格无关，又能节约计算成本，本文所有计算过程均在256×512网格密度下进行。此外，本文采用的所有物理量均已无量纲化。

本文涉及的无量纲参数：

式中：r为液滴的特征半径，本文中；ρ_l为液相密度。

2 结果与讨论 2.1 液滴在非均匀润湿壁面上的黏滑运动

Bo=0.352, γ=45°, θ₁=30°, θ₂=100°, μ_l=0.1, σ=0.1，疏水区宽度W=0.05时，液滴运动时前接触点位置及速度随时间变化曲线如图 3(a)所示。初始时液滴位于均匀亲水区，前接触点位于L=0.35处，在到达壁面缺陷前，即L=0.5之前充分发展至一定速度V_fcp=1.43后开始跨越非均匀润湿区。当液滴的前接触点到达第1个疏水区与亲水区的交界处时，速度开始迅速下降至接近为0，这是由于初始时前接触点位于亲水区，液固界面力小于气固界面力，使得壁面给予了接触点沿壁面前进的动力；而当前接触点到达疏水处时，液固及气固界面张力发生突变，最终体现为壁面突然给予前接触点一个沿壁面向上的阻力，因此前接触点的速度剧烈下降。而后由于重力作用依然大于前后接触角所产生的阻力作用，液滴仍作沿壁面下滑运动，前接触点速度又开始逐渐增加从而跨越疏水区域，待前接触点再次到达亲水区时，阻力又变为动力，前接触点的速度上升，随后将不断进行如此的周期性变化。如图 3(b)所示，后接触点的运动规律与前者相似，唯一不同点在于：当液滴位于亲水区域时，壁面给予后接触点的作用力与运动方向相反，成为液滴的阻力，而在疏水区为液滴的动力，并且后接触点的速度也随其所处位置呈周期性变化。最后液滴在非均匀润湿区域的运动表现为一种规律性的黏滑运动。接触点处的速度周期性变化特点及液滴运动过程并非在所有情况下都一致，下面根据不同的情形对液滴在非均匀润湿壁面上的运动过程进行研究。

2.2 黏滑运动受Bo数影响

重力是影响液滴在倾斜壁面上运动的一个重要因素，改变壁面的倾斜角度在某种意义上也是为了改变液滴运动方向上的重力作用。因此，本节通过仅改变重力来改变Bo数的大小，研究液滴运动随Bo数变化的规律。图 4为γ=45°, θ₁=30°, θ₂=100°, μ_l=0.1, σ=0.1, W=0.05, Bo=0.344~0.356时，液滴前接触点的位置随时间变化的曲线。可以发现，在前接触点到达非均匀润湿区之前，由于重力作用的时间较短，液滴运动过程差距不明显，位置随时间变化曲线几乎重叠；当其到达第1个疏水区时，液滴运动速度的差距开始显现但仍不明显；而液滴走到第1个亲水区时，液滴的运动受重力的影响开始变得明显，Bo数越大，液滴到达第1个亲水区域所需的时间越短，反之越长；当液滴在通过第3个疏水区时，不同Bo数下液滴的运动过程发生了较大的差异，Bo=0.344与Bo=0.348这2组液滴在跨越第3个疏水区域时失败，并稳定地停滞于第2个亲水区与第3个疏水区交界处，这是由于该情形下重力的作用无法克服由壁面润湿特性导致的阻力作用，在液滴前接触点通过第3个疏水区域的过程中，液滴的后接触点仍处于亲水区，前后接触点的阻力将重力所作功全部耗散，液滴停滞；而Bo=0.352与Bo=0.356两组液滴能够顺利地通过第3个疏水区域，并且在后续的运动过程中，虽液滴持续地进行黏滑运动，液滴的速度周期变化，但是整体趋势仍为加速运动，并且Bo数越大，液滴通过疏水区域所花费的时间越短，重力作用时间越长速度也越快。

2.3 黏滑运动受θ₁和θ₂影响

非均匀润湿壁面上，液滴发生黏滑运动是由于其受到因壁面周期性变化的润湿特性而带来的周期性变化的阻力作用而表现出来的运动形态。因此，研究液滴在非均匀润湿表面的黏滑运动，势必要将润湿差异的程度作为一个重要研究参数。定义液滴亲水接触角与疏水接触角之差为Δθ=θ₂-θ₁，Δθ代表润湿差异程度，Δθ越大，亲疏水壁面作用力差别越大。取Δθ=50°~90°这5种情况作为研究对象，具体的壁面接触角设置如图 5所示。图 5为Bo=0.360, σ=0.1, μ_l=0.1, γ=45°, W=0.01时，不同润湿差异下液滴前接触点随时间变化的曲线。初始时刻，亲水程度越高的表面上的液滴前接触点的移动速度越快，这是由于液滴的初始形态为半圆，受亲水壁面作用，液滴在下滑过程的同时也伴随着铺展过程，而壁面的亲水程度越高，液滴所受吸附力越大，铺展越快，因而前接触点的移动速度越快。而后，在跨越第1个疏水区域时，不同Δθ所体现的差别已经十分明显。30°~120°、35°~115°这2种情况在跨越第1个疏水区所耗费的时长远大于其他组，这是由于其前后接触角度差异较大分别为Δθ=90°、80°，此时当前后接触点分别位于疏水和亲水区时，由壁面润湿特性产生的阻力要大于Δθ=70°、60°、50°这3种情况，因此在跨越第2个疏水区时，30°~120°情况下液滴被阻碍无法前进，35°~115°在跨越第3个疏水条的过程中被阻碍无法继续前进；40°~110°情况下，液滴在前期能够快速跨越多个疏水区域，然而最终仍无法克服疏水表面阻力，因而在跨越第14个疏水区时停滞下来。不同于临界组40°~110°的是，当Δθ=60°、50°时，由壁面带来的阻力已经无法克服重力作用，随着时间的推移，前接触点速度变得越来越快。

2.4 黏滑运动受Oh数影响

液滴的流动性同样在液滴处于倾斜壁面上的运动过程起到关键作用，同时表面张力作用是决定液滴黏滑运动的主要因素，Oh数表征黏性力与表面张力相互关系，因此本节将Oh数作为研究对象。图 6为Bo=0.360，θ₁=40°，θ₂=110°，σ=0.1，γ=45°，W=0.01时，不同Oh数下液滴的前接触角位置随时间变化的曲线。可以看出，当液滴还未开始进入非均匀润湿区域前，液滴的运动速度随着Oh数的增加而减小，这是由于Oh数越大，流体受黏度的影响越大，流体的流动性越低，因此Oh数大的情况下液滴的移动速度会比Oh数小的时候低；当液滴开始进入非均匀润湿区后，情况发生反转，Oh数低的情况下，液滴的移动速度开始减小得较为明显，且Oh数越小，速度降低得越为明显，越难以顺利地通过非均匀润湿区，如Oh=0.631、0.673、0.715三组曲线所示，而Oh数大的情况下，液滴的速度开始追赶，并且保持整体做加速运动的趋势前进，且Oh数越大，加速趋势更为显著，通过非均匀润湿区的时间也越短，如Oh=0.757、0.799。这是由于Oh数越大，流体的黏性力与表面张力的比值越大，相较于表面张力的作用而言，黏性力在运动过程中占据了相对主导作用，接触点处不平衡的表面张力作用减弱，液滴沿壁面下滑运动过程受到壁面阻力的影响也就越来越小；反之，Oh数小，液滴受表面张力的作用就越为明显，受到的阻力作用也就越大，从而导致上述现象。

3 结论

本文采用FTM结合广义滑移边界，对液滴在倾斜45°的非均匀润湿壁面上的黏滑运动过程进行了数值模拟研究，研究过程通过改变Bo数、Oh数、亲水区接触角θ₁、疏水区接触角θ₂来观察液滴的前接触点的运动规律，研究结果如下：

1) 通过改变Bo数来改变液滴受重力影响，观察到当Bo数低于0.348时，液滴无法克服因壁面润湿特性的改变而带来的阻力，被壁面俘获。当Bo数大于0.352时，液滴在重力作用下能够通过非均匀润湿区域，并且Bo数越大，液滴通过疏水区域所花费的时间越短，且重力作用时间越长速度越快。

2) 通过改变θ₁与θ₂差值Δθ来表示非均匀润湿程度，从而观察其对液滴运动过程的影响。最终发现，当Δθ=70°、80°、90°时，壁面润湿差异带来的阻力作用大于重力，从而无法顺利通过非均匀润湿；而当Δθ=60°、50°时，由壁面带来的阻力已经无法克服重力作用，前接触点的速度变得越来越快，并且Δθ越小，阻力带来的效应越小，速度上升的也越快。

3) 通过改变Oh数来描述液滴黏滑运动受黏度和表面张力影响。发现当Oh数越大液滴受壁面阻力影响越小，黏滑运动带来的影响越不明显；反之液滴运动受壁面润湿特性变化带来的阻力影响越大，液滴黏滑运动越明显。当Oh数低于0.715时，液滴无法顺利通过非均匀润湿区域从而被壁面俘获，Oh数大于0.757时，液滴能够顺利地下滑，摆脱异质壁面。