Theoretical Investigation of the Intrinsic Mechanical Properties of

Article pubs.acs.org/JPCC

Theoretical Investigation of the Intrinsic Mechanical Properties of Single- and Double-Layer Graphene Balázs Hajgató,*,† Songül Güryel,† Yves Dauphin,‡ Jean-Marie Blairon,‡ Hans E. Miltner,‡ Gregory Van Lier,† Frank De Proft,† and Paul Geerlings† †

Free University of Brussels - Vrije Universiteit Brussel (VUB), Research Group General Chemistry (ALGC), Pleinlaan 2, B-1050 Brussels, Belgium ‡ SOLVAY S.A., Innovation Center, Ransbeekstraat 310, B-1120 Brussels, Belgium S Supporting Information *

ABSTRACT: In this study, the Young’s and flexural moduli of single- and double-layer graphene have been theoretically investigated using periodic boundary condition (PBC) density functional theory (DFT) with the PBE, HSE06H, and M06L functionals in conjunction with the 6-31G* and the 3-21G basis sets. The unit-cell size and shape dependence as well as the directional dependencies of the mechanical properties have also been investigated. Additionally, the calculated stretching and flexural strain−stress curves are reported. Finally, finiteelement simulations have been performed so as to find a homogeneous equivalent isotropic transverse material for single-layer graphene, in order to reproduce mechanical behavior in both tensile and bending sollicitations.

■

INTRODUCTION Until 1985, only two ordered forms of carbon were known to scientists, namely, diamond, with its perfect crystal structure, and graphite, also crystalline but black and flaky and not at all transparent. Besides those ordered forms, coal, coke, soot, lampblack, and the many kinds of charcoal were known. The graphite structure determines the properties, since it is made up of sheets of carbon atoms arranged in a hexagonal lattice, like a honeycomb of fused benzene rings, with weak bonding between adjacent sheets. This means that graphite easily forms flakes where the sheets can slide over each other, providing the use of graphite as a lubricant, and resulting in good electrical conductivity. However, it was only in 2007 that researchers in Manchester found a way to mechanically peel single two-dimensional sheets from three-dimensional graphite crystals.1 Graphene is the name given to this flat monolayer of carbon atoms tightly packed into a two-dimensional honeycomb lattice. Since the first experimental analysis,1 graphene has recently gained significant attention. In particular, its excellent mechanical properties are an important advantage for practical applications of graphene.2 These mechanical properties have extensively been investigated, and in particular, the most important elastic property, namely, the Young’s modulus, has been predicted using a range of experimental and theoretical approaches. On the experimental side, by ultrasonic, sonic resonance, and static test methods, Blakslee et al.3 reported a Young’s modulus of 1.06 TPa along the basal plane, for bulk pyrolytic graphite, © 2012 American Chemical Society

which has been highly ordered by annealing under compressive stress. It has to be stressed that this experiment was made decades before the discovery of graphene, and one has to assume that the interaction between graphite layers has no influence on the elastic constant of graphite along the basal plane (C11). If this latter statement is not true, then 1.06 TPa is only an approximation for the Young’s modulus of graphene. Frank et al.4 measured the Young’s modulus to be 0.5 TPa, using an atomic force microscope. The measured material was a stack of graphene sheets (less than five layers), fabricated by mechanical exfoliation from kish graphite, and was suspended over photolithographically defined trenches in silicon dioxide. More recently, by nanoindenting the center of a free-standing monolayer graphene membrane with an atomic force microscope, Lee et al.5 measured the Young’s modulus as 1.0 TPa, assuming the thickness of graphene to be 0.335 nm. Another measurement by Gómez-Navarro,6 using suspended graphene made from graphene oxide was determined the Young’s modulus of graphene to be 0.25 TPa in tip-induced deformation experiments. The Young’s moduli of the different samples were between ∼0.08 and ∼0.68 TPa. The above experimental results are rather different from each other. Since the four different experiments used four different materials (graphite, few-layer graphene, single-layer graphene, and reduced graphene oxide), it is not possible to determine the Received: July 27, 2012 Revised: September 17, 2012 Published: September 26, 2012 22608

dx.doi.org/10.1021/jp307469u | J. Phys. Chem. C 2012, 116, 22608−22618

The Journal of Physical Chemistry C

Article

accuracy of the different experimental techniques. However, evidence was found that the crystal defects of graphene have a possibly larger influence on the measured elastic properties than the experimental method itself.6 Another complication during the comparison of Young’s moduli of graphene obtained in different ways is that one needs to know the thickness of the graphene sample in order to get a Young’s modulus, which is often based on assumptions, since the thickness of single- or few-layer atomic structures is not well-defined. In a series of papers, the thickness of single-walled carbon nanotubes (SWCNs) was determined, based on the tip-end melting temperature of SWCNs, as well as their functional dependence on atomic coordination and bonding energy.7 On the basis of this thickness, the Young’s modulus of SWCN is 2.5585 TPa.7,8 With this determined thickness in combination with Raman studies, the Young’s modulus of graphene was determined to be 2.050 TPa.9,10 In most of the experimental studies and in recent theoretical studiess, the thickness of monolayer graphene is assumed to be equal to the interlayer distance of graphite. In this paper, we also follow this assumption. Many theoretical papers also discuss the Young’s modulus of graphene. The calculated values for Young’s modulus for singlelayer graphene range from 0.762 to 5.189 TPa (see Supporting Information). The methods of investigation are continuum mechanics on finite11−14/infinite systems15−17 combining the tight-binding (a.k.a. extended Hückel) method, bond-potential, or molecular mechanics/dynamics; tight-binding calculations on finite systems;18 analytical or numerical solutions of bond potential of finite19/infinite20−24 systems; molecular mechanics on finite systems;25 cellular material mechanics on finite and infinite systems;26 and density functional theory on finite27/ infinite28,29 systems. Besides calculations on pristine graphene, the Young’s modulus for graphene containing vacancies has also been predicted.30 The Young’s modulus calculations for multiple-layer systems are rather scarce; the only example known up to now is a molecular dynamics simulation of grafold, which is an architecture of folded graphene nanoribbon, for which the Young’s modulus is predicted to be between 0.71 and 0.85 TPa.31 Another intrinsic property of graphene is the bending elasticity of graphene. This property has been investigated less, and the published results are more confusing. In the literature, sometimes the bending modulus of a plate (or bending stiffness of a plate), flexural rigidity (or bending stiffness of a beam), and (three-point beam-bending) flexural modulus is often mixed. In this paper, we clarify these definitions.32,33 The formula for flexural rigidity or bending stiffness of a plate (units in Nm) is Dp =

