葛 宋 陳 民

(清華大學(xué)工程力學(xué)系,北京100084)

1 引言

當(dāng)熱流流過(guò)液固界面時(shí),界面處會(huì)存在溫度的不連續(xù).這種溫度的不連續(xù)是由液固兩相間的界面熱阻所引起.液固界面熱阻由Kapitza在1941年首先發(fā)現(xiàn),因此也被稱為Kapitza熱阻.1在常規(guī)條件下,界面熱阻往往可以忽略.但是近年來(lái),隨著微納米技術(shù)的發(fā)展,界面熱阻由于其在微納米尺度熱傳遞現(xiàn)象中的相對(duì)重要性而受到越來(lái)越多學(xué)者的關(guān)注.2界面熱阻的影響可由熱阻長(zhǎng)度來(lái)衡量.3實(shí)驗(yàn)測(cè)量和分子模擬的結(jié)果都表明,在液固相互作用較弱的界面,熱阻長(zhǎng)度可達(dá)數(shù)十納米.4,5因而,在微納米系統(tǒng)及含有大量液固界面的體系(如納米流體)中,界面熱阻均不可忽略.

使固體表面帶電荷及施加外電場(chǎng)是改變液固界面性質(zhì)的有效手段,有著廣泛的應(yīng)用前景.6-13朱素華等6的分子動(dòng)力學(xué)模擬研究顯示,外電場(chǎng)會(huì)影響水分子在膜材料中的電滲遷移運(yùn)動(dòng).李海蘭等7對(duì)不同電量修飾的碳納米管中水分子行為的模擬表明電荷對(duì)水分子的密度分布和擴(kuò)散速率以及水通量有顯著影響.在微納機(jī)電系統(tǒng)中,通常通過(guò)外加電場(chǎng)和通道表面的電荷來(lái)驅(qū)動(dòng)微通道內(nèi)液體流動(dòng).8Powell等9的研究顯示電場(chǎng)可以控制水分子在納米孔中的通過(guò)與阻斷.此外,通過(guò)調(diào)整施加在液體/固體間的電勢(shì)來(lái)改變兩者間接觸角的電潤(rùn)濕現(xiàn)象,10可以實(shí)現(xiàn)液體在表面上親水與疏水的可逆轉(zhuǎn)換,在顯示技術(shù)及微液滴操控和流動(dòng)控制等領(lǐng)域有廣泛的應(yīng)用.11-13表面電荷和電場(chǎng)能改變界面以及水分子體系本身的一些性質(zhì),14,15有可能對(duì)液固界面間的能量輸運(yùn)產(chǎn)生影響.Lee等16的研究表明納米顆粒的表面電荷對(duì)納米流體的導(dǎo)熱系數(shù)有重要影響,并認(rèn)為表面電荷能增強(qiáng)液固間的熱輸運(yùn).但目前對(duì)于表面電荷和電場(chǎng)對(duì)界面熱阻的影響還沒(méi)有直接的報(bào)道.本文利用分子動(dòng)力學(xué)方法模擬了固體表面電荷和外加電場(chǎng)對(duì)水與固體間界面熱阻的影響,在加深對(duì)界面熱阻影響機(jī)制理解的同時(shí),也能為實(shí)際應(yīng)用提供有益參考.

2 模擬方法與細(xì)節(jié)

本文所有模擬都采用了分子動(dòng)力學(xué)模擬軟件Large-Scale Atomic/MolecularMassively Parallel Simulator(LAMMPS).17模擬中,水分子采用常用的extended simple point charge(SPC/E)模型.18,19水分子間的相互作用可表述為

式(1)中前半部分為電荷間的靜電作用,后半部分為氧原子間的Lennard-Jones(LJ)作用.其中qα、qβ分別表示不同水分子中的點(diǎn)電荷,rαβ為兩電荷間的距離. εOO、σOO分別表示氧原子間LJ作用的能量參數(shù)和長(zhǎng)度參數(shù),rOO為氧原子間的距離.固定SPC/E水分子模型假設(shè)只有氧原子間存在LJ作用,且水分子的構(gòu)型在模擬中保持不變,即鍵長(zhǎng)和鍵角固定.固體采用以銅為原型的LJ固體,其勢(shì)能參數(shù)為ε=6.64× 10-20J,σ=2.277×10-10m.20固體原子按面心立方(face-centered cubic)結(jié)構(gòu)排列,晶格常數(shù)取a0= 1.56σ.21由于LJ固體熔點(diǎn)正比于能量參數(shù)ε,在本文模擬的溫度下此模型能很好地維持固體狀態(tài)而無(wú)需人為引入彈簧力來(lái)維持,并且能很好地描述固體的熱學(xué)性質(zhì).5水分子與固體間的相互作用由氧原子與固體原子間的LJ勢(shì)來(lái)描述,勢(shì)能參數(shù)εSO=2.78× 10-21J,σSO=(σOO+σSS)/2＝2.782×10-10m.

