微通道疏水表面滑移的耗散粒子动力学研究*

许少锋 楼应侯 吴尧锋 王向垟 何平

(浙江大学宁波理工学院,宁波 315000)

1 引 言

近年来,疏水表面在流动减阻、自清洁、防污等方面显示了广泛的应用前景[1,2].普遍认为,流体流经疏水表面时会发生滑移[3-5],尤其是在微观尺度下,滑移更加明显.经典流体力学中的无滑移假设适用于多数的宏观黏性流体流动问题,但对于疏水表面,特别是在微纳尺度下研究流体流动问题时,微小的滑移也会给流体的输运产生重要的影响,需要利用滑移边界才能准确获得流体的流动行为.

随着微纳机电系统及微流控芯片的快速发展,疏水表面的滑移研究受到越来越多学者的关注.疏水表面的滑移已得到大量实验研究的支持,如Choi等[6]以及Lee和Kim[7]利用流变仪研究了流体在疏水表面的滑移情况,测得在剪切速率为105s¯1时的滑移长度接近30 nm；Tretheway 和Meinhart[8]采用粒子图像测速技术(particle image velocimetry,PIV)测得含疏水涂层微管壁面的滑移长度为0.92 μm；Bhushan等[9]用原子力显微镜测得疏水和超疏水表面的滑移长度分别为43 和236 nm；Pit等[10]利用他们自己设计的实验装置直接观察到了疏水表面的滑移.尽管各种新的实验技术被用于疏水表面的滑移研究,但由于实验材料各异,实验条件也难以控制,而且微观尺度下的实验对精度要求苛刻,因此实验研究结果还难以揭示疏水表面的滑移规律及其机理.而计算机数值模拟可以直接研究系统的微观细节,能考虑到实验研究和理论分析所忽略的细节,这使得数值模拟在研究微尺度流动方面具有明显的优势.很多学者采用分子动力学方法、格子Boltzmann方法与耗散粒子动力学方法等数值模拟技术开展疏水表面的滑移研究,试图揭示其滑移规律(如滑移长度与接触角之间的关系)及其机理.曹炳阳等[11]采用分子动力学方法模拟了液态氩在铂纳米通道内的流动,通过改变壁面与流体之间的势能作用实现壁面的不同疏水性能,研究发现在疏水表面附近形成了一个低密度层,低密度层使流体在壁面发生了滑移.一般认为壁面接触角越大,滑移长度越大,但Voronov等[12,13]的分子动力学模拟得到了不同的结论,他们模拟了类石墨烯平板间LJ (Lennard-Jones)流体的Couette流,发现接触角在140°左右时滑移长度最大,并不是接触角越大滑移长度越大.黄桥高等[14]利用格子Boltzmann方法研究了疏水表面微通道内的流体流动,通过改变吸附参数(反映壁面对流体粒子的作用强度)得到疏水表面,结果表明疏水作用在近壁区诱导了一个低密度层,表观滑移则发生在低密度层,并发现滑移长度与吸附参数间存在近似的指数函数关系.Zhang等[15]也采用格子Boltzmann方法得到了与文献[14]类似的结果,滑移长度与接触角之间满足指数函数关系.Cupelli等[16]利用耗散粒子动力学研究了毛细管的润湿行为,改变固体壁面与流体粒子的相互作用参数得到了不同润湿性能的壁面,模拟结果显示接触角与固液相互作用参数呈近似线性关系.关于疏水表面滑移机理,Tretheway和Meinhart[17]提出的理论认为在疏水表面存在不连续的气体层或纳米气泡,这些不连续的气体区域是表面产生滑移的原因.

从以上研究可以看到,疏水表面会产生滑移已获得一致认同,但对其滑移机理以及滑移规律还没有得到一致的结论,有待进一步研究.本文采用耗散粒子动力学(DPD)方法模拟微通道Couette流动,通过改变壁面粒子与流体粒子之间排斥力系数获得不同疏水性能的壁面,探讨疏水表面的滑移问题,以期揭示疏水表面的滑移规律及其机理.

