韦 琦, 黄 君, 黄立新,2*, 吴 宇

(1. 广西大学土木建筑工程学院, 南宁 530004； 2. 广西大学工程防灾与结构安全教育部重点实验室, 南宁 530004)

石墨烯是单层碳原子面材料，其碳原子间相互围成蜂窝形的二维网状结构，被认为是一维结构的碳纳米管和三维结构的石墨片的最基本组成单元。石墨烯是已知的材料中最轻、最薄、最坚硬的材料，在力学、电学和热力学等方面都具有极其优异和独特的性质[1-4]。基于第一性原理[5]、分子动力学[6-8]和结构力学[9-11]等方法得到石墨烯的弹性模量约为1 TPa。石墨烯优异的力学性能使其成为纳米复合材料的理想纳米增强材料。

研究[12-16]表明，石墨烯材料增强的纳米复合材料的力学性能增强效果非常明显。随着石墨烯纳米复合材料市场需求的迅速增长，石墨烯性能参数的研究是必不可少的。石墨烯在不同尺度下往往会出现不同的力学特性，表现出明显的尺寸效应，这是石墨烯纳米复合材料性能增强的影响因素之一。Xu等[17]研究了石墨烯与聚对苯二甲酸乙二醇酯基之间切向界面尺寸产生的力学性能的变化和边界效应。通过原位拉曼光谱测定得出界面剪切应力具有尺寸效应，它的值随石墨烯尺寸的增大而显著减小，石墨烯达到宏观尺寸时趋于定值。在力学性能方面也存在着尺寸效应。Chu等[18]采用了分子动力学模拟的方法，利用单轴拉伸来研究有限尺寸下石墨烯纳米带的力学性能尺寸效应。研究发现，锯齿型石墨烯纳米带的尺寸效应显著而扶手型石墨烯纳米带尺寸效应可以忽略不计。值得注意的是，用有限元法、分子静力学方法、分子动力学方法对石墨烯尺寸效应的研究是无法比较的。最后，结果表明，方形锯齿型石墨烯尺寸增大，弹性模量增大，极限破坏应力和应变减小。对于长方形锯齿型石墨烯，宽度增加，弹性模量增大，极限破坏应力和应变减小。

采用杆件来模拟石墨烯之间碳-碳共价键键，杆件之间的连接为刚性连接。由此得到的杆件系统模拟单层石墨烯结构。现基于杆件系统的有限元方法，研究单层石墨烯的尺寸效应。首先构建8个不同尺寸的有限元计算模型，然后分析单层石墨烯力学性能随尺寸变化的情况。通过数值算例分析单层石墨烯力学性能的尺寸效应以及尺寸效应弱化的现象。

1 石墨烯有限元模型

1.1 石墨烯结构

石墨烯晶体属于复式晶格，单层石墨烯呈片状蜂窝型，如图1所示。两个碳原子组成一个元胞，每个碳原子会与3个相邻的碳原子形成共价键σ键，这是自然界中最强的化学键。研究表明，配位数为3的石墨烯两碳原子之间的距离为1.42110-10m。此外，每个键角为120°，因此几何形状上为正六边形。

图1 单层石墨烯结构示意图Fig.1 Schematic diagram of monolayer graphene structure

1.2 杆件系统的有限元计算模型

将石墨烯的碳碳共价键理解为杆件，碳原子理解为杆件之间的连接点，杆件的变形可用来反映实际共价键的拉伸、压缩以及面内键角的变化。在此基础上，得到杆件系统，如图2所示。用等效的杆件系统描述石墨烯结构，把石墨烯结构的受力行为用杆件系统进行模拟。基于最小势能原理，推导出杆件系统梁单元的刚度矩阵，由此得出有限元求解的基本步骤[19]。在大型有限元软件ABAQUS平台上，对杆件系统进行网格划分，定义梁单元的材料属性和尺寸参数，可以实现对石墨烯结构的力学性能分析。

图2 单层石墨烯有限元模型示意图Fig.2 Schematic diagram of finite element model of monolayer graphene

1.3 等效梁单元参数

