石墨炔納米帶熱傳導的分子動力學模擬

，，

(福州大學 電氣工程與自動化學院，福建　福州　350108)

1　引言

電子器件隨著制造工藝的不斷發展，微型化達到納米量級，但是急劇減小的尺寸引起器件內熱量快速增加，于是納米尺度下的器件能否穩定運行取決于其產生的熱量轉移效率的高低。只有熱量得到及時的排出，器件內部的局部溫度才不會超過其正常工作溫度，有利于延長器件的使用壽命，否則散熱問題將嚴重制約器件微型化的進一步發展。

Baughman等[1]預測了一種新型的碳結構-石墨炔。Li等[2]利用六炔基苯在銅片催化的作用下發生偶聯反應，在銅片表面上成功制備出大塊 graphdiyne 薄膜，引起科學界對 graphyne 的關注熱潮。Zhang等[3]通過分子動力學模擬方法研究了graphyne的導熱系數，結果表明graphyne存在乙炔鏈使得其導熱系數低于石墨烯的。歐陽等[4]通過格林方法對graphyne納米帶的熱輸運進行研究，結果發現graphyne納米帶的熱導率與其手性有關。對于小尺寸的石墨炔納米帶的尺寸效應及Airebo勢函數與Tersoff勢函數對熱導率的影響文獻比較少看到。

因為納米尺度下的石墨炔結構的變化會造成聲子在傳輸過程中的聲子移動速度和平均自由程發生變化，將影響聲子的傳熱效率[5]。為了更好地預計和控制石墨炔納米帶熱導率，本文主要研究的是長度和寬度及勢函數對石墨炔納米帶熱導率的影響。

2　分子動力學模擬

2.1　研究方法

分子動力學模擬方法是利用計算機進行模擬實驗，來研究和反應整個過程，這對于理論計算和實驗都具有重要意義。該法不僅可以“監控”原子，獲得原子整個過程的運動軌跡，也可以像作實驗一樣，對有興趣的量進行一系列的觀察。該法多用于多粒子體系中，其核心思想是將體系中的粒子當作經典粒子，使得粒子的運動過程遵循經典牛頓運動定律。本文通過分子動力學方法，利用經典牛頓力學可得到體系中所有粒子的運動軌跡，從而得到該體系隨時間變化的規律，再通過統計力學理論計算出體系的一些參數及熱輸運性質，然后通過這些數據計算得到溫度梯度，再帶入熱導率公式求出導熱系數，如圖1所示。

圖1　分子動力學基本流程圖

2.2　模型構建

石墨炔是由 1,3-二炔鍵將苯環共軛連接形成的二維平面網絡結構，它具有豐富的碳化學鍵，大的共軛體系、寬面間距、優良的化學穩定性以及半導體性能如圖2所示。利用MS軟件構建石墨炔納米帶模型，根據不同的切割方向可以切得扶手椅型(armchair)和鋸齒型(zigzag)石墨炔納米帶，如圖3所示。

2.3　勢函數的選取

在分子動力學模擬中，原子與原子之間的相互作用力由勢函數來體現。對于分子動力學模擬來說，勢能函數的選取決定著模擬計算的準度和精度。對于以碳為成分的物質和材料諸如石墨、金剛石、碳納米管、石墨炔、石墨烯等，采用的勢函數大多是Airebo勢函數[6]。Airebo勢函數是在Rebo勢函數的基礎上改進和推廣的，通過在Rebo勢函數中加入原子與原子之間的長程作用能和原子與原子之間及分子與原子之間的扭曲作用能，同時Airebo勢是一個依賴于周圍原子與原子之間的成鍵強度及鍵與鍵之間的耦合作用強度。Airebo勢函數的優勢在于它具有更高的精度，能精確的描述碳碳、碳氫之間的相互作用。因此，在本篇論文中，研究石墨炔納米帶的熱導率，C-C原子之間的相互作用采用Airebo勢函數來描述。

