An atomistic-based Föppl-von Kármán model for graphene
Cesare Davini Antonino Favata Roberto Paroni

	















 An atomistic-based 
Föppl–von Kármán model for graphene 
    
Cesare Davini^1 Antonino Favata^2 Roberto Paroni^3

    
==========================================================







	
	^1 Via Parenzo 17, 33100 Udine

		mailto:cesare.davini@uniud.itcesare.davini@uniud.it

		
		^2 Department of Structural and Geotechnical Engineering

		Sapienza University of Rome, Rome, Italy

		mailto:antonino.favata@uniroma1.itantonino.favata@uniroma1.it

		
		
		^3 Dipartimento di Ingegneria Civile e Industriale

	Università di Pisa, Pisa, Italy

		mailto:roberto.paroni@unipi.itroberto.paroni@unipi.it
	





myheadings
C. Davini, A. Favata, R. Paroni 
An atomistic-based Föppl–von Kármán model for graphene







§ ABSTRACT


We  deduce a non-linear continuum model of graphene  for the case of finite out-of-plane displacements and  small in-plane deformations. On assuming that the lattice interactions are governed by the Brenner's REBO potential of 2nd generation and that self-stress is present, we introduce discrete strain measures accounting for up-to-the-third neighbor interactions. The continuum limit turns out to depend on an average (macroscopic) displacement field and  a relative shift displacement  of the two Bravais lattices that give rise to the hexagonal periodicity. On minimizing the energy with respect to the shift variable, we formally determine a continuum model  of Föppl–von Kármán type, whose constitutive coefficients are given in terms of the atomistic interactions.






Keywords: Graphene, Continuum Modeling, Föppl–von Kármán equations.










§ INTRODUCTION


For its extraordinary mechanical, electrical, and thermal properties graphene is one of the most studied materials of the last two decades. Its discovery gave the occasion to renew a classical debate on the stability of 2D materials in nature <cit.>, and opened the way to the discussion of many theoretical issues. From the onset it was clear that the applications in various fields of technology might have been revolutionary. Though it's fair to say that graphene's potentialities are far from being fully explored and exploited, and remain the object of an intensive study, as much as it would be difficult to give an account of the huge literature on the subject. For a general picture of the subject and the new technological applications see the review by Ferrari et al. <cit.>.

The availability of macroscopic models is crucial to design applications and   experiments. The simplest models of structural mechanics such as membranes, plates, and shells have been often adopted in the past; in some cases, they have been assumed as a priori models and  the relevant constitutive constants have been estimated from ab initio or Molecular Dynamics simulations. Huang et al. <cit.> highlight how this has resorted to a bit of a stretch in some cases and led to paradoxical conclusions. More recent contributions start from atomistic analyses based on appropriate constitutive assumptions  on the interatomic potentials to obtain the continuum models of structural mechanics.

 Thus, Lu and Huang <cit.> estimated the elastic modulus and  bending stiffness of a graphene sheet from Molecular Mechanics calculations by  considering the one-dimensional stretching and the rolling on cylinders of various radii of a rectangular piece of graphene, and assuming that the interatomic forces are ruled by  Brenner's REBO potential of the second type. In this way they obtained values of the elastic constant  and bending stiffness that closely agree with those found by Kudin et al. <cit.> from ab initio calculations. These values are twice as much as  those found by Arroyo and Belitschko <cit.> in a paper where the dihedral contribution in Brenner's potential is taken into account; Arroyo and Belitschko also gave an atomistic-based membrane model for single layer crystalline films,  <cit.>. Finally, Davini  <cit.> deduced a 2D continuum model for the in-plane deformations of a graphene sheet within the framework of Γ-convergence. 

The out-of-plane deformations have been considered by various authors.  In particular,  by exploiting a formal analysis  we  deduced a continuum model of a graphene sheet <cit.>, and provided explicit expressions for the bending and Gaussian stiffness by starting from the study of the lattice kinematics and assuming the reactive empirical bond-order potential (REBO) of 2nd generation by Brenner et al. <cit.>. The approach takes into account the role of self-stress  and provides a quantitative estimate of the self-stress contribution to the overall bending and Gaussian stiffness. Indeed, the continuum model turns out to be the Γ-limit of the discrete graphene sheet, as proven in <cit.>. 

To understand the bending behavior of graphene is of the essence for several technological applications. The bending behavior controls the ripple formation and the performance of graphene nano-electro-mechanical devices <cit.>, and it is regarded as crucial in order to produce efficient hydrogen-storage devices <cit.>; moreover it can be instrumental to get inspiration for designing new metamaterials <cit.>.  Indeed, the intrinsic ripples are believed to be essential for the structural stability of the 2D graphene lattice and may have major impacts on the electronic and mechanical properties of graphene <cit.>.

In a recent review on Materials Today,  Deng and  Berry <cit.> give an overview of the hot problem of wrinkling, rippling, and crumpling, highlighting both formation mechanism and applications. The formation of these corrugations may have various explanations, see <cit.>. Basically, the out-of-plane deformations (wrinkles and ripples) can significantly reduce
the magnitude of in-plane stresses generated, for instance, by defects, <cit.>. Zhang et al. <cit.> adopted a generalized Föppl–von Kármán equation for a flexible solid membrane to describe ripples near defects such as disclinations (heptagons or pentagons) and dislocations (heptagon-pentagon dipoles) on graphene, and predicted the large scale graphene configurations under specific defect distributions. The paper closely follows a study of Seung and Nelson <cit.>. Comparison with atomistic simulations indicates that the proposed model is capable to predict the atomic scale wrinkles near disclination/dislocation cores. The analysis shows that considering the buckling into 3-dimensional deformations is energetically more favorable than restricting to the in-plane ones. Similar defect-guided ripples in graphene were also simulated and discussed in the work of Wang et al. <cit.>.

With an eye toward wrinkling and ripple formation, here we deduce a continuum model of graphene for the case of finite out-of-plane displacements and  small in-plane deformations. We consider an array of C-atoms sitting at the nodes of a hexagonal lattice, and assume that the lattice interactions are governed by the Brenner's REBO potential of 2nd generation and that self-stress is present. Thus, the starting point is the same as in <cit.>, but the changes of edge lengths, wedge angles and dihedral angles are calculated by keeping  the quadratic term in the out-of-plane displacements, according to the form of the in-plane Green-Lagrange strain used in Föppl–von Kármán plate theory.

