Elastic in-plane properties of 2D linearized models of graphene

Mechanics of Materials 62 (2013) 60–68

Contents lists available at SciVerse ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Elastic in-plane properties of 2D linearized models of graphene I.E. Berinskii a,b,⇑, F.M. Borodich a a b

School of Engineering, Cardiff University, Cardiff CF24 3AA, UK Institute for Problems in Mechanical Engineering of RAS, St. Petersburg 199178, Bol’shoy pr. V.O. 61, Russia

a r t i c l e

i n f o

Article history: Received 6 August 2012 Received in revised form 5 March 2013 Available online 22 March 2013 Keywords: Graphene Carbon nanomaterials Nanomechanics Molecular mechanics Elastic properties

a b s t r a c t Graphene is a monolayer of carbon atoms packed into a two-dimensional honeycomb lattice. This allotrope can be considered as mother of all graphitic forms of carbon. The elastic in-plane properties of graphene are studied and various existing linearized models of its elastic deformations are critically re-examined. Problems related to modelling of graphene by nonlinear multi-body potentials of interaction are also discussed. It is shown that experimental results for small deformations can be well described by both the two-parametric molecular mechanics model developed by Gillis in 1984, while some popular models have serious flaws and often the results obtained using these models do not have physical meaning. It is argued that in order to study elastic constants of linearized models of graphene layers, it is very convenient to use the four parameter molecular mechanics model. The advantages of this approach is demonstrated by its application to the Tersoff and Brenner nonlinear interaction potentials, and by its comparison with the Gillis two-parametric model. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction As it was noted in Geim and Novoselov (2007), graphene is a basic building block for graphitic materials of all other dimensionalities: it can be wrapped up into 0D buckyballs, rolled up into 1D nanotubes or stacked into 3D graphite. Thus, studies of mechanical properties of graphene are of enormous practical interest. Mechanical properties of graphite were studied for many years (Blakslee, 1970; Bosak et al., 2007). Carbon nanotubes were discovered experimentally in Institute of Physical Chemistry (USSR Academy of Sciences) in about 1950. That time Derjaguin and his co-workers started their intensive studies of various carbon structures, see e.g. their paper that described Derjaguin and Fedoseev results for 25years work in this area (Derjaguin and Fedoseev, 1975). In particular, Derjaguin and his group introduced chemical vapour deposition (CVD) technology. Although,

⇑ Corresponding author at: School of Engineering, Cardiff University, Cardiff CF24 3AA, UK. Tel.: +44 75515 77077. E-mail address: [email protected] (I.E. Berinskii). 0167-6636/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mechmat.2013.03.004

