石墨烯熱導率的分子動力學模擬研究

張正林 魏鑫 楊啟容 宮薛菲

摘要： 針對石墨烯所存在的熱導率高等優(yōu)點，本文采用非平衡態(tài)分子動力學方法，選用Airebo勢函數(shù)，建立了扶手型和鋸齒形石墨烯納米帶模型，研究了石墨烯模型不同特征尺寸（長度、寬度）、不同勢函數(shù)、不同溫度對石墨烯熱導率的影響。研究結(jié)果表明，特征尺寸的差異會導致石墨烯熱導率不同；選擇不同的勢函數(shù)對計算的石墨烯熱導率有顯著影響；設定相同的初始條件，石墨烯熱導率隨溫度的升高而明顯上升，且在320～340 K內(nèi)變化趨勢顯著。本文所研究的石墨烯相關(guān)熱性能在低維熱傳導領(lǐng)域具有一定的實際應用價值。

關(guān)鍵詞： 石墨烯； 分子動力學模擬； 熱導率

中圖分類號： TK124； O357.5+3文獻標識碼： A

收稿日期： 20170518； 修回日期： 20170829

作者簡介： 張正林（1994），男，碩士研究生，主要研究方向為材料導熱性能的分子動力學。

通訊作者： 楊啟容（1970），女，教授，主要研究方向為能源開發(fā)與利用中的傳熱傳質(zhì)。Email： luyingyi125@163.com石墨烯憑借高熱導率、高強度、高彈性的優(yōu)秀性能逐漸在復合材料領(lǐng)域得到研究和推廣，因測量尺度太小以及實驗條件受限，石墨烯的實驗研究進展較慢。石墨烯作為一種碳原子緊密排列的二維新型材料，由碳原子sp2雜化連接形成，擁有出色的化學物理性能，包括高強度、室溫下高速的電子遷移率、高模量和高熱導率等[1]。石墨烯厚度為034 nm，僅為一個碳原子層的厚度，使其在納米電子元器件[23]、化工[4]、復合材料[57]應用等領(lǐng)域擁有廣闊的應用前景。石墨烯熱導率的研究方法包括實驗計算和模擬分析，實驗方面包括石墨烯復合材料的性能測試[810]及石墨烯基納米流體的研究[1113]。微納米尺度加大了實驗難度，而基于馳豫時間的研究又受模型尺寸較大影響，分子動力學方法逐漸在石墨烯的研究中得到推廣。分子動力學模擬是在系統(tǒng)邊界條件的約束下，通過解原子間相互作用勢和在約束作用下的牛頓運動方程，模擬模型在微觀尺度各重要參數(shù)的變化量。自2004年英國科學家Geim和Novoselov通過膠帶反復粘離石墨片得到石墨烯后，國內(nèi)外學者在使用分子動力學研究石墨烯導熱性能方面做了很多工作。王建等人[14]研究了石墨烯層數(shù)對石墨烯/聚乙烯界面熱阻的影響，表明不同層數(shù)石墨烯使石墨烯/聚乙烯界面處出現(xiàn)明顯的溫度不連續(xù)現(xiàn)象；韓超群[15]研究了石墨烯在不平整表面上吸附的分子動力學，表明溫度對石墨烯剛性基底的吸附行為有較大影響；謝華清等人[16]用實驗測得了石墨烯的熱導率，指出自由懸架的石墨烯熱導率隨溫度升高而增大，在50 ℃時溫度梯度最高?；诖耍疚牟捎梅瞧胶鈶B(tài)分子動力學方法，以石墨烯納米帶為模型，選用Airebo勢函數(shù)，研究了不同石墨烯特征尺寸（長度、寬度）、選擇不同勢函數(shù)、設定不同初始溫度對石墨烯熱導率的影響，將模擬結(jié)果與石墨烯相關(guān)實驗結(jié)果進行對比。該研究對石墨烯熱導率的實驗具有更大的指導意義。

1模擬環(huán)境設定

本文采用Materials Studio，建立石墨烯納米帶模型，同時采用Lammps軟件設計程序并計算。在分子動力學模擬過程中，設置勢函數(shù)、邊界條件和模型等關(guān)鍵信息。

1.1勢函數(shù)選取

由于勢函數(shù)在分子動力學計算中起重要作用，選擇不同的勢函數(shù)，體系的勢能面會有不同的形狀，晶胞內(nèi)分子受力情況與運動狀況也有所區(qū)別，進而影響計算結(jié)果[17]。在N個粒子組成系統(tǒng)中，勢函數(shù)為

