刘乔 黄家宸 王昊 邓亚骏

1) (中国石化石油勘探开发研究院,北京 102206)

2) (北京大学工学院,北京 100871)

3) (深圳技术大学聚龙学院,深圳 518118)

1 引言

气-液-固三相接触线附近薄液膜是油气开采、相变换热和材料合成等系统的关键区域.在提高石油采收率技术中,润湿性接触角决定储层微纳米基质孔隙中毛细管力的大小,从而决定渗吸置换效率[1–3].热管中弯液面的蒸发和沸腾中气泡底部的蒸发过程中,此区域的薄液膜具有显著高于体相的传热效率[4,5].在天然气水合物的相变研究中发现,甲烷水合物优先在亲水石英表面的三相接触线区域生成[6].接触线作为液滴动态浸润的边界,此区域形成的动态接触角是模拟计算中的关键边界条件[5,7].目前虽然有众多关于接触线区域流体的实验研究,但只要涉及接触角的动态预测总会遇到困难,瓶颈在于对接触线区域薄液膜这一特殊结构认识不清[2,8–10].此区域具有从毫米到纳米的跨尺度特性,同时又受到多界面因素的共同影响,因此对其研究存在挑战.目前针对接触角有不同的计算模型[1,5,11,12],在流体动力学模型中,接触线处的摩擦耗散被认为可以忽略不计,动态微观接触角被认为可以采用静态接触角;而分子运动学模型考虑了移动接触线处局部耗散的重要性,认为动态微观接触角与速度相关[8].

对薄液膜微观结构和运移机制认识的缺乏是导致接触角普适预测模型没有形成的重要原因.不论是分子运动模型还是流体动力学模型,其关键争议点主要聚焦于微观接触角是否随速度变化、微观接触角与宏观静态接触角是否相同[7,13–16].因此,关键的一步是需要对接触线薄液膜区域微观接触角的动态过程进行监测.传统的测量手段是采用光学方法从液滴侧面对接触角区域液体进行观测,测量结果的分辨率取决于高速摄像机的放大倍数.光学方法获得的是杨氏方程所指出的宏观接触角,然而随着研究对象来到接触线的薄液膜区域,受界面分子间相互作用影响,杨氏方程无法预测微观接触角.Herminghaus 等[17]测量到己二醇液滴接触线区域20 nm 尺度的静态薄液膜形态.Yu 等[18]采用轻敲模式的原子力显微镜测量了丙三醇在云母和聚苯乙烯片上百纳米高度区域静态接触角,发现静态下液体微观接触角与宏观接触角相同.类似地,Deng 等[19]利用原子力显微镜测量了水滴在达到相变平衡状态下的100 nm 高度附近接触线形貌,获得了微观和宏观接触角一致的结论.那么,动态下接触线区域薄液膜微观形貌是否有所不同呢? Chen 等[20]针对移动接触线进行了实验探索,他们使用原子力显微镜成像系统揭示了前进接触线附近的纳米级形态,发现部分浸润前进接触线薄液膜呈现一个十几纳米的凸起结构,此凸起结构的存在使得移动接触线微观接触角不同于宏观接触角.进一步地,实验发现此介尺度凸起结构在部分浸润体系的前进接触线区域普遍存在,表明接触线处局部能量耗散的重要性[4,20,21].揭示接触线区域薄液膜结构与运移规律,对于理解接触线微观局部能量耗散以及进行接触角的预测具有重要意义.进一步地,需要对接触线区域的流体运动规律进行微观分析,作为纳米尺度动力学研究的重要手段,分子动力学模拟越来越多地被应用到微观流体的研究工作中.Yuan 和Zhao[16]研究了水滴在金固体基底上的铺展过程,发现接触线前端存在一层水分子网络,提出前驱液膜对于部分浸润体系同样具有普遍性.Blake 等[11]研究了有机物液滴在固体基底上的铺展细节,在此体系中没有揭示前驱液膜,也没有指出纳米凸起的存在.Lukyanov 等[15]模拟液体在平板间的运动认为动态接触角的变化受接触线区域高于体相区域的摩擦力影响.最近一项模拟工作表明,即使是纳米级的表面缺陷也可以对前进和后退的接触角产生不可忽略的影响[22].

