韩强 黄凌燕

( 华南理工大学 土木与交通学院，广东 广州510640)

自从碳纳米管［1］( CNTs) 被发现以来，由于其在力学、化学、热力学和电学等方面具有的独特性质，使得纳米研究成为当今世界的热点问题之一. 近年来，随着单层石墨烯［2］的制备成功，引发了新的碳质材料的研究热潮.单层石墨烯是以sp2轨道杂化的碳原子形成的蜂窝状六角平面晶体，厚度仅为0.335 nm，其中C—C 键键长约为0.142 nm，为目前世界上存在的最薄的材料. 理论计算［3］表明，单层石墨烯的弹性模量在700 ～1 000 GPa 之间，是普通钢材的5 倍.

关于石墨烯基础性能的研究是纳米材料的研究热点之一，辛浩［4］提出采用三层正交异性单向纤维的叠层作为石墨烯的等效模型，从理论上计算了石墨烯的弹性模量，为0.97 TPa.Gao 等［5］采用量子力学和分子动力学的方法得到锯齿型和扶手椅型石墨烯的杨氏模量，分别为0.6 TPa 和1.1 TPa. Ni 等［6］研究了室温下石墨烯弹性模量的尺寸效应，认为当尺寸变化时，锯齿型石墨烯的弹性模量值始终高于扶手椅型石墨烯的弹性模量值. 韩同伟等［7］通过分子动力学方法对单层石墨烯薄膜的拉伸变形进行了研究.现有的研究中，有关尺寸效应和层数效应的文章并不多见，还没有形成统一的结论.

文中基于分子动力学的方法，模拟不同尺寸单层石墨烯在两个方向上的拉伸变形，并研究拉伸过程中的应力－应变关系、力学性能以及尺寸变化对石墨烯拉伸性能的影响.在尺寸相同的条件下，通过改变石墨烯的层数，研究层数与石墨烯拉伸力学性能之间的关系.

1 分子动力学模拟

1.1 单层、多层石墨烯的物理模型

文中采用的单层石墨烯薄膜尺寸列于表1 中.其中，Lx表示单层石墨烯x 方向的长度，Ly表示单层石墨烯y 方向的长度.

图1 所示为单层石墨烯的原子构型图.可以看出:当对薄膜进行x 方向拉伸时，相当于分析扶手椅型石墨烯的拉伸性能;当对薄膜进行y 方向拉伸时，相当于分析锯齿型石墨烯的拉伸性能.

表1 单层石墨烯薄膜的尺寸Table 1 Size of monolayer graphene sheet

图1 单层石墨烯的原子构型图Fig.1 Atomic model of monolayer graphene sheet

多层石墨烯是由单层石墨烯在z 方向进行Bernal( AB)［8］堆叠后得到的，多层石墨烯的层间距为0.335 nm.

将表1 中模型编号为6 的单层石墨烯薄膜在z方向堆叠，得文中模拟的多层石墨烯模型尺寸，列于表2.

表2 多层石墨烯模型尺寸Table 2 Size of multilayer graphene sheet model

1.2 分子动力学模拟方法

分子动力学模拟的基本原理是通过原子间的相互作用势，求出每个原子所受到的力，在选定的时间步长、边界条件和初始条件下，对有限数目的分子( 原子) 建立其牛顿动力学方程组，用数值方法求解，得到这些原子的运动轨迹和速度分布，然后对足够长时间的结果求统计平均，从而得到所需要的宏观物理量和力学量.

势函数的选取是分子动力学模拟的关键，研究单层石墨烯时选用不考虑层间作用的Tersoff 势函数［9］.Tersoff 势函数是一种三体势函数，它可以很好地模拟C—C 共价键的各种特性，包括键长、键角、键能、晶格常数和键的断裂重组等动态行为，能够较真实地反应碳元素所构成的固态材料的物理性质.

原子间相互作用Tersoff 势函数为

式中，E 为体系的总能量; Ei为单个原子能量; Eij为i、j 原子之间的成键能量，分别为对势的吸引项和排斥项，fC是光滑截断函数，bij为键序函数，rij是i、j 原子之间的键长.

