黄 正 兴，　孙 童 童，　李 倩 倩

( 大连理工大学 电子科学与技术学院, 辽宁 大连　116024 )



纳米尺度水薄膜液-气相变分子动力学模拟

黄 正 兴*，孙 童 童，李 倩 倩

( 大连理工大学 电子科学与技术学院, 辽宁 大连116024 )

摘要：采用分子动力学模拟，研究了纳米尺度水薄膜厚度和系统温度对液-气相变蒸发率的影响．模拟了水薄膜厚度为2 nm时，温度(375～425 K)对相变蒸发率的影响；同时模拟了温度400 K条件下，水薄膜厚度(2、3、4 nm)对相变蒸发率的影响．结果表明：相同水薄膜厚度下蒸发率随着温度的升高而增大，相同温度下蒸发率随水薄膜厚度的增大而增大，蒸发率随时间呈指数形式递减．相关结果可为微热管及其他依靠内部液体相变传热的散热器件设计提供参考．

关键词：液膜厚度；液-气相变；分子动力学模拟

0引言

随着电子技术的发展，微电子设备的工作频率越来越高而其外形尺寸却越来越小，电子产品过热问题愈显突出[1]．这将在一定程度上影响电子产品的可靠性，导致产品失效[2]．确保微电子器件产生的热量能够快速散出成了微电子工业要研究的重要问题．

微槽道平板热管是一种高效的传热元件，它是利用工作液体在微小空间发生相变进行热量传递的[3]．纳米尺度空间液膜的相变蒸发率与微热管的传热效率息息相关，蒸发率的研究对微热管的小型化及性能改善有着举足轻重的作用．因此，要对微尺度空间液膜的相变蒸发率进行研究，以改善微热管的传热性能．

由于实验器材的限制，在纳米尺度下很难进行相变方面的实验，而且实验法很难预测达到相变平衡的时间，不能对相变后的微观现象和微观机理进行观察分析，因此有必要选取更有效的研究方法[4-5]．分子动力学模拟方法是探究纳米尺度相变传热特性的有效方法，对纳米尺度下相变机理的研究更具优势[6-8]．该方法可以从分子水平进行研究分析，模拟液膜厚度可低至纳米级，得到其他方法难以获得的微观信息，并且可以通过可视化软件观察到整个相变过程中分子的运动状态．因此,本文采用分子动力学模拟方法进行研究．由于LAMMPS提供并支持多种势函数，支持多种形态、多种系综下的模拟，并且具有良好的并行扩展性，本文的计算主体采用开源代码LAMMPS．

很早以前人们就开始对流体的基本性质及其相变现象进行研究，1997年Bertolini等采用平衡分子动力学的方法计算得到水的热导率[9]．1998年Matsumoto对流体的相变现象进行研究[10]．1999年Yang等对微尺度液滴的液-气界面进行研究[11]．2011年Cheng等对Lennard-Jones液体的蒸发进行研究[12]．2012年Yano研究了气-液界面的蒸发凝结现象[13]．而对于纳米空间中的薄液膜相变传热问题，尺度效应是研究的一个重点[14]．但目前对薄液膜厚度与相变蒸发率关系的研究还很少见，因此本文将主要进行这方面的模拟工作．

1水模型的确定

在进行液膜相变的模拟之前，首先对文献中使用的多种水模型进行对比分析．分别计算了温度T=298 K时3种常用水模型(TIP4P、SPC、SPC/E)的密度及热导率，并与实验值进行比较，3种模型使用的势能参数见表1．对比结果(见表2)表明：对于密度的计算，SPC/E模型和TIP4P模型与实验值较为接近，但是同时考虑热导率的比较，选择出最符合实际情况的模型TIP4P模型．这与文献[15-16]给出的结论是一致的．并且TIP4P模型在3个原子中间，H—O—H化学键的对角线上多了1个无质量的第4原子M，该原子为Lennard-Jones态．由于此模型采用4个粒子相互作用，模拟更加准确而被广泛应用[17]．

表1　TIP4P、SPC、SPC/E模型势能参数

