王宝和 ， 李瑞静 ， 王 维

(1.大连理工大学 化工学院 ， 辽宁 大连 116024 ； 2.大连理工大学 化工机械与安全学院 ， 辽宁 大连 116024)

•开发与研究•

纳米氩液膜蒸发行为的分子动力学模拟研究

王宝和1， 李瑞静1， 王 维2

(1.大连理工大学 化工学院 ， 辽宁 大连 116024 ； 2.大连理工大学 化工机械与安全学院 ， 辽宁 大连 116024)

采用非平衡分子动力学模拟技术，探讨模拟温度、能量系数、粗糙度因子、相面积分数等对纳米氩液膜蒸发行为的影响。模拟结果表明：纳米氩液膜在光滑壁面或纳米结构粗糙壁面上的蒸发，存在恒速蒸发和降速蒸发两个阶段；在各个蒸发阶段，蒸发通量相差不大；蒸发过程达到平衡后，纳米栏栅形粗糙壁面吸附的氩原子数最多，纳米方柱矩阵形的次之，光滑壁面的最少。对于光滑壁面上的纳米氩液膜蒸发，随着模拟温度的提高，恒速蒸发阶段的时间变短，蒸发通量迅速增大；在恒速蒸发阶段，能量系数对蒸发通量几乎没有影响；在降速蒸发阶段，随着能量系数的增加，固体壁面吸附的氩流体分子数增大。对于纳米栏栅形粗糙壁面上的纳米氩液膜蒸发，随着粗糙度因子或相面积分数的增大，恒速蒸发阶段的时间增加，蒸发通量减小；当蒸发过程达到平衡后，固体壁面吸附的氩原子数增多。

纳米氩液膜 ； 分子动力学 ； 模拟 ； 蒸发

众所周知，由于固体壁面上纳米薄液膜的厚度处于纳米尺度， 使得其蒸发过程的实验研究工作难以开展，但随着计算机技术的迅猛发展，采用分子动力学(MD)模拟技术，对这样一个传递过程进行研究已经成为可能。Yu等[1]利用分子动力学模拟技术，探讨了脱离压力对氩液膜蒸发—冷凝过程的影响，并获得了纳米薄液膜厚度的计算方法。Hens等[2]利用分子动力学模拟方法，研究了氩纳米液膜在固体铂壁面上的蒸发和沸腾过程，探究了温度和固体表面润湿性的影响。Shavik等[3]通过分子动力学模拟手段，研究了固—液之间的相互作用势能以及固—液界面润湿性(包括疏液性、亲液性和中性)对固体铂壁面上纳米氩液膜的蒸发和爆炸沸腾过程的影响规律。Wang等[4]借助于分子动力学模拟方法，探索了纳米氩液膜在铝纳米结构粗糙表面上的蒸发和爆炸沸腾过程。模拟结果表明：与光滑表面相比，在一定范围内，表面的粗糙纳米结构有利于提高蒸发速率。本文拟通过非平衡分子动力学模拟方法，研究纳米氩液膜在固体壁面上的蒸发特性，探讨模拟温度、能量系数、粗糙度因子、相面积分数等因素对纳米氩液膜蒸发行为的影响。

1 模拟方法

1.1 模拟体系的建立

模拟体系的初始构型如图1所示，模拟盒子尺寸为6.1 nm×6.1 nm×12.8 nm(包括2 nm的删除区)。3 136个氩原子随机分布成铺满壁面的纳米液膜，氩原子的初速度由随机数发生器确定[5]。液膜上方为真空区域。壁面采用Pt原子参数[ε(Pt)=8.35×10-20J；σ(Pt)=2.475×10-10m，讨论能量系数影响时除外]，以FCC的方式排布。壁面原子数量与粗糙度因子或相面积分数大小有关。

1.固体壁面 2.纳米氩液膜 3.删除区

1.2 势能模型的选取

氩流体分子之间的作用势能，采用Lennard-Jones(L-J)势能函数进行描述，分子i和j之间的相互作用势能函数U，如式(1)所示[6]。

(1)

式中：U为势能函数，N为模拟分子个数，rij为分子i和j之间的距离，ε为能量参数，σ为尺寸参数；对于氩流体，ε=1.655×10-21J，σ=3.405×10-10m。

