A finite deformation membrane based on inter-atomic potentials for the transverse mechanics of nanotubes

vectors. In addition to the extraction of elastic material tensors, these models have been used with the finite element method to solve boundary value problems, such as the nano-indentation of silicon (Tadmor et al., 1999). If the deformation is not uniform enough for the local theory to hold, mixed continuum atomistic approaches have been proposed to deal with inhomogeneities or defects (Tadmor et al., 1996; Shenoy et al., 1999). The appeal of the approach based on the Born rule stems from the fact that it gives rise to finite deformation constitutive models based on the nano-scale physics rather than phenomenological material models. The inter-atomic potentials, based on experimental data and quantum mechanical considerations or calculations, are thus at the core of the resulting strain energy density. A local quasicontinuum has also been developed based on the tight-binding method (Tadmor et al., 1999). This paper deals with the application of such a theory to the mechanics of carbon nanotubes. Since the discovery of these crystalline tubes in 1991, many studies have focused on their unique mechanical properties, through experiments (Yu et al., 2000; Chopra et al., 1995; Yu et al., 2001b), molecular dynamics (MD) and molecular mechanics (MM) simulations (Bernholc et al., 1998; Gao et al., 1998) and first-principles calculations (Zhou et al., 2001; Maiti, 2000). Although molecular simulations seem well suited to study these systems, they are not completely satisfactory. Indeed, they are very demanding from the computational point of view. The length scales ( AA) and time scales (ps) that must be resolved are often well below the scales of practical interest for a particular problem. In addition, although systems of over one million nuclei are currently being analyzed, one can always envision larger problems for which the computing capabilities do not suffice. An approach that can alleviate some of the drawbacks of molecular simulations is the use of continuum mechanics. The ability of continuum models to describe the mechanics of nanotubes has been noted by several authors. Cross-section continuum models have been used to explain experimental observations on the transverse stability of nanotubes (Chopra et al., 1995; Yu et al., 2001a). These extremely simplified models bring insights as well as quantitative information on the physical phenomena that govern the stability of the circular and the collapsed configurations observed in nanotubes. The vibrational properties of carbon nanotubes have been investigated through linear elasticity by Sohlberg et al. (1998). The elastic properties of carbon nanotubes as a continuum, neglecting all curvature effects, have investigated by Lu (1997). Zhong-Can et al. (1997) considered the nanotube to be an inextensible membrane, and obtained an expression of the elastic energy in terms of the curvature for a family of simple deformations. Yakobson et al. (1996) used the theory of elastic shells and linearized bifurcation analysis to study the buckling patterns of compressed carbon nanotubes observed in MD simulations. Qian et al. (2001) used a 3D continuum theory combined with a mesh-free approximation to study C60 molecules inside nanotubes. Nevertheless, the proposed models so far are either over-simplified, restricted to the linear regime, or to very particular situations, and do not constitute a systematic continuum approach to the mechanics of nanotubes. This is particularly true with regards to the large deformations. Indeed, experiments (Chopra et al., 1995; Falvo et al., 1997), MD/MM simulations (Bernholc et al., 1998) and first-principles calculations (Maiti, 2000) show that carbon nanotubes undergo very large deformations, with highly non-linear behavior and still remain elastic in the sense that the deformations are reversible, with stable bonds and intact bond topology. From these considerations, developing a finite deformation model based on the Born rule and on nano-scale physics applicable to nanotubes would be of great interest and would fill in a gap in the present use of continuum models to model crystalline films one atom thick. Apart from the physical insights that a continuum model brings, such a description could be the basis for an efficient numerical simulation methodology, in contrast with the sometimes too detailed molecular simulations. Efficiency becomes an issue when nano-ropes bundles of tens to hundreds of nanotubes or multi-walled nanotubes several microns long (Ruoff et al., 1993; Yu et al., 2000) are 2 to be analyzed. Furthermore, a continuum mechanics theory allows us to exploit the symmetry of certain situations explicitly, analogously to the plane strain and plane stress situations in 3D elasticity. Such reduced models for the transverse behavior of nanotubes are one of the topics presented in this paper. Unfortunately, the traditional approach based on the Born rule works for bulk materials but fails to extend directly to the case of crystalline films and ropes one atom thick deforming in higher dimensional spaces, i.e. 3D in the case of films, and 2D or 3D in the case of ropes. The present paper describes some mechanical effects that arise from a recent extension of the Born rule to membranes (Arroyo and Belytschko, 2002). The extension is based on the differential geometry concept of the exponential map, and is here called the exponential Born rule. An alternative approach has been reported by Friesecke and James (2000). The outline of this paper is as follows: we first present the Born rule for bulk materials, investigate its structure and explain why it breaks down for films (Section 2). Then, after some geometric preliminaries, Section 3 introduces the proposed exponential Born rule in an abstract and general way. This abstract presentation of the theory is complemented by its realization in the simplest, yet complete, situation, i.e. an atomic chain deforming in two dimensions. The formulation of the equivalent continuum rope-like object is detailed in Section 4, and a simple example illustrating the effectiveness of this model in mimicking the atomic chain is provided. Previous to the application of the theory to carbon nanotubes, their crystalline structure, as well as the instance of inter-atomic potential considered, are described in Section 5. Section 6 describes in detail the model for the transverse mechanics of carbon nanotubes. The implementation of the new theory to the arbitrary deformation of the continuum membrane in 3D is presented in Arroyo and Belytschko (2002). The Bravais multi-lattice nature of graphene requires the treatment of additional internal variables, the so-called inner displacements. Additionally, the van der Waals interactions are also accounted for in the continuum theory, and the continuum variational statement of the problem as well as the Lagrangian stress measures that naturally arise are described. Finally, a validation test comparing the proposed continuum model discretized with finite elements to molecular calculations is provided in Section 7. Several simulations highlighting the relevance of the van der Waals forces in the transverse configurations of carbon nanotubes and nano-ropes are also included in this section. 2. Breakdown of the Born rule for films The formulation of a finite deformation continuum model for space-filling defectless crystals based on the Born rule is relatively straightforward. The Born rule links the atomistic deformation to that of the continuum medium. Then, a representative crystallite is considered, and, for a given continuum deformation, the continuum strain energy density is defined to be the energy of the crystallite subject to the deformation divided by its volume. The details of the procedure are presented in several of the articles referenced in Section 1, and will also be briefly described later in the present paper. The focus of this section is on the fundamental kinematic assumption that links the continuum and the atomic deformations, i.e. the Born rule. The details of the atomic model are deliberately omitted. Later, an instance of an atomic model is adopted. 2.1. The standard Born rule Assume for the moment that we are dealing with space-filling continuum bodies, i.e. open subsets of the ambient Euclidean space. Let Ube the deformation that maps the undeformed body X0Rn, into Rn,nbeing either 1, 2 or 3. If X denotes a point in the undeformed body, its image after deformation is x¼UðXÞ. The deformed body is denoted as X¼UðX0Þand is an open set of Rn. The deformation gradient is the derivative of the vector-valued vectorial function U,F¼ DU¼oU=oX2Rnn. At each point X, the deformation gradient is a linear transformation from Rn into Rn, which maps ‘‘infinitesimal’’ material vectors, dx¼FdX(see Malvern, 1969, p. 156). 3 From a differential geometry point of view, the deformation gradient is called the tangent map of U, and is denoted as F¼TU. Let us call the infinitesimal neighborhoods of Xand xthe tangent spaces of the undeformed body and the deformed one, respectively denoted as TXX0and TxX(see Fig. 1 for an illustration). Then, using this language, the deformation map Umaps the undeformed body into the deformed one, and the tangent map F¼TUmaps the tangent space of the undeformed body into the tangent space of the deformed body. In the absence of slips, phase transitions and other special crystallographic phenomena, the Cauchy and Born hypothesis for crystals are equivalent for homogeneous deformations (Ericksen, 1984). What is referred to as the Born rule in some works is simply called the method of homogeneous deformations in others (Martin, 1975; Cousins, 1978). The Born hypothesis consists of assuming that the lattice vectors deform as would material line elements in a homogeneous deformation: a¼FA;ð1Þ where Adenotes an undeformed lattice vector and athe same vector in the deformed crystal. The geometry of the lattice vectors, that is their length and the angles they form with other lattice vectors in the deformed crystal, can therefore be extracted from the continuum deformation through the Green deformation tensor C¼FTFusing standard continuum mechanics relations: kak¼ ACA pand cos h¼ACB kakkbk;ð2Þ where Band brepresent another undeformed and deformed lattice vector and his the angle a and bform in the deformed crystal. Once the geometry of the deformed lattice vectors is linked to the continuum deformation, a constitutive model based on the atomic interactions can be constructed by identifying the continuum strain energy density with the potential energy of the atomic system for a representative cell divided by its volume. One could argue that the rule expressed by Eq. (1) is formally inconsistent, because the lattice vectors Aand a, each connecting two atomic positions, are physical entities that lie in the undeformed and deformed body respectively, while the tangent map F¼TUmaps elements of the tangent of the undeformed body into elements of the tangent of the deformed body. This inconsistency can also be viewed from a more classic standpoint: the lattice vectors have finite length while the deformation gradient maps ‘‘infinitesimal’’ material vectors, dx¼FdX. These objections are circumvented by noting that for homogeneous deformations, Eq. (1) holds exactly, even for material vectors of finite length. The Born rule assumes that, at least ‘‘locally’’, i.e. in the scale of the lattice vectors, the deformation is homogeneous. 2.2. Why the case of films is more difficult Consider now the case in which we have a single atom thick crystalline film (such as a graphene sheet) deforming arbitrarily in 3D. It is natural in this case to treat the continuum solid as a membrane without thickness. The sheet is then a twomanifold embedded in R3(a surface). It is assumed that the atoms lie on the surface (CauchyÕs hypothesis), and therefore the lattice vectors are chords of the surface. We would like to use the Born rule in order to express the geometry of the deformed lattice vectors in terms of the some continuum variable characterizing the deformation of the surface, such as the Green deformation tensor. Suppose that the undeformed body is planar, like a planar graphene sheet. In this case, X0an open set in R2. The deformation map transforms this originally planar body into a curved memFig. 1. Deformation map and its tangent map for space filling bodies. 4 F¼g1I1þg2I2;ð7Þ where fI1;I2gis the dual basis of fI1;I2g. Thus, a vector W¼WIII2TXX0is transformed by the deformation gradient into w¼FW ¼WIgI2 TUðXÞX. Therefore, the matrix representation of F in the Cartesian/convected basis fgIIJgis the 22 identity matrix, the information about the deformation being contained in the convected basis vectors. We can also define the Green deformation tensor C¼FTF. Its components in the Cartesian basis of X0coincide with these of the metric tensor in the convected basis presented in Eq. (6). The tensor Ccan be used to measure length, angle and area changes due to the deformation in terms of undeformed body quantities, i.e. Coperates in TX0. In particular the element of area of Xcan be written in terms of the element of area of X0as dX¼JdX0where the Jacobian is J¼ det gIJ p. Also, the stretch in the direction of a unit vector E2TXX0is K¼ECðXÞE p. 3.1.2. The second fundamental form: the curvature The unit normal to the surface Xcan be defined as n¼g1g2 kg1g2k;ð8Þ where kkdenotes the Euclidean norm. The second fundamental form of the deformed body kcan be expressed in the basis fgIgJgin terms of its components: kIJ ¼ngI;J;ð9Þ where gI;Jdenotes the derivative of gIwith respect to XJ. The normal curvature knat a point xof the surface Xand in a given direction defined by the unit vector v¼vIgI2TxX, is the minimum of the curvatures of all the curves of Xpassing through xtangent to v. It can be obtained as: knðxÞ¼kIJ vIvJ:ð10Þ Suppose the normal curvature of the deformed body is to be computed at a point x¼UðXÞ2Xin a given direction V¼VIII2TXX0of the undeformed body. This direction corresponds in the deformed body to v¼FV ¼VIgI, where Eq. (7) has been used. Therefore, after normalizing v, the resulting expression for the normal curvature is: knðxÞ¼ kIJ VIVJ gMN VMVN p;ð11Þ where the denominator corresponds to the Euclidean norm of v. 3.1.3. The exponential map A simple definition of the exponential map is given in Morgan (1993) for a manifold M: The exponential map exppat a point pin M maps the tangent space TpMinto Mby sending a vector vin TpMto the point in Ma distance jvjalong the geodesic from pin the direction v. The exponential map is invertible and differentiable in a neighborhood of each regular point pof the manifold. It can be defined because of the existence and uniqueness of geodesics at any point given a direction in the tangent space. The exponential map is defined here in abstract terms because its evaluation requires the knowledge of the geodesics. In general, obtaining the geodesics involves solving the geodesic differential equations. These equations are a system of non-linear ordinary differential equations whose unknowns are the parametric coordinates of the geodesic, and whose coefficients are the Christoffel symbols of the surface. Finding the geodesics, and thus the exponential map, is much simpler in some particular cases, as will be shown for the cylinder. More details about the exponential map for surfaces can be found in do Carmo (1976). Fig. 3 provides an illustration of how the exponential map brings a tangent vector to the surface. Fig. 3. Illustration of the exponential Born rule; the geodesic at xin the direction of wis represented by a dashed line. 6 3.2. Exponential Born rule In the present theory, the continuum solid equivalent to the original single layer crystalline film is a membrane without thickness. The nuclei of the atomistic system lie on this surface and consequently the lattice vectors are chords of the surface. Let Adenote an undeformed lattice vector. In the present setting, since the undeformed body is planar, X0and TX0can be identified, and consequently Acan be transformed through the deformation gradient. The result of the transformation w¼FA is simply the deformed lattice vector we would obtain through the standard Born rule. However the vector wis tangent to the deformed surface X, not a chord. Consider the following generalized kinematic rule, called the exponential Born rule in the following: a:¼expUðXÞFA:ð12Þ Described in words, this map takes a lattice vector in the undeformed body Aemanating from Xand transforms it into a vector win the tangent of the deformed body Xat x¼UðXÞ. Then, this vector is mapped from the tangent space to the deformed surface through the exponential map, which ‘‘brings’’ the result of the standard Born rule back to the surface, hence defining a chord (see Fig. 3 for an illustration of this procedure). Thus, the exponential Born rule links the deformation of the lattice vectors to the deformation of the continuum object, since both the deformation gradient and the exponential map are defined in terms of the deformation map U. This extended kinematic rule provides a theoretical framework for the application of crystal elasticity to curved crystals, by rectifying the shortcomings of the standard Born rule. Note however that its practical implementation is not straightforward, since the evaluation of the exponential map requires the determination of the geodesics, which in general entails the integration of a system of two non-linear differential equations. This results in a computationally very complex method that is necessarily non-local. Here we present approximations to the exponential map that render the model local and computationally feasible. 4. Atomic chain in 2D In this section, we illustrate the exponential Born rule for the simplest case, an atomic chain deforming in 2D. The resulting continuum model is a hyper-elastic rope whose strain energy density depends on the stretch and the curvature of the continuum object. This constitutive model is based exclusively on the atomistic description of the chain. In this case the exponential map is approximated at each point by the exponential map of the circle, for which a closed-form expression is straightforward. 4.1. Atomic model The strain energy of the atomic system is described by means of bond stretch Vsand bond angle Vhpotentials. The strain potential energy of the atomic chain can be written as a function of the nuclear positions xi: Pchainðx1;...;xnÞ¼X nB k1 VsðakÞþX mB l1 VhðhlÞ; ð13Þ where akdenotes the bond lengths, hldenotes the angle that adjacent bonds form, and nBand mBare the number of bonds and adjacent bonds, respectively. This particular atomistic model is chosen for simplicity, but the approach is not restricted to this structure of the inter-atomic potential by any means. The exponential Born rule provides a link between the atomistic and the continuum deformations, and can be combined with any atomistic model of choice, not restricted to closest-neighbor models. 4.2. Continuum model As illustrated in Fig. 4, the undeformed body is considered to be a 1D line segment that is allowed to deform in 2D. Therefore, the deformation map can be described as x¼UðXÞ¼U1ðXÞi1þU2ðXÞi2 with X2X0Rand fi1;i2gthe basis of R2. In this case, the components of the deformation gradient are ½F¼½U1 ;X;U2 ;XT, and the Green deformation tensor Cis a scalar, whose square root is the stretch Kof the deformed rope: 7 K¼C p¼ ðU1 ;XÞ2þðU2 ;XÞ2 q:ð14Þ The normal curvature knof the deformed rope can be written as: kn¼1 K3ðU2 ;XU1 ;XX U1 ;XU2 ;XX Þ;ð15Þ and can be interpreted geometrically as the inverse of the radius of curvature of the curve. As we mentioned in Section 3.2, in order to obtain a practical method the exponential Born rule needs to be approximated. It is desirable that the approximation of the exponential Born rule leads to a local model, i.e. one in which the strain energy depends on the local deformation of the rope. The strategy followed to obtain such an approximation is to perform the exponential map at each point, not of the original curve, but of a circle of radius r¼1=knwith the same normal as the original curve (see Fig. 4). Thus, locally, this circle replaces the original curve. The exponential map of the circle is readily available in closed form. The first part of the exponential Born rule maps the lattice vector Aof length Ainto a vector tangent to the curve whose components are ½w¼A½F. Therefore, its length is w¼C pA:ð16Þ The exponential map of the circle is illustrated in Fig. 5. The length of the tangent vector wis ‘‘walked’’ on the geodesic to obtain expUðXÞw, and therefore the chord a. Since the geodesic of the circle is trivially the circle itself, the length of the arc defined by the ends of ais w. Let hdenote the angle formed by two adjacent deformed lattice vectors. Consider the triangle formed by the ends of aand the center of the circle. This triangle is isosceles, and its equal angles are h=2. Therefore its third angle, the angle subtended by the arc of length w,isc¼ph. Consequently, we can relate the length of the arc, w,to the radius of the circle r¼1=knand the angle c: w¼cr¼ðphÞ=kn:ð17Þ Since the length of the unequal side of the triangle can be easily computed as a¼kak¼2rsin c 2;ð18Þ it follows that a¼2 kn sin knw 2and h¼pknw:ð19Þ Note from Eqs. (16) and (19) that the quantities a and h, which are the arguments of the atomistic energy (see Eq. (13)), are expressed in terms of the continuum deformation. The next step is to consider a representative crystallite of the atomistic system, which in this case is a cell of length Aincluding a single nucleus in the undeformed crystal. In a homogenization process, the energy of this deformed cell containing one bond and one angle between adjacent bonds is identified to the strain energy density of the continuum multiplied by the undeformed volume of the cell: AWðUÞ¼VsðaÞþVhðhÞ. Since our aim is to formulate a hyper-elastic continuum model, the elastic potential WðUÞis a strain energy per undeformed volume, in this case undeformed Fig. 5. The exponential map of the circle defined at each point of the curve by the unit normal and the normal curvature. Fig. 4. Illustration of the continuum rope like model for an atomic chain deforming in 2D. 8 length. The continuum strain energy density depends on the deformation map Uthrough the local strain measures Cand kn. Therefore, the hyperelastic potential of the continuum rope can be written as: WðC;knÞ¼1 AVs2=knsinðknC pA=2Þ hinþVhp hknC pAio:ð20Þ The total strain energy of the continuum system approximating the atomistic energy of Eq. (13) can then be written as: PropeðUÞ¼ZX0 WðC;knÞdX0:ð21Þ By taking derivatives of the hyper-elastic potential Wwith respect to the strain measures, Lagrangian stress measures arise: the work conjugate to Cis an axial stress analogous to the second Piola Kirchhoff stress tensor, and the conjugate of knis a bending moment-like stress. Second derivatives lead to the axial, bending and coupled axialbending elastic tangent moduli. 4.3. Example and discussion Suppose an initially rectilinear undeformed rope of length nA is bent into a circle of radius r with uniform stretch. The continuum stretch is K¼2pr=ðnAÞ, and the curvature is kn¼1=r. Since the atoms are postulated to lie on the continuum surface, the corresponding equispaced atomic chain containing nbonds is deformed into a regular polygon of nsides whose circumcircle has a radius r. The bond length and angle predicted by the continuum model (see Eqs. (16) and (19)) are: a¼2rsin p nand h¼p2p=n:ð22Þ It is easy to see that these predictions coincide exactly with the actual bond lengths and angles of the atomic chain deformed into a regular polygon. Therefore, the predicted energetics for this finite deformation are also exact. Of course, for a general deformation with non-constant stretch and curvature, the local approximation of the exponential Born rule will lead to approximate energetics. The examples presented later demonstrate, however, that this approximation is very accurate. The inadequacy of the standard Born rule can be illustrated easily in the present example. The standard Born rule corresponds to taking a¼ w¼FA. Suppose that our rectilinear 1D undeformed body is deformed into a circle without stretch, i.e. K¼C¼1. The application of the standard Born rule leads to deformed lattice vectors that are tangent to the rope. Consequently, two lattice vectors emanating from the same nucleus remain collinear after deformation, so the angle they form is unchanged irrespective of the bending of the rope. Furthermore, since the rope is bent without stretch, the length of the deformed lattice vectors also remains unchanged (see Eq. (16)). Therefore, the energy of such a model will remain unchanged, and the resulting rope has zero bending stiffness. However, the real lattice vectors do not remain coplanar and their length changes due to the curvature even if K¼C¼1, since from Eq. (22) it follows that for this isometric deformation a¼nA=psin p nand h¼p2p=n:ð23Þ Therefore, the energy of the atomic system will change when deformed in this fashion. Thus, a continuum model based on the standard Born rule is blind to the fact that the rope is being bent, and assigns zero energy change to the deformation, in sharp contrast with the exponential Born rule, which predicts the correct energetics. Although an intuitive approach would associate the continuum stretch to the stretch of the bonds, and the continuum curvature to changes in bond angles, the proposed model couples these deformation modes. Indeed, the continuum bond length a, which is the argument for the inter-atomic stretch potential, depends both on Cand knin a non-linear fashion. The same applies to the continuum bond angle h. This feature is essential and makes the continuum model exact for deformations that map an initially straight chain into a circular arc with constant stretch. Thus, as in the case of the standard Born rule for bulk crystalline materials, the resulting continuum model for the rope is exact for homogeneous deformations. 9 The continuum statement of the problem of finding stable equilibrium solutions is then given by: U¼arg inf W2CPðWÞ  ;ð59Þ where Cis the appropriate space of deformations or trial functions accounting for essential boundary conditions. According to the principle of stationary energy, the equilibrium solutions of the system are stationary points of the potential energy functional, and they verify the principle of virtual work: 0¼dPðUÞ ¼ZX0 ob WW oC:dCþob WW okn dkn!dX0 ZX0 BdUdX0þdPnb;ð60Þ where dUdenotes the virtual deformation. The variations of the non-bonded continuum potential can be written as: dPnb ¼1 2 4 S2 0ZX0ZX0BX V0 nb kUðXÞUðYÞk ½UðXÞUðYÞ ½dUðXÞdUðYÞdX0YdX0X:ð61Þ Let us also define the stress measures, always evaluated at the relaxed inner displacements ^ gg. Recalling the expression of the strain energy density in Eq. (48) and following a similar rationale to that used to obtain Eq. (52), we obtain: S¼2ob WW oC¼2oW oC ¼2 S0X 3 l1 V0 s okalk oC "þ2X 3 k1 oVh oh ohk oC  þoVh or1 okaik oCþoVh or2 okajk oC#;ð62Þ and m¼ob WW okn¼oW okn ¼1 S0X 3 l1 V0 s okalk okn "þ2X 3 k1 oVh oh ohk okn  þoVh or1 okaik oknþoVh or2 okajk okn#:ð63Þ The in-plane stress Scorresponds precisely to the Second Piola Kirchhoff stress, while mis a moment-like stress. Note that, because of the special form of C(see Eq. (33)), Shas only two non-zero components, which are related to the tractions in the axial and the circumferential directions. On the other hand, mis here a scalar. The general theory for arbitrary deformations is given in Arroyo and Belytschko (2002). Note that, since the membrane has no thickness, the units of Sare force divided by length, while mis expressed in units of force (bending moment divided by length). Using the Green strain tensor E¼1=2ðCIÞ, we can rewrite the principle of virtual work as: 0¼ZX0ðS:dEþmdknÞdX0 ZX0 BdUdX0þdPnb:ð64Þ Depending on the treatment of the axial stretch K1 (see Eq. (30)) different situations can be studied: Plane strain: We can consider the situation in which the value of K1is prescribed. In this case, the unknowns of the variational problem (64) are U2and U3(see Eq. (30)). If K1¼1, a deformation analogous to plane strain conditions is achieved. This applies to very long or axially constrained nanotubes. K1can also be prescribed an arbitrary value to study the transverse behavior of stretched or compressed nanotubes; the behavior will change due to the non-linearity of the model. In this situation, dK1¼0 and the axial component of the membrane stress does not appear in the variational principle. This means that the axial stress can be computed a posteriori, but does not play a role in the solution of the problem. Plane stress: Alternatively, the axial component of the membrane stress can be prescribed, for in16 stance, to be zero. This would be the case of axially unconstrained nanotubes. In this case, in addition to U2and U3, the axial stretch K1becomes an unknown of the problem. 7. Validation and representative simulations In this section, numerical simulations of stable configurations of carbon nanotubes in different situations are reported. The reduced continuum model described in the previous section is used and the variational principle described in Eq. (64) is discretized by Galerkin finite elements (FE). Thus, the original discrete molecular system is replaced by a continuum model which is subsequently transformed by the FE method into another discrete system. However, in principle we are free to design the FE discretization so that the FE model has fewer degrees of freedom than the original system. Furthermore, since the continuum model is 2D, while the full atomistic model is 3D, the computational cost is further reduced. First, the exponential Born rule-based continuum model is validated by comparing FE simulations based on it with full atomistic calculations. In these comparisons the inter-atomic potentials used in the MM simulations are used to construct the continuum constitutive equation, and analogous boundary conditions are considered in both calculations. Since the continuum model is intended to mimic the atomistic system, which is viewed as ‘‘true’’, the term error should be understood as deviation form the atomistic model. Simulations show that the agreement is excellent with regard to the energetics as well as to the stable configurations. Simulations of a model based on the standard Born rule are also provided, illustrating the deficiencies of such a model. The exponential Born rule simulations also show that, for the tested situations, the relaxation of the inner displacements greatly affects the energetics but has very little impact on the minimum energy configurations. The continuum model is then applied to simulate several situations where the transverse behavior of carbon nanotubes and the effect of van der Waals interactions are important. A final example of the generalization of the model to three dimensions is presented, with a twisting test of a [10,10] nanotube beyond the point of structural instabilities. The inter-atomic potentials fall into the general form described in Eq. (28). The two-body potential Vsis a Morse potential while the three-body potential depends only on the angle Vhand is harmonic with a sextic correction. The parameters are taken from the MM2 model. The non-bonded interactions are based on the classical LennardJones (6 12) potential. The variational principle in Eq. (64) imposes restrictions on the finite element interpolation spaces. The virtual internal work term involves variations on the curvature of the test functions, and therefore the finite element space needs to be H2, i.e. have up to second square integrable derivatives. This is why C1Hermite finite elements are chosen. Note that the discretization of the configuration described in Eq. (30) requires the approximation of the scalar functions U2ðÞ and U3ðÞ, i.e. the curve in R2described by these functions needs to be parametrized with respect to the finite element degrees of freedom. Each of these functions is approximated by piecewise C1 cubic polynomials, and therefore, each node I carries four degrees of freedom: U2 I,U3 I,ðU2Þ0 Iand ðU3Þ0 I. The internal and external work terms of the variational principle are integrated using 3 Gauss points per element, while the integration of the non-bonded interactions term may require more integration points depending on the size of the finite elements relative to the van der Waals equilibrium distance. Four integration points are required for this term in some of the simulations. The BFGS quasi-Newton technique is used both in the relaxation of the inner displacements and in the global energy minimization. This iterative method only requires gradients of the objective function and approximates the inverse of its Hessian using information from the previous iterations. For some of the larger examples involving more than one nanotube, and when the initial configuration is very far from equilibrium, dynamical relaxation is used to obtain a good first guess which is further refined with the BFGS minimization algorithm. 17 7.1. Validation test To validate the proposed reduced continuum model, a FE discretized version is compared to a MM model. A [32,0] zigzag carbon nanotube (the standard description of carbon nanotubes in terms of two integers is described by Saito et al. (1992)) is considered (see Fig. 9(a)). The molecular model used in the comparison has 384 nuclear positions, that is 1152 degrees of freedom, while the FE model has 20 nodes and consequently 80 degrees of freedom. Note that the discrete FE model reduces the computational cost, not only because large elements relative to the crystal cell size can be used, but also because of its reduced dimensionality. The first configuration studied consists of simply rolling a graphene sheet into a tube in an isometric transformation, without any kind of relaxation. This configuration is called Original tube in Table 1. The table shows the excellent agreement between the energy obtained with the molecular model and that obtained via the continuum model and FE. According to the last remark of Section 6.2, the energy of the continuum model should be exact in this situation. Note however that the continuum membrane is discretized using an approximation space that does not reproduce exactly a circle, and thus introduces discretization errors. Then several ‘‘plane strain’’ situations are considered. This condition is enforced in the molecular model by prescribing to zero the nuclear displacements in the direction of the axis of the tube at the nuclei located at both ends of the tube. In the continuum model, we simply enforce K1¼1. Also, two kinds of energy minimization are considered. The first one freezes the inner displacements to those of the graphene sheet in equilibrium, i.e. in the continuum model by prescribing g¼0. In this particular example, invoking symmetry considerations, this constrained minimization can be easily implemented in the molecFig. 9. (a) Actual molecular model used in comparison, (b) comparison of 20 element exponential Born rule continuum model with MM and (c) results obtained with a model constructed from the standard Born rule. Table 1 Comparison of 20 element model (CþFE) with MM: energy in J/mol g0 Relaxed g MM CþFE Error (%) MM CþFE Error (%) Original tube 14.58 14.69 0.81 Relaxed tube 10.26 10.28 0.22 6.338 6.324 0.22 Squeezed A 20.85 21.22 1.8 12.95 13.11 1.2 Squeezed B 48.56 49.17 1.3 30.68 30.46 0.75 18 ular model by prescribing to zero the displacements of all the nuclei in the direction of the tube axis. This incomplete analysis is performed to highlight the effect of the inner relaxation. The other analysis is an unconstrained structural optimization of all the nuclear positions. In the continuum, the inner displacements are relaxed in order to calculate b WW at each Gauss point. The situations considered are: Relaxed tube: The Original tube is relaxed without any constraint other than the plane strain conditions. Squeezed A: Displacements at the ends of one diameter of the tube are prescribed so that this diameter of is squeezed to 3=4 of its original size. Squeezed B: Displacements at the ends of one diameter of the tube are prescribed so that this diameter of is squeezed to 1=2 of its original size. Table 1 presents the equilibrium energies for both the MM and the Continuum FE simulations, as well as the relative error of the FE calculation with respect to MM. A positive value of error means that the MM energy is lower than the FE energy. Note that this error includes contributions not only from the modelling of the discrete atomic system as a membrane, but also from the FE discretization. Fig. 9(b) compares the equilibrium configurations for the continuum/FE model and the MM model in the Squeezed B situation. Despite the large deformations to which the tube is subjected, the agreement is excellent. Table 1 shows that the equilibrium energies obtained with the continuum model are in all the cases very accurate approximations of the MM energies. The discrepancies are in all the cases below 2%. The effect of the inner relaxation in the magnitude of the energies is very important. In this table 20 finite elements have been used, while 32 hexagonal cells span the same perimeter in the MM model. Therefore, we expect the FE model to be more constrained and therefore yield higher equilibrium energies. This can be noticed in the columns corresponding to frozen inner displacements. However, when those are relaxed, the FE model reaches lower energies than the molecular model, still remaining very accurate. Probably the continuum treatment of the inner displacements allows for this extra relaxation. Although the effect of the inner relaxation in the equilibrium energies is very important, in these simulations its effect on the stable configurations is negligible. This can be explained by noting that the in-plane behavior of the model is very stiff, while the flexural behavior is very compliant. Therefore, a slight perturbation of in-plane deformation (the inner rearrangements are an in-plane effect) has dramatic influence on energy, but not in these flexural-dominated optimal deformations. This suggests that in these examples, the inner relaxation is nearly uncoupled from the bending deformation. This is not the case for other types of deformation (Arroyo and Belytschko, 2002). Table 2 shows the results obtained with 36 finite elements. The errors obtained are smaller in all the cases except in the Squeezed B situation with inner relaxation. This indicates that in general the richer discretization decreases the overall error, but also that the finer mesh allows other modeling errors to manifest themselves. Indeed, the error probably increases in the last case because the continuum model is more compliant than the molecular one with regard to the inner displacements. However, simulations carried out with even finer meshes indicate that the results ‘‘converge’’ to a very accurate result. Thus, even if the FE model is refined beyond the unit cell size, the continuum model Table 2 Comparison of 36 element model (CþFE) with MM: energy in J/mol g0 Relaxed g MM CþFE Error (%) MM CþFE Error (%) Original tube 14.58 14.60 0.14 Relaxed tube 10.26 10.28 0.21 6.338 6.324 0.22 Squeezed A 20.85 21.05 0.96 12.95 12.99 0.31 Squeezed B 48.56 49.01 0.93 30.68 30.34 1.1 19 apparently does not exhibit fine features that cannot be present in the molecular model. This excellent behavior contrasts with the situation encountered when a continuum model for the membrane is directly constructed from the Born rule without the proposed exponential extension. In this case the resulting hyper-elastic potential is non-convex. Indeed, as discussed in Section 4.3, the energy of such a model is invariant under isometric deformations (bending without stretch), i.e. the model has zero bending stiffness. This reflects in a pathological mesh dependency in the numerical implementation of such a model: since the discrete FE space cannot represent all isometric deformations, the discrete problem can still be solved, but as the mesh is refined, the numerical method picks solutions with increasingly finer features. Fig. 9(c) illustrates this fact, and sharper kinks in the numerical solution are observed as the mesh is refined. The equilibrium energy of the FE solutions is almost zero, which is not realistic. This is reminiscent of the situation encountered in other materials, for which the Fig. 13. Equilibrium configuration of a bundle of seven closely packed [22,0] nanotubes. Fig. 10. Which is more stable, circular or collapsed? (Answer: for the [20,0] and [26,0] tubes, circular, and for the [32,0] and [40,0] tubes, collapsed.) Fig. 12. Equilibrium configurations for pairs of nanotubes in van der Waals contact. Fig. 11. Transverse stability of a multi walled nanotube. 20 strain energy density is physically con-convex, leading to non-unique solutions with increasingly fine features, as reported by Dacorogna (1989, p. 276) and references therein. 7.2. Transverse deformation simulations The next simulations illustrate the application of the continuum/FE model to the transverse mechanics of nanotubes in different situations. In these applications, the computational cost of analogous MM simulations would be much higher than the cost of the presented calculations. This is especially true with regards to the non-bonded interactions. The first example studies the stability of the circular and the collapsed configurations of carbon nanotubes. Because of the van der Waals attraction potential, the energy of the system is reduced when two walls adhere. On the other hand, for the wall of a nanotube to come in contact with itself, significant elastic energy is required. This tradeoff is probably responsible for the observation by Gao et al. (1998) that below a certain radius, only the circular configuration is stable. For greater radii, the collapsed configuration is at least meta-stable. Subsequently, another threshold radius separates the nanotubes for which the circular configuration is energetically favorable from those in which the collapsed configuration is. Fig. 10 shows the simulations performed for several nanotubes. In this and subsequent figures, the nodes shown are nodes of the finite element mesh; they are not atoms. In these calculations, the fully relaxed circular configuration is deformed so that the wall of the nanotube is brought in contact with itself at the van der Waals equilibrium distance, and then the energy is minimized. The sign of the difference in energy between the circular configuration and the relaxed configuration is also reported, i.e. a positive difference means that the energy of the configuration presented on the right is lower. In some cases, the nanotube goes back to the original configuration (this is the case of the [20,0] nanotube). This implies that the collapsed configuration is not stable. The collapsed configuration is stable for the [26,0] nanotube, but this only constitutes a local minimum of the energy since the circular configuration has lower energy. For the [32,0] and [40,0] nanotubes, the collapsed configuration is the energetically favorable structure. This is expected because larger nanotubes are more flexible and have more wall area to gain adhesion energy. Fig. 11 displays a similar analysis for a multi-walled nanotube for which the collapsed configuration yields lower energy . A similar competition of elastic and adhesion energy occurs when two nanotubes are brought to the van der Waals equilibrium distance. Fig. 12 shows the equilibrium configurations obtained when this numerical experiment is performed with nanotubes of different sizes. Again, the larger nanotubes have larger portions of flattened walls. We also report a simulation of a bundle of nanotubes under plane strain. Fig. 13 shows the equilibrium configuration of the system. A TEM image of such a nanorope has been reported by Salvetat et al. (1999). Carbon nanotubes tend to be closely packed in hexagonal lattices in the nanoropes and crystals of nanotubes (Thess et al., 1996; Schlittler et al., 2001). As can be seen from Fig. 13, the equilibrium configuration displays a flattening of the nanotube walls, or partial polygonalization. 7.3. Three dimensional simulation The theory presented has been used to construct a membrane applicable in the general 3D Fig. 14. Twisting of a [10,10] nanotube: deformed geometry for twisting angles of 38,210and 360, and cross section of the deformed membrane at the center of the tube for the above three configurations. 21 deformation of carbon nanotubes (Arroyo and Belytschko, 2002). This more general membrane can be discretized with subdivision finite elements, and the structural instabilities reported in experiments and atomistic simulations can be analyzed at very low computational cost. The analysis of twisting a [10,10] nanotube is provided in Fig. 14, for a Tersoff-Brenner potential. Note that, the deformed geometries have been post-processed, and the computational mesh has about 18 elements around the perimeter. Each end of the nanotube is incrementally rotated 360in opposite orientations. The first snapshot of the deformation shows the configuration when the first instability from a uniform twisting occurs, and the corresponding cross-section is shown at the bottom of the figure. Further twisting causes the wall of the nanotube to come in van der Waals contact with itself, as clearly shown in the cross-section in the bottom of Fig. 14. Beyond 210, a secondary instability develops, and the tube folds onto itself. From the cross-section it is apparent that the van der Waals interactions are responsible for this buckled morphology. In the absence of these long-range forces, the membrane inter-penetrates and the secondary structure is not observed. This 3D membrane has been shown to provide very accurate energetics and deformed geometries even for very large deformations (Arroyo and Belytschko, 2002). 8. Conclusions We have further explored a methodology to construct continuum models for one-atom thick crystalline films. The proposed model is a hyperelastic membrane whose elastic potential energy is written in closed-form exclusively in terms of the inter-atomic potentials that constitute the molecular description of the system. The analysis of the present work is based on the exponential the Born rule (Arroyo and Belytschko, 2002), a kinematic assumption linking the atomic and the continuum deformations when the crystal is a curved film. This extension is based on the exponential map. An illustrative example of an atomic chain deforming in two dimensions has been presented. The resulting simple rope-like continuum model encompasses all of the fundamental ideas. The general methodology then is particularized to analyze the transverse mechanics of carbon nanotubes. This model explicitly exploits the symmetry of such a deformation, and leads to a model of reduced dimensionality. The hyper-elastic potential, as well as strain and stress measures are provided, and a continuum formulation of the non-bonded interactions is derived. The proposed model is discretized using finite elements, yielding an alternative simulation method that is faster than atomistic calculations. Several simulations highlighting the relevance of van der Waals interactions in the transverse mechanics of nanotubes are reported. The results show that the continuum model based on the exponential Born rule very well approximates the stable configurations and energies of the corresponding MM model. Results agree with MM calculations within 2% in the equilibrium energies. This sharply contrasts with the non-physical results obtained from a model based on the standard Born rule. We also show the important effect of the inner rearrangements of the crystal structure on the equilibrium energies. A full 3D simulation illustrates the application of the present theory to analyze the structural instabilities of nanotubes observed in experiments and atomistic calculations. Acknowledgements The support of the ‘‘la Caixa’’ Graduate Program to M. Arroyo, and the National Science Foundation and the U.S. Army Research Office is gratefully acknowledged. References Arroyo, M., Belytschko, T., 2002. An atomistic based finite deformation membrane for single layer crystalline films. Journal of the Mechanics and Physics of Solids 50, 1941 1977. Bernholc, J., Brabec, C.J., Nardelli, M.B., Maiti, A., Roland, C., Yakobson, B.I., 1998. Theory of growth and mechanical properties of nanotubes. Applied Physics A 67, 39 46. 22 Brenner, D.W., 1990. Empirical potential for hydrocarbons for use in simulating chemical vapor deposition of diamond films. Physical Review B 42 (15), 9458 9471. 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. Cousins, C.S.G., 1978. Inner elasticity. Journal of Physics C, 4867 4879. Dacorogna, B., 1989. Direct methods in the calculus of variations. In: Applied Mathematical Sciences, vol. 78. Springer Verlag, Berlin. do Carmo, M.P., 1976. Differential geometry of curves and surfaces. Prentice Hall, Englewood Cliffs, NJ. Ericksen, J.L., 1984. Phase transformations and material instabilities in solids. In: Gurtin, M.E. (Ed.), The Cauchy and Born Hypotheses for Crystals. Academic Press, Lon don, pp. 61 77. Falvo, M.R., Clary, G.J., Taylor, R.M., Chi, V., Brooks, F.P., Washburn, S., Superfine, R., 1997. Bending and buckling of carbon nanotubes under large strain. Nature 389, 582 584. Friesecke, G., James, R.D., 2000. A scheme for the passage from atomic to continuum theory for thin films, nanotubes and nanorods. Journal of the Mechanics and Physics of Solids 48, 1519 1540. Gao, G., C ßa ggin, T., Goddard III, W., 1998. Energetics, structure, mechanical and vibrational properties of single walled carbon nanotubes. Nanotechnology 9, 184 191. Lu, J.P., 1997. Elastic properties of carbon nanotubes and nanoropes. Physical Review Letters 79 (7), 1297 1300. Maiti, A., 2000. Mechanical deformation in carbon nanotubes bent tubes vs tubes pushed by atomically sharp tips. Chemical Physical Letters 331, 21 25. Malvern, L.E., 1969. Introduction to the mechanics of a continuous medium. Prentice Hall, Englewood Cliffs, NJ. Marsden, J.E., Hughes, T.J., 1983. Mathematical foundations of elasticity. Prentice Hall, Englewood Cliffs, NJ. Martin, J.W., 1975. Many body forces in metals and the brug ger elastic constants. Journal of Physics C 8, 2837 2857. Morgan, F., 1993. Riemannian geometry, a beginnerÕs guide. Jones and Bartlett Publishers. Nevins, N., Chen, K., Allinger, N.L., 1996. Moleculas mechan ics (MM4) calculations on alkenes. Journal of Computa tional Chemistry 17 (5 6), 669 694. Qian, D., Liu, W.K., Ruoff, R.S., 2001. Mechanics of C60 in nanotubes. Journal of Physical Chemistry B 105, 10753 10758. Ruoff, R.S., Tersoff, J., Lorents, D.C., Subramoney, S., Chan, B., 1993. Radial deformation of carbon nanotubes by van der Waals forces. Nature 364, 514 516. Saito, R., Fujita, M., Dresselhaus, G., Dresselhaus, M.S., 1992. Electronic structure of chiral graphene tubules. Applied Physics Letters 60 (18), 2204 2206. Salvetat, J.P., Briggs, G.A.D., Bonard, J.M., Bacsa, R.R., Kulik, A.J., Stoockli, T., Burnham, N.A., Forr oo, L., 1999. Elastic and shear moduli of single walled carbon nanotube ropes. Physical Review Letters 82 (5), 944 947. Schlittler, R.R., Seo, J.W., Gimzewski, J.K., Durkan, C., Saifullah, M.S.M., Welland, M.E., 2001. Single crystals of single walled carbon nanotubes formed by self assembly. Science 292, 1136 1139. Shenoy, V.B., Miller, R., Tadmor, E.B., Rodney, D., Phillips, R., Ortiz, M., 1999. An adaptive finite element approach to atomic scale mechanics the quasicontinuum method. Journal of the Mechanics and Physics of Solids 47, 611 642. Sohlberg, K., Sumpter, B.G., Tuzun, R.E., Noid, D.W., 1998. Continuum methods of mechanics as a simplified approach to structural engineering of nanostructures. Nanotechnol ogy 9, 30 36. Tadmor, E.B., Ortiz, M., Phillips, R., 1996. Quasicontinuum analysis of defects in solids. Philosophical Magazine A 73 (6), 1529 1563. Tadmor, E.B., Smith, G.S., Bernstein, N., Kaxiras, E., 1999. Mixed finite element and atomistic formulation for complex crystals. Physical Review B 59 (1), 235 245. Thess, A., Lee, R., Nikolaev, P., Dia, H., Petit, P., Robert, J., Xu, C., Lee, Y.H., Kim, S.G., Rinzler, A.G., Colbert, D.T., Scuseria, G.E., Tom aanek, D., Fischer, J.E., Smalley, R.E., 1996. Crystalline ropes of metallic carbon nanotubes. Science 273, 483 487. Weiner, J.H., 1983. Statistical mechanics of elasticity. Wiley, New York. Yakobson, B.I., Brabec, C.J., Bernholc, J., 1996. Nanome chanics of carbon tubes: Instabilities beyond the linear response. Physical Review Letters 76 (14), 2511 2514. Yu, M., Dyer, M.J., Ruoff, R.S., 2001a. Structure and mechanical flexibility of carbon nanotube ribbons: An atomic force microscopy study. Journal of Applied Physics 89 (8), 4554 4557. Yu, M., Kowalewski, T., Ruoff, R.S., 2001b. Structural analysis of collapsed, and twisted and collapsed, multiwalled carbon nanotubes by atomic force microscopy. Physical Review Letters 86 (1), 87 90. Yu, M., Lourie, O., Dyer, M.J., Moloni, K., Kelly, T.F., Ruoff, R.S., 2000. Strength and breaking mechanism of multi walled carbon nanotubes under tensile load. Science 287, 637 640. Zanzotto, G., 1996. The Cauchy Born hypothesis, nonlinear elasticity and mechanical twinning in crystals. Acta Crys tallographica A 52, 839 849. Zhong Can, O. Y., Su, Z. B., Wang, C. L., 1997. Coil forma tion in multishell carbon nanotubes: competition between curvature elasticity and interlayer adhesion. Physical Re view Letters 78 (21), 4055 4058. Zhou, G., Duan, W., Gu, B., 2001. First principles study on morphology and mechanical properties of single walled carbon nanotube. Chemical Physical Letters 333, 344 349. 23