these studies were mainly focused on diamond growth on diamond specimens, Derjaguins group also got as a by-product various carbon whiskers including carbon nanotubes (Derjaguin and Fedoseev D.V., 1977). The electron microscopy studies in the Institute of Physical Chemistry were fulfilled mainly by Radushkevich and Lukyanovich from another group led by Dubinin. In 1952 Radushkevich and Lukyanovich (Radushkevich and Lukyanovich, 1952 presented experimental images of tubular carbon structures that later started to be called carbon nanotubes (CNTs). Derjaguin and Fedoseev called such structures as hollow whiskers or thread-like carbon crystals. Their diffraction studies showed: (i) these whiskers are made from pure carbon and there are no metallic inclusions; (ii) the graphite planes elongated along the whiskers and this can explain their high tensile strength (Derjaguin and Fedoseev D.V., 1977). The possibility of tubular structure of carbon allotropes was rediscovered many times (see e.g. Kornilov, 1985). Nesterenko et al. (1982) in Nesterenko et al. (1982) studied electron microscop images of CNTs, gave a mathematical description of their crystal structure, and presented images of

61

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

several variants of the closure of the monoatomic hexagonal carbon lattices that form cylindrical tube-like crystals. Nowadays these variants are called chirality of CNTs. Rediscovery of CNTs by Iijima (1991), was the trigger for intensive studies of properties of CNTs Dresselhaus et al., 1996. After publishing of Derjaguin and Fedoseev book (Derjaguin and Fedoseev D.V., 1977), it was clear that studies of mechanical properties of graphite and graphite layers are of interst for studies of mechanical properties of CNTs. However, studies of monoatomic structures were considered as auxiliary studies for real life problems because, it was assumed for many years that monoatomic carbon structure does not exist in the free state. Hence, the first mechanical models of graphene appeared as a tool to study properties of in-plane properties of graphite. Although the term ‘‘graphene’’ as a single carbon layer of the graphitic structure first appeared in 1987 Djurado et al. (1987) and became well-established since 1994 Boehm et al., 1994, the first isolated individual graphene planes were obtained by using adhesive tape only in 2004 Novoselov et al. (2004). In fact, Geim and Novoselov and their coworkers showed that graphene can exist in the free state and it is stable (Geim and Novoselov, 2007). Their further studies of graphene were devoted mainly to its remarkable electronic properties. The present paper is devoted to critical re-examination of linearized models of elastic deformations of graphene. 2. Molecular mechanics models of graphene lattices The interest to mechanical properties of honeycomb structures was connected with two main directions: (i) engineering two-dimensional cellular constructions used mainly in aircraft industry; and (ii) models of graphite and its derivatives (graphene, CNT, fullerens so on). 2.1. Structural mechanics approaches to engineering constructions and graphene lattices Honeycomb structures were used in engineering applications to make light and stiff constructions such as aircraft panels or skis. An analysis of in-plane properties of hexagonal structures was presented in El-Sayed et al. (1979). In particular, the effective in-plane elastic moduli and Poisson ratio were calculated for honeycomb structures loaded in two orthogonal directions and then some plastic properties of such constructions were investigated. These structures were also studied in Gibson et al. (1982) and Gibson and Ashby (1997) where linear and nonlinear elastic and plastic properties of honeycombs were considered and compared with experiments. In Masters and Evans (1996) previous models were developed by taking into account the deformation of honeycombs by stretching and hinging in addition to the flexure considered before. All these models based on using the flex-like beams (the Bernoulli–Euler type beams whose neutral fiber can elongate) showed the isotropy of hexagonal honeycomb structure in the case of linear elasticity. The same conclusion follows from the molecular mechanics model of Gillis

(1984). Note that that the conclusion about the isotropy of graphene lattice can be obtained without employment of any structural or spring model but using just the continuum models of graphene Berinskii and Borodich, 2013. In Li and Chou (2003), Tserpes and Papanikos (2005) it was suggested to combine the molecular mechanics (MM) and the structural mechanics approaches in order to model a single-walled carbon nanotube (SWCNT). First the MM approach was applied. For sake of simplicity, the simplest harmonic forms of potentials were adopted and the dihedral angle torsion and the improper torsion were merged into a single equivalent term, i.e.,

U r ¼ 12 kr ðr  r 0 Þ2 ¼ 12 kr ðDrÞ2 ;

U s ¼ 12 ks ðDuÞ2

2

U h ¼ 12 kh ðh  h0 Þ ¼ 12 kh ðDhÞ2 ;

ð1Þ

where kr ; kh and ks are the bond stretching force constant, bond angle bending force constant and torsional resistance respectively, the symbols Dr; Dh and Du represent the bond stretching increment, the bond angle change and the angle change of bond twisting respectively, and the symbols U r ; U h ; U s represent the corresponding harmonic potentials. Then it was noted that in a carbon nanotube, the carbon atoms are bonded to each other by covalent bonds and form hexagons on the wall of the tube. It was suggested to consider the covalent bonds as connecting elements between carbon atoms. According to the theory of classical structural mechanics (Timoshenko, 1955), the strain energy of a uniform beam of length L subjected to pure axial force N is

UA ¼

1 2

Z

N2 1 N2 L 1 EA ¼ ðDLÞ2 dx ¼ 2 EA 2 L EA

L 0

ð2Þ

where DL is the axial stretching deformation, A is a section area, and E is Young modulus. The strain energy of a uniform beam under pure bending moment M is

UM ¼

Z

1 2

L

0

M2 EI 1 EI ð2aÞ2 dx ¼ 2 a2 ¼ L 2 L EI

ð3Þ

where a denotes the rotational angle at the ends of the beam, I is a moment of inertia. The strain energy of a uniform beam under pure torsion T is

UT ¼

1 2

Z

L 0

T2 1 T 2 L 1 GJ ¼ ðDbÞ2 dx ¼ 2 GJ 2 L GJ

ð4Þ

where Db is the relative rotation between the ends of the beam, J is a polar moment of inertia, G is a shear modulus. Comparing the Eqs. (1)–(4), it was suggested the following direct relationship between the structural mechanics parameters EA, EI and GJ and the molecular mechanics parameters kr ; kh and ks

EA ¼ kr ; L

EI ¼ kh ; L

GJ ¼ ks L

ð5Þ

The beam is assumed to have a circular cross-section of radius R. It is known that the initial carbon–carbon bond length is L ¼ 0:142nm. The force constant values were chosen based upon the experience with graphite sheets (Cornell et al.,

62

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68 1

1

1995): kr =2 ¼ 469  kcal  mol  Å2, kh =2 ¼ 63  kcal  mol  2 1 2 rad , and ks =2 ¼ 20  kcal  mol  rad . The similar approaches are often used to describe the elastic properties of graphene and CNTs (Scarpa et al., 2009 and Scarpa et al., 2010). Unfortunately, this model has serious flaws. First of all, a covalent bond is modelled as a beam and the structural mechanics assumes that bended beams have both fibres that are in tension or compression. On the other hand, it is known from physics that a covalent bond involves an overlap of the electron clouds from each atom in order to produce a greater density of electron cloud along the line connecting the two atoms. In order to model the bond as an elastic beam, the electron cloud has to have simultaneously both regions that are in tension and ones in compression. It is rather difficult to imagine this occurs in nature. Second drawback is related to the radius of the ‘‘beam’’. Let us evaluate it for circular beam section with the following moment of inertia, section area, and shear modulus (Timoshenko, 1955)

J¼

pR4 2

;

A ¼ pR2 ;

G¼

E : 2ð1 þ mÞ

ð6Þ

dx ¼ dl sin a ¼

Using this, one can obtain

kr A 4 ¼ ¼ 2: J R kh

ð7Þ

Substituting the above values of the carbon–carbon force constants, one obtains

R2 ¼ 4

63 2 Å 469

ð8Þ

F F sin a cos a dy ¼ dl cos a ¼ cos2 a: ð10Þ c c

Now consider another problem. Let T be a torque applied to the angular spring with stiffness k as it is shown in Fig. 3. As T ¼ kda where da is an twist of the angular spring, then units of k are N  m. Let us connect the angular spring with a rigid rod with a length a. Then one can represent the torque through a force F applied along y direction:

T ¼ F c a ¼ aF cos b

or

R ¼ 0:73Å; D ¼ 2R ¼ 1:46Å:

applied to graphene (Davydov, 2010) in order to derive some estimations of the graphene constants. Apparently, the first two-parametric model (nuclei are connected by two types of linear elastic springs) was presented in 1984 by Gillis (1984) in order to describe the elastic properties of graphite. Namely, the in-plane interaction of carbon atoms were simulated with two independent force constants for interatomic bond stretching and valence angle deformation. Hexagonal ring shown at Fig. 1 simulates a part of graphene lattice. The interatomic bonds are simulated by linear and angular springs having the stiffnesses c and k respectively. The initial length of the hexagon side (carbon interatomic distance) is a ¼ 1:42Å and the initial angle between two any sides is a ¼ 120 . To explain the Gillis model let us consider two auxiliary problems. First problem is devoted to elongation of the linear spring fixed at one end by the external force F applied to its other end along the y direction (Fig. 2). The projection of the force F to the spring line is F c ¼ F cos a and F c ¼ cdl where c is a stiffness of the spring and dl is its elongation. Using this one can obtain the x and y components of displacement:

ð9Þ

Hence, the ‘‘beam’’ aspect ratio L=D  1 and the Bernoulli– Euler beam theory is not applicable for modelling a covalent bond in case of graphene parameters. One can see that the value of the Poisson ratio m does not affect the above arguments. This means that the above structural mechanics model should not be used to simulate bending and shear of graphene and single wall carbon nanotube membranes. 2.2. Linearized molecular mechanics models of graphene lattices Evidently, the graphene layer cannot be modelled solely by one-parameter model with elastic springs representing the central interactions. Indeed, the layer would not satisfy the restrictions of structural rigidity and the model structure would be a mechanism, i.e. under tension the model could have large displacements, while none of the model springs undergoes any elongation. Hence, two parameter models describing central and non-central interaction of atoms are widely accepted. It is difficult to trace back the origin of the two parameter models. In particular, a two parameter model introduced by Keating (1996) was used to describe diamond, silicon and germanium. Recently this model has been

ð11Þ

Hence, the displacement of the edge of the rod could be written as

dx ¼ adb sin b dy ¼ adb cos b:

ð12Þ

Consider a hexagonal ring that simulates a part of graphene lattice constructed with the elements described above. All linear and angular springs have the stiffnesses c and k respectively (Fig. 1) and and simulate the interatiomic bonds. The initial length of the hexagon side is a carbon interatomic distance and the initial angle between two sides is a ¼ 120 . When the ring is elongated (see Fig. 1) the total change of the angles is 4da  2ð2daÞ ¼ 0. It can be shown that in this configuration

da ¼

T : 6k

ð13Þ

To obtain elastic modulus in ‘‘y’’ direction consider the characteristic repeating unit as it is shown at Fig. 4. In response to the loading shown, the lower bond in ‘‘y’’ direction stretches by dy1 ¼ 2F=c According to the (10), due to linear springs elongation the other two bonds stretch as

dy2 ¼

F cos2 60 : c

ð14Þ

At the same time there is an additional elongation of these springs caused by change of the valence angle. Using (11),(13) one can obtain

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

63

Fig. 1. Two-parametric model of hexagon.

Fig. 2. The linear spring under the acting force.

Fig. 4. A repeating unit loading in uniaxial tension in the y direction.

rGillis ¼

Fig. 3. The angular spring under the acting force.

dy3 ¼

Fa2 cos2 30 : 6k

ð15Þ

yy

2F

ð16Þ

Now let us find the corresponding stress component ryy . Gillis determines the stress components in their classical meaning as a force divided on some area. Forpthe ffiffiffi situation considered above the unit area is A ¼ h=2  2 3=2a, where h=2 is a half lattice spacing considered by Gillis as a thickness of graphene layer. For graphite h ¼ 0:68 nm. Hence, one has

ð17Þ

Nowadays there is still a discussion what is better to use as thickness of the graphene layer (its dimension in z direction). Therefore, to avoid this problem, another approach is used. The elastic modulus for graphene and SWCNTs is defined as a product of the conventional Young’s modulus with the layer thickness (see e.g. Zhang et al., 2002; Chang and Gao, 2003; Arroyo and Belytschko, 2004; Berinskii and Krivtsov, 2010; Berinskii et al., 2009). Hence, the elastic modules Ex and Ey will have physical dimension force/ length, i.e. it is represented in the units N/m. For loading in y direction one has

The ‘‘y’’ strain component is

  ðdy1 þ dy2 þ dy3 Þ F 18 a2 : ¼ ¼ þ 3a=2 12a c k

pffiffiffi 2F 4 3F ¼ : A 3ah

ryy ¼ pffiffiffi ¼ 3a

pffiffiffi 2 3F : 3a

ð18Þ

Finally, let us obtain elastic modulus in y direction:

pffiffiffi

Ey ¼

ryy 8 3ck ¼ : yy 18k þ ca2

ð19Þ

The Poisson ratio may be found as

mxy ¼ 

11 dx  ¼ pffiffiffi1 : 22 11 3a

ð20Þ

64

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

Using 10,11, one can obtain

! pffiffiffi pffiffiffi 3 1 F 1 3 Fa2 dx1 ¼ 2  : 2 2 c 2 2 6k

ð21Þ

Further, from 20, 21 follows

mxy ¼

ca2  6k ca2 þ 18k

ð22Þ

Now let us consider another stress state shown at Fig. 5. In this case one need to change (10), (11) to similar formulas for horizontal direction of the applied force to obtain a displacement that corresponds to this stress state. Then one can obtain dx ¼ dx1 þ dx2 where

dx1 ¼ 2

  F Fa2 2 dx2 ¼ 2ada sin 30 ¼ cos2 30 sin 30 c 6k ð23Þ

Using the last formula with determination of the elastic characteristics

rxx ¼

2F 3a

dx

xx ¼ pffiffiffi

3a

¼

pffiffiffi   3 18 a2 : F þ c k 36

ð24Þ

It follows from the above that to obtain the elastic modulus in x direction is

pffiffiffi

Ex ¼

rxx 8 3ck ¼ : xx 18k þ ca2

ð25Þ

One can see that the last formula coincides with (19) that means Ex ¼ Ey . Using the same procedure as described above for lading in y direction, the following expression for the Poisson ratio can be derived

mxy ¼

ca2  6k : ca2 þ 18k

ð26Þ

One can see that it coincides with (22). This corresponds to the symmetry requirements of graphene crystal lattice and corresponds to the in-plane isotropy of graphene mechanical properties. Due to its simplicity such models are widely used to simulate nanostructures. Chang and Gao, 2003 used a similar molecular mechanics approach to find the size-dependent elastic properties of SWCNTs. They calculated

Fig. 5. A repeating unit loading in uniaxial tension in the x direction.

analytically the elastic properties of chiral nanotubes as a function of nanotube diameter. The same approach to determine the elastic properties of SWCNTs as a function of the nanotube size in terms of the chiral vector integers ðn; mÞ, was used in Natsuki et al. (2004). Being unaware of the Gillis model, another two-parametric model of graphene was presented in Berinskiy et al. (2008), Berinskii et al. (2009) the model is formally equivalent to the Gillis model, however, the obtained elastic characteristics were presented as characteristics attributed to the effective thickness of a graphene layer. Two-parametric model allowed to calculate the elastic modules of graphene and compare them with experimental data. However, these models is not applicable to describe large deformations of carbon atoms. Nowadays such effects are investigated with molecular dynamics methods using nonlinear interaction potentials. 2.3. Molecular mechanics models with nonlinear interaction potentials Transformation of a molecular mechanics model (or an atomistic model) to an equivalent continuum mechanics model is a usual approach for modelling of crystall systems, e.g. this idea was applied for studying elastic crystal in Born and Huang (1954) and it is usually referred to as the Cauchy-Born (CB) rule. The CB rule or CB approximation states that the positions of the atoms within the crystal lattice follow the overall strain of the medium for small strains applied to the lattice. The CB rule is used to link the continuum strain energy density and microscopic parameters of the crystal lattice such as its geometry and parameters of the interaction potentials. This approximation generally holds for centrosymmetric crystal systems. However, for complex non-centrosymmetric lattices such as graphene this rule has to be modified to allow for internal degrees of freedom between the sublattices. A model of graphene crystal is shown in Fig. 6. Its lattice is not close-packed, i.e. the average density of atoms location in crystal lattice is not maximum possible, while the triangular crystal lattice is the unique plane close-packed structure with the highest possible average density. One can see that graphene crystal lattice can be decomposed to two triangular centrosymmetrical sub-lattices (see dark and light atoms in Fig. 6). In an elementary cell of the crystal lattice, there are two atoms shown as rhombus. Crystal lattices having more than one atom in elementary cell are called as complex lattices.

Fig. 6. Graphene crystal lattice and an elementary cell with two atoms.

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

65

A feature of both simple and complex non-closepacked lattices, especially connected by covalent bonds, is the impossibility to simulate interaction of the atoms using pair potentials (such as Lennard–Jones). Pair potentials do not contain any possibility to keep the valent angle between the atomic bonds. As a result, the crystal lattice misses the stability and transform into the close–packed structure. One of the possible ways to solve this problem is taking the multi-body interaction into account.

decomposed in Taylor series then one can obtain represent the interaction energy per volume occupied by one atom as quadratic form containing four parameters Gk , k ¼ 1; . . . ; 4

2.3.1. The Tersoff and Brenner interaction potentials In 1988 Tersoff developed the non-linear multi-body potential to describe carbon structures Tersoff, 1988. This potential was developed later by Brenner and co-workers Brenner, 1990; Brenner et al., 2002. Its general form may be represented as

Here atoms and respective bonds are noted by indexes a; b and c; V 0 is a volume of elementary cell, ja jb are deformations of the bonds a and b; nab is a change of angle between bonds. The summation sign labeled with a prime sign means that only adjacent bonds are taken into account. It is assumed that a – b – c takes place for each individual term. This four-parameter model is the linearized four-parameter model with respect to forces. This approach can be applied to any nonlinear potential, in particular to the Tersoff and Brenner potentials (Tersoff, 1988; Brenner, 1990; Brenner et al., 2002). For each particular potential, one can find the appropriate values of Gk . On the other hand, from macroscopic point of view the energy per a unite volume (per unit area in 2D case) may be represented as quadratic form of the strain tensor and the discrepancy vector

Pðrij Þ ¼ V R ðr ij Þ  Bij V A ðr ij Þ

ð27Þ

for atoms i and j, where rij is the distance between the atoms i and j and V R and V A are the repulsive and attractive pair terms (note there is no summation over repeated indices). A coefficient Bij represents the many-body coupling between the bond from atom i to atom j and the local environment of atom i. This coefficient is a function of angle between bonds i  j and i  k. Here index k runs over two nearest bonds so that k – j. Tersoff, Brenner and other potentials (see e.g. Allinger et al., 1989; Cornell et al., 1995; Stuart et al., 2000) were widely used to calculate different physical and chemical properties of graphene, CNT, and diamond. At the same time, it was not clear whether they are applicable to describe the mechanical properties of such objects and how is it possible to connect the parameters of such potentials with elastic characteristics of investigated materials. The main problem to describe the complex lattices is impossibility to use the classical CB rule. However, a nanoscale continuum theory that connects directly the parameters of multi-body interatomic potentials into a continuum model of single-wall nanotube was suggested in Zhang et al. (2002). It was shown that even though the hexagonal lattice of graphene sheet does not possess centrosymmetry, the CB rule can be applied independently to centrosymmetric sub-lattices of graphene (Fig. 6) and and the general deformation consists of the uniform deformation of the sublattices and the rigid translation of the sub-lattices with respect to each other. The elastic moduli of nanotubes were obtained with Tersoff and first-generation Brenner potentials Tersoff, 1988; Brenner, 1990. Arroyo and Belytschko Arroyo and Belytschko, 2004 used a similar approach combined with a finite deformation theory taken into account the inner displacements. They introduced an extention of the Cauchy–Born rule called the exponential Cauchy–Born rule, which accounts for the curvature of the nanotube surface. Similar approaches were used in Zhou and Huang (2008) and Lu and Huang (2009) to investigate internal lattice relaxation and non-linear mechanics of graphene. 2.3.2. A four-parameter model Berinskii and Krivtsov Berinskii and Krivtsov, 2010 noted that if nonlinear interaction potentials are

0 0 X X X 1 G1 j2a þ G2 n2ab þ G3 ðja þ jb Þnab V0 a a;b a;b ! 0 X þG4 nab nac :

W¼

ð28Þ

a;b;c

W¼

¼

2 1 X C ijkl ij kl 2 i;j;k;l¼1 2 2 2 X 1 X 1X C ijkl ij kl þ þ C ij fi fj þ C ijk fi jk 2 i;j;k;l¼1 2 i;j¼1 i;j;k¼1

ð29Þ

From (29) it follows that the elasticity tensor C ijkl is represented by tensor C ijkl connected with the uniform deformation of the sublattices and two additional second and third rank tensors C ij and C ijk . A discrepancy vector f is connected with the inner degree of freedom caused by shift of sublattices. The elastic properties of the graphene lattice may be found from comparison of (28) and (29) as functions of coefficients Gk . In particular, the following elastic modulus and Poisson ratio were found (Berinskii and Krivtsov, 2010):

E¼

m¼

36G1 ð2G1 G2  G1 G4  G23 Þ pffiffiffi 2 V 0 ðG1 þ 18G1 G2  9G1 G4  6G23  2 3G1 G3 Þ pffiffiffi G21  6G1 G2 þ 3G1 G4 þ 6G23  2 3G1 G3 pffiffiffi : G21 þ 18G1 G2  9G1 G4  6G23  2 3G1 G3

ð30Þ

2.3.3. Linearization of the Tersoff and Brenner interaction potentials The parameters of four-parameter model calculated in Berinskii and Krivtsov (2010), are given in Table 1. There following notations have been used: TP, BP-1 and BP-2 are respectively the Tersoff potential, first Brenner potential (Brenner, 1990), second Brenner potential (Brenner et al., 2002). Using Gk parameters in (30) one can obtain the elastic modulus and the Poisson ratio corresponding to Tersoff and Brenner potentials (see Table 2). The obtained values coincided with the elastic characteristics

66

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

Table 1 Parameters of linearized model for different potentials of interaction. Potential

G1 ; N  m

G2 ; N  m

G3 ; N  m

G4 ; N  m

TP BP–1 BP–2

40.568 45.634 43.945

9.2607 1.5905 1.5601

3.2795 3.1089 3.6373

3.7687 0.13979 0.13773

obtained using different approaches (Zhang et al., 2002; Arroyo and Belytschko, 2004; Reddy et al., 2006). It follows from Table 2 that there is a very significant discrepancy between the experimental and theoretical data for both the elastic modulus of graphene and the Poisson ratio. One of the paradoxical results is that the Tersoff potential gives a negative value of Poisson ratio. Hence, the values of the elastic characteristics of the graphene layers estimated by linearization of some popular multi-body potentials, are in disagreement with the experimental values of the characteristics. One can conclude that the above mentioned nonlinear models represented in Table 2 give a rather poor describtion of the mechanical properties of graphene even for small deformations and there is a need in further studies in this field. 2.3.4. Experiments and the two-parameter model values Let us choose coefficients Gk in (28) as

G1 ¼

1 2 ca 2

G2 ¼

1 cc a2 2

G3 ¼ 0 G4 ¼ 0

ð31Þ

In this case, the four–parameter model is reduced to the Gillis model with bonds of length a simulated by linear springs of stiffnesses c and angular springs with stiffnesses k ¼ cc a2 . Then one has

pffiffiffi E¼8 3

ccc c þ 18cc

m¼

c  6cc c þ 18cc

ð32Þ

The same result was obtained in Section 2.2 by the direct use of the two-parameter model (Berinskii et al., 2009). For graphene it gives the following values of the parameters c ¼ 730 N/m, cc ¼ 67 N/m. This corresponds to the experimental data from (Blakslee, 1970), namely to E ¼ 350 N/m and m ¼ 0:17. 2.3.5. The discrete-continuum model of graphene In Odegard et al. (2002), a discrete-continuum approach (also called rod or structural approach) was suggested to

Fig. 7. Truss model. Different rods connect the nearest and far atoms.

describe the interaction between carbon atoms by the so-called truss model. The simples (linear) interaction law was used, when the interatomic bonds are replaced by linear springs with different spring constants. As it shown on the Fig. 7, the nearest neighbour atoms are connected with springs type A of stiffness cA , and far neighbour atoms are connected with springs type B of stiffness cB . In the above paper and in the subsequent literature (Goldstein et al., 2008; Scarpa et al., 2009), these springs are called rods, although the flexural rigidity typical of rods is not taken into account in these models. Krivtsov (2009) has shown that this model gives the following elastic constants

pffiffiffi cA þ 6cB E ¼ 2 3c B cA þ 9cB

m¼

cA þ 3cB : cA þ 9cB

ð33Þ

Using the above formulae and experimental data from (Blakslee, 1970) (see Table 2), one can obtain that cA ¼ 305:9 N/m and cB ¼ 172:7 N/m. This result means that the discrete-continuum model is in disagreement with experimental observations and it has confusing properties. Indeed, analysis of the formula for Poisson ratio shows that in 2D case it may be varied as 1=3 < m < 1 (Krivtsov, 2009). At the same time experimental value m for graphene is 0.17. Hence, this model does not permit to obtain the correct value of Poisson ratio for graphene. Further, the parameter cA denotes not only the connection between two neighbour carbon atoms but also the stiffness of the bond between two sub–lattices of graphene (see Fig. 6). If cA is negative, then even the infinitesimally small perturbation leads to the increasing displacement

Table 2 Comparison of the experimental elastic characteristics of graphene and their estimations obtained by linearization of the Tersoff and Brenner potentials. E, N/m

m

Approach

Source

407 236 236 236 227 243 243 360 350 371

0.158 0.412 0.412 0.416 0.397 0.397 0.249 0.170 0.125

TP, 1988 BP-1, 1990 B-1, 1990 BP-1, 1990 BP-2, 2002 BP-2, 2002 BP-2, 2002 Experiment Experiment Experiment

Berinskii et al., 2010 Berinskii and Krivtsov, 2010 Berinskii et al., 2010 Berinskii and Krivtsov, 2010 Zhang et al., 2002 Zhang et al., 2002 Arroyo et al., 2004 Arroyo and Belytschko, 2004 Reddy et al., 2006 Reddy et al., 2006 Berinskii et al., 2010 Berinskii and Krivtsov, 2010 Arroyo et al., 2004 Arroyo and Belytschko, 2004 Bowman et al., 1958 Bowman and Krumhansl, 1958 Blakslee et al., 1970 Blakslee, 1970 Bosak et al., 2007 Bosak et al., 2007

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

of the sublattices from the equilibrium. This means that crystal lattice with these parameters is unstable. Finally, the springs type B are just a convenient mathematical abstraction and their physical meaning is unclear from chemical point of view. It can be explained that the rod model is simple and it provides a descriptive mechanical analogue for carbon bonds. However, although the described discrete-continuum model is increasingly used, the above arguments show that the model should not be applied to graphene.

67

1992 and references therein), however these results are out the scope of this paper. Acknowledgments The authors acknowledge the financial support of MINILUBES (FP7 Marie Curie ITN network 216011-2) by the European Commission and Russian Foundation for Basic Research (grants 12-01-00521-a and 12-01-31459mol_a). The authors are grateful to Prof. A. M. Krivtsov for valuable discussions.