Eh3 12(1 − ν 2)

Figure 1. Three-point beam bending. F is the force, L is the length, w is the width, d is the thickness, and b is the deflection (see eq 2).

Bernoulli beam theory for beams with rectangular cross section and assuming a third-order polynomial for the curvature Ef =

Ed3w 12

3FL 2wd 2 6bd L2

=

L3F 4wd3b

(3)

where F is the force, L is the length, w is the width, d is the thickness, and b is the deflection. The bending stiffness of graphite is often derived from the phonon spectrum of pyrolytic graphite, amounting to 0.192 aNm,34 and a recent work measures the bending stiffness of double-layer graphene using atomic force microscopy, amounting to 5.688 aNm.35 The methods for investigation are periodic DFT calculations on SWNT,28,36,37 analytic solution of some interatomic potentials,38−42 molecular mechanics,43 continuum mechanics combined with molecular mechanics,44 tight-binding with spherical PBC,45 and molecular dynamics combined with density functional tight binding (DFTB).46 Some of the theoretical results for bending properties of graphene are summarized in the Supporting Information. In this paper, the Young’s and flexural moduli of single- and double-layer graphene have been theoretically investigated using PBC DFT. The connecting stretching and bending strain−stress curves are also reported. The unit-cell size and shape dependence as well as the directional dependencies of the mechanical properties have also been investigated.

■

COMPUTATIONAL METHODS AND THEORY All quantum-chemical calculations have been performed with DFT using the pure PBE,47,48 the short-range hybrid HSE06H49−55 developed for solids (sometimes referred as HSEh1PBE), and the local M06L56 functionals, in conjunction with 3-21G,57 and 6-31G*58 split-valence basis sets using PBC implemented in the Gaussian 0959 program package. To minimalize errors arising from numerical instabilities, the convergence on self consistent field was set to 10−10 a.u. on the energy and 10−8 a.u. on the electron density, ultrafine integration grid (99 radial shells and 590 angular points per radial shells), and tight optimization criteria (15 μa.u. on maximum force, 10 μa.u. on RMS force, 60 μa.u. on predicted maximum displacement, and 40 μa.u. on RMS predicted displacement) were used exclusively unless mentioned othervise. The Young’s modulus (E) was calculated according to

(1)

where E is the Young’s modulus of the plate, h is the thickness of the plate, and ν is the Poisson ratio of the plate. The definition for flexural rigidity or bending stifness of a beam (units in Nm2) is Db =

σ flexural stress = f = flexural strain εf

(2)

where E is the Young’s modulus of the beam, d is the thickness of the beam, and w is the width of the beam. The flexural modulus (units in Pa) can be calculated from a three-point beam-bending test (see Figure 1), derived from the Euler22609


The Journal of Physical Chemistry C E=

σ tensile stress = = ε tensile strain

F A ΔL L0

=

FL0 AΔL

Article

(4)

where F is the force, A is the area where the force is applied, L0 is the original length, and ΔL is the elongation/compression. Since the thickness of atomic-sized plates is not well-defined, the thickness of graphene has been approximated as the graphite interlayer distance (3.35 Å) in all calculations in this paper. In the case of double-layer graphene, the thickness is assumed to be double (6.70 Å) compared to the case of singlelayer graphene. The area was calculated based on the unit cell parameters (the width varies with strain) and the abovementioned uniform graphene thickness. Young’s modulus was calculated using the average force acting on unit cells at 0.5% compression and elongation (applied on the translational vector). Strain−stress curves were calculated between 1% and 20% elongation, using 1% step size. The flexural modulus (Ef) was calculated according to the three-point flexural test, using eq 3. The flexural modulus was calculated using the force action on the carbon atoms in the middle of the unit cell, setting the distance to 0.05 Å perpendicular to the basal plane for the two carbon atoms in the middle of the unit cell and the two carbon atoms on the left extreme of the unit cell. In the case of double-layer graphene, the forces were calculated on the two carbon atoms in the middle of the top layer, while pushing them down 0.05 Å perpendicularly to the basal plane (z direction), keeping the position projected to the bending vector of two carbon atoms in the bottom layer on the left extreme of the unit cell. For the strain−stress curves, the previously mentioned distance was set from 0.1 to 2.0 Å, in 0.1 Å or 0.2 Å step size. During all bending calculations, the unit cell angle was kept to 90°. For the finiteelement calculations, we have used the ABAQUS 6.11 program package from SIMULIA.

Figure 2. Different unit cells for PBC parameter check. Translational Vector (TV) 1 (horizontal) is marked with red and TV2 is marked with green.