2 DPD方法与壁面边界模型

2.1 DPD方法

DPD是由Hoogerbrugge和Koelman[18]于1992年提出的一种介观尺度的模拟方法,后来Espanol和Warren[19]以及Marsh[20]建立了DPD统计力学理论模型,奠定了其理论基础,Groot和Warren[21]的研究工作进一步推广了DPD的应用.DPD模型中,每个粒子是由一团分子组成的,是一种粗粒化的粒子,粒子之间采用软势能作用,因此DPD模型模拟的时间和空间尺度比分子动力学大得多.DPD粒子之间具有保守力、耗散力和随机力三种形式的作用力,三种作用力都是成对出现的,可以确保体系动量守恒,十分适合模拟流体动力学问题.DPD方法已被广泛应用于液滴、多相流体、高分子溶液、血液等复杂流体的动力学模拟[22-24].

DPD方法中,粒子的运动规律通过牛顿第二定律得到.具体地说,对i粒子,其运动方程为

式中ri,vi分别为粒子的位置矢量和速度矢量.注意(1)式中,每个粒子的质量相等,并且质量采用无量纲单位1.Fext为施加的外场力.Fij是粒子j对粒子i的作用力,由保守力,耗散力,随机力三部分组成：

粒子i和粒子j之间只在截断半径rc之内具有相互作用力,两粒子之间的距离大于rc时则没有相互作用力,三种作用力均沿粒子之间连线方向,都在rc范围内才有相互作用.

式中aij为粒子i和j之间的最大排斥力系数,rij=ri-rj,rij=|rij|,eij=rij/|rij| .耗散力是依赖于粒子位置和速度的力,亦称耗散阻力,其表达式为

式中γ为耗散力系数,wD(rij)为耗散力权重函数,vij=vi-vj.随机力的表达式为

式中σ为随机力系数,wR(rij)为随机力权重函数,ζij是均值为0方差为1的随机数,θij为高斯分布的随机数,具有对称性θij=θji,且具有以下随机特性,对i≠j,k≠l：

式中 〈···〉 表示求平均.

耗散力阻止粒子间相对运动,使粒子速度减慢,其宏观效应是使系统温度降低,而随机力宏观效应是使系统升温.这样,在耗散力和随机力的协同作用下,系统能保持温度恒定,即引入的耗散力和随机力相当于分子动力学模拟中的恒温热浴.

DPD模拟中,各物理量均采用无量纲单位,通常选取rc为特征长度,密度ρ为特征密度,kBT为特征能量,其他物理量单位由这个三个参数推导得到.在DPD模拟中,通常取无量纲截断半径rc=1,无量纲温度kBT=1,粒子无量纲质量m=1,无量纲密度ρ= 3 — 6.通过定义DPD粒子的质量、平均热能和截断半径,可以得到对应的实际物理单位.

2.2 壁面边界模型

由于耗散粒子动力学中粒子之间作用为软势能的形式,固体壁面粒子的软的排斥力不足以阻止流体粒子穿透壁面,导致施加固-液界面边界条件比较困难.此外,壁面边界结构处理不合理还会造成壁面附近一些假象的波动,如壁面附近流体密度和温度波动.DPD固体壁面边界条件是DPD模拟中的研究重点和热点问题,过去20年很多学者提出了不同的施加固体壁面边界条件的方法[25-31],其中一个最常见的方法是采用固定住的粒子并结合适当的反射机制构建固体壁面边界,如Revenga等[25]提出了三种反射机制,即镜像反射、向后反弹发射与Maxwellian反射,将穿过壁面的流体粒子重新移回流体区域,刘谋斌和常建忠[31]通过冻结随机分布并且达到平衡状态的DPD粒子,形成复杂的固体区域,结合反弹反射,该方法能有效地模拟复杂固体区域.本文模拟中采用固定住的粒子构建固体壁面,并通过一定的反射机制将穿透壁面的流体粒子重新移回流体区域.图1为固体壁面结构示意图,图中黑色圆圈代表壁面粒子,空白圆圈代表流体粒子,空白圆圈上的箭头表示粒子的速度方向.如图1所示,采用三层固定住的粒子构建固体壁面,由Duong-Hong等[28]的研究工作可知,合理地调整这几层固体壁面粒子与固-液界面之间的距离,可以消除壁面附近的密度波动并能阻止大多数流体粒子穿过固-液界面.本文模型中,三层固定住的粒子离固-液界面的距离分别为0.15rc,0.35rc,0.75rc,第一层壁面粒子离固-液界面的距离设置得比较小,能适当地增大壁面的排斥力,有助于阻止大部分流体粒子穿过固-液界面,另外两层壁面粒子主要是为了给壁面附近的流体粒子提供均匀的壁面排斥力,进而消除壁面附近流体密度波动.相邻粒子间的距离在x和y方向分别为0.5rc和rc.壁面粒子与壁面附近的流体粒子之间也受到(2)式中的三种作用力.

