石墨烯热导率的分子动力学模拟研究

张正林 魏鑫 杨启容 宫薛菲

摘要： 针对石墨烯所存在的热导率高等优点，本文采用非平衡态分子动力学方法，选用Airebo势函数，建立了扶手型和锯齿形石墨烯纳米带模型，研究了石墨烯模型不同特征尺寸（长度、宽度）、不同势函数、不同温度对石墨烯热导率的影响。研究结果表明，特征尺寸的差异会导致石墨烯热导率不同；选择不同的势函数对计算的石墨烯热导率有显著影响；设定相同的初始条件，石墨烯热导率随温度的升高而明显上升，且在320～340 K内变化趋势显著。本文所研究的石墨烯相关热性能在低维热传导领域具有一定的实际应用价值。

关键词： 石墨烯； 分子动力学模拟； 热导率

中图分类号： TK124； O357.5+3文献标识码： A

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

作者简介： 张正林（1994），男，硕士研究生，主要研究方向为材料导热性能的分子动力学。

通讯作者： 杨启容（1970），女，教授，主要研究方向为能源开发与利用中的传热传质。Email： luyingyi125@163.com石墨烯凭借高热导率、高强度、高弹性的优秀性能逐渐在复合材料领域得到研究和推广，因测量尺度太小以及实验条件受限，石墨烯的实验研究进展较慢。石墨烯作为一种碳原子紧密排列的二维新型材料，由碳原子sp2杂化连接形成，拥有出色的化学物理性能，包括高强度、室温下高速的电子迁移率、高模量和高热导率等[1]。石墨烯厚度为034 nm，仅为一个碳原子层的厚度，使其在纳米电子元器件[23]、化工[4]、复合材料[57]应用等领域拥有广阔的应用前景。石墨烯热导率的研究方法包括实验计算和模拟分析，实验方面包括石墨烯复合材料的性能测试[810]及石墨烯基纳米流体的研究[1113]。微纳米尺度加大了实验难度，而基于驰豫时间的研究又受模型尺寸较大影响，分子动力学方法逐渐在石墨烯的研究中得到推广。分子动力学模拟是在系统边界条件的约束下，通过解原子间相互作用势和在约束作用下的牛顿运动方程，模拟模型在微观尺度各重要参数的变化量。自2004年英国科学家Geim和Novoselov通过胶带反复粘离石墨片得到石墨烯后，国内外学者在使用分子动力学研究石墨烯导热性能方面做了很多工作。王建等人[14]研究了石墨烯层数对石墨烯/聚乙烯界面热阻的影响，表明不同层数石墨烯使石墨烯/聚乙烯界面处出现明显的温度不连续现象；韩超群[15]研究了石墨烯在不平整表面上吸附的分子动力学，表明温度对石墨烯刚性基底的吸附行为有较大影响；谢华清等人[16]用实验测得了石墨烯的热导率，指出自由悬架的石墨烯热导率随温度升高而增大，在50 ℃时温度梯度最高。基于此，本文采用非平衡态分子动力学方法，以石墨烯纳米带为模型，选用Airebo势函数，研究了不同石墨烯特征尺寸（长度、宽度）、选择不同势函数、设定不同初始温度对石墨烯热导率的影响，将模拟结果与石墨烯相关实验结果进行对比。该研究对石墨烯热导率的实验具有更大的指导意义。

1模拟环境设定

本文采用Materials Studio，建立石墨烯纳米带模型，同时采用Lammps软件设计程序并计算。在分子动力学模拟过程中，设置势函数、边界条件和模型等关键信息。

1.1势函数选取

由于势函数在分子动力学计算中起重要作用，选择不同的势函数，体系的势能面会有不同的形状，晶胞内分子受力情况与运动状况也有所区别，进而影响计算结果[17]。在N个粒子组成系统中，势函数为

φ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为模型总原子数。石墨烯的模型为单层碳原子层，以研究碳碳原子间热物性为主，除计算碳碳原子、碳氢原子间相互作用力外，考虑原子间非共价键作用及模型的修正，本次模拟选择Airebo势函数，势函数为

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

式中，EREBOij是REBO势函数；ELJij是LJ势函数；ETORSkijl是四体势Torsion。REBO、LJ和Torsion是3种用于计算不同原子间相互作用的势函数，具体介绍如下：

EREBOij取决于原子的位置和类别，其公式为

EREBOij=VRij+bijVAij（3）

式中，VRij和VAij分别表示原子i与j的排斥力与吸引力，主要由原子间距离rij决定。

LJ势是描述非共价键作用的势函数，是对范德瓦尔斯力的表现，即

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

式中，ε为势阱常数；σ为平衡常数；rij为i和j原子间距离。

Torsion是一种明确的四体势函数，在碳氢化合物结构中描述各种二面角参数，其公式为

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

式中，ω为原子k、i、j、l组成二面角的扭转系数；VTORS是k、i、j、l四原子的势函数；rki是原子k、i间距离；wki是原子k、i相互作用的作用项。

1.2设定初始条件