Tersoff 势函数形式上是一个二体势，实际上是一个多体势，因为系数bij并非一个常数，而是一个依赖于i、j 原子位置并与i 粒子周围其他近邻原子有关的多体函数项.

研究多层石墨烯时，考虑到层间作用力的影响，需要在Tersoff 势函数的基础上考虑Lennard-Jones( L-J) 势函数.

L-J 势函数形式如下:

式中:右边括号内第1 项是排斥项，第2 项为吸引项;参数r 是连接i 和j 原子的位置矢量的模;ε 和l分别表示能量标度参数和碰撞直径参数，ε 的意义等同于离解能，l 则可以理解为两个原子间相互作用势为零时的原子间距.

目前，运用最为广泛的积分算法是Verlet 算法［10］，但是速度处理不是很方便，可能导致不同物理量计算精度不同. 在LAMMPS 中，默认的积分算法是Velocity-Verlet，它可以解决Verlet 算法的不足，具体形式如下:

式中，t 为时间，δt 为时间变量;r 为原子位置，v 为原子的速度，m 为原子质量，F 表示原子所处的力场.这种算法可以同时获得相同精度的原子位置和速度，在每步积分中只要存储一个时刻的状态变量，模拟稳定性较好，同时允许较大的时间步长，因此在分子动力学积分运算中得到广泛的运用.

模拟过程中，选取C 原子质量为12. 控制x、y方向为自由边界条件，z 方向为周期性边界条件.采用Nose-Hoover 热浴法［11-12］进行温度调节，控制温度为0.01 K，设计时间步长为1 fs.

当进行x 方向拉伸模拟时，固定薄膜左端碳原子，先进行充分的弛豫，稳定后对右端碳原子施加x方向拉伸应变载荷，每步应变量为0.001，弛豫步长为3000 步.持续加载，直至薄膜被拉断.

当进行y 方向拉伸模拟时，固定薄膜下端碳原子，对上端原子施加y 方向的拉伸应变载荷，其他步骤与x 方向拉伸相同.

2 模拟结果

2.1 拉伸过程的应力－应变曲线

材料的应力－应变关系是研究材料特性的重点之一.对表1 中列出的不同尺寸单层石墨烯和表2中列出的不同厚度多层石墨烯进行拉伸模拟，得到相应的应力－应变关系.

在计算石墨烯拉伸变形性能时，需要采用系统的平均应力.由离散原子构成的纳米系统中应力的概念与连续介质中应力的概念有所不同，Jin 等［13］在研究纳米管力学性能时给出了纳米系统下平均应力的计算公式:

式中:σαβ表示在笛卡尔坐标下系统原子水平的平均应力;分别表示原子i 与原子j 之间的相互作用力和距离; V0为模型初始状态下的体积. 对于单层石墨烯薄膜有V0=LxLyh，其中，Lx和Ly为单层石墨烯的长度和宽度值，h 为单层石墨烯的厚度，取为0.335 nm［14-16］.

图2 为表1 中模型编号为6 的单层石墨烯薄膜在两个方向上的应力－ 应变曲线. 图中，x 方向的拉伸曲线对应扶手椅型石墨烯的拉伸性能，y 方向的拉伸曲线对应锯齿型石墨烯的拉伸性能.

图3 为1 层、2 层和4 层石墨烯拉伸过程中的应力－应变曲线.为了便于比较层数变化对应力－应变关系的影响，将1、2、4 层的拉伸应力－应变关系曲线画在同一幅图上.

图2 和图3 表明，多层石墨烯和单层石墨烯具有相同的应力－应变曲线形态，且尺寸变化和层数变化对曲线形态几乎不产生影响. 石墨烯在两个方向的拉伸过程中均经历了弹性变形、屈服阶段、强化阶段和局部变形4个阶段，但两个方向所表现出来的力学性能相差较大，x 方向( 扶手椅型) 的拉伸屈服强度、断裂强度以及极限应变值均高于y 方向( 锯齿型) ，表明了石墨烯各向异性的特征.

