魯濟(jì)豹 楊楠楠,2 孫 蓉 汪正平,3

1(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院 廣東省高密度電子封裝關(guān)鍵材料重點(diǎn)實(shí)驗(yàn)室 深圳 518055)

2(中國(guó)科學(xué)技術(shù)大學(xué)納米科學(xué)技術(shù)學(xué)院 蘇州 215123)

3(香港中文大學(xué)電子工程系 香港 999077)

1 引 言

隨著電子元器件的微型化和智能化發(fā)展，高密度的電子器件使得集成電路的散熱問(wèn)題成為制約其穩(wěn)定性、可靠性和壽命長(zhǎng)短的關(guān)鍵因素。為使電子器件發(fā)揮最佳性能，必須開(kāi)發(fā)具有優(yōu)良導(dǎo)熱性能及機(jī)械性能的熱界面材料，以確保其產(chǎn)生的熱量及時(shí)排出。廣泛應(yīng)用的聚合物基熱界面復(fù)合材料采用聚合物為基底，通過(guò)添加高導(dǎo)熱率的填料以增加其熱導(dǎo)率。在各種備選填料材料中，氮化硼(BN)納米材料(納米片、納米管)由于具有高耐熱性、高導(dǎo)熱性以及優(yōu)異的高溫絕緣性(滿足了電子封裝對(duì)材料介電性能的要求)，從而在電子封裝熱傳導(dǎo)材料中具有重要的應(yīng)用價(jià)值[1-8]。

復(fù)合材料的導(dǎo)熱率與多種因素相關(guān)。如以氮化硼納米材料為填料的聚合物基復(fù)合材料，其導(dǎo)熱率與填料的尺寸、分散程度、是否形成三維導(dǎo)熱網(wǎng)絡(luò)以及填料與聚合物之間的結(jié)合強(qiáng)度等密切相關(guān)[8-10]。實(shí)驗(yàn)中對(duì)于氮化硼-聚合物復(fù)合材料的熱性能研究主要通過(guò)表征復(fù)合材料宏觀熱導(dǎo)率來(lái)實(shí)現(xiàn)。但是，由于復(fù)合材料導(dǎo)熱性能與材料內(nèi)部微觀結(jié)構(gòu)密切相關(guān)，所以想要了解復(fù)合材料導(dǎo)熱機(jī)理，并對(duì)材料進(jìn)行深層次的結(jié)構(gòu)優(yōu)化設(shè)計(jì)，需要對(duì)界面熱阻等熱相關(guān)量的精準(zhǔn)測(cè)量以及理論模擬計(jì)算。實(shí)驗(yàn)中進(jìn)行界面熱阻納米尺度上的表征對(duì)實(shí)驗(yàn)設(shè)備的要求具有非常大的挑戰(zhàn)。另外，實(shí)驗(yàn)測(cè)量結(jié)果強(qiáng)烈依賴于樣品質(zhì)量、樣品尺寸、測(cè)量方法和測(cè)量溫度等多種因素，這使得不同的實(shí)驗(yàn)結(jié)果難以有效對(duì)比。而分子動(dòng)力學(xué)模擬由于直接對(duì)物質(zhì)微觀結(jié)構(gòu)進(jìn)行表征，成為研究熱界面復(fù)合材料導(dǎo)熱機(jī)理及微觀結(jié)構(gòu)設(shè)計(jì)的重要工具，并對(duì)實(shí)驗(yàn)測(cè)量形成重要補(bǔ)充。

