引言

纳米流动系统因其高效、经济等优势在医学、制药、国防以及环境监测等领域具有广泛的应用前景，如DNA分析和测序，药物输送，化学分析等^[1-3].在纳米流动中，尺度的减小导致固液界面滑移现象对流动的影响更显著^[4-7].因而通过调控纳米流动系统中壁面润湿性将能够产生流动滑移(flow slip)和可观减阻(significant drag reduction)^[8-9]，进而提高系统效率.为描述纳米通道中流体的流动特性，需要全面掌握固液界面的滑移特性和产生机理.

分子动力学(MD)能以原子级精度描述流体在紧邻壁面处的流动特性，被广泛应用于纳米流动的研究.现有MD研究表明，滑移依赖于固液润湿特性，流动剪切率，固液密度比以及系统的热力学条件. 1997年，Thompson和Troian ^[10]提出了一种普适滑移模型，定量描述了在不同固液相互作用强度以及壁面密度下，滑移长度(slip length)与剪切率之间统一的依赖关系. 2004---2011年，Priezjev等^{[2, 11-14]}在Thompson等的工作基础上，开展一系列关于简单和聚合物液体滑移特性的研究，给出了通道表面热粗糙度、几何粗糙度、分子链长对滑移的剪切依赖特性的影响规律.

调控固液界面润湿特性能方便地控制纳米流动的滑移量，因此，受到研究者的广泛关注. 1990至2015年，大量的研究^[15-20]均发现, 滑移长度随固液润湿性的降低而增大. 1999年，Barrat等^[15]基于滑移长度的Kubo公式和MD模拟给出了滑移长度与液体在固体表面静态特征量之间的定量关系. 2008年，Huang等^[19]利用MD模拟方法，以液态水为流动介质，给出了滑移与接触角的准普适关系： $b \propto (\cos \theta_{\rm c}+ 1)^{-2}$ ，其中 $b$ 为滑移长度， $\theta_{\rm c}$ 为水在壁面上润湿时的接触角. 2011年，Priezjev^[14]研究了周期性高低表面能相间壁面上Couette流动，给出了高表面能条带与流动方向之间的夹角对滑移特性的影响规律.同年，Wang等^[21]提出了滑移的分子动理论(molecular kinetics)模型，并给出从无滑移到有滑移的临界剪切率对固液相互作用强度的指数依赖关系.

上述研究主要集中于均质壁面上纳米流动的滑移特性.然而在实际化工环境中，材料在制备中经常会混入少量杂质或被轻微污染.因此，需开展非均质低表面能壁面上流动滑移特性的研究.本论文以Poiseuille流动为对象，采用非平衡MD模拟方法，研究了低表面能壁面表层混入少量亲液杂质对壁面滑移特性的影响.与Priezjev ^[14]的研究不同，这里重点关注亲液杂质占比及其分布形式(集中式和均布式)对流动滑移特性的影响规律，并基于分子动理论讨论了其产生机理.

1.   分子动力学模拟方法

前人研究表明，二维MD模拟不仅在简单流动中可反映流动的物理本质^[22]，而且对于热对流、钝体绕流等复杂流动中也能获得与理论与实验相吻合的结果^[23-24].因此，这里采用基于降维思想^[22]的二维MD模拟方法来模拟两无限大平行平板间充分发展Poiseuille流动.在MD模型中，液体原子之间的相互作用采用经典的Lennard-Jones(LJ)12/6势能函数描述.该势能模型被广泛用于纳米流动的MD模拟，其表达式为

其中， $\varepsilon $ 和σ 分别为液体原子特征能量和特征长度， $r$ 为两个液体原子之间距离.同时，为提高计算效率，设置LJ势能的截断半径 $r_{\rm c} =2.5 \sigma $ ，即当两个原子间距离大于截断半径时, 认为这两个原子之间的势能为零.

数值模拟的几何模型如图 1所示，其中，空心菱形为疏液原子，实心菱形为亲液原子，两壁面之间为液体原子.通道长度 $L=50.709 2\sigma $ ，上下壁面厚度 $D =3.457 4 \sigma $ ，通道宽度 $B =17.439 7 \sigma $ .上下壁面均由312个原子构成面心立方结构，晶格常数为 $1.152 5 \sigma $ ；液体原子共512个.固液之间相互作用也由LJ12/6模型描述，其特征长度 $\sigma_{\rm SL}=0.867 3\sigma $ .亲液杂质原子在壁面表层上集中和均布形式如图 1(a)和图 1 (b)所示.均布亲液杂质以单原子方式在壁面表层上均匀分散排列，该方式与纳米尺度电渗流(electroosmotic flow)MD模拟中在壁面均匀排布带电原子^[25]方法相同.亲液杂质原子占比α定义为亲液原子个数与其所在固体原子层总原子个数之比，图 1中 $\alpha \approx 0.24$ .