液固界面熱阻的模擬體系如圖1所示,由液體區(qū)域和固體區(qū)域兩部分構(gòu)成.模擬盒子的尺寸為15a0×15a0×32a0,其中液體區(qū)域尺寸為15a0×15a0× 20a0.水分子的密度設(shè)為994 kg·m-3.體系共含有6700個(gè)水分子和11700個(gè)固體原子.如圖1所示,在沿z方向的兩端分別設(shè)定兩個(gè)固定的原子層,以維持模擬體系穩(wěn)定.在緊鄰固定原子層的液體和固體圖示區(qū)域分別設(shè)置熱源和熱沉,以此在體系中形成熱流和溫度梯度.液體中熱源區(qū)域和固體中熱沉區(qū)域的大小均為15a0×15a0×a0.在x和y方向均采用周期性邊界條件.模擬中為了提高效率,只有當(dāng)原子間的距離小于截?cái)喟霃?1.0 nm)時(shí)才計(jì)算LJ相互作用.而庫(kù)侖力為長(zhǎng)程相互作用,在本文中采用particle-particle particle-mesh(PPPM)方法22來(lái)處理.模擬中采用SHAKE算法23來(lái)固定水分子的鍵長(zhǎng)和鍵角.粒子運(yùn)動(dòng)方程的積分采用速度Verlet算法24,積分步長(zhǎng)取為0.5 fs.模擬首先采用NVT系綜,利用Nose-Hoover控溫方法25使體系在300 K下演化50萬(wàn)步,隨后將控溫移除,使體系在NVE系綜中演化.在熱源中不斷輸入熱量,在熱沉中移走相等的熱量,其具體做法是在每一步先通過(guò)坐標(biāo)來(lái)確定所有處于熱源和熱沉中的原子,然后通過(guò)重新調(diào)整熱源和熱沉區(qū)域中原子的熱運(yùn)動(dòng)速度來(lái)分別輸入和移出能量.26體系整體將保持能量守恒.每一步加入和移出熱量ΔE,體系中的熱流將為J=ΔE/AΔt(其中A為界面的橫截面積,Δt為時(shí)間步長(zhǎng)).此過(guò)程持續(xù)80萬(wàn)步令體系達(dá)到穩(wěn)定,體系中建立起穩(wěn)定的溫度分布,界面溫度跳躍ΔT及熱流J.界面溫度跳躍ΔT定義為界面處液體層和固體層的溫度差.繼續(xù)計(jì)算200萬(wàn)步以完成統(tǒng)計(jì)采樣.將體系沿z方向劃分成條狀區(qū)域,其厚度為a0,每個(gè)條狀區(qū)域中含有超過(guò)300個(gè)原子以減小溫度統(tǒng)計(jì)的漲落.將總的模擬時(shí)間分成10個(gè)時(shí)間塊,利用塊平均的方法27來(lái)確定溫度的平均值及誤差.統(tǒng)計(jì)各區(qū)域中的溫度可獲得溫度分布和界面溫度跳躍.液固界面熱阻RK則可由其定義RK=ΔT/J來(lái)獲得.

圖1 液固界面熱阻的模擬體系Fig.1 Simulation system for the liquid/solid interfacial thermal resistance

3 結(jié)果與討論

首先我們討論無(wú)表面電荷和外電場(chǎng)的情況.模擬需要先設(shè)定體系的熱流密度.設(shè)定的熱流密度不能太大以保證體系維持在線性響應(yīng)區(qū)域,又要保證能產(chǎn)生足夠強(qiáng)的溫度信號(hào).Hu等28對(duì)水與硅膠界面熱阻的模擬結(jié)果顯示在熱流較大(＞6000 MW·m-2)時(shí),界面熱阻幾乎與設(shè)定的熱流大小無(wú)關(guān).因此本文中采用的熱流為J=7350 MW·m-2.穩(wěn)定后體系內(nèi)的溫度分布如圖2所示.在給定的熱流密度下,穩(wěn)定后固體和液體體系內(nèi)的溫度分布均較好地滿足線性分布,體系處于線性響應(yīng)區(qū),傅里葉導(dǎo)熱定律仍然成立.可以看到液體中的溫度梯度要明顯大于固體中的溫度梯度,這是由于固體熱導(dǎo)率遠(yuǎn)大于液體熱導(dǎo)率.液體與固體在界面處存在顯著的溫度不連續(xù).由熱流密度和界面溫度跳躍計(jì)算得到固體與水之間的界面熱阻RK=3.73×10-9K·m2·W-1,與之相對(duì)應(yīng)的界面熱導(dǎo)(界面熱阻的倒數(shù))為268 MW·m-2· K-1.此界面熱導(dǎo)與Ge等29實(shí)驗(yàn)測(cè)量得到的水與AuPd雙金屬納米粒子間的界面熱導(dǎo)(100-300 MW·m-2·K-1)處于同一量級(jí).

