硅鍺超晶格結構熱導率的分子動力學模擬

孫兆偉，張興麗

(哈爾濱工業大學衛星工程技術研究所，150001哈爾濱，sunzhaowei@sina.vip.com)

超晶格材料具有良好的熱電性質，在計算機芯片、MEMS、納米材料制造、航空航天等領域有廣泛的應用.由于在很多情況下，熱問題一直是這些技術發展的瓶頸問題，因此超晶格結構熱物性已成為熱傳導領域的1個新的研究熱點.分子動力學計算機模擬是研究微尺度材料熱傳輸特性的有力工具，特別是對于一些理論上難以說明或實驗中難以觀察的現象給出基于物質微觀結構及分子動力學關系的解釋.從統計物理的角度可以將分子動力學分為平衡分子動力學模擬(EMD)和非平衡分子動力學模擬(NEMD).前者計算平衡系統的熱流與時間的相關函數，然后通過Green-Kubo關系式得到熱導率;后者需要對系統施加能量產生熱流，得到系統的溫度梯度，根據Fourier定律計算熱導率.近年來，許多學者通過分子動力學方法對超晶格材料的熱導率進行了深入的研究.Chen等［1］利用NEMD模擬方法研究了晶格周期長度對超晶格結構熱導率的影響，得出當晶格的周期長度小于聲子的平均自由程時，熱導率會出現最小值的結果;Deyu等［2］對Si/SiGe納米線結構進行研究，表明其熱導率受邊界散射機制的影響隨著溫度的升高而有下降的趨勢.這些研究都表明微尺度下超晶格結構的熱導率明顯小于大體積情況下的熱導率值.本文在此研究的基礎上利用NEMD模擬方法研究了Si/Ge超晶格結構熱導率隨周期長度、周期數及溫度的變化規律，并且討論了界面熱阻對超晶格結構熱導率的影響特性.

1 分子動力學模擬方法

1.1 超晶格結構導熱模型

超晶格結構導熱模型如圖1所示.在Y、Z方向施加周期性邊界條件時，由于垂直熱流方向的橫截面積過小會對熱導率的計算結果產生誤差［3-4］，因此本文選取的YOZ橫截面積為4 UC× 4 UC(Unit Cell，晶胞長度).為建立熱流方向的溫度梯度，在X方向上布置隨機恒溫熱墻.模型中的高低溫熱墻材料同各自臨近的超晶格材料相同，并設定其厚度為3 UC;模型最外層設置厚度為2 UC的絕熱壁，作用是減少導熱層內的粒子蒸發，防止與外界產生熱量交換，并設定該區域粒子的速度為0.

圖1 硅鍺超晶格結構非平衡分子動力學模擬模型

1.2 作用勢能模型

Stillinger和Weber在1985年提出了穩定金剛石結構的多體勢能函數來描述超晶格結構分子之間的相互作用，N個粒子的總勢能為［5］

上式中兩體、三體作用式分別用下列公式表達:

其中θjik是rij和rik的夾角，式中的各物理參數值如表1所示［6］.

表1 硅鍺S－W勢能參數

計算分子間作用勢能和作用力需要花費大量的運算時間，因此本文采用Cell列表法為Velet列表法建立列表，用Velet列表法計算原子間勢能和作用力，這樣可以大大提高計算速度［7］.

1.3 局域溫度修正

由Boltzmann分布函數可以獲得導熱區域內每個粒子的動能，系統溫度可以表示為

式中，kB是Boltzmann常數，N為系統中的粒子數.

式(1)成立的前提條件之一是系統溫度遠高于材料的debye溫度θD，由于模擬過程中的平衡溫度低于Si的Debye溫度(645 K)，因此需對系統的局域穩定進行量子化修正，才能獲得超晶格結構熱導率的真實值［8］.假設系統總能量是動能的2倍，通過求解積分方程得到系統的真實溫度，如下所示:

式的右邊是系統中粒子的總能量，D(ω)為聲子密度分布函數;ω為聲子頻率;n為對應于熱平衡溫度T的聲子平均占有數，該占有數滿足Planck分布(即Bose-Einstein統計).

這樣Fourier定律中的溫度梯度可以修正為

1.4 熱導率計算

利用非平衡分子動力學(NEMD)方法，對系統施加1個給定的熱流，在系統內部自然形成1個溫度梯度.系統的熱流通過動能的變化量反應，X方向的熱流密度可以表示為

由Fourier定律得到熱導率為

其中:ΔKinH、ΔKinc分別是高溫端和低溫端在時間Δt內的動能變化;Ly·Lz為截面YOZ的面積;Δt為時間步長，Δt=0.005τ0(τ0為特征時間，且

2 模擬結果

2.1 導熱區域的溫度分布

本次對Si/Ge超晶格結構的NEMD模擬總的模擬步長數為7×107，其中前106步使系統趨于平衡.在不同的溫度下，高低溫熱墻的溫度設為Thot=(T+20)K及Tcold=(T-20)K.根據式(2)得到周期長度為10 UC，周期數為2的超晶格結構X方向導熱層的溫度分布，如圖2所示.

圖2 界面溫度的變化

從圖中可以看出在Si-Ge每個交界面處都存在溫度跳躍，在溫度較高的第1個界面處的突變要比溫度較低處更加明顯.超晶格結構界面處平衡溫度的主要影響因素是:1)高溫熱墻與其臨近的Si粒子層和低溫熱墻與其臨近的Ge粒子層之間的界面特性;2)高低溫熱墻的厚度.由于隨著界面周期數的增加，高低溫熱墻厚度對導熱層溫度的影響逐漸減弱［9］，因此本文中的溫度跳躍變化主要是受界面熱阻機制的影響，即靠近高溫熱墻處粒子的界面反射率增加使界面熱阻增加，從而引起溫度的下降［10］.另外Si和Ge粒子通過第1個界面的聲子振動波譜是連續，而在剩余的界面聲子振動波譜的連續性降低，這也導致第1個界面溫度下降得比其余界面更加明顯［11］.

2.2 周期長度對超晶格結構熱導率的影響

通過改變超晶格結構的周期長度來研究周期長度與熱導率之間的變化關系，如圖3所示.隨著周期長度的增加，在系統溫度為300 K和500 K時熱導率也隨之增大，并且兩者接近于線性關系.

超晶格結構的熱阻由兩層材料的固有熱阻和界面熱阻組成，因此本文結構的熱阻表達式為

式中的RBD為界面熱阻;RSi、RGe為硅鍺材料每層的熱阻，可近似為常數;n為界面層數，其值為n= 2N-1，N為超晶格的周期數.Chen等［10］通過理論公式推導得出固體超晶格結構有效熱導率表達式為

式中LP為超晶格結構的周期長度.上式說明超晶格結構的熱導率隨著周期長度的增大而增大，并且兩者呈線性關系，與本次的模擬結果吻合較好.

圖3 熱導率同周期長度的關系

2.3 周期數對超晶格結構熱導率的影響

超晶格結構中，周期數的變化也會使界面熱效應發生改變，從而導致熱導率發生變化.圖4是溫度為400 K時，周期長度分別為5 UC、10 UC的結構熱導率隨周期數的變化關系.從圖中可以很明顯看出2種結構的熱導率隨周期數的增加而增大，并且周期長度越大的結構，其周期數對熱導率的影響越明顯.這是因為周期數越大，Si/Ge交界面就越多，界面總熱阻隨之增加;同時周期數越大，也使導熱層總厚度增大，因此界面熱阻在總熱阻中所占的比例不斷下降，界面效應慢慢減弱，熱導率就逐漸增大.對于周期長度更大的結構，整個導熱層厚度的增加使溫度梯度下降更加劇烈，熱導率上升更快.