图2 单层石墨烯拉伸的应力－应变曲线Fig. 2 Tensile stress-strain curves of monolayer graphene sheet

图3 多层石墨烯拉伸的应力－应变曲线Fig.3 Tensile stress-strain curves of multilayer graphene sheets

2.2 单层石墨烯的弹性模量

这里，仅对图2 中应变在0 ～0.03 段的曲线进行最小二乘法拟合，计算单层石墨烯的弹性模量如图4 所示为拟合范围内的应力－应变关系.

图4 拟合范围内的应力－应变曲线Fig.4 Tensile stress-strain curves in the fitting range

用Ex表示x 方向拉伸时对应的扶手椅型石墨烯的弹性模量，用Ey表示y 方向拉伸时对应的锯齿型石墨烯的弹性模量.通过计算，将不同尺寸单层石墨烯薄膜的弹性模量列于表3.

表3 单层石墨烯的弹性模量Table 3 Elastic modulus of monolayer graphene sheet

从表3 可以看出，在单层石墨烯薄膜尺寸由小增大的过程中，其弹性模量并未发生太大的改变，Ex总是在1078.02 GPa 附近，而Ey总是在1041.53 GPa附近，且Ex和Ey两者的差距始终小于6%. 可见，尺寸变化对单层石墨烯弹性模量的影响很小，且在尺寸变化过程中Ex始终大于Ey.

结合图4，当拉伸应变在0 ～0.03 范围内时，单层石墨烯两个方向的应力－应变曲线几乎是线性重合的，且在该阶段两个方向的弹性模量相差不超过6%.于是可以认为在拉伸变形的线弹性阶段，石墨烯薄膜是各向同性材料.

2.3 多层石墨烯的弹性模量

与单层石墨烯计算弹性模量值的方法一致，对图3 中应变在0.00 ～0.03 段的曲线进行最小二乘法拟合，将计算结果列于表4 中.

表4 不同层数石墨烯的弹性模量值Table 4 Elastic modulus of graphene sheets with different layers

从表4 可以看出，多层石墨烯在两个方向上的弹性模量与单层石墨烯的接近，均在1 050 GPa 附近.单层石墨烯两个方向弹性模量的差距为1.35%，而多层石墨烯的弹性模量差距均小于0.4%.可见，石墨烯层数的增加有助于减小两个方向弹性模量的差距，在拉伸变形的线弹性阶段，多层石墨烯比单层石墨烯表现出更为强烈的各向同性.

2.4 拉伸变形破坏形态分析

图5 为单层石墨烯在拉伸过程中的几何构型图.在整个拉伸变形过程中，同一尺寸单层石墨烯两个方向的破坏形态是不同的;不同尺寸条件下，单层石墨烯在相同方向的拉伸破坏形态保持一致.

图5 模型编号为6 的单层石墨烯在两个方向上的拉伸变形Fig.5 Tensile deformation in both directions of the NO.6 monolayer graphene sheet

在x 方向拉伸变形过程中，当石墨烯处于弹性变形阶段时，模型在轴向载荷作用下作简单的伸长变形，碳环保持为六边形; 随着载荷的增大，边缘六边形碳环开始变得不规则起来，当载荷达到一定值时，石墨烯一侧边缘碳环出现明显的破坏迹象，同时形成一条近似成45°角的斜向破坏线，迅速向石墨烯内部延伸，直至石墨烯被拉断.

在y 方向拉伸变形过程中，线弹性阶段的拉伸变形与x 方向拉伸变形类似;随着载荷的增大，边缘六边形碳环变得不规则起来，当载荷达到一定值时，石墨烯两侧边缘有一部分碳环出现明显的破坏迹象，随着载荷的增加，这种破坏形式越来越明显，并对称地向石墨烯内部延伸，直至石墨烯被拉断，形成图中所示的破坏形态.