此外,当流体粒子穿过固-液界面时,采用如图2所示的修正过的向前反弹机制将穿透壁面的流体粒子重新移回流体区域.如图2所示,t时刻粒子i处于位置1处,在t+Δt时刻粒子i穿过了固-液界面运动到位置2处,此时将粒子i重新移回到流体区域的位置3处,位置2与位置3关于固-液界面对称,并且移回到位置3处的流体粒子速度方向与位置2处的速度方向相反(图中从粒子中心引出的箭头线表示速度分量).传统的向前反弹机制是移回到位置3处的流体粒子速度方向与位置1处的速度方向相反,也即将t+Δt时刻的速度与t时刻的速度反向,两个时刻的速度大小相同,而本文改进的反弹机制是将t+Δt时刻位置2处速度反向,得到位置3处的速度,这更符合物理运动规律.

图1 固体壁面结构示意图Fig.1.Sketch of the solid wall structure.

图2 修正的向前反弹机制Fig.2.Sketch of the modified bounce-forward reflection.

2.3 模拟系统与模拟参数

本文以平板间的Couette流为例,探讨微通道疏水表面的滑移问题.图3所示为微通道(两平板间)Couette流模拟系统示意图,x和y方向采用周期性边界条件,z方向采用上述壁面边界模型.模拟区域的实际尺寸在x,y,z方向为Lx×Ly×Lz=60rc×5rc×30rc.DPD单位中,截断半径rc的无量纲单位为1,则模拟区域无量纲尺寸为Lx×Ly×Lz=60×5×30.笛卡尔坐标系的原点设在模拟区域的中心,上下固-液界面分别位于z=15与z=-15 处.所有流体粒子和壁面粒子的质量均相等,每个流体粒子或壁面粒子的质量设为1(m=1 ).所有的壁面粒子按照上述壁面边界模型布置在上下壁面区域,则一共含有3600个壁面粒子,上下壁面各1800个粒子.所有流体粒子的初始位置按照面心立方体晶格布置在模拟区域,流体的密度取ρ=4.0,则一个单位立方体内包含4个流体粒子,流体粒子的总数为36000.所有的模拟中,系统的温度取kBT=1,流体粒子的初始速度根据此温度随机设定.取耗散力系数γ=4.5,随机力系数为σ=3.0 .

图3 模拟系统示意图Fig.3.Sketch of simulation system of Couette flow.

上下壁面沿x方向分别施加等值反向的速度v和—v,用于驱动形成Couette流.采用Groot和Warren[21]提出的修正Velocity-Verlet算法数值求解运动方程(1)—(9)式,获得流体粒子的运动轨迹.数值求解的时间步长取 Δt=0.02 .为减少计算时间,采用文献[32]描述的元胞分割法和临位列表法计算粒子之间的作用力.

2.4 壁面疏水表征

本文模拟系统包含两种粒子,即壁面粒子和流体粒子,壁面粒子间排斥力系数为aww,流体粒子间的排斥力系数为aff,其中下角标‘w’表示壁面粒子,‘f’表示流体粒子,根据文献[27]定义壁面粒子与流体粒子之间的排斥力系数：