表2　密度及热导率模拟值与文献值比较

对选取的模型进行水的基本性质以及液-气相变的模拟，图1为水薄膜的热导率随温度的变化曲线，从图中可以看出：随着温度的升高，热导率近似线性增加．这是因为，温度越高水分子具有的能量越高，热运动越剧烈，因此热导率越大．对于液-气相变的模拟统计了不同温度下密度的分布曲线如图2所示．以上模拟主要得到了以下结论：(1)温度越高水的热导率越大．(2)液-气相变平衡后，从气相到液相的密度随时间是连续变化的，气相区的密度随温度的升高而升高,液相区则相反．(3)温度升高液-气界面区的厚度线性增加．

图1　水薄膜热导率随温度的变化关系曲线

图2　不同温度下密度的分布曲线

(4)液-气相变过程不是一个连续变化的过程，而是一个突变的过程，与文献[18-19]得出的结论是一致的，因而验证了本文模拟方法的准确性．

2模拟方法

本文的模拟分为两部分：等厚度水薄膜不同温度的液-气相变和等温度不同厚度水薄膜的液-气相变．模拟系统包括衬底硅和液体水，液膜厚度为d的水薄膜置于硅上．模拟体系在3个方向都采用周期性边界条件，时间步长为1 fs，分子初始速度按Maxwell分布随机取值．模拟系综为微正则系综NVE，即在分子动力学模拟过程中保持系统分子数N、体积V和能量E不变．采用PPPM方法计算长程作用力，采用Shake算法固定水分子的键长与键角，确保模拟过程中水分子的构型不变．

初始时刻，水分子以饱和液相的平均距离均匀密集分布在体系内部，并且在模拟过程中调整液膜重心回到初始位置．体系Y方向为液膜厚度方向，此方向的一端为硅衬底，另一端留有很长的真空区以便相变后有足够的气相空间．X、Z方向尺寸始终保持不变，均为2 nm．

硅与硅之间采用Tersoff势函数，硅与水分子、水分子内部氧原子与氢原子之间势函数均采用Lennard-Jones双体势[20]，其势能函数形式为

(1)

式中：E是两水分子间的对势，ε和σ为Lennard-Jones参数，ε为能量单位，σ为距离单位，r为内部分子之间的距离．硅与水之间采用的参数为σSi-O=3.44×10-10m，εSi-O=1.76×10-21J．

模拟的过程分为3个阶段：温度控制阶段、弛豫阶段和加热阶段．温度控制阶段通过速度标定法，每一步都重新调整衬底硅和水分子的速度，使系统内所有原子的温度保持在275 K．弛豫阶段将衬底硅的温度保持在275 K，水分子可以自由移动．加热阶段控制衬底硅的温度为400 K，在这个阶段开始发生相变，一段时间后相变达到平衡．由于液体相变达到平衡需要很长时间，模拟时间达到纳秒级，以确保相变平衡．

3模拟结果

3.1水薄膜蒸发率

模拟了水薄膜厚度为2 nm时温度对相变蒸发率的影响，以及温度为400 K时水薄膜厚度对相变蒸发率的影响．待相变平衡后，统计气液共存系统中气态分子数．在本文中，根据固定半径(1.5σO-O) 内具有的氧原子个数判定气态分子，若在设定的半径内某个分子中的氧原子周围有两个以上的氧原子，则认为此分子为液态分子，否则为气态分子．

模拟的第3个阶段(加热阶段)系统温度由275 K开始升温，进行5 ns的分子动力学模拟．在温度升高的过程中，水分子开始发生相变．图3为模拟水薄膜厚度d=2 nm相变前后的水分子瞬态位型X-Y投影图，可以看出模拟进行到1 ns时，系统开始发生相变，并且在相变平衡后(3 ns)，仍然有一定厚度的水薄膜没有蒸发为气态．

图3　水分子瞬态位型X-Y投影图

相变开始时气相区的分子数迅速增加，达到平衡后，蒸发至气相区的水分子数基本保持不变．统计气液共存系统中气态分子数,其随时间的变化呈指数形式增长[21]：