采用分子動(dòng)力學(xué)對(duì)熱導(dǎo)率以及聲子性質(zhì)的計(jì)算通常有兩種方法——非平衡態(tài)分子動(dòng)力學(xué)(Non-Equilibrium Molecular Dynamics，NEMD)和平衡態(tài)分子動(dòng)力學(xué)(Equilibrium Molecular Dynamics，EMD)方法[11-16]。相對(duì)于 NEMD 方法來(lái)說(shuō)，EMD 方法可以采用一個(gè)平衡態(tài)模擬計(jì)算得到完整的熱導(dǎo)率張量，從而得到材料在各個(gè)方向上的熱導(dǎo)率。并且，EMD 方法中的周期性邊界條件允許聲子在邊界處穿梭而不受散射，從而計(jì)算熱導(dǎo)率受體系的邊界條件以及尺寸的影響要小很多。但由于計(jì)算資源限制，分子動(dòng)力學(xué)模擬所構(gòu)建的模型尺寸一般在納米量級(jí)，通常遠(yuǎn)小于實(shí)驗(yàn)中測(cè)量樣品的實(shí)際尺寸。所以對(duì)于 EMD 方法，仍然要求在計(jì)算模擬中對(duì)樣品尺寸與性質(zhì)之間的關(guān)系具有清晰的理解。Wang 和 Ruan[11]針對(duì)單晶硅、石墨烯以及硅烯，詳細(xì)研究了模擬尺寸對(duì)采用分子動(dòng)力學(xué)計(jì)算其熱相關(guān)性質(zhì)的影響，發(fā)現(xiàn)樣品尺寸對(duì) 3 種材料的熱導(dǎo)率計(jì)算具有不同的影響，并且這種影響隨所采用計(jì)算方法的(EMD 或 NEMD)不同而不同。對(duì)于氮化硼納米結(jié)構(gòu)，雖然文獻(xiàn)中有一些針對(duì)聲子色散及熱導(dǎo)率的模擬計(jì)算，但尚未有針對(duì)熱學(xué)相關(guān)性質(zhì)與樣品尺寸之間依賴關(guān)系及規(guī)律的明確探討。因此，在此工作中，采用 EMD 結(jié)合 Green-Kubo方法，對(duì)氮化硼單層結(jié)構(gòu)的熱導(dǎo)率、聲子色散曲線和聲子態(tài)密度(Density of States，DOS)與模擬樣品尺寸之間的關(guān)系及影響機(jī)理進(jìn)行了深入探討，為后續(xù)對(duì)氮化硼材料以及氮化硼-聚合物復(fù)合材料的熱相關(guān)性質(zhì)的計(jì)算提供重要參考。

2 模型、方法與計(jì)算細(xì)節(jié)

2.1 模型：Tersoff 力場(chǎng)

分子動(dòng)力學(xué)對(duì)材料性質(zhì)的表征精度依賴于所采用的力場(chǎng)是否能夠很好地再現(xiàn)材料的相關(guān)實(shí)驗(yàn)結(jié)果或者第一性原理計(jì)算結(jié)果。Sevik 等[17]基于六方氮化硼的結(jié)構(gòu)、力學(xué)和振動(dòng)性質(zhì)，優(yōu)化構(gòu)建了氮化硼結(jié)構(gòu)的 Tersoff 力場(chǎng)，該力場(chǎng)很好地再現(xiàn)了實(shí)驗(yàn)中氮化硼的聲子色散關(guān)系。隨后，Tersoff 力場(chǎng)被廣泛用于氮化硼結(jié)構(gòu)的分子動(dòng)力學(xué)模擬中。本文采用 Tersoff 力場(chǎng)描述氮化硼二維單層材料中氮(N)與硼(B)的相互作用[17]。在Tersoff 力場(chǎng)中，原子之間的相互作用采用勢(shì)函數(shù)Vij描述。

其中，rij為原子i與j之間的距離；和分別為吸引和排斥相互作用；為截?cái)嗑嚯x(cutoff)，截?cái)嗑嚯x保證了只有最緊鄰原子之間存在相互作用；aij是對(duì)排斥相互作用勢(shì)的限制項(xiàng)，在此aij＝1；A、B表示量綱為[能量]的系數(shù)；λ表示量綱為[1/距離]的系數(shù)；e為自然對(duì)數(shù)的底數(shù)；鍵長(zhǎng)和鍵角項(xiàng)bij取決于原子i周圍原子的配位數(shù)以及原子i、j和k之間的夾角。

2.2 計(jì)算方法

2.2.1 采用平衡態(tài)分子動(dòng)力學(xué)與 Green-Kubo 方法計(jì)算熱導(dǎo)率