φij=∑iU1（ri）+∑i∑j>iU2（ri，rj）+∑i∑j>i∑k>j>iU3（ri，rj，rk）+…+∑…∑Un（ri，rj，rk，…，rn）（1）

式中，ri表示第i個粒子的空間位置；Un表示體作用勢；i、j為模型內(nèi)任意兩原子；n為模型總原子數(shù)。石墨烯的模型為單層碳原子層，以研究碳碳原子間熱物性為主，除計算碳碳原子、碳氫原子間相互作用力外，考慮原子間非共價鍵作用及模型的修正，本次模擬選擇Airebo勢函數(shù)，勢函數(shù)為

E=12∑i∑j≠i[EREBOij+ELJij+∑k≠i、j∑l≠i、j、kETORSkijl]（2）

式中，EREBOij是REBO勢函數(shù)；ELJij是LJ勢函數(shù)；ETORSkijl是四體勢Torsion。REBO、LJ和Torsion是3種用于計算不同原子間相互作用的勢函數(shù)，具體介紹如下：

EREBOij取決于原子的位置和類別，其公式為

EREBOij=VRij+bijVAij（3）

式中，VRij和VAij分別表示原子i與j的排斥力與吸引力，主要由原子間距離rij決定。

LJ勢是描述非共價鍵作用的勢函數(shù)，是對范德瓦爾斯力的表現(xiàn)，即

VLJ=4εσrij12-σrij6（4）

式中，ε為勢阱常數(shù)；σ為平衡常數(shù)；rij為i和j原子間距離。

Torsion是一種明確的四體勢函數(shù)，在碳氫化合物結(jié)構(gòu)中描述各種二面角參數(shù)，其公式為

ETORSkijl=wki（rki）wij（rij）wjl（rjl）VTORS（ωkijl）（5）

式中，ω為原子k、i、j、l組成二面角的扭轉(zhuǎn)系數(shù)；VTORS是k、i、j、l四原子的勢函數(shù)；rki是原子k、i間距離；wki是原子k、i相互作用的作用項。

1.2設定初始條件

非平衡態(tài)分子動力學模擬是通過對系統(tǒng)施加外界條件，使其在勢差作用下發(fā)生非平衡傳熱過程，再通過公式計算其熱導率。本次模擬設固定邊界條件，即在模型中固定兩端原子，在兩端分別定義熱端和冷端原子的坐標。模型邊界條件的設定如圖1所示，圖1的中間位置為待測區(qū)域。

圖1模型邊界條件的設定系綜是被一定宏觀性質(zhì)限定的空間相點的集合，所有的模擬過程都要在各自的系綜下進行。本次模擬應用的系綜先后有網(wǎng)絡虛擬終端（network virtual terminal，NVT）（用于使模型穩(wěn)定）和網(wǎng)絡可視實體（network visible entities，NVE）（用于模擬溫度變化）。其中，NVE表示在晶胞內(nèi)有確定的粒子數(shù)N、體積V和總能量E，系綜的壓強P和溫度T可能在某一平均值附近起伏變化，系統(tǒng)為與外界無粒子進出能量交換的孤立系統(tǒng)。NVT系綜指在晶胞內(nèi)有確定的粒子數(shù)N、體積V和溫度T，系統(tǒng)為恒溫的封閉系統(tǒng)。由于系綜中物理量的變化范圍與其數(shù)值大小隨著體系的增大而減小，所以針對某一特定宏觀體系，不同系綜產(chǎn)生的計算結(jié)果近似相同，差異可忽略不計。

1.3建立模型

常見的石墨烯模型主要有扶手型和鋸齒型，扶手型石墨烯納米帶模型如圖2所示，鋸齒型石墨烯納米帶模型如圖3所示。為減少邊緣結(jié)構(gòu)對石墨烯內(nèi)部熱流傳遞的影響，本文采用扶手型石墨烯模型。在Materials Studio中導入碳原子參數(shù)，標準模型為852 nm×492 nm×142 nm，對模型添加cvff力場，并進行能量最小化處理。

在Lammps模擬過程中，通過程序?qū)S建立的模型進行固定邊界條件的定義。假設寬度和厚度方向的模型尺寸為492 nm×142 nm保持不變，長度方向分別設0～05 nm和82～852 nm為固定位置，設定05～1 nm為冷端，77～82 nm為熱端。設定溫度為500 K，在NVT系綜下模擬200萬步，使整個模型處于穩(wěn)定狀態(tài)，在NVE體系兩端增加或減少固定能量Q，再運行200萬步，根據(jù)模擬結(jié)果統(tǒng)計溫度分布，并通過將結(jié)果導入Origin軟件進行最小二乘法擬合，求得溫度梯度。