图6 为双层石墨烯在拉伸过程中的几何构型图.通过分析发现:多层石墨烯的拉伸破坏形态与单层石墨烯一致;多层石墨烯在拉伸过程中，各单层的变形破坏形态也是一致的.

图6 双层石墨烯在两个方向上的拉伸变形Fig.6 Tensile deformation in both directions of doublelayer graphene sheets

3 结论

文中采用分子动力学的方法，对不同尺寸、不同层数石墨烯的拉伸力学性能进行研究.所得结论如下:(1) 石墨烯在整个拉伸变形过程中，两个方向的应力－应变关系呈现不同的规律，表明石墨烯是各向异性的，且这一特性不受尺寸和层数的影响; (2) 对单层石墨烯，x 方向( 扶手椅型) 的弹性模量在1078.02 GPa附近，y 方向( 锯齿型) 的弹性模量在1041.53 GPa附近，在拉伸线弹性变形阶段，单层石墨烯是各向同性的;对多层石墨烯，在拉伸线弹性变形范围内，各向同性更为明显; (3) 在拉伸载荷作用下，扶手椅型石墨烯的破坏从一侧边缘开始，并沿45°方向向石墨烯内部延伸，锯齿型石墨烯的破坏从两侧边缘对称地向内部延伸，随着C—C 键的断裂，石墨烯被拉断.

［1］Iijima S. Helical microtubules of graphitic carbon ［J］.Nature，1991，354:56-58.

［2］Novoselov K S，Geim A K，Morozov S V，et al. Electric field effect in atomically thin carbon films［J］. Science，2004，306(22) :666-669.

［3］Zhou J，Huang R.Internal lattice relaxation of single-layer graphene under in-plane deformation［J］. Journal of the Mechanics and Physics of Solids，2008，56( 4) : 1609-1623.

［4］辛浩.石墨层/碳纳米管力学性能的研究［D］. 广州:华南理工大学土木与交通学院，2010.

［5］Gao Y W，Hao P. Mechanical properties of monolayer grapheme under tensile and compressive loading ［J］.Physica E，2009，41(8) :1561-1566.

［6］Ni Z H，Bu H，Zou M，et al. Anisotropic mechanical properties of graphene sheets from molecular dynamics［J］.Physica B，2010，405(5) :1301-1306.

［7］韩同伟，贺鹏飞，王建，等. 单层石墨烯薄膜拉伸变形的分子动力学模拟［J］. 新型碳材料，2010，25( 4) :261-266.Han Tong-wei，He Peng-fei，Wang Jian，et al. Molecular dynamics simulation of monolayer graphene under tensile deformation［J］. New Carbon Materials，2010，25( 4) :261-266.

［8］Bernal J D. The structure of graphite［J］. Proc Roy Soc A，1924，106(740) :749-773.

［9］Tersoff J. Modeling solid-state chemistry: interatomic potentials for multicomponent systems ［J］. Phys Rev B，1989，39(8) :5566-5568.

［10］Allen M P，Tildesley D J.Computer simulation of liquids［M］.Oxford:Clarendon Press，1989:385.

［11］Nose S.A unified formulation of the constant temperature molecular dynamics methods［J］.The Journal of Chemical Physics，1984，81(1) :511-519.

［12］Hoover W G. Canonical dynamics: equilibrium phasespace distributions ［J］. Physical Review A，1985，31(3) :1695-1702.

［13］Jin Y，Yuan F G.Simulation of elastic properties of single-walled carbon nanotubes ［J］. Composites Science and Technology，2003，63(11) :1507-1515.

［14］Lee C G，Wei X D，Kysar W J，et al.Measurement of the elastic properties and intrinsic strength of monolayer graphene［J］.Science，2008，321(5887) :385-388.

［15］Ni Z H，Wang H M，Kasim J，et al. Graphene thickness determination using reflection and contrast spectroscopy［J］.Nano Letters，2007，7(9) :2758-2763.

［16］Dresselhaus M S，Dresselhaus G，Eklund P C.Science of fullerenes and carbon nanotubes［M］. San Diego: Academic Press，1996:965.