魯濟豹 楊楠楠,2 孫 蓉 汪正平,3
1(中國科學(xué)院深圳先進技術(shù)研究院 廣東省高密度電子封裝關(guān)鍵材料重點實驗室 深圳 518055)
2(中國科學(xué)技術(shù)大學(xué)納米科學(xué)技術(shù)學(xué)院 蘇州 215123)
3(香港中文大學(xué)電子工程系 香港 999077)
隨著電子元器件的微型化和智能化發(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ì)的計算提供重要參考。
分子動力學(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之間的夾角。
相對于 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)率由和二者平均值得到。
氮化硼二維單層結(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ù)而得到。
分子動力學(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)率的平均值。
值得說明的是,文中所涉及的計算均在國家超級計算廣州中心完成。
在 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]。
有文獻報道,采用 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
圖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ù)雜,從而需要更多的工作進行進一步研究。
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
本工作采用 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.