3. Conclusions The elastic in-plane properties of two-dimensional graphene layer have been studied and the existing linearized models of its elastic deformations have been critically re-examined. It is shown that some very popular models of graphene, e.g. the structural mechanics model and the discretecontinuum model, are physically meaningless and they do not reflect the mechanical properties of graphene lattice. Hence these models should not be used. The two-parameter molecular mechanics model developed by Gillis (1984) has been considered in detail. It has been shown that the model can well describe the available experimental results for small deformations. In addition, this model has very attractive properties due to its simplicity and because its parameters have clear physical meaning. Since the Gillis model can be used only to model the linear behaviour of material, more sophisticated models, in particular, models based on the use of the multi-body interaction potentials, have to be used in non-linear case (large deformations). On the other hand, the values of the elastic characteristics of the graphene layers estimated by linearization of the currently used multi-body potentials (the Tersoff and the first and second Brenner potentials) with generally accepted parameters, are in disagreement with the experimental values of the elastic characteristics. Hence, there is a lack of proper nonlinear models for description of graphene mechanical properties. ‘‘Ab initio’’ methods were not considered in this paper, however they also can be used for determineation of the elastic properties of carbon nanostructures, see e.g. (Kudin et al., 2001; Bichoutskaia et al., 2006). For the above nonlinear potentials that are frequently used to model covalent structures and graphenes, it is shown by an expansion the potentials in Taylor series about an equilibrium point that they can be reduced to four-parameter models (Berinskii and Krivtsov, 2010). It is shown that the Gillis model is a particular case of the four-parameter models. The results presented are of importance mainly in application to the graphene samples in tension. If there are compressive regions then one will observe almost immediately buckling of the structure. Hence, some elements will have the out-of-plane deformations while just in-plane properties of graphenes have been considered. The out-of-plane deformations of honeycombs have been studied by several authors (see e.g. Zhang and Ashby,

References Geim, A.K., Novoselov, K.S., 2007. The rise of graphene. Nature Materials 6 (3), 183–191. Blakslee, O.L., 1970. Elastic constants of compression-annealed pyrolytic graphite. Journal of Applied Physics 41 (8), 3373–3383. Bosak, A., Krisch, M., Mohr, M., et al., 2007. Elasticity of single-crystalline graphite: inelastic x-ray scattering study. Physical Review B 75 (15), 153408. Derjaguin, B.V., Fedoseev, D.V., 1975. The synthesis of diamond at low pressure. Scientific American 233 (5), 102–109. Derjaguin, B.V., Fedoseev, D.V., 1977. Growth of Diamond and Graphite from the Gas Phase. Nauka, Moscow (in Russian). Radushkevich, L.V., Lukyanovich, V.M., 1952. Structure of the carbon produced in the thermal decomposition of carbon monoxide on an iron catalyst. Zhurnal Fizicheskoi Khimii 26, 88–95, in Russian. Kornilov, M.Y., 1985. There is a necessity in the tubular carbon. Chemistry and Life 8, pp. 22–23 (in Russian). Nesterenko, A.M., Kolesnik, N.F., Akhmatov, Yu. S., et al., 1982. Special features of the phase composition and structure of products of NiO and Fe2 O3 reaction with carbon monooxide. Russian Metallurgy (Metally) 3, 14–17. Iijima, S., 1991. Helical microtubules of graphitic carbon. Nature 354, 56– 58. Dresselhaus, S., Dresselhaus, G., Eklund, P.C., 1996. Science of Fullerenes and Carbon Nanotubes. Academic Press, San Diego. Djurado, D., Cousseins, J.-C., Mouras, S., et al., 1987. Synthesis of first stage graphite intercalation compounds with fluorides. Revue de Chimie Minerale 24 (5), 572–582. Boehm, H.P., Setton, R., Stumpp, E., 1994. Nomenclature and terminology of graphite intercalation compounds. Pure and Applied Chemistry 66 (9), 1893–1901. Novoselov, K.S., Geim, A.K., Morozov, S.V., et al., 2004. Electric field effect in atomically thin carbon films. Science 306 (5696), 666–669. Abd El-Sayed, F.K., Jones, R., Burgess, I.W., 1979. A theoretical approach to the deformation of honeycomb based composite materials. Composites 10 (4), 209–214. Gibson, L.J., Ashby, M.F., Schajer, G.S., 1982. The mechanics of twodimensional cellular materials . In: Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, vol. 382 (1782)pp. 25–42 Gibson, L.J., Ashby, M.F., 1997. Cellular Solids: Structure and Properties. Cambridge University Press. Masters, I.G., Evans, K.E., 1996. Models for the elastic deformation of honeycombs. Composite Structures 35 (4), 403–422. Gillis, P., 1984. Calculating the elastic constants of graphite. Carbon 22 (4– 5), 387–391. Berinskii, I.E., Borodich, F.M. 2013. On the isotropic elastic properties of graphene crystal lattice. in: H. Altenbach and N.F. Morozov (eds.). Surface effects in Solid Mechanics. Advanced Structured Mechanics, 30. Springer-Verlag Berlin Heidelberg, 33–42. Li, C.Y., Chou, T.W., 2003. A structural mechanics approach for the analysis of carbon nanotubes. International Journal of Solids and Structures 40 (10), 2487–2499. Tserpes, K.I., Papanikos, P., 2005. Finite element modeling of single-walled carbon nanotubes. Composites part B –Engineering 36 (5), 468–477. Timoshenko, S., 1955. Strength of Materials: Elementary Theory and Problems, v.1. Van Nostrand Company, Princeton, NJ. Cornell, W.D., Cieplak, P., Bayly, C.I., et al., 1995. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society 117 (19), 5179–5197.

68

I.E. Berinskii, F.M. Borodich / Mechanics of Materials 62 (2013) 60–68

Scarpa, F., Adhikari, S., Srikantha Phani, A., 2009. Effective elastic mechanical properties of single layer graphene sheets. Nanotechnology 20 (6), 065709. Scarpa, F., Adhikari, S., Gil, A.J., Remillat, C., 2010. The bending of single layer graphene sheets: the lattice versus continuum approach. Nanotechnology 21 (12), 125702. Keating, P.N., 1996. Effect of invariance requirements on the elastic strain energy of crystals with application to the diamond structure. Physical Review 145, 637–645. Davydov, S., 2010. Elastic properties of graphene: the Keating model. Physics of the Solid State 52, 810–812. Zhang, P., Huang, Y., Geubelle, P.H., 2002. The elastic modulus of singlewall carbon nanotubes: a continuum analysis incorporating interatomic potentials. International Journal of Solids and Structures 39 (13–14), 3893–3906. Chang, T., Gao, H., 2003. Size-dependent elastic properties of a singlewalled carbon nanotube via a molecular mechanics model. Journal of the Mechanics and Physics of Solids 51 (6), 1059–1074. Arroyo, M., Belytschko, T., 2004. Finite crystal elasticity of carbon nanotubes based on the exponential cauchy-born rule. Physical Review B 69, 115415. Berinskii, I.E., Krivtsov, A.M., 2010. On using many-particle interatomic potentials to compute elastic properties of graphene and diamond. Mechanics of Solids 45 (6:815–83), DEC 2010. Berinskii, I.E., Krivtsov, A.M., Kudarova, A.M. 2009. Two-parameric multibody model to describe elastic properties of graphene. In: Advances in Continuum Mechanics: to the 70th anniversary of acad. V.A. Levin. pp. 67–82. Dal’nauka, Vladivistok (in Russian). Natsuki, T., Tantrakarn, K., Endo, M., 2004. Prediction of elastic properties for single-walled carbon nanotubes. Carbon 42 (1), 39–45. Berinskiy, I.E., Krivtsov, A.M., Kudarova, A.M. 2008. Determination of macroscopic characteris-tics for graphene layer using angledepending atomic interactions. In: Proc. of XXXVI Summer School Advanced Problems in Mechanics, St.-Petersburg, Russia., pp. 122– 132. Born, M., Huang, K., 1954. Dynamical theory of crystal lattices. International Series of Monographs on Physics. Clarendon Press. Tersoff, J., 1988. Empirical interatomic potential for carbon, with applications to amorphous-carbon. Physical Review Letters 61 (25), 2879–2882.

Brenner, D., 1990. Empirical potential for hydrocarbons for use in simulating the chemical vapor-deposition of diamond films. Physical Review B 42 (15), 9458–9471. Brenner, D.W., Shenderova, O.A., Harrison, J.A., 2002. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. Journal of Physics–Condensed Matter 14 (4), 783–802. Allinger, N.L., Yuh, Y.H., Lii, J.H., 1989. Molecular mechanics. the MM3 force field for hydrocarbons. 1. Journal of the American Chemical Society 111 (23), 8551–8566. Stuart, S.J., Tutein, A.B., Harrison, J.A., 2000. A reactive potential for hydrocarbons with intermolecular interactions. The Journal of Chemical Physics 112 (14), 6472–6486. Zhou, J., Huang, R., 2008. Internal lattice relaxation of single-layer graphene under in–plane deformation. Journal of the Mechanics and Physics of Solids 56 (4), 1609–1623. Lu, Q., Huang, R., 2009. Nonlinear mechanics of single-atomic-layer graphene sheets. International Journal of Applied Mechanics 1 (3), 443–467. Reddy, C.D., Rajendran, S., Liew, K.M., 2006. Equilibrium configuration and continuum elastic properties of finite sized graphene. Nanotechnology 17 (3), 864–870. Bowman, J.C., Krumhansl, J.A., 1958. The low-temperature specific heat of graphite. Journal of Physics and Chemistry of Solids 6 (4), 367–368. Odegard, G.M., Gates, T.S., Nicholson, L.M., et al., 2002. Equivalentcontinuum modeling of nano-structured materials. Composites Science and Technology 62 (14), 1869–1880. Goldstein, R.V., Osipenko, N.M., Chentsov, A.V., 2008. To determination of the strength of nanodimensional objects. Mechanics of Solids 43 (3), :453–469. Krivtsov, A.M., 2009. Theoretical mechanics. Elastic properties of singleatomic and two-atomic crystals. SPbSTU Pub., St. Petersburg. Kudin, K.N., Scuseria, G.E., Yakobson, B.I., 2001. C 2 F, BN, and C nanoshell elasticity from ab initio computations. Physical Review B 64 (23), 235406. Bichoutskaia, E., Heggie, M.I., Popov, A.M., et al., 2006. Interwall interaction and elastic properties of carbon nanotubes. Physical Review B 73 (4), 045435. Zhang, J., Ashby, M.F., 1992. The out-of-plane properties of honeycombs. International Journal of Mechanical Sciences 34 (6), 475–489.