图  1  含亲液杂质原子的纳米通道流动系统示意图，图中右侧为纵向坐标，单位为σ

Figure  1.  Schematic diagram of flow in nanochannel with wetting impurity atoms. The vertical coordinate axes are shown in right side of figures and their units are σ

下载: 全尺寸图片 幻灯片

Gennes ^[26]基于杨氏方程给出液滴在固体表面上接触角θ与固液间势能特征能量 $\varepsilon_{\rm SL}$ 之间的定量关系

其中， $\varepsilon $ 为液体原子间势能特征能量.曹炳阳等^[27]已利用MD模拟方法验证了上述关系.因此，文中通过调整固液原子之间LJ势能的特征能量 $\varepsilon_{\rm SL}$ 来调控液体在固体表面的润湿性.设置疏液壁面与液体原子之间 $\varepsilon _{\rm SL}=0.25\varepsilon $ ，对应固液界面接触角为120 $^\circ$ ；亲液杂质原子与液体原子之间 $\varepsilon_{\rm SL}=1.0\varepsilon $ ，对应接触角为0 $^\circ$ .

模拟过程中，采用Nosé-Hoover温度控制方法将系统温度始终控制为0.8 $\varepsilon /k_{\rm B}$ .待系统从初始状态达到热力学平衡后，对每个液体原子施加一个大小为 $0.005 \varepsilon / \sigma $ 的驱动力，以模拟Poiseuille流动.另外，前人研究已证实, 上述量级的驱动力所产生的流动处于低剪切率区内^[28-30]，此时刚性壁面和柔性壁面产生的流动滑移和流场特性基本相同^[30].因此，为减小计算量，采用刚性壁面模型，即固体原子在晶格点处始终保持不动.文中所有算例均采用LAMMPS开源代码，计算时间步长 $\Delta t =0.002 \tau $ ，其中， $\tau =( m \sigma^{2}/\varepsilon)^{1/2}$ .施加驱动力后，系统继续迭代5×10 $^{5}$ 步使流动达到平衡状态，同时将通道沿垂直壁面方向分为872层，每层厚度为0.02σ .最后通过统计平均1×10 $^{6}$ 个时间步内每层内速度、数密度等物理量以获得通道内沿流向平均的流场信息.

文中采用滑移长度的概念表征滑移量的大小.滑移速度和滑移长度的计算方法为：通过抛物线方程拟合主流区的速度分布得到表达式 $v_{x}(y)$ ；基于此表达式得到主流区速度分布延伸至壁面处的速度 $v_{x}\vert _{\rm wall}$ 和壁面剪切率 $\gamma $ ；滑移速度 $v_{\rm S}=U-v_{x} \vert_{\rm wall}$ ( $U$ 为壁面速度，文中壁面静止)，而滑移长度 $L_{\rm S}=v_{\rm S}/\gamma $ .其中，壁面位置设定为靠近液体的第一层固体原子的位置处.上述滑移模型称为Navier模型，被广泛应用于界面滑移的实验和数值模拟研究中.由该模型计算得到的滑移量被称为表观滑移(effective slip length)^[29].

2.   模拟结果与讨论

2.1   液体的密度分布

亲液杂质均布时，不同亲液杂质占比(α)条件下液体的法向密度分布如图 2(a)所示.其中，横坐标表示在垂直固体壁面方向上的位置，坐标原点与通道中心重合；纵坐标为液体原子的数密度，即单位体积内液体原子的个数.由于密度分布沿通道中心线对称，图 2仅展示下壁面至通道中心的密度分布.下壁面靠近液体的第一层固体原子位于 $y=-9.219 9\sigma $ .

图  2  亲液杂质占比对通道内流体密度分布的影响

Figure  2.  The influence of wetting impurity on fluid density

下载: 全尺寸图片 幻灯片

不同α时，液体密度分布均先出现一个较高的峰值点，然后在约3~4个液体原子直径范围内振荡衰减至体相密度值.且随α增大，壁面对液体原子的平均势能约束增强，致使液体密度振荡程度逐渐增加，但振荡周期和相位并未变化.液体密度分布的第一个极大值点称为液体在固体表面的接触密度( $\rho_{\rm C}$ )^[28]，一般认为其与滑移密切相关.因此，本文考察了 $\rho_{\rm C}$ 对α的依赖关系(图 2(b)).从中发现, $\rho_{\rm C}$ 随α的增大量线性增加，即( $\rho_{\rm C}-\rho^{0}_{\rm C}) \propto \alpha $ ，这里 $\rho^{0}_{\rm C}$ 为 $\alpha=0$ 时的接触密度.但是， $\rho_{\rm C}$ 随α增加缓慢，α从0增至0.28时，接触密度仅增加约18%.这表明亲液杂质对密度振荡程度的影响并不明显，从图 2(a)也可以看出该规律.