and 6 atom (1×1b) unit cells, since the same are expected for the 2 × 2 and 3 × 3 supercell systems. The results are reported in Table 2. Inspecting the data reveals that the calculated values for the Young’s modulus are slightly larger using the shortrange hybrid functional (HSE06H) compared to the pure one (PBE). The calculated Young’s modulus (∼1.05 TPa) for single-layer graphene is in an excellent agreement with the experimental in-plane Young’s modulus (1.06 TPa) of graphite.3 In Figure 3, the strain−stress curves are depicted for only those methods giving the maximum (MAX) and minimum (MIN) stresses, and the average (AVE) and median (MED) of all calculations. The differences between the maximum and minimum stress values are also shown as a percentage (DIF). The parameters were calculated in the same way as the Young’s modulus. This implies that the area calculated takes into account the changes in the unit cell shape (width perpendicular to the vector of elongation), but the thickness of graphene was assumed to be constant (3.35 Å). In the case of higher strain, some calculations are not converged (geometry, wave function, etc.); those data points were left out. The convergence problems were not fixed, since they would be computationally demanding and the results are apparently similar to each other. The number of converged data points (#DATA) is also depicted at each strain value in Figure 3. Furthermore, a linear-regression polynomial fit was employed to the strain−stress curves, and the derivative at zero strain was calculated, which gives the Young’s modulus. Different orders of polynomials were used, from 2 to rank 6. From third- and higher-order polynomial fits, rather similar results were obtained, so only the third-order polynomial fits are reported. If the number of data points were less than 8, then the polynomial fits were not performed, since the numerical noise arising from the uncertainties of our calculations would not be

■

RESULTS AND DISCUSSION The PBC calculations for Young’s modulus were carried out using three different unit cells for both single- and double-layer graphene. For each single-layer unit cell, two additional supercells were calculated, in order to check the default PBC parameters and the uniformity of the graphene. The supercells were designed by repeating the basic unit cells in both translational vector directions two and three times. In the case of all nine unit cells (see Figure 2), a full optimization (using default criteria, 0.00095 maximum and 0.00064 root-meansquare step-size) was performed without any restrictions other than those arising from the PBC. The total energy per carbon atom was always found to be −38.08030 and −38.07453 hartree using HSE06H/6-31G* and PBE/6-31G* levels of theory, respectively, with default PBC parameters. The average bond lengths along with the standard deviations of the optimized structures are listed in Table 1. On the basis of the uniform energies and geometries, we found the default PBC parameters to be suitable, and we also found single-layer graphene to be uniform within the boundaries of the applied accuracy of geometry optimization with the restrictions arising from PBC. Since the PBC parameters are expected to be less of an influence on the results of double-layer graphene compared to the single-layer one, a parameter check was not performed on double-layer graphene. The Young’s moduli were numerically calculated for singlelayer graphene, using only the basic cells 2 (1×1h), 4 (1×1c), 22610



Article

Table 1. Statistical Parameters of Optimized Bond Lengths with Different Unit Cells HSE06H/6-31G*

1×1h

1×1c

1×1b

2×2h

2×2c

2×2b

3×3h

3×3c

3×3b

average st. deviation PBE/6-31G*

1.4193 0.0001 1×1h

1.4193 0.0002 1 ×1c

1.4193 0.0005 1×1b

1.4193 0.0005 2×2h

1.4192 0.0004 2×2c

1.4193 0.0004 2×2b

1.4193 0.0004 3×3h

1.4193 0.0004 3×3c

1.4193 0.0004 3×3b

average st. deviation

1.4292 0.0002

1.4291 0.0003

1.4292 0.0003

1.4292 0.0004

1.4293 0.0005

1.4292 0.0003

1.4292 0.0003

1.4292 0.0003

1.4292 0.0003

smoothed. The Young’s moduli calculated from the polynomial fits are reported in Table 2. The virtually matching values to averages of the expansion and compression Young’s moduli shows the validity of the parameter choice in the numerical derivation. In the case of double-layer graphene, dispersion effects might be more important than in the single-layer case. Therefore, an additional local DFT functional was employed, namely, M06L, because it describes π−π stacking better than PBE (about 1 order of magnitude in energy), keeping the computational cost about the same, which is considerably lower than using HSE06H hybrid functional. The basis set dependence of the results was also checked here, since it was expected to be more important than in the case of single-layer graphene. Similar to the single-layer graphene, three different unit cells were considered (Figure 4), namely, the minimal 4-carbon unit cell (d1×1h), the 8-carbon atom cell (d1×1c), and the 12-carbon atom one (d1×1b). All double-layer unit cells represent AB stacking, since it is more stable than the AA stacking. Young’s moduli were again numerically calculated, similar to the singlelayer case, which implies that the two layers were stretched identically. The results are summarized in Tables 3 and 4. Inspecting the data reveals that the calculated values of Young’s modulus are slightly larger using the short-range hybrid functional (HSE06H) compared to the pure one (PBE), and the M06L local functional gives rather similar results compared to the HSE06H functional. The results are slightly larger in the case of smaller basis set. The calculated Young’s modulus (∼1.09−1.10TPa) for double-layer graphene is again in an excellent agreement with the experimental in-plane Young’s modulus of graphite (1.06 TPa). The predicted accuracy of our results is around 10% (±5%), based on numerical noise,

Table 2. Calculated Young’s Moduli of Single-Layer Graphene 0.5% change in length

E (TPa) expansion

E (TPa) compression

E (TPa) average

E (TPa) extrapolatedb

PBE/6-31G* 1×1h TV1a PBE/6-31G* 1×1h TV2 HSE06H/6-31G* 1×1h TV1 HSE06H/6-31G* 1×1h TV2 PBE/6-31G* 1×1c TV1 PBE/6-31G* 1×1c TV2 HSE06H/6-31G* 1×1c TV1 HSE06H/6-31G* 1×1c TV2 PBE/6-31G* 1×1b TV1 PBE/6-31G* 1×1b TV2 HSE06H/6-31G* 1×1b TV1 HSE06H/6-31G* 1×1b TV2 average minimum maximum

1.01

1.04

1.02

1.02

1.01

1.03

1.02

1.02

1.07

1.10

1.09

1.08

1.07

1.10

1.09

1.08

0.991

1.06

1.03

1.02

1.03

1.07

1.05

1.02

1.08

1.08

1.08

1.08

1.10

1.12

1.11

1.08

0.994

0.995

0.994

1.01

1.00

0.984

0.993

1.03

1.06

1.04

1.05

1.07

1.06

1.04

1.05

1.09

1.04 0.991 1.10

1.06 0.984 1.12

1.05 0.993 1.11

1.05 1.01 1.09

The format of label is “level of theory” “unit-cell” (see Figure 2) “direction of length change” (see Figure 2). bBy fitting a third-order polynomial to the strain−stress curves and calculating the derivative at 0 strain. a

Figure 3. Strain−stress in curves for single-layer graphene, using all methods listed in Table 2. Both the difference between the minimum and maximum in percentage (DIF) and the number of calculation methods converged (#DATA) are displayed on the right axis. 22611



Article

Figure 4. Unit cells for double-layer calculations (AB stacking). Bottom layer atoms are in gold. Translational vector (TV) 1 (horizontal) is marked with red and TV2 is marked with green.

Table 3. Calculated Young’s Moduli of Double-Layer Graphene Using 6-31G* Basis Set 0.5% change in length PBE/6-31G* d1×1h TV1a PBE/6-31G* d1×1h TV2 HSE06H/6-31G* d1×1h TV1 HSE06H/6-31G* d1×1h TV2 M06L/6-31G* d1×1h TV1 M06L/6-31G* d1×1h TV2 PBE/6-31G* d1×1c TV1 PBE/6-31G* d1×1c TV2 HSE06H/6-31G* d1×1c TV1 HSE06H/6-31G* d1×1c TV2 M06L/6-31G* d1×1c TV1 M06L/6-31G* d1×1c TV2 PBE/6-31G* d1×1b TV1 PBE/6-31G* d1×1b TV2 HSE06H/6-31G* d1×1b TV1 HSE06H/6-31G* d1×1b TV2 M06L/6-31G* d1×1b TV1 M06L/6-31G* d1×1b TV2 average min max

E (TPa) expansion

E (TPa) compression

E (TPa) average

1.01

1.05

1.03

c

Table 4. Calculated Young’s Moduli of Double-Layer Graphene Using 3-21G Basis Set


1.02

1.07

1.05

c

1.15

1.14

1.15

1.08

1.10

1.12

1.11

c

1.13

1.16

1.14

1.01

1.10

1.12

1.11

1.00

1.03

1.05

1.04

c

1.03

1.06

1.05

1.01

1.10

1.11

1.11

c

1.10

1.11

1.11

1.07

1.03

1.17

1.10

c

1.05

1.15

1.10

c

1.01

1.06

1.04

1.00

1.01

1.06

1.03

0.997

1.12

1.07

1.09

1.07

1.14

1.05

1.10

1.07

1.01

1.18

1.09

c

1.00

1.18

1.09

1.07

1.06 1.00 1.15

1.11 1.05 1.18

1.09 1.03 1.15

1.04 0.997 1.08

0.5% change in length PBE/3-21G d1×1h TV1a PBE/3-21G d1×1h TV2 HSE06H/3-21G d1×1h TV1 HSE06H/3-21G d1×1h TV2 M06L/3-21G d1×1h TV1 M06L/3-21G d1×1h TV2 PBE/3-21G d1×1c TV1 PBE/3-21G d1×1c TV2 HSE06H/3-21G d1×1c TV1 HSE06H/3-21G d1×1c TV2 M06L/3-21G d1×1c TV1 M06L/3-21G d1×1c TV2 PBE/3-21G d1×1b TV1 PBE/3-21G d1×1b TV2 HSE06H/3-21G d1×1b TV1 HSE06H/3-21G d1×1b TV2 M06L/3-21G d1×1b TV1 M06L/3-21G d1×1b TV2 average min max

E (TPa) expansion

E (TPa) compression

E (TPa) average


1.00

1.08

1.04

1.03

1.03

1.09

1.06

1.01

1.14

1.16

1.15

1.11

1.15

1.17

1.16

c

1.11

1.18

1.15

1.10

1.15

1.22

1.18

1.09

0.950

1.15

1.05

c

0.977

1.12

1.05

1.04

1.08

1.15

1.12

1.09

1.12

1.13

1.12

1.07

1.09

1.13

1.11

c

1.09

1.11

1.10

1.12

1.03

1.06

1.04

0.993

1.02

1.06

1.04

c

1.09

1.14

1.11

c

1.06

1.17

1.11

c

1.06

1.14

1.10

1.06

1.07

1.13

1.10

1.08

1.07 0.950 1.15

1.13 1.06 1.22

1.10 1.04 1.18

1.07 0.993 1.12

a

The format of label is “level of theory” “unit-cell” (see Figure 4) “direction of length change” (see Figure 4). bBy fitting a third-order polynomial to the strain−stress curves, and calculating the derivative at 0 strain. cLess than 8 points converged.

The format of label is “level of theory” “unit-cell” (see Figure 4) “direction of length change” (see Figure 4). bBy fitting a third-order polynomial to the strain−stress curves, and calculating the derivative at 0 strain. cLess than 8 points converged.

deviations from isotropic values, and the differences between the applied levels of theory. Since all three employed DFT

functionals give rather similar results, it is worth it to consider the computational demands for the different levels of theory.

a

22612



Article

Figure 5. Expansion strain−stress curves for double-layer graphene, using all methods listed in Table 3. (6-31G*) Both the difference between the minimum and maximum in percentage (DIF) and the number of calculation methods converged (#DATA) are displayed on the right axis.

Figure 6. Expansion strain−stress curves for double-layer graphene, using all methods listed in Table 4. (3-21G) Both the difference between the minimum and maximum in percentage (DIF) and the number of calculation methods converged (#DATA) are displayed on the right axis.

The pure (local) functionals, by nature, are faster compared to hybrid ones, and should therefore be preferred in terms of resources. In general, it is hard to make any assumption on the speed ratio of pure and hybrid calculations, but in our case, the hybrid calculations usually take about 1.5 to 2.5 more time than the calculations employing pure (local) DFT functionals. Using smaller basis sets also speeds up the calculations. The 3-21G basis set has 9 basis functions per carbon atom, while 6-31G* has 15 basis functions per carbon. The DFT calculations are scaling formally with the fourth power of the number of basis functions; therefore, the theoretical speed ratio between the two basis sets are 94/154 = ∼0.13. Nonetheless, for practical reasons this speed ratio is always larger, so using small basis sets has clear advantages in computational times. On the basis of all this findings, the best method to calculate the Young’s modulus in pristine graphene is at the M06L/3-21G level of theory. As in the case of single-layer graphene, a linear-regression polynomial fit was employed to the strain−stress curves. The Young’s moduli calculated from the polynomial fits are reported in Tables 3 and 4. Again, the matching values show the appropriate choice of our parameters during the numerical derivation. In Figures 5 and 6, the strain−stress curves are again depicted for only those methods giving the maximum and minimum stresses, together with the average and median for all calculations. The differences between the maximum and minimum stress values are also shown in percentage (DIF). The parameters were calculated in the same way as in the case of the Young’s modulus. This implies that the area calculated taking into account the changes in the unit cell shape (width perpendicular to the vector of elongation), but the thickness of both layers of graphene, was assumed to remain constant (6.70

Å). In contrast to the Young’s modulus, the smaller basis set calculations are giving smaller stresses for the same strain. More detailed investigation reveals that the stresses are slightly smaller in the case of lower strains calculated with larger basis sets. Due to the decrease of data points for higher strains, significant conclusions cannot be drawn for higher strains. The number of data points is decreasing with applied strain, due to convergence problems. These convergence problems might be an indication of entering into the plastic regime in the strain−stress curves, since this type of convergence error is typical in our type of calculations if bond-reordering happens and straightforward DFT is not capable of correctly describing both bond-breaking and bond-formation. These bond-reorderings are clearly present in the (irreversible) plastic regime, so our static calculations give questionable results here. The convergence problems are seen to be more severe than in the case of single-layer graphene, which probably is an indication that the elastic regime of multilayer graphene is smaller than for single-layer graphene. This point of view is also supported by the strain−stress curve for the HSE06H/3-21G d1×1h TV2 case (Figure 7). During these calculations, the plastic regime was clearly achieved, the only case in all of our 36 calculated double-layer graphene strain−stress curves. The geometry in the plastic regime is reorganized, and the carbon−carbon bonds in both layers are broken with interlayer bonds being formed. The interlayer distance decreases from 3.15 Å to 1.60 Å, at 6% and 9% elongation, respectively (see Figure 8). This means that single-layer (pristine) graphene is embedded in a medium and freestanding multiple-layer (pristine) graphene will have a smaller elastic regime than freestanding (pristine) single-layer graphene. This is understandable, as one considers that during the elongation of the graphene sheet the carbon−carbon bonds 22613



Article

rearrangements) in the strain−stress curve than in the case of any static calculation. In a dynamic calculation, the change would be smoother, as the carbon atoms might have different (kinetic) energy; therefore, bond reordering would not occur suddenly. Finally, it has to be mentioned that the strain corresponding to the elastic−plastic change is not conclusive in the case of any static calculation using a method that is not capable of describing correctly bond breaking and bond forming. Performing (molecular) dynamic calculations combined with quantum chemical methods capable to describe bond-forming and bond-breaking (for example, multiconfigurational or multireference methods) is clearly out of the scope at the current state of computational chemistry. The PBC calculations for flexural modulus were carried out using five different unit cells in the case of single-layer graphene (Figure 9). On the basis of the findings during the Young’s modulus calculations, only the M06L functional was used. The interlayer distance in the case of double-layer calculations seems to be too short using the 3-21G basis set; therefore, we used only the 6-31G* basis set. It is seen to increase the calculation time, but the interlayer interactions are thought to be much more important in the case of flexural modulus compared to the Young’s modulus. The results of the calculations are summarized in Table 5. We note that in the three-point beam bending test the flexural modulus is derived assuming that the curvature of the beam can be described with a third-order polynomial. In our case, the curvature of graphene is not necessarily the same as in a three-point beam-bending test. The error on the flexural modulus of single-layer graphene is estimated to be about 15% (±7.5%), based on the differences between different calculations. Only the smallest unit cell was used to calculate the strain−stress curves, since the unit-cell size has no significant influence on both the calculated flexural. The flexural strain−stress curve for single-layer graphene is depicted in Figure 10. Similarly to the Young’s modulus of graphene, a third-order polynomial was fit to the strain−stress curve, and the derivative was calculated at zero strain. The bending modulus calculated in this way is 314 GPa, which is again in

Figure 7. Expansion strain−stress curve for double-layer graphene at the HSE06H/3-21G d1×1h TV2 level of theory.

become electron-poor and unstable. One possible way to decrease the electron deficit of the carbon atoms in the case of double- or multiple-layer graphene is to form bonds perpendicular to the elongation, which is clearly more desirable than in-plane bond reorganizations due to the already electrondepleted carbon−carbon bonds. The errors arising from our static calculations might be further decreased with the use of a supercell model, allowing the system to be less periodic. We experienced this elastic−plastic transition in the case of the smaller unit cell (d1×1h), and all the perpendicular carbon− carbon bond breaking and formation has to occur at the same time, because this setup does not allow decreasing the periodic symmetry. In the case of a larger unit cell, we expect that the elastic−plastic transition would occur at lower strain and with a softer transition, due to the lower periodic symmetry. In our calculations, the stresses are clearly not valid in the plastic regime, since the uniform thickness of double-layer graphene is not even valid as an approximation, as a consequence of the interlayer carbon−carbon bonds. Additionally, it has to be mentioned that in (molecular) dynamic calculations one might not expect such an abrupt change (that originated from bond

Figure 8. Top and side views of 6% and 9% elongation for double-layer graphene at the HSE06H/3-21G d1×1h TV2 level of theory. Note that the unit cell is replicated for easier view, and the vector of elongation is horizontal. 22614



Article

single-layer graphene using finite-element simulations. The model for the simulation is in Figure 11. For the simulations, the experimental3 linear elastic constant of graphite were used. Then, we modified only the out-of-plane shear modulus in order to reproduce the same deflection obtained from quantum computation in a bending test. The (out-of-plane) shear modulus of the equivalent homogeneous material of singlelayer graphene was found to be 745 MPa. As a comparison, the same value for graphite is 4 GPa, determined experimentally.3 We have found that deflection in bending was only slightly dependent on the Young’s modulus in the direction perpendicular to the basal plane (out-of-plane) of graphite. Overall, the three-point bending tests that have been used in our quantum mechanical simulations were a rather accurate measurement of shear moduli. From the values obtained in the case of single-layer graphene, it is clear that the unit cell size has no significant influence on the calculated flexural properties. Therefore, strain−stress curves and double-layer flexural modulus were calculated using only the smallest unit cell. The unit cell for the calculation of flexural modulus of double-layer graphene is given in Figure 12. The flexural modulus of double-layer graphene obtained using 0.05 Å deflection (1.34% strain) is 64.8 GPa at M06L/6-31G* level of theory. The flexural strain− stress curve for double-layer graphene is given in Figure 13. The character of the strain−stress curve for double-layer graphene is completely different from that of the single-layer one, since it is not monotonous, and the stress is negative around 5% strain. From analysis of the forces in our calculations, it turns out that the residual forces are not (close to) zero, and some carbon atoms undergo upward forces, while other carbon atoms undergo downward forces (Figure 14). This is not entirely surprising because the carbon atoms in the double-layer graphene are no longer uniform due to the AB stacking. There are two different types of carbons: one has a carbon atom on the top; the other has a middle-of-a-ring on the top. The up- and downward forces would wrinkle the surface; however, we kept the atoms in line, to be able to calculate the bending strain. The residual forces of the nonsupported atoms are decreasing with increasing bending strain; meanwhile, the difference between the forces of the supported atoms (Figure 14) remains about the same order of magnitude during the entire bending strain range of our calculations. Calculating the flexural modulus using numerical differentiation from the 5− 12% bending strain calculations is in the range 95.2−106 GPa. Due to the complications arising in the case of calculation, the flexural modulus of double-layer graphene, the accuracy of the results is significantly worse than that of the rest of our calculations, and we estimated around 20−30% (±10−15%) for

Figure 9. PBC unit cells for calculation of the flexural modulus for single-layer graphene. Translational Vector (TV) 1 (horizontal) is marked with red and TV2 is marked with green.

Table 5. Calculated Flexural Moduli of Single-Layer Graphene Using M06L/6-31G* Level of Theory Using PBC unit cella

strain

Ef (GPa)

Ef (GPa) extrapolatedb

S5 S6 S7 S10 SD5

0.66% 0.46% 0.68% 0.66% 0.66%

315 322 313 344 315

314 c c c c

For unit cell definitions, see Figure 9. By fitting a third-order polynomial to the strain−stress curve, and calculating the derivative at 0 strain. cNot calculated. a

b

excellent agreement with the 315 GPa value obtained from numerical differentiation above. We have calculated a homogeneous equivalent isotropic transverse material for the

Figure 10. Flexural strain−stress curve for single-layer graphene, using M06L/6-31G* level of theory and S5 unit cell (see Figure 9) for the PBC. 22615



Article

Figure 11. Finite element model used to fit a homogeneous equivalent isotropic transverse material reproducing the mechanical behavior of singlelayer graphene.

graphene, no finite element simulations have been made in order to find a homogeneous equivalent material for doublelayer graphene.

■

COMMENT ON THE ORIGIN OF ERRORS In our calculations, the errors of different moduli are arising from different origins. The most apparent errors arise from the applied quantum chemical methods, as they only approximate the full solution of Schrödinger equation. Another source of errors is the numerical differentiation of strain−stress curves in order to obtain the moduli. The moduli are the derivatives of the strain−stress curves at zero strain, and since analytical solutions of this problem are not known in the cases of quantum mechanical methods, we have to perform numerical derivations. This type of error is most apparent when we investigate the difference between the compression and expansion Young’s moduli. In theory, the compression Young’s moduli should be larger than the expansion Young’s moduli, if it was obtained by numerical differentiation, due to the asymmetric nature of the compression/expansion potentials and due to bond-vibration anharmonicities. Numerical errors are typically smaller in the cases of more symmetric compression/expansion potentials and in the cases of less anharmonic bond-vibrations. The parameter choices in our numerical differentiations are sufficient, as the derivatives of polynomial fits and the numerical differentiation of the strain− stress curves are practically the same. Further errors arise from the fact that our calculations were performed practically only with periodic symmetry, e.g., no symmetry constraints were applied within the unit-cell. This type of error is most apparent in the cases of 1×1h and 1×1b unit cells, since the two 120°

Figure 12. PBC unit cell (D5) for calculation of the flexural modulus for double-layer graphene. Translational vector (TV) 1 (horizontal) is marked with red and TV2 is marked with green.

the flexural modulus and at least 30−40% (±15−20%) for the flexural strain−stress curves. It has to be mentioned that negative stresses are sometimes measured in real three-point beam-bending tests as well. This usually happens if the material of the bent beam has a high Poisson ratio, and the beam takes a saddle-point shape. The Poisson ratio of graphite3 is only 0.165, but this value alone does not mean that a “graphene-beam” would not take a saddle-point shape during bending. One has to bear in mind that our PBC approach would not allow the double-layer graphene to take a saddle-point shape, so one should use nonperiodic calculations. All our trials to inspect the flexural modulus in a finite system failed using hydrogenterminated double-layer graphite either due to wave function convergence problems or due to the bonding of edge-carbon atoms. Due to the structural instability of double-layer

Figure 13. Flexural strain−stress curve for double-layer graphene, using M06L/6-31G* level of theory and D5 unit cell (see Figure 12) for the PBC. 22616



Article

Figure 14. Forces of the atoms at 0.1 Å (a) and 0.4 Å (b) deflection. The atoms that would move upward are colored in red, and the atoms that would move downward are colored in green. Atoms in black are atoms that are supported.

graphene. This material is available free of charge via the Internet at http://pubs.acs.org.

compression/expansion directions should be the same due to the hexagonal symmetry of the graphene lattice. Inspecting our data in Tables 2, 3, and 4, we can conclude that our parameter choice of geometry optimization is fine, and we can confirm the (in-plane) isotropy of graphene. Some errors are specific for the calculation of flexural moduli. First, the Euler-Bernoulli equations have been solved assuming a third-order polynomial curvature for rectangular beams, and not to our actual curvatures and cross sections. The errors arising from this issue are still acceptable, as the differences between the calculated flexural moduli using different unit cells are about in the same range as in the case of Young’s moduli calculated in different ways. Second, the wrinkling of the double-layer during the flexural modulus calculations causes further errors.

■

■

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support of SOLVAY, S. A. The authors also wish to acknowledge support from the COST Materials, Physical and Nanosciences (MPNS) Action MP0901: “Designing Novel Materials for Nanodevices - from Theory to Practice (NanoTP)”.

CONCLUSIONS

■

In this work, we have predicted the intrinsic mechanical properties of graphene using periodic boundary condition (PBC) calculations. Using different sized unit cells, containing from 2 to 54 carbon atoms, we have found that graphene is uniform (bond lengths are equal) within the boundaries of the employed geometry-convergence criteria. The Young’s modulus and flexural modulus in our periodic calculations are found to be unit-cell and direction independent within the expected accuracy (±5%). The Young’s modulus of single-layer graphene is found to be 1.05 TPa using both numerical differentiation and third-order polynomial fits on the strain−stress curves. The flexural modulus of single-layer graphene is 315 GPa using numerical differentiation and 314 GPa using the derivative of the fitted third-order polynomial on strain−stress curve at zero stress. Finite-element simulations have been performed to find a homogeneous equivalent isotropic transverse material for single-layer graphene, and the out-of-plane shear modulus is found to be 745 MPa. The Young’s modulus of double-layer graphene is calculated to be 1.10 TPa, and 1.07 TPa using numerical differentiation and third-order polynomial fits on strain−stress curves, respectively. The flexural modulus of double-layer graphene is 64.8 GPa, using numerical differentiation. During the bending of double-layer graphene, we have found an unexpected wrinkling of both layers. Furthermore, we found that the strain−stress curves can be well approximated using a third-order polynomial in the investigated range, except in the case of flexural strain−stress curves of double-layer graphene.

■

AUTHOR INFORMATION

REFERENCES

(1) Geim, A. K.; Novoselov, K. S. Nat. Mater. 2007, 6, 183−191. (2) Van Noorden, R. Nature 2006, 442, 228−229. (3) Blakslee, O. L.; Proctor, D. G.; Seldin, E. J.; Spence, G. B.; Weng, T. J. Appl. Phys. 1970, 41, 3373−3382. (4) Frank, I. W.; Tanenbaum, D. M.; van der Zande, A. M.; McEuen, P. L. J. Vac. Sci. Technol., B 2007, 25, 2558−2561. (5) Lee, C.; Wei, X.; Kysar, J. W.; Hone, J. Science 2008, 321, 385− 388. (6) Gómez-Navarro, C.; Burghard, M.; Kern, K. Nano Lett. 2008, 8, 2045−2049. (7) Sun, C. Q.; Bai, H. L.; Tay, B. K.; Li, S.; Jiang, E. Y. J. Phys. Chem. B 2003, 107, 7544−7546. (8) Zheng, W.-T.; Sun, C. Q. Energy Environ. Sci. 2011, 4, 627−655. (9) Sun, C. Q.; Nie, Y.; Pan, J.; Zhang, X.; Ma, S. Z.; Wang, Y.; Zheng, W. RSC Advances 2012, 2, 2377−2383. (10) Yang, X. X.; Li, J. W.; Zhou, Z. F.; Wang, Y.; Yang, L. W.; Zheng, W. T.; Sun, C. Q. Nanoscale 2012, 4, 502−510. (11) Natsuki, T.; Tantrakarn, K.; Endo, M. Carbon 2004, 42, 39−45. (12) Reddy, C. D.; Rajendran, S.; Liew, K. M. Int. J. Nanosci. 2005, 4, 631−636. (13) Reddy, C. D.; Rajendran, S.; Liew, K. M. Nanotechnology 2006, 17, 864−870. (14) Meo, M.; Rossi, M. Compos. Sci. Technol. 2006, 66, 1597−1605. (15) Wu, Y.; Zhang, X.; Leung, A. Y. T.; Zhong, W. Thin-Walled Structures 2006, 44, 667−676. (16) Hemmasizadeh, A.; Mahzoon, M.; Hadi, E.; Khandan, R. Thin Solid Films 2008, 516, 7636−7640. (17) Shokrieh, M. M.; Rafiee, R. Mater. Des. 2010, 31, 790−795. (18) Hernández, E.; Goze, C.; Bernier, P.; Rubio, A. Phys. Rev. Lett. 1998, 80, 4502−4505. (19) Zhao, H.; Min, K.; Aluru, N. R. Nano Lett. 2009, 9, 3012−3015. (20) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. J. Phys.: Condens. Matter 2002, 14, 783−802. (21) Gupta, S.; Dharamvir, K.; Jindal, V. K. Phys. Rev. B 2005, 72, 165428.

ASSOCIATED CONTENT

S Supporting Information *

A short review of theoretically obtained mechanical properties (Young’s, flexural, and bending moduli) of single and few layer 22617



Article

(22) Xiao, J. R.; Gama, B. A.; Gillespie, J. W., Jr. Int. J. Sol. Struct. 2005, 42, 3075−3092. (23) Huang, Y.; Wu, J.; Hwang, K. C. Phys. Rev. B 2006, 74, 245413. (24) Terdalkar, S. S.; Huang, S.; Yuan, H.; Rencis, J. J.; Zhu, T.; Zhang, S. Chem. Phys. Lett. 2010, 494, 218−222. (25) Li, C.; Chou, T.-W. Int. J. Sol. Struct. 2003, 40, 2487−2499. (26) Scarpa, F.; Adhikari, S.; Srikantha Phani, A. Nanotechnology 2009, 20, 065709. (27) Van Lier, G.; Van Alsenoy, C.; Van Doren, V.; Geerlings, P. Chem. Phys. Lett. 2000, 326, 181−185. (28) Kudin, K. N.; Scuseria, G. E.; Yakobson, B. I. Phys. Rev. B 2001, 64, 235406. (29) Liu, F.; Ming, P.; Li, J. Phys. Rev. B 2007, 76, 064120. (30) Güryel, S.; Hajgató, B.; Dauphin, Y.; Blairon, J.-M.; Miltner, H. E.; De Proft, F.; Geerlings, P.; Van Lier, G. Effect of Structural Defects and Chemical Functionalisation on the Intrinsic Mechanical Properties of Graphene. Submitted, 2012. (31) Zheng, Y.; Wei, N.; Fan, Z.; Xu, L.; Huang, Z. Nanotechnology 2011, 22, 405701. (32) Landau, L. D.; Lifchitz, E. M. Theory of Elasticity, 2nd (revised) ed.; Pergamon Press: Oxford, 1970. (33) Young, W. C. ROARK’S Formulas for Stress & Strain, 6th ed.; McGraw Hill International: New York, 1989. (34) Nicklow, R.; Wakabayashi, N.; Smith, H. G. Phys. Rev. B 1972, 5, 4951−4962. (35) Lindahl, N.; Midtvedt, D.; Svensson, J.; Nerushev, O. A.; Lindvall, N.; Isacsson, A.; Campbell, E. E. B. Nano Lett. 2012, 12, 3526−3531. (36) Kürti, J.; Kresse, G.; Kuzmany, H. Phys. Rev. B: Condens. Matter 1998, 58, R8869−R8872. (37) Sánchez-Portal, D.; Artacho, E.; Soler, J. M. Phys. Rev. B 1999, 59, 12678−12688. (38) Tu, Z.-c.; Ou-Yang, Z.-c. Phys. Rev. B 2002, 65, 233407. (39) Arroyo, M.; Belytschko, T. Phys. Rev. B 2004, 69, 115415. (40) Behfar, K.; Seifi, P.; Naghdabadi, R.; Ghanbari, J. Thin Solid Films 2006, 496, 475−480. (41) Lu, Q.; Arroyo, M.; Huang, R. J. Phys. D: Appl. Phys. 2009, 42, 102002. (42) Cranford, S.; Sen, D.; Buehler, M. J. Appl. Phys. Lett. 2009, 95, 123121. (43) Duan, W. H.; Wang, C. M. Nanotechnology 2009, 20, 075702. (44) Lu, Q.; Huang, R. Int. J. Appl. Mech. 2009, 1, 443−467. (45) Koskinen, P.; Kit, O. O. Phys. Rev. B 2010, 82, 235420. (46) Zhang, D.-B.; E., A.; Dumitrică, T. Phys. Rev. Lett. 2011, 106, 255503. (47) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (48) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (49) Henderson, T. M.; Izmaylov, A. F.; Scalmani, G.; Scuseria, G. E. J. Chem. Phys. 2009, 131, 044108. (50) Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. J. Phys. Chem. 2005, 123, 174101. (51) Heyd, J.; Scuseria, G. E. J. Chem. Phys. 2004, 121, 1187−1192. (52) Heyd, J.; Scuseria, G. E. J. Chem. Phys. 2004, 120, 7274−7281. (53) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. J. Chem. Phys. 2006, 124, 2199906. (54) Izmaylov, A. F.; Scuseria, G. E.; Frisch, M. J. J. Chem. Phys. 2006, 125, 104103. (55) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 224106. (56) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101. (57) Gordon, M. S.; Binkley, J. S.; Pople, J. A.; Pietro, W. J.; Hehre, W. J. J. Am. Chem. Soc. 1980, 102, 939−947. (58) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654−3666. (59) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci,

B.; Petersson, G. A. et al., Gaussian 09, rev B.01; Gaussian, Inc.; Wallingford, CT, 2009.

22618


Theoretical Investigation of the Intrinsic Mechanical Properties of

Recommend Documents