本文通过构建大规模液滴的动态润湿体系,通过分子动力学模拟研究薄液膜的微观结构,获取不同润湿特性下接触线薄液膜结构的差异,揭示微观接触角与宏观接触角的演变规律.针对纳米区域的特殊结构,进一步采用粒子追踪的方法,探究不同接触线速度下薄液膜流体的运移机制.研究结果有望为构建普适的接触角预测模型提供关键物理机理.

2 模拟方法

如图1 所示,液滴润湿体系由固体基底和液滴组成.模拟过程分两个步骤.第1 步通过能量最小化过程使位于固体基底上方的液体达到平衡.对于初始构建的液滴,其中每个原子的初始给定位置对于整个液滴体系来说不是处于平衡状态,因此需要通过这一步使原子在相互作用下进行位置调整,逐渐达到平衡.在能量最小化过程之后,对液滴施加一定的速度使其接近固体基底.之后,外部施加的此速度被取消.第2 步即液滴铺展,液滴在液体和固体粒子间的相互作用下在固体基底上进行铺展.模拟采用三维直角坐标系,使用周期性边界条件,体系大小(x,y,z)=(300,6,240) nm.所有的模拟过程均由NVT 系综控制(即正则系综,指系统原子数量,系统体积和系统温度在模拟过程中保持恒定),系统温度控制在82 K,时间步长设置为2 fs.关于模拟体系的材料,选用目前势能参数最成熟、可靠性最佳的铂-氩体系,组成气态氩-液态氩-固态铂的三相系统[9].液滴是由相对原子质量为39.95 g/mol 的氩原子构建,直径为50 nm,包含319520 个氩原子.固体使用相对原子质量为195.11 g/mol 的铂原子构建.分子动力学计算的迭代过程通过LAMMPS 软件进行.

图1 液滴在固体基底表面铺展的分子动力学模型,液滴直径为50 nm.构建了两种固体基底,分别是绝对光滑固体基底A,表面无摩擦;另一种为晶格结构的固体基底B,单晶Pt(110)表面,原子尺度粗糙,表面有摩擦Fig.1.Molecular dynamics simulation model of droplet spreading on solid substrates,with a droplet diameter of 50 nm.Two kinds of solid substrates were constructed,one was an absolutely smooth solid substrate A with no friction on the surface;the other was lattice structure solid substrate B with friction,which has atomic scale roughness,single crystal Pt(110).

为了探究固体表面对于液滴铺展的影响,构建了不同结构的固体表面.首先构建绝对光滑的理想表面A,这是一个绝对光滑表面,没有任何缺陷和晶格间隙,在平行于表面的方向上没有摩擦.液态氩原子与表面A 之间的相互作用为

其中,ε 和σ 是的Lennard-Jones 参数.下标a和b分别表示不同种类的原子,它们可分别由l (液体)和s (固体)替代.同时,构建了具有晶格结构的表面B,即面心立方晶格结构的表面,晶格常数为0.392 nm.在B 表面上原子间相互作用为

对于氩原子间相互作用,εll和σll分别为0.2404 kcal/mol 和0.3405 nm.对于铂固体原子间相互作用,εss和σss分别为12.0181 kcal/mol和0.2475 nm[9].为了控制不同表面上液滴浸润性相同,表面A,B 上液-固原子间相互作用参数见表1.在探讨薄液膜液体分子的运移机制部分,研究表明接触线在高速运动时,气压对接触线区域气体裹挟的发生存在一定影响[23],因此在体系中加入氮气气氛,组成氮气+气态氩-液态氩-固态铂的三相系统,其中,氮气分子与铂、氩原子间相互参数见表2.

表1 不同固体表面粗糙度、液固原子间(液态氩、固态铂)相互作用参数和对应的体系平衡接触角Table 1.Surface roughness of different solids,interaction parameters between liquid-solid atoms(argon in liquid state,platinum in solid state) and corresponding equilibrium contact angles of systems.