y=y0+C1eτ t

(2)

其中y0表示相变平衡后，蒸发至气相区的分子总数，它受系统温度和水薄膜厚度d的影响．图4为模拟温度为400 K，水分子数为380，水薄膜厚度d=3 nm时，蒸发分子数Nvapor随时间的变化曲线．

图4　气态水分子数随时间的变化曲线

蒸发率的计算方法[21]是：首先对不同情况得到的气态分子数随时间的变化曲线1 ns以后的数据进行最小二乘拟合，利用式(2)，得到参数y0、C1和τ的值，然后将得到的函数对时间进行求导，如式(3)所示．进而求出蒸发率α的表达式(4)，其中M为水的摩尔质量，A为液体薄膜的表面积．

(3)

(4)

3.2温度对蒸发率的影响

对温度为375～425 K，厚度为2 nm的水薄膜进行相变模拟，模拟方法如上所述．分别计算不同温度条件下的蒸发率，图5给出了蒸发率随时间的变化关系曲线．从图中可以看出：蒸发率随温度的升高呈增加趋势，温度为425 K时的蒸发率最大，并且其随温度的升高增加幅度逐渐变小．温度低于400 K时，温度升高蒸发率迅速上升；温度高于400 K时，蒸发率随温度的增加上升幅度变小，达到饱和．

图5　不同温度条件下蒸发率随时间的变化曲线

以上模拟情况均有一部分液体没有变成气态, 未蒸发的水薄膜厚度d随温度升高而变小．5种温度情况下未蒸发水薄膜厚度分别为d375 K=1.59 nm，d387 K=1.49 nm，d400 K=1.29 nm，d412 K=1.19 nm，d425 K=1.00 nm．原因是在纳米尺度下表面力不可忽略，衬底硅与水分子之间存在范德华力，由于范德华力的作用，衬底附近一定厚度范围内的水分子很难脱离硅的引力束缚而蒸发．但是温度升高，水分子运动加剧，具有更高的动能，更容易脱离衬底硅及周围水分子的束缚变为气态，因此随着温度的升高，未蒸发水薄膜厚度减小，蒸发率增加．

3.3水薄膜厚度对蒸发率的影响

用同样的方法，对水薄膜厚度为2、3、4 nm，水分子数分别为250、380、504的薄液膜分别进行相变模拟．图6为不同水薄膜厚度下蒸发率随时间的变化曲线．从图中可以看出，3种厚度水薄膜的蒸发率随时间呈指数衰减．水薄膜厚度越大蒸发率越大，衰减速度越慢，达到平衡所需要的时间越长．由于范德华力的作用以上情况也会产生一定厚度未蒸发水薄膜，这说明范德华力对一定半径内的分子有较大的束缚力．

图6　不同水薄膜厚度下蒸发率随时间的变化曲线

在宏观情况下，由于较厚的液膜升温速度慢，开始发生相变时，表面温度较低，因此蒸发率小，但是纳米尺度液膜呈现出了不同特性．统计3种水薄膜厚度液-气界面处的温度随时间的变化曲线如图7所示，发现液-气界面处的温度波动较大，但升温速率并没有随水薄膜厚度的增加而变慢．因此，在相变开始时水薄膜表面的温度基本相同，而对于厚度较大的水薄膜，为蒸发提供的水分子数更多，使得水薄膜越厚蒸发率越大．这说明不只有表面的水分子发生相变，其内部一部分水分子也参与了相变．因此，对于纳米尺度的相变过程，增加厚度可以提高蒸发率．

图7不同水薄膜厚度下液-气界面处温度随时间的变化曲线

Fig.7Changecurveoftemperatureofliquid-vaporinterfacewithtimeforvaryingwaterfilmthickness

4结语

本文采用分子动力学模拟的方法研究了系统温度和水薄膜厚度对相变的影响，从分子的角度分析了其动力学特性，研究结果表明：

蒸发率随时间呈指数形式递减，其影响因素有系统温度及水薄膜厚度．模拟系统温度越高，水薄膜的厚度越大，蒸发率越大．并且当温度很高(400K)时，升高温度对蒸发率的影响减弱，因此想进一步增加蒸发率就要适当增加水薄膜厚度．

