TiO2/ZnO納米薄膜界面熱導率的分子動力學模擬*

楊平吳勇勝許海鋒許鮮欣張立強李培

TiO2/ZnO納米薄膜界面熱導率的分子動力學模擬*

楊平吳勇勝 許海鋒 許鮮欣 張立強 李培

(江蘇大學機械工程學院，鎮江212013)
(2010年7月5日收到;2010年9月28日收到修改稿)

采用平衡分子動力學方法及Buckingham勢研究了金紅石型TiO2薄膜與閃鋅礦型ZnO薄膜構筑的納米薄膜界面沿晶面［0001］(z軸方向)的熱導率．通過優化分子模擬初始條件中的截斷半徑rc和時間步后，計算并分析了平衡溫度、薄膜厚度、薄膜截面大小對熱導率的影響．研究表明，薄膜熱導率受薄膜溫度和厚度的影響很大，當溫度由300 K升高600 K時，薄膜的熱導率逐漸減小;當薄膜厚度由1.8 nm增大到5 nm時，熱導率會逐漸增大;并在此基礎上對更大厚度范圍內TiO2/ZnO界面熱導率可能變化性態進行了探討．

熱導率，分子動力學，TiO2/ZnO納米薄膜界面，數值模擬

PACS:66.70.Df，83.10.Rs，68.35.－p，02.60.－x

1.引言

近年來，納米科學特別是納米材料的制備和性能研究引起研究者的廣泛關注．這是由于當材料小到納米尺度時，處于材料邊界處的原子在組成整個納米材料的原子中占有很大的比例，這樣的結果就體現為納米材料的尺寸效應．另一方面，界面作為微納制造中最普通的結構，對材料的力學、電學、熱傳導等物理和化學特性具有非常重要的影響．在界面附近，材料的物理和化學性質有可能發生很大的變化，這就體現為納米材料的界面效應．因此，界面結構及界面效應研究的重要性是不能忽略的，界面科學在材料和物理研究領域中已成為一個新的分支［1］．由于尺寸效應和界面效應的影響;納米薄膜界面傳熱等相關問題甚至需要從分子/原子的層面來進行研究［2］．而分子動力學方法(molecular dynamics，簡記為MD)就是通過求解有相互作用的各個粒子的運動方程，得到每個粒子的空間位置、運動狀態隨時間的演化狀況，從而統計出材料的宏觀行為特性的方法．前期的研究顯示了MD在研究微/納尺度傳熱特性方面的可行性;如吳國強等［3］采用平衡分子動力學方法(EMD)和非平衡分子動力學方法(NEMD)模擬了氬晶體薄膜的法向熱導率，發現熱導率隨著薄膜厚度的增加而呈線性增加，賈明等［4］則采用改進分析型嵌入原子模型(MAEAM)，通過MD來模擬了Mo薄膜的熱力學性質．侯泉文等［5］通過NEMD研究了碳納米管熱導率與長度的關系，發現碳納米管熱導率隨其長度的增加而增加;并詳細研究了導熱中彈道運輸階段、擴散運輸階段與碳納米管長度的關系．Kulkarni和Zhou［6］利用MD模擬了長度從19—41(1= 0.1 nm)的ZnO納米帶在溫度范圍為500—1500 K，沿［0110］方向上的熱導率，研究顯示ZnO納米帶具有明顯的尺寸效應和溫度效應．Hegedus［7］應用NEMD討論了異質系統界面熱傳導，并與基于聲子散射理論的界面熱阻模型進行了比較，兩者取得了良好匹配．Liang和Sun［8］采用NEMD發現，界面熱阻在納米界面結構熱傳遞中非常重要，并詳細討論了界面特性對熱導率的影響．Volz等［9］采用NEMD模擬研究了界面、晶格應變對整個Si/Ge超晶格熱導率的影響，并且模擬結果與實驗數據相符;另外，Yang等［10，11］則采用MD與有限元結合的方法研究了兩種材料所形成的界面對熱傳輸的影響．