固体壁面原子与氩流体分子间的相互作用势能，依然采用L-J的12-6势能形式，如式(2)所示。相应的尺度参数采用Lorentz-Berthelot混合规则来确定，如式(3)所示[7]，相应的能量参数如式(4)所示[8]。

(2)

σLS=(σL+σS)/2

(3)

εLS=ErεL

(4)

式中：rLS为固体壁面原子与氩流体分子之间的距离，σS为固体壁面原子之间的尺度参数，σL为氩流体分子之间的尺度参数，σLS为固体壁面原子与氩流体分子之间的尺度参数，εL为氩流体分子之间的能量参数，εLS为固体壁面原子与氩流体分子之间的能量参数，Er为能量系数。由式(4)可见，可以通过调整能量系数Er的大小，来考察不同表面能材料对润湿性的影响[8]。

1.3 模拟细节

模拟在x和y方向均采用周期性边界条件，z方向采用固壁边界条件。系统采用NVT系综进行模拟。采用速度Verlet法积分运动方程，时间步长为1 fs。并用Woodcock变标度恒温法进行恒温控制。截断半径为10。模拟时间为3 ns，前1 ns为弛豫时间，确保1 ns后系统达到平衡状态。采用删除分子法，实现非平衡蒸发过程，即液相氩原子受热气化脱离液膜之后，一旦进入删除区(如图1所示)，模拟系统自动将其删除。LAMMPS程序设置每隔0.002 ns，最多删除20个氩原子[7]。

1.4 纳米粗糙壁面的构建

本文构建的纳米栏栅形和纳米方柱矩阵形粗糙壁面的结构参数，分别如图2和图3所示。其中，凸起部分的宽度为a，凸起高度为h，凹陷部分的宽度为b。

纳米栏栅形和纳米方柱矩阵形粗糙壁面的粗糙度因子r1、r2与结构参数的关系，如式(5)和(6)所示；其相面积分数f1、f2与结构参数的关系，如方程(7)和(8)所示[7]。

图2 纳米栏栅形粗糙壁面的结构参数示意图

图3 纳米方柱矩阵形粗糙壁面的结构参数示意图

(5)

(6)

(7)

(8)

1.5 蒸发通量的确定

蒸发通量是指单位时间内离开固体表面吸附区域的原子数除以固体表面的面积。Amsterdam规则认为，对于氩原子i，若以r=1.5 为半径的球体范围内存在至少4个氩原子，则可以认为氩原子i属于液相[9]。在模拟计算过程中，对模拟得到的轨迹文件每隔0.4 ns统计汽、液相原子数，就可以计算出1 ns离开液膜的氩流体分子总数，再根据固体表面积，就可以得到蒸发通量的统计值。

2 结果与讨论

2.1 纳米氩液膜在光滑壁面上的蒸发

2.1.1 能量系数的影响

当氩流体分子数为3 136，模拟温度为100 K，采用MD模拟计算得到的能量系数Er=0.4～10时，纳米氩液膜在光滑壁面上的液相氩原子数与时间的关系，如图4所示；蒸发通量与时间的关系如图5所示。由图4和图5可见，纳米氩液膜的蒸发过程存在恒速蒸发和降速蒸发两个阶段。在过程的初期，蒸发通量近似不变，即为恒速蒸发阶段；然后，随着过程的继续进行，蒸发通量不断下降，即为降速蒸发阶段，直至蒸发通量趋近于零的平衡过程。在恒速蒸发阶段，能量系数对蒸发通量几乎没有影响；在降速蒸发阶段，随着能量系数的增加，固体壁面吸附的氩流体分子数增多。

图4 不同能量系数下的液相氩原子数随时间的变化

图5 不同能量系数下的蒸发通量随时间的变化

2.1.2 温度的影响

当氩流体分子数为3 136，能量系数Er=2时，采用分子动力学模拟方法，计算得到的不同模拟温度下，纳米氩液膜在光滑壁面上的液相氩原子数与时间的关系，见图6；蒸发通量与时间的关系见图7。由图6和图7可见，蒸发过程可以分为恒速蒸发和降速蒸发两个阶段；且随着模拟温度的提高，恒速蒸发阶段的时间变短，蒸发通量迅速增大。

