用分子動力學方法預測硅納米薄膜的熱導率

白素媛， 李晶晶， 王立凡， 支明蕾， 王 斌

(遼寧師范大學 物理與電子技術學院，遼寧 大連 116029)

用分子動力學方法預測硅納米薄膜的熱導率

白素媛， 李晶晶， 王立凡， 支明蕾， 王 斌

(遼寧師范大學 物理與電子技術學院，遼寧 大連 116029)

采用非平衡分子動力學方法，基于優化的集成勢函數COMPASS力場，預測了室溫下(300 K)硅納米薄膜的熱導率.模擬結果表明：厚度約4～10 nm的硅薄膜的熱導率在3.06～7.28 W/(m·K)范圍，并且隨著膜厚的增加而增大，表現出明顯的尺度效應.在所計算的薄膜厚度范圍內，硅納米薄膜的熱導率與薄膜厚度呈現近似線性變化的關系.應用氣動理論對產生的尺度效應進行了初步的理論分析，當薄膜厚度在幾納米到十幾納米時，有效聲子平均自由程與膜厚有關，不再等于體材料的平均自由程.同時也將本文的預測結果與其他研究者采用Stillinger-Weber勢所進行的模擬結果進行了比較.為分子動力學方法在低維材料熱物性方面的研究提供了有益的參考.

分子動力學方法；熱導率；硅納米薄膜；尺寸效應

薄膜作為二維低維材料廣泛應用于微電子器件及微機電系統，伴隨著器件特征尺寸的不斷減小，薄膜厚度已經進入到納米量級.薄膜熱物性通常不同于體材料.熱物性，特別是熱導率直接影響器件和系統的穩定性和可靠性.由于現階段商用的熱物性測試儀器對于微納米薄膜熱導率測試還存在非常大的困難，采用分子動力學方法對納米尺度薄膜的熱導率進行預測成為行之有效的方法.該方法按照分子系統的時間演化來進行,并由此產生被研究體系內部相互作用分子的詳細運動軌跡，進而可以求出平衡或偏離平衡情況下分子體系的動力學性質.在分子動力學方法中，研究的分子體系是確定的，分子體系又嚴格遵守物理定律.由于該方法不需要對幾何對稱、輸運性行為、熱力學行為等作先驗假設，熱力學及輸運性質直接就是計算結果，因此亦被看作“計算機實驗”[1-2].

硅是微電子器件和微機電系統中最重要的薄膜材料之一，因此用分子動力學方法對硅納米薄膜熱導率的研究成為當前一個研究熱點[3-6].但由于人們選用的方法存在不同、采用的力場也不相同,研究系統的尺寸以及溫度范圍不同，這些都導致當前得到的硅納米薄膜熱導率的相關結果存在差異性及數據的零散性，還需要進一步驗證和明確.

針對厚度范圍約在4～10 nm的硅薄膜在室溫下(300 K)的熱導率進行了分子動力學模擬計算，采用的是非平衡動力學方法,在作用勢的選擇上，采用了優化的集成的勢函數COMPASS力場[7].

1 非平衡分子動力學方法

依據模擬過程中體系所處的熱力學狀態，分子動力學方法分為兩大類別：第一種是平衡分子動力學方法；第二種是非平衡分子動力學方法.平衡分子動力學方法主要根據線性響應理論中Green-Kubo公式所建立的輸運系數與平衡時間相關聯來獲得模擬體系的熱導率；非平衡分子動力學方法主要根據體系內的溫度分布和熱流密度情況由Fourier定律獲得熱導率.由于本文選擇了第二種方法，所以下面將重點介紹它.

在非平衡分子動力學方法中，需要建立起一個非平衡輸運模型.這個非平衡輸運模型通過對系統施加擾動來形成，比較常用的擾動方式是加入溫度梯度或者加入熱流.圖1即為本文模擬體系所建立起來的模型.模擬體系用三維空間坐標X、Y、Z表示，3個坐標方向上都取周期性邊界條件，在待考查的薄膜厚度Z方向上建立起穩態的非平衡熱輸運過程.實現的具體方法是：沿Z方向將體系均勻地分成多個層，使每層中所包含的原子數目相等或接近相等，將處于Z坐標1/4和3/4處的層分別作高溫控制層和低溫控制層(控溫層可以是1個或多個).當然，由于采用了周期性邊界條件，控溫層不是必須選在前面所述的位置上的，只需要選擇滿足2個控溫層的中心處相距等于模擬體系厚度的一半就可以.本文所提到的薄膜厚度即為模擬體系厚度的一半.

圖1 非平衡輸運模型

2 COMPASS力場