過渡金屬氧化物因其具有良好的光電磁等特性而成為了當今研究的熱點．ZnO和TiO2作為過渡金屬氧化物中很常見的兩種物質，在太陽電池和光催化劑中有著廣泛的應用，但對這兩種物質特別是兩者構筑的界面結構在納米尺度上的傳熱特性和機理還缺乏足夠的了解．本文采用EMD及優化的Buckingham勢函數來研究金紅石型TiO2與閃鋅礦型ZnO構筑的納米薄膜界面特性，研究在xy平面上晶胞數目、整個薄膜系統的厚度以及系統溫度與熱導率的關系．

2 ．建模與計算

在建模時，首先必須設置原子的初始位型，原子的初始位型包括在t0時刻速度以及其他運動信息，這些運動屬性可描述為

原子i上受到的合力為MD勢函數負的梯度，即

通過數值求解牛頓方程得到一個時間序列，即可得到任意物理量的系綜平均為

因此，從計算數學的角度來看，MD是微分方程的數值解問題．

對于系統溫度保持不變的正則系統(NVT)，采用Nose-Hoover熱浴進行溫度調節，描述該系統的運動方程為

其中mi，ri，fi，pi和U(r1，r2，…，rn)分別是第i個粒子的質量、坐標、所受到的合力、動量和其他粒子對其的作用勢，ξ是約束溫度不變的參數．由于溫度不變，系統動能=0的關系，所以參數

2.1.局域溫度的量子修正

在經典分子動力學中，系統溫度由Boltzmann能量均分定理得到

其中N為模擬區域的粒子總數，kB為Boltzmann常數，TMD為系統模擬溫度，mi和vi分別為粒子的質量和速度．當溫度低于材料Debye溫度時，由于系統的比熱容與溫度有關，(8)式不再成立，即系統實際溫度T≠TMD．為了獲得正確的仿真結果，必須對模擬溫度TMD、模擬結果kMD進行量子化修正．采用Debye模型進行溫度統計的量子修正，通過求解下列積分方程進行修正

求解此方程可以得到與MD的局域溫度Tj，MD相對應的真實晶格溫度Tj，而求解方程(9)須引入晶格格波譜，即格波波矢q與頻率ω的關系．本文采用近似的Debye模型

Debye模型給出的晶格熱容為

其中C'v為質量定容熱容，ΘD為Debye溫度．在Debye近似下，能量方程(9)轉化為

本文通過求解方程(12)，實現了對MD得到的局域溫度的量子修正．

2.2 ．EMD中熱導率的計算

EMD意味著通過模擬系統的平衡態求出導熱系數．從溫度梯度加上到形成穩態的過程中，熱流的變化有可采用如下公式描述［12］:

其中τ是聲子的弛豫時間，q是穩態的熱流，qz是到達穩態前在z軸方向隨時間變化的熱流，且

其中，εi是包含勢能和動能的總能量．(14)式只適用于兩體勢．

采用平衡態方法時，沒有宏觀的溫度梯度，所以q=0．根據G-K公式，熱導率可以由熱流qz(t)的自相關函數積分得到［13］

其中kMD為熱導率，T為溫度，V為模擬區域的體積．

2.3.原子之間的作用勢函數

在MD模擬中，材料系統內粒子間相互作用是通過勢函數來實現，勢函數的選擇直接關系到模擬結果的準確性．氧化物與氧化物組成的薄膜界面原子之間的作用力，還沒有成熟的勢函數來描述［14］．本文中的TiO2/ZnO薄膜界面系統，原子作用對存在10種類型．這里將采用Buckingham勢函數的兩體勢形式來描述TiO2及ZnO薄膜內部原子之間的作用力，Buckingham函數形式［15］如下:

其中Eij為原子i和j的相互作用能．對于ZnO和TiO2中原子的A，C，ρ的取值見文獻［15，16］;對于ZnO薄膜和TiO2薄膜界面上的原子作用對如Ti-Zn，Ti-OZn，Zn-OTi，OZn-OTi，本文將采用文獻［17，18］中方法來處理．

2.4.模型的建立

本文研究的TiO2，ZnO的晶體類型分別為金紅石型和閃鋅礦型，晶體結構如圖1所示，它們的參數如下:TiO2的參數為密度ρ=4.25×103kg/m3，晶胞參數a=4.59，c=2.96;Zn O的參數為密度ρ=5.45×103kg/m3，晶胞參數a=4.63．

模型建立過程如下:本文模擬的是TiO2/ZnO沿晶面［0001］的熱導率，在xy平面上將采用4×4個晶胞，且xy平面平行于晶面［0001］，在z方向上對TiO2/ZnO薄膜界面建立了1.8，2.9，3.9，5 nm的4個厚度．建立的模型如圖2所示．

圖1 晶體結構示意圖(a)TiO2晶體晶胞結構;(b)ZnO晶體晶胞結構

圖2 TiO2/ZnO薄膜界面模型

2.5.MD的初始條件

在x，y方向上采用周期性邊界條件，z軸為熱流方向，采用Nose-Hoover熱浴進行溫度調節，模擬在正則系統(NVT)中進行，采用Gear算法進行運動方程積分，庫侖力的計算采用Ewald方法，時間步長為0.1 fs，TiO2中原子Ti和O的價位數分別為+2.196和－1.098，ZnO中原子Zn和O的價位數分別為+2和－2．

MD初始條件中的另外兩個參數截斷半徑rc和模擬時間步數是兩個很重要的初始條件，對模擬結果和計算時間有著直接的影響．為了保證模擬結果的精度和盡量減少計算時間，下面通過計算來確定這兩個參數的大小．

2.5.1 ．模擬時間的確定

由于文中是通過EMD來計算1.8，2.9，3.9，5 nm這4個厚度的薄膜在溫度為300，400，500，600 K下的熱導率，所以可以通過系統平衡溫度與模擬時間的關系來確定MD模擬時間的大小，即當系統溫度平衡時所對應的時間即為我們所需要的MD初始時間．圖3中(a)—(d)中的曲線分別表示的是TiO2/ZnO薄膜界面厚度為3.9 nm時，在模擬溫度為300，400，500，600 K下薄膜界面的平衡溫度與模擬時間的關系，此處分子模擬的初始條件中的rc= 10，模擬時間為20 ps(時間步長為0.1 fs，步數為2×105)．從圖3可以看出，不管在哪個模擬溫度下，系統溫度平衡的時間大約都在10 ps(1×105步)左右，在10 ps后，系統溫度基本不變，即整個薄膜系統達到平衡．因此，溫度對平衡時間沒有太大影響，我們就把10 ps(1×105步)當做MD的初始時間．

2.5.2 ．截斷半徑rc的確定

在確定模擬時間為10 ps的情況下，再來考慮rc大小的選擇，為了消除其他因素的影響，此處選擇的薄膜系統及初始條件與上述相同，只是此時的變量為rc．計算結果如圖4所示，模擬結果顯示出同一薄膜系統不同溫度下rc的選取對模擬結果影響各不相同，如模擬溫度為300 K時，可選取0.4 nm及以上的長度為rc，而600 K時則要選取0.6 nm及以上的長度為rc才能保證模擬結果的精度．但不管在哪個模擬溫度下，當rc為0.6 nm時，計算出來的熱導率已具有足夠的準確性．為了保證模擬結果的精度和減少計算時間以及滿足不同的薄膜系統，本文將設rc=0.6 nm作為MD的初始條件．

圖3 薄膜界面平衡溫度與時間關系曲線(a)T=300 K;(b)T=400 K;(c)T=500 K;(d)T=600 K

