劉英光 邊永慶 韓中合

(華北電力大學(xué)能源動(dòng)力與機(jī)械工程學(xué)院, 保定 071003)

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

1 引 言

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

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

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

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

2 模型與計(jì)算方法

2.1 模 型

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

表1 對(duì)稱傾斜角度與對(duì)應(yīng)的Σ值Table 1.List of the symmetrical tilt angles and the corresponding Σ values.

2.2 計(jì)算方法

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

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

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

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

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

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

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

基于傅里葉定律的有效熱導(dǎo)率 κeff的表達(dá)式為

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

3 模擬結(jié)果和討論

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

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

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

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

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

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

式中, 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)出很好的一致性, 進(jìn)一步驗(yàn)證模擬結(jié)果的正確性.

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

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

圖5 300 K時(shí)卡皮查熱阻與樣本長(zhǎng)度的關(guān)系Fig.5.Kapitza resistance as a function of sample length at 300 K.

圖6 熱導(dǎo)率與樣本長(zhǎng)度的關(guān)系Fig.6.Thermal conductivity as a function of sample length.

其中 c 為 比熱 容, v 為 聲 子 群速 度, Λb為 整體 平均自由程.由于材料結(jié)構(gòu)變化對(duì)比熱容和群速度的影響很小, 晶界結(jié)構(gòu)主要影響聲子的散射, 我們專注于平均自由程對(duì)尺寸的依賴性.當(dāng)樣本長(zhǎng)度遠(yuǎn)小于聲子平均自由程時(shí), 此時(shí)聲子在雙晶中處于彈道輸運(yùn)模式, 聲子?晶界散射作為主要散射機(jī)制, 對(duì)聲子的抑制作用明顯, 散射強(qiáng)度隨長(zhǎng)度的增加而減小;當(dāng)樣本長(zhǎng)度遠(yuǎn)大于聲子平均自由程時(shí), 此時(shí)聲子處于擴(kuò)散輸運(yùn)模式, 聲子?聲子散射起主導(dǎo)作用, 但由于其散射強(qiáng)度較弱, 對(duì)聲子傳輸幾乎沒有阻礙作用.

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

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

圖8 5.45°和38.95°晶界結(jié)構(gòu)的熱導(dǎo)率隨溫度的變化Fig.8.Temperature?dependent thermal conductivity for the 5.45° and 38.94° GBs.

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

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

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

4 結(jié) 論

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

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

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

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

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