The computation of the approximated measures of strain is done in Section 3. Unlike Zhang et al., <cit.>, that assume a triangular lattice for the continuum analysis as  done by Seung and Nelson in <cit.>, here we use the real geometry of a hexagonal lattice. With due modifications, we adopt a harmonic approximation of the interatomic potential, which yields a splitting of the energy into membrane and bending parts, see Section 4. It follows that the bending part keeps the form already discussed in <cit.>, while the membrane part turns out to be affected by the non-linearity of the assumed in-plane Green-Lagrange strain. The continuum limit of the membrane energy is computed in Section 5 according to the formal approach followed in <cit.>. The founding assumption is that, to within a 
remainder tending to zero with the lattice size, the nodal displacements can be described by an average (macroscopic) displacement  and a relative shift displacement  of the two Bravais lattices that give rise to the hexagonal periodicity. On minimizing the energy with respect to the shift variable, we formally determine a continuum model  of Föppl–von Kármán type, whose constitutive coefficients are given in terms of the atomistic interactions.




A full validation of the obtained continuum limit within the scheme of Γ-convergence is left for future work.





§ KINEMATICS AND ENERGETICS 
 


As reference configuration we use the 2–lattice generated by two simple Bravais lattices 

    [ L_1(ℓ) = {𝐱∈ℝ^2:
    𝐱 = n^1ℓ_1 + n^2ℓ_2     (n^1,
    n^2) ∈ℤ^2 },;                                       L_2(ℓ) = ℓ𝐩 + L_1(ℓ), ]

simply shifted with respect to one another by ℓ, see
Fig. <ref>.



 
In (<ref>), ℓ denotes
the lattice size (the reference interatomic distance), while
ℓ_α and ℓ𝐩  respectively are the lattice
	vectors and the shift vector, with


    _1 =√(3)_1,   _2 = √(3)/2_1+
    3/2_2     𝐩 =
    √(3)/2_1+ 1/2_2.


The sides of the hexagonal cells in Figure <ref> stand for the bonds
between pairs of next nearest neighbor atoms and are represented by the
vectors

    _α = _α -    (α = 1, 2)     _3 = -  .

For convenience we also set 
    _3=_2-_1.


In what follows we denote by

    ^ℓ=n^1ℓ_1+n^2ℓ_2+mℓ,   (n^1,n^2,m)∈ℤ^2×{0,1}

the lattice points and label them by the triplets (n^1,n^2,m): the points with m=0 belong to L_1(ℓ), while those in L_2(ℓ) correspond to m=1.



Graphene energetics depends on the description chosen to mimic atomic  interactions. Our model is based on the 2nd-generation Brenner potential <cit.>, which is one of the most used in molecular dynamics simulations of graphene. Accordingly, the  binding energy V of an atomic aggregate is given as a sum over nearest neighbors:


    V=∑_i∑_j<i V_ij ,    V_ij=V_R(l_ij)+b_ij(ϑ_hij,Θ_hijk)V_A(l_ij),


where the individual effects of the repulsion and attraction functions V_R(l_ij) and V_A(l_ij), which model pair-wise interactions of the atoms i and j depending on their distance l_ij, are modulated by the bond-order function b_ij; for a given bond chain  h,i,j,k the function b_ij depends in a complex manner on the angle between the edges hi and ij and on the dihedral angle between the planes spanned by (hi,ij) and (ij,jk). This potential reveals that,  in order to properly account  for the mechanical behavior of  graphene, it is necessary to consider three types of energetic contributions: 

	
  * binary interactions between next nearest atoms  (edge bonds),
	
  * three-body interactions between consecutive pairs of next nearest atoms  (wedge bonds),
	
  * and four-body interactions between three consecutive pairs of next nearest atoms (dihedral bonds).  There are two types of relevant dihedral bonds: the Z-dihedra, in which the edges connecting the four atoms form a Z-shape, and the C-dihedra, in which the edges form a C-shape (see Fig. <ref>).

Moreover, it is possible to show  <cit.> that  the angle  at ease between consecutive edges is greater than 2/3π: this means that in the flat reference configuration  the graphene sheet is not  stress-free, and we will proper account for this feature.






§ APPROXIMATED STRAIN MEASURES

In this section we calculate the strain measures associated to a change of configuration described by a displacement field ^ℓ : (L_1(ℓ)∪ L_2(ℓ))∩Ω→ℝ^3. Since we have in mind to deduce a model with the same non-linearities as the Föppl–von Kármán one, we write the displacements of the nodes in the form

    ^ℓ(^ℓ)= ξ(^ℓ)+ξ^1/2(^ℓ)_3,

where ξ is a positive scalar measuring smallness, :=ξ^-1^ℓ := ξ^-1(_1⊗_1+_2⊗_2)^ℓ and := ξ^-1/2:=w^ℓξ^-1/2^ℓ·_3 stand for the in-plane and out-of-plane normalized displacements, respectively. 




 §.§ Change of the edge lengths

With δ l_i(^ℓ) we denote the change in length of the edge parallel to _i and starting from the lattice point ^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω. We fix our attention to lattice points in L_2(ℓ). Thus,

    δ l_i(^ℓ)   =|(^ℓ + ℓ_i + ^ℓ(^ℓ+ℓ_i)) - (^ℓ +^ℓ(^ℓ))|-ℓ
       =ℓ|_i + ^ℓ(^ℓ+ℓ_i)- ^ℓ(^ℓ)/ℓ |-ℓ.

On introducing the notation

    if():=f(+ℓ_i)-f()/ℓ,

the axial strain measure can be recast as

    δℓ_i/ℓ_i=|_i+i^ℓ(^ℓ)|-1=|_i+ξi(^ℓ)+ξ^1/2i(^ℓ)_3|-1,

where we have made use of (<ref>). The   expansion up to the first order in ξ of the non-linear strain measure (<ref>) is

    δℓ_i/ℓ_i=ξ(i(^ℓ)·_i+1/2 |i(^ℓ)|^2  )+o(ξ),

which allows to define the  edge strain measure, once rescaled-back by ξ:

    ε_i(^ℓ):=i^ℓ(^ℓ)·_i+1/2 |iw^ℓ(^ℓ)|^2      i=1,2,3,

and δℓ_i=ε_iℓ_i.



 §.§ Change of the wedge angles