圖4 熱導率隨著截斷半徑的變化

由于上述兩個參數在確定過程中只考慮了溫度因素，而本文在模擬過程中存在4種不同厚度，因此在進行計算前先對上述確定的兩個參數進行論證以確定它們的正確性．圖5計算的是在系統溫度為400 K，TiO2/ZnO薄膜界面厚度為1.8，2.9，3.9， 5 nm情況下熱導率與時間的關系曲線．從圖5可以看出，不管哪一個厚度下，當模擬時間大于10 ps時，熱導率的變化已經很小．這與上面用溫度來確定的模擬時間步是一致的，并且時間步數的選擇與薄膜界面厚度沒有太大的關系．因此，上述確定的兩個參數是可靠的．

3 ．模擬結果與討論

3.1 ．TiO2/ZnO薄膜界面水平方向尺寸與熱導率的關系

圖6顯示的是TiO2/ZnO薄膜界面在厚度為5 nm時水平方向尺寸對熱導的影響．從圖6可以看出，隨著水平尺寸的增大(原子數的增多)，薄膜界面熱導率也隨之增大．這說明薄膜在水平方向也具有明顯的尺寸效應．可能的原因是水平尺寸越大，薄膜的熱阻越小，而熱導率與熱阻之間存在反比關系，因而水平尺寸越大，熱導率也越大．這與文獻［19］的結果相符．另外，從圖6還可以看出，原子數增多時，計算結果的波動越來越小，即精度得到很大的提高．這是因為截面積的增加，模擬的原子數相應增加，這就更接近于實際情況，因而模擬得到的結果也就更接近．由此可知，要想得到精度高的模擬結果，選擇合適的水平尺寸是必要的，這也是本文中選擇4×4個原胞進行模擬的原因．

圖5 不同厚度下熱導率與模擬時間的關系

圖6 水平方向模型尺寸對熱導率的影響

3.2 ．TiO2/Zn O薄膜界面厚度及溫度與熱導率的關系

圖7顯示了不同溫度下TiO2/ZnO薄膜界面沿z軸方向的尺寸對薄膜熱導率的影響．從圖7可以看出，TiO2/ZnO薄膜界面熱導率隨厚度的增加而增加，并且隨著厚度增加，熱導率增加趨勢變得平緩．模擬得到的界面熱導率處于1.95—2.78 W/mk的范圍內，這不僅遠小于對應的單層TiO2薄膜或ZnO薄膜的熱導率(TiO2薄膜熱導率為4 W·m－1·K－1左右［20］，ZnO薄膜熱導率為10 W·m－1·K－1左右［6］)，而且也小于這兩種薄膜材料熱導平均值的一半，這一結論與文獻［21］相符．這種薄膜尺寸效應以及小于兩種薄膜熱導率平均值一半的原因可作如下解釋．

圖7 薄膜界面厚度對熱導率的影響

對于半導體材料而言，聲子是主要的熱載流子．組成TiO2/Zn O薄膜界面的單層薄膜中，根據固體物理理論中聲子氣的動力論，電介質固體的晶格導熱系數可表示為

其中v為聲速，l為聲子的平均自由程．經典固體物理理論中認為:尺寸效應一般在薄膜厚度小于聲子平均自由程的情況下發生．在本文模擬的溫度范圍，聲子平均自由程可近似表示為

其中γ為Gruneisen常數，通常γ≈2，Tm為材料熔點溫度，a為晶格常數．對于TiO2和Zn O兩種材料而言，通過(18)式可以估算出它們的平均聲子自由程為10 nm以上．這要大于本文研究的薄膜厚度．因而整個薄膜界面顯示出尺寸效應．

另外，在組成TiO2/Zn O薄膜界面中的單層薄膜內部，薄膜邊界的散射對熱導率也有很大影響，根據晶格導熱系數理論，可得到晶格導熱系數為