圖2　石墨炔原胞

圖3　石墨炔納米帶結構

在Rebo勢函數的基礎上，增加并加以改進了長程相互作用項和扭曲項，得到了我們所需要的Airebo勢函數，它是一個依賴于周圍原子配置的鍵合強度的勢函數。其表達式為：

(1)

(2)

(3)

其中

VR(r)=fc(r)(1+Q/rAe-αr)

(4)

(5)

2.4　非平衡分子動力學方法

目前常見的關于用分子動力學方法模擬材料的導熱性能有以下幾種方法,分別是Green-Kubo 平衡態分子動力學方法(EMD)[7-8]、正向非平衡態分子動力學方法(Non-Equilibrium Molecular Dynamics，NEMD)[9]和Muller-Plathe 逆向非平衡態分子動力學方法(RNEMD)[10]。由于EMD 方法難以處理單質以外的物質,所以通常在微正則或正則系綜下進行模擬。

本文采用的方法是非平衡分子動力學模擬方法(NEMD)，該法是由 Muller-Plathe[11]提出的，是通過交換粒子的速度實現外加熱流的算法，該算法所描述的是每隔N個時間步長交換任何兩個不同區域內的兩粒子間的速度，從而產生溫度梯度，利用公式得到熱導率。位于納米帶中間的是熱層，位于帶兩端的是冷層，且在軸向設置周期邊界條件，如圖4所示。通過交換熱層運動最慢的原子的速度和冷層運動最快的原子，實現熱層和冷層間的能量交換，結果是熱層更熱、冷層更冷，從而實現熱層到冷層方向的熱傳導。產生的熱傳導方向與粒子速度交換引起的能量傳遞方向相反[12]。

某一時刻第k層的溫度Tk為：

(6)

其中，N為層中的原子數，mi表示i原子的質量，vi為對應原子的速度，kB=1.3806505×10-23J×K-1為Boltzmann常量。

通過得到的動能差就能求得體系所施加的熱流值，進一步求得石墨炔納米帶體系的導熱系數為：

(7)

圖4　交換速度方法示意圖

3　模擬結果分析

在非平衡態分子動力學模擬過程中,判斷模擬結果可靠性的一個重要原則就是沿著熱流方向建立的溫度梯度應該呈線性分布。整個模擬結構的溫度曲線除了在熱域和冷域附近區域出現較大波動外，在模擬結構的中間區域均表現出了一致的線性度。熱域和冷域附近區域出現較大波動在非平衡態分子動力學模擬中總是存在的，這是由于聲子在熱域和冷域發生了強烈的散射，導致溫度曲線出現了非線性。通過模擬結構的溫度曲線獲取溫度梯度時，本文只考慮了模擬結構中間的線性區域，所以使用Fourier定律計算熱導率是可行的[14]。

為了加快計算速率和通過計算精度，使用Verlet算法求解運動方程，積分時間步長為0.4fs。整個模擬過程分為3步：首先，在等溫等壓正則(NPT)系綜(模擬過程中保持N-原子數量，P-壓力，T-溫度不變的系綜)下模擬4×105step；再轉到正則(NVT)系綜(模擬過程中保持N-原子數量、V-體積、T-溫度不變的系綜)下繼續運行4×105step；然后在微正則(NVE)系綜下運行2×105step，通過交換冷熱原子使系統內部形成一定的溫差，并統計系統產生的熱流和響應的溫差，最后通過熱流與溫度梯度的比值求得納米帶的熱導率。

3.1　軸向長度、寬度對AGYNRs熱導率的影響