从不同亲液杂质分布形式下模拟结果发现，集中和均匀分布对应的液体密度分布几乎重合，如图 3(a)所示，其中 $\alpha =0.28$ .为定量衡量二者的差异，这里给出了层内平均密度差 $\overline {\Delta \rho } $ 随α值的变化关系，如图 3 (b)所示，其中纵坐标取为平均密度差与液体体相密度的比值，误差线为 $\Delta \rho_{i}$ 的标准差.平均密度差 $\overline {\Delta \rho } $ 定义为两种杂质分布条件下对应片层内密度差的统计平均值.

图  3  两种分布形式密度分布比较

Figure  3.  The comparison between fluid density profiles at different distributions of wetting impurity

下载: 全尺寸图片 幻灯片

其中， $\rho _i^1 $ 和 $\rho _i^2 $ 分别为亲液杂质集中和均布条件下第 $i$ 层液体的密度.由于 $\rho _i^1 $ 和 $\rho _i^2$ 的差别主要体现在前两个密度振荡区域，这里统计时 $i$ 的范围仅覆盖从壁面至密度分布的第2个极小值点.从图 3(b)可见， $\overline {\Delta \rho }$ 在不同α取值时均接近为零，说明相同亲液杂质占比时，杂质集中和均匀分布形式下，液体密度分布完全一致，即杂质分布形式不影响流体密度分布规律.

2.2   液体速度分布和滑移

图 4为亲液杂质占比及分布形式不同时液体在通道内的速度分布，其中纵坐标为速度，单位为 $\varepsilon/\sigma $ ，横坐标为通道内的位置，坐标原点在通道中心.由于液体的速度分布沿通道中心对称，为增强对比性，图 4的左右半边同时给出亲液杂质集中和均匀分布时纳米通道内的速度分布曲线.因 $\alpha = 0.02$ 代表仅有一个亲液原子，没有集中和均布分布之别，因而该条件下左右两边速度分布曲线关于通道中心线对称.且此时由于亲液杂质占比很小，对流场几乎无影响，因此速度分布与 $\alpha = 0$ (单纯疏液壁面)基本重合.

图  4  亲液杂质占比及分布对液体速度分布的影响

Figure  4.  The influence of percent of wetting impurity atoms and its distribution on velocity of liquid

下载: 全尺寸图片 幻灯片

随着α增加，两种杂质分布形式下主流区液体的速度分布曲线逐渐降低，表明α的增加会减弱固液界面滑移效果.且相同α时，均匀分布的亲液杂质会导致液体速度分布曲线比集中分布时更低.由此可见，滑移和法向密度分布特性之间并不总是存在直接的对应关系，或者滑移并不单独依赖于液体的法向密度分布特性.

为进一步定量考察亲液杂质对滑移量的影响，图 5又给出了不同亲液杂质占比及分布形式时液体的滑移长度 $(L_{\rm S} )$ .其中，散点为MD模拟结果，虚线通过方程 $L_{\rm S} = B\alpha + L_{\rm S}^0 $ 拟合得到，这里 $L_{\rm S}^0 $ 为 $\alpha = 0$ 时的滑移长度.从图 5中可见，不同亲液杂质分布形式时， $L_{\rm S}$ 均随α以线性规律迅速减小，且在亲液杂质均匀分布时滑移长度减小得更快.如在 $\alpha = 0.28$ 时，集中和均匀分布的亲液杂质使 $L_{\rm S} $ 比单纯疏液壁面分别降低约50\%和56\%.这表明少量亲液杂质就会对疏液表面的滑移量产生显著影响，且亲液杂质均匀分布后滑移量的减小程度更大.

图  5  亲液杂质占比及分布对滑移长度的影响

Figure  5.  The influence of percent of wetting impurity atoms and its distribution on slip length

下载: 全尺寸图片 幻灯片

2.3   亲液杂质对滑移的影响机理

Lichter等^[31]的研究证明, 滑移的本质是第一液体层中原子的比率过程.第一液体层定义为壁面至密度分布的第一个极小值点之间的区域.第一层液体原子由于受到壁面周期性势能的约束而优先占据势能最低点.部分第一层液体原子会由于受到主流区原子的扰动将越过势能的最高点跃迁至相邻的势能最低点.当原子在各个方向上的跃迁概率相同时，第一液体层的统计平均速度为零，即无滑移.当沿流向的跃迁概率大于逆流向的概率时，第一层液体原子统计平均速度不为零，即出现滑移. Wang等^[21]基于该原理给出了滑移速度的定量表达式