其中n為原子的數密度，τ為聲子平均弛豫時間．由于薄膜邊界對聲子的散射，薄膜中有效的聲子平均弛豫時間可根據Matthiessen法則確定，

其中τb為塊體材料中的聲子平均弛豫時間，p為薄膜邊界對聲子的鏡反射率，d為薄膜的厚度．根據(19)和(20)式可以得出，d和p越大，有效弛豫時間τ也越大，從而薄膜熱導率就越大，反之就越小．這就是熱導率隨著薄膜厚度的增加而變大的原因．

對于TiO2/ZnO薄膜界面而言，聲子在傳熱過程中除了在薄膜內部和邊界發生散射外，界面處對聲子的影響也非常大．聲失配理論認為，聲子像平面波一樣是連續體，當通過兩種不同的材料時聲子會在界面處發生反射、散射以及模轉化，從而導致聲子傳熱能力的降低．因此，整個薄膜界面的熱導率不但小于單層薄膜的值，而且要小于兩層薄膜熱導率平均值的一半．

圖8 TiO2/ZnO薄膜界面熱導率與溫度的關系

圖8 所示為TiO2/ZnO2薄膜界面熱導率與溫度的關系．從圖8可看出，當溫度由300 K上升到600 K時，熱導率降低了15%左右，顯示出薄膜熱導率具有明顯的溫度效應．這主要是由于溫度升高時，材料的熱軟化以及原子振動幅度變大，以致高頻聲子在熱傳遞中起的作用越來越大，而高頻聲子的散射更強，聲子自由程更短，從而導致整個薄膜界面熱導率降低．據文獻［20］中對TiO2單層薄膜和文獻［6］中對ZnO單層薄膜熱導率的實驗研究可知，在一定的溫度范圍內，單層TiO2薄膜和單層ZnO薄膜熱導率隨著溫度的增加而減小．本文的研究結果也表明，TiO2/Zn O薄膜界面熱導率隨溫度的變化趨勢和其單層薄膜隨溫度的變化趨勢具有相似性．

4 ．結論與討論

通過計算，選擇合理的截斷半徑rc=0.6 nm和模擬時間t=10 ps，用EMD對TiO2/ZnO薄膜界面熱導率進行了模擬．模擬結果表明:薄膜界面熱導率在水平方向和沿晶面［0001］方向都存在明顯的尺寸效應，在水平方向上晶胞數位4×4時的熱導率幾乎為1×1晶胞熱導率的2.5倍．晶面［0001］方向上，薄膜界面厚度為5 nm時的熱導率比1.9 nm時高出約20%，并且整個TiO2/ZnO薄膜系統熱導率值要小于這兩種材料熱導率平均值的一半．其次，晶面［0001］方向的熱導率還具有明顯的溫度效應，熱導率隨著溫度的升高而減小，減小幅度達到15%．

對于更大范圍的TiO2/Zn O薄膜界面厚度與熱導率的的關系，用于計算容量的限制，本文未作進一步的研究，但進行以下兩種推測:1)隨著TiO2/ZnO薄膜界面厚度的增加，其熱導率也隨著增加，最終直至宏觀TiO2與宏觀ZnO所成界面的熱導率．2)隨著薄膜厚度的增加，熱導率變化有不可預測性，如出現反向非線性等．若為第二種推測，則說明TiO2/Zn O薄膜界面熱導率具有特殊性．這一特殊性是普遍存在于氧化物界面之間，還是只有TiO2/Zn O薄膜界面所特有?這有待進一步探索．

［1］Allara D L 2005 Nature 437 638

［2］Chou F C，Lukes J R，Liang X G 1999 J．Heat Transfer．10 141

［3］Wu G Q，Kong X R，Sun Z W，Wang Y H 2006 Acta Phys．Sin．55 1(in Chinese)［吳國強、孔憲仁、孫兆偉、王亞輝2006物理學報55 1］

