基于分子動力學對超晶格結構界面熱阻的模擬研究

張興麗， 孫兆偉

(哈爾濱工業大學衛星技術工程研究所，哈爾濱 150001)

超晶格結構材料是兩種或兩種以上材料按照周期性結構排列而形成的復合材料，具有良好的熱電性質，能顯著提高熱電轉換設備效率，在計算機芯片、MEMS器件、航空航天等領域有廣泛應用。對超晶格材料進行傳熱分析是當今的熱點問題，因為超晶格材料的尺寸達到了微/納米量級，經典熱傳導理論已經不能正確解析體系內非常規傳熱特性［1］。傳熱學及傳熱分析方法正面臨著從宏觀向微觀理論和方法過渡，許多理論及研究方法急需從更高層次和深度來觀察與解決。分子動力學模擬方法就是探求微/納尺度條件下的熱現象的規律和內在機制的有效方法，它可以從微觀細節著手，研究熱載流子(如聲子和電子)的行為，并依據統計力學原理得到系統的宏觀性質。從統計物理的角度可以將分子動力學分為平衡分子動力學模擬(EMD)和非平衡分子動力學模擬(NEMD)兩種模擬方法。前者是計算平衡系統的熱流與時間的相關函數，然后通過Green-Kubo關系式得到熱導率;后者需要對系統施加能量產生熱流，得到系統的溫度梯度，根據Fourier定律計算熱導率［2～4］。

界面熱阻是熱載流子在兩接觸固體界面層相互作用的結果。在微觀傳熱領域內，不同材料間界面熱阻的計算是當前研究的熱點問題，因為界面熱阻直接影響到材料的熱傳導性能，從而對微/納米器件的設計和熱優化產生影響。利用分子動力學方法對界面熱阻進行計算機模擬，可以從不同角度對界面熱阻的內在機理進行探討。分子動力學(MD)模擬通過求解牛頓運動方程得到每個粒子空間位置和運動狀態隨時間的演進狀況，在計算界面熱阻時，它不需要考慮每個粒子本身的散射特性，而只需根據勢能函數確定粒子之間的相互作用規律［5～8］。

本工作以Si/Ge超晶格結構為例，利用非平衡態分子動力學模擬方法從微觀機制出發研究超晶格結構界面熱阻的一些變化趨勢。

1 分子動力學模型

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

圖1 硅鍺超晶格結構非平衡分子動力學模擬模型Fig.1 The NEMD simulation model of Si/Ge superlattice structure

利用非平衡分子動力學方法模擬了平均溫度為400K時不同周期長度的硅鍺超晶格薄膜熱導率。模擬中采用Stillinger-Webber多體勢能函數來描述硅、鍺分子之間的相互作用［10，11］;采用 Verlet推導的Leap-frog算法進行粒子運動方程的數值積分。由于模擬過程中的平衡溫度低于Si的Debye溫度(645K)，因此需對系統的局域穩定進行量子化修正，才能獲得超晶格結構界面熱阻的真實值。根據經典Boltzmann統計可以獲得第j層中局域溫度為:

式中，[]...表示在總的模擬時間內的統計平均;kB是Boltzmann常數;Nj為第j層的粒子數。依據經典統計下的能量均分定理，第j層中的能量等式可以表示為:

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

在Debye近似下，能量等式(2)可以轉化為:

通過數值求解式(3)，得到與分子動力學模擬(MD)的局域溫度Tj，MD相對應的真實晶格溫度Tj。

當高溫熱墻或低溫熱墻的粒子與其他粒子作用時，溫度會發生改變。為使在導熱層區域形成穩定的溫度梯度，需要通過改變粒子的速率來增加高溫熱墻一定數量的動能，同時從低溫熱墻移出同等數量的動能。高、低溫熱墻能量的變化量可以表示為:

通過溫度梯度方向(X方向)的熱流計算公式為:

式中A為熱流方向的橫截面積;τ為模擬的時間步長。

超晶格結構兩種材料界面的示意圖如圖2所示。假設在連接區域的原子散射都是彈性的，這樣可忽略連接區域長度的影響，界面熱阻近似等于連接區域的熱阻［12］，界面熱阻可定義為:

式中TR，TL分別為界面左右兩端的溫度。

圖2 硅鍺超晶格結構界面示意圖Fig.2 Schematic diagram of the Si/Ge superlattice interface

3 結果和討論

非平衡分子動力學模擬在微正則(NVE)系統條件下進行，模擬的時間步長為1fs，總的模擬步長數為5×106，其中前2×106步使系統平衡。在不同的模擬溫度下，高低溫熱墻的溫度設為Thot=T+20 K及Tcold=T－20 K以形成溫度梯度。