由(3)式可知,排斥力系数代表相互作用粒子间的最大排斥力,即awf为壁面粒子与流体粒子间的最大排斥力.因此,awf的大小可以反映出壁面与流体的相互作用强度,通过调整awf的大小即可实现壁面不同的亲/疏水性能,awf越大,壁面与流体的排斥作用越强,疏水性能越强.文献[16]和文献[33]就是通过调整壁面粒子与流体粒子之间的排斥力系数awf,从而实现壁面的不同润湿性能,文献[16]研究结果表明,排斥力系数awf与接触角呈近似线性关系,壁面与流体的排斥作用越强,流体的接触角越大,调整awf大小可以实现从亲水到疏水的转变.本文模拟中取流体粒子间的排斥力系数为aff=18.75,为使壁面具有不同的亲/疏水性能,实现无滑移到滑移的变化,设计了不同壁面粒子间的排斥力系数aww=5,10,15,20,25,对应的壁面粒子与流体粒子之间的排斥力系数为awf=9.68,13.69,16.77,19.36,21.65.

3 微通道疏水表面滑移

3.1 壁面无滑移

在研究壁面滑移之前,先讨论壁面无滑移的情况.一方面有助于与壁面滑移进行比较,另一方面可以验证上述壁面边界模型的有效性以及编制的DPD程序的正确性.本文研究的表面滑移均为部分滑移(partial slip),图4为壁面无滑移和部分滑移示意图(壁面是静止不动的).部分滑移边界条件可以用Navier滑移边界模型描述如下：

式中us为滑移速度,b为滑移长度,ux为x方向的速度,滑移速度和滑移长度的定义如图4所示.

图4 壁面边界的无滑移和部分滑移示意图Fig.4.No-slip and partial slip status at a solid-fluid interface.

在壁面无滑移模拟中,取壁面粒子间的排斥力系数aww=5,对应的壁面粒子与流体粒子间的排斥力系数为awf=9.68 .模拟先运行一段时间,使系统达到热力学平衡状态,然后给上壁面施加沿x正方向的速度=3.0,同时给下壁面施加沿x负方向的速度=-3.0,驱动形成Couette流动.模拟再运行一段时间,使系统达到稳定的Couette流后开始统计相关物理量.为获得速度、密度、温度以及剪切应力在通道z方向上的分布,通道在z方向上被划分为200层,每 2×104个时间步长对每层的数据进行统计平均得到相关物理量.

图5 DPD模拟的速度分布与Navier-Stokes(NS)分析解对比Fig.5.The velocity profile of DPD simulation result compares to the analytical solution.

图6 DPD模拟的密度、温度分布与Navier-Stokes(NS)分析解对比Fig.6.The density and temperature profiles of DPD simulation results compare to the analytical solutions.

图7 DPD模拟的剪切应力分布与分析解对比Fig.7.The shear stress profile of DPD simulation result compares to the analytical solution.

图5为DPD模拟的Couette流的速度分布与NS方程得到的分析解对比,图6为DPD模拟得到的密度、温度曲线与NS分析解比较,图7为模拟的剪切应力分布与分析解对比.从图5—图7可以看到,DPD模拟得到速度、密度、温度、剪切应力分布与理论分析解符合得很好,验证了本文编写的DPD程序的正确性.更重要的是,从图5可知壁面处流体的速度与壁面的速度相同,即壁面无滑移.从图6可以看到,模拟得到的密度与温度分布比较均匀,密度和温度曲线与理论曲线几乎重合,在壁面附近密度、温度几乎没有波动,图7显示模拟得到的剪切应力分布与分析解(Szx=µ∂vx/∂z,µ为动力黏度)符合得也很好,验证了上述边界模型的有效性.

由此可见,壁面粒子与流体粒子间的排斥力系数为awf=9.68 时,壁面实现了无滑移边界条件,壁面附近密度分布均匀无波动.由模拟结果也可以看到,模拟得到的速度、温度、密度以及剪切应力等物理量的分布与Navier-Stokes理论分析解符合得很好,这也说明了本文固体壁面边界模型以及DPD模型具有很好的收敛性.