无论水薄膜厚度为2、3或4nm，都有一部分水薄膜没有发生相变，这是由于衬底硅与水分子之间存在范德华力，抑制靠近衬底硅的部分水分子蒸发．未发生相变的水薄膜厚度随着温度的升高而变小．

与宏观情况相比，微尺度下液膜具有更大的蒸发率，因此对于微热管及其他依靠内部液体相变传热的传热器件来说，增加微尺度沟道数量、增加液膜厚度和适当升高温度能够有效改善传热性能．

参考文献：

[1] 孙中原. 微小型槽道平板热管传热特性研究[D]. 南京:南京航空航天大学, 2006.

SUNZhong-yuan.Researchontheheattransfercharacteristicsofflatminiaturegroovedheatpipe[D].Nanjing:NanjingUniversityofAeronauticsandAstronautics, 2006. (inChinese)

[2]吕永超,杨双根. 电子设备热分析、热设计及热测试技术综述及最新进展[J]. 电子机械工程, 2007, 23(1):5-10.

LYUYong-chao,YANGShuang-gen.Areviewofthermalanalysis,thermaldesignandthermaltesttechnologyandtheirrecentdevelopment[J].Electro-MechanicalEngineering, 2007, 23(1):5-10. (inChinese)

[3]张丽春. 微小型热管的传热特性及蒸发薄液膜的研究[D]. 合肥:中国科学技术大学, 2003.

ZHANGLi-chun.Flatminiatureaxiallygroovedheatpipesandthemicroscaleevaporativethinliquidfilms[D].Hefei:UniversityofScienceandTechnologyofChina, 2003. (inChinese)

[4]王遵敬. 蒸发与凝结现象的分子动力学研究及实验[D]. 北京:清华大学, 2002.

WANGZun-jing.Moleculardynamicsresearchandexperimentofevaporationandcondensation[D].Beijing:TsinghuaUniversity, 2002. (inChinese)

[5]YiP,PoulikakosD,WaltherJ, et al.Moleculardynamicssimulationofvaporizationofanultra-thinliquidargonlayeronasurface[J].InternationalJournalofHeatandMassTransfer, 2002, 45(10):2087-2100.

[6]WangCS,ChenJS,ShiomiJ, et al.Astudyonthethermalresistanceoversolid-liquid-vaporinterfacesinafinite-spacebyamoleculardynamicsmethod[J].InternationalJournalofThermalSciences, 2007, 46(12):1203-1210.

[7]NagayamaG,KawagoeM,TokunagaA, et al.Ontheevaporationrateofultra-thinliquidfilmatthenanostructuredsurface:Amoleculardynamicsstudy[J].InternationalJournalofThermalSciences, 2010, 49(1):59-66.

[8]PoulikakosD,ArcidiaconoS,MaruyamaS.Moleculardynamicssimulationinnanoscaleheattransfer:Areview[J].MicroscaleThermophysicalEngineering, 2003, 7(3):181-206.

[9]BertoliniD,TaniA.Thermalconductivityofwater:moleculardynamicsandgeneralizedhydrodynamicsresults[J].PhysicalReviewE.StatisticalPhysics,Plasmas,Fluids,andRelatedInterdisciplinaryTopics, 1997, 56(4):4135-4151.

[10]MatsumotoM.Moleculardynamicsoffluidphasechange[J].FluidPhaseEquilibria, 1998, 144(1-2):307-314.

[11]YANGChun,CHENMin,GUOZeng-yuan.Moleculardynamicssimulationoninterfacecharacteristicsofmicrodroplets[J].ChinesePhysicsLetters, 1999, 16(11):803-804.

[12]CHENGSheng-feng,LechmanJB,PlimptonSJ, et al.EvaporationofLennard-Jonesfluids[J].JournalofChemicalPhysics, 2011, 134(22):224704.

