氮化硼納米管熱輸運(yùn)性能的分子動(dòng)力學(xué)模擬

高宇飛 孟慶元 張 璐 劉甲秋 荊宇航,*

(1哈爾濱工業(yè)大學(xué)航天學(xué)院,哈爾濱150001;2中國(guó)醫(yī)科大學(xué),沈陽(yáng)110001;3哈爾濱玻璃鋼研究院,哈爾濱150001)

氮化硼納米管熱輸運(yùn)性能的分子動(dòng)力學(xué)模擬

高宇飛1孟慶元1張 璐2劉甲秋3荊宇航1,*

(1哈爾濱工業(yè)大學(xué)航天學(xué)院,哈爾濱150001;2中國(guó)醫(yī)科大學(xué),沈陽(yáng)110001;3哈爾濱玻璃鋼研究院,哈爾濱150001)

采用基于聲子散射理論的Boltzmann-Peierls聲子傳輸方程(BTE)和非平衡態(tài)分子動(dòng)力學(xué)模擬(NEMD)方法研究了氮化硼納米管(BNNT)的熱輸運(yùn)性能.分析了BNNT的熱力耦合效應(yīng),通過(guò)BTE與NEMD兩種方法相結(jié)合,分析了溫度和長(zhǎng)度對(duì)BNNT熱輸運(yùn)性能的影響,并應(yīng)用量子修正擴(kuò)大了NEMD的研究范圍.結(jié)果表明:隨著拉伸或壓縮應(yīng)變的增加,BNNT熱輸運(yùn)性能均呈降低的趨勢(shì).通過(guò)計(jì)算聲子態(tài)密度(PDOS)在理論上分析了以上結(jié)果,發(fā)現(xiàn)在拉伸狀態(tài)下,聲子模式的變化是決定BNNT熱輸運(yùn)性能變化的主要因素;在壓縮狀態(tài)下,熱導(dǎo)率變化是由于模型發(fā)生明顯的屈曲變形引起的.在低溫段,BNNT的熱輸運(yùn)性能受量子效應(yīng)影響最初有一個(gè)線性增加的過(guò)程,當(dāng)溫度超過(guò)一定值時(shí),其開(kāi)始顯著地降低;當(dāng)BNNT長(zhǎng)度小于120 nm時(shí),隨著長(zhǎng)度的增加,其彈道性能逐漸減弱,但仍主要體現(xiàn)為彈道特征,其熱導(dǎo)率(κ)與長(zhǎng)度(L)基本滿足κ∝Lα這一關(guān)系.

氮化硼納米管;熱輸運(yùn)性能;Boltzmann-Peierls聲子傳輸方程;非平衡態(tài)分子動(dòng)力學(xué)模擬;聲子態(tài)密度;量子修正

1 引言

納米材料因其獨(dú)特的結(jié)構(gòu)和優(yōu)良的力學(xué)、電學(xué)、熱學(xué)、光學(xué)等性能,表現(xiàn)出廣闊的應(yīng)用前景.至1991年發(fā)現(xiàn)碳納米管(CNT)以來(lái),1各種不同種類的納米管相繼被發(fā)現(xiàn),對(duì)它們的研究也隨之展開(kāi).

1994年Rubio等2預(yù)測(cè)到一定存在氮化硼納米管(BNNT).之后人們?cè)趯?shí)驗(yàn)中得到了多壁BNNT. 1996年Loiseau等3采用電弧放電法得到了單壁BNNT.有關(guān)BNNT及氮硼化合物的研究,如電子結(jié)構(gòu)和力學(xué)性能的研究已經(jīng)展開(kāi),4但有關(guān)其熱學(xué)性能的研究卻較少,在力學(xué)和熱學(xué)交叉領(lǐng)域,如熱力耦合效應(yīng)方面還未見(jiàn)報(bào)道.BNNT作為一種典型的納米材料,它的很多結(jié)構(gòu)參數(shù)都對(duì)其熱導(dǎo)率,尤其是軸向熱導(dǎo)率有著顯著的影響,5因此,可以利用BNNT的這種優(yōu)良特性制成可控制的導(dǎo)熱設(shè)備,這與傳統(tǒng)的通過(guò)選擇不同的材料來(lái)實(shí)現(xiàn)對(duì)器件導(dǎo)熱性能控制的方法相比有了革命性的突破.因此有必要對(duì)BNNT的熱傳導(dǎo)機(jī)理和導(dǎo)熱性能進(jìn)行系統(tǒng)的研究.