其中， $k_{\rm B} $ 为玻尔兹曼常量， $h$ 为普朗克常量， $T$ 为绝对温度， $F_ + $ 和 $F_0 $ 分别为激发态和平衡态配分函数， $S$ 为沿流向第一液体层单个分子的平均有效面积，λ为相邻固体原子平衡位置的距离， $\widetilde{\tau }$ 为施加的剪切应力. $\Delta E$ 为不存在流动时，原子跃迁所需克服的能垒.剪切应力的存在使第一液体层原子沿流向跃迁所需克服的能垒降低了 $\widetilde{\tau }\lambda S / 2$ ，而使原子逆流向跃迁的能垒增加了 $\widetilde{\tau }\lambda S / 2$ ，从而导致原子沿流向跃迁的概率大于逆流向的概率.这一过程由方程(4) 中 $v'_{\rm S} $ 的表达式描述.当 $v'_{\rm S}$ 大于第一液体层原子的热运动速度的涨落( $v_{\rm f}$ )时，滑移出现，否则施加的剪切力只会增加液体的热运动.

通常 $F_ + $ 和 $F_0 $ 相等，且文中温度和驱动力均保持不变，因此， $k_{\rm B} TF_ + / hF_0 $ 和 $\widetilde{\tau }\lambda S / 2$ 为常数.由此可得滑移的大小取决于不存在流动时原子跃迁所需克服的能垒 $\Delta E$ . 图 6为固定法向位置在密度分布的第一个极大值点处，壁面势能沿平行通道方向的变化曲线，其中纵坐标为势能，横坐标为沿通道方向的位置，空心圆为疏液原子，实心圆为亲液原子.在该位置处的势能曲线已可反映壁面势能分布沿平行壁面方向的变化特点，且第一液体层原子也主要集中于第一个极大值点附近.由于壁面势能具有周期性，因此对于集中分布的亲液杂质，需计算含有3个亲液杂质原子的情况，如图 6(a)所示.而亲液杂质在均布形式下，相邻两个亲液杂质原子之间的距离较远，因此仅需计算含有一个亲液杂质原子的势能曲线，如图 6(b)所示.

图  6  亲液杂质和集中分布时壁面第一层原子产生的势能分布

Figure  6.  The potential distribution of first layer of wall as the wetting impurity atoms were distributed uniformly and concentrated

下载: 全尺寸图片 幻灯片

当壁面含有亲液杂质时，由于亲液杂质与液体之间的相互作用强度较大，因此第一液体层原子在跨越亲液原子时需克服更大的能垒，致使原子发生跃迁的概率下降从而滑移减小.而随着亲液原子占比增加，较大的能垒在壁面产生势能中的占比增加，因此，滑移量会逐渐减小.

虽然亲液杂质集中和均匀分布时，由亲液原子产生的能垒基本相同, 但均匀分布时，由于与液体相互作用强度不同，致使亲液原子与相邻两个疏液原子产生的势能极小值点之间的距离略大于相邻两个固体原子之间的距离λ, 这等效于降低了液体与固体之间的非公度性(incommensurability)，从而使得滑移量减小(文献[10, 20]均证实较低的非公度性会导致较小的滑移量); 而亲液杂质集中分布时，亲液原子相邻排列，对固液之间的非公度性影响很小(仅在亲液杂质最边缘处会有影响), 因此，均匀分布的亲液杂质相较于集中分布会进一步降低流动系统的滑移量.

3.   结论

采用非平衡MD模拟方法，研究了含有少量亲液杂质的疏液表面对纳米通道内液体流场和滑移特性的影响.

由于亲液杂质与液体之间相互作用强度较大，亲液杂质的存在会增强液体在固体表面的类固体现象；同时随着亲液杂质占比增加，液体的类固体现象更明显，且液体在固体表面的接触密度与亲液杂质占比呈线性关系; 而相同亲液杂质占比条件下，不同亲液杂质分布形式对液体的密度分布并无影响.

亲液杂质的占比和分布形式对液体的速度分布和滑移特性均有着显著的影响.在亲液杂质集中和均匀分布条件下，滑移长度均随亲液杂质占比增加而线性减小，且亲液杂质均匀分布时滑移长度减小的更快.少量亲液杂质就会对壁面滑移产生显著的影响，从而减弱疏液表面降低流动阻力的能力.

亲液杂质在壁面局部会增加第一液体层原子在发生跃迁行为时的能垒，从而降低壁面滑移.而相较于集中分布，均匀分布的亲液杂质由于降低固液之间的非公度性，致使其对滑移特性的破坏更严重.