2模擬過程與結(jié)果

2.1石墨烯模型特征尺寸（長度、寬度）對熱導率的影響

本文采用非平衡態(tài)分子動力學模擬，通過對系統(tǒng)施加外界擾動，建立非平衡導熱過程。對模型熱端和冷端每隔一段時間增加或減去固定能量，一段時間后形成穩(wěn)定的溫度梯度。采用0142 nm作為單層石墨烯模型的厚度，計算熱導率的方向是冷端向熱端方向。熱導率為

K=QS×ΔT（6）

式中，Q為單位時間內(nèi)增加或減去的熱流能量；S為模型橫截面積；ΔT為溫度梯度。

由于石墨烯是碳原子sp2雜化連接形成的單層碳原子層，其模型厚度為0142 nm，在三維模型中特征尺寸的變化有模型長度和模型寬度。首先對石墨烯納米帶模型進行等比例擴展，在固定模型寬度和厚度為492nm×142 nm后，模型長度分別設置特征尺寸為426，639，852，1065，1278 nm，建立模型后，導出Lammps模擬計算，帶入公式求熱導率。模型不同長度對熱導率影響如圖4所示。同理，在石墨烯納米帶模型中，固定模型長度和厚度分別為852 nm×142 nm，模型寬度分別設置特征尺寸為12，246，369，492，615 nm，建立模型后，導出Lammps模擬計算，帶入公式求熱導率。模型不同寬度對熱導率影響如圖5所示。

由圖4可以看出，隨著模型長度的增大，石墨烯熱導率顯著升高，整個變化未出現(xiàn)極值和收斂，呈現(xiàn)出線性變化趨勢。利用最小二乘法進行擬合，得石墨烯熱導率變化的線性方程為

K=-653+17387L（7）

式中，L為模型A方向尺寸，nm。

造成這一現(xiàn)象的原因是分子動力學模擬尺寸有限，石墨烯納米帶的分子平均自由程遠大于模擬尺寸，可達775 nm，在模擬尺寸范圍內(nèi)得到石墨烯熱導率比實驗所得小，且石墨烯屬于低維材料，低維材料的熱導率依賴于體系的尺寸。文獻[18]通過化學氣相沉積法制備了單層石墨烯，通過測試石墨烯拉曼光譜2 D帶的溫度敏感性變化，發(fā)現(xiàn)單層石墨烯尺寸對熱導率測試有較大影響，模擬結(jié)果也證實了這點。

由圖5可以看出，隨著模型寬度的增大，石墨烯熱導率明顯升高，但趨勢與A方向相比較平緩，造成與長度方向變化趨勢的差別在于石墨烯納米帶在熱流傳遞過程中，聲子在不同方向上的傳遞速率不同，而且石墨烯模型本身長度與寬度方向上結(jié)構(gòu)的差異，造成石墨烯材料在熱量傳遞上的各向異性。同時，在研究不同影響因素對石墨烯熱導率的變化時，應建立統(tǒng)一模型，以排除特征尺度對實驗結(jié)果的干擾。

2.2不同石墨烯勢函數(shù)對熱導率的影響

建立石墨烯納米帶模型時選擇的勢函數(shù)為Airebo，包括共價鍵作用勢REBO、LJ勢及四體勢Torsion。在石墨烯納米帶模型中，兩原子間作用勢與兩原子間距離有關(guān)。在Airebo勢函數(shù)中，LJ勢和四體勢Torsion屬于可選擇項，取決于兩原子間是否受非共價鍵作用。本次用于對比模擬的是REBO勢函數(shù)，REBO勢函數(shù)用于研究短程碳碳原子，碳氫原子，氫氫原子的相互作用（r<02 nm），這些相互作用通過鍵階參數(shù)，具有很強的配位依賴

圖6不同勢函數(shù)對石墨烯熱導率的影響性，但與Airebo勢函數(shù)相比，REBO勢函數(shù)計算過程中，原子缺乏長程原子的相互作用，缺乏非共價鍵的作用項，缺乏由于鍵角引起的勢函數(shù)修正，本次模擬計算后，易忽略其他條件對模擬結(jié)果的影響。