表2 气液固原子间相互作用参数Table 2.Parameters of gas-liquid-solid atomic interaction.

采用均方根粗糙度对不同固体表面的粗糙情况进行描述,与Semal 等[24]接触线微观实验中原子力显微镜获取固体表面的均方根粗糙度原理一致.其中,hn为n点的表面高度,〈h〉为表面平均高度值,N为采样点总数.微观接触角与宏观接触角的计算与Santiso 等[25]采用的方法相一致,在模拟中提取液滴密度分布,在接触线区域做切线获取接触角,宏观接触角与微观接触角获取方式一致.

值得指出的是,本模拟体系具有两点优于以往分子动力学模拟体系的优势.1)构建了大规模分子模拟体系.以往的模拟体系在综合考虑模拟时间和体系合理性的基础上,普遍集中在十几纳米甚至几纳米范围.而纳米凸起结构在实验中被首次观测到的高度约为20 nm,因此在模拟中显著观测到凸起结构及微观接触角演变,需要采用大规模的模拟体系,从而避免小尺寸体系对结果的掩盖.2)采用不同类型固体表面来研究液滴前进接触线,实现了有无摩擦因素的模拟对比.

3 液滴前进接触线薄液膜结构

完全浸润液滴在固体B 表面上薄液膜结构对比如图2(a)所示.对于动态铺展过程中的完全浸润液滴,其前进接触线薄液膜为液滴主体平滑过渡至2—3 个分子层厚度的前驱液膜,前驱液膜覆盖了固体表面,改变了固体表面属性.该结果验证了前人对于完全浸润体系实验获得的结论[26].对于完全浸润体系,前驱液膜距离接触线处越近液膜越厚.本文所采用的大规模分子动力学模拟在揭示前驱液膜分子结构的同时,也证明了模拟方法的有效性.

图2 固体A,B 表面上不同浸润性液滴在1600 ps 时刻动态铺展结构对比,白色拟合虚线根据液滴球冠状理论所做,红色为接触线真实形貌,左上角分别对应各体系平衡态 (a)完全浸润液滴前进接触线区域存在数层分子厚度的前驱液膜;(b)部分浸润液滴前进接触线区域存在纳米凸起结构;(c)部分浸润液滴在绝对光滑固体A 表面前进接触线区域没有纳米凸起结构Fig.2.Comparison of dynamic spreading structures of liquid droplets with different wettability on the surface B at 1600 ps,the white fitted line was based on the droplet spherical crown theory,and the red was the real morphology of the contact line,the upper left corner corresponds to the equilibrium state of each system: (a) There was several molecular-layer thick precursor films in the vicinity of advancing contact line of complete wetting droplet;(b) there was convex nanobending in the vicinity of the advancing contact line of the partially wetting droplet;(c) no convex nanobending structure appeared in the vicinity of the advancing contact line on absolutely smooth surface A.

部分浸润液滴在固体B 基底上的铺展如图2(b)所示.从液滴整体来看,根据杨氏方程依据的液滴球冠状理论所做的圆形白色拟合虚线在接触线附近与实际形貌发生偏离,由此形成了接触线区域的纳米凸起结构,这是部分浸润体系与完全浸润体系不同的特殊结构.图中白色拟合虚线与固体基底在相交点M处形成宏观接触角,红色虚线与固体基底在相交点N处形成微观接触角.点M与点N之间的长度,来自于根据液滴球冠状理论所做的拟合虚线与实际形貌之间的偏离量,称为纳米凸起的长度.此模拟体系中的纳米凸起高度可达到18 nm,这也验证了采用大规模模拟体系能够成功获得此纳米结构的猜想.