正確選擇被模擬材料系統的相互作用勢，是使用分子動力學方法進行準確模擬的前提.相互作用勢一般由勢函數來描述，它實際上描述的是電子云重疊的量子力學作用結果.由于材料的復雜性，大多數原子的排列的精確勢函數目前仍不明確，因此經驗或半經驗的勢函數常被采用，比如應用最為廣泛的對勢L-J勢[8]、包括兩體項及三體項的Stillinger-Weber勢[9]和Tersoff勢[10]等.另外還有集成的勢函數文件，被叫作力場，比如UNIVERSAL[11]、DERIDING[12]、COMPASS力場[7]等.由于COMPASS是新近發展起來的適合凝聚態應用的一個全新力場，具有更多的優點和先進性，因此本文選擇了該力場.

COMPASS力場是采用從頭計算方法并結合實驗數據經驗參數化技巧而研發出來的.分子內的鍵參數主要由從頭計算方法得到；范德華非鍵合參數通過分子動力學模擬結果和實驗數據相結合方法而獲得，也就是：除了像一般的分子力場通常所要考慮的兩個實驗數據(X光和振動光譜)作為參照之外，COMPASS力場還采用對液態分子動力學模擬得到的結合能、液態密度等同時作為參數化的實驗標準，從而得到了更加優化的力場性能.COMPASS力場的勢函數形式如下：

函數展開式由13個統計求和項構成，其內容可劃分為兩種項目：鍵合項和非鍵合項.鍵合項由展開式的前11個統計求和項構成，包括對角和非對角的交叉耦合項，主要由鍵長(b)、鍵角(θ)、二面角(φ)、面外角(χ)及含有上述參數兩個以上的耦合項構成.非鍵合項對應展開式中最后兩個統計求和項，包括范德華能和庫侖能.非鍵合項中，q為電荷量，rij為i原子與j原子間的距離.函數展開式中其他參數的意義為：k、H、V、A、B為力常數，ε為能量參數.

3 結果與討論

采用非平衡分子動力學方法模擬計算了厚度分別為4.344、5.430、6.516、7.602、8.688、9.744 nm(分別對應在Z方向上取16、20、24、28、32、36個基本單元的超級晶胞)的硅薄膜的熱導率，系統溫度為300 K，時間步1.5 fs，計算步數300萬步.

以厚度為8.688 nm的硅薄膜為例，經300萬步計算后得到各層的溫度分布曲線如圖2所示.由圖2可以看到，各層的溫度分布在高低控溫層之間所形成的曲線是比較光滑的.淺色圓點分布區域即為控溫區，控溫區內存在溫差，但這是合理的，因為有熱流從該區域通過，因此必然也有溫差.

在非平衡動力學方法中,穩態的非平衡過程的建立十分關鍵，它也可以通過檢查得到的薄膜熱導率隨計算步數的變化而趨于不變得以確認.仍以8.688 nm的硅薄膜為例，圖3所示的是它的熱導率隨計算總步數的變化情況.可以看到，80萬步(橫坐標刻度乘以200×103等于總計算步數)之前，熱導率波動較大，而到了約160萬步之后，熱導率變化已經非常微小，因此300萬步的計算時間是滿足建立穩態非平衡過程的.

圖2 厚度為8.688 nm的硅薄膜中各層的溫度分布情況

溫度300 K下硅納米薄膜的熱導率的模擬結果列在表1中.厚度范圍在4.344～9.744 nm 的硅薄膜的熱導率在3.06～7.28 W/(m·K)范圍.熱導率隨膜厚的變化關系如圖4所示.由圖4可以看到，在所計算的厚度范圍內，硅薄膜的熱導率隨著薄膜厚度的增加而增加，接近于線性變化，表現出明顯的尺度效應，且遠低于體材料值.

表1 硅納米薄膜的熱導率的計算結果

圖4 熱導率隨膜厚的變化

由氣體動力學理論導出的介質導熱系數表達式為[13]

(2)

其中，c是體積比熱容，v是聲子的平均速度，leff是材料的有效聲子平均自由程.當薄膜的厚度d遠大于體材料的聲子平均自由程l∞時，leff=l∞；當薄膜厚度d只有幾個納米到十幾個納米時，leff遠小于聲子平均自由程l∞，此時，有效聲子平均自由程leff由膜厚d所決定，因此，熱導率與膜厚呈現出近似的線性關系.在聲子氣動理論中，影響晶格熱導率的決定因素是聲子的有效平均自由程leff，而影響有效聲子平均自由程的因素很多，如邊界對聲子的散射、缺陷和雜質對聲子的散射、聲子與聲子間的散射等.

文獻[6]中利用Stillinger-Weber勢計算了300K下厚度在2.715～54.30nm的硅薄膜的熱導率，他們在低于20nm的厚度范圍內，觀察到了熱導率與膜厚呈近似線性關系，20nm以上厚度的薄膜的熱導率仍隨厚度增加而增加，但增加趨勢變緩.在與本文研究相近的厚度范圍內，他們得到3～10nm的硅薄膜的熱導率約在3.5～12W/(m·K),本文采用COMPASS力場得到的結果與之相接近,并且同樣得到了熱導率與膜厚之間近似線性變化的關系，可見本文的預測是合理的和可信的，至于還存在差異性，是由于選擇不同力場的原因.后繼工作將對更大厚度范圍內硅薄膜的熱導率開展進一步的研究.

