陳 慧 魏志勇 陳偉宇 劉晨晗 畢可東 陳云飛

(1東南大學機械工程學院, 南京 211189)(2江蘇海事職業(yè)技術學院船舶與海洋工程學院, 南京 211170)

摻雜硅納米薄膜法向熱導率的分子動力學模擬

陳 慧1,2魏志勇1陳偉宇1劉晨晗1畢可東1陳云飛1

(1東南大學機械工程學院, 南京 211189)(2江蘇海事職業(yè)技術學院船舶與海洋工程學院, 南京 211170)

采用非平衡態(tài)分子動力學方法計算了平均溫度為300 K時膜厚2.2～104.4 nm的硅納米薄膜以及摻鍺的硅納米薄膜的法向熱導率.采用隨機摻雜方式在硅納米薄膜中摻入鍺原子,摻雜濃度分別為5%和50%.模擬結果表明,相同膜厚的摻鍺硅薄膜的法向熱導率比純硅薄膜的法向熱導率小很多,摻雜前后的硅薄膜法向熱導率均隨著膜厚的增大而增大.通過計算體態(tài)硅熱導率關于聲子平均自由程的累積函數(shù),發(fā)現(xiàn)對體態(tài)硅熱導率主要貢獻是聲子平均自由程為2～2 000 nm的聲子,而摻鍺的體態(tài)硅中對熱導率貢獻約占80%的聲子平均自由程為0.1～30 nm,遠小于純硅中對熱導率主要貢獻的聲子平均自由程.

硅納米薄膜;熱導率;摻雜;分子動力學;聲子平均自由程;熱導率累積函數(shù)

微電子器件的日益小型化使得其特征尺寸已縮小至納米量級.一方面,對高密度熱量的有效轉移已成為影響微電子器件工作性能和可靠性的關鍵因素,因而期望構成微電子器件的納米材料具有較高的熱導率;另一方面,提升熱電器件制冷效率的一個重要途徑是在不影響熱電材料電導率和Seebeck系數(shù)的情況下降低熱導率,將顯著提高熱電材料的品質因數(shù)[1-2],大幅提升熱電器件的性能.目前熱電設備中使用較多的硅材料,其熱電性能較差.因此,有效降低硅材料的熱導率以提高其品質因數(shù)一直是該領域的研究熱點之一.研究發(fā)現(xiàn),當納米器件或結構的特征尺寸降至與材料中傳遞熱量的聲子平均自由程(MFP)相當時,納米結構的熱傳導將出現(xiàn)顯著的尺寸效應[3].此時,材料的熱導率主要取決于結構尺寸.由于電子的MFP比聲子的MFP小得多,當材料尺度接近聲子的MFP時,聲子熱導率會降低,而電子的電導率則影響較小.因此,利用納米結構對聲子輸運的影響,降低材料的熱導率,可有效提高熱電材料的能量轉換效率.納米材料導熱出現(xiàn)尺寸效應的主要原因在于,當材料的特征尺寸遠小于聲子平均自由程時,基本不會發(fā)生聲子-聲子間散射,導熱處于彈道輸運階段,熱導率隨著尺寸的增加而增大.例如,Feng等[4]采用非平衡態(tài)分子動力學(NEMD)方法模擬了硅薄膜法向熱導率,發(fā)現(xiàn)在膜厚為2～10 nm范圍內,硅薄膜熱導率顯著低于體態(tài)實驗值,且隨著薄膜厚度的增大而增大.Gomes等[5]采用平衡態(tài)分子動力學(EMD)方法模擬了薄膜厚度為2～217 nm的硅薄膜面向和法向熱導率.蘇高輝等[6]采用蒙特卡洛方法研究了硅薄膜法向熱導率出現(xiàn)明顯尺度效應的臨界尺寸,當薄膜厚度等于1 μm 時,開始呈現(xiàn)明顯的尺度效應,此時擴散輸運和彈道輸運共存.Sun等[7]采用分子動力學方法研究了硅中聲子輸運的尺寸效應,認為硅薄膜熱導率下降的主要原因是由于聲子的限制和界面熱阻的影響.Yang等[8]采用NEMD方法研究了硅納米結構熱導率的尺寸效應,發(fā)現(xiàn)硅納米結構的長度對其熱導率具有顯著影響,并將其歸因于硅納米結構中聲子有效平均自由程的顯著減小以及納米結構尺寸的限制導致對熱傳導作貢獻的聲子數(shù)量的減少.此外,硅的晶格缺陷以及摻雜等因素也會增加對聲子的散射幾率而使聲子平均自由程減小,從而降低材料的熱導率.Asheghi等[9]采用實驗方法測量了溫度為15～300 K時分別摻有硼和磷的單晶硅薄膜的熱導率,結果表明雜質的存在顯著降低了材料的熱導率.Bi等[10]采用分子動力學模擬方法研究了在氬晶體中摻雜氪后對其熱導率的影響,結果表明氬晶體的熱導率降低了近1.9倍.Garg等[11]的研究表明硅和鍺之間的質量差引起的質量失配對于抑制硅鍺合金中的熱傳導起到了重要作用.目前關于摻鍺的硅納米薄膜的導熱性能研究較少,且對熱導率做貢獻的聲子平均自由程的大小尚未獲得一致的結論.

