Graphene as a hexagonal 2-lattice: evaluation of the in-plane material constants for the linear theory. A multiscale approach
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n This document is the Accepted Manuscript version of a Published Work that appearedin final form in Journal of Applied Physics, copyright American Institute of Physics afterpeer review and technical editing by the publisher. To access the final edited and publishedwork seehttp://aip.scitation.org/doi/full/10.1063/1.49284641 raphene as a hexagonal 2-lattice: evaluation of thein-plane material constants for the linear theory. Amultiscale approach.
D. Sfyris, E.N. Koukaras, N. Pugno, C. GaliotisOctober 17, 2018
Abstract
Continuum modeling of free-standing graphene monolayer, viewed as a two di-mensional 2-lattice, requires specification of the components of the shift vector thatacts as an auxiliary variable. If only in-plane motions are considered the energydepends on an in-plane strain measure and the shift vector. The assumption of ge-ometrical and material linearity leads to quadratic energy terms with respect to theshift vector, the strain tensor, and their combinations. Graphene’s hexagonal sym-metry reduces the number of independent moduli then to four. We evaluate thesefour material parameters using molecular calculations and the AIREBO potentialand compare them with standard linear elastic constitutive modeling. The resultsof our calculations show that the predicted values are in reasonable agreement withthose obtained solely from our molecular calculations as well as those from litera-ture. To the best of our knowledge, this is the first attempt to measure mechanicalproperties when graphene is modeled as a hexagonal 2-lattice. This work targets atthe continuum scale when the insight measurements comes from finer scales usingatomistic simulations.
Keywords: graphene; hexagonal 2-lattice; molecular dynamics; AIREBO potential;material modulus; linear elasticity.
Ever since its discovery ([8]) graphene attracted significant attention in the mechanics lit-erature. Many works are devoted on evaluating graphene’s Young’s modulus and Poissonratio, either by experimental or computational means. Lee et al ([14]) conduct nanoiden-tation measurements using an atomic force microscope and measure Young modulus, E, of340 ±
40 N m − , or of 1 ± − (0.93 TPa) and Poissonratio ν = 0 .
31. Other tight binding claculations ([10]) report an E of 1.21 TPa.Zhou and Huang ([31]) utilize molecular dynamics and employ the Tersoff–Brennerpotential to evaluate E=235 N m − (0.70 TPa), and ν = 0 . ν =0.147, while Reddy et al ([23]) use the Tersoff–Brenner potential to arrive at E=0.67TPa, ν = 0 . − (1.15 TPa), and ν =0 . − (1.02 TPa), and ν =0 . − (1.04TPA), and ν = 0 . − (0.95 TPa).A density functional theory ([12]) render E=1.24 TPa.Arroyo and Belytschko ([2]) use a finite deformation continuum theory derived frominteratomic potentials to derive E=235 N m − (0.70 TPa), and ν = 0 . ν = 0 .
45. The braced truss model ([24]) using theAMBER force field result at E=1.22 TPa, while when the Morse force field is used theyrender E=1.91 TPa. In a recent review paper ([7]) we summarize the relevant literatureon graphene mechanics as probed by deformation and spectroscopic measurements and ascalculated by ab-initio, molecular simulations and continuum mechanics methods.The present work is a continuation of our previous efforts ([25, 26, 27, 28]) to properlymodel graphene at the continuum level. The starting point is the modeling of graphene asa hexagonal 2-lattice ([6]), in line with well established theories of multilattices ([4, 17, 18,19, 20]). By making appropriate hypothesis (see [25]) one works with an energy dependingon an in-plane strain measure, the curvature tensor and the shift vector. Graphene’ssymmetry is taken into account, in this framework, by adding the structural tensor tothe list of independent variables of the energy. This way the complete and irreduciblerepresentation of the energy is evaluated and from it, the stress tensor, the couple stresstensor as well as the driving force related with the shift vector.Simple closed form solutions for this genuinely geometrically and materially nonlineartheory are reported in [26]. The geometrically and materially linear counterpart of theabove theory is given in [27]. There, graphene’s energy is assumed to have quadraticdependence on the in-plane strain measure, the curvature tensor, the shift vector as well asto quadratic combinations of them. Hexagonal symmetry reduces then the overall numberof moduli to nine. If in-plane motions are considered, only four material parameters shouldbe determined; these are the constants c , c , c , c in the terminology of [27].The present work is concerned with the evaluation of these four material parametersusing molecular mechanics with the AIREBO potential. So, while the overall theory appliesto the continuum scale, the calculations come from atomistic insights in finer scales. We3orrespond to the material parameters at the continuum level, four well defined measuresfrom molecular considerations. The strategy for doing that is non-standard and goes asfollows: we start at the discrete level where we focus on graphene’s unit cell and distinguishbetween the measured length of the shift vector and the length of the lattice vectors. Then,we apply a tensile strain up to 6% along the armchair direction for a graphene monolayerthat contains 31600 carbon atoms. Then, we evaluate the radial distribution diagramdescribing length change due to loading for carbon–carbon connections.At zero strain level, we find two peaks on the radial distribution diagrams: one cor-responding to the equal length of the lattice vectors (approximately 0,242 nm), while theother peak corresponds to the shift vector (approximately 0,140 nm). As strain is gradu-ally applied, we find that these peaks split into 2 new peaks each. These peaks measurechanges that happen to graphene’s unit cell due to applied strain. To these four peakswe correspond, at the continuum level, the four required material parameters c , c , c , c .To do this we first define four strain measures as the differences between the length atthe peak point for strain level of 6%, minus the initial length corresponding to the peakat zero strain, divided by the initial length. We plot the applied stress versus these fournewly defined strain quantities. Slopes of these four diagrams correspond to the materialparameters c , c , c , c .As a minimum validation/calibration of our approach we compare with reported valuesof E, and ν from the literature. To obtain the relation of ( E, ν ) with ( c , c , c , c ) we solvethe equations ruling the shift vector, to express the components of the shift vector as afunction of the strain components. Since the problem is geometrically and materially linearthis is feasible; for the nonlinear case this would have been cumbersome, if not non-solvableexplicitly. Having these expressions at hand, we can invert the stress–strain relations toobtain the required expression of ( E, ν ) as function of ( c , c , c , c ). Having ( c , c , c , c )evaluated from molecular calculations, we can then evaluate ( E, ν )=(1.37 TPa, 0.41) forour framework.Values for E and ν from the reported literature (see the first four paragraphs of thisSection) range as E=0.67-1.91 TPa and ν =0.14-0.45 depending on the methodology used.The central tendency of these values for E is the value 1 ± E, ν )=(1.37 TPa, 0.41) overestimates these quantities but stillremain within the range of acceptable values. From the literature cited, the continuummethods (i.e. the finite element approaches of [2, 22, 24]) tend to have greater discrepancyfrom the value 1 ± E, ν )=(0.95TPa, 0.20). These values are obtained using the definition of the AIREBo manual. But, thevalues (
E, ν )=(1.37 TPa, 0.41) are based on a different definition of (
E, ν ): they are basedon a genuinely continuous definition which is non-standard since it uses c , c , c , c . Cer-tainly, the two definitions (the discrete and the continous one) measure the same quantitiesin a different way. So, the discrepancy in their reported values is based on the differentdefinition but still remains in the range of an admissible difference.The paper is organized as follows. Section 2 presents compactly the key findings of4ur previous works ([25, 26, 27]), to which we refer for more information. In Section 3 wepresent the core of our caclulations. We lay down the strategy for obtaining/defining therequired material parameters at the continuum level starting from discrete pictures andmeasurements using the AIREBO potential. Section 4 gives the minimum validation bycorrelating with standard results. The article ends up in Section 5 with some concludingremarks as well as future directions. As far as notation is concerned, we use tensor notationin component form throughout the paper. All indices are assumed to refer to the sameCartesian coordinate system and range from 1 to 2. We begin by presenting the main findings of our relevant previous works ([25, 26, 27])that constitute the theoretical backbone of this work. Some more detailed informationregarding the continuum modeling of crystalline materials can be found in the excellentbook by Pitteri and Zanzotto ([20]).Graphene is modeled as a hexagonal 2-lattice ([6]) at the crystalline level. One arrivesat the continuum level by assuming validity of the Cauchy–Born rule ([5]). Validity of thisrule, together with confinement to weak transformation neighborhoods ([18, 19]) enablesone to work with the symmetry groups classical elasticity uses; for the case of graphenethese are rotations by 60 .For the geometrically and materially linear case energy depends on the in-plane straintensor e αβ = 12 ( u α,β + u β,α ) , (1) u being the in-plane displacement, the curvature tensor b as well as the shift vector p .Explicitly, energy has the form ([27]) W ( e , b , p ) = 12 C ijkl e ij e kl + 12 C ij p i p j + 12 C ijk e ij p k + 12 C ijkl b ij b kl + 12 C ijkl e ij b kl + 12 C ijk b ij p k . (2)Tensors C , ..., C are tensors of material moduli. When out-of-plane motions are ne-glected, terms related with the curvature should be set equal to zero. In this case thestress–strain relations read ([27]) σ = c e + c e − c p , (3) σ = c e + c e + c p , (4) σ = c − c e − c p , (5)steming from the expression for the stress tensor σ = ∂W∂ e = C ijkl e kl + C ijk p k . (6)5or the components related with the shift vector we have ∂W∂p i = C ij p j + C ijk e jk . (7)So, we finally take ∂W∂p = c p − c e , (8) ∂W∂p = c p − c e + c e . (9)The field equations for such a model are the momentum equation in the absence ofbody forces and inertial terms σ ij,j = 0 , (10)and the equation ruling the shift vector ∂W∂p i = 0 . (11)From the physical point of view the momentum equation is the force balance for the surface.The equation ruling the shift vector express that the shift vector adjusts so as equilibriumis reached.So, for the above framework we should evaluate material parameters ( c , c , c , c ).These are material moduli at the continuum level which are present in the constitutive law,thereby characterizing graphene’s mechanical properties in the small strain regime. Nextsection describes how these material parameters can be obtained/defined using molecularcalculations. c , c , c , c from molecular pictures For our purposes we load a graphene sheet along the armchair direction as Figure 1 shows.In this figure we depict the unit cell at ease, as well as the loading direction of the uni-axially strain which applies. Within the unit cell, we differentiate with numbers 1, 2, 3,4 pair of carbon–carbon distances. At ease, pair of distances no. 1 have equal length ofapproximately 0.242 nm. The same holds true for the pair of distances denoted by no. 2.The pair of distances denoted by no. 4 have equal lengths of 0.140 nm, approximately.The same length is shared by distance denoted by no. 3 in Figure 1. Inspecting the unitcell of graphene when it is modeled as a 2-lattice (see figures in [25, 26, 27]) it appearsthat lenghts no. 1, 2 correspond to the lattice vectors of graphene. Lengths numbered asno. 3, 4 correspond to graphene’s shift vector.Now we apply an axial strain along the direction shown in Figure 1. This is done byemploying the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) ([29])6igure 1: Unit cell and loading direction for graphene. Numbers 1 and 2 pertain to pairsof lattice vectors having equal length changes due to loading, while no. 4 pertains to apair of the shift vector components having the same length change due to loading. No. 3is another component of the shift vector which has different length change from no. 4.potential to describe the carbon–carbon interatomic forces. An orthogonal periodic com-putation cell is used with dimension 42.6 nm × x -axis along the armchair direction (the direction of loading). The computationalcell is initially relaxed, leading to an equilibrium structure for the given potential. Thein-plane symmetry of the structure is broken by assigning randomized velocities with aGaussian distribution to all of the atoms corresponding to a temperature of T = 40 K.The large number of atoms considered in the computational cell is needed to properlycapture the distribution of the nearest and next nearest neighbor distances, in what fol-lows. An energy equilibriation is performed within the microcanonical ensemble (NVE).The structure and unit cell are further relaxed by a follow up equilibriation within theisothermal–isobaric ensemble (NPT) at the same temperature. Uniaxial tensile strain ap-plies then by a deformation control method. The strain applies every 200 time steps on the x -axis homogeneously with a strain rate of 0.0005 ps − . The strain rate is very small so wedisregard viscous response, i.e. the system is given ample time to respond to the applieddeformations. The Poisson’s effect is accounted for by allowing the y -axis to relax duringthe slow elongation process. All of the MD simulations are performed using the LAMMPS([21]) software package.As an outcome of this loading process, lenghts denoted by no. 1, 2, 3, 4 within the unitcell change. This change will not be same for all of them, due to their different position7ith respect to the loading direction. Pair of distances denoted by no. 1 tend to shortensince they are perpendicular to the loading direction; they are also affected equally due toloading, this is why we group them together. Pair of distances denoted by no. 2 tend toextent since they are inclined with respect to the extension direction. Certainly, due tosymmetry and homogeneity of the applied extension, lengths denoted by no. 2 experiencethe same change. Pair of distances denoted by no. 4 tend to elongate by the same amountbetween them. For the distance denoted by no. 3 one expects elongation as well, since itis parallel to the direction of loading.Now, the idea is to correlate the continuum material parameters c , c , c , c with dis-tances denoted by no. 1, 2, 3, and 4 in the unit cell. This can be done using the followingprocedure: at zero strain level, we plot the radial distribution diagram (see Figure 2). On
0% strain6% strain R ad i a l D i s t r i bu t i on ( a r b i t r a r y un i t s ) R C − C ( Å ) Figure 2: The radial distribution diagram at ease (above) and at strain level 6% (below).The horizontal axis measures carbon–carbon distances in Angstrom.the horizontal axis of this diagram we have carbon–carbon distances measured in Angstrom.The radial distribution (or pair correlation) diagram describes the probability of encoun-tering a carbon atom at any given distance from another carbon atom. To produce theradial distribution diagrams we extract atomic configuration at regular time step intervals.These are the same configuration that we later use for the calculation of the Poisson’s ra-tio and Young modulus from the molecular modeling itself. For each atom we identify thefirst and second neighbors and calculate the corresponding distances. Thermal fluctuationsalter these distances in a canonical (random) manner around a central value. Each of thedistances is accounted for uniquely. In Figure 2 the formation of distinct density bandscan be seen, one near the radial distance of 0.145 nm and another around 0.25 nm. Thesebands correspond to first and second nearest neighbors. At ease these band are Gaussianpeaks, one centered at 0.242 nm and the second at 0.140 nm. These correspond to lengthsof distances denoted by no. 1 and 2 and no. 3 and 4 at Figure 1, respectively. We remindhere that at ease, pairs no. 1 and 2 have equal length of 0.242 nm. Also, pairs no. 3 and4 have an equal length of 0.140 nm at ease. 8s loading applies these peaks split into two new peaks each (see Figure 2, bottom).The bottom plot of Figure 2 corresponds to the radial distribution at strain level 6%. Thesplitting of each peak into two new ones is apparent. Essentially, these peaks describe thebehavior of lengths no. 1–4 upon loading. Inspecting Figure 2 we see that at 6% the peakcentered at 0.242 splits into two new peaks: one at 0.240 nm and a second at 0.252 nm.These most probable values can be found by fitting two Gaussian functions to the band.As noted earlier, fitting with Gaussian functions is most appropriate due to the stochasticnature of the bond length variation which is a direct result of thermal movement. Eachof the Gaussian functions correspond to specific lengths no. 1 and no. 2, and the centers(positions) of which provide the length values. Similarly, the peak centered at 0.140 nmat ease, upon strain splits into two peaks, one at 0.142 nm and one at 0.146 nm. Overall,this provides the evolution of lengths no. 3–4 upon application of strain.Now, using these peaks we want to define the material parameters c , c , c , c . Werepeat that the radial distribution diagram render the most probable lengths for no. 1–4of Figure 1, after loading applies. We first define strain measures from the evaluated peaksas strain measure for each peak = ( final - initial ) value of the peakfinal value of the peak . (12)Namely, for each peak value we subtract and divide by its initial value (namely the peakvalue for strain level zero). We then plot the applied stress (in absolute value) as a functionof the above defined strain measures and obtain Figures 3–6. So, Figures 3–6 render P xx ( T P a ) reduced shift vector | ℘ R | Figure 3: Applied stress vs the strain defined by the left peak which initially was at 0.140nm. The slope of the fitted line is 1.534.how carbon–carbon distances denoted by no. 1–4 in Figure 1 change with applied loadingwhen one can do calculations at the molecular level. Figure 3 relate the applied stresswith the increment of the left peak when initially the peak is at 0.140 nm; we call this leftpeak strain measure as the reduced shift vector P R . Figure 4 plot the increment in theright peak when initially the peak is at 0.140 nm with applied stress; we call the strainmeasure calculated from the right peak as reduced shift vector P L . Figures 5 and 6 are9 .000 0.004 0.008 0.012 0.016 0.0200.000.010.020.030.040.050.06 P xx ( T P a ) reduced shift vector | ℘ L | Figure 4: Aplied stress vs the strain defined by the right peak which initially was at 0.140nm. The slope of the fitted line is 3.117. P xx ( T P a ) reduced vector | R | Figure 5: Applied stress vs strain defined by the right peak which initially was at 0.240nm. The slope of the fitted line is 1.102. -0.0060 -0.0055 -0.0050 -0.0045 -0.0040 -0.0035 -0.00300.0050.0100.0150.0200.0250.0300.0350.0400.045 P xx ( T P a ) reduced vector | L | Figure 6: Applied stress vs strain defined by the left peak which initially was at 0.240 nm.The dependence exhibits a parabolic form, nevertheless, not taking into account a smallregion near zero strain, the behavior is near linear with a slope of -12.7.10nalogous for the peak which initially was at 0.240 nm; the corresponding strain measuresare denoted as reduced vectors | R | and | L | , respectively.One remark is in order here regarding the definition of the applied stress P xx . Thewhole process is deformation controlled. Nevertheless, one may convert the applied strainto stress using the manual convertion of the LAMMPS programs, which is based on theuse of the virial theorem. This renders a three dimensional definition of applied stress. Toprojected this quantity to graphene’s sheet one has to divide this three dimensional stessby graphene’s thickness, taken to be approximately 0.335 nm.We now define the material parameters c , c , c , c as the slope of diagrams 3-6. Theoutcome regarding the sorresponding slopes render values 3 . , . , . , − .
07 inTPa. Inspecting these diagrams one can see the different behaviour of Figure 6 from therest figures. The slope of diagram 6 is negative, while for the rest figures render a positiveslope. This is physically reasonable since we expect length no. 1 in Figure 2 to shorten.Also, it has a parabolic character, nevertheless, not taking into account a small region nearzero strain, the behaviour is near linear, as is seen in Figure 6.Now, what it remains to be done is to juxtapose lengths no. 1–4 to material parameters c , c , c , c . At a first look it seems reasonable to associate c , c with 1 and 2 and c , c with 3 and 4, namely to have the fourtuple ( c , c , c , c ) = ( − . , , , . , . c , c , c , c pertain. Further discussion on the subject is presented in the following section wherevalidation comparing with standard measurement is presented. There we see that the ex-pectation of corresponding lenghts no. 1, 2, 3, 4 of Figure 1 to c , c , c , c is not borneout, when we take as a minimum requirement to have values for ( E, ν ) comparable to theone’s from literature.To sum up, we define material parameters ( c , c , c , c ) as the slope in the diagramsof applied stress versus the four newly defined strain measures. These new strain mea-sures pertain to length changes denoted as numbers 1–4 in Figure 1. The values are then − . , . , . , .
534 in TPa. 11igure 7: On the left we see the disrete view within the unit cell of graphene. On theright we have the continuum analog, which is a continuum point with material parameters c , c , c , c . As a method of validation/calibration of our theory, we compare our calculated values withsome of the most well-accepted measurements from the vast literature on the topic (see theIntroduction section for relevant citations). This can also help as a guideline for makingthe correct association between material parameters c , c , c , c and the evaluated slopes.To obtain the relation between ( E, ν ) and ( c , c , c , c ) we start by solving the equationsruling the shift vector: these are eqs. (11). From eq. (8), setting the right hand side equalto zero, we obtain 2 c p = 4 c e c . (13)Eq. (9), with a zero on the right hand, solves to give c p = c ( e − e ) c . (14)Substituting these expressions to eqs. (3–5) we obtain in matrix form σ σ σ = c − c c c + c c c + c c c − c c
00 0 c ( c − c ) − c c e e e . (15)Inversion of the above relations render in matrix form e e e = αα − β − βα − β − βα − β αα − β
00 0 c σ σ σ , (16)where α = c − c c , β = c + c c and c = c ( c − c ) − c c . Setting all stress tensor componentsequal to zero except σ we obtain that e = αα − β σ . (17)12ince for the linear case Young’s modulus is defined as e = E σ , we obtain for ourcase E = α − β α = [ c − c c ] − [ c + c c ] c − c c . (18)Poisson ration is defined for our framework as ν = − e e = − − βα − β αα − β = βα = c + c c c − c c . (19)The last two equations are the connection between our mechanical approach with thestandard material parameters of the linear modeling of graphene, at small strains. It isobvious that material parameters c , c , c , c work synergetically to produce the Youngmodulus and the Poisson ratio. One cannot attribute changes of length to the directionof loading to only one of the material parameters c , c , c , c . The same holds true forchanges in length along the direction perpendicular to loading. Also, it is obvious that onecannot invert eqs. (18, 19) to solve in a unique way for c , c , c , c as functions of ( E, ν ).This is expected since the present framework has four material parameters while classicallinear elasticity has only two.This should not be confused with the inversion procedure for the passage from eq. (15)to eq. (16). There, we solve the equation ruling the shift vector, eq. (11), to obtain theshift vector as a function of the strain components (see eqs. (13, 14)). Due to the fact thatwe use a linear theory, the relation between the shift vector components and the straincomponents is linear. Having the shift vector components as linear functions of the straincomponents we substitute them to the constitutive law (eqs. (3-5)). This way we ontaineq. (15), which is a non-standard constitutive law (since it contains c , c , c , c ) relatingstresses to strains. Even though it is non-standard it is linear. Thus, it can be invertedgiving eq. (16).To find the appropriate values of c , c , c , c we choose from the pool of values − . , . , . , .
534 in TPa and substitute them to eqs. (18, 19). The choice ( c , c , c , c ) =(1 . , . , . , − .
07) is the optimum, in the sense of having values of (
E, ν ), cal-culated from eqs. (18, 19), as close to the one reported in literature as possible. It isclear that the outcome values of (
E, ν ) = (1 .
37 TPa , .
41) overestimate both measures,but nevertheless, remains within the range of reasonable values for these measures.On the other hand, for the calculation of the Poisson’s ratio and Young modulusfrom the molecular dynamics simulations solely (namely, without the need of introducing c , c , c , c ), atomic configurations are extracted from the simulation trajectory at regulartime step intervals that correspond to regular increases in applied strain. The calculatedvalues are ( E, ν )=(0.95 TPa, 0.20). Care is taken so that the structure has sufficient timeto equilibriate following the latest deformation (strain increase). The Poisson’s ratio iscalculated by measuring the changes of both the lateral and transverse dimensions of thecomputational cell. Using the same information the corresponding strain levels are recordedas well. For the Young modulus, in addition to the previous, the corresponding applied13ressures P xx at the graphene edges are also needed, of course suitably scaled to accountfor the difference between height of the computational cell and thickness of graphene. TheYoung modulus is then calculated as the slope of the pressure–strain diagram at smallstrains ( < ν from the reported literature (see the first four paragraphs of thisSection) range as E=0.67-1.91 TPa and ν =0.14-0.45 depending on the methodology used.The central tendency of these values for E is the value 1 ± E, ν )=(1.37 TPa, 0.41) overestimates these quantities but stillremain within the range of acceptable values. From the literature cited, the continuummethods (i.e. the finite element approaches of [2, 22, 24]) tend to have greater discrepancyfrom the value 1 ± E, ν )=(0.95TPa, 0.20). These values are obtained using the definition of the AIREBO manual. Onthe other hand, the values (
E, ν )=(1.37 TPa, 0.41) are based on a different definition of(
E, ν ): they are based on a genuinely continuous definition which is non-standard since ituses c , c , c , c . Certainly, the two definitions (the discrete and the continous one) measurethe same quantities in a different way. So, the discrepancy in their reported values is basedon the different definition but still remains in the range of an admissible difference.So, all in all, from molecular calculations we determine/define the following valuesfor the material parameters ( c , c , c , c ) = (1 . , . , . , − .
07) in TPa. Theseare the material parameters needed when graphene is modeled as a hexagonal 2-latticeat the continuum level. They appear to the non-standard constitutive law (eqs. (3-5))and characterize the stress-strain response in this case. Ultimately, they lead to values(
E, ν ) = (1 .
37 TPa , .
41) through eqs. (18, 19); these are slightly overestimated valueswhich nevertheless, remain in the range of reasonably accepted values.
The present work involves a molecular study with the purpose of measuring in-plane ma-terial moduli for graphene at the continuum level. The theoretical framework adopted isrestricted to material and geometrical linearities. Graphene is modeled as a hexagonal2-lattice, so for the linear regime there are four material parameters for in-plane motions.We evaluate these material parameters using molecular mechanics and the AIREBO po-tential. The material parameters are defined as the slopes of stress–strain diagrams ofsuitably defined strain measures from graphene’s unit cell at the discrete level. The finalvalues evaluated are ( c , c , c , c ) = (1 . , . , . , − .
07) in TPa and correspond toYoung modulus and Poisson ratio (
E, ν ) = (1 .
37 TPa , . Acknowledgments
We thank G.I. Sfyris (Athens, Greece) for reading the manuscript throughout and fornumerous most valuable addittions in the draft. This research has been co-financed bythe European Union (European Social Fund – ESF) and Greek national funds throughthe Operational Program ”Education and Lifelong Learning” of the National StrategicReference Framework (NSRF) Research Funding Program: ERC-10 ”Deformation, Yieldand Failure of Graphene and Graphene- based Nanocomposites”. The financial support ofthe European Research Council through the projects ERC AdG 2013 (Tailor Graphene)is greatfully acknowledged. The authors would like to acknowledge the financial sup-port of Graphene FET Flagship (Graphene-Based Revolutions in ICT and Beyond, GrantNo. 604391N.M.P. is supported by the European Research Council (ERC StG Ideas 2011BIHSNAM no. 279985 on Bio-inspired hierarchical supernanomaterials, ERC PoC 2013-1REPLICA2 no. 619448 on Large-area replication of biological anti-adhesive nanosurfaces,ERC PoC 2013-2 KNOTOUGH no. 632277 on Super-tough knotted fibres), by the Euro-pean Commission under the Graphene Flagship (WP10 Nanocomposites, no. 604391) andby the Provincia Autonoma di Trento (Graphene nanocomposites, no. S116/2012-242637and reg. delib. no. 2266).
References [1] M. Arroyo, T. Belytscko, An atomistic-based finite deformation membrane for singlelayer crystalline films. J. Mech. Phys. Sol. 50, 1941-1977 (2002).152] M. Arroyo, T. Belytschko, Finite crystal elasticity of carbon nanotubes based on theexponential Cauchy-Born rule. Phys. Rev. B, 69, 115415 (2004).[3] E. Cadelano, P.L. Palla, S. Giordano, L. Colombo, Nonlinear elasticity of monolayergraphene. Phys. Rev. Lett. 102, 235502 (2009).[4] J.L. Ericksen, On the symmetry of deformable crystals. Arch. Rat. Mech. Anal. 72,1-19 (1979).[5] J.L. Ericksen, On the Cauchy-Born Rule, Math. Mech. Sol. 13, 199-220 (2008).[6] G. Fadda, G. Zanzotto, The arithmetic symmetry of monoatomic 2-nets. Acta Cryst.(2000) A56, 36-48.[7] C. Galiotis, O. Frank, E.N. Koukaras, D. Sfyris, Graphene mechanics: current statusand perspectives. Annu. Rev. Chem. Biomol. Eng., to appear (2015).[8] A.K. Geim, K.S. Novoselov, The rise of graphene, Nature Materials 6, 183-191 (2007).[9] S. Gupta, K. Dharamvir, V.K. Jindal, Elastic moduli of single-walled carbon nanotubesand their ropes. Phys. Rev. B 72, 165428 (2005).[10] E. Hernandez, C. Goze, P. Bernier, A.; Rubio, Elastic properties of C and bxcynzcomposite nanotubes. Phys. Rev. Lett. 80, 45024505 (1998).[11] G. Kalosakas, N.N. Lathiotakis, C. Galiotis, K. Papagelis, In plane force fields andelastic properties of graphene. J. Appl. Phys. 113, 134307 (2013).[12] E. Konstantinova, S.O. Dantas, P.M.V.B. Barone, Electronic and elastic properties oftwo-dimensional carbon planes. Phys. Rev. B 74, 035417 (2006).[13] K.N. Kudin, E. Scusersia, B.I. Yakobson, C F, BN, and C nanoshell elasticity fromab initio computations. Phys. Rev. B, 64, 235406 (2001).[14] C. Lee, X.D. Wei, J.W. Kysar, J. Hone, Measurement of the elastic properties andinstrinsic sterngth of monolayer graphene. Science 321, 385-388 (2008).[15] F. Liu, P. Ming, J. Li, Ab initio calculation of ideal strength and phonon instabilityof graphene under tension. Phys. Rev. B, 76, 064120 (2007).[16] K.H. Michel, B. Verberck, Theory of the elastic constants of graphite and graphene.Phys. Stat. Sol. B, 245 2177-2180 (2008).[17] G.P. Parry, On diatomic crystals. Int. J. Sol. Struct. 14, 283-287 (1978).[18] M. Pitteri, Reconciliation of local and global symmetries of crystals. J. Elast. 14,175-190 (1984). 1619] M. Pitteri, On ν +1 lattices. J. Elast. 15, 3-25 (1985).[20] M. Pitteri, G. Zanzotto, Continuum models for phase transition and twinning incrystals, Chapman and Hall, Boca Raton (2003).[21] S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Com-put. Phys. 117, 1-19 (1995).[22] C. D. Reddy, S. Rajendran and K. M. Liew, Equivalent continuum modeling ofgraphene sheets. Int. J. Nanosci. 4, 631-636 (2005).[23] C.D. Reddy, S. Rajendran, K.M. Liew, Equilibrium configuration and continuum elas-tic properties of finite sized graphene. Nanotechnology 17, 864 (2006).[24] F. Scarpa, S. Adhikari, A. Srikantha Phani,Effective elastic mechanical properties ofsingle layer graphene sheets Nanotechnology 20 (2009) 065709.[25] D. Sfyris, C. Galiotis, Curvature dependent surface energy for free standing monolayergraphene. Math. Mech. Sol, in press, doi: 10.1177/108128651453667.[26] D. Sfyris, G.I. Sfyris, C. Galiotis, Curvature dependent surface energy for free standingmonolayer graphene. Some closed form solutions of the nonlinear theory. Int. J. NonlinearMech. 67, 186-197 (2014)[27] D. Sfyris, G.I. Sfyris, C. Galiotis, Curvature dependent surface energy for free standingmonolayer graphene. Geometrical and material linearization with closed form solutions.Int. J. Engng. Sci. 85, 224-233 (2014).[28] D. Sfyris, G.I. Sfyris, C. Galiotis, Constitutive modeling of some 2D crystals:graphene, hexagonal BN, MoS , WSe and NbSe2