纳米凹槽表面结冰的分子动力学模拟研究

王军,张悦,夏国栋

(北京工业大学 传热强化与过程节能教育部重点实验室;传热与能源利用北京市重点实验室,北京 100124)

固体表面上过冷水结冰的现象广泛存在于航空[1]、气象[2]、能源[3]等领域.结冰现象实际上是液态水分子在冷表面上发生相变的过程,在相变过程中大量水分子冻结并逐渐累积形成一定厚度的冰.在微观上,结冰过程的开始通常以临界冰核的形成作为标志,随后发生冰生长过程.在自然界中,冰晶的结构多种多样[4].这是因为冰成核过程和冰生长过程总是发生在很宽的温度和压力条件下,导致冰的形成过程具有复杂化和多样化的特点,甚至当温度低至-40 ℃,纯水依然可能维持过冷液态[5].然而在实际生活中,室外温度并没有达到如此低温,水就已经发生相变.这是因为日常生产生活中的水大多含有杂质,如:固体烟灰颗粒[6]、金属表面[7-8]和表面空隙[9]等.发生在固体或者杂质表面的结冰过程,称之为异质结冰.水分子与冷表面之间的相互作用和表面结构会影响界面水的性质,进而影响水结冰过程.扫描隧道显微镜实验揭示了金属表面上覆盖的第一层水具有许多复杂的结构[10].分子动力学模拟研究发现,石墨表面结冰的分层程度取决于界面水的结构[11].尽管如此,学者们在微观层面上对固体表面影响结冰的潜在机制仍知之甚少.因此,继续探索异质结冰的机理具有十分重要的意义.

对于大多数社会活动和工业生产领域而言,结冰是不利的[12],严重时甚至会对公共设施和能源系统带来不必要的危害;而对于低温保存等特殊过程却需要制冰,结冰是有利的[13].因此,合理的控制结冰过程(控制成核位点和冰生长方向等)也是十分必要的.由于冰成核是结冰过程开始的关键步骤,因此许多学者通过改变固体表面与水分子之间的相互作用和表面微观形态来调整异质冰成核过程.研究发现,表面润湿性是影响结冰的重要因素之一.当水在硅表面结冰时,润湿性较好的表面能够显著地促进冰成核,而水在润湿性较差的壁面上成核率则较低[14].晶体结构与冰相匹配的固体表面也可以促进冰成核,但晶格匹配并不是促进冰成核的必要条件[15-17].相关文献研究表明,改变表面粗糙度可以为冰成核提供更多的成核位点并且可以促进或抑制冰成核[18].LUPI等[19]研究发现,水分子更易在曲率较大的碳表面的凹面而不是凸面结冰,其主要原因是凹面更有利于界面处的水分子形成有序结构,从而提高成核率.经典成核理论指出,与凸表面相比,凹表面可以减小临界核的尺寸,因而有利于冰核的形成[20].LI等[21]在石墨烯表面引入纳米凹槽,当凹槽尺寸与冰分子的结构尺寸相匹配时,水分子在纳米凹槽结构形成的受限空间内整齐排列,因而更容易结冰成核,相对于平面石墨烯结构,带有纳米凹槽结构的石墨烯表面的成核率提高了2～3倍.

由上述可知,通过在表面上引入纳米凹槽或者改变表面的粗糙度,可以提高其冰成核率.一般而言,只有当纳米凹槽结构与冰分子的结构尺寸相匹配时,才能获得比较高的成核率.对于表面润湿性较差的表面,其结冰成核率较低,甚至低于均相成核率.通过在疏水表面引入纳米凹槽,则难以提高其成核率.这是因为水分子难以浸润疏水表面的纳米凹槽,形成Cassie态,难以提高其表面异质成核率[21].

基于上述原因,本文设计出一种新型的纳米凹槽表面,在疏水表面上引入亲水纳米凹槽.分子动力学模拟结果表明,新型纳米凹槽表面能够显著提高疏水表面的成核率(接近于亲水表面的成核率),其主要原因是水分子可以被亲水纳米凹槽吸附,在纳米受限空间内形成稳定的冰核结构,从而提高了冰成核率.研究结果为疏水表面结冰过程的微观调控提供了新的思路.

