基于分子動(dòng)力學(xué)對(duì)超晶格結(jié)構(gòu)界面熱阻的模擬研究

張興麗， 孫兆偉

(哈爾濱工業(yè)大學(xué)衛(wèi)星技術(shù)工程研究所，哈爾濱 150001)

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

界面熱阻是熱載流子在兩接觸固體界面層相互作用的結(jié)果。在微觀傳熱領(lǐng)域內(nèi)，不同材料間界面熱阻的計(jì)算是當(dāng)前研究的熱點(diǎn)問(wèn)題，因?yàn)榻缑鏌嶙柚苯佑绊懙讲牧系臒醾鲗?dǎo)性能，從而對(duì)微/納米器件的設(shè)計(jì)和熱優(yōu)化產(chǎn)生影響。利用分子動(dòng)力學(xué)方法對(duì)界面熱阻進(jìn)行計(jì)算機(jī)模擬，可以從不同角度對(duì)界面熱阻的內(nèi)在機(jī)理進(jìn)行探討。分子動(dòng)力學(xué)(MD)模擬通過(guò)求解牛頓運(yùn)動(dòng)方程得到每個(gè)粒子空間位置和運(yùn)動(dòng)狀態(tài)隨時(shí)間的演進(jìn)狀況，在計(jì)算界面熱阻時(shí)，它不需要考慮每個(gè)粒子本身的散射特性，而只需根據(jù)勢(shì)能函數(shù)確定粒子之間的相互作用規(guī)律［5～8］。

本工作以Si/Ge超晶格結(jié)構(gòu)為例，利用非平衡態(tài)分子動(dòng)力學(xué)模擬方法從微觀機(jī)制出發(fā)研究超晶格結(jié)構(gòu)界面熱阻的一些變化趨勢(shì)。

1 分子動(dòng)力學(xué)模型

本工作所采用的超晶格結(jié)構(gòu)導(dǎo)熱模型如圖1所示。模擬對(duì)應(yīng)的硅的晶格常數(shù)為0.543nm，鍺的晶格常數(shù)為0.5657nm。在X方向上布置隨機(jī)恒溫?zé)釅σ越崃鞣较虻臏囟忍荻?，高溫?zé)釅偷蜏責(zé)釅Φ牟牧暇c各自臨近的超晶格材料相同，并設(shè)定其厚度為3UC(UC，晶格長(zhǎng)度);在Y，Z方向施加周期性邊界條件，由于垂直熱流方向的橫截面積過(guò)小會(huì)對(duì)熱導(dǎo)率的計(jì)算結(jié)果產(chǎn)生誤差［9］，因此設(shè)定YOZ橫截面積為4UC×4UC;模型的最外層設(shè)置厚度為2UC的絕熱壁，它的作用是減少導(dǎo)熱層內(nèi)的粒子蒸發(fā)，防止與外界產(chǎn)生熱量交換，并且設(shè)定該區(qū)域粒子的速率為0。

圖1 硅鍺超晶格結(jié)構(gòu)非平衡分子動(dòng)力學(xué)模擬模型Fig.1 The NEMD simulation model of Si/Ge superlattice structure

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

式中，[]...表示在總的模擬時(shí)間內(nèi)的統(tǒng)計(jì)平均;kB是Boltzmann常數(shù);Nj為第j層的粒子數(shù)。依據(jù)經(jīng)典統(tǒng)計(jì)下的能量均分定理，第j層中的能量等式可以表示為:

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

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

通過(guò)數(shù)值求解式(3)，得到與分子動(dòng)力學(xué)模擬(MD)的局域溫度Tj，MD相對(duì)應(yīng)的真實(shí)晶格溫度Tj。

當(dāng)高溫?zé)釅虻蜏責(zé)釅Φ牧Ｗ优c其他粒子作用時(shí)，溫度會(huì)發(fā)生改變。為使在導(dǎo)熱層區(qū)域形成穩(wěn)定的溫度梯度，需要通過(guò)改變粒子的速率來(lái)增加高溫?zé)釅σ欢〝?shù)量的動(dòng)能，同時(shí)從低溫?zé)釅σ瞥鐾葦?shù)量的動(dòng)能。高、低溫?zé)釅δ芰康淖兓靠梢员硎緸?

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

式中A為熱流方向的橫截面積;τ為模擬的時(shí)間步長(zhǎng)。

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

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

圖2 硅鍺超晶格結(jié)構(gòu)界面示意圖Fig.2 Schematic diagram of the Si/Ge superlattice interface

3 結(jié)果和討論