相對(duì)于 NEMD 模擬，EMD 模擬計(jì)算熱導(dǎo)率的方法受體系的邊界條件以及尺寸的影響要小很多。在此工作中，采用 Green-Kubo 方法[18]結(jié)合EMD，計(jì)算氮化硼單層結(jié)構(gòu)的熱導(dǎo)率

其中，kB、T分別為玻爾茲曼常數(shù)和溫度；V為體系體積；τm為時(shí)間；為熱流密度在時(shí)間t時(shí)的α項(xiàng)；分別為沿α和β方向的熱導(dǎo)率。氮化硼單層的面內(nèi)熱導(dǎo)率由和二者平均值得到。

2.2.2 聲子色散關(guān)系及聲子態(tài)密度的計(jì)算

氮化硼二維單層結(jié)構(gòu)的聲子色散關(guān)系以及聲子態(tài)密度(DOS)可以由體系的動(dòng)態(tài)矩陣得出。在分子動(dòng)力學(xué)中，動(dòng)態(tài)矩陣可以采用漲落耗散理論計(jì)算得到。假設(shè)晶體中包含N個(gè)原包，標(biāo)記為li＝(l1,l2,l3)，其中l(wèi)i為整數(shù)。每個(gè)原包定義在由 3 個(gè)線性獨(dú)立的基矢a1、a2、a3組成的平行六面體中，并且每個(gè)原包包含的基元原子數(shù)目為K，基元原子標(biāo)記為k?；跐q落耗散理論，晶體體系在倒格空間中的力常數(shù)系數(shù)為[15,19]：

其中，G為格林函數(shù)系數(shù)[13,19]。

其中，R是原子的瞬時(shí)位置；代表原子位置的平均值。這種方法可以給出與位移方法同樣的結(jié)果，并且便于在分子動(dòng)力學(xué)代碼中實(shí)現(xiàn)。

一旦得到了力常數(shù)矩陣，則動(dòng)態(tài)矩陣D可由如下公式得出：此矩陣的本征值即為q處的聲子頻率。所以，通過(guò)求解D的本征值，可以得到所有聲子模式的頻率。

其中，v為聲子頻率；e為極化向量。由此，q和v之間的關(guān)系得以求出。聲子態(tài)密度可由聲子色散關(guān)系除以頻率范圍以得到多個(gè)小片段，然后計(jì)算每個(gè)片段內(nèi)的態(tài)數(shù)而得到。

2.2.3 分子動(dòng)力學(xué)模擬細(xì)節(jié)

分子動(dòng)力學(xué)模擬采用 LAMMPS[20]軟件包進(jìn)行。如前所述，N 和 B 之間的相互作用采用Tersoff-2010 力場(chǎng)進(jìn)行描述。所有計(jì)算均采用x、y、z方向的周期性邊界條件。為了研究模擬尺寸對(duì)氮化硼單層結(jié)構(gòu)熱導(dǎo)率、聲子色散關(guān)系以及聲子態(tài)密度計(jì)算的影響，構(gòu)建了不同尺寸的氮化硼單層結(jié)構(gòu)，分別為 10×10、18×18、35×35、52×52、70×70。其中，10×10 代表在x、y方向上分別包含 10 個(gè)氮化硼原包，其他依此類推。這些結(jié)構(gòu)分別對(duì)應(yīng)的面積為 5.24 nm2、16.99 nm2、64.20 nm2、141.71 nm2、256.80 nm2。圖 1 為18×18(16.99 nm2)氮化硼單層結(jié)構(gòu)模型。分子動(dòng)力學(xué)步長(zhǎng)為 0.5 fs[11,12,16]。對(duì)于每個(gè)尺寸的結(jié)構(gòu)，首先采用 Nose-Hoover 恒溫恒壓方法[21]，在NPT系綜(體系原子數(shù)目、壓強(qiáng)、溫度恒定)中，在 0 bar 壓強(qiáng)、300 K 溫度下平衡 1 ns，以確保體系完全達(dá)到平衡態(tài)。