对比部分浸润液滴在固体A,B 两种不同结构表面的动态铺展过程.A,B 两个系统唯一的不同是固体表面是否具有粗糙结构带来的水平摩擦.如图2(c)所示,不同于表面B 上纳米凸起结构的产生,在无摩擦的表面A 上,根据液滴球冠状理论所做的圆形白色拟合虚线在接触线附近与实际形貌发生向内的偏离,没有纳米凸起结构形成.此对比直接揭示了部分浸润液滴纳米凸起结构产生的原因是固-液界面动摩擦.如图3 所示,分析液滴在表面A,B 上的接触角变化全过程可知,液滴在理想无摩擦的表面A 上铺展时,始终保持接触角恒定,即等于静态接触角,无动态接触角变化现象.而在表面B 上,微观接触角随铺展过程而变化,始终大于静态接触角,当液滴达到平衡时纳米凸起消失,宏观接触角、微观接触角与静态接触角相等.

图3 部分浸润液滴在固体A,B 表面铺展过程接触角随时间变化.在绝对光滑A 表面,微观接触角始终保持与静态接触角相同(64°),不随时间发生变化.在单晶B 表面,宏观接触角与微观接触角均随时间发生变化,微观接触角大于宏观接触角,达到平衡时,接触角相同Fig.3.Contact angle changed with time during the spreading process of partial wetting droplets on the surface of solid A and B.On an absolutely smooth A surface,the microscopic contact angle was always the same as the static contact angle (64°) and did not change with time.On the surface of single crystal B,both the macro contact angle and the micro contact angle change with time,and the micro contact angle was larger than the macro contact angle,and the contact angles were similar to equilibrium angle when the equilibrium was reached.

4 薄液膜液体分子的运移

为了进一步探究纳米凸起结构的存在对于薄液膜液体分子运动的影响,从微观层面揭示其运动轨迹.采用粒子标记追踪的方法,通过对接触线区域液体分子的染色标记,追踪此部分液体分子的位置,并统计粒子最终分布数据.如图4 所示,液滴主体粒子标记为蓝色,接触线区域液体分子标记为红色,固体标记为褐色,氮气分子标记为粉色.模拟分为3 组,如图4(a)所示,液滴为自发铺展,定义从图中t0时刻到(t0+7320) ps 时刻液滴前进接触线铺展位移与时间之比为平均铺展速度,则此时液滴铺展过程的接触线平均速度约为0.007 nm/ps;如图4(b)所示,液滴通过施加恒定速度场,实现从t0时刻到(t0+1660) ps 时刻接触线平均速度约为0.03 nm/ps;如图4(c)所示,进一步提高液滴的铺展速度,实现从t0时刻到(t0+240) ps 时刻铺展过程中接触线平均速度约为0.21 nm/ps.整个铺展过程中,选取5 个主要代表时刻,对红色标记粒子的位置和分布进行追踪统计.

图4 液滴在固体B 表面不同速度下接触线区域红色粒子运动模式,液滴主体粒子标记为蓝色,接触线区域液体分子标记为红色,固体标记为褐色,氮气分子标记为粉色 (a)自发铺展液滴,从t0 时刻到(t0+7320) ps 时刻平均接触线移动速度0.007 nm/ps;(b)从t0 时刻到(t0+1660) ps 时刻平均接触线移动速度0.03 nm/ps;(c)从t0 时刻到(t0+240) ps 时刻平均接触线移动速度0.21 nm/ps.接触线区域液体分子的运动模式随速度提高,出现从滑动为主向滚动为主的转变模式,进一步发生气体入侵导致接触线滑移Fig.4.Droplet contact line particles movement mode at different velocities on the surface of solid B,the droplet body particles are marked in blue,the liquid molecules in the contact line area are marked in red,the solids are marked in brown,and the nitrogen molecules are marked in pink: (a) Spontaneous spreading droplet with an average contact line moving speed of 0.007 nm/ps from time t0 to time (t0+7320) ps;(b) average contact line moving speed is 0.03 nm/ps from time t0 to time (t0+1660) ps;(c) average contact line moving speed is 0.21 nm/ps from time t0 to time (t0+240) ps.The liquid sliding and rolling modes at the advancing contact line coexist,and tend to be the rolling mode as the speed increases.As the speed increases,air entrainment at the bottom of the contact line causes slippage.