目前針對(duì)BNNT導(dǎo)熱性能的研究成果指出:對(duì)BNNT的熱輸運(yùn)起主要貢獻(xiàn)的是聲子,這與CNT的熱輸運(yùn)機(jī)理不同(對(duì)CNT熱導(dǎo)率起主要作用的是聲子和電子).BNNT的熱導(dǎo)率低于CNT,但BNNT比CNT有著更強(qiáng)的熱穩(wěn)定性和化學(xué)穩(wěn)定性;Liew和Yuan6通過(guò)分析CNT和BNNT混合物發(fā)現(xiàn),BNNT具有更好的穩(wěn)定性,較之CNT更適宜做管狀結(jié)構(gòu)的外壁;然而也有部分研究7指出,在低溫段,BNNT的熱導(dǎo)率要高于CNT.此外,同位素對(duì)BNNT熱導(dǎo)率有著顯著的影響,Zettl等8,9的研究指出,當(dāng)BNNT中含有高濃度的B同位素時(shí),其熱導(dǎo)率可與CNT相比;純的11B同位素可將BNNT熱導(dǎo)率提升50%;10而B(niǎo)arman11的研究表明,同位素10B會(huì)削弱氮化硼晶體的熱導(dǎo)率.BNNT具有較好的導(dǎo)熱性能,在室溫下,外直徑為30-40 nm的BNNT的熱導(dǎo)率可達(dá)到350 W·m-1·K-1;Stewart等12通過(guò)計(jì)算表明(8,0) BNNT的熱導(dǎo)率為200 W·m-1·K-1.Nakayama等13發(fā)現(xiàn)氮化硼納米帶可以增強(qiáng)部分金屬氧化物的導(dǎo)熱性能.Mingo和Broido14通過(guò)研究輸運(yùn)方程和彈道函數(shù)發(fā)現(xiàn)納米管的熱傳導(dǎo)性能表現(xiàn)出彈道特征.以上結(jié)果表明:BNNT在掃描探針顯微鏡、原子力顯微鏡的針尖,高溫陶瓷以及納米尺度熱量計(jì)和微電子處理器等方面有著很大的潛在應(yīng)用.但由于實(shí)驗(yàn)測(cè)量中的實(shí)際困難和材料在納米尺度下的誤差,目前理論和實(shí)驗(yàn)針對(duì)納米管導(dǎo)熱性能的定量研究結(jié)果具有較大的差別,如有部分研究結(jié)果15,16表明BNNT具有極低的熱導(dǎo)率(熱導(dǎo)率約1.5-14 W·m-1·K-1).總體來(lái)說(shuō),目前對(duì)于BNNT熱輸運(yùn)機(jī)理和性能的定量研究成果仍很缺乏,而且隨著聲子散射理論的逐步完善,一些對(duì)BNNT導(dǎo)熱性質(zhì)產(chǎn)生重要影響的較復(fù)雜的因素,如納米管的彈道傳輸與擴(kuò)散傳輸?shù)姆纸缧阅?、?yīng)變作用下動(dòng)態(tài)納米管結(jié)構(gòu)的熱輸運(yùn)性能等逐漸被提出,并逐漸成為研究的熱點(diǎn).

本文首先采用非平衡態(tài)分子動(dòng)力學(xué)模擬(NEMD)方法研究了BNNT在軸向應(yīng)變作用下的熱輸運(yùn)性能,并通過(guò)計(jì)算聲子態(tài)密度(PDOS)的方法對(duì)模擬結(jié)果進(jìn)行了理論上的補(bǔ)充和驗(yàn)證.進(jìn)一步通過(guò)求解Boltzmann-Peierls聲子傳輸方程(BTE)和NEMD相結(jié)合的方法分析了溫度和長(zhǎng)度對(duì)BNNT熱輸運(yùn)性能的影響,通過(guò)量子修正擴(kuò)大了NEMD的研究范圍.為研究BNNT彈道傳輸與擴(kuò)散傳輸?shù)姆纸缫约皠?dòng)態(tài)納米管結(jié)構(gòu)的熱輸運(yùn)性能提供了一些理論與計(jì)算模擬的參考數(shù)據(jù).

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

2.1 非平衡態(tài)分子動(dòng)力學(xué)模擬方法

本文在模擬過(guò)程中應(yīng)用Tersoff多體勢(shì)17來(lái)描述BNNT中原子間的相互作用.求解原子的運(yùn)動(dòng)方程時(shí)采用預(yù)測(cè)—校正算法,積分步長(zhǎng)為0.5 fs.邊界條件選擇為自由邊界條件.構(gòu)建BNNT模型時(shí),首先根據(jù)BNNT的空間結(jié)構(gòu)特點(diǎn)得到其初始構(gòu)型的空間排列,然后將BNNT的初始構(gòu)型在0 K下弛豫50000步,隨后逐漸升溫,每次升溫100 K,并弛豫相同的步數(shù),直到溫度升到1600 K,最后進(jìn)行退火處理來(lái)得到其穩(wěn)定結(jié)構(gòu).

計(jì)算模型的熱導(dǎo)率采用NEMD,18其基本思想為:將BNNT沿軸向等分為N層,設(shè)0層和N/2層為“冷區(qū)”和“熱區(qū)”(見(jiàn)圖1).每隔一定步數(shù)提取0層中動(dòng)能最大的原子和N/2層中動(dòng)能最小的原子,并交換兩者速度.一般情況下冷區(qū)中動(dòng)能最大的原子的動(dòng)能大于熱區(qū)中動(dòng)能最小的原子的動(dòng)能,這樣便產(chǎn)生了一個(gè)從0層到N/2層的熱流.又因?yàn)闊釁^(qū)的溫度高于冷區(qū)的溫度,因此便存在一個(gè)與熱流方向相反的溫差.基于上述模型,經(jīng)過(guò)統(tǒng)計(jì)可獲得軸向各層的溫度,進(jìn)一步計(jì)算可得溫度梯度(▽T),再結(jié)合能量交換過(guò)程中統(tǒng)計(jì)得到的熱流密度(J),由Fourier導(dǎo)熱定律J=-κ?T,可計(jì)算得到BNNT的熱導(dǎo)率(κ).整個(gè)非平衡態(tài)分子動(dòng)力學(xué)模擬過(guò)程是通過(guò)Lammps開(kāi)源軟件包19實(shí)現(xiàn)的.