圖2 穩(wěn)定后體系內(nèi)的溫度分布Fig.2 Temperature distribution of the system after equilibrium

為考慮固體表面電荷對(duì)界面熱阻的影響,模擬中設(shè)定與水分子接觸的最外層原子分別帶0.01e-0.05e的正負(fù)電荷,對(duì)應(yīng)的表面電荷密度為1.6-8.0 mC·cm-2.體系中除增加了水分子上的電荷與表面電荷間的靜電作用外,其它勢(shì)能均保持不變.圖3所示為表面電荷對(duì)界面熱阻的影響.采用簡(jiǎn)單的函數(shù)形式對(duì)此關(guān)系進(jìn)行擬合,發(fā)現(xiàn)二次函數(shù)能較好地描述表面電荷密度與界面熱阻的關(guān)系,擬合曲線如圖3所示.

可以看到,無(wú)論表面帶正電荷還是負(fù)電荷,界面熱阻都隨著表面電荷密度的增加而減小.表面電荷可有效減小界面熱阻,促進(jìn)界面處的能量輸運(yùn).

當(dāng)考慮外加電場(chǎng)對(duì)界面熱阻的影響時(shí),由于水分子模型中點(diǎn)電荷的存在,電荷將受到電場(chǎng)力F=qE的作用,其中q為點(diǎn)電荷電量,E為電場(chǎng)強(qiáng)度.我們考慮電場(chǎng)的兩種不同布置方式,即垂直于界面與平行于界面,其中垂直于界面還可以分為垂直于界面指向液體和垂直于界面指向固體.為了表述方便,我們將指向液體記為正方向,指向固體記為負(fù)方向.圖4所示為電場(chǎng)垂直于界面時(shí)對(duì)界面熱阻的影響.同樣采用二次函數(shù)擬合能較好地描述電場(chǎng)強(qiáng)度與界面熱阻的關(guān)系,擬合曲線如圖4所示.

圖3 界面熱阻隨表面電荷的變化Fig.3 Interfacial thermal resistance as a function of surface charge

由圖4可見(jiàn),電場(chǎng)垂直于界面時(shí),無(wú)論電場(chǎng)方向指向液體還是固體,界面熱阻都隨著電場(chǎng)強(qiáng)度的增加而減小.電場(chǎng)強(qiáng)度與界面熱阻的關(guān)系也同樣能通過(guò)二次函數(shù)來(lái)較好地?cái)M合,電場(chǎng)和表面電荷對(duì)界面熱阻的影響存在類似的效果.值得一提的是,在電潤(rùn)濕現(xiàn)象中,接觸角的余弦與外加電壓之間也存在二次函數(shù)的關(guān)系.10由于界面熱阻和表面濕潤(rùn)性間存在一定的聯(lián)系,30電場(chǎng)和表面電荷對(duì)界面熱阻和濕潤(rùn)性兩種界面現(xiàn)象的影響可能存在相似的機(jī)理,還有待進(jìn)一步探究.

當(dāng)電場(chǎng)平行于界面方向時(shí)界面熱阻的變化情況則有所不同.圖5所示為電場(chǎng)平行于界面時(shí)界面熱阻的變化,界面熱阻幾乎不受外加電場(chǎng)的影響.

圖4 電場(chǎng)垂直于界面時(shí)界面熱阻隨電場(chǎng)強(qiáng)度的變化Fig.4 Interfacial thermal resistance as a function of electric field strength when the electric field is perpendicular to the interface

圖5 電場(chǎng)平行于界面時(shí)界面熱阻隨電場(chǎng)強(qiáng)度的變化Fig.5 Interfacial thermal resistance as a function of electric field strength when the electric field is parallel to the interface