在图4(a)中,对于自发铺展的平均速度约为0.007 nm/ps 液滴,红色标记粒子首先从液体主体右下侧运动至前进接触线区域,此时,粒子比较集中地聚集在前进接触线前沿区域,保持集中向前滑动,在接触线前进过程中,此部分粒子始终保持在前沿位置齐头并进,进一步追踪,待液滴铺展结束后,红色粒子水平移动至静态接触角的接触线区域,呈现为聚集在静态接触角附近,有一定的扩散.对比图4(b),平均速度约为0.03 nm/ps,红色粒子运动至接触线前沿区域,顶端粒子开始向下方运动至固体表面,并停留在固体基底表面,粒子不断同上方运动至固体表面呈滚动模式,待液滴主体铺展结束,红色标记粒子继续滚动前进,最终呈现为水平铺展在固体基底表面的状态.进一步提高平均速度至0.21 nm/ps,如图4(c)所示,气体侵入接触线下方,导致接触线前沿发生滑移,此时即发生溅射现象,即高速下气体入侵导致接触线处发生滑移.

为进一步证实液体分子运移模式的转变,对图4(a),(b)中铺展结束的红色标记粒子的水平分布进行统计.平行于固体基底方向,以2 nm 为单位,统计红色粒子的粒子数分布百分比.如图5(a)所示,对于接触线平均移动速度为0.007 nm/ps 的液滴,红色标记粒子的分布较为集中,百分比数据方差为34.72%.正如图4(a)所示,从中间时刻至铺展结束时刻,粒子主要以滑动的方式前进,最终粒子主要聚集在静态接触角附近区域.而对于接触线平均速度为0.03 nm/ps 的液滴,红色标记粒子的分布较为分散,每个单位间隔内的粒子数占粒子总数的百分比均低于13%,数据方差为15.21%.如图4(b)结果所示,从中间时刻至铺展结束时刻,接触线区域上方粒子不断向下滚动,粒子主要以滚动前进的方式运动,最终粒子水平铺开在固体基底表面.可以明确,从0.007 nm/ps 到0.03 nm/ps 的前进接触线区域液体分子随速度增大,逐渐从滑动为主转变为滚动前进为主.综上分析可知,在固体B表面由于摩擦的存在产生了接触线纳米凸起结构,使得接触线薄液膜液体随速度增大从滑动为主转变为滚动模式为主,进而接触线底部出现气体裹挟导致滑移、发生溅射.

图5 (a)平均接触线移动速度0.007 nm/ps 时,接触线处标注红色原子分布百分比方差为34.72%;(b)平均接触线移动速度0.03 nm/ps 时,接触线处标注红色原子分布百分比方差为15.21%Fig.5.(a) When the average contact line moving speed is 0.007 nm/ps,the variance of the distribution of atoms marked in red at the contact line is 34.72%;(b) when the average contact line moving speed is 0.03 nm/ps,the variance of the distribution of atoms marked in red at the contact line is 15.21%.

5 结论

本文通过大规模分子动力学模拟液滴接触线薄液膜结构与运移机制,结论如下.

1)通过模拟直径为50 nm 液滴的动态铺展过程,复现了完全浸润体系的前驱液膜结构和部分浸润体系的纳米凸起结构.

2)通过对比绝对光滑和实际固体表面动态润湿过程的不同,首次在模拟中获得了微观接触角与宏观接触角各自随时间的动态演变规律.

3)采用粒子示踪的方式对接触线区域液体分子进行追踪和统计,揭示了接触线薄液膜液体在固体表面摩擦作用下随速度增大从滑动为主转变为滚动模式为主,进而接触线底部出现气体裹挟导致滑移、发生溅射的机制.

研究结果有望为以下3 类方向提供理论指导:1)基于润湿性调控提高微纳米器件传热效率;2)基于润湿性调控提高页岩油储层基质渗吸置换效率;3)构建普适的接触角预测模型.