包含傾斜晶界的雙晶ZnO的熱輸運行為*

劉英光 邊永慶 韓中合

(華北電力大學能源動力與機械工程學院, 保定 071003)

用嵌入位錯線法和重合位置點陣法構建含有小角度和大角度傾斜角的雙晶氧化鋅納米結構.用非平衡分子動力學方法模擬雙晶氧化鋅在不同傾斜角度下的晶界能、卡皮查熱阻, 并研究了樣本長度和溫度對卡皮查熱阻和熱導率的影響.模擬結果表明, 晶界能在小角度區(qū)域隨傾斜角線性增加, 而在大角度區(qū)域達到穩(wěn)定,與卡皮查熱阻的變化趨勢一致.熱導率隨樣本長度的增加而增加, 卡皮查熱阻表現(xiàn)出相反的趨勢.然而隨著溫度的增加, 熱導率和卡皮查熱阻都減小.通過比較含5.45°和38.94°晶界樣本的聲子態(tài)密度, 發(fā)現(xiàn)聲子光學支對熱傳導的影響不大, 主要由聲子聲學支貢獻, 大角度晶界對聲子散射作用更強, 聲學支波峰向低頻率移動.

1 引 言

氧化鋅(ZnO)作為寬禁帶半導體材料, 具有優(yōu)異的壓電性、光電性和氣敏性等[1?3].隨著納米材料制備技術的提高, 微納尺度的ZnO在諸多領域有較為廣泛和有效的應用[4,5].由于較高的缺陷密度和表面積?體積比, ZnO納米裝置的熱效率和熱穩(wěn)定性已經(jīng)成為眾多科研工作者關注的主要問題.納米結構ZnO中存在大量的晶界[6], 當聲子的平均自由程大于平均晶粒尺寸時, 晶界對聲子的散射成為主要的散射機制.關于熱輸運的許多基本問題還沒有得到充分解釋, 仍然需要發(fā)展不同方法來深入理解熱輸運細節(jié).

目前, 國內(nèi)外研究者主要對不同類型納米結構ZnO的熱性質(zhì)開展研究.Zhang和Chen[7]討論了ZnO納米線的直徑、長度對熱導率的影響, 發(fā)現(xiàn)其具有明顯的尺寸效應, 同時引入合金、超晶格以及多孔結構等控制聲子的熱輸運行為.El?Brolossy等[8]采用溶劑熱法制備了摻雜鋁的ZnO納米顆粒, 當顆粒尺寸減小時, 比熱容明顯提高, 而鋁摻雜對比熱容的影響較小; 由于存在幾種嵌套的減熱機制, 所測得的熱導率比單晶ZnO小一個數(shù)量級左右.Giri等[9]合成了氧化鋅/對苯二酚(ZnO/HQ)納米結構材料, 研究晶界密度對聲子散射和熱傳導的影響, 結果表明晶界散射是導致熱導率顯著降低的主要原因.Liang和Shen[10]探索多晶ZnO的聲子和電子輸運, 并確定卡皮查熱阻和有效勢壘高度與晶界間隙有關; 當材料中加入硫化鋅(ZnS)納米薄膜時, 卡皮查熱阻增加了三倍, 勢壘區(qū)域擴大了兩倍.這些研究表明納米結構ZnO的熱導率主要受到晶界對聲子的散射作用的影響, 因此有必要對晶界如何影響熱性質(zhì)進行深層次, 多方面探討.

聲學失配模型[11]和擴散失配模型[12]通常被用于評估界面熱阻.然而, 這兩個理論都有特定的應用范圍并且不能考慮晶界的原子結構對納米材料熱傳導的影響.同時實驗上需要花費大量的時間和精力得到需要的晶界結構, 例如具有特定傾斜角度的晶界, 使得定性地研究晶界對熱傳導的影響變得非常困難.分子動力學(MD)模擬克服了傳統(tǒng)的理論和實驗方法的缺陷, 廣泛地應用于設計不同晶界結構和研究晶界上的聲子散射[13].我們嘗試采用這種方法去控制一定形式的晶界并得到相應的晶界性質(zhì)與卡皮查熱阻和熱導率之間的關系.

本文采用非平衡分子動力學(NEMD)模擬研究不同傾斜角度的晶界對納米結構ZnO的熱性質(zhì)的影響.基于NEMD模擬和Frank?Bilby公式分別計算晶界能和位錯密度.選擇幾組小角度和大角度晶界作為對比, 研究在不同樣本長度和溫度下卡皮查熱阻和熱導率的變化.最后利用聲子態(tài)密度深入理解聲子?晶界相互作用, 進一步解釋傾斜角度對熱阻的影響.通過設計特定的晶界結構為獲得想要的熱性質(zhì)提供可能.

2 模型與計算方法

2.1 模 型