［4］Jia M，Lai Y Q，Tian Z L，Liu Y X 2009 Acta Phys．Sin．58 1139(in Chinese)［賈明、賴延清、田忠良、劉業翔2009物理學報58 1139］

［5］Hou Q W，Cao B Y，Guo Z Y 2009 Acta Phys．Sin．58 7809 (in Chinese)［侯泉文、曹炳陽、過增元2009物理學報58 7809］

［6］Kulkarni A J，Zhou M 2006 Appl．Phys．Lett．88 141921

［7］Hegedus P J，Abramson A R 2006 Heat Mass Transfer．49 4921

［8］Liang X G，Sun L 2005 Microscale Thermophys．Eng．9 295

［9］Volz S G，Saulnier J B 2000 Microelectron．J．31 815

［10］Yang P，Liao N B 2008 Appl．Phys．A 92 329

［11］Yang P，Liao N B 2008 J．Thermophys Heat Transfer 22 581

［12］Zhong Z R，Wang X W，Xu J 2004 Numer．Heat Transfer B 45 429

［13］Schelling P K，Phillpot S R，Keblinski P K 2002 Phys．Rev．B 65 144306

［14］Sinnott S B，Dickey E C 2003 Mater．Sci．Eng．R 43 1

［15］Wunderlich W 1998 Phys．Stat．Sol．A 170 99

［16］Naicker P K，Cummings P T，Zhang H Z，Banfield J F 2005 J．Phys．Chem．B 109 15243

［17］Sergey V D，Nobuhiro Y，Yutaka K 2004 Acta Mater．52 1959

［18］White A 2000 Intermolecular Potentials of Mixed System:Testing the Lorentz-Berthelot Mixing Rules With Ab Initio calculations (Melbourne Victoria:DSTO Aeronautical and Maritime Research Laboratory)p1

［19］Stevens R J，Zhigilei L V，Norris P M 2007 Int．J．Heat Mass Transfer．50 3977

［20］Kim D J，Kim D S，Cho S 2004 Int．J．Thermophys．25 281

［21］Abramson A R，Tien C L，Majumdar A 2002 J．Heat Transfer．124 963

PACS:66.70.Df，83.10.Rs，68.35.－p，02.60.－x

*Project supported by the National Natural Science Foundation of China(Grant Nos．61076098，50875115)，the Natural Science Foundation of

Jiangsu Province，China(Grant No．BK2008227)，and the Graduate Innovative Program of Jiangsu Province，China(Grant No．CX10 B-252 Z)． E-mail:yangping1964@163.com

Molecular dynamics simulation of thermalconductivity for the TiO2/ZnO nano-film interface*

Yang PingWu Yong-Sheng Xu Hai-Feng Xu Xian-Xin Zhang Li-Qiang Li Pei
(School of Mechanical Engineering，Jiangsu University，Zhenjiang 212013，China)
(Received 5 July 2010;revised manuscript received 28 September 2010)

In the paper，the equilibrium molecular dynamics and Buckingham potential function are used to investigate the thermal conductivity of TiO2/ZnO nano-film interface along to［0001］(z-axis)．The effects of the equilibrium temperature，the thickness and the cross section of the nano-film interface on the thermal conductivity of TiO2/ZnO are investigated by optimizing the cut-off radius(rc)and the time step for initial condition of molecular dynamics．The results indicate that the thermal conductivity of TiO2/ZnO nano-film interface decreases with temperature increasing from 300 K to 600 K，and increases with film thickness increasing from 1.8 to 5 nm．Finally，the relationship between the thermal conductivity and the thickness of TiO2/ZnO nano-film interface is discussed．

thermal conductivity，molecular dynamics，TiO2/ZnO nano-film interface，numerical simulation

*國家自然科學基金(批準號:61076098，50875115)、江蘇省自然科學基金(批準號:BK2008227)和江蘇省研究生創新計劃(批準號: CX10 B-252 Z)資助的課題．

E-mail:yangping1964@163.com