3.2 疏水表面滑移

由前面分析知道,壁面粒子与流体粒子间的排斥力系数为awf=9.68 时壁面为无滑移壁面,这也能反映出壁面为亲水性壁面.现增大awf,即增强壁面与流体的排斥作用,使壁面实现从亲水到疏水的转变,分析疏水表面的滑移机理,讨论滑移长度与awf的关系.

图8为壁面粒子与流体粒子间的排斥力系数awf分别为9.68,13.69,16.77,19.36和21.65时的速度分布.由于速度分布是关于坐标原点对称的,为了显示更清晰些,图8只显示了速度分布的一半,即z=-15—0 之间的速度曲线.awf=9.68 时,壁面满足无滑移条件,而当增大awf,也即增强壁面与流体之间的排斥作用,使壁面向疏水转变,速度分布发生了变化,也即当awf＞9.68 时,壁面产生了滑移.根据Navier滑移边界的定义,将主流区的速度曲线线性延伸至壁面,得到的速度若小于壁面的速度则为滑移边界.由图8可以看到,awf=13.69,16.77,19.36,21.65时,壁面均产生了滑移,并且随着awf的增大,即疏水性增强,壁面滑移越明显.

图8 排斥力系数 awf分别为9.68,13.69,16.77,19.36,21.65时的速度分布Fig.8.Velocity distributions due to different values of awf=9.68,13.69,16.77,19.36,21.65.

为探索排斥力系数awf(疏水强度)对壁面滑移的影响,下面分析滑移速度和滑移长度与排斥力系数awf的关系.根据Navier滑移边界的定义,滑移速度和滑移长度可以通过如下方法获得：对主流区的速度线性拟合得到直线形的速度分布,拟合的速度直线延伸至壁面处时得到速度与壁面速度的差值即为滑移速度,拟合的速度直线向外延伸出实际流场区域直至满足无滑移边界条件时离壁面的距离即为滑移长度.表1为awf= 13.69,16.77,19.36,21.65时通过线性拟合得到的速度分布以及对应的滑移速度和滑移长度,表中还列出无滑移时的情况(滑移速度和滑移长度均为0).根据表1中的数据,得到了图9所示滑移速度us与排斥力系数为awf的关系,以及图10所示滑移长度b与排斥力系数awf的关系.由图9和图10可以看到,随着壁面粒子与流体粒子间的排斥力系数awf的增大,即壁面与流体之间排斥作用增强(疏水性增强),滑移速度和滑移长度均增大.分别对滑移速度us和排斥力系数awf、滑移长度b和排斥力系数awf进行曲线拟合,如图9和图10所示,拟合结果显示滑移速度us或滑移长度b与排斥力系数awf之间存在近似的二次函数关系,函数表达式分别为：

表1 不同排斥力系数 awf时拟合的速度分布及对应的滑移速度和滑移长度Table 1.The linear fit velocity profiles,slip velocity and slip length with respect to different awf.

图9 滑移速度 us与排斥力系数 awf的关系Fig.9.Slip velocity versus awf.

图10 滑移长度b与排斥力系数 awf的关系Fig.10.Slip length versus awf.

表面的亲/疏水性一般是通过接触角表征的,表面疏水性越强,接触角越大,表面滑移也越明显,但目前关于滑移长度和接触角之间还没有明确的关系式.根据文献[16]的DPD模拟结果可以看到,壁面与流体之间作用参数awf与接触角呈近似线性关系,结合本文研究结果(13)式,可以认为滑移长度与接触角之间存在近似的二次函数关系.接触角是衡量液体对固体材料表面润湿性能的重要参数,反映了液体与固体表面亲和作用的强弱,也即接触角的大小取决于液体与固体表面的相互作用强度.而疏水表面的滑移本质上是由液体与固体表面相互作用引起的,滑移长度的大小同样取决于液体与固体表面的相互作用强度.因此,滑移长度与接触角之间必然存在着某种函数对应关系.本文模拟中,壁面粒子与流体粒子之间的排斥力系数awf反映了流体与固体壁面的相互作用强度,由(13)式可以看到,滑移长度与接触角之间满足近似二次函数关系.