本文采用NEMD方法研究了摻鍺對硅納米薄膜熱導率的影響,討論了熱導率隨膜厚的變化關系,計算了熱導率關于聲子平均自由程的累積函數(shù)[12],得出了對熱導率做貢獻的聲子平均自由程的范圍.

1 分子動力學模型

在硅納米薄膜中以隨機方式摻入鍺(Ge)原子,原子的初始位置按照金剛石晶格結構排列在格點上,采用SW作用勢[13]描述硅原子之間、鍺原子之間以及硅原子與鍺原子之間的相互作用力.圖1為隨機摻鍺的硅納米薄膜法向熱導率的NEMD模型示意圖,模擬區(qū)域大小用Nxa×Nya×Nza表示,其中Nx,Ny和Nz分別為x,y和z方向上的晶胞(unit cell, UC)數(shù),a為硅的晶格常數(shù).x方向為熱流方向,分別在y和z方向上施加周期性邊界條件,以實現(xiàn)對無限大平面納米薄膜的模擬.模擬發(fā)現(xiàn),垂直于熱流方向的橫截面積大小對模擬結果影響很小,本文中模擬區(qū)域橫截面積均為10 UC×10 UC.在模擬區(qū)域的兩端分別設置厚度為2 UC的固定絕熱墻,熱、冷域的厚度也均為2 UC,模擬樣品的膜厚均不計算固定墻和熱、冷域的厚度.

圖1 摻雜硅納米薄膜法向熱導率的NEMD模型

模擬系統(tǒng)內原子速度的初始化遵循Maxwell分布,采用Velocity-Verlet算法對運動方程進行時間積分,積分步長Δt=0.5 fs,整個MD過程共模擬1.05×107步(5.25 ns)以上.最初的2.0×106步使模擬系統(tǒng)在NVT系綜下演進,通過調整原子的速度使模擬系統(tǒng)的溫度為恒定值.隨后,系統(tǒng)在NVE系綜下運行5.0×105步,保持系統(tǒng)能量守恒.此后至少再模擬8.0×106步,通過加減能量法在模擬區(qū)域內形成溫度梯度,即每次給熱域內增加能量Δε,同時從冷域內減少能量Δε.系統(tǒng)達到穩(wěn)定狀態(tài)時再運行6.0×106步以統(tǒng)計薄膜法向(x方向)的溫度梯度.將整個模擬區(qū)域沿x方向劃分為2Nx層,根據(jù)能量均分定理可計算每一層的區(qū)域溫度Tj為

(1)

式中,m為原子質量;Nj為第j層所含的原子數(shù);vi為原子i的速度;kB為Boltzmann常數(shù).每隔1 000步統(tǒng)計一次各模擬層的溫度,將最后運行的6.0×106步內統(tǒng)計的溫度取平均值作為各層在熱平衡后的平均溫度,溫度分布曲線如圖2所示.根據(jù)分布曲線可計算出溫度梯度▽T.模擬區(qū)域的法向熱流密度Jx為

(2)

式中,A為薄膜的橫截面積;Δt為積分步長.依據(jù)Fourier定律,摻鍺硅納米薄膜的法向熱導率k為

(3)

圖2 摻鍺硅納米薄膜的溫度分布曲線

2 模擬結果與分析

本文分別模擬了300 K時鍺原子的摻雜濃度為5%(Si0.95Ge0.05)、50%(Si0.5Ge0.5)以及不摻雜時的硅納米薄膜在不同膜厚(2.2～104.4 nm)時的法向熱導率,模擬結果如圖3所示.在相同膜厚時,摻鍺硅薄膜的法向熱導率均顯著低于純硅薄膜的法向熱導率.Yang等[14]采用NEMD方法對同位素摻雜的硅納米線的熱導率研究表明,隨摻雜濃度的增加,熱導率呈指數(shù)變化,極小的摻雜濃度即會引起熱導率顯著降低,而當熱導率下降到一個最低值后,會隨著摻雜濃度的增大而增大.