1 模型系统

1.1 分子动力学模型

本文建立的分子动力学模拟系统由底部的固体表面和液体水构成,见图1.采用的分子动力学模拟软件为Large-scale Atomic/Molecular Massively Parallel Simulator(LAMMPS)[22].图1中,固体由具有面心立方结构的金原子(20层)组成,晶格常数为0.408 nm.系统中共有52 88个水分子,采用单原子水分子模型[23-24].z方向上的总长度约为190 nm,考虑到水结冰之后的体积膨胀,系统顶部留有厚度为95 nm的真空层[25],本文只显示了部分真空层(见图1).x和y方向上的尺寸分别为81.6 nm和40.8 nm.系统x和y方向上设置为周期性边界条件,z方向设置为固定边界条件.

如图1(a)所示,本文首先构建了固体平面上液态水冻结模型,固体原子数为8 000个.在图1(b)中,去掉固体表面的部分固体原子,在固体平面上布置了纳米凹槽.本文设计了10个具有不同凹槽尺寸的表面,并进行编号,见表1.各个凹槽的槽深h均为2.04 nm,通过改变凹槽宽度w以获得不同的深宽比α(α=h/w).各表面凹槽宽w及深宽比 详见表1.图1(b)中,凹槽宽为2.448 nm,对应的表面编号为10.

表1 纳米凹槽尺寸Tab. 1 The groove size of each surface

1.2 分子间相互作用势

水分子间相互作用采用Stillinger-Weber(SW)势函数

(1)

(2)

(3)

式中,E代表总能量,而Ф2和Ф3分别代表总能量中的二体项和三体项.总能量包括在截断半径内所有相邻原子l、m和n的能量和.其中,参数Alm=7.05,Blm=0.60,Plm=4.00,qlm=0.00,αlm=1.80,λlmn=23.15,γlm=1.20,cosθlm=-0.33.ε和σ指水分子之间的碰撞直径和势阱深度,详见表2.

表2 LJ势函数参数Tab. 2 Parameter setting value in LJ potential function

固-固和固-液原子之间的相互作用采用Lennard-Jones(LJ)势函数

(4)

其中,参数ε和σ分别代表分子/原子间碰撞直径和势阱深度,rij表示为原子i和j之间的距离.其中Au-water原子之间的势能参数遵循Lorentz-Berthelot法则,即

(5)

(6)

为了模拟亲水及疏水表面,采用了3种不同的固-液相互作用参数ε,分别为Au-water、AuQ-water和AuS-water,见表2.

1.3 接触角

为了直观地描述表面的润湿性强弱,需模拟计算液滴在金表面上的形状和接触角.如图2(a)所示,模拟系统的初始模型由3 000个水分子和金固体平壁面组成.首先,采用Nose-Hoover热浴使系统的温度达到280 K.系统中所有原子维持在此温度下运行至少5 ns以确保液滴在表面上充分铺展开,收集最后1 ns的模拟数据用于接触角的测量和计算.本文至少取5次独立模拟数据的平均值作为最终的接触角.

在不改变固体原子排列结构的条件下,对比分析具有不同亲疏水特性的固体表面异质结冰过程.考虑3种不同的固液界面相互作用势(具有不同的Au-water界面势参数ε),因此,表面与水之间的接触角有所不同.如图2(c)所示,对于Au-water 界面而言,其接触角约为93°,此时,认为表面既不亲水也不疏水.此外,AuQ-water和AuS-water分别表示亲水和疏水的Au原子与水发生相互作用.水分子在疏水表面(AuS)上的接触角较大,约为140°,见图2(b);如图2(d)所示,亲水表面(AuQ)上形成两层类似固体的液体水膜,本文在计算液滴在亲水表面(AuQ)上的接触角时忽略了表面液膜的影响[26-27],进而得到接触角约为45°.

1.4 成核率

统计了各表面系统中液态水分子数与水分子总数的比

(7)

系统中冰的成核率