[13]YanoT.Moleculardynamicsstudyofnonequilibriumprocessesofevaporationandcondensationatavapor-liquidinterface[C] // 28thInternationalSymposiumonRarefiedGasDynamics2012.CollegePark:AmericanInstituteofPhysics, 2013.

[14]曲 伟,马同泽. 微小空间薄液膜相变传热的微尺度效应[J]. 航天器工程, 2004, 13(2):36-45.

QUWei,MATong-ze.Microscaleeffectsofheattransferduringphasechangeforthinliquidfilminmicrospace[J].SpacecraftEngineering, 2004, 13(2):36-45. (inChinese)

[15]吴方棣,郑辉东,吴燕翔. 水性质的分子动力学模拟研究[J]. 武夷学院学报, 2011, 30(5):61-66.

WUFang-di,ZHENGHui-dong,WUYan-xiang.Moleculardynamicssimulationonphysicalpropertiesofwater[J].JournalofWuyiUniversity, 2011, 30(5):61-66. (inChinese)

[16]刘娟芳,曾丹苓,刘 朝,等. 水导热系数的分子动力学模拟[J]. 工程热物理学报, 2007, 28(2):196-198.

LIUJuan-fang,ZENGDan-ling,LIUChao, et al.Moleculardynamicssimulationofconductivity[J].JournalofEngineeringThermophysics, 2007, 28(2):196-198. (inChinese)

[17]RömerF,LervikA,BresmeF.Nonequilibriummoleculardynamicssimulationsofthethermalconductivityofwater:AsystematicinvestigationoftheSPC/EandTIP4P/2005models[J].JournalofChemicalPhysics, 2012, 137(7):074503.

[18]YangTH,PanC.Moleculardynamicssimulationofathinwaterlayerevaporationandevaporationcoefficient[J].InternationalJournalofHeatandMassTransfer, 2005, 48(17):3516-3526.

[19]王发刚,李瑞阳,刘永启,等. 相变传热相关问题的分子动力学模拟及研究[J]. 上海理工大学学报, 2003, 25(2):130-138.

WANGFa-gang,LIRui-yang,LIUYong-qi, et al.Moleculardynamicssimulationandstudyonsomerelatedtopicsofphasechangeheattransfer[J].JournalofUniversityofShanghaiforScienceandTechnology, 2003, 25(2):130-138. (inChinese)

[20]LandryES,MikkilineniS,PahariaM, et al.Dropletevaporation:Amoleculardynamicsinvestigation[J].JournalofAppliedPhysics, 2007, 102(12):124301.

[21]MarooSC,ChungJN.Nanoscaleliquid-vaporphase-changephysicsinnonevaporatingregionatthethree-phasecontactline[J].JournalofAppliedPhysics, 2009, 106(6):064911.

Molecular dynamics simulation of liquid-vapor phase-change for nanoscale water films

HUANGZheng-xing*,SUNTong-tong,LIQian-qian

( School of Electronic Science and Technology, Dalian University of Technology, Dalian 116024, China )

Abstract：The influence of the thickness of nanoscale water film and the system temperature on the evaporation rate of the liquid-vapor phase-change is studied using the molecular dynamics (MD) simulation. The influence of temperature (375-425 K) on the evaporation rate of the phase-change is simulated with the water film thickness of 2 nm. In addition, at the temperature of 400 K, the influence of the water film thickness (2, 3, 4 nm) on the evaporation rate of the phase-change is simulated. The simulation results indicate that the evaporation rate increases with temperature under the condition of the same water film thickness, the evaporation rate increases with the thickness of the water film at the same temperature and exponentially decreases with time. The related results can provide a reference for the design of heat transfer devices whose principle is based on the liquid phase-change, such as micro heat pipe.

Key words：thickness of liquid films; liquid-vapor phase-change; molecular dynamics simulation

中图分类号：TK124

文献标识码：A

doi:10.7511/dllgxb201603002

作者简介:黄正兴*(1975-)，男，副教授，E-mail:huangzx@dlut.edu.cn.

基金项目:国家自然科学基金资助项目(61376115).

收稿日期：2015-12-07；修回日期: 2016-03-14．

文章编号：1000-8608(2016)03-0230-06