For each fixed node ^ℓ∈ L_2(ℓ))∩Ω we denote by _i(^ℓ) the angle of the wedge delimited by the edges _i+1 and _i+2; that is, the wedge angle opposite to the i-th edge (see Fig. <ref>). Here, i, i+1, and i+2 take values in {1,2,3} and the sums should be interpreted mod 3: for instance, if i=2 then i+1=3 and i+2=1.

To keep the notation compact, we set

    ^ℓ_i:=^ℓ(^ℓ+ℓ_i),    ^ℓ_0:=^ℓ(^ℓ).

Let 

    _i+1=ℓ(_i+1+ξi+1(^ℓ)+ξ^1/2i+1 (^ℓ)_3      )

and 

    _i+2=ℓ(_i+2+ξi+2(^ℓ)+ξ^1/2i+2 (^ℓ)_3      )


be the images of the edges parallel to _i+1 and _i+2 and starting at ^ℓ. Then, the angle _i(ξ)= is given by

    _i(ξ)=arccos(_i+1·_i+2/|_i+1||_i+2|).
 

Now, in the light of (<ref>) and (<ref>), we get

    1/ℓ^2_i+1·_i+2=_i+1·_i+2+ξ(_i+1·i+2(^ℓ)+_i+2·i+1(^ℓ)+ i+1(^ℓ)i+2(^ℓ)      ) +o(ξ),

and

    1/ℓ^2|_i+1|^2=1+ξ(|i+1(^ℓ)|^2 +2_i+1·i+1(^ℓ)  )+o(ξ),
       1/ℓ^2|_i+2|^2= 1+ξ(|i+2(^ℓ)|^2 +2_i+2·i+2(^ℓ)  ) +o(ξ),

so that the expansion up to the first order in ξ of _i(ξ) yields

    _i(ξ)=arccos(_i+1·_i+2+aξ+o(ξ)/√(1+bξ+o(ξ)) √(1+cξ+o(ξ)))= 2/3π - 1/2√(3) (4a+b+c) ξ+o(ξ)  ,

with

    a:=_i+1·i+2(^ℓ)+_i+2·i+1(^ℓ)+ i+1(^ℓ)i+2(^ℓ),
       b:=|i+1(^ℓ)|^2 +2_i+1·i+1(^ℓ),
       c:=|i+2(^ℓ)|^2 +2_i+2·i+2(^ℓ).

This allows to define the wedge strain measure, once the displacement is  rescaled-back by ξ:

    _i(^ℓ)=_i-2/3π=
    
    -1/√(3)( i+1^ℓ(^ℓ)·_i+1 +i+2^ℓ(^ℓ)·_i+2  +2( i+1^ℓ(^ℓ)·_i+2 +i+2^ℓ(^ℓ)·_i+1  )   +
    1	/2( |i+1w^ℓ(^ℓ)|^2+ |i+2w^ℓ(^ℓ)|^2)     +2i+1w^ℓ(^ℓ) i+2w^ℓ(^ℓ)) ,

for  ^ℓ∈ L_2(ℓ))∩Ω. 
In particular, algebraic manipulations allow to conclude that

    ∑_i=1^3_i(^ℓ)=-3√(3)(1/3∑_i=1^3iw^ℓ(^ℓ))^2.


If we consider a lattice point belonging to L_1(ℓ), it is not difficult to see that

    (1)ψ_i(^ℓ)=
    1/√(3)( i+1^ℓ(^ℓ)·_i+1 +i+2^ℓ(^ℓ)·_i+2  +2( i+1^ℓ(^ℓ)·_i+2 +i+2^ℓ(^ℓ)·_i+1  )   +
    
    -1	/2( |i+1w^ℓ(^ℓ)|^2+ |i+2w^ℓ(^ℓ)|^2)     -2i+1w^ℓ(^ℓ) i+2w^ℓ(^ℓ)) .







 §.§ Change of the dihedral angles

For each fixed node ^ℓ∈ L_2(ℓ)∩Ω and  for each edge parallel to _i
and starting at ^ℓ we need to define four types of dihedral angles __i^+(^ℓ), __i^-(^ℓ), __i_i+1(^ℓ) and __i_i+2(^ℓ):


    cos__i^+=(_i×_i+1)·(_i×_i^+)/|_i×_i+1||_i×_i^+|,
       cos__i^-=(_i+2×_i)·(_i^-×_i)/|_i+2×_i||_i^-×_i)|,
       cos__i_i+1=(_i×_i+1)·(_i^-×_i)/|_i×_i+1||_i^-×_i|,
       cos__i_i+2=(_i+2×_i)·(_i×_i^+)/|_i+2×_i||_i×_i^+|,

where	

    _i^+=   ^ℓ+ℓ_i-ℓ_i+2+^ℓ(^ℓ+ℓ_i-ℓ_i+2)-( ^ℓ+ℓ_i+^ℓ(^ℓ+ℓ_i) )
    
    =   -ℓ_i+2+^ℓ_i^+-^ℓ_i,         ^ℓ_i^+:=^ℓ(^ℓ+ℓ_i-ℓ_i+2),
    _i^-=   ^ℓ+ℓ_i-ℓ_i+1+^ℓ(^ℓ+ℓ_i-ℓ_i+1)-( ^ℓ+ℓ_i+^ℓ(^ℓ+ℓ_i) )
    
    =   -ℓ_i+1+^ℓ_i^--^ℓ_i,         ^ℓ_i^-:=^ℓ(^ℓ+ℓ_i-ℓ_i+1)

are the images of vectors ℓ_i^+ and ℓ_i^- (see Fig. <ref>, for i=1), parallel to _i+2 and _i+1 and starting at the image of the point ^ℓ+ℓ_i.


Also here, i, i+1, and i+2 take values in {1,2,3} and the sums should be interpreted mod 3: for instance, if i=3 then i+1=1 and i+2=2.

The  C-dihedral angle __i^+(^ℓ) is the angle corresponding to the C-dihedron with middle edge ℓ_i and oriented as _i^⊥,  while __i^-(^ℓ) is the angle corresponding to the C-dihedron oriented opposite to _i^⊥ (see Fig. <ref> for i=1). 
The Z-dihedral angle __i_i+1(^ℓ) 
corresponds to the Z-dihedron with middle edge ℓ_i and the other two edges parallel to _i+1
(see Fig. <ref> for i=1).




To fix the ideas, we focus on the  dihedral angle __1_2, sketched in Fig. <ref>; the other strains can be obtained in analogous manner.  