在分子力学中，分子的总能量由分子动能及势能叠加而成。分子势能其表达式为

U=Unb+Ur+Uθ+Uφ+Uτ+Uel

(1)

式(1)中：Unb为范德华势；Ur为共价键的拉伸势；Uθ为共价键的键角弯曲势；Uφ为二面角的扭转势；Uτ为平面震动势；Uel为库仑静电势。而在势能计算中，因为非键势的作用相对于其他势能影响较小，因此有时可以忽略其作用，Cornell等[20]得出石墨烯的拉伸势与键角变化势等势能函数。为了将梁单元等效碳-碳连接键，Li等[11]做出假定：用圆形截面梁单元模拟碳-碳键，从而得到梁单元的应变能函数。根据分子势能函数以及结构应变能函数，可以推导出等效梁单元的抗拉刚度EA、抗弯刚度EI以及抗扭刚度GJ。图3为圆形截面梁单元示意图。

F为受力，下标x、y、z为受力方向；M为弯矩， 下标x、y、z为弯矩作用方向图3 圆形截面梁单元示意图Fig.3 Schematic diagram of beam element with circular section

在自定义圆形截面梁时，根据在分子模型下的能量变化的拉伸常数kγ，键角弯曲常数kθ，抗扭常数kτ和经典力学中能量变化的等价关系，由式(2)、式(3)和式(4)确定单元的相关参数,包括弹性模量E，泊松比μ，直径d。

(2)

(3)

(4)

Li等[11]得出参数kγ/2=469 kcal/(mol·Å2),kθ/2=63 kcal/(mol·rad2),kτ/2=20 kcal/(mol·rad2)，其中，1 cal=4.186 J,从而得出弹性模量E，泊松比μ，直径d。Li等指出式(4)经过计算无法得出使用的泊松比，但因为kτ对模型的弹性模量影响较小，因此欧拉梁的μ取0。

2 有限元建模

2.1 自定义梁单元建模

碳原子相对位置变化引起的势能等于空间梁单元变形引起的势能，由此可以导出分子力学的力常数与结构构件单元刚度系数的联系。假设杆件的截面为圆截面，杆件之间为刚性连接，如图4所示。

基于有限元法导出梁单元的抗拉刚度矩阵，抗弯刚度矩阵以及抗扭刚度矩阵。根据式(2)、式(3)和式(4)，从而得出梁单元的材料参数E、μ以及几何参数d，有关参数值如表1所示。在得到一系列的杆件参数后，采用Fortran语言编写有限元软件ABAQUS的圆截面梁单元的UEL子程序,以备后续中导入自定义单元进行建模。

在编辑UEL子程序时，对圆截面梁单元的结点坐标、自由度以及单元刚度进行定义。因为最后是通过有限元软件ABAQUS进行相应的力学模型分析，因此需要完成UEL子程序与ABAQUS语言的对接。

图4 石墨烯有限元建模结点连接Fig.4 Node connection of finite element modeling of graphene

表1 梁单元材料参数和几何参数

2.2 单层石墨烯模型建立

单层石墨烯呈正六边形集合的蜂窝型结构，键角均为120°。通过AutoCAD辅助设计软件建立单层石墨烯模型，设置正六边形的边长为0.142 1，这是单层石墨烯碳-碳键的长度，即0.142 1 nm。

在CAD建立完单层石墨烯的几何模型后，导入ABAQUS有限元计算软件进行相应的前处理。在前处理中，需要对单元进行截面属性和材料属性的定义，而这些工作在UEL子程序已经考虑。在划分网格时，杆件的尺寸设置为0.142 1，这与碳-碳键的尺寸一致。对于载荷以及边界条件，均在inp文件里进行相应的设置，并输入单元编号和节点编号。在ABAQUS中建立作业后，通过编写好的inp文件进行计算得出结果，并在ABAQUS中对计算结果进行后处理。

在建立模型时，假设单层石墨烯(正六边形结构)无团聚，无褶皱。如图5所示，单层石墨烯的长和宽分别为l和h。为了研究单层石墨烯力学性能的尺寸效应，建立有限元计算模型时，考虑模型尺寸的8个工况，表2是8个工况的模型尺寸参数。

