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

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

用分子動(dòng)力學(xué)方法預(yù)測硅納米薄膜的熱導(dǎo)率

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

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

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

分子動(dòng)力學(xué)方法；熱導(dǎo)率；硅納米薄膜；尺寸效應(yīng)

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

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

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

1 非平衡分子動(dòng)力學(xué)方法

依據(jù)模擬過程中體系所處的熱力學(xué)狀態(tài)，分子動(dòng)力學(xué)方法分為兩大類別：第一種是平衡分子動(dòng)力學(xué)方法；第二種是非平衡分子動(dòng)力學(xué)方法.平衡分子動(dòng)力學(xué)方法主要根據(jù)線性響應(yīng)理論中Green-Kubo公式所建立的輸運(yùn)系數(shù)與平衡時(shí)間相關(guān)聯(lián)來獲得模擬體系的熱導(dǎo)率；非平衡分子動(dòng)力學(xué)方法主要根據(jù)體系內(nèi)的溫度分布和熱流密度情況由Fourier定律獲得熱導(dǎo)率.由于本文選擇了第二種方法，所以下面將重點(diǎn)介紹它.

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

圖1 非平衡輸運(yùn)模型

2 COMPASS力場

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

COMPASS力場是采用從頭計(jì)算方法并結(jié)合實(shí)驗(yàn)數(shù)據(jù)經(jīng)驗(yàn)參數(shù)化技巧而研發(fā)出來的.分子內(nèi)的鍵參數(shù)主要由從頭計(jì)算方法得到；范德華非鍵合參數(shù)通過分子動(dòng)力學(xué)模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)相結(jié)合方法而獲得，也就是：除了像一般的分子力場通常所要考慮的兩個(gè)實(shí)驗(yàn)數(shù)據(jù)(X光和振動(dòng)光譜)作為參照之外，COMPASS力場還采用對液態(tài)分子動(dòng)力學(xué)模擬得到的結(jié)合能、液態(tài)密度等同時(shí)作為參數(shù)化的實(shí)驗(yàn)標(biāo)準(zhǔn)，從而得到了更加優(yōu)化的力場性能.COMPASS力場的勢函數(shù)形式如下：

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

3 結(jié)果與討論

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

以厚度為8.688 nm的硅薄膜為例，經(jīng)300萬步計(jì)算后得到各層的溫度分布曲線如圖2所示.由圖2可以看到，各層的溫度分布在高低控溫層之間所形成的曲線是比較光滑的.淺色圓點(diǎn)分布區(qū)域即為控溫區(qū)，控溫區(qū)內(nèi)存在溫差，但這是合理的，因?yàn)橛袩崃鲝脑搮^(qū)域通過，因此必然也有溫差.

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

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

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

表1 硅納米薄膜的熱導(dǎo)率的計(jì)算結(jié)果

圖4 熱導(dǎo)率隨膜厚的變化

由氣體動(dòng)力學(xué)理論導(dǎo)出的介質(zhì)導(dǎo)熱系數(shù)表達(dá)式為[13]

(2)

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

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

4 結(jié) 論

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

[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].北京：科學(xué)出版社，2002：2-5.

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

[4] 吳國強(qiáng),孔憲仁,孫兆偉,等. 單晶硅薄膜法向熱導(dǎo)率的分子動(dòng)力學(xué)模擬[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào),2007,39(9):1366-1369.

[5] 華鈺超,董源,曹炳陽. 硅納米薄膜中聲子彈道擴(kuò)散導(dǎo)熱的蒙特卡羅模擬[J].物理學(xué)報(bào),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 基金項(xiàng)目：遼寧省自然科學(xué)基金資助項(xiàng)目(2015020079) 作者簡介：白素媛(1971-)，女，遼寧海城人，遼寧師范大學(xué)副教授，博士.Email:dl_bsy@163.com

1000-1735(2017)01-0030-05

10.11679/lsxblk2017010030

O484；TK124

A