圖1 EMD 模擬中構(gòu)建的尺寸為 16.99 nm2 的氮化硼單層結(jié)構(gòu)Fig. 1 The sturcture of the boron nitride monolayer with the size of 16.99 nm2 built in EMD simulations

對(duì)于動(dòng)態(tài)矩陣的計(jì)算，體系在NPT中平衡后，采用 LAMMPS[20]中的 fix phonon 命令在NVE系綜(體系原子數(shù)目、體積以及總能量保持恒定)中進(jìn)行，并采用輔助代碼 phana 進(jìn)行數(shù)據(jù)后處理得到聲子色散譜和態(tài)密度[13,15,22]。對(duì)于熱導(dǎo)率的計(jì)算，體系在NPT中平衡后，采用NVE系綜運(yùn)行 6 ns 以進(jìn)行數(shù)據(jù)收集。其中，原子位置軌跡每隔 10 個(gè)步長(zhǎng)記錄一次。采用 Green-Kubo方法計(jì)算熱導(dǎo)率，其中對(duì)熱流密度自相關(guān)函數(shù)積分的時(shí)間上限設(shè)置為 400 ps。我們分別測(cè)試了更長(zhǎng)時(shí)間上限(如 500 ps、1 ns)，證實(shí) 400 ps 的時(shí)間積分上限可以得到收斂的熱導(dǎo)率值。為了計(jì)算氮化硼單層的體系，氮化硼單層的厚度選取為0.335 nm[11]。熱導(dǎo)率的計(jì)算沒(méi)有包含量子修正。為了得到更為精確的熱導(dǎo)率值及其誤差棒，對(duì)每個(gè)尺寸的樣品，分別進(jìn)行 5 次獨(dú)立的計(jì)算，其中每次計(jì)算具有獨(dú)立的原子初始位置及初始速度。熱導(dǎo)率的最終值以及誤差棒由 5 次計(jì)算的平均值和標(biāo)準(zhǔn)偏差得出。值得注意的是，氮化硼單層結(jié)構(gòu)的熱導(dǎo)率具有輕微的各向異性，即在 armchair和 zigzag 方向熱導(dǎo)率的值略有不同。本文報(bào)道的熱導(dǎo)率值為兩個(gè)方向熱導(dǎo)率的平均值。

值得說(shuō)明的是，文中所涉及的計(jì)算均在國(guó)家超級(jí)計(jì)算廣州中心完成。

3 結(jié)果和討論

3.1 熱流密度自相關(guān)函數(shù)

在 EMD 計(jì)算熱導(dǎo)率的方法中，熱導(dǎo)率由體系中熱流密度自相關(guān)函數(shù)(HCACF)對(duì)時(shí)間的積分求出(公式 7)。理論上，對(duì)時(shí)間的積分應(yīng)該至無(wú)窮大，即至 HCACF 衰減為 0 時(shí)。但在實(shí)際計(jì)算中，采用積分至?xí)r間τm代替時(shí)間無(wú)窮大，并近似認(rèn)為 HCACF 在時(shí)間τm已衰減至 0。在本工作的 EMD 模擬中，氮化硼單層模型為無(wú)缺陷結(jié)構(gòu)，并且采用了周期性邊界條件，所以不存在聲子-邊界以及聲子-缺陷散射效應(yīng)。因此，HCACF的衰減由非簡(jiǎn)諧聲子-聲子散射(Umpklapp 散射)相互作用導(dǎo)致。在本工作中，τm選取為 400 ps。圖 2 為 HCACF 隨時(shí)間的變化情況。由圖 2 可看出，在τm＝400 ps 時(shí)，HCACF 已經(jīng)近似衰減至0。值得注意的是，圖 2 所示 HCACF 隨時(shí)間表現(xiàn)出了劇烈的震蕩。這種劇烈震蕩是由于氮化硼單層結(jié)構(gòu)的高頻光學(xué)聲子所致，并且可以采用傅里葉低通濾波的方法將高頻光學(xué)聲子對(duì) HCACF的貢獻(xiàn)濾去[23]。