图5 单层石墨烯计算模型工况Fig.5 Monolayer graphene computational model

表2 8个工况的尺寸参数

如图6所示，考虑扶手型与锯齿型的情况，单层石墨烯的扶手型与锯齿型的规定是根据石墨烯的边界来确定的。在固定好坐标轴后，垂直于x轴的边界形状为锯齿形状，则研究该x轴方向性质的模型表述为锯齿型石墨烯。垂直于x轴的边界形状为扶手型，则研究该x轴方向性质的模型表述为扶手型石墨烯。扶手型与锯齿型的边界情况并不相同，从而在inp文件中编辑的扶手型与锯齿型的边界条件、位移条件等也不同。

图6 单层石墨烯片的手性定义Fig.6 Definition of chirality of monolayer graphene

在单层石墨烯模型中，碳-碳共价键理解为圆截面梁，碳原子分布在各个结点上。因此忽略碳原子的大小，仅以圆截面梁刚性连接。在对单层石墨烯进行尺寸效应的分析中，将会对不同类型的单层石墨烯从几纳米到几十纳米情况下进行研究，并且石墨烯的厚度取0.34 nm[21]。单层石墨烯模型理论上为正方形薄片，但因其扶手型和锯齿型的尺寸并不相等，且不一定恰好等于整数。因此，未得到整数尺寸下的模型力学性能。在这种情况下，取能完整显示六边形的模型，并取最接近该尺寸的模型来研究。同时考虑石墨烯在锯齿型方向与扶手型方向上力学性质的不同，从而全面地对单层石墨烯的尺寸效应进行分析。

2.3 力学性能分析

在小变形的基础上，通过有限元计算，将研究单层石墨烯的弹性模量、剪切模量和泊松比。在单层石墨烯的边界设置约束，使其位移在某些方向为零。同时在一条边界施加适当的小变形。根据有限元计算求得单层石墨烯的支座反力，从而求得石墨烯的材料参数。

如图7所示，针对锯齿型和扶手型石墨烯，在相应边上施加拉伸小变形。由胡克定律得到应变，即

(5)

应力σ和应变ε分别为

(6)

(7)

式中：Fy为y方向支座反力的合力(锯齿型)；Fx为x方向支座反力的合力(扶手型)；Lx和Ly分别为石墨烯x和y方向的初始长度；t为石墨烯的厚度；ΔLy为石墨烯y方向的变形(锯齿型)；ΔLx为石墨烯x方向的变形(扶手型)。

将式(6)和式(7)代入式(5),求得弹性模量为

(8)

根据横向线应变与纵向线应变的比值求得泊松比为

(9)

图7 单层石墨烯受拉伸变形Fig.7 Monolayer graphene under tensile deformation

如图8所示，针对锯齿型和扶手型石墨烯，在相应边上施加剪切小变形。由剪切胡克定律求得剪切模量为

(10)

剪应力和剪应变分别为

(11)

(12)

式中：Fx为x方向支座反力的合力(锯齿型)；Fy为y方向支座反力的合力(扶手型)；Lx和Ly分别为石墨烯x和y方向的初始长度；t为石墨烯的厚度；ΔLx为石墨烯x方向的变形(锯齿型)；ΔLy为石墨烯y方向的变形(扶手型)。式(11)和式(12)代入式(10),求得剪切模量为

(13)

3 结果与分析

在单层石墨烯的边界设置约束，针对锯齿型和扶手型石墨烯，在相应边上施加小变形。借助于大型有限元计算软件ABAQUS，计算单层石墨烯的变形和支座反力，根据上述公式计算出单层石墨烯的性能参数。

