水滴撞击结冰过程的分子动力学模拟∗

董琪琪 胡海豹2) 陈少强 何强 鲍路瑶

1)(西北工业大学航海学院,西安 710072)

2)(西北工业大学深圳研究院,深圳 518057)

3)(中国船舶重工集团公司第七〇五研究所,西安 710077)

(2017年10月6日收到;2017年12月11日收到修改稿)

1 引 言

水滴结冰是常见的自然现象,也给人类生活带来无数的挑战.过度积冰会对基础设施造成严重破坏,如建筑物、输电线路、风力涡轮机、太阳能电池板、船舶和气动机翼结构等[1−9].这些结冰现象实际为水滴撞击到冷表面的动态结冰问题.因此,为真实反映壁面水滴结冰的规律,亟需深入研究水滴撞击结冰过程.

在实验研究方面,Quero等[10]在研究过冷水滴撞击水膜表面的冻结过程时发现,与静态结冰相似,水滴撞击冷表面的冻结过程受表面温度和环境温度的影响最大.Mishchenko等[11]通过水滴撞击超疏水表面结冰实验发现,表面温度在−25—−30°C之间的超疏水表面可使水滴在形核之前回弹,能抑制结冰.Jung等[12]发现材料的抑冰性能同时受接触角和表面粗糙度的影响,同等粗糙度情况下,超疏水材料比亲水材料具有更强的抑冰性能.Li等[13]研究了冻结对铝板表面水滴撞击物理过程的影响,结果表明,冻结不影响水滴的扩散过程,但可以减小回缩过程.Yang等[14]对过冷水滴撞击不同材质金属管表面的冻结机理进行了研究,指出除环境温度和冷表面温度外,水滴过冷度、表面特性对结冰也有影响,且表面温度越低撞击速度对结冰的影响越弱.

在数值模拟方面,张大林等[15]模拟发现,随着过冷水滴平均直径的增大,结冰区的面积和水收集系数越大.杨倩等[16]数值模拟给出了飞行高度、速度及水滴半径对撞击特性的影响规律.陈科和曹义华[17]则发现飞机发生结冰时,冷空气会先在飞机表面凝华成霜,之后再诱导结冰.盛强等[18]对机翼流场进行数值模拟,在计算机翼表面的结冰量的同时,分析了结冰对机翼阻力、升力和压力系数的影响.

需要指出的是,现有研究报道主要关注了宏观尺度下水滴撞击冷表面的结冰过程,对微观尺度下冷壁面上水滴撞击结冰过程仍缺乏精细的探究手段.为此,本文采用三维分子动力学模拟方法,从纳米尺度研究水滴撞击冷壁面的结冰过程,并在分析水滴结冰判断方法的基础上,初步探究了温度和表面能对水滴结冰的影响规律.

2 分子动力学模拟方法

2.1 模拟系统

建立了一个包含纳米级水滴和固体壁面的三维分子动力学模型,如图1所示.其中,固体壁面长度L=262.0875 Å,壁厚D=30.7 Å,z方向长度为289.2 Å.

壁面采用fcc晶格结构,晶格常数a=3.615 Å,整个壁面分为9层,共有196667个固体原子组成.其中,位于壁面最下端的3层原子为固定原子,作用是阻止模型中的壁面原子和水滴原子下移;中间4层原子为温度控制原子组,用于调控整个模型,使其保持在一个稳定的温度;最上端的2层原子为自由原子组,该组原子的厚度为7.23 Å.

直径为几个纳米的水滴就可以准确反映水滴特性[19−21],因此,为减少计算量,选取水滴直径为9 nm.水分子共15625个,由15625个氧原子和31250个氢原子组成.模型中的气相成分则是由水分子在相应模拟条件下自由扩散形成,且所有气相的密度、原子质量以及势能参数均接近真实的气相环境.

图1 模拟系统示意图Fig.1.Schematic of simulation system.

2.2 数学模型

这里使用TIP4P/ice模型来模拟氧原子、氢原子和水分子,具体势能函数形式如下:

其中roo为两氧原子之间的距离,ε和σ是LJ势能的特征能量和特征长度,e为质子电荷,ε0为介电常数,a和b为分子的带电位置.水分子的键长和键角则使用SHAKE算法来进行固定.

温度低于熔点时,水滴处于固态,水分子会表现出近程和长程均有序的特性;温度高于熔点时,水滴处于液态,水分子会表现出近程有序、长程无序的特性.利用TIP4P/ice模型,能够自动捕捉到水滴结冰和融化过程.在其他特性不受影响的情况下,该模型预测的六边形冰的最优熔点为−0.95°C/1 bar(1 bar=0.1 MPa)[22].因此,该模型非常适合研究水滴撞击结冰.值得注意的是水滴的表面张力(大量水分子间实际存在的相互作用的统计结果)已经包含在了TIP4P/ice模型中.

另外参考文献[23],采用LJ/126模型模拟水分子和固体原子之间的相互作用,通过调控势能参数来模拟不同润湿性壁面.

2.3 模拟算法