The first order approximation of the dihedral angle is all we need to evaluate the corresponding energy contribution. 

Let us introduce the vector 

    _5=ℓ_5+(_5-_1)ξ+_5ξ^1/2,

image of _5 under the deformation. We have that 

    __1_2=arccos(   _1×_2/|_1×_2|·_5×_1/|_5×_1|).

Cumbersome computations yield:

    _1×_2/|_1×_2|·_5×_1/|_5×_1|  =1-2/3ℓ^2(_5-_1+_2-_0)^2ξ+o(ξ).

The expansion up to the first order in ξ of __1_2 yields

    __1_2(ξ)=arccos(1-a^2ξ+o(ξ))=√(2 )a ξ^1/2+O(ξ^3/2),    a^2:=2/3ℓ^2(_5-_1+_2-_0)^2.

This allows to define the Z-dihedron strain measure, once the displacement be  rescaled-back by ξ:

    __1_2=__1_2(ξ)-arccos(   _1×_2/|_1×_2|·_5×_1/|_5×_1|)=2√(3)/3ℓ(w^ℓ_5-w^ℓ_1+w^ℓ_2-w^ℓ_0  ).

For a generic Z-dihedral angle centered in _i, we get

    __i_i+1(^ℓ)=2√(3)/3ℓ[w^ℓ(^ℓ+ℓ_i-ℓ_i+1)-w^ℓ(^ℓ+ℓ_i)+w^ℓ(^ℓ+ℓ_i+1)-w^ℓ(^ℓ)],
       __i_i+2(^ℓ)=2√(3)/3ℓ[w^ℓ(^ℓ+ℓ_i-ℓ_i+2)-w^ℓ(^ℓ+ℓ_i)+w^ℓ(^ℓ+ℓ_i+2)-w^ℓ(^ℓ)].

Analogous computations allows to determine the C-dihedron strain measure

    __i^+(^ℓ)=2√(3)/3ℓ[2w^ℓ(^ℓ)-w^ℓ(^ℓ+ℓ_i+1)+w^ℓ(^ℓ+ℓ_i-ℓ_i+2)-2w^ℓ(^ℓ+ℓ_i)],
       __i^-(^ℓ)=-2√(3)/3ℓ[2w^ℓ(^ℓ)-w^ℓ(^ℓ+ℓ_i+2)+w^ℓ(^ℓ+ℓ_i-ℓ_i+1)-2w^ℓ(^ℓ+ℓ_i)].
		






§ MEMBRANE AND BENDING ENERGY


The above calculations show that ε and  ψ depend upon both the in-plane and the out-of-plane  components of ^ℓ, while , and  depend upon the out-of-plane component of ^ℓ; moreover, (<ref>) shows that the sum of all ψ depends on  the out-of-plane component of ^ℓ. We introduce the following splitting of the energy into membrane and bending parts:

    𝒰_ℓ=𝒰_ℓ^(m)+𝒰_ℓ^(b),     𝒰_ℓ^(m):=  𝒰_ℓ^(e)+ 𝒰_ℓ^(w)  𝒰_ℓ^(b):=_ℓ^(s) +_ℓ^(d)
   
defined by

    𝒰_ℓ^(e):=1/2 ∑_ℰ k^l (δℓ)^2=1/2 ∑_ℰ k^l ℓ^2ε^2 
       𝒰_ℓ^(w):=1/2 ∑_𝒲 k^ϑℓ^2  ψ^2 
       _ℓ^(s) :=τ_0∑_𝒲ψ, 
       _ℓ^(d) :=1/2 ∑_𝒵 k^ ^2+1/2 ∑_𝒞 k^ ^2,

where 𝒰_ℓ^(e),  𝒰_ℓ^(w),  _ℓ^(s) and _ℓ^(d) are the 
edge, wedge, self-, and dihedral energy, respectively.

The self-stress term is the outcome of the fact  that  the angle  at ease between consecutive edges is greater than 2/3π: this means that in the flat reference configuration  the graphene sheet is not  stress-free, and τ_0 represents the pre-stress couple (see <cit.>).  The constants k^l, k^ϑ, k^, and k^ can be deduced by means of the 2nd-generation Brenner potential. In the wedge energy we may interpret the presence of ℓ^2 as a scaling of the constant 
k^ϑ,  introduced to keep the energy finite as the lattice size ℓ goes to zero. We notice that the multiplication factor ℓ^2 appears in all contributions to the membrane energy.


 

With the notation introduced in Section <ref> we now write the  energy more explicitly.
The edge energy can be written as 

    𝒰_ℓ^(e)=1/2  k^lℓ^2 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3   ε_i^2(^ℓ),

while the wedge energy reads:

    𝒰_ℓ^(w)=1/2  k^ϑℓ^2∑_^ℓ∈ L_1(ℓ)∩Ω∑_i=1^3  (1)ψ_i^2(^ℓ)+ 1/2  k^ϑℓ^2 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3  (2)ψ_i^2(^ℓ).

Analogously, the self-energy becomes

    𝒰^(s)_ℓ
    = ∑_^ℓ∈ L_1(ℓ)∩Ωτ_0 ∑_i = 
    	1^3(1)ψ_i (^ℓ)+∑_^ℓ∈ L_2(ℓ)∩Ωτ_0 ∑_i = 
    	1^3(2)ψ_i (^ℓ).

We further split the dihedral energy  _ℓ^(d) in

    _ℓ^(d):=𝒰^𝒵_ℓ+𝒰^𝒞_ℓ,

where

    𝒰^𝒵_ℓ=1/2  k^𝒵∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (__i_i+2(^ℓ))^2+ (__i_i+1(^ℓ))^2

is the contribution of the Z-dihedra, and

    𝒰^𝒞_ℓ=1/2 k^𝒞 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (__i^+(^ℓ))^2+ (__i^-(^ℓ))^2

is the contribution of the C-dihedra. 

We notice that the bending energy here obtained is the same as <cit.>. In the following, we then focus on the membrane energy.












§ THE CONTINUUM LIMIT OF THE MEMBRANE ENERGY