非平衡态分子动力学模拟是通过对系统施加外界条件，使其在势差作用下发生非平衡传热过程，再通过公式计算其热导率。本次模拟设固定边界条件，即在模型中固定两端原子，在两端分别定义热端和冷端原子的坐标。模型边界条件的设定如图1所示，图1的中间位置为待测区域。

图1模型边界条件的设定系综是被一定宏观性质限定的空间相点的集合，所有的模拟过程都要在各自的系综下进行。本次模拟应用的系综先后有网络虚拟终端（network virtual terminal，NVT）（用于使模型稳定）和网络可视实体（network visible entities，NVE）（用于模拟温度变化）。其中，NVE表示在晶胞内有确定的粒子数N、体积V和总能量E，系综的压强P和温度T可能在某一平均值附近起伏变化，系统为与外界无粒子进出能量交换的孤立系统。NVT系综指在晶胞内有确定的粒子数N、体积V和温度T，系统为恒温的封闭系统。由于系综中物理量的变化范围与其数值大小随着体系的增大而减小，所以针对某一特定宏观体系，不同系综产生的计算结果近似相同，差异可忽略不计。

1.3建立模型

常见的石墨烯模型主要有扶手型和锯齿型，扶手型石墨烯纳米带模型如图2所示，锯齿型石墨烯纳米带模型如图3所示。为减少边缘结构对石墨烯内部热流传递的影响，本文采用扶手型石墨烯模型。在Materials Studio中导入碳原子参数，标准模型为852 nm×492 nm×142 nm，对模型添加cvff力场，并进行能量最小化处理。

在Lammps模拟过程中，通过程序对MS建立的模型进行固定边界条件的定义。假设宽度和厚度方向的模型尺寸为492 nm×142 nm保持不变，长度方向分别设0～05 nm和82～852 nm为固定位置，设定05～1 nm为冷端，77～82 nm为热端。设定温度为500 K，在NVT系综下模拟200万步，使整个模型处于稳定状态，在NVE体系两端增加或减少固定能量Q，再运行200万步，根据模拟结果统计温度分布，并通过将结果导入Origin软件进行最小二乘法拟合，求得温度梯度。

2模拟过程与结果

2.1石墨烯模型特征尺寸（长度、宽度）对热导率的影响

本文采用非平衡态分子动力学模拟，通过对系统施加外界扰动，建立非平衡导热过程。对模型热端和冷端每隔一段时间增加或减去固定能量，一段时间后形成稳定的温度梯度。采用0142 nm作为单层石墨烯模型的厚度，计算热导率的方向是冷端向热端方向。热导率为

K=QS×ΔT（6）

式中，Q为单位时间内增加或减去的热流能量；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可以看出，随着模型长度的增大，石墨烯热导率显著升高，整个变化未出现极值和收敛，呈现出线性变化趋势。利用最小二乘法进行拟合，得石墨烯热导率变化的线性方程为

K=-653+17387L（7）

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

造成这一现象的原因是分子动力学模拟尺寸有限，石墨烯纳米带的分子平均自由程远大于模拟尺寸，可达775 nm，在模拟尺寸范围内得到石墨烯热导率比实验所得小，且石墨烯属于低维材料，低维材料的热导率依赖于体系的尺寸。文献[18]通过化学气相沉积法制备了单层石墨烯，通过测试石墨烯拉曼光谱2 D带的温度敏感性变化，发现单层石墨烯尺寸对热导率测试有较大影响，模拟结果也证实了这点。

由图5可以看出，随着模型宽度的增大，石墨烯热导率明显升高，但趋势与A方向相比较平缓，造成与长度方向变化趋势的差别在于石墨烯纳米带在热流传递过程中，声子在不同方向上的传递速率不同，而且石墨烯模型本身长度与宽度方向上结构的差异，造成石墨烯材料在热量传递上的各向异性。同时，在研究不同影响因素对石墨烯热导率的变化时，应建立统一模型，以排除特征尺度对实验结果的干扰。

2.2不同石墨烯势函数对热导率的影响

建立石墨烯纳米带模型时选择的势函数为Airebo，包括共价键作用势REBO、LJ势及四体势Torsion。在石墨烯纳米带模型中，两原子间作用势与两原子间距离有关。在Airebo势函数中，LJ势和四体势Torsion属于可选择项，取决于两原子间是否受非共价键作用。本次用于对比模拟的是REBO势函数，REBO势函数用于研究短程碳碳原子，碳氢原子，氢氢原子的相互作用（r<02 nm），这些相互作用通过键阶参数，具有很强的配位依赖

图6不同势函数对石墨烯热导率的影响性，但与Airebo势函数相比，REBO势函数计算过程中，原子缺乏长程原子的相互作用，缺乏非共价键的作用项，缺乏由于键角引起的势函数修正，本次模拟计算后，易忽略其他条件对模拟结果的影响。