模拟中,统计系统均采用微正则系综(NVE),温度控制策略采用速度定标法.而在速度定标法中,采用场基温控法(profile-unbiased thermostat)[24]来计算温度.系统中所有水滴分子运动均满足牛顿运动方程:

这里采用文莱特算法来迭代求解这些运动方程,

模拟时间步设定为0.002 ps,运行总步数为400000步.

2.4 结冰判断方法

相比于固体分子,液体分子之间空隙较大,分子不会长时间地在同一个固定位置上振动,而是振动一段时间后,转到另一个平衡位置上振动,即液体分子可以在液体内移动.因此,若水分子中氧原子位置保持不变,则可判定其所在位置已从液态转变为固态.图2为采用上述模型模拟出的−40°C冷壁面上10260号氧原子的运动轨迹,可以看出,约75 ps之后氧原子三个方向的坐标基本保持不变(仅在该位置振动),即此时该氧原子所属水分子被冻结.表明通过此方法理论上可以分析出每个水分子的结冰时间.

图2 10260号氧原子的轨迹变化坐标图Fig.2.The coordinate graph of trajectory variation of oxygen atom No.10260.

但由于水滴中含水分子数众多,对不同时刻下各分子位置的统计过于复杂,致使计算量在现有硬件条件下无法承受.为此,这里利用TIP4P/ice模型中水的凝固温度恒为−0.95°C的特点,提出了一种通过统计不同时刻水滴内水分子的温度来判断水滴结冰的方法.考虑到水滴结冰自下而上的规律,这里采用垂直分层统计温度的方法来简化计算量.当水滴顶部温度也达到结冰点(−0.95°C)时,判定水滴完全结冰.

3 模拟结果与分析

3.1 温度对水滴结冰的影响

模拟4种壁面温度下的水滴撞击结冰过程.模拟参数:固-液/液-液势能比ε=11.5,水滴起始温度20°C,撞击速度为300 m/s,壁面温度与空间温度(环境温度)相同,分别为20,−40,−50,−100°C.

水滴撞击结冰过程如图3所示,可以看出,冷壁面上水滴最大铺展直径明显小于常温,且水滴最大铺展直径随着壁面温度的降低而减小.同时,图3中水滴撞击20°C壁面时,壁面上固/液接触线在740 ps趋于稳定;而撞击−40,−50,−100°C壁面时,接触线分别在95,81和60 ps后基本保持不变,且壁面温度越低,水滴接触线振荡时间越短.胡海豹等[25]的实验结果表明,同种材料的壁面上,温度T越低,最大铺展系数βmax越小,而且,低温壁面上固/液接触线在振荡过程中始终保持不变.这与图3中的仿真结果相一致,证明此模型是可行的.

图3 光滑壁面上水滴撞击结冰时序图Fig.3.Sequence diagram of freezing behavior of water droplet on the smooth wall.

水滴撞击−50°C和−40°C壁面时,垂直方向各层水分子的温度变化情况如图4所示.分析−40°C壁面水滴撞击结冰过程(图4(a))可知,120 ps时水滴底部(靠近固体层)的水分子温度已低于结冰点,而水滴顶部的水分子温度仍高于结冰点,这说明在水滴由底部到顶端的结冰过程中存在着较大的温度梯度136 ps时水滴顶部的温度达到结冰点,水滴完全结冰;330 ps时水滴完全降至壁面温度.而水滴撞击−50°C壁面结冰过程中(图4(b)),132 ps时水滴顶部的温度就已达到结冰点,471 ps后水滴才完全降至壁面温度.说明壁面温度越低,水滴完全结冰所需的时间越少,但水滴降至壁面温度所需的时间却越多.这是由于随着壁面温度的降低,固/液间热流密度增大,单位时间内水滴被吸收的热量增多,从而使水滴整体温度快速降低,减少了完全结冰的时间;但水滴需要降低较多的温度和进行较多的热量传递才能降至较低的壁面温度,并且越接近壁面温度时,热量传递的速率就会越慢,所以壁面温度越低,水滴降至壁面温度所需的时间越多.

图4 不同温度壁面水滴垂直方向各层水分子的温度分布(a)T=−40°C;(b)T=−50°CFig.4. Thetemperaturedistributionofwater molecules in the vertical direction of different temperature wall:(a)T= −40°C;(b)T=−50°C.

3.2 壁面表面能属性对水滴结冰的影响

固/液间作用势能可以影响壁面的表面能属性,作用势越大,壁面越亲水,反之越疏水.本文模拟了不同表面能壁面上水滴撞击结冰过程,其中,水滴起始温度为20°C,撞击速度为300 m/s,壁面温度和环境温度为−40°C,固-液/液-液作用势能比(ε)依次为8,6,2(壁面亲水性越来越弱).

图5 不同表面能壁面水滴垂直方向各层水分子的温度分布 (a) ε =8;(b) ε =6;(c) ε =2Fig.5. Thetemperaturedistributionofwater molecules in the vertical direction of different surface energy wall:(a) ε =8;(b) ε =6;(c) ε =2.