圖3為周期長度為10UC的硅鍺超晶格薄膜在平均溫度為500K時從高溫熱墻到低溫熱墻導熱層溫度示意圖，從圖中可以明顯看出受界面熱阻的影響在硅鍺界面處溫度存在明顯的跳躍。三個硅鍺交界面的溫度變化值ΔT分別約為19K，10K，8K，由式(6)可知靠近高溫熱浴的第一個界面熱阻要遠遠大于其余界面，這個模擬結果與Abramson等［13］的研究成果相一致。他指出界面熱阻與界面數之間變化的非線性行為，是由于不同界面間傳遞熱流的聲子振動頻率不同致使熱流穿過每個界面產生的熱阻也不同。因此，最靠近高溫熱浴的界面熱阻比其余界面大，這也與傳遞熱流的聲子類型不同有關。

圖3 導熱層長度方向的溫度分布圖Fig.3 The temperature profile along the length of simulation cell for Si/Ge superlattice system

圖4是溫度為400K和500K時利用分子動力學模擬計算的最靠近高溫熱浴的Si/Ge界面熱阻隨周期長度的變化。由圖4可知隨著周期長度的增大，該界面熱阻逐漸減小，因此界面熱阻在總熱阻中所占的比例逐漸減小，界面效應減弱，Si/Ge超晶格薄膜的熱傳導性能不斷提高(如圖5所示)。

圖4 最靠近高溫熱浴的Si/Ge界面熱阻與周期長度的變化關系Fig.4 Thermal boundary resistance as a function period length for Si/Ge interface nearest to the hot bath nearest to the hot bath

平均溫度在300K到800K之間時，周期長度為10UC的Si/Ge超晶格薄膜界面熱阻模擬結果如圖6所示。從圖中可以明顯看出不同溫度下薄膜的界面熱阻也是不同的，隨著溫度的逐漸升高，界面熱阻越來越小。圖中比較了分子動力學模擬、散射失配理論(Diffuse Mismateh Model，DMM)研究［14，15］以及實驗數據［16］三種不同方法的計算結果，可以看出他們符合得比較好，變化趨勢大致相同。因為實驗過程中試樣的誤差以及DMM理論本身的適用條件限制，二者得到的界面熱阻結果比本工作利用分子動力學模擬結果略大。界面熱阻隨溫度升高而減小的原因主要是因為溫度的升高使界面處發生非彈性散射的概率增加。在高溫下界面處的高頻聲子只能發生非彈性散射，轉化為低頻聲子才可以傳遞能量。因此聲子的非彈性散射增加了界面的熱傳導能力，界面熱阻也隨之下降［7］。

4 結論

(1)受界面熱阻機制的影響，在超晶格導熱區域的界面處會發生溫度的突變，并且在最靠近高溫熱浴的界面溫度突變最為明顯。因此，最靠近高溫熱浴的界面熱阻對整個結構的熱傳導能力起著決定性作用。

(2)界面熱阻會隨超晶格結構周期長度的增大而逐漸減小，超晶格結構的導熱能力會隨之相應提高。

(3)隨著溫度的逐漸升高，界面熱阻會越來越小。這與聲子在高溫下發生非彈性散射概率的增加有關。

［1］吳勇華，楊決寬，陳云飛，等.超晶格薄膜熱傳導的分子動力學模擬［J］.東南大學學報，2003，33(4):468 －470.

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

［3］SHIGEO M.Molecular dynamic method for microscale heat transfer［J］.Advances in Numerical Heat Transfer，2000，2(6):189－226.

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

［5］MAITI A，MAHAN G D，PANTELIDES S T.Dynamical simulations of nonequilibrium processes—heat flow and the kapitza resistance across grain boundaries［J］.Solid Communications，1997，102(7):17－21.

［6］MARUYAMA S，KIMURA T.A study on thermal resistance over a solid-liquid interface by molecular dynamics method［J］.Thermal Science Engineering，1999，7:63－68

［7］TWU C J，HO J R.Molecular dynamics study of energy flow an the kapitza conductance across an interface with imperfection formed by two dielectric thin films［J］.Phys Rev(B)，2003，67(20):205422

［8］李博翰，江建軍，朱玲，等.硅鍺超晶格薄膜界面熱傳導的分子動力學模擬［J］.功能材料與器件學報，2007，13(3):293－296.

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

［10］STILLINGER F，WEBER T.Computer simulation of local order in con densed phases of Silicon［J］.Phy Rev(B)，1985，31:5262－5271.

［11］DING K，ANDERSEN H C.Molecular-dynamics simulation of amorphous germanium［J］.Phys Rev(B)，1986，34:6987－6991.

［12］LANDRY E S，MCGAUGHEY A J H.Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations［J］.Phys Rev(B)，2009，80:165304.

［13］ABRAMSON A R，TIEN C L，MAJUMDAR A.Interface and strain effects on the thermal conductivity of heterostructures:a molecular dynamics study［J］.J Heat Transfer，2002，124:963 －967.

［14］SWARTZ E T，POHLR O，Thermal boundary resistance［J］.Rev Modern Phys，1989，61:605 －668.

［15］PRASHER R S，PHELAN P E.A scattering-mediated acoustic mismatch model for the prediction of thermal boundary resistance［J］.J Heat Transfer，2001，123:105－112.

［16］BORCA T，LIU W，LIU J，et al.Thermal conductivity of symmetrically strained Si/Ge superlattices［J］.Superlattices Microstruct，2000，28:199 －206.