The Gaussian Stiffness of Graphene deduced from a Continuum Model based on Molecular Dynamics Potentials
Cesare Davini, Antonino Favata, Roberto Paroni
TL;DR
This paper derives an analytical expression for the Gaussian stiffness of graphene using a continuum model based on molecular dynamics potentials, clarifying its atomic-scale origins and resolving previous conflicting values.
Contribution
It introduces the first analytical formula for graphene's Gaussian stiffness derived from a continuum limit of a molecular dynamics model.
Findings
Analytical expression for Gaussian stiffness of graphene.
Quantitative evaluation of atomic-scale sources of stiffness.
Clarification of previously conflicting values in literature.
Abstract
We consider a discrete model of a graphene sheet with atomic interactions governed by a harmonic approximation of the 2nd-generation Brenner potential that depends on bond lengths, bond angles, and two types of dihedral angles. A continuum limit is then deduced that fully describes the bending behavior. In particular, we deduce for the first time an analytical expression of the Gaussian stiffness, a scarcely investigated parameter ruling the rippling of graphene, for which contradictory values have been proposed in the literature. We disclose the atomic-scale sources of both bending and Gaussian stiffnesses and provide for them quantitative evaluations.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
**The Gaussian Stiffness of Graphene
deduced from a Continuum Model
based on Molecular Dynamics Potentials**
Cesare Davini1
Antonino Favata2
Roberto Paroni3
1 Via Parenzo 17, 33100 Udine
2 Department of Structural and Geotechnical Engineering
Sapienza University of Rome, Rome, Italy
3 Dipartimento di Architettura, Design e Urbanistica
University of Sassari, Alghero (SS), Italy
Abstract
We consider a discrete model of a graphene sheet with atomic interactions governed by a harmonic approximation of the 2nd-generation Brenner potential that depends on bond lengths, bond angles, and two types of dihedral angles. A continuum limit is then deduced that fully describes the bending behavior. In particular, we deduce for the first time an analytical expression of the Gaussian stiffness, a scarcely investigated parameter ruling the rippling of graphene, for which contradictory values have been proposed in the literature. We disclose the atomic-scale sources of both bending and Gaussian stiffnesses and provide for them quantitative evaluations.
Keywords: Graphene, Continuum Modeling, Gaussian stiffness.
Contents
1 Introduction
Graphene has attracted increasing interest during the past few years, and is nowadays used in a great variety of applications, taking advantage of its extraordinary mechanical, electrical and thermal conductivity properties. Nevertheless, its potentialities, and those of graphene-based materials, are far from being fully explored and exploited, and many studies are carried out by the scientific community in order to develop new technological applications [17].
The understanding of the bending behavior of graphene is of paramount importance in several technological applications. It is exploited, for example, to predict the performance of graphene nano-electro-mechanical devices and ripple formation [21, 28, 41, 25, 35, 22, 34, 20, 19, 38, 30, 15], and it is proposed to be the key point to produce efficient hydrogen-storage devices [36, 18, 37]. A very recent review in Materials Today by Deng & Berry [12] gives an overview on the hot problem of wrinkling, rippling and crumpling of graphene, highlighting formation mechanism and applications. Indeed, these corrugations can modify its electronic structure, create polarized carrier puddles, induce pseudo-magnetic field in bilayers and alter surface properties. Although a great effort has been done on the experimental side, predictive models are still wanted. They are of crucial importance when these phenomena need to be controlled and designed.
In particular, since the bending stiffness and the Gaussian stiffness —that is the reluctance to form non-null Gaussian curvatures— are the two crucial parameters governing the rippling of graphene, it is necessary to accurately determine them for both the design and the manipulation of graphene morphology. Although several evaluations of the bending stiffness have been proposed in the literature, the Gaussian stiffness has not been object of an intensive study. Indeed, as pointed out in a very recent review on mechanical properties of graphene [1], only two conflicting evaluations have been proposed. In [23], periodic boundary conditions have been used within a quantum-mechanical framework, and the value eV has been found. While in [38] the estimate of eV has been obtained by combining the configurational energy of membranes determined by Helfrich Hamiltonian with energies of fullerenes and single wall carbon nanotubes calculated by Density Functional Theory (DFT). At a discrete level there are two main difficulties in the evaluation of the Gaussian stiffness: on the one hand, controlling a discrete double curvature surfaces is problematic, and on the other hand, a suitable notion of Gaussian curvature at the discrete level should be introduced. Instead, when well established continuum models are adopted, such as plate theory, one has the problem of determining the equivalent stiffnesses, letting alone the conceptual crux of giving a meaning to the notion of thickness (see [21], [5] and references therein).
In this paper we deduce a continuum 2-dimensional model of a graphene sheet inferred from Molecular Dynamics (MD). In particular, looking at the 2nd-generation reactive empirical bond-order (REBO) potential [6], we give a nano-scale description of the atomic interactions and then we deduce the continuum limit, avoiding the problem of postulating an “equivalent thickness” and circumventing artificial procedures to identify the material parameters that describe the mechanical response of a plate within the classical theory.
Our analysis of the atomic-scale interaction relies on the discrete mechanical model proposed in [14] and exploited in [16, 15, 2], whose results are also based on the 2nd-generation Brenner potential. This potential is largely used in MD simulations for carbon allotropes; for a detailed description of its general form and that adopted in our theory we refer the reader to Appendix B of [14]. Here, we recall the key ingredients needed:
- (i)
the kinematic variables associated with the interatomic bonds involve first, second and third nearest neighbors of any given atom. In particular, the kinematical variables we consider are bond lengths, bond angles, and dihedral angles; from [6] it results that these latter are of two kinds, that we here term C and Z, as carefully described in Sec. 2. 2. (ii)
graphene suffers an angular self-stress, and the self-energy associated with the self-stress (sometimes called cohesive energy in the literature) is quantitatively relevant; 3. (iii)
the energetic contribution of dihedral interaction is very relevant in bending.
For the first time, we propose a continuum model able to predict both the bending and the Gaussian stiffnesses. The analytical formula we obtain for the former predicts exactly the same value as that computed with MD simulations of the last generation. The value of the Gaussian stiffness we obtain is in very good agreement with DFT computations proposed in [38].
For the modeling of graphene many different approaches at different scales can be found in the literature, ranging from first principle calculations [24, 26], atomistic calculations [40, 42, 31] and continuum mechanics [7, 39, 29, 33, 32, 9, 8]. Furthermore, mixed atomistic formulations with finite elements have been reported for graphene [3, 4].
The paper is organized as follows. In Sec. 2, we describe the kinematics and the energetic of the graphene sheet at the nano-scale. In Sec. 3, we deduce the strain measures for the change of edge lengths, wedge angles and dihedral angles, approximated to the lowest order that makes the energy quadratic in the displacement. In Sec. 4, the total energy is split in its in-plane and out-of-plane contributions, and focus is set on the latter, having the first already been considered in [10]. In Sec. 5, we deduce a continuous energy that approximates the discrete energy and in Sec. 6 the limit energy is rearranged in a more amenable form, able to put in evidence the equivalence with plate theory. In Sec. 7, quantitative results for the continuum material parameters are deduced, by means of the 2nd-generation Brenner potential, and compared with the literature. Appendix A, containing some computations ancillary to Sec. 3, completes the paper.
2 Description of kinematics and energetics of the graphene sheet
At the nano-scale a graphene sheet is a discrete set of carbon atoms that, in the absence of external forces, sit at the vertices of a periodic array of hexagonal cells. More specifically, atoms occupy the nodes of the –lattice, see Figure 1, generated by two simple Bravais lattices
[TABLE]
simply shifted with respect to one another. In (1), denotes the lattice size (the reference interatomic distance), while \ell\mathbf{d}_{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 267\relax}}{\mbox{\boldmath\textstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptstyle\mathchar 267\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 267\relax}}}} and respectively are the lattice vectors and the shift vector, whose Cartesian components are given by
[TABLE]
The sides of the hexagonal cells in Figure 1 stand for the bonds between pairs of next nearest neighbor atoms and are represented by the vectors
[TABLE]
For convenience we also set
[TABLE]
As reference configuration we take the set of points contained in a bounded open set .
Graphene mechanics is ruled by the interactions between the carbon atoms given by some suitable potential. According to the 2nd-generation Brenner potential [6], as detailed in [14, 16], in order to account properly for the mechanical behavior of a bended graphene sheet it is necessary to consider three types of energetic contributions, respectively coming from: binary interactions between next nearest atoms (edge bonds), three-bodies interactions between consecutive pairs of next nearest atoms (wedge bonds) and four-bodies 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. 2).
We consider a harmonic approximation of the stored energy and assume that it is given by the sum of the following terms:
[TABLE]
, \mathcal{U}_{\ell}^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}} and are the energies of the edge bonds, the wedge bonds and the dihedral bonds, respectively; denotes the distance between nearest neighbor atoms, the angle between pairs of edges having a lattice point in common and and the Z- and C-dihedral angles between two consecutive wedges, to be defined later (see Fig. 3); is the edge length at ease, {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}^{\mbox{\scriptsize nat}} the angle at ease between consecutive edges and the dihedral angle at ease.
The sums extend to all edges, , all wedges, , all Z-dihedra, , and all C-dihedra, , contained in the set . The bond constants , k^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}, , and will be deduced by making use of the 2nd-generation Brenner potential.
The graphene sheet does not have a configuration at ease (i.e. stress-free). Indeed, in [14] it has been shown that
[TABLE]
where {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{0}\neq 0. We set {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\Theta:=\Theta, l=\ell+{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}l and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}=\frac{2}{3}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}+{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}} and write (4) as
[TABLE]
In particular, up to a constant, the wedge energy takes the form
[TABLE]
with
[TABLE]
the angle self-stress. The dihedral bonds play an important role because they contribute to the stored energy by about 50%, see [14, 15], the rest is due to the angle self-stress {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 284\relax}}{\mbox{\boldmath\textstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 284\relax}}}_{0} associated to the wedge bonds.
The energy decomposition (5) is based on the choice of the set of kinematical variables \{l,{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}},\operatorname{\overset{(z)}{\Theta}},\operatorname{\overset{(c)}{\Theta}}\}. This choice is the most natural, if one considers the 2nd-generation Brenner potential, where all those variables appear in explicit manner. A harmonic approximation in each of those parameters is of course unique.
In the next section we shall make explicit the change of length {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}l, the change of wedge angle \textstyle\mathchar 270\relax$$\textstyle\mathchar 291\relax, and the changes of the Z- and C-dihedral angles {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(z)}{\Theta}} and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(c)}{\Theta}}. In Section 4, with the notation introduced in the next section, we shall write the energies (5) more explicitly.
3 Approximated strain measures
In this section we calculate the strain measures associated to a change of configuration described by a displacement field , approximated to the lowest order that makes the energy quadratic in .
3.1 Change of the edge lengths
With {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}l_{i}(\mathbf{x}^{\ell}) we denote the change in length of the edge parallel to and starting from the lattice point . Thus,
[TABLE]
and up to terms can be rewritten as
[TABLE]
In particular, the first order changes are determined by the in-plane components of only.
3.2 Change of the wedge angles
For each fixed node we denote by {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}(\mathbf{x}^{\ell}) the angle of the wedge delimited by the edges and ; that is, the wedge angle opposite to the -th edge (see Fig. 4). Here, , and take values in and the sums should be interpreted mod 3: for instance, if then and .
From (6) we see that the change in the wedge angle enters into the energy not just quadratically but also linearly, therefore the variations of the wedge angle should be computed up to the second order approximation. To keep the notation compact, we set
[TABLE]
Let
[TABLE]
and
[TABLE]
be the images of the edges parallel to and and starting at . Then, the angle {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}(\mathbf{x}^{\ell}) is given by
[TABLE]
Calculations given in Appendix A.1 yield that
[TABLE]
where {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{(1)} and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{(2)} are the first order and the second order variation, respectively, of the wedge angle with respect to the reference angle \frac{2}{3}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}. Therefore, keeping up to second order terms one has that
[TABLE]
It turns out that the first order variation takes the form
[TABLE]
with defined by
[TABLE]
that is, the unit vector orthogonal to , (cf. equation (19) in [13]).
Figure 5 illustrates the geometrical meaning of formula (12). In particular,
[TABLE]
as it could have been deduced from geometrical considerations.
The second order variation is given by, see Appendix A.1,
[TABLE]
By algebric manipulation one finds that
[TABLE]
where denotes the out-of-plane component of the displacement, that is
[TABLE]
where is the unit vector perpendicular to the undeformed sheet. Note that, by (15), the \sum_{i}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}(\mathbf{x}^{\ell}) is non-positive and hence the contribution of the self-stress to the strain energy is non-negative for {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 284\relax}}{\mbox{\boldmath\textstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 284\relax}}}_{0}<0, i.e., for {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{0}>0, see (7).
3.3 Change of the dihedral angles
For each fixed node and for each edge parallel to and starting at we need to define four types of dihedral angles and :
[TABLE]
where
[TABLE]
are the images of vectors and (see Fig. 6, for ), parallel to and and starting at the image of the point .
Also here, , and take values in and the sums should be interpreted mod 3: for instance, if then and .
The C-dihedral angle is the angle corresponding to the C-dihedron with middle edge and oriented as , while is the angle corresponding to the C-dihedron oriented opposite to (see Fig. 6 for ). The Z-dihedral angle corresponds to the Z-dihedron with middle edge and the other two edges parallel to (see Fig. 6 for ).
Then, recalling that {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\Theta=\Theta, calculations in Appendix A.2 yield that
[TABLE]
and
[TABLE]
Analogous formulas hold for {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{i}^{-}} and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{i}\mathbf{p}_{i+2}}:
[TABLE]
4 Splitting of the energy
The above calculations show that {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{(1)} as well as {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}l_{i} depend upon the in-plane components of , cf. (8) and (12), while {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{(2)}, {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(c)}{\Theta}}, and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(z)}{\Theta}} depend upon the out-of-plane component of , cf. (15), (18), and (19). This yields a splitting of the energy into membrane and bending parts
[TABLE]
defined by
[TABLE]
where is the self-energy (corresponding to the so-called cohesive energy in the literature) and is the dihedral energy. The analysis in a paper by Davini [10] applies here to the in-plane deformations, providing a continuum model of the graphene sheet within the framework of -convergence theory. Hereafter, we concentrate on the out-of-plane deformations. With the notation introduced in Section 3 we now write the bending energy more explicitly. The self-energy can be written as
[TABLE]
where \sum_{i=1}^{3}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{(2)}(\mathbf{x}^{\ell}) is given in (15) in terms of the out-of-plane component of the displacement . We further split the dihedral energy in
[TABLE]
where
[TABLE]
is the contribution of the Z-dihedra, and
[TABLE]
is the contribution of the C-dihedra. The Z- and C-dihedral angles appearing in (38) and (37) are given in terms of in (18)-(20).
In the next section we deduce, by means of a formal analysis, a continuous version of the discrete bending energy
[TABLE]
from which we shall deduce expressions for the sheet’s bending stiffnesses. A rigorous analysis based on -convergence theory will be done in a forthcoming paper [11].
5 The continuum limit
In this section we find a continuous energy, defined over the domain , that approximates the discrete bending energy defined over the lattice . This is achieved by letting the lattice size go to zero so that invades . With this in mind, in place of a function , we consider a twice continuously differentiable function .
Given two vectors and , with we denote the second partial derivative of in the directions and , that is
[TABLE]
where denotes the Hessian of . Clearly, we also have
[TABLE]
The change of the Z-dihedra, see (20)2, can be rewritten, after setting
[TABLE]
as
[TABLE]
where the third equality follows from (26). Similarly, setting
[TABLE]
we have that
[TABLE]
Taking (3) into account, we may rewrite the vectors and , defined in (27) and (29), in terms of , for instance and , and then rewrite the Z-dihedral energy, see (38) and rewritten below for the reader convenience, as
[TABLE]
where is the area of the hexagon of side centred at (see Fig. 7).
Let {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 287\relax}}{\mbox{\boldmath\textstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptstyle\mathchar 287\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 287\relax}}}_{E^{\ell}(\mathbf{x}^{\ell})}(\mathbf{x}) be the characteristic function of , i.e., it is equal to if and [math] otherwise, and let
[TABLE]
Then, we may simply write
[TABLE]
and since converges, as goes to zero, to
[TABLE]
we deduce that
[TABLE]
The functional , defined in (31), is the continuum limit of the Z-dihedral energy.
Working in a similar manner, we find the continuum limit of the C-dihedral energy:
[TABLE]
and the continuum limit of the self-energy:
[TABLE]
Detailed calculations leading to (32) and (33) are found in Appendix A.3.
The total bending limit energy is therefore
[TABLE]
6 The equivalent plate equation
In this section we rewrite the limit energies in a more amenable form. We start by manipulating the limit C-dihedral energy. We first note that
[TABLE]
and similarly
[TABLE]
from which we find that
[TABLE]
where denotes the Laplacian of . Hence, the C-dihedral energy defined in (32) rewrites as
[TABLE]
We now tackle the Z-dihedral energy. Recalling (31), by a calculation similar to that carried on in (35) we find
[TABLE]
where , and from these equations we deduce that
[TABLE]
Similar identities hold for and . Thence, we find that the Z-dihedral energy takes the form
[TABLE]
The second sum is equal, as it can be checked, to the last line of (36), and with a similar calculation we also find that
[TABLE]
and hence
[TABLE]
We now deal with the self-energy. Again with a calculation similar to that carried on in (35) we find
[TABLE]
and hence
[TABLE]
By summing (37), (38), and (39), we find the total bending energy, defined in (34):
[TABLE]
where
[TABLE]
are the bending and the Gaussian stiffnesses, respectively. is called Gaussian because it multiplies the Gaussian curvature , while is twice the mean curvature. The analytical expression for the bending stiffness (42)1 coincides with the one deduced in [28, 15], within a discrete mechanical framework, if one assumes that . It clearly shows that the origin of the bending stiffness is twofold: a part depends on the dihedral contribution, and a part on the self-stress. Quite surprisingly, the self-stress has no role in the Gaussian stiffness. It is worth noticing that in the above approach there is no need of introducing any questionable effective thickness parameter.
7 Numerical results
In this section, we adopt the 2nd-generation Brenner potential [6] to obtain quantitative results for the continuum material parameters deduced in Sec. 6.
The 2nd-generation REBO potentials developed for hydrocarbons by Brenner et al. in [6] accommodate up to third-nearest-neighbor interactions through a bond-order function depending, in particular, on dihedral angles. Following Appendix B of [14], we here give a short account of the form of this potential.
The binding energy of an atomic aggregate is given as a sum over nearest neighbors:
[TABLE]
the interatomic potential is given by
[TABLE]
where the individual effects of the repulsion and attraction functions and , which model pair-wise interactions of atoms and depending on their distance , are modulated by the bond-order function . The repulsion and attraction functions have the following forms:
[TABLE]
where is a cutoff function limiting the range of covalent interactions, and where , , , , and , are parameters to be chosen fit to some material-specific dataset. The remaining ingredient in (44) is the bond-order function:
[TABLE]
where apexes and refer to two types of bonds: the strong covalent -bonds between atoms in one and the same given plane, and the -bonds responsible for interlayer interactions, which are perpendicular to the plane of -bonds. The role of function b_{ij}^{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 283\relax}}{\mbox{\boldmath\textstyle\mathchar 283\relax}}{\mbox{\boldmath\scriptstyle\mathchar 283\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 283\relax}}}-{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}} is to account for the local coordination of, and the bond angles relative to, atoms and ; its form is:
[TABLE]
Here, for each fixed pair of indices , (a) the cutoff function limits the interactions of atom to those with its nearest neighbors; (b) {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 277\relax}}{\mbox{\boldmath\textstyle\mathchar 277\relax}}{\mbox{\boldmath\scriptstyle\mathchar 277\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 277\relax}}}_{ijk} is a string of parameters designed to prevent attraction in some specific situations; (c) function depends on and , the numbers of and atoms that are nearest neighbors of atom ; it is meant to adjust the bond-order function according to the environment of the C atoms in one or another molecule; (d) for solid-state carbon, the values of both the string {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 277\relax}}{\mbox{\boldmath\textstyle\mathchar 277\relax}}{\mbox{\boldmath\scriptstyle\mathchar 277\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 277\relax}}}_{ijk} and the function are taken null; (e) function modulates the contribution of each nearest neighbour of atom in terms of the cosine of the angle between the and bonds; its analytic form is given by three sixth-order polynomial splines. Function b_{ij}^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}} is given a split representation:
[TABLE]
where the first addendum depends on whether the bond between atoms and has a radical character and on whether it is part of a conjugated system, while the second addendum depends on dihedral angles and has the following form:
[TABLE]
where function is a tricubic spline depending on , , and , a function of local conjugation, and the dihedral angle is defined as
[TABLE]
The values of the constant and can be deduced by deriving twice of the potential, and computing the result in the ground state (GS): , {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 274\relax}}{\mbox{\boldmath\textstyle\mathchar 274\relax}}{\mbox{\boldmath\scriptstyle\mathchar 274\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 274\relax}}}_{ijk}=2/3{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 281\relax}}{\mbox{\boldmath\textstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptstyle\mathchar 281\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 281\relax}}}, . In particular, we find:
[TABLE]
where is the value of in the GS.
Remark 1
With this notation the bending stiffness becomes:
[TABLE]
This expression coincides with that given in [28]:
[TABLE]
after noticing that
[TABLE]
Instead, in references [4] and [21] the dihedral energies are not contemplated and the bending stiffness found, up to notational differences, coincides with ours after setting .
With the values reported in [6], we get:
[TABLE]
From [14], we take the value of the selfstress {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 284\relax}}{\mbox{\boldmath\textstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptstyle\mathchar 284\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 284\relax}}}_{0}:
[TABLE]
With (52) and (53), we obtain:
[TABLE]
The value of is in complete agreement with the literature [27, 38]; from (52) and (53) it is possible to check that the contribution of the self-stress and the dihedral stiffness amounts to about % and % of the total, respectively. Neither analytical evaluation of , nor MD computations, have been proposed so far. The value we obtain is in good agreement with the value of eV, reported in [38] and determined by means of DFT.
8 Conclusions
Starting form a discrete model inferred from MD, we have deduced a continuum theory describing the bending behavior of a graphene sheet. Atomic interactions have been modeled by exploiting the main features of the 2nd-generation Brenner potential and adopting a quadratic approximation of the energy. The deduced continuum limit fully describes the bending behavior of graphene. To our knowledge, it is the first time that an analytical expression of the Gaussian stiffness is given and an explanation of its origins at the atomistic scale is provided. We also derived a quantitative evaluation of the related constitutive parameters.
Acknowledgments
AF acknowledges the financial support of Sapienza University of Rome (Progetto d’Ateneo 2016 — “Multiscale Mechanics of 2D Materials: Modeling and Applications”).
Appendix A Appendix
In order to compute the strain measures for small changes of configuration of the graphene foil, we write the displacements of the nodes in the form
[TABLE]
where is a positive scalar measuring smallness and \operatorname{\boldsymbol{\mathfrak{u}}}:=\frac{\bf u}{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}} stands for the displacement distribution normalized accordingly.
A.1 Change of the bond angle
Let us define the bond angle as
[TABLE]
with
[TABLE]
where we have set
[TABLE]
Then, from Taylor’s expansion we get
[TABLE]
where the various terms can be calculated by successive differentiations of Eq. (55). Thus,
[TABLE]
which yields
[TABLE]
where we take into account that \sin({{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}(0))=\frac{\sqrt{3}}{2}, , and .
Moreover, by differentiating Eq. (55) twice, we get
[TABLE]
which gives
[TABLE]
Computations yield that
[TABLE]
Since and \left(\frac{\mathbf{m}\cdot\mathbf{n}}{|\mathbf{m}||\mathbf{n}|}\right)^{\prime}\Big{|}_{{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}=0}=-\sin{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}(0){{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{\prime}(0), we finally have that
[TABLE]
All in all, we have that:
[TABLE]
Recalling that {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{\prime}(0)\,{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}+\frac{1}{2}\,{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{i}^{\prime\prime}(0){{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{2}+O({{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{3}) and that \mathbf{u}={{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}\operatorname{\boldsymbol{\mathfrak{u}}}, we get:
[TABLE]
with
[TABLE]
and
[TABLE]
We write (66) in the simpler form
[TABLE]
where the unit vectors ’s are defined by (13).
A.2 Change of the dihedral angle
To fix the ideas, we focus on the dihedral angles {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(c)}{\Theta}}_{\mathbf{p}_{1}^{+}}, and {{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 270\relax}}{\mbox{\boldmath\textstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptstyle\mathchar 270\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 270\relax}}}\operatorname{\overset{(z)}{\Theta}}_{\mathbf{p}_{1}\mathbf{p}_{2}}, sketched in Fig. 8; 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 denote by
[TABLE]
the edge vectors after the deformation. We have that
[TABLE]
it is easy to see that
[TABLE]
whence
[TABLE]
On recalling the identity
[TABLE]
we obtain:
[TABLE]
where we have used the fact that and set . Now,
[TABLE]
whence
[TABLE]
Thus, on recalling that w=\mathfrak{w}{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}, we conclude that
[TABLE]
For a generic C-dihedral angle centered in , we get
[TABLE]
For the Z-dihedral angle , we introduce the vector
[TABLE]
image of under the deformation. We have that
[TABLE]
and
[TABLE]
whence
[TABLE]
Again, on making use of (73) and recalling that , we obtain:
[TABLE]
On noticing that |\mathbf{q}_{5}\times\mathbf{q}_{1}|=\frac{\sqrt{3}}{2}\ell^{2}+O({{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 280\relax}}{\mbox{\boldmath\textstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptstyle\mathchar 280\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 280\relax}}}^{2}), we get
[TABLE]
whence
[TABLE]
For a generic Z-dihedral angle centered in , we get
[TABLE]
A.3 Deduction of the continuum limits and
In this Appendix we give a justification of (32) and (33) following the lines outlined in Section 5 for the derivation of the continuum limit of the Z-dihedral energy.
As in Section 5, in place of a function , we consider a twice continuously differentiable function . We first consider the contribution due to the C-dihedra. Momentarily, to keep the notation compact, we set
[TABLE]
By Taylor’s expansion we have
[TABLE]
Thence, since , we find from (18) that
[TABLE]
and, by substituting , we eventually get
[TABLE]
where is defined in (13). Similarly, we also find that
[TABLE]
With the same steps taken in the study of the Z-dihedral energy we now derive the limit of the C-dihedral energy. With the expressions of the change of the C-dihedra (87) and (88), we can rewrite the C-dihedral energy, see (37), as
[TABLE]
where the function
[TABLE]
converges to \sum_{i=1}^{3}\Big{(}\partial^{2}_{\mathbf{p}_{i}\mathbf{p}_{i}^{\perp}}w\Big{)}^{2}, as goes to zero. Thus,
[TABLE]
which is (32).
We now compute the limit of the self-energy. By Taylor’s expansion, with the notation introduced in (86), and taking into account that , we find
[TABLE]
With this expression the self-energy (22) takes the form
[TABLE]
Since the sum is over the points of both lattices, whereas in the previous cases the sum was only over the nodes of , we cannot use the hexagons earlier introduced. Let be the triangle centered at of side as depicted in Figure 9.
Let
[TABLE]
and note that the area of is . The self-energy rewrites as
[TABLE]
Since W^{{}\mathchoice{\mbox{\boldmath\displaystyle\mathchar 291\relax}}{\mbox{\boldmath\textstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptstyle\mathchar 291\relax}}{\mbox{\boldmath\scriptscriptstyle\mathchar 291\relax}}}_{\ell} converges to (\partial^{2}_{\mathbf{p}_{1}\mathbf{p}_{1}}w+\partial^{2}_{\mathbf{p}_{2}\mathbf{p}_{2}}w+\partial^{2}_{\mathbf{p}_{1}\mathbf{p}_{2}}w\Big{)}^{2} as goes to zero, we find
[TABLE]
which is (33).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] D. Akinwande, C.J. Brennan, J.S. Bunch, P. Egberts, J.R. Felts, H. Gao, R. Huang, J.-S. Kim, T. Li, Y. Li, K.M. Liechti, N. Lu, H.S. Park, E.J. Reed, P. Wang, B.I. Yakobson, T. Zhang, Y.-W. Zhang, Y. Zhou, and Zhu Y. A review on mechanics and mechanical properties of 2d materials - graphene and beyond. Preprint at ar Xiv:1609.07187 .
- 2[2] R. Alessi, A. Favata, and A. Micheletti. Pressurized CN Ts under tension: A finite-deformation lattice model. Compos. Part B Eng. , in press, http://dx.doi.org/10.1016/j.compositesb.2016.10.006.
- 3[3] M. Arroyo and T. Belytschko. An atomistic-based finite deformation membrane for single layer crystalline films. J. Mech. Phys. Solids , 50(9):1941 – 1977, 2002.
- 4[4] M. Arroyo and T. Belytschko. Finite crystal elasticity of carbon nanotubes based on the exponential cauchy-born rule. Phys. Rev. B , 69:115415, 2004.
- 5[5] C. Bajaj, A. Favata, and P. Podio-Guidugli. On a nanoscopically-informed shell theory of carbon nanotubes. Europ. J. Mech. A/Solids , 42:137–157, 2013.
- 6[6] 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):783, 2002.
- 7[7] E. Cadelano, P.L. Palla, S. Giordano, and L. Colombo. Nonlinear elasticity of monolayer graphene. Phys. Rev. Lett. , 102:235502, 2009.
- 8[8] C. Galiotis D. Sfyris, G.I. Sfyris. Curvature dependent surface energy for a free standing monolayer graphene: Some closed form solutions of the non-linear theory. Int. J. nonl. Mech. , 67:186–197, 2014.