3.2 模擬尺寸對(duì)氮化硼單層熱導(dǎo)率的影響

有文獻(xiàn)報(bào)道，采用 EMD 方法，硅(Si)晶體熱導(dǎo)率隨尺寸的增加而增加，而石墨烯熱導(dǎo)率則表現(xiàn)出相反的性質(zhì)[11]?？紤]到 Si 晶體為三維結(jié)構(gòu)，而石墨烯為二維，那么對(duì)于二維材料來(lái)說(shuō)，由 EMD 方法計(jì)算出的熱導(dǎo)率是否總是隨尺寸的增加而減小呢？以往文獻(xiàn)尚沒(méi)有關(guān)于采用EMD 結(jié)合 Green-Kubo 方法計(jì)算氮化硼熱導(dǎo)率隨模擬尺寸變化規(guī)律的明確報(bào)道。為此，本文采用EMD 結(jié)合 Green-Kubo 方法，計(jì)算了廣泛應(yīng)用于熱傳導(dǎo)材料的氮化硼二維材料的熱導(dǎo)率。圖 3 為本文在 300 K 溫度下計(jì)算的氮化硼單層結(jié)構(gòu)熱導(dǎo)率隨模擬尺寸的變化。由圖 3 可以看出，隨著模擬尺寸從 5.241 nm2開(kāi)始增加，氮化硼單層熱導(dǎo)率迅速下降，并從面積約 17 nm2開(kāi)始，熱導(dǎo)率達(dá)到一個(gè)相對(duì)穩(wěn)定的值；繼續(xù)增加模擬尺寸不再明顯影響熱導(dǎo)率的計(jì)算。對(duì)于這種現(xiàn)象的解釋為：隨著模擬尺寸的增大，聲子數(shù)量增多，從而聲子之間散射過(guò)程增多，這一散射過(guò)程是不利于熱傳導(dǎo)的。在氮化硼單層或石墨烯中，增加的聲子之間的散射過(guò)程對(duì)熱導(dǎo)率的阻礙作用，超過(guò)了由于尺寸的增加而增加的聲子模式對(duì)于熱導(dǎo)率的貢獻(xiàn)，從而使熱導(dǎo)率總體上表現(xiàn)出隨模擬尺寸的增大而減小的現(xiàn)象。收斂后的氮化硼熱導(dǎo)率值為(349±22)W/(m·K)。其中，zigzag 方向?yàn)?285±46)W/(m·K)，armchair 方向?yàn)?413±67)W/(m·K)，與之前文獻(xiàn)報(bào)道的符合良好[23,24]。相比于 EMD 模擬中石墨烯熱導(dǎo)率收斂尺寸(10 nm×10 nm)而言[11]，氮化硼熱導(dǎo)率在約 4.1 nm×4.1 nm 即開(kāi)始達(dá)到收斂(圖 3)。這從側(cè)面表明，在氮化硼單層結(jié)構(gòu)中，聲子之間的散射作用遠(yuǎn)大于石墨烯中聲子之間的散射。

圖2 氮化硼單層結(jié)構(gòu)在 EMD 模擬中的熱流自相關(guān)函數(shù)(HCACF)隨時(shí)間的變化(內(nèi)部插圖更詳細(xì)地顯示了 HCACF 的震蕩行為)Fig. 2 The heat fl ux auto-correlation function (HCACF) as a function of time in the boron nitride monolayer structure( The inset shows more details of the strong oscillations of the HCACF)

圖3 EMD 方法計(jì)算得到的氮化硼單層結(jié)構(gòu)熱導(dǎo)率(σ)隨模擬體系面積(S)的變化Fig. 3 Thermal conductivity (σ) of the boron nitride monolayer as a function of the area (S) of the structure

3.3 模擬尺寸對(duì)氮化硼單層聲子色散關(guān)系及聲子態(tài)密度的影響