圖1 非平衡態(tài)分子動(dòng)力學(xué)模擬過(guò)程示意圖Fig.1 Schematic diagram of NEMD process

圖2為BNNT在300 K溫度下軸向方向的溫度分布圖,其中各層溫度是通過(guò)統(tǒng)計(jì)各層中所有原子的動(dòng)能,通過(guò)公式(1)獲得

其中N為層中的原子數(shù),T代表溫度,kB為Boltzmann常量,mi和vi分別為模型中i原子的質(zhì)量和速度.從圖中可發(fā)現(xiàn)中間各層的溫度呈線性分布,而在冷層和熱層附近區(qū)域,則有明顯的非線性溫度響應(yīng),這是由于熱源和熱槽的強(qiáng)相互作用產(chǎn)生的.20因此在計(jì)算BNNT熱導(dǎo)率過(guò)程中,需要統(tǒng)計(jì)中間各層的溫度,進(jìn)而獲得BNNT軸向方向的溫度梯度以及熱導(dǎo)率.

圖2 BNNT軸向(L)的溫度分布Fig.2 Temperature distribution of BNNT along axial direction(L)

模擬的整個(gè)過(guò)程分為五步:首先用Conjugate Gradient(CG)方法弛豫整個(gè)構(gòu)型,以獲得BNNT的穩(wěn)定結(jié)構(gòu);然后將納米管兩端固定,在NPT系綜下運(yùn)行10000步(5 ps)標(biāo)定壓力使其處于最小狀態(tài);再在最小壓力狀態(tài)下運(yùn)行100000步(50 ps),使系統(tǒng)達(dá)到穩(wěn)定狀態(tài);然后在NVT系綜下運(yùn)行200000步(100 ps)來(lái)弛豫整個(gè)系統(tǒng);最后在NVT系綜下運(yùn)行1000000步(500 ps)統(tǒng)計(jì)熱流密度和溫度梯度,計(jì)算得出BNNT的熱導(dǎo)率(在最后一步過(guò)程中每100000步統(tǒng)計(jì)一次熱導(dǎo)率,通過(guò)計(jì)算10個(gè)熱導(dǎo)率數(shù)值的偏差分布可獲得模型熱導(dǎo)率的不確定度).

2.2 Boltzmann-Peierls聲子傳輸方程方法

微觀固體中的熱量主要是靠晶格振動(dòng)產(chǎn)生的聲子傳輸?shù)?聲子是一些被抽象出的離散粒子,單一聲子的能量為?ω,?=h/2π,h為普朗克常量,ω為頻率.Boltzmann-Peierls聲子傳輸方程14,21(BTE)可以完整地描述固體中聲子間的相互作用.但獲得BTE的精確解是困難的,目前可以通過(guò)一些近似的方法來(lái)獲得不同模型對(duì)應(yīng)的BTE的近似解.求解BTE的基本思想是通過(guò)不同聲子的散射時(shí)間和處于非平衡態(tài)聲子的分布,并結(jié)合模型的邊界條件來(lái)獲得模型熱流與溫度梯度,最后求得模型的熱導(dǎo)率.

2.3 低溫段量子修正

一般情況下,經(jīng)典分子動(dòng)力學(xué)模擬(MD)方法在溫度接近和高于德拜溫度時(shí)被認(rèn)為是有效的,因?yàn)樵诖藴囟葏^(qū)域內(nèi)所有的聲子模式都被完全激發(fā)出來(lái).與此相對(duì)的,當(dāng)溫度顯著地低于德拜溫度時(shí),大量的聲子沒(méi)有被激發(fā),模型的量子效應(yīng)不能被忽略,此時(shí)經(jīng)典的MD方法無(wú)效.為突破這一限制,量子修正方法(QCs)22-25被提出,目前它被廣泛應(yīng)用于熱力學(xué)和力學(xué)的研究中.

QCs的基本思想是由經(jīng)典力學(xué)系統(tǒng)中的模擬溫度來(lái)計(jì)算量子系統(tǒng)中對(duì)應(yīng)的修正溫度,根據(jù)經(jīng)典力學(xué)系統(tǒng)與量子系統(tǒng)能量相等,具體為系統(tǒng)的動(dòng)能等于聲子能量的一半這一原則.模擬溫度(TMD)與修正溫度的關(guān)系可表示為

其中,N0為平衡態(tài)聲子分布函數(shù),Tq為量子修正溫度,g(ω)為聲子態(tài)密度.結(jié)合模擬中獲得的TMD可計(jì)算Tq.通過(guò)統(tǒng)計(jì)不同時(shí)間的軸向熱流密度和溫度梯度,可獲得修正后的熱導(dǎo)率.