本文首先在室溫(300K)下模擬相同寬度的AGYNRs熱導率與其長度間的關系，分別計算了寬固定為2.4699，3.7409，4.9398nm，長10.5，21，30.4，42,63和84nm的扶手椅型石墨炔納米帶的熱導率，計算結果如圖5所示。由圖5可以得到：同一寬度下，石墨炔納米帶熱導率隨著AGYNRs長度的增加而增加，并且直接到84nm時，熱導率為26.7701W/m·K，未收斂，說明石墨炔的分子平均自由程(polyatomic molecule mean free path)遠大于84nm。本文模擬的最大尺寸長為84nm，可見熱導率的大小與尺寸有關[15]，由于分子動力學模擬的尺寸是有限的，當模擬的尺寸小于其分子的平均自由程時會出現Casimir限制效應[7]，使得模擬得到的熱導率要比實驗值小，隨著長度的繼續增加，熱導率會趨于一個穩定值。不同長度下的AGYNRs，長度大的橫截面積大。在AGYNRs中，聲子的散射主要受邊界的限制，根據聲子氣動模型和彈射聲子模型，相同條件下熱量的輸運主要依靠原子的振動，橫截面積越大，其對聲子面內傳輸的影響越小，故橫截面積大的熱導率較大。由圖5還可以看出，寬度為2.4699，3.7409，4.9398nm 三個尺寸下，寬度對熱導率的影響不大。

圖5　完整石墨炔熱導率和乙炔鏈上單空位缺陷對比

3.2　勢函數對ZGYNRs熱導率的影響

圖6是鋸齒型石墨炔納米帶在不同勢函數時計算得到的熱導率,其中鋸齒型納米帶的寬度為63nm。從圖中可以看出,該尺寸下，采用Airebo勢函數計算得到的熱導率高于采用Tersoff勢函數計算得到的熱導率。

4　結語

本文使用非平衡態分子動力學方法計算了AGYNRs的熱導率。在室溫下，對于相同寬度的同種類型石墨炔納米帶，熱導率隨著納米帶長度的增加而逐漸增大。采用Airebo勢函數計算得到的熱導率高于采用Tersoff計算得到的熱導率。

圖6　利用Airebo勢函數與Tersoff勢函數得到石墨炔納米帶熱導率對比

[1]Baughman R H,Eckhardt H,Kertesz M.Structure-property predictions for new planar forms of carbon:Layered phases containing sp2 and sp atoms[J].The Journal of chemical physics,1987,87(11):6687-6699.

[2]Li J,Porter L,Yip S 1998 J.Nucl.Mater.255 139

[3]Zhang Y Y,Pei Q X,Wang C M.A molecular dynamics investigation on thermal conductivity of graphynes[J].Computational Materials Science,2012,65:406-410.

[4]Ouyang T,Chen Y,Liu L M,et al.Thermal transport in graphyne nanoribbons[J].Physical Review B,2012,85(23):235436.

[5]溫志宏.石墨炔納米管熱導率的分子動力學模擬[D].湘潭大學,2014.

[6]Stuart S J,Tutein A B,Harrison J A.The J Chem Phys,2000,112(14):6472-6486.

[7]Schelling P K,Phillpot S R,Keblinski P.Phys Rev B,2002,65(14):144306.

[9]Berber S,Kwon Y K,Tománek D 2000 Phys Rev Lett.84,4613

[10]Muller-Plathe 1999 Phys.Rev.E 59 4894.

[11]Florian Muller-Plathe 1997 J.Chem.Phys.106 6082.

[12]趙晗,周麗娜,魏東山,等.石墨炔納米管熱導率的分子動力學模擬研究[C]//中國化學會第 29 屆學術年會摘要集——第27 分會:多孔功能材料,2014.

[13]楊平，王曉亮，李培，等.氮摻雜和空位對石墨烯納米帶熱導率影響的分子動力學模擬[J].物理學報，2012,61(7)：76501-076501.

[14]魏志勇,畢可東,陳云飛.石墨烯納米帶熱導率的分子動力學模擬[J].東南大學學報(自然科學版),2010,40(2):306-310.

[15]姚承軍,汪曉明,李瑩瑩,等.一維到二維石墨烯聲子熱傳導的分子動力學模擬[J].揚州大學學報(自然科學版),2013(1):22-26.