圖4 為不同模擬尺寸的氮化硼單層結(jié)構(gòu)的聲子色散沿著Γ-M-K-Γ方向的曲線。由圖 4 可以看出，模擬尺寸對(duì)氮化硼單層聲子色散關(guān)系幾乎沒(méi)有影響。進(jìn)一步對(duì)聲子態(tài)密度(DOS)的計(jì)算表明，不同模擬尺寸情況下，計(jì)算得到的氮化硼單層結(jié)構(gòu)的聲子態(tài)密度很好地重合在一起，表明了聲子態(tài)密度分布對(duì)模擬尺寸的不敏感性(圖 5)。值得注意的是，本工作中考慮的氮化硼單層結(jié)構(gòu)為無(wú)缺陷晶體結(jié)構(gòu)。對(duì)于含有缺陷的結(jié)構(gòu)，如摻雜、空位以及表面化學(xué)吸附等的氮化硼結(jié)構(gòu)，模擬尺寸對(duì)熱導(dǎo)率、聲子色散和態(tài)密度的影響可能會(huì)更加復(fù)雜，從而需要更多的工作進(jìn)行進(jìn)一步研究。

3.4 與國(guó)內(nèi)外相似研究的對(duì)比分析

Zhu 和 Ertekin[25]采用 NEMD 方法計(jì)算得到氮化硼單層結(jié)構(gòu)熱導(dǎo)率值為(600±30)W/(m·K)。而 Li 和 Chen[24]采用第一性原理分子動(dòng)力學(xué)求解線性聲子玻爾茲曼方程(LBTE)的方法，計(jì)算得到的氮化硼單層結(jié)構(gòu)熱導(dǎo)率約為 220 W/(m·K)，遠(yuǎn)小于 Zhu 和 Ertekin[25]得到的數(shù)值。Zhou 等[26]和 Jo 等[27]對(duì)氮化硼結(jié)構(gòu)的熱導(dǎo)率值進(jìn)行了實(shí)驗(yàn)測(cè)量，測(cè)得的氮化硼結(jié)構(gòu)的熱導(dǎo)率值在 227～280 W/(m·K)，與 Li 和 Chen[24]得到的數(shù)值接近。本文采用 EMD 結(jié)合 Green-Kubo 方法計(jì)算得到的氮化硼單層結(jié)構(gòu)熱導(dǎo)率為(349±22)W/(m·K)，和文獻(xiàn)報(bào)道的實(shí)驗(yàn)測(cè)量[26,27]和理論計(jì)算[24]值有比較好的吻合。本文得到的熱導(dǎo)率與實(shí)驗(yàn)測(cè)量和其他理論計(jì)算結(jié)果在定量上的差別主要是由于：(1)Tersoff 力場(chǎng)模擬氮化硼實(shí)際上采用的是擬合的經(jīng)典力學(xué)方程，而非從第一性原理出發(fā)；(2)由于受計(jì)算體系以及采樣數(shù)目的限制，模擬計(jì)算過(guò)程中不可避免地受到熱漲落影響而引入誤差。本文計(jì)算的氮化硼聲子色散曲線與文獻(xiàn)中采用 Tersoff 力場(chǎng)報(bào)道的計(jì)算結(jié)果吻合非常好[24,28]，并與 Serrano 等[29]實(shí)驗(yàn)測(cè)量的聲子散射有很好的吻合度，說(shuō)明了本文計(jì)算的可靠性。文獻(xiàn)尚無(wú)明確地針對(duì) EMD 結(jié)合 Green-Kubo 方法研究氮化硼單層結(jié)構(gòu)熱導(dǎo)率隨模擬尺寸變化規(guī)律的討論，只有關(guān)于石墨烯的相關(guān)詳細(xì)報(bào)道[11]。相比于 EMD 模擬中石墨烯熱導(dǎo)率收斂尺寸(10 nm×10 nm)而言[11]，氮化硼熱導(dǎo)率在約4.1 nm×4.1 nm 即開(kāi)始達(dá)到收斂(圖 3)。這從側(cè)面表明，在氮化硼單層結(jié)構(gòu)中，聲子之間的散射作用大于石墨烯中聲子之間的散射。