由上面的結(jié)果可以看出通過(guò)施加表面電荷和特定方向上的外電場(chǎng)確實(shí)可以改變液固間的界面熱阻.由于水分子本身帶有等量的正負(fù)電荷,呈電中性,在外電場(chǎng)或電荷的作用下整體受力為零,因此不會(huì)產(chǎn)生水分子在界面處的聚集等現(xiàn)象.在靜電力不是非常強(qiáng)的情況下也不會(huì)對(duì)水分子的整體結(jié)構(gòu)性質(zhì)(如長(zhǎng)程有序度等)產(chǎn)生影響.我們分析不同電場(chǎng)強(qiáng)度下液體水中氧原子的徑向分布函數(shù)gOO(r) (徑向分布函數(shù)表征的是距離目標(biāo)原子為r的球殼內(nèi)另一種原子的密度與平均密度的比值,其可以反映液體體系的結(jié)構(gòu)特征).如圖6所示,電場(chǎng)對(duì)氧氧原子的徑向分布函數(shù)幾乎沒(méi)有產(chǎn)生可見(jiàn)的影響,這也與Evans等14的結(jié)論一致.因此可排除界面水分子聚集和長(zhǎng)程結(jié)構(gòu)變化的影響.考慮到水分子雖然呈電中性,但由于水分子的結(jié)構(gòu)特征,水分子中存在偶極矩,從而使得水分子的取向會(huì)受到表面電荷和外電場(chǎng)的影響.在無(wú)表面電荷和外電場(chǎng)的條件下,由于水分子取向的隨機(jī)排列,水分子體系的平均偶極矩為零,而正負(fù)表面電荷或垂直于界面的正負(fù)電場(chǎng)將導(dǎo)致水分子垂直于界面向上或向下排列.此現(xiàn)象已得到了模擬的驗(yàn)證,如Evans等14分析了電場(chǎng)作用下水分子的結(jié)構(gòu),表明水分子會(huì)發(fā)生沿電場(chǎng)方向的定向排列,平均偶極矩會(huì)隨電場(chǎng)的增大而增加.考慮液固間能量傳遞的機(jī)理,液體中的能量傳遞與固體有所不同,由于液體不能承受剪切,液體中的能量傳遞主要以縱向聲波為主.1水分子在電場(chǎng)和表面電荷影響下的排列方式可能會(huì)明顯增強(qiáng)液體水中的縱向聲波的能量輸運(yùn),從而促進(jìn)水分子與固體間的能量交換.而平行于界面施加的電場(chǎng)引起的偶極矩排列可能對(duì)熱量輸運(yùn)不產(chǎn)生明顯影響.但其中的機(jī)理還有待進(jìn)一步研究.

圖6 垂直界面施加電場(chǎng)時(shí)水分子中氧氧徑向分布函數(shù)Fig.6 Radial distribution function for oxygen atoms in water when the external electric field is perpendicular to the interface

4 結(jié)論

采用非平衡分子動(dòng)力學(xué)模擬分析了固體表面電荷及外電場(chǎng)對(duì)水與固體間界面熱阻的影響.模擬結(jié)果表明,壁面帶正電荷或負(fù)電荷都將使得界面熱阻減小,且界面熱阻與表面電荷密度滿足二次函數(shù)關(guān)系;外加電場(chǎng)平行于界面時(shí),其對(duì)界面熱阻幾乎沒(méi)有產(chǎn)生可見(jiàn)的影響,而垂直于界面時(shí),界面熱阻將隨著電場(chǎng)強(qiáng)度的增大而減小.界面熱阻與電場(chǎng)強(qiáng)度也滿足二次函數(shù)關(guān)系.界面熱阻對(duì)表面電荷和電場(chǎng)的影響可能與水分子偶極矩的定向排列有關(guān).本文的模擬結(jié)果表明改變電場(chǎng)和施加表面電荷是控制界面熱阻的有效方法.

(1) Swartz,E.T.;Pohl,R.O.Rev.Mod.Phys.1989,61,605.doi: 10.1103/RevModPhys.61.605

(2)Cahill,D.G.;Ford,W.K.;Goodson,K.E.;Mahan G.D.; Majumdar,A.;Maris,H.J.;Merlin,R.;Phillpot,S.R.J.Appl. Phys.2003,93,793.doi:10.1063/1.1524305

(3) Barrat,J.L.;Chiaruttini,F.Mol.Phys.2003,101,1605.doi: 10.1080/0026897031000068578

(4) Ge,Z.B.;Cahill,D.G.;Braun,V.Phys.Rev.Lett.2006,96, 186101.doi:10.1103/PhysRevLett.96.186101