2.2 纳米氩液膜在纳米粗糙壁面上的蒸发

2.2.1 固体壁面形貌的影响

由式(5)和(6)可见，当取a=b=h=2个晶格长度时，纳米栏栅形和纳米方柱矩阵形两种粗糙壁面的粗糙度因子相等(r=2)。

图6 不同温度下的液相氩原子数与时间的关系

图7 不同温度下的蒸发通量与时间的关系

当氩流体分子数为3 136，模拟温度为100 K，能量系数Er=2，纳米栏栅形和纳米方柱矩阵形两种粗糙壁面的粗糙度因子都等于2时，采用分子动力学模拟方法，计算得到的纳米氩液膜在两种形貌壁面上的液相氩原子数及蒸发通量随时间的变化，分别见图8和图9。图中还与光滑壁面的MD模拟结果进行了比较，可见，纳米氩液膜在光滑壁面、纳米栏栅形粗糙壁面及纳米方柱矩阵形粗糙壁面上的蒸发过程，都存在恒速蒸发和降速蒸发两个阶段；在各个蒸发阶段，蒸发通量相差不大；但蒸发过程达到平衡后，纳米栏栅形粗糙壁面吸附的氩原子数最多，纳米方柱矩阵形的次之，光滑壁面的最少。

图8 不同形貌壁面上的液相氩原子数随时间的变化

图9 不同形貌壁面上的蒸发通量随时间的变化

2.2.2 粗糙度因子的影响

对于纳米栏栅形粗糙壁面，当取a=b时，由方程(7)计算得到，相面积分数f=0.5；在保持(a+b)=4个晶格长度的情况下，通过改变壁面突起高度h的大小分别为1、2个晶格长度，壁面粗糙度因子r分别为1.5、2.0。

对于纳米栏栅形粗糙壁面，当氩流体分子数为3 136，模拟温度为100 K，相面积分数f=0.50，能量系数Er=2时，采用分子动力学模拟技术，计算得到的纳米氩液膜在粗糙度因子r分别为1.5、2.0时的液相氩原子数及蒸发通量随时间的变化，分别如图10和11所示。可见，随着粗糙度因子的增大，恒速蒸发阶段的时间增加，蒸发通量减小；当蒸发过程达到平衡后，壁面吸附的氩原子数增多。

图10 不同粗糙度因子壁面上的液相氩原子数随时间的变化

图11 不同粗糙度因子壁面上的蒸发通量随时间的变化

2.2.3 相面积分数的影响

对于纳米栏栅形粗糙壁面，当取h=2.5个晶格长度，(a+b)=5个晶格长度时，由式(5)计算得到，粗糙度因子r=2；在保持(a+b)=4个晶格不变的情况下，变化壁面凸起宽度a分别为1、1.5个晶格长度，使相面积分数f分别为0.25、0.5。

对于纳米栏栅形粗糙壁面，当氩流体分子数为3 136，模拟温度为100 K，粗糙度因子r1=2，能量系数Er=2 时，采用分子动力学模拟技术，计算得到的纳米氩液膜在相面积分数f分别为0.25、0.5时的液相氩原子数及蒸发通量随时间的变化，分别如图12和13所示。可见，当粗糙度因子一定时，随着相面积分数的增大，恒速蒸发阶段的时间增加，蒸发通量减小；当蒸发过程达到平衡后，壁面吸附的氩原子数增多。

图12 不同相面积分数壁面上的液相氩原子数随时间的变化

图13 不同相面积分数壁面上的蒸发通量随时间的变化

3 结论