雙晶模型通常被用于研究單一晶界的相關特征.兩個晶粒平行放置在界面兩邊, 將左邊的晶粒沿Z軸順時針旋轉 θL角度, 將右邊的晶粒逆時針旋轉 θR角度, 定義晶界傾斜角 θM= θL+ θR.理論上傾斜角可以設定為任意值, 但是會產(chǎn)生很多原子排列混亂的缺陷晶界, 因此采用重合位置點陣方法構建合理的晶界來滿足模擬要求[14].在經(jīng)過旋轉后, 原始的兩個晶粒在晶界處具有一定數(shù)量的重合位置點, 重合位置點與晶格點比例的倒數(shù)用Σ表示.根據(jù)傾斜角, 晶界分為小角度和大角度晶界.小角度晶界由沿[001]方向的邊界位錯構成, 構建了5個從 5.45°到30.51°小角度晶界的雙晶ZnO.小角度晶界處的位錯排列不適用于大角度晶界, 基于重合位置點陣法構建了5個從36.86°到67.38°大角度晶界的雙晶ZnO.所選晶界傾斜角與對應的Σ值見表1.此外, 根據(jù)非線性共軛梯度法對樣本進行能量最小化, 得到最穩(wěn)定的結構.選擇5.45°小角度晶界和38.94°大角度晶界作為對比, 采用OVITO軟件將其可視化[15], 如圖1所示.可以清楚地看出5.45°晶界由位錯線構成, 而38.94°晶界含有大量缺陷.

表1 對稱傾斜角度與對應的Σ值Table 1.List of the symmetrical tilt angles and the corresponding Σ values.

2.2 計算方法

NEMD方法的基本思想是在材料體系上施加熱流.本文采用的納米尺度傳熱試樣的原理模型如圖2所示.兩個晶粒以不同的取向放置在模擬盒子中, 中間位置包含晶界.為了避免邊界原子的相互作用, 在熱流傳輸方向上選擇兩層原子固定在模擬盒子兩端.包含五層原子的熱源和熱匯區(qū)域分別在固定原子層內(nèi)側.熱源區(qū)的溫度要高于熱匯區(qū)溫度, 使體系中產(chǎn)生從左到右的溫度梯度.在熱流傳輸方向采用固定邊界條件, 通過設定不同的樣本長度來研究卡皮查熱阻和熱導率的尺寸效應.橫截面方向采用周期性邊界條件, 防止模擬過程中原子的運動超出邊界而丟失, 保持系統(tǒng)中原子數(shù)量恒定,同時消除橫截面邊界對原子受力的影響.為了避免橫截面尺寸對結果產(chǎn)生干擾, 將其大小定為10 ×10 單位晶胞 (46.3 ? × 46.3 ?).

所有的模擬過程都是基于LAMMPS軟件進行的[16], 利用麥克斯韋?玻爾茲曼速度分布對體系進行初始化.Buckingham勢函數(shù)用于描述Zn—O和O—O的原子間相互作用, 設置參數(shù)下的模擬值與ZnO的實驗結果有很好的一致性[17], 其表達式為

其中對于 Zn2+—O2—化合鍵, A = 529.7 eV, B =0.03581 nm, C = 0; 對于 O2——O2—化合鍵, A =9547.96 eV, B = 0.021916 nm, C = 32 eV·A6.

圖2 NEMD模擬計算熱性質(zhì)的示意圖Fig.2.Schematic diagram of the NEMD model for calculat?ing the thermal properties.

首先在零溫條件下對體系進行能量最小化, 優(yōu)化原子位置; 然后在正則系綜(NVT)下進行1 ns的平衡態(tài)模擬來達到控制的溫度和體積; 最后在微正則系綜(NVE)下施加溫度梯度并測量熱流值.根據(jù)雙晶樣本的長度將其分割為一定數(shù)目的原子層.選擇熱匯區(qū)中最熱的原子和熱源區(qū)中最冷的原子進行動能交換獲得每個原子層中的平均溫度.動能交換的過程重復2 ns, 保證達到穩(wěn)定熱流, 其表達式為

式中t為模擬時間, A為樣本的橫截面積, n為動能交換次數(shù), m為Zn和O原子的平均質(zhì)量, vh和vc分別表示最熱和最冷原子的速度.

在穩(wěn)定狀態(tài)下晶界的兩側建立線性溫度梯度,然而在體系的兩端存在非線性區(qū)域.在晶界處溫度不連續(xù), 具有非常明顯的溫度跳躍 Δ T .卡皮查熱阻R為晶界處的溫度跳躍 Δ T 與熱流 J 之比,

基于傅里葉定律的有效熱導率 κeff的表達式為

式中 ? T/?x 為溫度梯度.

3 模擬結果和討論