非平衡分子動(dòng)力學(xué)模擬在微正則(NVE)系統(tǒng)條件下進(jìn)行，模擬的時(shí)間步長(zhǎng)為1fs，總的模擬步長(zhǎng)數(shù)為5×106，其中前2×106步使系統(tǒng)平衡。在不同的模擬溫度下，高低溫?zé)釅Φ臏囟仍O(shè)為Thot=T+20 K及Tcold=T－20 K以形成溫度梯度。

圖3為周期長(zhǎng)度為10UC的硅鍺超晶格薄膜在平均溫度為500K時(shí)從高溫?zé)釅Φ降蜏責(zé)釅?dǎo)熱層溫度示意圖，從圖中可以明顯看出受界面熱阻的影響在硅鍺界面處溫度存在明顯的跳躍。三個(gè)硅鍺交界面的溫度變化值ΔT分別約為19K，10K，8K，由式(6)可知靠近高溫?zé)嵩〉牡谝粋€(gè)界面熱阻要遠(yuǎn)遠(yuǎn)大于其余界面，這個(gè)模擬結(jié)果與Abramson等［13］的研究成果相一致。他指出界面熱阻與界面數(shù)之間變化的非線性行為，是由于不同界面間傳遞熱流的聲子振動(dòng)頻率不同致使熱流穿過(guò)每個(gè)界面產(chǎn)生的熱阻也不同。因此，最靠近高溫?zé)嵩〉慕缑鏌嶙璞绕溆嘟缑娲?，這也與傳遞熱流的聲子類型不同有關(guān)。

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

圖4是溫度為400K和500K時(shí)利用分子動(dòng)力學(xué)模擬計(jì)算的最靠近高溫?zé)嵩〉腟i/Ge界面熱阻隨周期長(zhǎng)度的變化。由圖4可知隨著周期長(zhǎng)度的增大，該界面熱阻逐漸減小，因此界面熱阻在總熱阻中所占的比例逐漸減小，界面效應(yīng)減弱，Si/Ge超晶格薄膜的熱傳導(dǎo)性能不斷提高(如圖5所示)。

圖4 最靠近高溫?zé)嵩〉腟i/Ge界面熱阻與周期長(zhǎng)度的變化關(guān)系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之間時(shí)，周期長(zhǎng)度為10UC的Si/Ge超晶格薄膜界面熱阻模擬結(jié)果如圖6所示。從圖中可以明顯看出不同溫度下薄膜的界面熱阻也是不同的，隨著溫度的逐漸升高，界面熱阻越來(lái)越小。圖中比較了分子動(dòng)力學(xué)模擬、散射失配理論(Diffuse Mismateh Model，DMM)研究［14，15］以及實(shí)驗(yàn)數(shù)據(jù)［16］三種不同方法的計(jì)算結(jié)果，可以看出他們符合得比較好，變化趨勢(shì)大致相同。因?yàn)閷?shí)驗(yàn)過(guò)程中試樣的誤差以及DMM理論本身的適用條件限制，二者得到的界面熱阻結(jié)果比本工作利用分子動(dòng)力學(xué)模擬結(jié)果略大。界面熱阻隨溫度升高而減小的原因主要是因?yàn)闇囟鹊纳呤菇缑嫣幇l(fā)生非彈性散射的概率增加。在高溫下界面處的高頻聲子只能發(fā)生非彈性散射，轉(zhuǎn)化為低頻聲子才可以傳遞能量。因此聲子的非彈性散射增加了界面的熱傳導(dǎo)能力，界面熱阻也隨之下降［7］。

4 結(jié)論

(1)受界面熱阻機(jī)制的影響，在超晶格導(dǎo)熱區(qū)域的界面處會(huì)發(fā)生溫度的突變，并且在最靠近高溫?zé)嵩〉慕缑鏈囟韧蛔冏顬槊黠@。因此，最靠近高溫?zé)嵩〉慕缑鏌嶙鑼?duì)整個(gè)結(jié)構(gòu)的熱傳導(dǎo)能力起著決定性作用。

(2)界面熱阻會(huì)隨超晶格結(jié)構(gòu)周期長(zhǎng)度的增大而逐漸減小，超晶格結(jié)構(gòu)的導(dǎo)熱能力會(huì)隨之相應(yīng)提高。

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

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

［2］吳國(guó)強(qiáng)，孔憲仁，孫兆偉，等.單晶硅薄膜法向熱導(dǎo)率的分子動(dòng)力學(xué)模擬［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào).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］李博翰，江建軍，朱玲，等.硅鍺超晶格薄膜界面熱傳導(dǎo)的分子動(dòng)力學(xué)模擬［J］.功能材料與器件學(xué)報(bào)，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.