魯濟豹 楊楠楠,2 孫 蓉 汪正平,3

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

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

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

1 引 言

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

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

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

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

2.1 模型：Tersoff 力場

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

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

2.2 計算方法

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

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

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

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

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

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

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

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

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

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

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

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

值得說明的是，文中所涉及的計算均在國家超級計算廣州中心完成。

3 結(jié)果和討論

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

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

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

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

圖2 氮化硼單層結(jié)構(gòu)在 EMD 模擬中的熱流自相關(guān)函數(shù)(HCACF)隨時間的變化(內(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 方法計算得到的氮化硼單層結(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 模擬尺寸對氮化硼單層聲子色散關(guān)系及聲子態(tài)密度的影響

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

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

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

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