為研究不同勢函數(shù)對石墨烯熱導率的影響，在Lammps中設定Airebo勢函數(shù)時，將原先的pairstyle airebo 20 1 1改為pairstyle airebo 20 0 0，即可將勢函數(shù)改為單純的REBO勢函數(shù)，其他條件不變。同樣，針對石墨烯納米帶A方向的特征尺寸不同，不同勢函數(shù)對石墨烯熱導率的影響如圖6所示。由圖6可以看出，在確定模型、初始條件和程序設定的情況下，單純運用REBO勢函數(shù)求得的石墨烯熱導率，隨A方向模型尺寸變化與Airebo勢函數(shù)（REBO、LJ勢及四體勢Torsion）有較大差別，在A方向模型尺寸為4～6 nm時，兩種勢函數(shù)求得熱導率基本相同，但在A方向模型尺寸大于6 nm之后，相同尺寸的模型REBO勢函數(shù)求得的石墨烯熱導率要比Airebo勢函數(shù)求得的石墨烯熱導率大，出現(xiàn)這種差異的原因在于REBO勢函數(shù)是描述原子間共價鍵的作用勢，而附加的LJ勢是用來描述原子間的非共價鍵作用，當兩原子間距離小于02 nm時，原子間作用勢為REBO勢函數(shù)，如果兩原子間距離大于02 nm，但小于設定的截斷距離cutoff（碳碳原子之間為034 nm）兩原子間為非共價鍵作用，即存在LJ勢。在本次研究建立的模型中，兩原子間最小距離為0246 nm，即兩原子間存在描述非共價鍵作用的LJ勢，因此需要設定勢函數(shù)為Airebo勢函數(shù)。在分子動力學的研究中，勢函數(shù)的選擇對研究結(jié)果起決定性影響，應根據(jù)具體研究選擇合適勢函數(shù)。

2.3不同溫度對石墨烯熱導率的影響

在石墨烯的分子動力學模擬過程中，溫度是重要的設定參數(shù)，而在石墨烯的實際應用過程中，溫度也是影響石墨烯導熱性能的重要因素。許多研究人員用實驗方法測定了溫度對石墨烯導熱性能的影響，文獻[19]用實驗測試了有懸架的石墨烯和有支撐的石墨烯熱導率分別為2 500 W/（m·K）和370 W/（m·K）；文獻[20]采用四點測量方法，得到了還原氧化石墨烯的熱導率，發(fā)現(xiàn)石墨烯的氧化水平極大影響了其熱傳輸特性；文獻[16]通過實驗研究了自由懸架的石墨烯熱傳輸特性，表明石墨烯納米帶熱導率與溫度存在很強的相關(guān)性，當溫度從-75 ℃升高至100 ℃時，其熱導率從12621 W/（m·K）增至87732 W/（m·K），且在50 ℃附近石墨烯熱導率隨溫度變化趨勢最明顯。

通過改變Lammps中的模型初始條件后，運用分子動力學原理計算出石墨烯熱導率隨溫度變化趨勢。選用尺寸為852 nm×492 nm×142 nm的扶手型石墨烯納米帶模型，分別在Lammps的Velocity和fix模塊設定不同的初始溫度，利用最小二乘擬合計算結(jié)果求出溫度梯度，帶入方程，求得不同初始溫度對石墨烯熱導率的影響

由圖7可以看出，在300～600 K的溫度變化范圍內(nèi)，石墨烯熱導率隨溫度升高而增大，且在320～340 K附近，橫縱坐標的斜率最大，即熱導率隨溫度變化最顯著，這與文獻[16]的結(jié)論相符。

在文獻[16]中，單層石墨烯納米帶尺寸為508 μm×385 μm，模擬建立模型尺寸為852 nm×492 nm，數(shù)量級相差巨大。根據(jù)21，石墨烯尺寸對熱導率的測定有較大影響。因計算機計算條件的限制，把納米尺度改為微米尺度，計算量增加較大。將模擬結(jié)論與實驗所得熱導率相比，熱導率隨不同溫度變化趨勢相同。說明分子動力學模擬是可信的，而且可節(jié)約實驗成本。

3結(jié)束語

本文主要研究了不同石墨烯特征尺寸（長度、寬度）、選擇不同勢函數(shù)、設定不同初始溫度對石墨烯熱導率的影響。石墨烯為低維材料，其熱傳導性能為各向異性，尺度加大，尤其長度方向加大尺度，有利于增強導熱性能；與大多數(shù)常規(guī)材料不同，石墨烯熱導率隨溫度的增加表現(xiàn)為非線性增長，在320～340 K附近增速最大；本次研究中模擬模型尺寸與實驗材料尺寸相差數(shù)量級巨大，所得到的石墨烯熱導率隨溫度變化趨勢相同，但具體數(shù)值不同。今后應在建模方向繼續(xù)研究，使模擬結(jié)果更精準，對實驗具有較大的指導意義；由于石墨烯獨特的導熱性能，可以與其他材料復合，有效提高材料導熱性能。