圖4 由格林函數(shù)法計(jì)算得到的不同模擬尺寸(面積)氮化硼單層結(jié)構(gòu)的聲子色散曲線Fig. 4 The phonon dispersion of the boron nitride monolayer with different sizes calculated with the Green-Kubo method

圖5 不同模擬尺寸下氮化硼單層結(jié)構(gòu)的聲子態(tài)密度Fig. 5 The phonon density of states of the boron nitride monolayer with different sizes

4 結(jié) 論

本工作采用 EMD 結(jié)合 Green-Kubo 方法，對(duì)純凈氮化硼單層結(jié)構(gòu)的熱導(dǎo)率、聲子色散、聲子態(tài)密度進(jìn)行了計(jì)算，研究了模擬尺寸對(duì)這些性質(zhì)的影響規(guī)律。研究表明，氮化硼單層結(jié)構(gòu)的聲子色散曲線以及聲子態(tài)密度并不受模擬尺寸的影響。不同的是，氮化硼單層熱導(dǎo)率隨模擬尺寸的增大而減小，至單層結(jié)構(gòu)面積達(dá)到 17 nm2時(shí)，熱導(dǎo)率開(kāi)始達(dá)到收斂值(349±22)W/(m·K)。這不同于硅晶體熱導(dǎo)率在 EMD 方法中隨著尺寸的增大而增大的現(xiàn)象，但與石墨烯、硅烯等單層結(jié)構(gòu)類似。這是由于在單層結(jié)構(gòu)中，隨模擬尺寸的增大，增加的聲子之間的散射對(duì)熱導(dǎo)率的抑制效應(yīng)超過(guò)了其對(duì)于熱導(dǎo)率的貢獻(xiàn)。氮化硼單層熱導(dǎo)率收斂尺寸 17 nm2小于石墨烯的收斂尺寸(～100 nm2)，這說(shuō)明了氮化硼單層結(jié)構(gòu)中的聲子-聲子散射要大于石墨烯，與氮化硼單層熱導(dǎo)率小于石墨烯相符。

[1]孫娜, 曾小亮, 陳鵬, 等. 高導(dǎo)熱氮化硼/聚碳酸酯復(fù)合材料的制備與性能研 [J]. 集成技術(shù), 2017,6(6): 15-23.

[2]Hu JT, Huang Y, Yao YM, et al. Polymer composite with improved thermal conductivity by constructing a hierarchically ordered three-dimensional interconnected network of BN [J]. ACS Applied Materials & Interfaces, 2017, 9(15): 13544-13553.

[3]Zhang J, Wang XN, Yu CP, et al. A facile method to prepare flexible boron nitride/poly(vinyl alcohol)composites with enhanced thermal conductivity[J]. Composites Science and Technology, 2017,149(Supplement C): 41-47.

[4]Huang T, Zeng XL, Yao XM, et al. A novel h-BNRGO hybrids for epoxy resin composites achieving enhanced high thermal conductivity and energy density [J]. RSC Advances, 2017, 7(38): 23355-23362.

[5]Wu K, Fang JC, Ma JR, et al. Achieving a collapsible, strong, and highly thermally conductive fi lm based on oriented functionalized boron nitride nanosheets and cellulose nano fi ber [J]. ACS Applied Materials & Interfaces, 2017, 9(35): 30035-30045.

[6]Li LH, Chen YL. Atomically thin boron nitride:unique properties and applications [J]. Advanced Functional Materials, 2016, 26(16): 2594-2608.

[7]Chen J, Huang XY, Zhu YK, et al. Cellulose nanofiber supported 3D interconnected BN nanosheets for epoxy nanocomposites with ultrahigh thermal management capability [J].Advanced Functional Materials, 2017, 27(5):1604754.

[8]Sun JJ, Yao YM, Zeng XL, et al. Preparation of boron nitride nanosheet/nanofibrillated cellulose nanocomposites with ultrahigh thermal conductivity via eengineering interfacial thermal resistance[J]. Advanced Materials Interfaces, 2017, 4(17):1700563.

[9]Zeng XL, Sun JJ, Yao YM, et al. A combination of boron nitride nanotubes and cellulose nanofibers for the preparation of a nanocomposite with high thermal conductivity [J]. ACS Nano, 2017, 11(5):5167-5178.