除了基于MD的QCs修正外,還存在兩種基于勢(shì)函數(shù)的量子修正方法,Feynman-Hibbs勢(shì)函數(shù)(FH勢(shì))和Wigner-Kirkwood勢(shì)函數(shù)(WK勢(shì))方法.FH勢(shì)是由量子均分方程的Feynman-Hibbs估算變量獲得的,是一種依賴于溫度的二體勢(shì),可以方便地植入分子動(dòng)力學(xué)或蒙特卡羅算法的計(jì)算程序中.應(yīng)用這種方法,Sese26計(jì)算了Lennard-Jones系統(tǒng)的熱力學(xué)和結(jié)構(gòu)特性,Tchouar等27應(yīng)用由FH方法演化的Lennard-Jones勢(shì)研究了氖的熱力學(xué)、結(jié)構(gòu)和傳輸特性,Goharshadi等28應(yīng)用Hartree-Fock dispersion近似勢(shì)函數(shù)和FH勢(shì)研究了氦的熱力學(xué)性質(zhì).WK勢(shì)是由對(duì)均分方程中德布羅意波長(zhǎng)的擴(kuò)展獲得的,可以被用來(lái)研究不同物理特性系統(tǒng)的量子效應(yīng).

3 計(jì)算結(jié)果與討論

3.1 熱力耦合效應(yīng)

本文選取20 nm的20_20扶手椅型BNNT(原子數(shù)為6320個(gè))為研究對(duì)象,分別對(duì)模型施加在0-0.1范圍的軸向拉伸和壓縮應(yīng)變.施加位移載荷過(guò)程中,應(yīng)變是通過(guò)重新標(biāo)定模型軸向原子坐標(biāo)實(shí)現(xiàn)的,具體為首先選擇BNNT軸向方向最外部?jī)蓪釉?分別對(duì)其施加反方向的軸向位移載荷,每施加完一次位移載荷讓系統(tǒng)弛豫一定步數(shù),使系統(tǒng)達(dá)到平衡狀態(tài).加載速率為每5000個(gè)步長(zhǎng)施加0.0025 nm的位移.在此期間統(tǒng)計(jì)不同時(shí)間BNNT軸向方向的熱流密度和溫度梯度,最終計(jì)算得到熱導(dǎo)率隨應(yīng)變變化趨勢(shì)如圖3所示.由圖3可發(fā)現(xiàn),在拉伸和壓縮應(yīng)變作用下,BNNT的軸向熱導(dǎo)率隨著應(yīng)變的增加而逐漸減小.

圖3 BNNT軸向熱輸運(yùn)性能隨軸向應(yīng)變的變化曲線Fig.3 Variation curve of axial thermal conductivity of BNNT under axial strain

圖4 BNNT在軸向拉伸應(yīng)變作用下的聲子態(tài)密度圖Fig.4 PDOS of BNNT under axial tensile strain

聲子態(tài)密度(PDOS)可以完整地描述各種模式聲子的分布,可用于研究動(dòng)態(tài)結(jié)構(gòu)材料的熱輸運(yùn)性能.PDOS是通過(guò)對(duì)速度自相關(guān)函數(shù)進(jìn)行傅里葉變換得到的.29-32速度自相關(guān)函數(shù)可以表示為

其中

對(duì)(4)式進(jìn)行傅里葉變換可得到PDOS的表達(dá)形式為

其中t代表時(shí)間,k=x,y,z代表模型的三個(gè)方向.通過(guò)上式,結(jié)合MD中統(tǒng)計(jì)的原子速度,可以獲得不同應(yīng)變作用下BNNT的聲子態(tài)密度分布圖.BNNT在拉伸和壓縮應(yīng)變作用下的聲子態(tài)密度分布如圖4、 5所示.

圖5 BNNT在軸向壓縮應(yīng)變作用下的聲子態(tài)密度高頻波峰圖Fig.5 High frequency peak in PDOS of BNNT under axial compressive strain

由圖4可見(jiàn),PDOS中有三個(gè)較高的波峰,其中前兩個(gè)位于低頻段,可被視為代表低頻聲子模式的波峰,后一個(gè)位于高頻段,被認(rèn)為是代表高頻聲子模式的波峰.隨著應(yīng)變的變化,相比于高頻波峰,兩個(gè)低頻波峰的變化較小.這意味著高頻聲子對(duì)應(yīng)變更加敏感.由于波峰的變化直接反映了聲子數(shù)量的增加或減少,以及振動(dòng)模式的增強(qiáng)或減弱,對(duì)比圖4和圖3我們可以發(fā)現(xiàn),高頻波峰的左移趨勢(shì)伴隨著模型熱導(dǎo)率的下降,因此可以認(rèn)為高頻聲子對(duì)模型熱輸運(yùn)性能有著主要影響.