參考文獻：

[1]方明. 石墨烯基納米復合材料的制備及性能[D]. 上海： 復旦大學， 2011.

[2]聞雷， 陳靜， 羅洪澤， 等. 石墨烯在柔性鋰離子電池中的應用及前景[J]. 科學通報， 2015， 60（7）： 630644.

[3]汪建德， 彭同江， 鮮海洋， 等. 三維還原氧化石墨烯/聚苯胺復合材料的制備及其超級電容性能[J]. 物理化學學報， 2015， 31（1）： 9098.

[4]何大方， 吳健， 劉戰(zhàn)劍， 等. 面向應用的石墨烯制備研究進展[J]. 化工學報， 2015， 66（8）： 28882894.

[5]陳中華， 王建川， 余飛， 等. 氧化石墨烯/密胺樹脂復合材料的制備及其熱性能研究[J]. 功能材料， 2015， 46（1）： 11251128.

[6]沈小軍， 孟令軒， 付紹云. 石墨烯多壁碳納米管協(xié)同增強環(huán)氧樹脂復合材料的低溫力學性能[J]. 復合材料學報， 2015， 32（1）： 2126.

[7]何飛雄， 卞軍， 藺海蘭， 等. 功能化納米石墨烯片/PPPPgMAH復合材料的制備與表征[J]. 復合材料學報， 2015， 32（1）： 4753.

[8]獨濤， 張洪迪， 范同祥. 石墨烯/金屬復合材料的研究進展[J]. 材料導報， 2015， 29（3）： 121129.

[9]Chen W Y， Tao X Q， Wang H S， et al. Research Progress on Application of Graphene and its Composite Material[J]. Material Sciences， 2015， 5（3）： 6271.

[10]于偉， 謝華清， 陳立飛， 等. 石墨烯制備方法及粒徑對復合材料熱導率的影響[J]. 工程熱物理學報， 2016， 37（1）： 168172.

[11]Lee G J， Chang K R. Enhanced Thermal Conductivity of Nanofluids Containing Graphenenanoplatelets Prepared by Ultrasound Irradiation[J]. Journal of Materials Science， 2014， 49（4）： 15061511.

[12]Mehrali M， Sadeghinezhad E， Latibari S T， et al. Preparation， Characterization， Viscosity， and Thermal Conductivity of NitrogenDoped Graphene Aqueous Nanofluids[J]. Journal of Materials Science， 2014， 49（20）： 71567171.

[13]Park S S， Kim N J. Influence of the Oxidation Treatment and the Average Particle Diameter of Graphene for Thermal Conductivity Enhancement[J]. Journal of Industrial & Engineering Chemistry， 2014， 20（4）： 19111915.

[14]王建， 苗瑞珍， 劉偉帥， 等. 石墨烯/聚乙烯界面導熱性能的分子動力學模擬[J]. 分子科學學報， 2016， 32（2）： 156161.

[15]韓超群. 石墨烯在不平整表面上吸附的分子動力學研究[D]. 上海： 華東師范大學， 2011.

[16]Xie H， Chen L， Yu W， et al. Temperature Dependent Thermal Conductivity of a FreeStanding Graphenenanoribbon[J]. Applied Physics Letters， 2013， 102（11）： 6661979.

[17]Lukes J R， Liang X G， Tien C L. Molecular Dynamics Study of Solid ThinFilm Thermal Conductivity[J]. Journal of Heat Transfer， 2000， 122（3）： 536543.

[18]Chen S， Moore A， Cai W， et al. Raman Measurements of Thermal Transport in Suspended Monolayer Graphene of Variable Sizes in Vacuum and Gaseous Environments[J]. ACS NaNo， 2011， 5（1）： 321328.

[19]Cai W， Moore A， Chen S， et al. Thermal Transport in Suspended and Supported Monolayer Graphene Grown by Chemical Vapor Deposition[J]. Nano Letters， 2010， 10（5）： 16451651.

[20]Schwamb T， Burg B R， Schirmer N C， et al. An Electrical Method for the Measurement of the Thermal and Electrical Conductivity of Reduced Graphene Oxide Nanostructures[J]. Nanotechnology， 2009， 20（40）： 15.