(5) Xue,L.;Keblinski,P.;Phillipot,S.R.;Choi,S.U.S.;Eastman, J.A.J.Chem.Phys.2003,118,337.doi:10.1063/1.1525806

(6) Zhu,S.H.;Yan,L.M.;Ji,X.B.;Shao,C.L.;Lu,W.C.Acta Phys.-Chim.Sin.2010,26,2659. [朱素華,嚴(yán)六明,紀(jì)曉波,邵長(zhǎng)樂(lè),陸文聰.物理化學(xué)學(xué)報(bào),2010,26,2659.]doi:10.3866/ PKU.WHXB20100934

(7) Li,H.L.;Jia,Y.X.;Hu,Y.D.Acta Phys.-Chim.Sin.2012,28, 573.[李海蘭,賈玉香,胡仰棟.物理化學(xué)學(xué)報(bào),2012,28, 573.]doi:10.3866/PKU.WHXB201112191

(8) Li,Q.L.;Chen,L.X.Prog.Chem.2008,20,1406.[李清嶺,陳令新.化學(xué)進(jìn)展,2008,20,1406.]

(9) Powell,M.R.;Cleary,L.;Davenport,M.;Shea,K.J.;Siwy,Z. S.Nat.Nanotech.2011,6,798.doi:10.1038/nnano.2011.189

(10) Mugele,F.;Baret,J.C.J.Phys.:Condes.Matter 2005,17,R705.

(11) Hayes,R.A.;Feenstra,B.J.Nature 2003,425,383.doi: 10.1038/nature01988

(12) Bhushan,B.;Pan,Y.Langmuir 2011,27,9425.

(13) Wang,Y.L.;Bhushan,B.Langmuir 2010,26,4013.doi: 10.1021/la903460a

(14) Evans,W.;Fish,J.;Keblinski,P.J.Chem.Phys.2007,126, 154504.doi:10.1063/1.2723071

(15) Xu,G.L.;Xiao,X.H.;Liu,Y.F.;Sun,J.F.;Zhu,Z.H.Acta Phys.-Chim.Sin.2007,23,746. [徐國(guó)亮,肖小紅,劉玉芳,孫金鋒,朱正和.物理化學(xué)學(xué)報(bào),2007,23,746.]doi:10.3866/ PKU.WHXB20070524

(16)Lee,D.;Kim,J.W.;Kim,B.G.J.Phys.Chem.B 2006,110, 4323.doi:10.1021/jp057225m

(17) Plimpton,S.J.Comp.Phys.1995,117,1.doi:10.1006/jcph. 1995.1039

(18)Ge,S.;Zhang,X.X.;Chen,M.J.Chem.Eng.Data 2011,56, 1299.doi:10.1021/je101138g

(19) Ge,S.;Zhang,X.X.;Chen,M.Int.J.Thermophys.2010,31, 2176.doi:10.1007/s10765-010-0863-0

(20)Agrawal,P.M.;Rice,B.M.;Thompson,D.L.Surf.Sci.2002, 515,21.doi:10.1016/S0039-6028(02)01916-7

(21) Chen,G.Nanoscale Energy Transport and Conversion;Oxford University Press:New York,2005;pp 85-86.

(22) Deserno,M.;Holm,C.J.Chem.Phys.1998,109(18),7678. doi:10.1063/1.477414

(23) Ryckaert,J.P.;Ciccotti,G.;Berendsen,H.J.C.J.Comp.Phys. 1977,23,327.

(24)Allen,M.P.;Tildesley,D.J.Computer Simulation of Liquids; Clarendon Press:Oxford,1997;pp 87-89.

(25) Evans,D.J.;Holian,B.L.J.Chem.Phys.1985,83(8),4069. doi:10.1063/1.449071

(26) Jund,P.;Jullien,R.Phys.Rev.B 1999,59(21),13707.doi: 10.1103/PhysRevB.59.13707

(27) Frenkel,D.;Smit,B.Understanding Molecular Simulation: From Algorithms to Application;Academic Press:New York, 1996;pp 529-532.

(28) Hu,M.;Goicochea,J.V.;Michel,B.;Poulikakos,D.Appl.Phys. Lett.2009,95,151903.doi:10.1063/1.3247882

(29) Ge,Z.B.;Cahill,D.G.;Braun,P.V.J.Phys.Chem.B 2004, 108,18870.doi:10.1021/jp048375k

(30) Shenogina,N.;Godawat,R.;Keblinski,P.;Garde,S.Phys.Rev. Lett.2009,102,156101.doi:10.1103/PhysRevLett.102.156101