在0-0.1的拉伸應(yīng)變作用下,隨著拉伸應(yīng)變的增加,從圖4可以看出高頻聲子波峰逐漸左移.這意味著高頻聲子數(shù)量逐漸減少,部分高頻聲子逐漸轉(zhuǎn)化為低頻聲子,它將導(dǎo)致聲子的平均頻率降低,因此此過(guò)程可以視為一個(gè)紅移過(guò)程.Xu和Buehler33研究CNT在軸向應(yīng)變作用下的熱輸運(yùn)性能,發(fā)現(xiàn)紅移過(guò)程將會(huì)減小聲子的速度.在聲子的彈道輸運(yùn)過(guò)程中,模型的熱導(dǎo)率(κ)與聲子速度有以下關(guān)系

其中c為比熱,l為平均自由程,v為聲子速度.因此圖4中高頻聲子的變化趨勢(shì)將導(dǎo)致BNNT的熱導(dǎo)率隨著拉伸應(yīng)變的增加而逐漸減小.觀察圖5可以發(fā)現(xiàn),在壓縮應(yīng)變作用下,不同應(yīng)變的PDOS變化較小,幾乎沒(méi)有明顯的變化.結(jié)合圖6中0.1的拉伸和壓縮應(yīng)變作用下BNNT的構(gòu)型圖可以發(fā)現(xiàn),在0.1的拉伸應(yīng)變作用下BNNT被逐漸拉長(zhǎng),但其沒(méi)有顯著的塑性變形,而在0.1的壓縮應(yīng)變作用下,BNNT發(fā)生了嚴(yán)重的塑性變形.隨著明顯的塑性變形的產(chǎn)生,系統(tǒng)的能量將顯著降低,進(jìn)而導(dǎo)致模型中聲子平均速度的降低,根據(jù)公式(6),在壓縮應(yīng)變下, BNNT的熱導(dǎo)率隨著應(yīng)變的增加而逐漸減小.這也說(shuō)明了在壓縮應(yīng)變下BNNT的熱導(dǎo)率發(fā)生變化的機(jī)理不同于拉伸應(yīng)變下的情況:在拉伸應(yīng)變下,是聲子模式的變化導(dǎo)致了BNNT熱導(dǎo)率的變化;在壓縮應(yīng)變下,是BNNT結(jié)構(gòu)的顯著變化導(dǎo)致了其熱導(dǎo)率的變化.

圖6 BNNT在(a)0.1的拉伸應(yīng)變作用和(b)0.1的壓縮應(yīng)變下的構(gòu)型圖Fig.6 Structures of BNNT under(a)0.1 tensile strain and (b)0.1 compressive strain

Buehler33和Zhang32等分別研究了CNT在應(yīng)變作用下導(dǎo)熱性能,本文在拉伸應(yīng)變作用下得到的PDOS與Buehler和Zhang的結(jié)果有著相同的變化趨勢(shì),高頻聲子波峰均隨著拉伸應(yīng)變的增加而逐漸左移.在壓縮應(yīng)變的作用下,Zhang計(jì)算了-0.02-0壓縮應(yīng)變下CNT的PDOS,發(fā)現(xiàn)其高頻聲子波峰有部分右移的趨勢(shì).本文認(rèn)為此處的不同是由于CNT的抗壓縮性能強(qiáng)于BNNT導(dǎo)致的.在-0.02的小壓縮應(yīng)變作用下CNT較之BNNT更不容易產(chǎn)生較大的塑性變形,因此CNT的能量減少將體現(xiàn)在聲子模式的變化上,而B(niǎo)NNT的能量減少體現(xiàn)在結(jié)構(gòu)的變化上.

3.2 溫度影響

本文應(yīng)用NEMD,結(jié)合低溫段溫度量子修正方法研究了10-1600 K溫度范圍內(nèi)長(zhǎng)度為20 nm的10_10扶手椅型BNNT(原子數(shù)為3160個(gè))的導(dǎo)熱性能.首先統(tǒng)計(jì)模擬過(guò)程中BNNT的溫度TMD,結(jié)合公式(2),可獲得量子修正溫度Tq.量子修正溫度Tq與模擬溫度TMD的關(guān)系見(jiàn)圖7.其中實(shí)線為T(mén)q隨TMD的變化曲線,虛線用來(lái)對(duì)比Tq隨TMD的變化率,斜率為1.從圖中可見(jiàn)在低溫段,?Tq/?TMD遠(yuǎn)大于1,隨著溫度的升高,其值逐漸減小,當(dāng)溫度達(dá)到1000 K時(shí),?Tq/?TMD≈1,這說(shuō)明此時(shí)Tq與TMD的變化率一致.結(jié)合Tq,可獲得BNNT的熱導(dǎo)率隨溫度的變化趨勢(shì)如圖8所示.

從圖8可以看出,BNNT的軸向熱導(dǎo)率在低溫段(10-70 K)隨著溫度的升高近似呈線性增長(zhǎng)趨勢(shì);當(dāng)溫度超過(guò)一定值(近似為100 K)時(shí),其熱導(dǎo)率隨著溫度的升高逐漸降低.