為了研究晶界結構對納米ZnO熱輸運性質(zhì)的影響, 首先計算了晶界能和位錯密度, 它們是表征晶界特性的重要參數(shù).其中晶界能 EGB的表達式為

式中A為特征參數(shù), 表示為 A =1+ln(b/2πr0) , 其中 r0為位錯核半徑, b 為伯格斯矢量; E0為體系初始能量, 表示為 E0=Gb/4π(1-v) , 其中 G 為剪切模量, v 為泊松比, ZnO的剪切模量 G =0.75 , 泊松比 v =0.337 .根據(jù)此模型計算得到的晶界能滿足MD的模擬值, 并且已經(jīng)在Al2O3等材料中得到證實[18].

位錯密度 ρ 為特定傾斜角晶界單位長度的位錯數(shù)量, 其表達式為[19]

圖3 晶界能和位錯密度與傾斜角的關系Fig.3.The GB energy and dislocation density as a func?tion of misorientation angle.

式中D為位錯間距.對于ZnO的六角纖鋅礦結構,其晶格常數(shù)為a = 3.28 ?, 伯格斯矢量的大小與晶格常數(shù)相等.圖3顯示位錯密度隨傾斜角的變化與晶界能具有相同的趨勢.當傾斜角大于過渡角時,近似地假設位錯密度為固定值1.28 nm—1.這是由于在小角度區(qū)域, 每個晶界由一系列獨立的位錯組成.位錯之間的間隙隨著傾斜角的增加而減小, 位錯間距減小.然而在大角度區(qū)域時, 位錯之間相互重疊, 沒有間隙, 位錯間距不隨傾斜角的改變而改變, 可以認為位錯密度保持不變.

雙晶ZnO的卡皮查熱阻隨傾斜角的變化關系如圖4所示.卡皮查熱阻在小角度區(qū)域中幾乎線性增加, 而在大角度區(qū)域趨于穩(wěn)定.這與晶界能隨傾斜角的變化非常類似, 表明卡皮查熱阻與晶界能之間有緊密的關聯(lián).這一發(fā)現(xiàn)與之前UO2[20]和金剛石[21]的研究結果一致.Read?Shockley的擴展模型能用于描述卡皮查熱阻與傾斜角之間的關系, 其表達式為

式中, Rc和Rs是可調(diào)參數(shù), 通過擬合MD數(shù)據(jù)可以確定= 3.474 m2·K/GW,= 0.348 m2·K/GW,= —0.214 m2K/GW,= 0.146 m2K/GW.圖4顯示擬合曲線與模擬值表現(xiàn)出很好的一致性, 進一步驗證模擬結果的正確性.

圖4 卡皮查熱阻與傾斜角的關系Fig.4.Kapitza resistance as a function of misorientation angle.

在不同樣本長度條件下研究雙晶ZnO的卡皮查熱阻和熱導率的尺寸效應.作為對比, 分別選擇包含 5.45°, 20.02°, 38.94°和 60.0°晶界的樣本, 樣本平均長度在23.2—185.2 nm之間, 平均模擬溫度設定為300 K.從圖5可以看出, 當樣本長度在23.2—92.6 nm之間時, 卡皮查熱阻隨長度的增加急劇減小, 具有明顯的尺寸效應.結合圖6可以看出, 熱導率在此長度范圍內(nèi)顯著增加, 而當長度大于92.6 nm時, 卡皮查熱阻和熱導率的變化趨于穩(wěn)定.通過線性外推法計算得到單晶ZnO的熱導率為98.5 W/mK, 含有晶界樣本的熱導率的最大值仍小于單晶熱導率.聲子是ZnO中主要的熱載流子, 基于經(jīng)典動力學理論, 晶格熱導率的表達式為

圖5 300 K時卡皮查熱阻與樣本長度的關系Fig.5.Kapitza resistance as a function of sample length at 300 K.

圖6 熱導率與樣本長度的關系Fig.6.Thermal conductivity as a function of sample length.

其中 c 為 比熱 容, v 為 聲 子 群速 度, Λb為 整體 平均自由程.由于材料結構變化對比熱容和群速度的影響很小, 晶界結構主要影響聲子的散射, 我們專注于平均自由程對尺寸的依賴性.當樣本長度遠小于聲子平均自由程時, 此時聲子在雙晶中處于彈道輸運模式, 聲子?晶界散射作為主要散射機制, 對聲子的抑制作用明顯, 散射強度隨長度的增加而減小;當樣本長度遠大于聲子平均自由程時, 此時聲子處于擴散輸運模式, 聲子?聲子散射起主導作用, 但由于其散射強度較弱, 對聲子傳輸幾乎沒有阻礙作用.