In the previous sections the discrete  energy 𝒰_ℓ was defined over the lattice (L_1(ℓ)∪ L_2(ℓ))∩Ω. By letting the lattice size ℓ go to zero the discrete set 
(L_1(ℓ)∪ L_2(ℓ))∩Ω invades the domain Ω, and the displacement functions
w^ℓ and ^ℓ  will approach two functions w and  defined over Ω. 
To  derive a continuous energy, defined over the domain Ω, from
the discrete  energy 𝒰_ℓ
we need to specify the relations between w^ℓ and w and between ^ℓ and .

We assume  w:Ω→ℝ to be twice continuously differentiable and to, almost, coincide with
 over the lattice (L_1(ℓ)∪ L_2(ℓ))∩Ω. More precisely, we assume that

    (^ℓ)=w(^ℓ)+o(ℓ),     ∀^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω.

If we think of w  and  as the macroscopic and microscopic displacements, respectively, then (<ref>) can be
thought as a Cauchy–Born rule. The assumption  (<ref>) is motivated and essentially justified in <cit.>.

For the in-plane displacement we assume :Ω→ℝ^2 to be continuously differentiable 
and 

    (^ℓ)   =(^ℓ)+o(ℓ),      ∀   ^ℓ∈ L_1(ℓ)∩Ω,
    (^ℓ)   =(^ℓ)-ℓ(^ℓ)+o(ℓ),        ∀   ^ℓ∈ L_2(ℓ)∩Ω.


Thus, over the lattice L_1(ℓ) we make the Cauchy–Born assumption for the in-plane displacement, while
over the lattice L_2(ℓ) this assumption is relaxed by introducing a “shift displacement” 
:Ω→ℝ^2 that we assume to be continuously differentiable. Clearly, if the shift displacement is set equal to zero we have the Cauchy–Born rule over both lattices, but energetically it might be convenient to have a shift displacement  different from zero. The minus sign in front of the term containing  is introduced simply for later convenience. The assumption (<ref>) is motivated and essentially justified in <cit.>. 

We note that for ^ℓ∈ L_2(ℓ) we have that ^ℓ+ℓ_i∈ L_1(ℓ), for i=1,2,3, and 

    (^ℓ+ℓ_i)=(^ℓ+ℓ_i)+o(ℓ)=(^ℓ)+ℓ∇(^ℓ)_i+o(ℓ)=(^ℓ)+ℓ(^ℓ)+ℓ∇(^ℓ)_i+o(ℓ),

while for ^ℓ∈ L_1(ℓ) we have that ^ℓ-ℓ_i∈ L_1(ℓ), for i=1,2,3, and 

    (^ℓ-ℓ_i)   =(^ℓ-ℓ_i)-ℓ(^ℓ-ℓ_i)+o(ℓ)=(^ℓ)-ℓ∇(^ℓ)_i-ℓ(^ℓ)+o(ℓ)
       =(^ℓ)-ℓ∇(^ℓ)_i-ℓ(^ℓ)+o(ℓ).

Thus,  

    i(^ℓ)   =(^ℓ+ℓ_i)-(^ℓ)/ℓ=∇(^ℓ)_i+(^ℓ)+o(1),    ^ℓ∈ L_2(ℓ),
    i(^ℓ)   =(^ℓ-ℓ_i)-(^ℓ)/ℓ=-∇(^ℓ)_i-(^ℓ)+o(1),    ^ℓ∈ L_1(ℓ),

and similarly for ^ℓ∈ L_2(ℓ) we have that

    i(^ℓ)=(^ℓ+ℓ_i)-(^ℓ)/ℓ=∇ w(^ℓ)·_i+o(1).


For i=1,2,3 and for ^ℓ∈ L_2(ℓ) the edge strain measure defined in (<ref>) writes as

    ε_i(^ℓ)   =_i·i(^ℓ)+1/2 |iw(^ℓ)|^2 
       =∇(^ℓ)_i·_i+(^ℓ)·_i+1/2 |∇ w(^ℓ)·_i|^2+o(1)
       =(^ℓ)_i·_i+(^ℓ)·_i+1/2 |∇ w(^ℓ)·_i|^2+o(1)
       =[(^ℓ)+1/2∇ w(^ℓ)⊗∇ w(^ℓ)]·_i⊗_i+(^ℓ)·_i+o(1)
       =[,w](^ℓ)·_i⊗_i+(^ℓ)·_i+o(1)

where

    :=∇ +(∇)^T/2,     [,w]:=+1/2∇ w⊗∇ w,

are the linearized and the von Kármán strain tensors, respectively.
Similarly, the wedge strain measure defined in (<ref>) rewrites as

    _i(^ℓ)   =-1/√(3)((∇(^ℓ)_i+1+(^ℓ))·_i+1 +(∇(^ℓ)_i+2+(^ℓ))·_i+2
       +2(∇(^ℓ)_i+1+(^ℓ))·_i+2 +2(∇(^ℓ)_i+2+(^ℓ))·_i+1
       +1	/2( |∇ w (^ℓ)·_i+1|^2+ |∇ w (^ℓ)·_i+2|^2)     +2∇ w (^ℓ)·_i+1 ∇ w (^ℓ)·_i+2)+o(1)
       =-1/√(3)((^ℓ)_i+1·_i+1 +(^ℓ)_i+2·_i+2+ 4 (^ℓ)_i+1·_i+2
       +3(^ℓ))·(_i+1 +_i+2)
       +1/2∇ w (^ℓ)⊗∇ w (^ℓ) ·( 
     _i+1⊗_i+1+ _i+2⊗_i+2+4 _i+1⊗_i+2)
    )+o(1)
       =-1/√(3)([,w](^ℓ)·_i -3(^ℓ)·_i)+o(1),

where we set

    _i:=_i+1⊗_i+1+ _i+2⊗_i+2+2 _i+1⊗_i+2+2 _i+2⊗_i+1,

and where we used the fact that _i+1 +_i+2=-_i.
A similar computation shows that, see (<ref>),

    (1)ψ_i(^ℓ)=-1/√(3)([,w](^ℓ)·_i -3(^ℓ)·_i)+o(1).


We now compute the energies. The edge energy (<ref>) becomes

    𝒰_ℓ^(e)   =1/2  k^lℓ^2 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3   ([,w](^ℓ)·_i⊗_i+(^ℓ)·_i+o(1))^2
       =o(1)+1/2  k^lℓ^2 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3   ([,w](^ℓ)·_i⊗_i+(^ℓ)·_i)^2,

 where the second equality is obtained by noticing that the number of points ^ℓ in L_2(ℓ)∩Ω
