An atomistic-based foliation model for multilayer graphene materials and nanotubes

An atomistic-based foliation model for multilayer graphene materials and nanotubes Susanta Ghosh 1 , Marino Arroyo n Laboratori de C alcul Num eric (LaC aN), Departament de Matem atica Aplicada III, Universitat Polit ecnica de Catalunya-BarcelonaTech, Campus Nord UPC-C2, E-08034 Barcelona, Spain abstract We present a three dimensional continuum model for layered crystalline materials made out of weakly interacting two dimensional crystalline sheets. We specialize the model to multilayer graphene materials, including multi walled carbon nanotubes (MWCNTs). We view the material as a foliation, partitioning of space into a continuous stack of leaves, thus loosing track of the location of the individual graphene layers. The constitutive model for the bulk is derived from the atomistic interactions by appro priate kinematic assumptions, adapted to the foliation structure and mechanics. In particular, the elastic energy along the leaves of the foliation results from the bonded interactions, while the interaction energy between the walls, resulting from van der Waals forces, is parametrized with a stretch transversal to the foliation. The resulting theory is distinct from conventional anisotropic models, and can be readily discretized with finite elements. The discretization is not tied to the individual walls and allows us to coarse grain the system in all directions. Furthermore, the evaluation of the non bonded interactions becomes local. We test the accuracy of the foliation model against a previously proposed atomistic based continuum model that explicitly describes each and every wall. We find that the new model is very efficient and accurate. Furthermore, it allows us to rationalize the rippling deformation modes characteristic of thick MWCNTs, highlighting the role of the van der Waals forces and the sliding between the walls. By exercising the model with very large systems of hollow MWCNTs and suspended multilayer graphene, containing up to 10 9 atoms, we find new complex post buckling deformation patterns. 1. Introduction The mechanical behavior of carbon nanotube (CNT) systems has attracted great attention over the last decade because of their importance in nanoscience and nanotechnology. Recently, graphene in the form of single sheets or multiple layer structures has taken over much of the research activity. Both families of systems exhibit phenomenal and tightly coupled mechanical, thermal and electronic properties with unique geometries and a high degree of crystalline order. Their slenderness, together with the strength and perfection of the graphene wall, allows these structures to undergo very large deformations in a reversible manner, without alterations of the lattice structure (Chopra et al., 1995;Poncharal et al., 1999; n Corresponding author. Tel.: þ34 934011805; fax: þ34 934011825. E-mail address: [email protected] (M. Arroyo). 1 Present address: Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI 48109, USA. Yu et al., 2000a;Tombler et al., 2000;Yap et al., 2007;Yakobson et al., 1996;Maiti, 2000). This resilience to undergo severe deformations without significant damage is the basis of CNT materials (Zhang et al., 2004;Cao et al., 2005). The basic element of these materials is the graphene sheet, a two dimensional arrangement of carbon atoms forming a strong and stable honeycomb lattice, flexible in bending yet very stiff in plane, which interacts with neighboring graphene sheets and substrates through weak van der Waals interactions. The interplay between these different mechanical aspects of MWCNTs results in complex deformation morphologies (Poncharal et al., 1999;Kuzumaki et al., 1998;Lourie et al., 1998), a characteristic power law nonlinear elastic response, and a strong size effect in the post buckling regime (Arias and Arroyo, 2008). The van der Waals interactions have been shown to result in a very smooth effective potential between two neighboring walls, which offers nearly no resistance in sliding (Kolmogorov and Crespi, 2000;Yu et al., 2000b). Although the effect of registry between adhered graphene walls has been experimentally observed (Yu et al., 2001), it is in general negligible (Lu et al., 2007). The mechanical coherence, hence the overall stiffness, of multiple wall structures can be significantly enhanced by creating covalent bridges under moderate irradiation (Huhtala et al., 2004;Kis et al., 2004;Peng et al., 2008;Locascio et al., 2009). Given the phenomenal computational requirements of all atom models of CNT materials and devices, and given their shell like behavior beyond the linear elastic regime first noted by Yakobson et al. (1996), these systems have been often understood, e.g. in interpreting experimental observations, in terms of continuum mechanics models. Continuum shell theories have provided valuable predictions (Pantano et al., 2004b), and also prompted long debates about the value of phenomenological parameters, such as the thickness of a thin shell modeling a two dimensional arrangement of atoms (see Huang et al., 2006 for a recent overview on the topic and Arroyo and Belytschko, 2004a;Lu et al., 2009 for derivations of the graphene in plane and bending stiffness from atomistics, without invoking any thickness). On the other hand, a number of continuum models constructed from atomistic models have been proposed for the graphene wall, a few of which are reviewed next. In Arroyo and Belytschko (2002), a surface model (i.e. without the notion of a shell thickness) based on interatomic potentials and on an extended Cauchy Born rule was presented, which accounts for the bending stiffness of the graphene wall and includes a continuum formulation of the van der Waals interactions. This theory has been systematically tested against all atom simulations in the full nonlinear regime, beyond geometric instabilities (Arroyo and Belytschko, 2004b), and has made possible multimillion atom simulations with physically sound atomistic models of MWCNTs at an affordable computational cost (Arroyo and Belytschko, 2003;Arias and Arroyo, 2008;Zou et al., 2009). Huang and co workers presented a Cauchy Born continuum model accounting for the in plane behavior (Zhang et al., 2002b) with applications to defect nucleation (Zhang et al., 2002a) amongst others, and more recently proposed a model also accounting for the bending stiffness (Wu et al., 2008). A mixed continuum atomistic approach has been proposed by Park et al. (2006).Yang and Weinan (2006),Guo et al. (2006) and Sun and Liew (2008) proposed alternative modifications of the Cauchy Born rule to account for curvature. In a different line of thought, Dumitrica and James (2007) coarsened all atom descriptions of objective structures, such as nanotubes, through parameterizations of molecular models in terms of group transformations. See Ericksen (2008) for a review of these and other related approaches. In atomistic models of MWCNTs and multiple layer graphene using empirical potentials, the evaluation of the van der Waals energy and forces takes by far most of the computational effort. While the bonded interactions are short ranged and very stable, the non bonded interactions are longer ranged and require the frequent update of the neighbor lists. In our previous atomistic based continuum approach for the van der Waals interactions (Arroyo and Belytschko, 2002,2004b), the resulting model for the wall to wall interactions was non local, i.e. expressed in terms of a double surface integral over the interacting walls, or twice over a given wall to capture self contact. Other contact models with a molecular basis (Sauer and Li, 2007) suffer from the same difficulty. When discretized by numerical quadrature, these models still require frequent neighbor searches and the density of the quadrature points needs to resolve the length scale of the van der Waals potential (about 1 nm) (Arroyo and Belytschko, 2004b). While this continuum theory allows us to safely reduce the number of degrees of freedom by two orders of magnitude when discretized with finite elements, the number of quadrature points needed to sample the continuum version of the van der Waals interactions is only marginally smaller than the total number of atoms. This limitation in the degree of coarse graining is the bottleneck of large scale simulations with this method, which have reached MWCNTs containing nominally some 10 7 atoms on supercomputing platforms (Arias and Arroyo, 2008). While considering a phenomenological continuum model for the bonded interactions, Pantano et al. (2004a) partially alleviated this aspect in using an effective pressure separation cohesive law to model the inter wall interactions derived from a Lennard Jones atom atom potential. A similar cohesive approach has been proposed by Lu et al. (2007). Yet, in all these approaches, each and every wall needs to be explicitly described. Bulk elastic orthotropic models for MWCNTs or multilayer graphene sheets avoid the need for an individual treatment of each wall and have qualitatively reproduced deformation patterns of MWCNTs observed in experiments, but at the expense of crude approximations (Liu et al., 2001;Garg et al., 2007). Here, we present a bulk model for multilayer graphene materials and nanotubes, built from standard interatomic potentials for the bonded and the van der Waals interactions between carbon atoms. It fills the gap between detailed descriptions such as all atom or explicit layer models and effective rod like mesoscopic models for CNTs (Zhigilei et al., 2005;Buehler, 2006;Arroyo and Arias, 2008). The proposed model has no constitutive phenomenological input, beyond the particular interatomic potential of choice, and only relies on phenomenology with regards to the kinematics. Indeed, it is motivated by the observation that the spacing between neighboring walls, even in highly deformed MWCNTs, varies smoothly in space, which suggests that the density of the graphene walls in the direction perpendicular to the graphene walls can serve as a meaningful continuum strain measure (see Fig. 1, top and middle). We call it a foliation model because this mathematical construction provides a natural partition of the three dimensional region of space occupied by a multilayer graphene system, both from a structural and from a mechanical point of view. Indeed, the leaves of the two dimensional foliation represent the graphene sheets, their stretch and bending kinematics can be related to the deformation of the covalent bonds within each sheet, and ultimately the energy of the bonded interactions can be evaluated. On the other hand, the kinematics of the co dimension 2 (i.e. one dimensional) transversal sections of the foliation measure the local density of the graphene sheets, and therefore can serve to evaluate the van der Waals interactions with an adequately parametrized potential. As it will become clear later, the resulting model fundamentally differs from standard anisotropic constitutive models, in that it accounts for the bending elasticity of the leaves, and observes exactly the invariance of the energy when sheets slide relative to each other without changing shape or transversal density. As illustrated in careful numerical tests, the proposed model very accurately mimics our previous surface model, which requires the explicit description of each wall and has been quantitatively validated against atomistic calculations (Arroyo and Belytschko, 2004b). The model we present here is very efficient computationally, in that the level of coarse graining is oblivious to the inter wall spacing or the length scale of the van der Waals potential, and the evaluation of the van der Waals energy is local at each point of the foliation, i.e. neighbor searches are not needed. Furthermore, the natural correspondence between the foliation geometry and the structure of multi layer graphene systems allows us to easily incorporate other ingredients such as the mechanical enhancement due to a distribution of covalent bonds bridging neighboring walls (Peng et al., 2008;Locascio et al., 2009;Huang et al., 2010). We briefly discuss preliminary calculations including this effect. Our model is applicable to other multi layer systems formed by elastic layers with bending elasticity that interact weakly, such as graphene oxide paper (Dikin et al., 2007). See also Novoselov (2011) for a catalog of new 2D crystals, which can form by stacking 3D heterostructures. The model looses track of the precise location of the individual graphene layers, and rather Fig. 1. Experimental evidence of the smooth variation of the inter-wall spacing for a severely deformed 45-walled bent CNT (Poncharal et al., 1999) (top), cross-sections from a numerical simulation of the rippling phenomenon with the model proposed in Arroyo and Belytschko (2002,2004b) for a 40-walled CNT (middle), and computational example of a local delamination in a bent 4-walled CNT (bottom). From P. Poncharal, Z.L. Wang, D. Ugarte, W.A. de Heer, Electrostatic de ections and electromechanical resonances of carbon nanotubes, Science 283 (1999) 1513–1516. Reprinted with permission from AAAS. considers a continuous distribution of leaves. Consequently, it may not be able to describe in detail situations in which the spacing between walls does not vary smoothly, such as local delaminations (see Fig. 1, bottom). Section 2 describes the kinematics of the foliated continuum, with a special attention to the incompatible reference configuration required to describe MWCNTs, and the a non standard strain measure, the transversal stretch, measuring the packing of the walls. Section 3 describes the formulation of the constitutive model from the bonded and non bonded interactions and Section 4 provides a concise description of the numerical implementation of the theory. A number of validation numerical examples are shown in Section 5. This section also provides new insight on the mechanics of rippling in MWCNTs, and presents applications of the model in graphene and MWCNT systems containing hundreds of millions of atoms. A summary and conclusions are given in Section 6. 2. Kinematics of the foliation model A foliation is a geometric construction used to study manifolds (Cohen, 1974;Rovenski, 1995;Candell and Conlon, 1999), particularly convenient to describe layered materials. A pdimensional foliation of a ndimensional manifold Mcan be thought of as a collection of disjoint connected pdimensional submanifolds M a (the leaves of the foliation) of M. The parameter labeling the leaves parameterizes the transversal sections of the foliation. Thus, locally, a foliation looks like a decomposition of a manifold in a set of parallel submanifolds of smaller dimension. In our context, the graphene walls in a multi layer graphene or CNTs can be understood as a particular finite collection of leaves of a two dimensional foliation of the region in three space occupied by the material. These leaves are indeed connected (through the network of atoms) and disjoint (by the Pauli exclusion principle), but are discretely placed in space with a specific spacing in a certain relaxed state. The foliated continuum is viewed as a continuous distribution of graphene layers in the transversal direction. The partition of the solid provided by the foliation structure allows us to incorporate the model in Arroyo and Belytschko (2002) to describe mechanics of the bonded interactions within the leaves, and the interlayer coarse grained model for the van der Waals interactions from Girifalco et al. (2000). 2.1. Reference and deformed configurations Let Vbe the parametric body, described by Cartesian coordinates f x 1 , x 2 , x 3 g. The coordinate x 3 parameterizes the transversal sections of the foliation. It is natural to require that the parametric body is delimited by two planar faces along this direction, for instance at x 3 ¼0 and x 3 ¼1. To model multilayered graphene, we will typically consider domains of the type V¼Sð0,1Þ, where Sis a two dimensional referential domain for the leaves. For MWCNTs, Vwill be the unit cube ð0,1Þ 3 periodic in the x 1 direction, corresponding to the circumferential coordinate. A configuration map j is an injective mapping from Vto Euclidean space R 3 , that we describe with Cartesian coordinates fx 1 ,x 2 ,x 3 ,g. It is customary to introduce a reference configuration map j 0 in order to define the deformation map F ¼ jJj 1 0 , and ultimately evaluate the strain measures, e.g. the deformation gradient F¼Dx j Dx j 1 0 . When possible, it is convenient to choose j 0 such that it describes a stress free state of the body. When considering multilayer graphene systems, there is no difficulty in proceeding this way. For MWCNTs, the choice of a reference configuration is not so straightforward. A straight cylindrical configuration seems a natural choice, but it is certainly not stress free. We find it more convenient to define a fully relaxed reference configuration, which is nevertheless not compatible, as described below. Note that the reference configuration is in fact a local concept, since the map j 0 only enters the formulation through its Jacobian matrix Dx j 0 . One can thus ignore the map altogether and define a tensor field T 0 ð x Þthat takes the role of Dx j 0 , which we simply call reference configuration. We choose Cartesian basis vectors fI 1 ,I 2 ,I 3 gto describe the local reference body at each material point, where again I 3 points in the transversal direction, I 1 in the circumferential direction, and I 2 in the longitudinal direction. The atomic structure of MWCNTs is not known in precise terms, and probably depends on the method of synthesis. It is typically assumed that adjacent tubes have nominal radii that differ by a length close to the equilibrium spacing of graphene sheets t 0 0:3415 nm. This is accomplished for instance by choosing arm chair tubes (Saito et al., 1992) with indices (n,n)and (nþ5,nþ5), since the radius difference is D R¼15a 0 =2 p ¼0:339 nm for a value of the graphene equilibrium bond length of a 0 ¼0:142 nm. MWCNTs defined in this manner can in principle be placed in an ideal configuration, in which each graphene sheet is uniformly bent into a cylinder, but not sheared or stretched in plane, and the spacing between adjacent tubes is nearly t 0 . Hence, such a configuration stores bending energy, but nearly no van der Walls or in plane energy. This is not exactly the case due to the small coupling between bending and in plane deformations (Arroyo, 2003) and the dependence of the effective van der Waals interactions on curvature (Lu et al., 2007). Both effects are only noticeable for tubes with very small radius, and we will ignore them in our definition of the reference configuration T 0 ð x Þ. Let us consider a n wall walled CNT with nominal length Land tube indices ðnþ5i,nþ5iÞ,i¼0,...,n wall 1. For the purpose of homogenizing the discrete set of layers into a continuum, we assign a thickness t 0 to each graphene wall (it is not the thickness of a thin shell). Then, the nominal inner and outer radii of the continuum are obtained as R inn ¼ð3a 0 =2 p Þn t 0 =2and R out ¼ð3a 0 =2 p Þ½nþ5ðn wall 1Þþt 0 =2, respectively. See Fig. 2 for an illustration. For a fixed value of x 3 , i.e. a given leaf of the foliation, the domain ð0,1Þð0,1Þis mapped into a cylinder of perimeter 2 p ½R inn x 3 þR out ð1 x 3 Þ and length Lin the ideal configuration of the MWCNT. On the other hand, the fibers along x 3 ,ofunitlengthinV, are mapped into radial line segments of length R out R inn in the ideal configuration. Hence, distributing these deformations homogeneously, the reference configuration characterizing the stretch and shear free states of the system is T 0 ð x Þ¼ 2 p ½R inn x 3 þR out ð1 x 3 Þ 00 0L0 00R out R inn 2 6 6 43 7 7 5 :ð1Þ This tensor field provides a reference to properly define the metric deformations, e.g. the right Cauchy Green tensor, of a configuration map j :V/ R 3 . See Fig. 2 for an illustration. It is clear that this field is not the Jacobian matrix of any configuration map, since r xT 0 ¼ 02 p ðR inn R out Þ0 000 000 2 6 43 7 5a0: For this reason, it can be viewed as an incompatible reference configuration, similar to what is done in other contexts, such as plasticity (Acharya and Bassani, 2000) or tissue growth (Grikipati, 2009). Despite we cannot define a reference body V 0 ,itis possible to define the volume element of the reference configuration as dV 0 ¼det T 0 d x 1 d x 2 d x 3 :ð2Þ 2.2. Strain measures We can now define the deformation gradient of a deformation map as F¼Dx j T 1 0 :ð3Þ The right Cauchy Green tensor follows simply from C¼F T F. We define next the specific strain measures required for our purposes, quantifying the in plane deformation of the leaves, their curvature, and the density of the stack of graphene walls. The tangent vectors to the leaves of the foliation are given by g a ¼@ j @ xa , a ¼1,2:ð4Þ From here on, Greek indices will run from 1 to 2. The metric tensor of each leaf can be expressed in its natural basis fg 1 ,g 2 g as g a b¼/g a ,gbS, where /,Sdenotes the Euclidean scalar product in R 3 , and JJits associated norm. The field of normal unit vectors to the leaves of the foliation is given by n¼g 1 g 2 Jg 1 g 2 J:ð5Þ The curvature of each leaf is characterized by the components of the Weingarten map k a b¼/n,gb , a S, where the comma denotes partial differentiation with respect to xa . We define the component of the reference configuration tensor field tangent to the leaves as T l 0 ð x Þ¼ 2 p ½R inn x 3 þR out ð1 x 3 Þ 0 0L "# :ð6Þ Fig. 2. Illustration of the incompatible local reference configuration of the foliation model applied to MWCNTs. We can then define the strain measures required to evaluate the continuum bonded potential of the graphene walls, i.e. the right Cauchy Green tensor of the leaves C l ¼ðT l 0 Þ T gðT l 0 Þ 1 ð7Þ and the pull back of the Weingarten map, measuring curvature K¼ðT l 0 Þ T kðT l 0 Þ 1 :ð8Þ See Arroyo and Belytschko (2002,2004b) for details. At each material point, labelled by its coordinates in the parametric body, x , the stretch of material fibers aligned along the transversal direction can be computed as l 3 ¼I 3 CI 3 p¼C 33 p:ð9Þ Note that, while this measure of deformation may be adequate to characterize the stretch of covalent bridges between adjacent walls, it does not measure the distance between the graphene sheets after deformation, which is the figure of merit in evaluating the van der Waals interactions. See Fig. 3 for an illustration, where l 3 41 (a transversal material fiber is stretched) but l t o1 (the leaves are compacted). The transversal stretch characterizing the compaction of the graphene wall is instead l t ¼/n,FI 3 S:ð10Þ This strain measure is crucial in the current formulation and significantly differs from the strain measures used in standard anisotropic elastic models. In particular, l t is invariant with respect to the shear occurring in the continuum when the graphene sheets slide relative to each other without deforming or changing their spacing. For post processing, it is useful to define the shear deformations measuring the relative sliding of the walls. We first note that the component of FI 3 tangential to the leaves is t¼FI 3 l t n. The (nonstandard) shear deformations along the leaves of the foliation are then defined as g l 1 ¼t, 1 l 1 FI 1  , g l 2 ¼t, 1 l 2 FI 2  measuring the relative sliding of the graphene walls along the circumferential and longitudinal directions, respectively. 3. Atomistic-based strain energy for the foliated bulk continuum The constitutive model for the foliated bulk continuum results from two ingredients: (1) the surface model for the leaves, accounting for the bonded interactions and (2) an effective inter surface (rather than inter atomic) model for the van der Waals interactions. We describe below each of these ingredients, the resulting model, and briefly outline the derived stress measures appearing in the first variation of the total energy, required in the numerical implementation in Section 4. 3.1. Constitutive model for the bonded interactions of graphene Graphene is a multi lattice, whose in plane elastic behavior follows from standard finite crystal elasticity based on the Cauchy Born rule, with the local relaxation of the inner displacements (see Arroyo and Belytschko, 2002;Zhang et al., 2002b and references therein). Such continuum models are written explicitly in terms of the atomistic analytical potential of choice and inherit the anisotropy of the underlying lattice. As outlined in the Introduction, a number of approaches have been proposed to extend the Cauchy Born method to curved low dimensional lattices such as graphene. We follow here the method based on the exponential Cauchy Born rule, which links the atomic deformation measures bond lengths, bond bending angles, or possibly dihedral angles (Lu et al., 2009) and the continuum strain measures of the equivalent continuum surface model, here C l and Kin each leaf. Then, by averaging over the unit cell the discrete energy of the inequivalent bonds and angles computed with an interatomic potential, e.g. the potential for hydrocarbons in Brenner (1990), and minimizing out the relative shift between the two simple lattices of graphene, an expression is obtained for the strain energy density per unit reference surface area W ECB ðC l ,KÞ. We refer to Arroyo and Belytschko (2002,2004b) for details. A much simpler surface model was proposed in Arroyo and Belytschko (2005), which was shown to provide reasonably accurate results as compared to all atom simulations even for very large deformations. This model is based on the Kirchhoff Saint Venant strain energy density and on the Helfrich curvature energy W KSVH ðC l ,KÞ¼ 1 2 ½2 m E l :E l þ l ðtr E l Þ 2 þ 1 2 ½c a ð2HÞ 2 þc b K,ð11Þ t0 λ3t0λtt0 Fig. 3. Illustration of the standard stretch in the transversal direction l 3 and the unconventional stretch l t measuring the compaction of the leaves of the foliation. where E l ¼1=2ðC l IdÞ,H¼1=2ðC l Þ 1 :Kis the mean curvature, and K¼det½ðC l Þ 1 Kis the Gaussian curvature. As for the elastic constants, m and l are the two dimensional Lame ´coefficients (Arroyo and Belytschko, 2004a), and c a and c b are bending elastic moduli. By considering the uniform bending of a graphene sheet into a cylindrical surface, the modulus c a can be derived from atomistics (Lu et al., 2009). Although graphene materials are generally not closed surfaces, hence the integral of the Gaussian curvature is in general not zero, the term involving c b was ignored in Arroyo and Belytschko (2005). This model captures well the elastic moduli close to equilibrium and the geometric nonlinearities, which play an important role in the CNT mechanics. However, the material nonlinearities stemming from the interatomic potentials are ignored in the Kirchhoff Saint Venant Helfrich model. We present this model because in many situations its simplicity may outweigh the missing physics, although in this work the model based on the exponential Cauchy Born rule, W ECB ðC l ,KÞ, is adopted. In the following, we simply denote this potential by WðC l ,KÞ. Note that, in the foliated bulk continuum, the surface bonded energy density WðC l ,KÞcan be evaluated at any point in the body, which is made out of a continuous distribution of leaves. We can then distribute the bonded energy density per unit reference area in the thickness assigned to each sheet t 0 (the equilibrium spacing of graphene sheets) to obtain the volumetric bonded energy density WðC l ,KÞ¼ 1 t 0 WðC l ,KÞ:ð12Þ It is certainly possible to calculate the usual stress measures, such as the second Piola Kirchhoff stress tensor S¼2@W=@C, which can be assigned a precise mechanical meaning. However, in the calculation of the first variation of the energy, and in the calculation of the out of balance forces of the discretized model, the relevant stress measures are a leaf second Piola Kirchhoff stress tensor S l ¼2@W=@C l and a couple stress associated to the bending of the leaves m¼@W=@K. The calculation of these stresses in terms of the inter atomic bond order potential used here is quite involved. We refer the reader to Arroyo and Belytschko (2004b, Appendix). However, for the potential W KSVH , these can be explicitly computed as t 0 S l ¼2 m E l þ l trðE l ÞId 4c a HðC l Þ 1 KðC l Þ 1 c b KðC l Þ 1 ð13Þ and t 0 m¼2c a HðC l Þ 1 þc b 2KK 1 :ð14Þ 3.2. Constitutive model for the van der Waals interactions The van der Waals interactions give rise to very weak attractive forces when two interacting atoms are far apart, and to strong repulsive forces as the atoms are brought together. These interactions are important in the mechanics of MWCNTs, or single walled CNTs in their collapsed configurations, as well as in bundles, forests or ropes of CNTs. These interactions are generally modeled with a pair wise potential between non bonded atoms, such as the classical Lennard Jones (LJ) potential (Lennard Jones, 1924;Girifalco et al., 2000) given by V LJ ðrÞ¼4 E 2 1=6 s r ! 12 2 1=6 s r ! 6 2 43 5,ð15Þ where s denotes the distance at which the potential attains its minimum, E is the potential well depth, and ris the distance between the interacting atoms. The values of the constants for graphitic systems are E ¼0:153 10 2 aJ and s ¼0:3834 nm (Girifalco and Lad, 1956;Qian et al., 2002). In a full atom simulation, the non bonded energy follows from the sum over all the atom pairs of the system (excluding bonded pairs of atoms) P nb ¼X i X j4i V LJ ðr ij Þ,ð16Þ where r ij is the distance between the ith and the jthatoms.Thisdoublesumisaveryexpensivepartofthecomputation,and is in practice alleviated by truncating the LJ potential beyond a cut off distance r cut . To avoid a discontinuity in the potential and its derivative, the truncated potential is shifted linearly as V LJ,cut ðrÞ¼ 0ifr4r cut , V LJ ðrÞV LJ ðr cut ÞV 0 LJ ðr cut Þðr r cut Þotherwise: (ð17Þ The second sum in Eq. (16) is then performed only for neighboring atoms within a distance r cut to the ith atom. Furthermore, the neighbor list is updated periodically but not in every evaluation of the energy. In spite of these techniques, these interactions remain non local in nature, and there has been an interest for effective potentials between graphitic systems, e.g. binding energies between pairs of objects such as bucky balls, straight infinitely long CNTs and planar graphene sheets, as a function of the distance between the objects and possibly their relative pose (Girifalco et al., 2000). Here, we are interested in the effective potential (energy per unit area) between two parallel infinite graphene sheets interacting with the potential V LJ as a function of their distance t,whichisfoundtobe V g ðtÞ¼ 9V g ðt 0 Þ9 0:6 t 0 t  4 0:4t 0 t  10 "# ,ð18Þ where t 0 is the equilibrium spacing between the graphene walls and 9V g ðt 0 Þ9denotes the well depth. Their numerical values are 3.415 ˚ Aand15:36 meV=˚ A 2 , respectively. The well depth is given by 9V g ðt 0 Þ9¼2 pr 2 c E ð2 1=6 s Þ 2 ,where r c is the area density of the carbon atoms in the equilibrium configuration of graphene. For later reference, we provide also the effective potential when the atoms in each graphene sheet interact through V LJ,cut . Such a potential is not provided in the literature. Using a cut off for the LJ potential does not change in any way the computational cost of the foliation model. Nevertheless, the effective potential V g,cut is needed for a quantitative comparison with the explicit all layer model, which uses a cut off to alleviate the computational cost. Consider two infinite parallel graphene sheets at a distance t. The energy per unit area is V g cut ðtÞ¼ r c Z A V LJ,cut ðrÞ r c dA ¼2 pr 2 c Z 1 0 V LJ,cut ðR 2 þt 2 qÞRdR,ð19Þ which, after a change of integration variable, becomes V g cut ðtÞ¼2 pr 2 c Z 1 t V LJ,cut ðxÞdx ¼2 pr 2 c Z r cut t V LJ,cut ðxÞdx:ð20Þ Recalling V LJ,cut from Eq. (17), we obtain V g cut ðtÞ¼ 0ift4r cut , V g ðtÞV g ðr cut Þþ2 pr 2 c V LJ ðr cut Þðt 2 r 2 cut Þ=2 otherwise: þ2 pr 2 c V 0 LJ ðr cut Þðt 3 =3r cut t 2 =2þr 3 cut =6Þ 8 > > < > > : ð21Þ Note that this formula is not the standard cut off of Eq. (18) We provide a numerical validation of the effective graphene graphene van der Waals potential, by comparing it with an all atom calculation with the pair wise potential. In the discrete calculations, two large parallel graphene sheets in their ground state are considered, separated by a variable distance. Then, the pair wise interaction energy between one atom in one of the graphene sheets and all the atoms in the other sheet is computed, and divided by the area per atom. This energy depends on the particular registry between the two graphene sheets. For each separation distance, this effect is removed by averaging the position of the atom of the first sheet over the triangle formed by the high symmetry points of the unit cell. The excellent agreement between discrete and the continuum potentials is shown in Fig. 4. We now define the volumetric potential for the van der Waals interactions. For this, we note the following assumptions and remarks:  We assume that, locally, the graphene sheets remain nearly parallel. Consequently, the local van der Waals energy density can be computed assuming the graphene sheets are infinitely large and planar using the potential V g . The distance between a nominal pair of graphene sheets at each point of the foliation is taken as l t t 0 .  When homogenizing the graphene pair potential V g , we should account for the fact that a n wall walled CNT has in fact n wall 1 interactions between adjacent walls. 0.3 0.4 0.5 0.6 0.7 0.8 −0.2 −0.15 −0.1 −0.05 Graphene−graphene seperation (nm) vdW energy density (aJ/nm 2 ) Discrete Potential Continuum Potential Fig. 4. Graphene–graphene interaction potential densities obtained with the discrete pairwise potential and with the continuum model by Girifalco and Lad (1956).  The potential V g measures the van der Waals energy per unit area for graphene sheets in their ground state, i.e. with the equilibrium density of atoms per unit area. If the deformation stretches the graphene walls, then the number of atoms per unit area changes, and therefore, the effective potential should be corrected accordingly. With these considerations in mind, the volumetric effective local potential for the van der Waals interactions is V vdw ð l t ,C l Þ¼ ðn wall 1Þ n wall 1 t 0 J l V g ð l t t 0 Þ,ð22Þ where J l denotes the leaf Jacobian determinant J l ¼det C l q¼det g p det T l 0 ð23Þ measuring the changes in local area in the leaves. In practice, changes in area are energetically very unfavorable, and therefore, the correction that this term introduces is insignificant in many cases. This model disregards the effect on the effective van der Waals interactions of the curvature of the graphene sheets, the registry of the graphene walls, and the axial finiteness of actual CNTs. These effects have been shown to be insignificant by Lu et al. (2007). Note also that, contrary to the assumption that adjacent walls remain parallel, any interesting deformation morphology such as rippling necessarily requires the separation between adjacent walls to be non uniform, although smoothly varying (see Fig. 1). Our model tacitly assumes that the nonlocal effects in the van der Waals interactions, e.g. arising from gradients in the separation, are not very important. We check the validity of this assumption a posteriori,by direct comparison of the equilibria of the foliation model and those of the explicit layer model that treats the van der Waals interactions in a fully non local manner. It is conceivable that for other multilayer materials this assumption may lead to inaccuracies, but for multilayer graphene and MWCNTs, we find it is adequate. Further theoretical work estimating the modeling error would be very interesting. Other non local effects such as the edge effects in the tube tube interactions, responsible for MWCNT based oscillators (Zheng and Jiang, 2002) are also ignored by the model. Since this model accounts only for the interaction between neighboring walls, it does not capture the self interaction of the inner and outer shells of a MWCNT. This effect can be easily included by considering the explicit layer model in Arroyo and Belytschko (2004b) for the outer surfaces of the foliation model. The explicit layer and the foliation models can also be combined to study inhomogeneous delamination or exfoliation. The relevant stress measures in the first variation of the total energy are, one the one hand, the pressure conjugate to the transversal stretch P t,vdw ¼@V vdw @ l t ¼ðn wall 1Þ n wall 1 J l V 0 g ð l t t 0 Þ ð24Þ and on the other hand, the leaf second Piola Kirchhoff stress tensor due to the van der Waals interactions S l,vdw ¼2@V vdw @C l ¼ðn wall 1Þ n wall 1 t 0 J l V g ð l t t 0 ÞðC l Þ 1 ,ð25Þ where we have used the fact that @1=J l @C l ¼1 2J l ðC l Þ 1 : 3.3. Total potential energy and its variation The total strain energy per unit reference volume for the foliated continuum is simply WðC l ,KÞþV vdw ð l t ,C l Þ, and therefore, the total internal energy of the system is P int ½ j ¼Z V ½WðC l ,KÞþV vdw ð l t ,C l Þ dV 0 ,ð26Þ where Eq. (2) needs to be kept in mind. The application of external forces on the nuclei, for instance electrostatic forces, gives rise to a body force in the continuum. The corresponding potential is given by P ext ½ j ¼Z V /B, j SdV 0 ,ð27Þ where Bis the body force per unit reference volume. The total potential energy is then P ½ j ¼ P int ½ j  P ext ½ j :ð28Þ The (meta )stable equilibrium deformation configurations are the (local) minimizers of the total energy j ¼arg inf c 2 G P ½ c ,ð29Þ 5.4. Hollow MWCNTs Depending on the synthesis method and the conditions, MWCNTs with a large hollow core can be produced, and are common in applications. It is reasonable to expect that their mechanical response to deformation will lie between the behavior of single walled CNTs and that of thick MWCNTs without a seizable empty core. For instance, in bending, single walled CNTs have been shown to develop sharp kinks that absorb in a single location most of the deformation, whereas we have seen in previous examples that thick MWCNTs, which lack internal space to accommodate a single deep buckle, present a distributed buckling pattern. Hollow MWCNTs are quite floppy, which makes it difficult to simulate them. For instance, it is difficult to obtain converged solutions with the surface model describing each wall. The proposed foliation model can access these systems very easily. We report here the four point bending and twisting of micron sized tubes, see Fig. 11. The MWCNT subject to bending contains nominally 1.23 million atoms, while the one in torsion contains a little over 10 9 atoms. In bending, it can be seen as initially, the system develops two kinks in the central region. As bending becomes more pronounced, these kinks lock and more buckles develop. The resulting deformation morphology does not follow a well defined pattern, and different simulations with slightly different parameters can produce morphologically distinct deformations, yet with similar energy (Arroyo and Belytschko, 2004b). In torsion, we find the system switches between modes, from a fourfold symmetric cross section, to a threefold mode, and finally a twofold mode. The walls are quite floppy at the length scale of the large diameter of this tube, which leaves space for a secondary bucking pattern in the helical furrows, highlighted in the color plot of the longitudinal shear tangential to the leaves. The frequency of the secondary buckling motif is highest for the four petals mode and lowest for the two petals mode. Deformation patterns in MWCNTs arise from a competition between in plane, bending and van der Waals energies, and a full understanding has been missing. From a general consideration, we expect that the secondary motif allows the system to further relax in plane stretching or shearing at the expense of a small amount of curvature and van der Waals energy. In this particular example, we check that the ridges are under slight longitudinal tension, while the furrows are under slight longitudinal compression. This slack may explain the secondary buckling motif. 5.5. Nanoindentation of multi layer graphene The mechanical behavior of graphene has been investigated through atomic force microscopy (AFM) (Lee et al., 2008, 2009). In these experiments, flakes of monolayer to up to three layer graphene were deposited over a Si substrate with x1 x2 t l Fig. 12. Wrinkling of a 10-walled graphene square sample with a side length of 1 m m suspended over a circular hole of 0:5 m m in diameter and indented. The system contains about 380 million atoms. Top: deformation morphology at several indentation depths. Bottom: color maps of the shear tangential to the graphene walls along the x 1 direction (left) and transverse stretch, measuring the packing of the sheets (right), both at the top surface of the 10-walled graphene sample. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) circular holes from 1 to 1:5 m m in diameter. Subsequently, the suspended graphene membranes were indented. Similar experiments were performed by Poot and van der Zant (2008) for thicker multi walled (8 100 layers) graphene suspended over circular holes of radii from 84 nm to 0:54 m m. However, little is known about the detailed deformation patterns and the likely geometric instabilities in such experiments. We exercise here the foliation model to study the indentation of multilayer graphene of experimentally accessible dimensions. Accurate simulations tend to consider much smaller samples (Wang et al., 2009), a compromising assumption for the size dependent mechanics of such systems. The foliation model allows us to easily access the actual scales of experiments. We model the substrate as a flat surface interacting with graphene via a Lennard Jones potential. As proposed by Aitken and Huang (2010), the interaction between graphene and a flat SiO 2 substrate surface is modeled by the following 9 3 Lennard Jones energy area density V gsi ðzÞ¼ E gsi 3 2 h 0 z  3 1 2 h 0 z  9 "# ,ð32Þ where zis the distance between SiO 2 substrate surface and its nearest graphene mono layer. h 0 and E gsi are the equilibrium separation and equilibrium energy area density, respectively. We take h 0 ¼0:6 nm following Aitken and Huang (2010) and E gsi ¼0:096 J=m 2 , as suggested by Ishigami et al. (2007). The indenter is modeled as a flat disc like surface with 10 nm radius interacting with the uppermost graphene layer through a Lennard Jones 9 3 potential. We consider a 10 walled square graphene sample of lateral dimension 1 m m. The sample is suspended over a hole of radius 0:25 m m, and then deformed with an indenter vertically displaced at its center. For small indentation depths, the graphene sample is stretched and bent. Up to this point the membrane tension in the central part is not sufficient to overcome the adhesion between graphene and the substrate. The deformed shape is very smooth, except below the indenter. As the indentation depth increases, a critical point is reached and the graphene sample wrinkles as material slips from the adhered part into the free standing region. The amplitude of wrinkles increases with the indentation depth. See Fig. 12 for an illustration. Graphene wrinkling under indentation has been reported in molecular simulations for much smaller systems (Wang et al., 2009). Wrinkling of graphene sheets under indentation is not surprising, and similar phenomena have been observed in graphene in other situations, such as uniaxial tension (Bao et al., 2009). Given the fact that wrinkling of mono and multi layered graphene strongly influences the electronic and transport properties (Low et al., 2011;Katsnelson and Geim, 2008;Morozov et al., 2008), this simulation suggests that nano electro mechanical devices may exploit the controlled rippling of graphene systems under indentation. Our preliminary observation suggests that the onset and morphology of wrinkling is sensitive to the strength of adhesion between graphene and SiO 2 , the number of layers of graphene, and radius of the hole. A systematic study is the subject of current research. The color plot of the shear tangential to the walls g l 1 at the uppermost surface shows that, again, these wrinkling deformations produce significant sliding between the graphene walls, and may be interpreted as an observable manifestation of the smooth tangential potential between the walls. A similar pattern of wall sliding is observed at internal graphene walls, or at the bottom surface. These results also suggest that moderate irradiation, producing covalent bridges between the walls, may delay the emergence of wrinkling. The color plot of l t at the topmost surface shows that, except in the vicinity of the indenter, the separation between the walls remains close to the equilibrium separation. At internal planes or at the bottom surface, l t is nearly one everywhere. Thus, contrary to the wrinkling patterns in bending and twisting, this pattern is very relaxed from the point of view of the van der Waals forces. 6. Conclusions We have presented a new bulk atomistic based continuum model for layered crystalline materials made out of two dimensional crystalline sheets. Such systems are emerging as a new family of materials with tunable and exceptional properties (Novoselov, 2011), but here we particularize the model to multi layer graphene systems, including multi walled carbon nanotubes. The proposed model exploits the foliation structure and mechanics of these systems, which allows us to express the strain energy of the leaves (walls) in terms of the bonded atomistic interactions and the interaction energy between the leaves in terms of the van der Waals potential. The key feature of the model is that, contrary to previous atomistic based continuum models of graphene, it deals with a continuous distribution of walls in the transverse direction, rather than tracking each individual wall. This allows us to coarse grain the system in all directions in the finite element discretization of the theory. Another salient feature of the foliation model is that the non bonded interactions are parametrized in terms of a local strain measure, thus avoiding expensive neighbor searches and summations. We have tested the accuracy and efficiency of the proposed foliation model against simulations based on an explicit surface model previously validated against full atomistics, obtaining excellent results. We have exercised the model in large hollow MWCNTs and graphene samples containing nominally up to 10 9 atoms, showing the ability of the proposed model to access the length scales of actual experiments. The foliation model provides insight on the commonly observed rippling deformations of MWCNTs, highlighting the importance of a network of compressive struts due to van der Waals forces, which stabilize the post buckling motifs. We have also shown the very large relative sliding between the graphene walls that occur as a result of wrinkling in MWCNTs and graphene systems. The wrinkling morphologies can thus be interpreted as an observable signature of the smooth inter layer potential against relative sliding, and can presumably be delayed by covalent bridges between the walls. Current research includes the systematic treatment of these covalent bridges, including the possibility of bond breaking. Acknowledgment We acknowledge the support of the European Research Council (FP7/2007 2013)/ERC grant Agreement No. 240487. S.G. acknowledges the support of the Spanish Ministry of Science and Innovation through the Juan de la Cierva program. M.A. acknowledges the support received through the prize ‘‘ICREA Academia’’ for excellence in research, funded by the Generalitat de Catalunya. References Acharya, A., Bassani, J., 2000. Lattice incompatibility and a gradient theory of crystal plasticity. J. Mech. Phys. Solids 48, 1565–1595. Aitken, Z.H., Huang, R., 2010. Effects of mismatch strain and substrate surface corrugation on morphology of supported monolayer graphene. J. Appl. Phys. 107 (12), 123531. Arias, I., Arroyo, M., 2008. Size-dependent nonlinear elastic scaling of multiwalled carbon nanotubes. Phys. Rev. Lett. 100 (085503). Arroyo, M., 2003. Finite Crystal Elasticity of Curved Monolayer Lattices: Applications to Carbon Nanotubes. Ph.D. Thesis, Northwestern University. Arroyo, M., Arias, I., 2008. Rippling and a phase-transforming mesoscopic model for multiwalled carbon nanotubes. J. Mech. Phys. Solids 56, 1224–1244. Arroyo, M., Belytschko, T., 2002. An atomistic-based finite deformation membrane for single layer crystalline films. J. Mech. Phys. Solids 50 (9), 1941–1977. Arroyo, M., Belytschko, T., 2003. Nonlinear mechanical response and rippling of thick multi-walled carbon nanotubes. Phys. Rev. Lett. 91 (21), 215505. Arroyo, M., Belytschko, T., 2004a. Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy–Born rule. Phys. Rev. B 69, 115415. Arroyo, M., Belytschko, T., 2004b. Finite element methods for the non-linear mechanics of crystalline sheets and nanotubes. Int. J. Numer. Methods Eng. 59 (3), 419–456. Arroyo, M., Belytschko, T., 2005. Continuum mechanics modeling and simulation of carbon nanotubes. Meccanica 40, 455–469. Bao, W., Miao, F., Chen, Z., Zhang, H., Jang, W., Dames, C., Lau, C.N., 2009. Controlled ripple texturing of suspended graphene and ultrathin graphite membranes. Nat. Nanotechnol. 4 (9), 562–566. Brenner, D.W., 1990. Empirical potential for hydrocarbons for use in simulating chemical vapor deposition of diamond films. Phys. Rev. B 42 (15), 9458–9471. Buehler, M., 2006. Mesoscale modeling of mechanics of carbon nanotubes: self-assembly, self-folding and fracture. J. Mater. Res. 21 (11), 2855–2869. Candell, A., Conlon, L., 1999. Foliations I (Graduate Studies in Mathematics). American Mathematical Society. Cao, A., Dickrell, P.L., Sawyer, W.G., Ghasemi-Nejhad, M.N., Ajayan, P.M., 2005. Super-compressible foamlike carbon nanotube films. Science 310, 1307–1310. Chang, Y.-C., Liaw, Y.-H., Huang, Y.-S., Hsu, T., Chang, C.-S., Tsong, T.-T., 2008. In situ tailoring and manipulation of carbon nanotubes. Small 4, 2195–2198. http://dx.doi.org/10.1002/smll.200800563. Chopra, N.G., Benedict, L.X., Crespi, V.H., Cohen, M.L., Louie, S.G., Zettl, A., 1995. Fully collapsed carbon nanotubes. Nature 377, 135–138. Cohen, M., 1974. Foliations of 3-manifolds. Am. Math. Mon. 81 (5), 462–473. Dikin, D., Stankovich, S., Zimney, E., Piner, R., Dommett, G., Evmenenko, G., Nguyen, S., Ruoff, R., 2007. Preparation and characterization of graphene oxide paper. Nature 448, 457–460. Dumitrica, T., James, R., 2007. Objective molecular dynamics. J. Mech. Phys. Solids 55, 2206–2236. Ericksen, J., 2008. On the Cauchy–Born rule. Math. Mech. Solids 13, 199–220. Garg, M., Pantano, A., Boyce, M., 2007. An equivalent orthotropic representation of the nonlinear elastic behavior of multiwalled carbon nanotubes. J. Eng. Mater. Technol. 129, 431. Gilbert, J., Nocedal, J., 1992. Global convergence properties of conjugate gradient methods for optimization. SIAM J. Optim. 2 (1). Girifalco, L., Lad, R., 1956. Energy of cohesion, compressibility and the potential energy functions of the graphite system. J. Chem. Phys. 25 (4), 693–697. Girifalco, L.A., Hodak, M., Lee, R.S., 2000. Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential. Phys. Rev. B 62 (19), 13104–13110. Grikipati, K., 2009. The kinematics of biological growth. Appl. Mech. Rev. 62, 030801. Guo, X., Wang, J., Zhang, H., 2006. Mechanical properties of single-walled carbon nanotubes based on higher order Cauchy–Born rule mechanical properties of single-walled carbon nanotubes based on higher order Cauchy–Born rule. Int. J. Solids Struct. 43, 1276–1290. Huang, X., Yuan, H., Liang, W., Zhang, S., 2010. Mechanical properties and deformation morphologies of covalently bridged multi-walled carbon nanotubes: multiscale modeling. J. Mech. Phys. Solids 58, 1847–1862. Huang, Y., Wu, J., Hwang, K., 2006. Thickness of graphene and single-wall carbon nanotubes. Phys. Rev. B 74, 245413. Huhtala, M., Krasheninnikov, A., Aittoniemi, J., Stuart, S., Nordlund, K., Kaski, K., 2004. Improved mechanical load transfer between shells of multiwalled carbon nanotubes. Phys. Rev. B 70, 045404. Ishigami, M., Chen, J.H., Cullen, W.G., Fuhrer, M.S., Williams, E.D., 2007. Atomic structure of graphene on SiO 2 . Nano Lett. 7 (6), 1643–1648. Katsnelson, M., Geim, A., 2008. Electron scattering on microscopic corrugations in graphene. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 366 (1863), 195–204. Kis, A., Csa ´nyi, G., Salvetat, J.P., Lee, T.-N., Couteau, E., Kulik, A.J., Benoit, W., Brugger, J., Forro ´, L., 2004. Reinforcement of single-walled carbon nanotube bundles by intertube bridging. Nat. Mater. 3, 153–157. Kolmogorov, A.N., Crespi, V.H., 2000. Smoothest bearings: interlayer sliding in multiwalled carbon nanotubes. Phys. Rev. Lett. 85 (22), 4727–4730. Kuzumaki, T., Hayashi, T., Ichinose, H., Miyazawa, K., Ito, K., Ishida, Y., 1998. In-situ observed deformation of carbon nanotubes. Philos. Mag. A 77 (6), 1461–1469. Lee, C., Wei, X., Kysar, J.W., Hone, J., 2008. Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science 321 (5887), 385–388. Lee, C., Wei, X., Li, Q., Carpick, R., Kysar, J.W., Hone, J., 2009. Elastic and frictional properties of graphene. Phys. Status Solidi (B) 246 (11–12), 2562–2567. Lennard-Jones, J., 1924. On the determination of molecular fields. II. From the equation of state of a gas. Proc. R. Soc. Lond. A 106 (738), 463–477. Liu, D., Nocedal, J., 1989. On the limited memory method for large scale optimization. Math. Program. B 45 (3), 503–528. Liu, J., Zheng, Q., Jiang, Q., 2001. Effect of a rippling mode on resonances of carbon nanotubes. Phys. Rev. Lett. 86 (21), 4843–4846. Locascio, M., Peng, B., Zapol, P., Zhu, Y., Li, S., Belytschko, T., Espinosa, H., 2009. Tailoring the load carrying capacity of MWCNTs through inter-shell atomic bridging. Exp. Mech. 49 (2), 169–182. Lourie, O., Cox, D.M., Wagner, H.D., 1998. Buckling and collapse of embedded carbon nanotubes. Phys. Rev. Lett. 81 (8), 1638–1641. Low, T., Guinea, F., Katsnelson, M.I., 2011. Gaps tunable by electrostatic gates in strained graphene. Phys. Rev. B 83, 195436. Lu, Q., Arroyo, M., Huang, R., 2009. Elastic bending modulus of monolayer graphene. J. Phys. D 42 (10), 102002. Lu, W., Wu, J., Jiang, L., Huang, Y., Hwang, K., Liu, B., 2007. A cohesive law for multi-wall carbon nanotubes. Philos. Mag. 87 (14), 2221–2232. Maiti, A., 2000. Mechanical deformation in carbon nanotubes-bent tubes vs tubes pushed by atomically sharp tips. Chem. Phys. Lett. 331, 21–25. Morozov, S.V., Novoselov, K.S., Katsnelson, M.I., Schedin, F., Elias, D.C., Jaszczak, J.A., Geim, A.K., 2008. Giant intrinsic carrier mobilities in graphene and its bilayer. Phys. Rev. Lett. 100, 016602. Novoselov, K., 2011. Nobel lecture: Graphene: materials in the flatland. Rev. Mod. Phys. 83, 837–849. Pantano, A., Parks, D., Boyce, M., 2004a. Mechanics of deformation of single and multi-wall carbon nanotubes. J. Mech. Phys. Solids 52 (4), 789–821. Pantano, A., Parks, D., Boyce, M., Nardelli, M., 2004b. Mixed finite element-tight-binding electromechanical analysis of carbon nanotubes. J. Appl. Phys. 96 (11), 6756–6760. Park, J., Cho, Y., Kim, S., Jun, S., Im, S., 2006. A quasicontinuum method for deformations of carbon nanotubes. Comput. Model. Eng. Sci. 11, 61–72. Peng, B., Locascio, M., Zapol, P., Li, S., Mielke, S., Schatz, G., Espinosa, H., 2008. Measurements of near-ultimate strength for multiwalled carbon nanotubes and irradiation-induced crosslinking improvements. Nat. Nanotechnol. 3, 626–631. Piegl, L., Tiller, W., 1997. The NURBS Book. Springer Verlag. Poncharal, P., Wang, Z.L., Ugarte, D., de Heer, W.A., 1999. Electrostatic deflections and electromechanical resonances of carbon nanotubes. Science 283, 1513–1516. Poot, M., van der Zant, H.S.J., 2008. Nanomechanical properties of few-layer graphene membranes. Appl. Phys. Lett. 92 (6), 063111. Qian, D., Wagner, G., Liu, W., Yu, M., Ruoff, R., 2002. Mechanics of carbon nanotubes. Appl. Mech. Rev. 55 (6), 495–553. Rovenski, V., 1995. Foliations on Riemannian Manifolds and Submanifolds. Birkh¨ auser. ISBN: 978-0-8176-3806-1. Saito, R., Fujita, M., Dresselhaus, G., Dresselhaus, M., 1992. Electronic structure of chiral graphene tubules. Appl. Phys. Lett. 60 (18), 2204–2206. Sauer, R., Li, S., 2007. A contact mechanics model for quasi-continua. Int. J. Numer. Methods Eng. 71, 931–962. Sen, D., Novoselov, K., Reis, P., Buehler, M., 2010. Tearing graphene sheets from adhesive substrates produces tapered nanoribbons. Small 6, 1108–1116. Sun, Y., Liew, K., 2008. Application of the higher-order Cauchy–Born rule in mesh-free continuum and multiscale simulation of carbon nanotubes. Int. J. Numer. Methods Eng. 75, 1238–1258. Tombler, T., Zhou, C., Alexseyev, L., Kong, J., Dai, H., Liu, L., Jayanthi, C., Tang, M., Wu, S., 2000. Reversible electromechanical characteristics of carbon nanotubes under local-probe manipulation. Nature 405, 769–772. Wang, C.Y., Mylvaganam, K., Zhang, L.C., 2009. Wrinkling of monolayer graphene: a study by molecular dynamics and continuum plate theory. Phys. Rev. B 80 (15), 155445. Wu, J., Hwang, K., Huang, Y., 2008. An atomistic-based finite-deformation shell theory for single-wall carbon nanotubes. J. Mech. Phys. Solids 56, 279–292. Yakobson, B.I., Brabec, C.J., Bernholc, J., 1996. Nanomechanics of carbon tubes: instabilities beyond the linear response. Phys. Rev. Lett. 76 (14), 2511–2514. Yang, J., Weinan, E, 2006. Generalized Cauchy–Born rules for elastic deformation of sheets, plates, and rods: derivation of continuum models from atomistic models. Phys. Rev. B 74, 184110. Yap, H., Lakes, R., Carpick, R., 2007. Mechanical instabilities of individual multiwalled carbon nanotubes under cyclic axial compression. Nano Lett. 7 (5), 1149–1154. Yu, M., Dyer, M., Chen, J., Qian, D., Liu, W., Ruoff, R., 2001. Locked twist in multiwalled carbon-nanotube ribbons. Phys. Rev. B 64 (24). Yu, M., Lourie, O., Dyer, M., Moloni, K., Kelly, T., Ruoff, R., 2000a. Strength and breaking mechanism of multiwalled carbon nanotubes under tensile load. Science 287, 637–640. Yu, M.-F., Yakobson, B.I., Ruoff, R.S., 2000b. Controlled sliding and pullout of nested shells in individual multiwalled carbon nanotubes. J. Phys. Chem. B 104 (37), 8764–8767. Zhang, D., Akatyeva, E., Dumitrica, T., 2011. Bending ultrathin graphene at the margins of continuum mechanics. Phys. Rev. Lett. 106, 255503. Zhang, M., Atkinson, K.R., Baughman, R.H., 2004. Multifunctional carbon nanotube yarns by downsizing an ancient technology. Science 306 (5700), 1358–1361. Zhang, P., Huang, Y., Gao, H., Hwang, K., 2002a. Fracture nucleation in single-wall carbon nanotubes under tension: a continuum analysis incorporating interatomic potentials. J. Appl. Mech. 69, 454–458. Zhang, P., Huang, Y., Geubelle, P., Klein, P., Hwang, K., 2002b. The elastic modulus of single-wall carbon nanotubes: a continuum analysis incorporating interatomic potentials. Int. J. Solids Struct. 39 (13–14), 3893–3906. Zheng, Q., Jiang, Q., 2002. Multiwalled carbon nanotubes as gigahertz oscillators. Phys. Rev. Lett. 88, 045503. Zhigilei, L., Wei, C., Srivastava, D., 2005. Mesoscopic model for dynamic simulations of carbon nanotubes. Phys. Rev. B 71 (16), 165417. Zou, J., Huang, X., Arroyo, M., Zhang, S.L., 2009. Effective coarse-grained simulations of super-thick multi-walled carbon nanotubes under torsion. J. Appl. Phys. 105, 033516.