在低溫段,當(dāng)系統(tǒng)溫度遠(yuǎn)小于德拜溫度時(shí),由于溫度較低,晶格振動(dòng)較弱,被激發(fā)出的聲子較少.因此在低溫段,材料的邊界對(duì)長(zhǎng)波聲子的散射(也即邊界散射)是影響材料熱導(dǎo)率的主要因素.隨著系統(tǒng)溫度的降低,具有低頻率、長(zhǎng)波長(zhǎng)的聲子在聲子中的比重增加.根據(jù)聲子理論,長(zhǎng)波聲子較易與邊界發(fā)生散射,這就導(dǎo)致了聲子平均自由程減小,由(6)式可知,熱導(dǎo)率也隨之減小,因此在低溫段,BNNT的熱導(dǎo)率隨著溫度的增加而增加.在高溫段,晶格振動(dòng)劇烈,散發(fā)出大量的聲子,因此聲子與聲子間的相互碰撞即聲子間的散射取代邊界對(duì)長(zhǎng)波聲子的散射成為影響材料熱導(dǎo)率的主要因素.隨著溫度的增高,聲子數(shù)增多,聲子間的碰撞概率增大,聲子平均自由程減小,因此在高溫段,BNNT的熱導(dǎo)率隨著溫度的增加而逐漸減小.

圖7 模擬溫度(TMD)與量子修正溫度(Tq)的關(guān)系Fig.7 Relationship between simulated temperature(TMD) and quantum corrected temperature(Tq)The slope of dotted line is 1 and can be used to contrast the change rate of corrected temperature with simulated temperature.

圖8 BNNT軸向熱輸運(yùn)性能隨溫度的變化曲線Fig.8 Variation of axial thermal conductivity of BNNT under different temperatures

由圖8還可觀察到,在低溫段,修正后的熱導(dǎo)率與分子模擬直接獲得的熱導(dǎo)率有較大差別,隨著溫度的增加,這種差距逐漸減小.文獻(xiàn)23對(duì)比了量子修正前后硅晶體的熱導(dǎo)率變化曲線,當(dāng)溫度在800 K左右時(shí),兩者誤差約為10%.本文用兩種方法獲得的數(shù)據(jù)誤差為8%.文獻(xiàn)34用MD計(jì)算了1200-2000 K下(10,10)BNNT的熱導(dǎo)率,其熱導(dǎo)率隨著溫度升高逐漸降低,且1200 K時(shí)BNNT熱導(dǎo)率近似為48 W·m-1·K-1,1400 K時(shí)近似為43 W·m-1·K-1,1600 K時(shí)近似為38 W·m-1·K-1.計(jì)算得到上述三個(gè)溫度的熱導(dǎo)率分別為46、40、37 W·m-1·K-1,與文獻(xiàn)34結(jié)果較為接近.文獻(xiàn)35中通過(guò)理論分析BNNT的熱導(dǎo)率隨溫度變化規(guī)律,發(fā)現(xiàn)幾種直徑(3,3)(5,0)和幾種軸向長(zhǎng)度(1,5,10 μm)的BNNT的熱導(dǎo)率發(fā)生轉(zhuǎn)折時(shí)系統(tǒng)的溫度均在70-100 K范圍內(nèi),而本文熱導(dǎo)率發(fā)生轉(zhuǎn)折的溫度在70-90 K之間,這可從一個(gè)角度印證本文結(jié)果,可證明本文方法的正確性.

本文通過(guò)理論求解Boltzmann-Peierls聲子傳輸方程得出了在低溫段BNNT的導(dǎo)熱性能與外界溫度的關(guān)系.

BTE聲子傳輸方程14可以表示為

左項(xiàng)描述聲子漂移部分,右項(xiàng)描述聲子散射部分,兩項(xiàng)可分別表示為36,37

其中N為聲子分布函數(shù),N0為平衡態(tài)聲子分布函數(shù), n可表示為n=N-N0,它用來(lái)描述實(shí)際聲子分布偏離平衡態(tài)的部分,τtot為總體弛豫時(shí)間,所以(7)式可以轉(zhuǎn)化為

結(jié)合公式(10),熱流可表示為

其中s為聲子模式,q為波數(shù),ω為頻率.考慮熱流和溫度梯度的關(guān)系38

熱導(dǎo)率可表示為

其中V為體積.在低溫段,影響熱導(dǎo)率的主要因素是聲子與邊界的散射,近似認(rèn)為聲子與邊界的散射時(shí)間是散射時(shí)間τtot的主要部分.邊界散射時(shí)間可表示為

其中p為描述邊界粗糙程度的參數(shù),Λ為描述材料橫截面積的參數(shù).

本文主要考慮邊界散射對(duì)熱導(dǎo)率的影響,可近似認(rèn)為在低溫狀態(tài)下,τtot≈τb.并將公式(13)轉(zhuǎn)化為相對(duì)于波數(shù)q的積分形式

所以在低溫段納米管的熱導(dǎo)率與溫度符合κ∝T關(guān)系.

3.3 尺寸效應(yīng)

根據(jù)聲子散射理論,模型軸向長(zhǎng)度的變化將因?yàn)槠渲饾u接近或者遠(yuǎn)離聲子平均自由程而對(duì)聲子在熱輸運(yùn)過(guò)程中的散射產(chǎn)生影響,進(jìn)而影響模型的熱輸運(yùn)性能.本文建立了長(zhǎng)度為10-120 nm的10_ 10扶手椅型BNNT,模擬溫度選擇為300 K,獲得結(jié)果如圖9所示.