圖3 硅納米薄膜法向熱導率隨膜厚的變化(T=300 K)

硅薄膜中對熱傳導的貢獻主要來自于聲子的輸運,忽略自由電子的作用,薄膜內除了聲子與聲子間散射外,薄膜邊界以及晶格缺陷和雜質等對聲子的散射占主導地位.在摻鍺硅薄膜中,由于鍺原子和硅原子質量以及直徑的差異引發(fā)了聲子-雜質間散射,而且,隨機分布的大質量原子增加了硅薄膜中的晶格缺陷,使聲子群速率降低,聲子自由度減小,因此,摻鍺大大降低了硅薄膜的熱導率.

從圖3可看出,當溫度為300 K時,2種摻鍺濃度以及不摻鍺的硅納米薄膜的法向熱導率均隨著膜厚的增加而增大,呈現(xiàn)出明顯的尺寸效應.

如圖4所示,采用基于Matthiessen的倒數(shù)擬合方法[15]外推,有1/k=A/t+B,其中,k為熱導率，t為硅薄膜的法向厚度,當1/t=0時,1/B即為體態(tài)硅的熱導率kbulk,A/B為聲子平均自由程Λbulk,具體數(shù)據(jù)見表1.當溫度為300 K時,體態(tài)硅的熱導率為153.14 W/(m·K),聲子平均自由程為110.7 nm.摻鍺濃度為5%和50%的體態(tài)硅熱導率分別為5.84和1.62 W/(m·K).

關于硅的聲子平均自由程有多種不同的研究結果.根據(jù)聲子氣動理論計算出室溫下硅晶體的聲子平均自由程為41 nm[16];Alvarez等[17]根據(jù)傳統(tǒng)的Debye模型,計算出室溫時硅的聲子有效平均自由程大約為40～43 nm;而Dames等[18]研究認為,硅的聲子平均自由程為210 nm左右.為了明確不同平均自由程的聲子對熱導率的貢獻,最近研究者[12, 19-21]提出了熱導率關于聲子平均自由程的累積函數(shù)計算方法，并進行了深入研究.

圖4 硅納米薄膜法向熱導率倒數(shù)擬合(T=300 K)

物理量不摻鍺摻雜鍺濃度/%550kbulk/(W·(m·K)-1)153.145.841.62Λbulk/nm110.708.8010.50

(4)

式中,Λ為體態(tài)硅中三聲子散射平均自由程.根據(jù)文獻[19]可計算出

(5)

從而得到需求解的Fredholm積分方程為

(6)

式中,kfilm(t)為膜厚為t的硅薄膜法向熱導率,由NEMD模擬獲得;α(Λ)為熱導率關于聲子平均自由程的累積函數(shù),表示聲子平均自由程小于Λ的聲子對熱導率貢獻的百分數(shù).采用類似文獻[20]的方法,先給定一個初值,然后采用一種優(yōu)化方法求解方程(6)[12],得到

(7)

式中,Λ0可用表1中的Λbulk代替.圖5給出了300 K時不摻鍺、摻鍺濃度為5%和50%的體態(tài)硅的熱導率累積函數(shù)與聲子平均自由程的關系.

由圖5可看出,對于不摻鍺的體態(tài)硅,熱導率主要取決于聲子平均自由程為2～2 000 nm的聲子,其對熱導率的貢獻約占80%,聲子平均自由程大于40 nm的聲子對熱導率的貢獻超過50%,這與Regner等[21]采用測量方法獲得晶體硅的熱導率累積函數(shù)的研究結果較為一致.由圖5可知,聲子平均自由程大于100 nm的聲子對熱導率的貢獻仍約占40%.因此,當硅薄膜的模擬厚度小于100 nm時,大多數(shù)聲子處于彈道熱輸運機制,薄膜邊界對聲子的散射起主導作用,聲子的有效平均自由程與膜厚相關,導致硅薄膜法向熱導率的NEMD模擬值呈現(xiàn)明顯的尺寸效應.

圖5 熱導率對聲子平均自由程的累積函數(shù)

對于摻鍺的體態(tài)硅,摻鍺濃度為5%和50%的熱導率累積函數(shù)隨聲子平均自由程的變化曲線非常接近.圖5表明,摻鍺的體態(tài)硅中對熱導率的貢獻80%均來自于聲子平均自由程為0.1～30 nm的聲子.對摻鍺硅熱導率作主要貢獻的聲子平均自由程遠小于不摻鍺硅的聲子平均自由程,這是由于聲子-雜質間的散射引起的.