卡皮查熱阻和熱導率還會受到溫度的影響.選擇含有5.45°和38.94°晶界的雙晶樣本, 設定長度為23.2 nm.在300—700 K的溫度范圍內(nèi)計算得到的結果如圖7和圖8所示.從圖7中看到, 38.94°晶界樣本的卡皮查熱阻隨溫度的增加顯著減小, 而對于5.45°晶界樣本其變化曲線相對平坦.同時在不同溫度下, 38.94°晶界樣本的卡皮查熱阻總是大于5.45°晶界樣本對應值.卡皮查熱阻主要來源于聲子在晶界處的散射, 當溫度升高時穿過晶界的聲子數(shù)量增加, 有更多的聲子參與熱傳導過程, 導致卡皮查熱阻減小.圖8顯示, 對于5.45°晶界樣本其熱導率同樣隨溫度的增加而減小, 而對于38.94°晶界樣本, 當溫度小于400 K時, 熱導率呈增加趨勢,然后迅速減小.這是由于在低溫區(qū)域, 被激發(fā)的聲子數(shù)量較少, 晶格的振動受到強烈限制, 熱導率主要受晶界對長波聲子散射的影響.隨著溫度的增加, 具有長波長、低頻率的聲子在晶界處的穿透率大于散射率, 表現(xiàn)為熱導率增加.在高溫區(qū)域, 材料發(fā)生熱軟化并且晶格劇烈振動, 激發(fā)出大量的低波長、高頻率聲子.聲子振動幅度隨溫度的增加而增加, 高頻率聲子之間的散射作用占據(jù)主導, 使聲子的平均自由程減小, 從而導致熱導率進一步降低.

圖7 5.45°和38.95°晶界結構的卡皮查熱阻隨溫度的變化Fig.7.Temperature?dependent Kapitza resistance for the 5.45° and 38.94° GBs.

圖8 5.45°和38.95°晶界結構的熱導率隨溫度的變化Fig.8.Temperature?dependent thermal conductivity for the 5.45° and 38.94° GBs.

晶界處的原子排列對聲子模式有重要的影響.為了深層次理解晶界對卡皮查熱阻和熱導率的影響, 基于速度自相關函數(shù)的傅里葉轉換計算得到聲子態(tài)密度, 其表達式為[22]

式中 v (t) 為原子在時間t時的速度, v (0) 為原子初速度, ω 為聲子振動頻率.比較5.45°和38.94°晶界樣本的聲子態(tài)密度, 結果如圖9所示, 設定固定長度為92.6 nm, 溫度為300 K.在ZnO中存在4個聲子振動模式, 分別是低頻率的橫向聲學分支(TA)、縱向聲學分支(LA)和高頻率的橫向光學分支(TO)、縱向光學分支(LO).兩個樣本的聲子態(tài)密度在大于9 THz時幾乎是重合的, 因為光學聲子的群速度非常小, 對晶體內(nèi)部的熱傳導影響可以忽略.然而在1—7 THz頻率范圍, 聲學模式存在顯著差異.相比于5.45°晶界樣本, 38.94°晶界樣本的TA和LA波峰移動到更低頻率的位置, 表明具有相同頻率的聲子能夠通過5.45°晶界, 而在38.94°晶界處發(fā)生散射, 轉變?yōu)楦皖l率的聲子.聲子模式失配現(xiàn)象進一步說明大角度晶界的原子排列比小角度晶界更加混亂, 對聲子散射作用更強, 最終造成大角度晶界樣本的卡皮查熱阻大于小角度晶界的卡皮查熱阻.

圖9 5.45°和38.94°晶界結構的聲子態(tài)密度比較Fig.9.Phonon density of states for the 5.45° and 38.94°GBs.

4 結 論

本文利用NEMD方法對含有傾斜晶界的雙晶ZnO的熱輸運性質(zhì)進行模擬, 系統(tǒng)地分析了晶界傾斜角、樣本長度以及溫度對卡皮查熱阻和有效熱導率的影響, 研究發(fā)現(xiàn)：

1)在傾斜角小于36.86°時, 晶界能和位錯密度隨傾斜角線性增加, 當傾斜角大于36.86°時, 晶界能和位錯密度幾乎穩(wěn)定不變;

2)卡皮查熱阻隨傾斜角的變化趨勢與晶界能和位錯密度的變化趨勢相同, 證明它們之間有較強的相關性, 基于Read?Shockley的擴展模型的計算值與模擬值擬合良好;

3)卡皮查熱阻和有效熱導率具有明顯的尺寸效應和溫度效應, 在小角度區(qū)域受傾斜角的影響較大, 而在大角度區(qū)域受傾斜角的影響較小;

4)通過計算聲子態(tài)密度得到光學分支幾乎對熱輸運沒有影響, 主要由聲學分支貢獻.大角度晶界的聲子?晶面散射作用更強, 聲子頻率降低, 聲學分支的波峰向低頻移動.