有限元计算采用以下数据和条件：单层石墨烯的厚度t=0.34 nm；单层石墨烯的长和宽尺寸考虑8个工况，如表2所示。图7和图8所示的计算模型中，8个工况的Lx和Ly分别取表2的l和h的尺寸。图7(a)锯齿型石墨烯的y方向的拉伸变形量ΔLy=0.001 nm;
图7(b)扶手型石墨烯的x方向的拉伸变形量ΔLx=0.001 nm；图8(a)锯齿型石墨烯的x方向的剪切变形量ΔLx=0.001 nm;
图8(b)扶手型石墨烯的y方向的剪切变形量ΔLy=0.001 nm。统一采用国际单位，根据公式计算得到的弹性模量和剪切模量的单位为Pa，然后换算成GPa。

3.1 弹性模量

为了考察弹性模量随尺寸变化而变化，定义如下量：

(14)

基于图7所示模型的有限元软件ABAQUS计算，分别得到y方向支座反力的合力Fy(锯齿型)和x方向支座反力的合力Fx(扶手型)，代入式(8)得出弹性模量，其结果如表3所示。

表3 单层石墨烯弹性模量

3.2 泊松比

定义

(15)

基于图7所示模型的有限元软件ABAQUS计算，分别得到锯齿型石墨烯x方向的变形ΔLx和扶手型石墨烯y方向的变形ΔLy，代入式(9)得出泊松比，其结果如表4和图10所示。

表4 单层石墨烯泊松比

图10 单层石墨烯泊松比Fig.10 Poisson’s ratio of monolayer graphene

3.3 剪切模量

定义

(16)

基于图8所示模型的有限元软件ABAQUS计算，分别得到x方向支座反力的合力Fx(锯齿型)和y方向支座反力的合力Fy(扶手型)，代入式(13)得出剪切模量，其结果如表5和图11所示。

表5 单层石墨烯剪切模量

图11 单层石墨烯剪切模量Fig.11 Shear modulus of monolayer graphene

无论是弹性模量、剪切模量，还是泊松比，在小尺度下都会存在显著的尺寸效应，但尺寸大于一定值时，尺寸效应消失，数值将会收敛到一定的值。为了进一步研究单层石墨烯力学性能的尺寸效应，基于图7所示的模型，尺寸从大约1 nm变化到80 nm，根据有限元计算得到弹性模量随尺寸变化的情况，如图12所示。从图12中可以看出，小尺度下尺寸效应十分明显，并且扶手型和锯齿型在小尺寸下的弹性模量相差较大。随着单层石墨烯尺寸逐渐增大，扶手型和锯齿型单层石墨烯的弹性模量逐渐下降并且趋近，最后趋近于同一个值。

图12 单层石墨烯弹性模量变化曲线Fig.12 Modulus of elasticity curve of monolayer graphene

4 结论

在片状蜂窝型的单层石墨烯中，将碳-碳共价键视为杆件，碳原子视为杆件之间的连接点，由此得到杆件系统，其变形可用来反映实际共价键的拉伸、压缩、以及面内键角的变化。基于杆件系统的有限元方法，计算出单层石墨烯的弹性模量、剪切模量和泊松比等力学性能参数。针对有限元计算模型，构建了8个从小到大的尺寸工况，求出单层石墨烯的力学性能参数，然后研究这些参数的尺寸效应。通过有限元计算分析，得到以下结论。

(1)小尺寸(0.923 nm×0.984 nm)的计算模型，单层石墨烯的弹性模量和剪切模量均达到最大值。当计算模型尺寸增大到1.969 0 nm×1.989 4 nm时，弹性模量和剪切模量均明显下降。可见，在小尺寸的情况下，单层石墨烯的力学参数变化幅度大，尺寸效应很显著。

(2)随着计算模型尺寸的增大，单层石墨烯的弹性模量和剪切模量逐渐下降，并且下降幅度逐渐减缓，最终弹性模量和剪切模量分别趋近于相应的值，尺寸效应出现弱化；而单层石墨烯的泊松比的数值迅速趋于某一区间范围之内，并且出现上下波动的变化。

(3)相同尺寸的情况下，扶手型单层石墨烯的弹性模量和剪切模量总是大于锯齿型单层石墨烯的相应量。除了小尺寸(0.923 nm×0.984 nm)的一个工况外，扶手型单层石墨烯泊松比的值大于锯齿型的值。