3 結論

1) 采用NEMD方法研究了300 K時硅納米薄膜以及摻鍺的硅納米薄膜的法向熱導率隨膜厚的變化關系，相同膜厚的摻鍺硅薄膜的法向熱導率比純硅薄膜的小很多.

2) 不論摻鍺與否，硅薄膜的法向熱導率均隨著膜厚的增大而增大，呈現(xiàn)出明顯的尺寸效應.

3) 采用倒數(shù)擬合方法得到體態(tài)硅的熱導率為153.14 W/(m·K)，摻鍺濃度為5%和50%的體態(tài)硅熱導率分別為5.84和1.62 W/(m·K).

4) 通過計算熱導率關于聲子平均自由程的累積函數(shù)，分析了對熱導率做主要貢獻的聲子平均自由程.當溫度為300 K時，對體態(tài)硅熱導率作主要貢獻來自于聲子平均自由程為2～2 000 nm的聲子，聲子平均自由程大于100 nm的聲子對熱導率的貢獻約占40%.這表明，傳統(tǒng)觀點嚴重低估了對體態(tài)硅熱導率貢獻的聲子平均自由程的大小.摻鍺濃度為5%和50%時，對熱導率貢獻約占80%的聲子平均自由程均為0.1～30 nm，遠小于純硅中對熱導率作主要貢獻的聲子平均自由程.

References)

[1]Liu W S, Yan X, Chen G, et al. Recent advances in thermoelectric nanocomposites[J].NanoEnergy, 2012, 1(1): 42-56. DOI:10.1016/j.nanoen.2011.10.001.

[2]Hochbaum A I, Chen R K, Delgado R D, et al. Enhanced thermoelectric performance of rough silicon nanowires[J].Nature, 2008, 451(7175): 163-167. DOI:10.1038/nature06381.

[3]Chen G.Nanoscaleenergytransportandconversion[M]. Oxford,UK: Oxford University Press, 2005: 283-291.

[4]Feng X L, Li Z X,Guo Z Y. Molecular dynamics simulation of thermal conductivity of nanoscale thin silicon films[J].MicroscaleThermophysicalEngineering, 2003, 7(2): 153-161. DOI:10.1080/10893950390203332.

[5]Gomes C J, Madrid M, Goicochea J V, et al. In-plane and out-of-plane thermal conductivity of silicon thin films predicted by molecular dynamics[J].JournalofHeatTransfer, 2006, 128(11): 1114-1121. DOI:10.1115/1.2352781.

[6]蘇高輝,楊自春,孫豐瑞. 硅薄膜導熱系數(shù)微尺度效應的臨界尺寸[J]. 納米技術與精密工程, 2014,12 (3): 176-181. DOI:10.13494/j.npe.20130086. Su Gaohui, Yang Zichun, Sun Fengrui. Critical size for microscale effect of silicon film thermal conductivity[J].NanotechnologyandPrecisionEngineering, 2014, 12(3): 176-181. DOI:10.13494/j.npe.20130086. (in Chinese)

[7]Sun L, Murthy J Y. Domain size effects in molecular dynamics simulation of phonon transport in silicon[J].AppliedPhysicsLetters, 2006, 89(17): 171919. DOI:10.1063/1.2364062.

[8]Yang Y W, Liu X J,Yang J P. Nonequilibrium molecular dynamics simulation for size effects on thermal conductivity of Si nanostructures[J].MolecularSimulation, 2008, 34(1): 51-55. DOI:10.1080/08927020701730419.

[9]Asheghi M, Kurabayashi K, Kasnavi R, et al. Thermal conduction in doped single-crystal silicon films[J].JournalofAppliedPhysics, 2002, 91(8): 5079-5088. DOI:10.1063/1.1458057.

[10]Bi K D, Zhao Y Y, Chen Y F, et al. The effects of different doping patterns on the lattice thermal conductivity of solid Ar[J].JournalofPhysicsandChemistryofSolids, 2012, 73(2): 204-208. DOI:10.1016/j.jpcs.2011.11.001.

[11]Garg J, Bonini N, Kozinsky B, et al. Role of disorder and anharmonicity in the thermal conductivity of silicon-germanium alloys: A first-principles study[J].PhysicalReviewLetters, 2011, 106(4): 045901. DOI:10.1103/PhysRevLett.106.045901.