(8)

式中,V表示水的体积.

1.5 模拟计算过程

如图1所示,将金表面分成紫色、黄色和红色.其中,为了防止表面在模拟过程中发生变形,将表面红色层的金原子固定.紫色和黄色的原子为能够与水分子相接触并发生相互作用的“真实”的金原子.时间积分算法选择Velocity-Verlet算法,时间步长设置为10 fs,截断半径设为1 nm.在模拟结冰的过程中,系统首先在280 K弛豫2 ns,然后在200 K运行至完全冻结(总模拟时长不超过202 ns).采用Nose-Hoover热浴对固体原子(除表面红色层原子外)和液体分子进行控温,并基于遍历原理对系统的温度和能量等物理量做长时间的统计平均计算.为了确定模拟结果的误差,各系统在不同的初始条件下(即各模拟过程中原子的初始速度有所不同,但符合同一个Maxwell速度分布函数)进行至少10次的独立模拟,将独立模拟结果的平均值作为本文的最终结果.

1.6 纳米凹槽表面

纳米凹槽表面主要由两部分组成,即图1(b)中黄色的凹槽内面和紫色的槽外固体原子.本文设计了一种综合考虑表面润湿性和表面粗糙度的新型纳米凹槽表面,在疏水表面上布置不同宽度的亲水纳米凹槽,即新型表面(Case Ⅰ).此凹槽内面为亲水Au原子(AuQ),槽外的上表面为疏水Au原子(AuS),两种表面与水分子之间的相互作用势参数见表2.此凹槽内部表面润湿较好,能够使水分子完全浸润凹槽.新型纳米凹槽表面旨在诱导界面水分子在凹槽内部形成冰核.作为对比,同时模拟了凹槽内和凹槽外表面均为疏水的凹槽表面(Case Ⅱ)、亲水平面结构(Case Ⅲ)以及疏水平面结构(Case Ⅳ)的结冰成核过程.

2 结果与讨论

2.1 表面成核率

图3显示了不同表面上水结冰过程中的成核率.为了进行比较,还计算了平坦Au表面上(见图3中紫色点线)的冰成核率J0为1.71×1033m-3·s-1,该值大约比均质成核率高3个数量级[5].由图3可知,与平坦表面对比,亲水平面(Case Ⅲ,蓝色虚线)能够促进结冰,而疏水平面(Case Ⅳ,绿色点划线)则抑制冰成核.亲水表面容易能够诱导界面水分子在表面形成冰核,而疏水表面不会诱导界面水形成冰核,疏水表面系统中的冰核大多出现在水膜表面或水体当中.与疏水程度相同的平坦表面(Case Ⅳ)相比,各疏水凹槽表面(Case Ⅱ)的表面成核率更低.对于Case Ⅰ而言,其表面成核率随凹槽尺寸的变化而改变,对于某些尺寸的纳米凹槽,其表面成核率较高,甚至接近于亲水表面的成核率.

2.2 新型纳米凹槽结构表面结冰过程

亲水凹槽可以吸附水分子,而疏水凹槽不能吸附水分子.因此新型纳米凹槽结构表面可以通过调节表面润湿性和表面粗糙度的协同作用来控制界面水的性质,并通过选择性吸附水分子来调节冰成核的位置和冰生长的方向,提高水分子在疏水表面结冰的成核率.

基于选择性吸附表面上(Case Ⅰ)结冰的分子动力学模拟研究,计算了水在不同深宽比的凹槽表面上冻结的成核率.如图3所示,随凹槽宽度增大,冰成核率的变化较为复杂,但遵循一定的变化规律:成核率先增大后减小,随后再增大.其中,当表面编号低于3时,成核率随凹槽宽度的增加而增加;在表面编号为3时,达到最大值1.84×1033m-3·s-1.随着凹槽宽度的进一步增大,成核率明显降低.当表面编号高于6时,成核率又迅速增大至1.83×1033m-3·s-1.其中,水在凹槽尺寸w3=1.020 nm、w7=1.836 nm、w8=2.040 nm、w9=2.244 nm和w10=2.448 nm的表面结冰的成核率显著增大,甚至冰成核率接近亲水平面(Case Ⅲ)的冰成核率,也就是说,特定凹槽宽度的新型表面能达到促进疏水表面冰成核的效果.