is of order 1/ℓ^2.
 Let |E^ℓ(^ℓ)| = ℓ^2 3√(3)/2 be the area of the hexagon E^ℓ(^ℓ) of side ℓ centred at ^ℓ, see Figure <ref>, and  let χ_E^ℓ(^ℓ)() be the characteristic function of E^ℓ(^ℓ), i.e., the function equal to 1 if ∈ E^ℓ(^ℓ) and 0 otherwise. The energy 𝒰_ℓ^(e) may be rewritten as
 
    𝒰_ℓ^(e)   =o(1)+k^l/3√(3) ∑_i=1^3 ∫_Ω∑_^ℓ∈ L_2(ℓ)∩Ω([,w](^ℓ)·_i⊗_i+(^ℓ)·_i)^2χ_E^ℓ(^ℓ)() d

and since the function ∑_^ℓ∈ L_2(ℓ)∩Ω([,w](^ℓ)·_i⊗_i+(^ℓ)·_i)^2χ_E^ℓ(^ℓ)() converges, as ℓ goes to zero, to ([,w]()·_i⊗_i+()·_i)^2
   we have that

    lim_ℓ→ 0𝒰_ℓ^(e)=k^l/3√(3) ∑_i=1^3 ∫_Ω([,w]()·_i⊗_i+()·_i)^2 d=: 𝒰_0^(e)(,w,).

 
By taking into account that the expression for  (1)ψ_i and
(2)ψ_i are identical, the wedge energy (<ref>) rewrites:


    𝒰_ℓ^(w)=o(1)+1/2  k^ϑℓ^2 ∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω∑_i=1^3  1/3([,w](^ℓ)·_i -3(^ℓ)·_i)^2





 

By introducing the triangles  T^ℓ(^ℓ) centered at ^ℓ of area 3√(3)ℓ^2/4 as depicted in
Figure <ref> and proceeding as above, we deduce that

    lim_ℓ→ 0𝒰_ℓ^(w)=2k^ϑ/9√(3) ∑_i=1^3 ∫_Ω([,w]()·_i-3()·_i)^2 d=:𝒰_0^(w)(,w,).

The limit of the membrane energy is

    𝒰_0^(m)(,w,)   :=lim_ℓ→ 0𝒰_ℓ^(m)=𝒰_0^(e)(,w,)+𝒰_0^(w)(,w,)
       =k^l/3√(3) ∑_i=1^3 ∫_Ω([,w]()·_i⊗_i+()·_i)^2 d
       +2k^ϑ/9√(3) ∑_i=1^3 ∫_Ω([,w]()·_i-3()·_i)^2 d
    
       = ∫_Ω∑_i=1^3 (k^l/3√(3)([,w]()·_i⊗_i)^2+
    2k^ϑ/9√(3)([,w]()·_i)^2)
    
       +·∑_i=1^3 (
    2k^l/3√(3)[,w]()·_i⊗_i-4k^ϑ/3√(3)[,w]()·_i
    )_i
       +  (k^l/3√(3)+2k^ϑ/√(3)) ∑_i=1^3(()·_i)^2 d

By using the relations

    _1=√(3)/2_1-1/2_2,     _2=_2,     _3=-√(3)/2_1-1/2_2

we find that

    _1·_1   =3/4 ()_11+1/4 ()_22-√(3)/2 ()_12,
     _2·_2   =()_22,
     _3·_3   =3/4 ()_11+1/4 ()_22+√(3)/2 ()_12,

where ()_αβ are the components of =[,w] with respect to the basis {_i}.
A simple computation then shows that

    ∑_i=1^3 ( _i·_i)^2   =1/8(9()_11^2+9()_22^2+6()_11()_22+12()_12^2),
       =3/8(3()^2-4).

Similarly we find 

    ∑_i=1^3 (·_i)^2=27/8(()^2-4),
    ∑_i=1^3 (·_i⊗_i)_i,=-3/2 ()_12_1-3/4(()_11-()_22)_2, 
    ∑_i=1^3 (·_i)_i=-9/2 ()_12_1-9/4(()_11-()_22)_2,
    ∑_i=1^3 (·_i)^2=3/2 ||^2.


With these identities it follows that

    𝒰_0^(m)(,w,)   = ∫_Ω3(k^l+2k^ϑ)/8√(3)([,w])^2
    -k^l+6k^ϑ/2√(3)[,w]
    
       -k^l-6k^ϑ/2√(3)·(
     2([,w])_12_1+
    ([,w])_11-([,w])_22)_2)
    
       +  k^l+6k^ϑ/2√(3)||^2 d.

The shift displacement  enters into the membranal energy without any derivatives and therefore
it can be minimized once and for all. We readily see that the shift displacement that minimizes the energy is

    =k^l-6k^ϑ/2(k^l+6k^ϑ)(
     2([,w])_12_1+
    ([,w])_11-([,w])_22)_2)

and that

    𝒰_0^(m)(,w)   :=min_𝒰_0^(m)(,w,)
    
       = ∫_Ω3(k^l+2k^ϑ)/8√(3)([,w])^2
    -k^l+6k^ϑ/2√(3)[,w]
    
       -(k^l-6k^ϑ)^2/8√(3)(k^l+6k^ϑ)|
     2([,w])_12_1+
    ([,w])_11-([,w])_22)_2|^2 d.

By expanding the squares and reorganizing the terms we find:

    𝒰_0^(m)(,w)   =∫_Ω3(k^l+6k^ϑ)/8√(3)([,w])^2
    -k^l+6k^ϑ/2√(3)[,w]
    
       -(k^l-6k^ϑ)^2/8√(3)(k^l+6k^ϑ)(([,w])^2-4[,w]) d,
    
       = ∫_Ωk^l(k^l+18k^ϑ)/4√(3)(k^l+6k^ϑ)([,w])^2
    -4√(3)k^lk^ϑ/k^l+6k^ϑ[,w]
     d,

which is clearly an isotropic energy.
By means of the relation =(()^2-||^2)/2, which holds for every two by two matrix , we may write

    𝒰_0^(m)(,w)   =
     ∫_Ω1/2λ([,w])^2
    +μ|[,w]|^2
     d,

with


    λ:=k^l(k^l-6k^ϑ)/2√(3)(k^l+6k^ϑ),     μ:=2√(3)k^lk^ϑ/k^l+6k^ϑ.





§ CONCLUSIONS

