納米尺度水薄膜液-氣相變分子動力學模擬

黃 正 興，　孫 童 童，　李 倩 倩

( 大連理工大學 電子科學與技術學院, 遼寧 大連　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時的蒸發率最大，并且其隨溫度的升高增加幅度逐漸變?。疁囟鹊陀?00 K時，溫度升高蒸發率迅速上升；溫度高于400 K時，蒸發率隨溫度的增加上升幅度變小，達到飽和．

圖5　不同溫度條件下蒸發率隨時間的變化曲線

以上模擬情況均有一部分液體沒有變成氣態, 未蒸發的水薄膜厚度d隨溫度升高而變?。?種溫度情況下未蒸發水薄膜厚度分別為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，都有一部分水薄膜沒有發生相變，這是由于襯底硅與水分子之間存在范德華力，抑制靠近襯底硅的部分水分子蒸發．未發生相變的水薄膜厚度隨著溫度的升高而變?。?/p>

與宏觀情況相比，微尺度下液膜具有更大的蒸發率，因此對于微熱管及其他依靠內部液體相變傳熱的傳熱器件來說，增加微尺度溝道數量、增加液膜厚度和適當升高溫度能夠有效改善傳熱性能．

參考文獻：

[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