如圖9可知,BNNT的軸向熱導(dǎo)率隨著軸向長(zhǎng)度的增加而增加,且增加的速率逐漸減小.

圖9 BNNT軸向熱輸運(yùn)性能與長(zhǎng)度的關(guān)系Fig.9 Relationship between axial thermal conductivity and length

根據(jù)聲子散射理論,隨著B(niǎo)NNT軸向長(zhǎng)度的增加,聲子在傳輸過(guò)程中與邊界的碰撞概率逐漸減小,這將增加聲子平均自由程,根據(jù)公式(6),BNNT的熱輸運(yùn)性能逐漸增加;隨著B(niǎo)NNT軸向長(zhǎng)度逐漸接近聲子平均自由程,聲子間相互碰撞概率逐漸增大,這對(duì)隨著長(zhǎng)度增加而逐漸減小的聲子與邊界的碰撞概率起到一個(gè)抵消作用,根據(jù)上文分析,BNNT的熱輸運(yùn)性能增加趨勢(shì)逐漸減緩.Mingo和Broido14通過(guò)計(jì)算CNT的BTE得到熱導(dǎo)率與長(zhǎng)度呈κ∝Lα這一關(guān)系,且隨著L的增加,α值(描述模型彈道特性的參數(shù))逐漸減小.本文通過(guò)MD獲得的結(jié)果如擬合為κ∝Lα形式,其α值大致介于0.1989-0.4899之間,且隨著長(zhǎng)度的增加,α值逐漸減小,這與Mingo的理論結(jié)果相吻合.同時(shí)這也意味著隨著B(niǎo)NNT長(zhǎng)度的增加,模型的彈道性能逐漸減弱,但其仍主要體現(xiàn)彈道特征.

4 結(jié)論

運(yùn)用NEMD研究了BNNT的熱輸運(yùn)性能,分別研究了長(zhǎng)度、溫度和應(yīng)變對(duì)BNNT的熱輸運(yùn)性能的影響,并得到了以下結(jié)論:在拉伸應(yīng)變作用下,隨著應(yīng)變的增加,受聲子模式變化的影響,BNNT的熱導(dǎo)率呈逐漸降低的趨勢(shì);在壓縮應(yīng)變作用下,隨著應(yīng)變的增加,受BNNT發(fā)生顯著塑性變形的影響,其熱導(dǎo)率逐漸減小.當(dāng)外界溫度低于70 K時(shí),受量子效應(yīng)的影響,BNNT的熱導(dǎo)率隨著溫度的增加呈線性增加趨勢(shì);當(dāng)溫度高于100 K時(shí),其熱導(dǎo)率隨著溫度的增加而逐漸減小.當(dāng)BNNT軸向長(zhǎng)度小于120 nm時(shí),隨著長(zhǎng)度的增加,其彈道性能逐漸減弱,但仍主要體現(xiàn)為彈道特征,其熱導(dǎo)率與長(zhǎng)度基本滿足κ∝Lα這一關(guān)系.

(1) Iijima,S.Nature 1991,354,56.

(2) Rubio,A.;Corkill,J.L.;Cohen,M.L.Phys.Rev.B 1994,49 (7),5081.

(3) Loiseau,A.;Willaime,F.;Demoncy,N.;Hug,G.;Pascard,H. Phys.Rev.Lett.1996,76(25),4737.

(4) Wirtz,L.;Rubio,A.Phys.Rev.B 2003,68(4),045425.

(5) Zhang,G.;Li,B.W.The Journal of Physical Chemistry B 2005, 109(50),23823.

(6) Liew,K.M.;Yuan,J.H.Nanotechnology 2011,22,085701.

(7) Xiao,Y.;Yan,H.;Cao,J.X.;Ding,J.W.;Mao,Y.L.;Xiang,J. Phys.Rev.B 2004,69(20),205415.

(8) Chang,C.W.;Fennimore,A.M.;Afanasiev,A.;Okawa,D.; Ikuno,T.;Garcia,H.;Li,D.Y.;Majumdar,A.;Zettl,A.Phys. Rev.Lett.2006,97,085901.

(9) Chang,C.W.;Okawa,D.;Majumdar,A.;Zettl,A.Science 2006,314,1112.

(10) Zhi,C.Y.;Bando,Y.;Tang,C.C.;Golberg,D.Materials Science and Engineering R 2010,70,92.

(11) Barman,S.Europhysics Letters 2011,96,16004.

(12) Stewart,D.A.;Savic,N.;Mingo,N.Nano Lett.2009,9,81.

(13) Cho,H.B.;Tokoi,Y.;Tanaka,S.;Suematsu,H.;Suzuki,T.; Jiang,W.H.;Niihara,K.;Nakayama,T.Composites Science and Technology 2011,71,1046.

(14) Mingo,N.;Broido,D.A.Nano Lett.2005,5(7),1221.

(15) Chang,C.W.;Han,W.Q.;Zettl,A.Journal of Vacuum Science &Technology B 2005,23,1883.