2.4 溫度對超晶格結構熱導率的影響

圖5是周期長度為10 UC的Si/Ge超晶格結構溫度與熱導率的變化關系，從圖中可以看出熱導率隨著溫度的升高而略微增大.由于溫度越高引起的熱流波動越大，使振動頻率大的聲子比振動頻率低的聲子要多，從而增加了粒子在界面處的非彈性散射幾率，界面熱阻變小，熱導率也隨之逐漸增大.與文獻［12］中的相應的SiGe合金結構相比，超晶格結構可以有效的減小材料熱導率，從而提高材料的熱電品質指數.

圖4 熱導率同周期數的關系

圖5 熱導率同溫度的變化關系

3 結論

利用非平衡分子動力學方法分析了硅鍺超晶格結構的熱傳導性能，仿真結果表明:

1)溫度在靠近高溫熱墻的第1個交界面發生明顯突變，這種現象應歸結于界面幾何特性的不同，導致聲子反射率不同.第1個交界面的反射率增加，界面熱阻因此增加，使溫度發生跳躍.

2)超晶格結構的法向熱導率隨著周期長度及周期數的增加都呈近似線性地增加.周期長度越大的結構，周期數對熱導率產生的影響越大.界面熱阻效應隨著周期數的增加而慢慢減弱.

3)隨著溫度的升高Si/Ge超晶格結構的熱導率也隨之略微增大，這與高溫下分子的界面散射機制有關，同時模擬結果明顯小于相應的合金結構熱導率.

［1］CHEN Yunfei，DEYU L，JENNIFER R.Minimum superlattice thermal conductivity from molecular dynamics［J］.Physical Review B，2005，72(174302):1-6.

［2］DEYU L，YIYING W，ARUN M.Thermal conductivity of Si/SiGe superlattice nanowires［J］.Applied Physics Letters，2003，15(83):3186-3188.

［3］孔憲仁，吳國強，孫兆偉，等.單晶碳和鍺薄膜熱導率的分子動力學模擬［J］.哈爾濱工業大學學報(自然科學版)，2006，38(4):517-519.

［4］SCHELLING P K，PHILLPOT S R，KEBLINSKI P. Comparison of atomic-level simulation methods for computing thermal conductivity［J］.Physical Review B，2002，65(14):144306.

［5］STILLINGER F，WEBER T.Computer simulation of local order in condensed phases of silicon［J］.Physical Review B，1985，31(8):5262-5271.

［6］SRINIVASAN S，MILLER R S.On parallel nonequilibrium molecular dynamics simulations of heat conduction in heterogeneous materials with three-body potentials: Si/Ge superlattice［J］.Numerical Heat Transfer B，2007，52(4):297-321.

［7］弗蘭克.分子模擬——從算法到應用［M］.汪文川，譯.北京:化學工業出版社，2002:351-354.

［8］吳國強，孔憲仁，孫兆偉，等.單晶硅薄膜法向熱導率的分子動力學模擬［J］.哈爾濱工業大學學報(自然科學版)，2007，39(9):1366-1369.

［9］LANDRY E S，MCGAUHEY A，HUSSEIN M.Molecular dynamics prediction of the thermal conductivity of Si/ Ge superlattices［C］//Proceedings of the HT2007 ASME-JSME Thermal Engineering Summer Heat Transfer.Vancouver:［s.n.］，2007:664-671.

［10］CHEN Yunfei，LI D，YANG J.Molecular dynamics simulation of Ar/Kr superlattice nanowires［J］.Physical Review B，2004，349(1/2/3/4):270-280.

［11］SAMVEDI V，TOMAR V.The role of interface thermal boundary resistance in the overall thermal conductivity of Si/Ge multilayered structures［J］.Nanotechnology. 2009，20(36):365701.

［12］ASHTON S，PATRICK K.Thermal resistivity of Si-Ge alloys by molecular dynamic simulation［J］.Journal of Applied Physics，2008，103(11):1-6.