本文采用非平衡分子动力学模拟方法，探讨了模拟温度、能量系数、粗糙度因子、相面积分数等因素对纳米氩液膜蒸发行为的影响，得到如下结论：①纳米氩液膜在光滑壁面或纳米结构粗糙壁面上的蒸发过程，都存在恒速蒸发和降速蒸发两个阶段；在各个蒸发阶段，蒸发通量相差不大；蒸发过程达到平衡后，纳米栏栅形粗糙壁面吸附的氩原子数最多，纳米方柱矩阵形的次之，光滑壁面的最少。②对于光滑壁面上的纳米氩液膜蒸发过程，随着模拟温度的提高，恒速蒸发阶段的时间变短，蒸发通量迅速增大；在恒速蒸发阶段，能量系数对蒸发通量几乎没有影响；在降速蒸发阶段，随着能量系数的增加，固体壁面吸附的氩流体分子数增大。③对于纳米栏栅形粗糙壁面上的纳米氩液膜蒸发，随着粗糙度因子的增大，恒速蒸发阶段的时间增加，蒸发通量减小；当蒸发过程达到平衡后，固体壁面吸附的氩原子数增多。④对于纳米栏栅形粗糙壁面上的纳米氩液膜蒸发，随着相面积分数的增大，恒速蒸发阶段的时间增加，蒸发通量减小；当蒸发过程达到平衡后，固体壁面吸附的氩原子数增多。

[1] Yu J,Wang H.Thin-film phase change driven by disjoining pressure difference[J].Heat and Mass Transfer,2012,48(7):1135-1140.

[2] Hens A,Agarwal R,Biswas G.Nanoscale study of boiling and evaporation in a liquid Ar film on a Pt heater using molecular dynamics simulation[J].International Journal of Heat and Mass Transfer,2014,71:303-312.

[3] Shavik S M,Hasan M N,Morshed A K M M,et al.Molecular dynamics study of effect of different wetting conditions on evaporation and rapid boiling of ultra-thin argon layer over platinum surface[J].Procedia Engineering,2015,105:446-451.

[4] Wang W D,Zhang H Y,Tian C H,et al.Numerical experiments on evaporation and explosive boiling of ultra-thin liquid argon film on aluminum nanostructure substrate[J].Nanoscale Research Letters,2015,10(1):1-14.

[5] Heermann D W.Computer-simulation methods[M].Springer Berlin Heidelberg,1990.

[6] Long L N,Micci M M,Wong B C.Molecular dynamics simulations of evaporation[J].Computer Physics Communications,1996,96:167-172.

[7] 李瑞静.纳米液膜传递特性的分子动力学模拟[D].大连:大连理工大学,2016.

[8] 蒋大林.固液界面润湿的分子动力学研究及实验[D].苏州:江苏大学,2007.

[9] Sumardiono S,Fischer J.Molecular simulations of droplet evaporation by heat transfer[J].Microfluidics and Nanofluidics,2007,3(2):127-140.

Study on Molecular Dynamics Simulation of Evaporation Behavior of Argon Nanofilm

WANG Bao-he1, LI Rui-jing1, WANG Wei2

(1.School of Chemical Engineering , Dalian University of Technology , Dalian 116024 , China ; 2.School of Chemical Machinery and Safety , Dalian University of Technology , Dalian 116024 , China)

Non-equilibrium molecular dynamics simulation technology is used to research influence of factors such as simulated temperature,energy factor,roughness factor,phase area fraction on evaporation behavior of argon nanofilm on solid surfaces.Simulation results show that there are two stages of a constant rate evaporating period and a falling rate evaporating period for evaporation process of argon nanofilm on smooth or nanotextured rough surfaces;evaporation flux has little change during the constant rate evaporating period or the falling rate evaporating period,but there is difference in number of argon molecules absorbed on the solid surfaces and the ranking (large to small) is as follows:fencing form surface,square column matrix form surface,smooth surface after evaporation process reaches equilibrium state.For evaporation of argon nanofilm on the smooth surface,the constant rate evaporating period becomes short and evaporation flux increases quickly with the increasing of simulated temperature,and with the increase of energy factor,evaporation flux changes hardly at the constant rate evaporating stage;however,number of argon molecules absorbed on the solid surface increases at the falling rate evaporating stage.For evaporation argon nanofilm on the fencing form surface,with increase of roughness factor or phase area fraction,the constant rate evaporating period becomes long and evaporation flux decreases;but number of argon molecules absorbed on the solid surface increases after evaporation process reaches equilibrium state.

argon nanofilm ; molecular dynamics ; simulation ; evaporation

2016-10-28

国家自然科学基金(21676042)

王宝和(1959-)，男，副教授，从事不同形貌微纳结构的制备、干燥及分子动力学模拟研究，电话：0411-84986167，E-mail：wbaohe@163.com。

TB383；O561

A

1003-3467(2017)01-0016-05