为研究不同势函数对石墨烯热导率的影响，在Lammps中设定Airebo势函数时，将原先的pairstyle airebo 20 1 1改为pairstyle airebo 20 0 0，即可将势函数改为单纯的REBO势函数，其他条件不变。同样，针对石墨烯纳米带A方向的特征尺寸不同，不同势函数对石墨烯热导率的影响如图6所示。由图6可以看出，在确定模型、初始条件和程序设定的情况下，单纯运用REBO势函数求得的石墨烯热导率，随A方向模型尺寸变化与Airebo势函数（REBO、LJ势及四体势Torsion）有较大差别，在A方向模型尺寸为4～6 nm时，两种势函数求得热导率基本相同，但在A方向模型尺寸大于6 nm之后，相同尺寸的模型REBO势函数求得的石墨烯热导率要比Airebo势函数求得的石墨烯熱导率大，出现这种差异的原因在于REBO势函数是描述原子间共价键的作用势，而附加的LJ势是用来描述原子间的非共价键作用，当两原子间距离小于02 nm时，原子间作用势为REBO势函数，如果两原子间距离大于02 nm，但小于设定的截断距离cutoff（碳碳原子之间为034 nm）两原子间为非共价键作用，即存在LJ势。在本次研究建立的模型中，两原子间最小距离为0246 nm，即两原子间存在描述非共价键作用的LJ势，因此需要设定势函数为Airebo势函数。在分子动力学的研究中，势函数的选择对研究结果起决定性影响，应根据具体研究选择合适势函数。

2.3不同溫度对石墨烯热导率的影响

在石墨烯的分子动力学模拟过程中，温度是重要的设定参数，而在石墨烯的实际应用过程中，温度也是影响石墨烯导热性能的重要因素。许多研究人员用实验方法测定了温度对石墨烯导热性能的影响，文献[19]用实验测试了有悬架的石墨烯和有支撑的石墨烯热导率分别为2 500 W/（m·K）和370 W/（m·K）；文献[20]采用四点测量方法，得到了还原氧化石墨烯的热导率，发现石墨烯的氧化水平极大影响了其热传输特性；文献[16]通过实验研究了自由悬架的石墨烯热传输特性，表明石墨烯纳米带热导率与温度存在很强的相关性，当温度从-75 ℃升高至100 ℃时，其热导率从12621 W/（m·K）增至87732 W/（m·K），且在50 ℃附近石墨烯热导率随温度变化趋势最明显。

通过改变Lammps中的模型初始条件后，运用分子动力学原理计算出石墨烯热导率随温度变化趋势。选用尺寸为852 nm×492 nm×142 nm的扶手型石墨烯纳米带模型，分别在Lammps的Velocity和fix模块设定不同的初始温度，利用最小二乘拟合计算结果求出温度梯度，带入方程，求得不同初始温度对石墨烯热导率的影响

由图7可以看出，在300～600 K的温度变化范围内，石墨烯热导率随温度升高而增大，且在320～340 K附近，横纵坐标的斜率最大，即热导率随温度变化最显著，这与文献[16]的结论相符。

在文献[16]中，单层石墨烯纳米带尺寸为508 μm×385 μm，模拟建立模型尺寸为852 nm×492 nm，数量级相差巨大。根据21，石墨烯尺寸对热导率的测定有较大影响。因计算机计算条件的限制，把纳米尺度改为微米尺度，计算量增加较大。将模拟结论与实验所得热导率相比，热导率随不同温度变化趋势相同。说明分子动力学模拟是可信的，而且可节约实验成本。

3结束语

本文主要研究了不同石墨烯特征尺寸（长度、宽度）、选择不同势函数、设定不同初始温度对石墨烯热导率的影响。石墨烯为低维材料，其热传导性能为各向异性，尺度加大，尤其长度方向加大尺度，有利于增强导热性能；与大多数常规材料不同，石墨烯热导率随温度的增加表现为非线性增长，在320～340 K附近增速最大；本次研究中模拟模型尺寸与实验材料尺寸相差数量级巨大，所得到的石墨烯热导率随温度变化趋势相同，但具体数值不同。今后应在建模方向继续研究，使模拟结果更精准，对实验具有较大的指导意义；由于石墨烯独特的导热性能，可以与其他材料复合，有效提高材料导热性能。

参考文献：

[1]方明. 石墨烯基纳米复合材料的制备及性能[D]. 上海： 复旦大学， 2011.

[2]闻雷， 陈静， 罗洪泽， 等. 石墨烯在柔性锂离子电池中的应用及前景[J]. 科学通报， 2015， 60（7）： 630644.

[3]汪建德， 彭同江， 鲜海洋， 等. 三维还原氧化石墨烯/聚苯胺复合材料的制备及其超级电容性能[J]. 物理化学学报， 2015， 31（1）： 9098.

[4]何大方， 吴健， 刘战剑， 等. 面向应用的石墨烯制备研究进展[J]. 化工学报， 2015， 66（8）： 28882894.

[5]陈中华， 王建川， 余飞， 等. 氧化石墨烯/密胺树脂复合材料的制备及其热性能研究[J]. 功能材料， 2015， 46（1）： 11251128.

[6]沈小军， 孟令轩， 付绍云. 石墨烯多壁碳纳米管协同增强环氧树脂复合材料的低温力学性能[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.