We have proposed a non-linear continuum model for the mechanical behavior of graphene inferred from Molecular Dynamics potentials. Starting from a harmonic approximation of the energy as depicted by the 2nd-generation Brenner potential, we have found a discrete energy, depending on the displacement of each atom, which is sensitive to change of (i) the distance between two atoms (edge energy), (ii) the angle spanned by three subsequent atoms (wedge energy), (iii) two types of dihedral angles generated by the plane spanned by four subsequent atoms (C- and Z-dihedral energy). Thus, up-to-third neighbors interactions have been considered. Moreover, we have taken into account the presence of the self-stress, as predicted by the 2nd-generation Brenner potential (self-energy).

We have introduced a different scaling for in-plane and out-of-plane components, and this feature has produced a  coupling in the  strain measures at the discrete level. In particular, while the dihedral and the self-stress energies depend just on  the out-of-plane components, the edge and the wedges energies depend on both. These two latter contributions determine the membrane energy, while the former two are part of the bending energy. 

With the scaling here adopted, the discrete bending energy  turns out to be the same as that already considered in <cit.>; for this reason, we focused on the membrane energy.

The deduced discrete energy  is defined over the two Bravais lattices  generating the graphene sheet. By letting the size of the lattices to zero, they invade a continuum domain Ω and the  discrete displacement functions approach two continuous functions w and , representing the  in-plane and the out-of-plane continuum displacements, defined over Ω. To obtain the continuum energy, it has been necessary to specify the relation between the discrete displacements and the corresponding continuum functions. To this end, motivated by <cit.>, we have made the Cauchy–Born assumption for the in-plane displacement defined over one of the Bravais lattices, while we have relaxed this assumption for the second lattice, by introducing a shift displacement . 

With these assumptions, we have found a continuum membrane energy depending on w,  and . Since  the shift displacement  enters into the membranal energy without any derivatives, we have minimized it once and for all; this has lead to a membrane energy depending just on w and  . 

Thus, on considering the bending energy already deduced in <cit.> and the membrane energy here found, we can state that the total continuum energy of graphene reads:

    𝒰_0(,w)=𝒰_0^(m)(,w)+𝒰_0^(b)(w),

where the membrane energy is given by

    𝒰_0^(m)(,w)   =
    ∫_Ω1/2λ([,w])^2
    +μ|[,w]|^2
     d,

and the bending energy is

    𝒰^(b)_0(w)
    =1/2∫_Ω ( Δ w)^2+_G∇^2 w d.

We find that 

    λ=k^l(k^l-6k^ϑ)/2√(3)(k^l+6k^ϑ),   μ=2√(3)k^lk^ϑ/k^l+6k^ϑ,   = 5√(3)/3 k^𝒵+2√(3)/3 k^𝒞-τ_0/2,   _G=- 8/55√(3)/3 k^𝒵-42√(3)/3 k^𝒞,

where  k^l, k^ϑ, k^𝒞, k^𝒵 and τ_0 are the  constants entering the edge, wedge, C-dihedral, Z-dihedral, and self-, atomistic energies.

Within the limits of this formal deduction, we have  found that graphene can be  modeled as a classical  Föppl–von Kármán plate, where the constitutive constants depend on  the atomistic interactions, as described by the 2nd-generation Brenner potential.







§ ACKNOWLEDGMENTS

A.F. acknowledges support from  Sapienza University of Rome through the   projects RP116154C92AF8A4  
“Multiscale Mechanics of 2D Materials: Modeling and Applications' and  RM11715C7F61C3E8 “Shape morphing. From advanced differential geometry to applications in engineering and architecture”.

R.P.  acknowledges support from the Università di Pisa through the project PRA_2018_61 “Modellazione multi-scala in ingegneria strutturale”.

	
	






10



arroyo2002
 M. Arroyo and T. Belytschko. An atomistic-based finite deformation membrane for single layer crystalline films, J. Mech. Phys. Solids, 50 (9) (2002), 1941 – 1977.

arroyo2004
 M. Arroyo and T. Belytschko. Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy–Born rule, Phys. Rev. B, 69 (11) (2004): 115415.

Brenner_2002
 D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, and S.B. Sinnott. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J. Phys. Cond. Mat., 14 (4) (2002), 783.

davini2014
 C. Davini, Homogenization of a graphene sheet.Cont. Mech. Thermod., 26 (1) (2014), 95–113. 

davini2017
 C. Davini, A. Favata, and R. Paroni. The Gaussian stiffness of graphene deduced from a continuum model based on Molecular Dynamics potentials, J. Mech. Phys. Solids, 104 (2017), 96 –114.



davini2017_1
 C. Davini, A. Favata, and R. Paroni.
A new material property of graphene: The bending Poisson coefficient.
EPL, 118 (2017), 26001.

davini2017_2
 C. Davini, A. Favata, A. Micheletti, and R. Paroni.
A 2D microstructure with auxetic out-of- plane behavior and non-auxetic in-plane behavior.
Smart Mater. Struct., 26 (2017) 125007.

davini2018
 C. Davini, A. Favata, and R. Paroni. A REBO-Potential-Based Model for Graphene Bending by Γ-Convergence, Arch. Rational Mech. Anal., 229 (3)(2018), 1153 – 1195.


Deng_2016
S. Deng and V. Berry.
Wrinkled, rippled and crumpled graphene: an overview of formation
mechanism, electronic properties, and applications.
Materials Today, 19(4):197 – 212, 2016.


fasolino2007
  Fasolino A., Los J.H., Katsnelson M.I., Intrinsic ripples in graphene, Nature Materials. Letters, 6 (2007), 858–861.

	Favata_2016
	A. Favata, A. Micheletti, P. Podio-Guidugli, and N.M. Pugno.
	Geometry and self-stress of single-wall carbon nanotubes and graphene
	via a discrete model based on a 2nd-generation REBO potential.
	J. Elasticity, 125:1–37, 2016.

Favata_2016b
A. Favata, A. Micheletti, P. Podio-Guidugli, and N.M. Pugno.
How graphene flexes and stretches under concomitant bending couples
and tractions.
Meccanica, (2017) 52:1601–1624