目前关于疏水表面滑移现象的物理机理还没有定论,比较常见的两种解释为：一是不连续气体层或纳米气泡理论,认为在表面存在不连续气体区域,因此固体表面同时存在液-固界面和液-气界面,液-气界面流体流动阻力小,产生了滑移；二是低密度层理论,由于排斥作用,疏水表面附近存在低密度区域,阻止了动量的传递,壁面发生了滑移.为探讨疏水表面的滑移机理,图11给出了排斥力系数awf=13.69,16.77,19.36,21.65时的密度分布和无滑移时的密度分布,由于远离壁面区域的密度几乎一样,并且密度分布关于z=0 对称,为了显示更清晰,图11只显示了下壁面z=-15 至z=-10 之间的密度曲线.无滑移时理论密度分布是均匀的,当awf=13.69,16.77,19.36,21.65 时,即壁面产生滑移时,在紧邻壁面附近流体密度很低,存在低密度区域,这些低密度区域阻碍了流体的动量传递,致使壁面产生滑移.从图11可以看到,在紧邻壁面处密度低,即存在一个低密度区域,紧邻低密度区域又存在一个密度大于4.0的高密度区域.这是由于壁面附近存在低密度区域,导致高密度区域处的流体在垂直于壁面方向上受到的排斥力不平衡,进而出现了高密度区域.根据文献[34,35]的研究可以知道,在DPD方法中,流体粒子受到的排斥力在垂直于壁面方向上的分力控制着密度的波动.本文中,在紧邻壁面区域存在低密度区域,密度远小于4.0,即流体粒子数目较少,而在远离壁面区域密度分布均匀,密度为4.0.因此,低密度区域流体对高密度区域流体的排斥力小于密度均匀区域流体对高密度区域流体的排斥力,也即高密度区域流体粒子受到朝着壁面方向上的净排斥力,该净排斥力使流体粒子向低密度区域方向运动,流体粒子往该方向运动,使其与低密度区域一侧流体粒子距离减小,由(3)式可知,粒子之间的距离减小,相互之间的排斥力增大,因此净排斥力使流体粒子往低密度区域一侧逐渐靠近,低密度区域一侧流体粒子对该流体粒子的排斥力逐渐增大,必然存在一个平衡位置,当流体粒子运动到该平衡位置时,增大的排斥力与净排斥力平衡,流体粒子会在该平衡位置堆积,故在紧邻壁面低密度区域又会存在高密度区域.由图11还可以看到,排斥力系数awf越大,即疏水性越强,紧邻壁面附近的密度越低,阻碍动量传递的能力越强,滑移越明显.

图11 不同排斥力系数时和无滑移时的密度分布Fig.11.The density profiles for different awfvalues and for no slip condition.

4 结 论

本文利用耗散粒子动力学模拟了平板间的Couette流,研究了微通道疏水表面的滑移现象.采用三层固定住的粒子并结合修正的向前反弹机制构建了固体壁面边界模型,通过调整壁面粒子与流体粒子之间的排斥力系数awf大小,壁面可以实现从亲水到疏水的转变,实现从无滑移到滑移的变化.当awf=9.68 时,壁面满足无滑移边界条件,壁面附近密度分布均匀；当awf＞9.68 时,壁面产生了滑移,并随着排斥力系数awf增大(疏水性增强),壁面滑移越明显,滑移速度和滑移长度越大.模拟结果还表明滑移速度和滑移长度与排斥力系数awf存在近似二次函数关系,也即与接触角之间存在近似二次函数关系.模拟得到的密度分布显示,无滑移时壁面附近的密度分布均匀,而疏水表面附近存在低密度区域,低密度区域阻碍了流体的动量传递,导致壁面产生了滑移.