(16) Tang,C.C.;Bando,Y.;Liu,C.H.;Fan,S.S.;Zhang,J.;Ding, X.X.;Golberg,D.The Journal of Physical Chemistry B 2006, 110(21),10354.

(17) Tersoff,J.Phys.Rev.B 1988,37(12),6991.

(18)Wang,X.W.;Huang,Z.;Wang,T.;Tang,Y.W.;Zeng,X.C. Physica B 2008,403,2021.

(19) Sandia National Laboratories.LAMMPS,Lammps-12Oct10; GNU:USA,2010.

(20) Schelling,P.K.;Phillpot,S.R.;Keblinski,P.Phys.Rev.B 2002, 65(14),144306.

(21) Lindsay,L.;Broido,D.A.;Mingo,N.Phys.Rev.B 2009,80, 125407.

(22) Wang,S.C.;Liang,X.G.;Xu,X.H.;Ohara,T.J.Appl.Phys. 2009,105,014316.

(23)Turney,J.E.;McGaughey,A.J.H.;Amon,C.H.Phys.Rev.B 2009,79(22),224305.

(24) Munoz,E.;Lu,J.X.;Yakobson,B.I.Nano Lett.2010,10(5), 1652.

(25) Hu,J.N.;Ruan,X.L.;Chen,Y.P.Nano Lett.2009,9(7),2730.

(26) Sese L.M.Mol.Phys.1993,78,1167.

(27) Tchouar,N.;Ould,K.F.;Levesque,D.J.Chem.Phys.2004, 121,7326.

(28) Goharshadi,E.K.;Abbaspour,M.;Kashani,H.;Baherololoom, M.Theor.Chem.Acc.2008,119,355.

(29) Yang,N.;Zhang,G.;Li,B.W.Nano Lett.2008,8(1),276.

(30) Kang,D.D.;Hou,Y.;Dai,J.Y.;Yuan,J.M.Phys.Rev.A 2009, 79,063202.

(31) Huang,M.J.;Weng,C.C.;Chang,T.M.International Journal of Thermal Sciences 2009,49,1095.

(32) Ren,C.L.;Zhang,W.;Xu,Z.J.;Zhu,Z.Y.;Huai,P. The Journal of Physical Chemistry C 2009,114(13),5786.

(33) Xu,Z.P.;Buehler,M.J.Nanotechnology 2009,20,185701.

(34) Shen,H.J.Computational Materials Science 2009,47,220.

(35) Jiang,J.W.;Wang,J.S.Phys.Rev.B 2011,84(8),085439.

(36) Goicochea,J.V.;Madrid,M.;Amon,C.ASME J.Heat Transfer 2010,132,012401.

(37)Chen,Y.F.;Li,D.Y.;Lukes,J.R.;Majumdar,A.ASME J.Heat Transfer 2005,127,1129.

(38) Nika,D.L.;Pokatilov,E.P.;Askerov,A.S.;Balandin,A.A. Phys.Rev.B 2009,79(15),155413.

November 30,2011;Revised:January 28,2012;Published on Web:February 27,2012.

Molecular Dynamics Simulation of Thermal Transport Properties for Boron Nitride Nanotubes

GAO Yu-Fei1MENG Qing-Yuan1ZHANG Lu2LIU Jia-Qiu3JING Yu-Hang1,*
(1School of Astronautics,Harbin Institute of Technology,Harbin 150001,P.R.China;2China Medical University, Shenyang 110001,Liaoning Province,P.R.China;3Harbin FRP Research Institute,Harbin 150001,P.R.China)

The Boltzmann-Peierls phonon transport equation(BTE)and non-equilibrium molecular dynamics simulation(NEMD)are used to investigate the thermal transport properties of boron nitride nanotubes(BNNTs).First,the thermal-mechanical coupling is explored using NEMD.Then,by combining BTE and NEMD,the influence of temperature and length is investigated.Quantum correction is used to extend the range over which NEMD can be used.The results demonstrate that under low-strain conditions,the thermal conductivity decreases with increasing tensile or compressive strain.Then the phonon density of state(PDOS)is used to analyze the trends in thermal transport properties theoretically; it is found that the variations in thermal transport properties under tension are caused by changes in the phonon modes,and that under compression changes are induced by the flection of the BNNT structure. The BNNT thermal conductivity increases linearly with increasing temperature because of the quantum effect at low temperatures,and it decreases significantly as the temperature reaches a certain value. When the BNNT length is less than 120 nm,the BNNT?s ballistic characteristics weaken with increasing length,but it also performs ballistic characteristic mainly,and thermal conductivity(κ)and length(L)obey the relationshipκ∝Lα.

Boron nitride nanotubes;Thermal transport property;Phonon Boltzmann-Peierls transport equation;Non-equilibrium molecular dynamics simulation;Phonon density of state;Quantum correction

10.3866/PKU.WHXB201202273

?Corresponding author.Email:jingyuhang@gmail.com;Tel:+86-451-86414143.

The project was supported by the National Natural Science Foundation of China(10772062).

國(guó)家自然科學(xué)基金(10772062)資助項(xiàng)目

O641;O642