图5为不同表面能壁面水滴垂直方向各层水分子的温度分布.可以明显看出,同一时刻下,ε=2时的温度分布明显高于ε=8,6时的温度分布;且ε=2时液滴各分子层间温差较小,内部温度相对均匀.ε=8,6,2时,水滴完全结冰时间分别为133,136,198 ps,水滴完全降至壁面温度所需时间分别为462,560,755 ps;即随着亲水性的减弱,水滴完全结冰和降至壁面温度的时间都逐渐增多.这是由于随着壁面亲水性降低,水滴内部的热量传递的速率会减慢(特别是冷壁面与近壁面的水分子层间的热传递),延长了结冰时间.

4 结 论

采用分子动力学方法模拟研究了水滴撞击冷壁面的结冰过程,结果发现:

1)通过统计水滴垂直方向水分子的温度判定水滴是否结冰比通过统计微观原子的位置坐标来判定更为简洁;

2)壁面温度是水滴结冰现象的重要影响因素之一,壁面温度越低,水滴完全结冰所需的时间越少,但降至壁面温度所需的时间却越大;

3)壁面亲水性是水滴结冰现象的另一重要影响因素,随着壁面亲水性降低,水滴内部的热传递速度减慢(尤其是冷壁面与水滴底端分子层间),内部温度趋于均匀,延长了水滴结冰时间.

另外,水滴晶格与表面晶格的匹配度同样会影响水滴的结冰现象[26,27].因此,后续仍需要进一步研究晶格结构对结冰特性的影响,揭示其对结冰过程的影响机制.

[1]Jung S,Tiwari M K,Doan N V,Poulikakos D 2012Nat.Commun.3 615

[2]Jin Z,Wang Z,Sui D 2016Int.J.Heat.Mass.Trans.97 211

[3]Wang Y,Orol D,Owens J,Simpson K,Lee H J 2013Mater.Sci.Appl.04 347

[4]Dalili N,Edrisy A,Carriveau R 2009Renew Sust.Energ.Rev.13 428

[5]Zou L,Xu H J,Gong S K,Li D W 2010China Safety Sci.J.20 105(in Chinese)[周莉,徐浩军,龚胜科,李大伟2010中国安全科学学报20 105]

[6]Xiao S,He J,Zhang Z 2017Acta Mech.Solida Sin.30 224

[7]Yao Y,Li C,Zhang H,Yang R 2017Appl.Surf.Sci.419 52

[8]Zou M,Beckford S,Wei R,Ellis C,Hattonc G,Millerb M A 2011Appl.Surf.Sci.257 3786

[9]Zhang C,Liu H 2016Phys.Fluids28 260

[10]Quero M,Hammond D W,Purvis R,Smith F T 2006AIAA466

[11]Mishchenko L,Hatton B,Bahadur V,Taylor J A,Krupenkin T,Aizenberg J 2010ACS Nano4 7699

[12]Jung S,Dorrestijn M,Raps D,Das A,Megaridis C M,Poulikakos M 2011Langmuir27 3059

[13]Li H,Roisman I V,Tropea C 2011Proceeding of the Sixth International Conference on Fluid Mechanics1376 451

[14]Yang G,Guo K,Li N 2011Int.J.Refrig.34 2007

[15]Zhang D L,Yang X,Ang H S 2003J.Propul.Power18 87(in Chinese)[张大林,杨曦,昂海松2003航空动力学报18 87]

[16]Yang Q,Chang S N,Yuan X G 2002Acta Aeronaut.Astronaut.Sin.23 173(in Chinese)[杨倩,常士楠,袁修干2002航空学报23 173]

[17]Chen K,Cao Y H 2008Aeronaut.Comput.Tech.38 36(in Chinese)[陈科,曹义华2008航空计算技术38 36]

[18]Sheng Q,Xing Y M,He C 2009Aeronaut.Comput.Tech.39 37(in Chinese)[盛强,邢玉明,何超 2009航空计算技术39 37]

[19]Yuan Q Z,Zhao Y P 2010Phys.Rev.Lett.104 246101

[20]Xiao S,He J Y,Zhang Z X 2016Nanoscale8 14625

[21]Bi Y,Cao B,Li T 2017Nat.Commun.8 15372

[22]Abascal J L,Sanz E,García F R,Vega C 2005J.Chem.Phys.122 234511

[23]Hong S D,Ha M Y,Balachandar S 2009J.Colloid Interf.Sci.339 187

[24]Evans D J,Morriss G P 1986Phys.Rev.Lett.56 2172

[25]Hu H B,He Q,Yu S X,Zhang Z Z,Song D 2016Acta Phys.Sin.65 104703(in Chinese)[胡海豹,何强,余思潇,张招柱,宋东2016物理学报65 104703]

[26]Fitzner M,Sosso G C,Cox S J,Michaelides A 2015J.Am.Chem.Soc.137 13658

[27]Liu K,Wang C,Ma J,Shi G,Yao X,Fang H,Song Y L,Wang J J 2016Proc.Natl.Acad.Sci.USA113 14739