Ferrari2014
A.C. Ferrari, F. Bonaccorso, V. Fal'ko, K.S. Novoselov, S. Roche,
P. Bøggild, S. Borini, F.H.L. Koppens, V. Palermo, N.M. Pugno, J.A.
Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis,
J. Kivioja, C. Marinelli, T. Ryhänen, A. Morpurgo, J.N. Coleman,
V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G.F.
Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A.N.
Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen,
A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini,
A. Tredicucci, G.M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath,
B.J. van Wees, H. van der Zant, L. Occhipinti, A. Di Matteo, I.A. Kinloch,
T. Seyller, E. Quesnel, K. Feng, X.and Teo, N. Rupesinghe, P. Hakonen,
S. R.T. Neil, Q. Tannock, T. Löfwander, and J. Kinaret.
Science and technology roadmap for graphene, related two-dimensional
	crystals, and hybrid systems.
Nanoscale, 7(11):4587–5062, 2015.

Genoese_2019
A. Genoese, A. Genoese, N.L. Rizzi, G. Salerno.
Buckling analysis of single-layer graphene sheets using molecular mechanics.
Front. Mater., in press.

Goler_2013
S. Goler, C. Coletti, V. Tozzini, V. Piazza, T. Mashoff, F. Beltram,
V. Pellegrini, and S. Heun.
Influence of graphene curvature on hydrogen adsorption: Toward
hydrogen storage devices.
Phys. Chem. C, 117(22):11506–11513, 2013.

Hajgato_2012
B. Hajgató, S. Güryel, Y. Dauphin, J.-M. Blairon, H.E. Miltner,
G. Van Lier, F. De Proft, and P. Geerlings.
Theoretical investigation of the intrinsic mechanical properties of
single- and double-layer graphene.
J. Phys. Chem. C, 116(42):22608–22618, 2012.


Hartmann_2013
M.A. Hartmann, M. Todt, F.G. Rammerstorfer, F.D. Fischer, and O. Paris.
Elastic properties of graphene obtained by computational mechanical
tests.
Europhys. Lett., 103(6):68004, 2013.

huang2006
 Y. Huang,  J. Wu, and K. C. Hwang, Thickness of graphene and single-wall carbon nanotubes,Phys. Rev. B, 74 (1): 245413.

Kim_2012
S.M. Kim, E.B. Song, S. Lee, J. Zhu, D.H. Seo, M. Mecklenburg, S. Seo, and K.L.
Wang.
Transparent and flexible graphene charge-trap memory.
ACS Nano, 6(9):7879–7884, 2012.

kudin2001
 K.N. Kudin, G.E. Scuseria, and B.I. Yakobson. C_2F,  BN, and C nanoshell elasticity from ab initio computations, Phys. Rev. B, 64 (23) (2001): 235406.

landau1980
 Landau, L. D., Lifshits, E. M., Pitaevskii, L. P., Sykes, J. B. and Kearsley, M. J.,Statistical Physics, Part I, Pergamon, Oxford, (1980). 







	Lindahl_2012
	N. Lindahl, D. Midtvedt, J. Svensson, O.A. Nerushev, N. Lindvall, A. Isacsson,
	and E.E. B. Campbell.
	Determination of the bending rigidity of graphene via electrostatic
	actuation of buckled membranes.
	Nano Lett., 12(7):3526–3531, 2012.

luHuang2009
 Q. Lu and R. Huang. Nonlinear mechanics of single-atomic-layer graphene sheets, Int. J. Appl. Mech., 1(03) (2009), 443–467.

Lu_2009
 Q. Lu, M. Arroyo, and R. Huang. Elastic bending modulus of monolayer graphene. J. Phys. D, 42 (10) (2009),102002.

mermin1968
 Mermin, N. D., Crystalline order in two dimensions,Phys. Rev., 176 (1), 250–254. 

meyer2007
  Meyer J.C., Geim A.K., Katsenelson M.I., Novoselov K.S., Booth T.J., Roth S., The structure of suspended graphene sheets, Nature, 446 (2007), 60–63.

nelsonPeliti1987
  Nelson D.R. and Peliti L., Fluctuations in membranes with crystalline and hexatic order, J. Physique 48 (1987), 1085-1092.


Pacheco_2014
A.A. Pacheco Sanjuan, Z. Wang, H.P. Imani, M. Vanević, and
S. Barraza-Lopez.
Graphene's morphology and electronic properties from discrete
differential geometry.
Phys. Rev. B, 89:121403, 2014.

seungNelson1988
 Seung H. S.  and Nelson D. R., Defects in flexible membranes with crystalline order. Phys. Rev. A, 38 (1988), 1005–1018.


Shi_2012
X. Shi, B. Peng, N.M. Pugno, and H. Gao.
Stretch-induced softening of bending rigidity in graphene.
Appl. Phys. Let., 100(19), 2012.


Tapaszto_2012
L. Tapaszto, T. Dumitrica, S.J. Kim, P. Nemes-Incze, C. Hwang, and L.P. Biro.
Breakdown of continuum mechanics for nanometre-wavelength rippling
	of graphene.
Nat. Phys., 8(10):739–742, 2012.

	Tozzini_2011
	V. Tozzini and V. Pellegrini.
	Reversible hydrogen storage by controlled buckling of graphene
	layers.
	Phys. Chem. C, 115(51):25523–25528, 2011.
	
	Tozzini_2013
	V. Tozzini and V. Pellegrini.
	Prospects for hydrogen storage in graphene.
	Phys. Chem., 15:80–89, 2013.
	
	wang2013
	  Wang C.G. , Lan L. , Liu Y.P., Tan H.F., Defect-guided wrinkling in graphene, Comput. Mater. Sci., 77 (2013), 250–253.
	
	
	Wei_2013
	Y. Wei, B. Wang, J. Wu, R. Yang, and M.L. Dunn.
	Bending rigidity and Gaussian bending stiffness of
	single-layered graphene.
	Nano Lett., 13(1):26–30, 2013.

Zelisko_2017
M. Zelisko, F. Ahmadpoor, H. Gao, and P. Sharma
Determining the Gaussian Modulus and Edge Properties of 2D Materials: From Graphene to Lipid Bilayers
Phys. Rev. Lett., 119, 068002 (2017).

Zhang_2011
D.-B. Zhang, E. Akatyeva, and T. Dumitric ăă.
Bending ultrathin graphene at the margins of continuum mechanics.
Phys. Rev. Lett., 106:255503, 2011.

zhang2014
 Zhang T., Li X., Gao H., Defects controlled wrinkling and topological design in graphene. J. Mech. Phys. Solids, 67 (2014), 2–13.











	