4 結 論

采用COMPASS力場，運用非平衡分子動力學方法計算了室溫下硅納米薄膜的熱導率.穩態的非平衡過程的建立十分關鍵，通過選擇合適的時間步和計算總步數，通過觀察模擬系統各層溫度分布情況以及熱導率隨計算時間趨于穩定來確保非平衡過程進入穩態.模擬結果表明，厚度為4.344～9.774nm范圍的硅薄膜的熱導率在3.06～7.28W/(m·K)范圍，并且隨著膜厚的增加而增大，呈近似線性變化的關系.對于觀察到的尺度效應，用氣動理論進行了初步的理論分析.研究結果與采用Stillinger-Webber勢研究硅薄膜熱導率的文獻中的結果基本一致，為分子動力學方法在低維材料熱物性研究的應用提供了有益的參考.

[1] KOTAKE S,WAKURI S. Molecular dynamics study of heat conduction in solid materials[J].Fluids and Thermal Engineering,1994,37(1):103-108．

[2] 楊小震.分子模擬與高分子材料[M].北京：科學出版社，2002：2-5.

[3] 肖鵬，馮曉利，李志信．單晶硅薄膜法向熱導率分子動力學研究[J]．工程熱物理學報，2002，23(6)：195-199.

[4] 吳國強,孔憲仁,孫兆偉,等. 單晶硅薄膜法向熱導率的分子動力學模擬[J].哈爾濱工業大學學報,2007,39(9):1366-1369.

[5] 華鈺超,董源,曹炳陽. 硅納米薄膜中聲子彈道擴散導熱的蒙特卡羅模擬[J].物理學報,2013(24):178-186.

[6] WANG Z H，LI Z X. Research on the out-of-plane thermal conductivity of nanometer silicon film[J].Thin Solid Films,2006,515:2203-2206.

[7] SUN H. COMPASS:An ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds[J].J Phys Chem B,1998,102(38):7338-7364.

[8] LENNARD-JONES J E. On the determination of molecular fields[J]．Proc R Soc A,1924，106(738):463-477．

[9] STILLINGER F H,WEBER T A. Computer simulation of local order in condensed phases of silicon[J].Physical Review B,1985,31(8):5262-5271.

[10] TERSOFF J. Empirical interatomic potential for silicon with improved elastic properties[J].Physical Review B,1988,38(14):9902-9905.

[11] DAUBER-OSGUTHORPE P,ROBERTERS V A,OSGUTHORPE D J, et al. Structure and energetics of ligand binding to proteins:Escherichia coli dihydrofolate reductase-trimethoprim,a drug-receptor system[J].Proteins:Structure,Function,and Genetics,1988,4(1):31-47.

[12] MAYO S L,OLAFSON B D，GODDARD W A. DREIDING:A generic force field for molecular simulations[J].J Phys Chem,1990,94(26):8897-8909.

[13] TAMMA K K，ZHOU X. Macroscale and microscale thermal transport and thermo-mechanical interactions-some noteworthy perspectives[J].Journal of Thermal Stress,1998,21(3/4):405-449.

Prediction of the thermal conductivity of nanometer silicon films by molecular dynamics simulations

BAISuyuan,LIJingjing,WANGLifan,ZHIMinglei,WANGBin

(School of Physics and Electronic Technology， Liaoning Normal University， Dalian 116029， China)

Thermal conductivities of a serial of nanometer silicon films with thickness of 4～10 nm are calculated by no-equilibrium molecular dynamics method and an ab-initio force field, COMPASS. At room temperature (300 K), the thermal conductivities are 3.06～7.28 W/(m·K). The calculated results show that the thermal conductivity of nanometer silicon film decreases with the decrement of film thickness and is significantly lower than the bulk value. A notable size effect is observed due to the approximate linear relationship between the thermal conductivity and the film thickness. Phonon aerodynamic theory is employed to explain the size effect. The predicted results agree well with the results of other researcher’s obtained with Stillinger-Weber potential model. This work will be helpful in applying molecular dynamic simulations to study the thermal properties of low dimension materials.

molecular dynamics simulation;thermal conductivity;nanometer silicon film;size effect

2016-11-20 基金項目：遼寧省自然科學基金資助項目(2015020079) 作者簡介：白素媛(1971-)，女，遼寧海城人，遼寧師范大學副教授，博士.Email:dl_bsy@163.com

1000-1735(2017)01-0030-05

10.11679/lsxblk2017010030

O484；TK124

A