[10]Yao YM, Zeng XL, Wang FF, et al. Significant enhancement of thermalconductivity in bioinspired freestanding boron nitride papersfilled with graphene oxide [J]. Chemistry of Materials, 2016,28(4): 1049-1057.

[11]Wang ZY, Ruan XL. On the domain size effect of thermal conductivities from equilibrium and nonequilibrium molecular dynamics simulations[J]. Journal of Applied Physics, 2017, 121(4):044301.

[12]Singh SK, Neek-Amal M, Costamagna S, et al.Thermomechanical properties of a single hexagonal boron nitride sheet [J]. Physical Review B, 2013,87(18): 184106.

[13]Kong LT. Phonon dispersion measured directly from molecular dynamics simulations [J]. Computer Physics Communications, 2011, 182(10): 2201-2207.

[14]Lindsay L, Broido DA. Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene [J]. Physical Review B,2010, 81(20): 262-265.

[15]Kong LT, Bartels G, Campa?á C, et al.Implementation of green’s function molecular dynamics: an extension to LAMMPS [J]. Computer Physics Communications, 2009, 180(6): 1004-1010.

[16]Tabarraei A,Wang XN. Anomalous thermal conductivity of monolayer boron nitride [J].Applied Physics Letters, 2016, 108(18): 181904.

[17]Sevik C, Kinaci A, Haskins JB, et al. Characterization of thermal transport in low-dimensional boron nitride nanostructures [J]. Physical Review B,2011, 84(8): 085409.

[18]Kubo R. Statistical-mechanical theory of irreversible processes. I. general theory and simple applications to magnetic and conduction problems[J]. Journal of the Physical Society of Japan, 1957,12(6): 570-586.

[19]Campa?á C, Müser MH. Practical Green’s function approach to the simulation of elastic semi-infinite solids [J]. Physical Review B, 2006, 74(7): 075420.

[20]Plimpton S. Fast parallel algorithms for short-range molecular dynamics [J]. Journal of Computational Physics, 1995, 117(1): 1-19.

[21]Nosé S. A unified formulation of the constant temperature molecular dynamics methods [J]. The Journal of Chemical Physics, 1984, 81(1): 511-519.

[22]Kong LT, Denniston C, Müser MT. An improved version of the Green’s function molecular dynamics method [J]. Computer Physics Communications,2011, 182(2): 540-541.

[23]Mahdizadeh SJ, Goharshadi EK, Akhlamadi G.Thermo-mechanical properties of boron nitride nanoribbons: a molecular dynamics simulation study [J]. Journal of Molecular Graphics and Modelling, 2016, 68(Supplement C): 1-13.

[24]Li S, Chen Y. Thermal transport and anharmonic phonons in strained monolayer hexagonal boron nitride [J]. Scienti fi c Reports, 2017, 7: 43956.

[25]Zhu TS, Ertekin E. Resolving anomalous strain effects on two-dimensional phonon fl ows: the cases of graphene, boron nitride, and planar superlattices[J]. Physical Review B, 2015, 91(20): 205429.

[26]Zhou HQ, Zhu JX, Liu Z, et al. High thermal conductivity of suspended few-layer hexagonal boron nitride sheets [J]. Nano Research, 2014, 7(8):1232-1240.

[27]Jo I, Pettes MT, Kim J, et al. Thermal conductivity and phonon transport in suspended few-layer hexagonal boron nitride [J]. Nano Letters, 2013,13(2): 550-554.

[28]Anees P, Valsakumar MC, Panigrahi BK. Effect of strong phonon-phonon coupling on the temperature dependent structural stability and frequency shift of 2D hexagonal boron nitride [J]. Physical Chemistry Chemical Physics, 2016, 18(4): 2672-2681.

[29]Serrano J, Bosak A, Arenal R, et al. Vibrational properties of hexagonal boron nitride: inelastic X-ray scattering and ab initio calculations [J].Physical Review Letters, 2007, 98(9): 095503.