[12]Wei Z, Yang J, Chen W, et al. Phonon mean free path of graphite along the c-axis[J].AppliedPhysicsLetters, 2014, 104(8): 081903. DOI:10.1063/1.4866416.

[13]Stillinger F H,Weber T A. Computer simulation of local order in condensed phases of silicon[J].PhysicalReviewB, 1985, 31(8): 5262-5271. DOI:10.1103/physrevb.31.5262.

[14]Yang N, Zhang G, Li B W. Ultralow thermal conductivity of isotope-doped silicon nanowires[J].NanoLetters, 2008, 8(1): 276-280. DOI:10.1021/nl0725998.

[15]Schelling P K, Phillpot S R,Keblinski P. Comparison of atomic-level simulation methods for computing thermal conductivity[J].PhysicalReviewB, 2002, 65(14): 144306. DOI:10.1103/physrevb.65.144306.

[16]Henry A S,Chen G. Spectral phonon transport properties of silicon based on molecular dynamics simulations and lattice dynamics[J].JournalofComputationalandTheoreticalNanoscience, 2008, 5(2): 141-152. DOI:10.1166/jctn.2008.2454.

[17]Alvarez F X,Jou D. Memory and nonlocal effects in heat transport: From diffusive to ballistic regimes[J].AppliedPhysicsLetters, 2007, 90(8): 083109. DOI:10.1063/1.2645110.

[18]Dames C, Chen G. Theoretical phonon thermal conductivity of Si/Ge superlattice nanowires[J].JournalofAppliedPhysics, 2004, 95(2): 682-693.

[19]Yang F,Dames C. Mean free path spectra as a tool to understand thermal conductivity in bulk and nanostructures[J].PhysicalReviewB, 2013, 87(3): 035437. DOI:10.1103/physrevb.87.035437.

[20]Minnich A J. Determining phonon mean free paths from observations of quasiballistic thermal transport[J].PhysicalReviewLetters, 2012, 109(20): 205901. DOI:10.1103/PhysRevLett.109.205901.

[21]Regner K T, Sellan D P, Su Z, et al. Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance[J].NatureCommunications, 2013, 4: 1640-01-1640-07. DOI:10.1038/ncomms2630.

[22]Majumdar A. Microscale heat-conduction in dielectric thin-films[J].JournalofHeatTransfer, 1993, 115(1): 7-16. DOI:10.1115/1.2910673.

Molecular dynamics simulation of out-of-plane thermal conductivity of doped silicon nanofilms

Chen Hui1, 2Wei Zhiyong1Chen Weiyu1Liu Chenhan1Bi Kedong1Chen Yunfei1

(1School of Mechanical Engineering, Southeast University, Nanjing 211189, China)(2Department of Ship and Ocean Engineering, Jiangsu Maritime Institute, Nanjing 211170, China)

The out-of-plane thermal conductivities of both pure silicon nanofilm and that doped with germanium (Ge) with film thicknesses ranging from 2.2 to 104.4 nm were calculated by non-equilibrium molecular dynamics (NEMD) simulation at an average temperature of 300 K. Ge atoms were doped in a random pattern with a doping density of 5% and 50%, respectively. The results show that there is a large reduction in thermal conductivities of the silicon nanofilm after doping with the same thickness. The out-of-plane thermal conductivity of the Ge-doping silicon nanofilm as well as that of a pure one increases with the increase of the film thickness. By calculating the thermal conductivity accumulation function in terms of phonon mean free path (MFP), it indicates that at 300 K the dominant contribution to the bulk thermal conductivity of pure silicon comes from phonons with MFPs between 2 and 2 000 nm, while the MFPs of phonons that contribute about 80% to the bulk thermal conductivity of Ge-doping silicon are between 0.1 and 30 nm, thus they are much less than those of the pure silicon.

silicon nanofilm; thermal conductivity; doping; molecular dynamics; phonon mean free path; thermal conductivity accumulation function

10.3969/j.issn.1001-0505.2017.03.013

2016-09-16. 作者簡介: 陳慧(1982—)，女，博士生；陳云飛(聯(lián)系人)，男，博士，教授，博士生導師，yunfeichen@seu.edu.cn.

國家自然科學基金資助項目(51205061, 51375092).

陳慧，魏志勇，陳偉宇，等.摻雜硅納米薄膜法向熱導率的分子動力學模擬[J].東南大學學報(自然科學版),2017,47(3):490-494.

10.3969/j.issn.1001-0505.2017.03.013.

TK124

A

1001-0505(2017)03-0490-05