模拟结冰的初始时刻(即系统的温度达到200 K时)到系统中出现临界成核的这段时间定义为成核时间.基于模拟的结冰过程,计算并统计了水在新型凹槽表面(Case Ⅰ)上结冰的成核时间,见图4.结果表明,当表面编号低于3时,成核时间随凹槽宽度的增加而减小,在表面编号为3 时达到最小值.随着凹槽宽度的进一步增大,等待成核时间逐渐增大.当表面编号大于6时,成核率又迅速减小.因此,系统等待成核的时间较短时,对应系统发生成核的速率较大,反映出系统中的水容易结冰;反之,系统等待成核的时间较长时,对应系统发生成核的速率较小,反映出系统中的水不易结冰.

为深入理解成核率与成核时间变化的机制,进一步分析了水/冰分子在凹槽内的分布情况.图5显示了Case Ⅰ中部分凹槽内部水/冰分子的质量密度分布和结冰情况,图5中ρ和m分别为凹槽内部水/冰分子的密度和质量.由密度分布和系统模拟快照可知,由于凹槽亲水,凹槽上紧密吸附了一层水分子,凹槽空隙内也会积聚水分子.其中,黄色表面接触层的水分子密度最大,而未与表面直接接触的水分子密度较低.

如图5(a)所示,当凹槽宽度为1.02 nm(表面编号为3)时,凹槽内中间两层原子以双层二维冰的形式存在,因此成核率较高[29];当水分子突破凹槽约束时,在凹槽正上方形成多元冰环,多元冰环周围的水分子通过氢键与其建立连接,迅速引发冰生长过程.众所周知,冰结构大多为规则且立体的六方冰.图5的结果表明,异质结冰过程会受到纳米凹槽结构的影响,六方冰结构与宽度小于1.836 nm(即表面编号大于7)的凹槽结构不匹配,因此表面上结冰过程的成核率较低.根据系统模拟快照分析冰分子的位置:当纳米凹槽宽度大于1.836 nm时,凹槽内部的冰分子排列规则,因而成核率逐渐增大.这是由于凹槽宽度增大,凹槽对水分子的纳米约束效应减小,有利于槽内冰核形成并突破凹槽结构的限制,从而加快了结冰进程.此时,凹槽内的六方冰作为初始冰核向z方向生长.

图6展示了冰在凹槽(表面编号为10)内形成并向上生长的过程.模拟发现,当凹槽尺寸大于1.836 nm并且可以容纳整数个冰结构时,临界冰核往往出现在凹槽底部,而不是上表面或者液体水中.因此,随后的冰生长过程也具有普遍的方向性,即沿着z轴的正方向定向生长.冰晶从凹槽内呈树枝状迅速地生长,一旦冰晶突破纳米凹槽的限制,便向周围漫延.

3 结 论

亲水表面容易能够诱导界面水分子在表面形成冰核,而疏水表面不会诱导界面水形成冰核,疏水表面系统中的冰核大多出现在水膜表面或水体当中.通过在疏水表面设置亲水凹坑设计出新型纳米凹槽结构表面,此表面可以通过调节表面润湿性和表面粗糙度的协同作用来控制界面水的性质.当凹槽宽度为1.02 nm时,凹槽内中间两层原子以双层二维冰的形式存在,其成核率接近于亲水表面的成核率; 当凹槽宽度>1.836 nm时,纳米凹槽结构与六方冰结构相匹配,凹槽内部的冰分子排列规则,有利于槽内冰核形成并突破凹槽结构的限制,因而成核率较高.当凹槽宽度<1.836 nm时,纳米凹槽结构与六边形冰结构不匹配,凹槽内不易形成冰核,结冰过程受到抑制.研究结果为控制异质结冰的成核位点和冰生长方向,特别是为疏水表面结冰过程的微观调控提供了新的思路.