A new shell formulation for graphene structures based on
existing ab-initio data
Reza Ghaffari1 , Thang X. Duong, Roger A. Sauer2
Aachen Institute for Advanced Study in Computational Engineering Science (AICES),
RWTH Aachen University, Templergraben 55, 52056 Aachen, Germany

arXiv:1612.08965v2 [physics.comp-ph] 1 Dec 2017

Published3 in International Journal of Solids and Structures, DOI: 10.1016/j.ijsolstr.2017.11.008
Submitted on 28. December 2016, Revised on 1. November 2017, Accepted on 10. November 2017

Abstract
An existing hyperelastic membrane model for graphene calibrated from ab-initio data (Kumar
and Parks, 2014) is adapted to curvilinear coordinates and extended to a rotation-free shell
formulation based on isogeometric finite elements. Therefore, the membrane model is extended
by a hyperelastic bending model that reflects the ab-inito data of Kudin et al. (2001). The
proposed formulation can be implemented straight-forwardly into an existing finite element
package, since it does not require the description of molecular interactions. It thus circumvents
the use of interatomic potentials that tend to be less accurate than ab-initio data. The proposed
shell formulation is verified and analyzed by a set of simple test cases. The results are in
agreement to analytical solutions and satisfy the FE patch test. The performance of the shell
formulation for graphene structures is illustrated by several numerical examples. The considered
examples are indentation and peeling of graphene and torsion, bending and axial stretch of
carbon nanotubes. Adhesive substrates are modeled by the Lennard-Jones potential and a
coarse grained contact model. In principle, the proposed formulation can be extended to other
2D materials.
Keywords: adhesive contact; carbon nanotubes; graphene; hyperelasticity; isogeometric finite
elements; rotation-free shell.

1

Introduction

Due to their extraordinary mechanical, thermal and electrical properties, carbon based structures (like graphene, carbon nanotubes (CNT), fullerene) are of high interest in many industrial
applications (Balandin, 2011; Marinho et al., 2012). Graphene can be considered as the basic
structure to build other carbon based structures. So, the first step to analyze and obtain material properties of these structures is the investigation of graphene.
Several methods are used to model and simulate graphene structures. Among them, first principle simulations are the most exact numerical methods. They simulate interaction of electrons,
but they are limited to small time and length scales. The tight-binding (TB) method and density functional theory (DFT) can solve larger problems. However, they are still restricted to a
few hundred atoms. The high computational cost of the mentioned methods motivates the usage
of empirical or semi-empirical potentials in the framework of molecular dynamics (MD) simulations. These potentials concentrate on atomic interactions and ignore electron interactions
1

email: ghaffari@aices.rwth-aachen.de
corresponding author, email: sauer@aices.rwth-aachen.de
3
This pdf is the personal version of an article whose final publication is available at http://sciencedirect.com
2

1

(Ansari et al., 2012). When the size of a problem reaches a few hundreds of nano-meters, MD
simulations become impractical and other methods should be used. Continuum based methods
are another option, however classical continuum mechanics does not capture size and boundary
effects at the nano scale (Cauchy, 1851). These issues can be resolved by the surface CauchyBorn rule (Park et al., 2006), boundary Cauchy-Born rule (Qomi et al., 2011) and modified
boundary Cauchy-Born rule (Ghaffari et al., 2015).
The time and size restriction of quantum mechanics (QM) and MD methods and the inaccuracy
of classical continuum mechanical theories motivate the use of multiscale methods. The most
well-known multiscale method to simulate graphene and CNTs is the exponential Cauchy-Born
(ECB) rule developed by Arroyo and Belytschko (2002, 2004a,b, 2005). These works employ
the first generation Brenner potential (FGBP) (Brenner, 1990) and second generation Brenner
potential (SGBP) (Brenner et al., 2002). These potentials and the multiscale methods based on
them, underestimate the elastic modulus and bending stiffness compared to QM simulations.
Hence, the usage of QM becomes necessary (Arroyo and Belytschko, 2004a). The elastic properties obtained by different methods are compared by Cao (2014). QM data can be either used to
enrich the atomistic potential (Lindsay and Broido, 2010), or calibrate the strain energy density
(Xu et al., 2012).
Classical continuum theories have been extended to lattice structures by Ericksen (1979), Pitteri
(1985) and Fadda and Zanzotto (2000). Based on this, Sfyris et al. (2014b) use Taylor expansion of the strain energy in order to describe graphene. Sfyris et al. (2014a) obtain a set of
invariants based on the right Cauchy-Green strain tensor, curvature tensor and shift vector and
propose a strain energy functional based on a linear combination of those invariants. Delfani
et al. (2013) and Delfani and Shodja (2015, 2016) use Taylor expansion for the strain energy
and apply the symmetry operators to the elasticity tensors in order to reduce the number of
independent variables. These works consider a finite thickness for graphene and use through
the thickness integration of the strain energy in order to obtain the bending stiffness. However,
the thickness of graphene is a controversial quantity. On the other hand, the bending modulus
can be obtained directly from the change of the bond and dihedral angles (Lu et al., 2009)
without introducing a thickness. Kumar and Parks (2014) proposed a membrane surface strain
energy per unit area that does not account for bending resistance. This is included here by
using the Canham bending strain energy (Canham, 1970) in conjunction with the nonlinear
shell formulation of Duong et al. (2017), which does not require the notion of a thickness.
Experimental tests can be considered as an alternative method to capture the properties of
graphene. Such tests are very important to calibrate and validate numerical simulations. The
most common experimental tests are uniaxial stretching, pure shear, simple shear, torsion and
indentation. The uniaxial and biaxial tests require grips to hold a specimen and it is difficult to
build nanoscale grips. The experimental indentation test setup needs no grips and it can be built
by laminating a graphene sheet on a substrate, e.g. SiO2 , via exfoliation. Force-displacement
curves can then be obtained by using an atomic-force microscope (AFM). It should be noted
that initial stress within the graphene and adhesion to the substrate play an important role in
experimental results (Lee et al., 2008).
The adhesion energy between a substrate and graphene can be obtained from blister tests
(Koenig et al., 2011; Boddeti et al., 2013). The atomic structure of graphene on a substrate can
be investigated experimentally (Ishigami et al., 2007) and theoretically (Hossain, 2009; Bellido
and Seminario, 2010; Lu and Dunn, 2010; He et al., 2013; Gao et al., 2014; Sfyris and Sfyris,
2017). The corrugation of substrates is investigated by Li and Zhang (2010) and Aitken and
Huang (2010). Adhesion effects on micro-mechanical structures are investigated by Maboudian
and Howe (1997). Recent development on the adhesion of graphene membranes is reviewed by
Bunch and Dunn (2012). A continuum model is developed for multi-buckling of graphene on a
substrate by Gao et al. (2016). The effects of substrate morphology on graphene is investigated

2

by a continuum model and verified by an atomistic simulation (Boddeti et al., 2016). The effects
of substrate rippling on the bending stiffness of a graphene are considered by Jomehzadeh and
Pugno (2015).
Results of the indentation test can be used to obtain graphene properties. However, connecting
in-plane material properties with the force-displacement curve of the indentation test is a challenging issue. The connection between the indentation results and in-plane material properties
is investigated by Han et al. (2015). Moreover, this test does not yield any direct information
about the constitutive model, like the functional form of the strain energy density. It has been
verified that graphene behaves isotropically in infinitesimal strain Cao (2014). But under large
deformations the material response depends on the relative direction of loading and the chirality. Therefore, graphene can not be considered isotropic in general (Larsson and Samadikhah,
2011; Cao, 2014). The anisotropic response of graphene can be simulated with the structural
tensor obtained from the symmetry groups of the lattice structure (Kumar and Parks, 2014).
In this approach, a set of virtual DFT experiments are conducted and the coefficients of the
strain energy are calibrated. Hence, there is no need to study the deformation of the lattice as
it is done in the Cauchy-Born rule.
In the research group of the authors, a novel curvilinear finite element (FE) formulation has
been developed and applied to computations of liquid and solid membranes (Sauer et al., 2014).
It has been extended to rotation-free shells (Sauer and Duong, 2017; Duong et al., 2017) and
employed in the modeling of anisotropic materials (Roohbakhshan et al., 2016). The new formulation is based on isogeometric FE (Cottrell et al., 2009) that provides higher accuracy than
Lagrange-based FE. It is extended here to model and simulate graphene based structures.
The highlights and novelties of the proposed new shell model are
∙ It combines the existing ab-initio data of Kudin et al. (2001) and Kumar and Parks (2014)
with the rotation-free isogeometric shell model of Duong et al. (2017).
∙ It can be implemented straightforwardly within the finite element method, which is more
efficient than molecular dynamics for large systems.
∙ It is closer to experimental and quantum results than shell models based on interatomic
potentials, since those often underestimate the elastic modulus of graphene by about one
third.
∙ It is directly based on a surface strain energy defined in closed form that avoids the
evaluation of molecular interactions as they appear, e.g. in the Cauchy-Born rule.
∙ It accounts for contact and adhesion, allowing to study indentation and peeling of graphene.
∙ It is fully nonlinear and can be used to accurately analyze the buckling behavior of
graphene structures.
∙ It can be extended to other 2D material models, like those of Sfyris et al. (2015).
The remainder of this paper is organized as follows: In Sec. 2, the kinematics of deforming
surfaces are presented and based on that, in Secs. 3.1 and 3.2, the membrane and bending
constitutive laws are developed. Details of the considered FE formulation are provided in Sec. 4.
In Sec. 5, the proposed constitutive law and FE formulation are verified by several elementary
benchmark tests. Further, numerical examples are presented in Sec. 6 to demonstrate the
capability of the proposed model for the simulation of graphene based specimens. The paper is
concluded in Sec. 7.

3

2

Kinematics

In this section, the kinematics of deforming surfaces are described. This description is then
used in the constitutive modeling and FE formulation.

2.1

Kinematics of deforming surfaces

The surface description is essential for membrane and shell formulations. Sauer and Duong
(2017) proposed a rotation-free curvilinear framework for membranes and shells. The parametric
description of the surface in the reference and the current configuration can be written as
𝑋 = 𝑋(𝜉 1 , 𝜉 2 ) ,

(1)

𝑥 = 𝑥(𝜉 1 , 𝜉 2 ) ,

(2)

and the corresponding tangent vectors are
𝐴𝛼 =

𝜕𝑋
,
𝜕𝜉 𝛼

(3)

𝑎𝛼 =

𝜕𝑥
.
𝜕𝜉 𝛼

(4)

Next, the co-variant components of the metric tensors are defined from the inner product as
𝐴𝛼𝛽 = 𝐴𝛼 · 𝐴𝛽 ,

(5)

𝑎𝛼𝛽 = 𝑎𝛼 · 𝑎𝛽 ,

(6)

and the contra-variant components of the metric tensors are defined as the inverse of the covariant metric tensors, i.e.
[𝐴𝛼𝛽 ] = [𝐴𝛼𝛽 ]−1 ,
(7)
[𝑎𝛼𝛽 ] = [𝑎𝛼𝛽 ]−1 .

(8)

The dual tangent vectors can then be defined as
𝐴𝛼 := 𝐴𝛼𝛽 𝐴𝛽 ,

(9)

𝑎𝛼 := 𝑎𝛼𝛽 𝑎𝛽 .

(10)

The normal unit vector of the surface in its reference and current configuration can then be
written as
𝐴1 × 𝐴2
𝑁=
,
(11)
‖𝐴1 × 𝐴2 ‖
𝑛=

𝑎1 × 𝑎2
.
‖𝑎1 × 𝑎2 ‖

(12)

The 3D identity tensor 1 can be written in terms of the surface identity tensors in reference
configuration, 𝐼 , and current configuration, 𝑖 , as
1=𝐼 +𝑁 ⊗𝑁 =𝑖+𝑛⊗𝑛 ,

(13)

where 𝐼 and 𝑖 can be written as
𝐼 = 𝐴𝛼 ⊗ 𝐴𝛼 = 𝐴𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 = 𝐴𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ,
4

(14)

𝑖 = 𝑎𝛼 ⊗ 𝑎𝛼 = 𝑎𝛼𝛽 𝑎𝛼 ⊗ 𝑎𝛽 = 𝑎𝛼𝛽 𝑎𝛼 ⊗ 𝑎𝛽 .

(15)

Next, the curvature tensor can be written as
𝑏 = 𝑏𝛼𝛽 𝑎𝛼 ⊗ 𝑎𝛽 ,

(16)

where 𝑏𝛼𝛽 are the co-variant components of the curvature tensor defined as
𝑏𝛼𝛽 := 𝑛 · 𝑎𝛼,𝛽 = 𝑛 · 𝑎𝛼;𝛽 ,

(17)

and 𝑎𝛼,𝛽 and 𝑎𝛼;𝛽 are the parametric and co-variant derivatives of the tangent vectors. They
are connected by
𝑎𝛼;𝛽 = 𝑎𝛼,𝛽 − Γ𝛾𝛼𝛽 𝑎𝛾 ,
(18)
where Γ𝛾𝛼𝛽 is the Christoffel symbol of the second kind, which is defined as
Γ𝛾𝛼𝛽 := 𝑎𝛼,𝛽 · 𝑎𝛾 .

(19)

The contra-variant components of 𝑏 are defined as
𝑏𝛼𝛽 := 𝑎𝛼𝛾 𝑏𝛾𝛿 𝑎𝛿𝛽 .

(20)

The frame indifference of the material model is an important requirement in continuum mechanics. It is satisfied if models are written in term of tensor invariants. The mean and Gaussian
curvatures are a good set of invariants to characterize the bending energy. They are defined as
1
1
𝐻 := 𝑏 : 𝑖 = 𝑏𝛼𝛼 ,
2
2
𝜅 := det(𝑏) =

det[𝑏𝛼𝛽 ]
.
det[𝑎𝛾𝛿 ]

(21)
(22)

The definition of the determinant of surface tensors can be found in Javili et al. (2014). Alternatively, the invariants can be defined based on the principal curvatures 𝜅𝛼 as
𝐻 :=

1
(𝜅1 + 𝜅2 ) ,
2

𝜅 := 𝜅1 𝜅2 .

(23)
(24)

The principal curvatures are the eigenvalues of matrix [𝑏𝛼𝛽 𝑎𝛽𝛾 ].

2.2

Kinematics of deformation

The first step to develop a material model is the introduction of stress and strain measures. The
logarithmic strain 𝐸 (0) is sensitive to material nonlinearities even for small strains. Furthermore,
it can be additively decomposed into volumetric and shear parts. These features make it an
excellent candidate for the development of material models (Neff et al., 2013; Kumar and Parks,
2014; Neff et al., 2015a,b; Ghiba et al., 2015; Montella et al., 2016; Neff and Ghiba, 2016). The
Lagrangian logarithmic strain is defined as the logarithm of the right stretch tensor 𝑈 . The
standard method to calculate non-integer powers and the logarithm of a tensor is based on the
spectral decomposition. The surface deformation gradient can be written as
𝐹 = 𝑅𝑈 ,
5

(25)

where 𝑅 is the rotation tensor. Using the spectral decomposition, 𝑈 can be written based on
its eigenvalues 𝜆𝑖 and eigenvectors 𝑌 𝑖 , which are the principal stretches and their directions,
respectively. So, the spectral decomposition of 𝑈 is
∑︁
𝑈=
𝜆𝑖 𝑌 𝑖 ⊗ 𝑌 𝑖 .
(26)
𝑖=1,2

Using the same decomposition, the right Cauchy-Green tensor 𝐶 can be written based on its
eigenvalues Λ𝑖 and eigenvectors 𝑌 𝑖 . It should be noted that 𝑈 and 𝐶 have the same principal
directions, i.e. eigenvectors, and Λ𝑖 = 𝜆2𝑖 . So, 𝐶 can be written as
∑︁
𝐶=
Λ𝑖 𝑌 𝑖 ⊗ 𝑌 𝑖 .
(27)
𝑖=1,2

The eigenprojection tensors 𝑌 𝑖 ⊗ 𝑌 𝑖 can be analytically obtained from Sylvester’s formula
(Itskov, 2015)
𝐶 − Λ2 𝐼
1
P1 := 𝑌 1 ⊗ 𝑌 1 =
=: 𝑃𝛼𝛽
𝐴𝛼 ⊗ 𝐴𝛽 ,
Λ1 − Λ2
(28)
𝐶 − Λ1 𝐼
2
𝛼
𝛽
=: 𝑃𝛼𝛽 𝐴 ⊗ 𝐴 ,
P2 := 𝑌 2 ⊗ 𝑌 2 =
Λ2 − Λ1
𝑖 are defined as
where 𝑃𝛼𝛽

1
(𝑎𝛼𝛽 − Λ2 𝐴𝛼𝛽 ) ,
Λ1 − Λ2
1
:=
(𝑎𝛼𝛽 − Λ1 𝐴𝛼𝛽 ) .
Λ2 − Λ1

1 :=
𝑃𝛼𝛽
2
𝑃𝛼𝛽

(29)

Using the spectral decomposition of 𝑈 , the logarithmic strain can be written as
𝐸 (0) = ln(𝑈 ) =

1
ln(𝐽) 𝐼 + ln 𝜆 (𝑌 1 ⊗ 𝑌 1 − 𝑌 2 ⊗ 𝑌 2 ) ,
2

(30)

where the two terms on the right hand side are the aeral (surface dilatation) and deviatoric
parts of the logarithmic strain, and 𝜆 and 𝐽 are defined as
√︂
𝜆1
(31)
𝜆 :=
,
𝜆2
𝐽 := 𝜆1 𝜆2 .

(32)

They are tensor invariants, just like 𝐻 and 𝜅.
Based on the presented kinematics, a hyperelastic strain energy function will be introduced for
curvilinear coordinates in the next section.

3

Material model

In this section, the membrane constitutive law of Kumar and Parks (2014) is extended to the
curvilinear shell theory of Sauer and Duong (2017) and its computational counterpart (Duong
et al., 2017).
For hyperelastic materials, it is assumed that all energy is stored as elastic energy and the
dissipation is zero. The strain energy density is defined directly on the surface (as energy per
reference area). It should not be confused with strain energy per volume, which is not needed
here. For Kirchhoff-Love kinematics, the strain energy can be written as a function of the metric
6

and curvature tensor. Here, the strain energy is defined analytically in closed form such that
the evaluation of molecular interactions is avoided. In the framework of the shell model of Sauer
and Duong (2017), the derivative of the strain energy density with respect to 𝑎𝛼𝛽 and 𝑏𝛼𝛽 gives
the stress and the bending moment as
𝜏 𝛼𝛽 = 2

𝜕𝑊
,
𝜕𝑎𝛼𝛽

(33)

𝑀0𝛼𝛽 =

𝜕𝑊
.
𝜕𝑏𝛼𝛽

(34)

In addition, the elasticity tensors are defined as
𝑐𝛼𝛽𝛾𝛿 := 4

𝜕2𝑊
,
𝜕𝑎𝛼𝛽 𝜕𝑎𝛾𝛿

(35)

𝑑𝛼𝛽𝛾𝛿 := 2

𝜕2𝑊
,
𝜕𝑎𝛼𝛽 𝜕𝑏𝛾𝛿

(36)

𝑒𝛼𝛽𝛾𝛿 := 2

𝜕2𝑊
,
𝜕𝑏𝛼𝛽 𝜕𝑎𝛾𝛿

(37)

𝜕2𝑊
,
𝜕𝑏𝛼𝛽 𝜕𝑏𝛾𝛿

(38)

𝑓 𝛼𝛽𝛾𝛿 :=

where 𝑊 can be written as summation of membrane and bending parts,
𝑊 = 𝑊m + 𝑊b .

(39)

Given 𝜏 𝛼𝛽 and 𝑀0𝛼𝛽 , the Cauchy stress tensor can be written as
𝜎 = 𝑁 𝛼𝛽 𝑎𝛼 ⊗ 𝑎𝛽 + 𝑆 𝛼 𝑎𝛼 ⊗ 𝑛 ,

(40)

𝑁 𝛼𝛽 = 𝜎 𝛼𝛽 + 𝑏𝛼𝛾 𝑀 𝛾𝛽 ,

(41)

𝛽𝛼
,
𝑆 𝛼 = −𝑀;𝛽

(42)

where

and 𝜎 𝛼𝛽 = 𝜏 𝛼𝛽 /𝐽 and 𝑀 𝛼𝛽 = 𝑀0𝛼𝛽 /𝐽 . The mixed in-plane components of 𝜎 are defined as
𝑁𝛽𝛼 := 𝑁 𝛼𝛾 𝑎𝛾𝛽 .

3.1

(43)

Membrane energy

Graphene is a thin 2D structure with a thickness of only one atom. This motivates the development of material models based on membrane and shell theories. A suitable strain energy
function should be selected to capture the anisotropic behavior of graphene. Based on the
isotropization theorem, an anisotropic tensorial functional can be written as an isotropic tensorial functional by utilizing structural tensors (Zheng, 1994). So, the strain energy density can
be isotropized by using the invariants obtained based on isotropization. Based on the symmetry group of graphene and using the logarithmic strain, the following three invariants can be
introduced (Kumar and Parks, 2014)
𝒥1 := 𝜖a = ln 𝐽 ;

𝐽 := 𝜆1 𝜆2 ,

1 2
1 (0)
(0)
𝛾𝑖 := 𝐸 dev : 𝐸 dev = (ln 𝜆)2 ; 𝜆1 ≥ 𝜆2 ,
4
2
[︂
)︁
(︁
)︁ (︁
)︁ ]︂
1 3
1 (︁ ^
(0) 3
(0)
(0) 2
^
^
𝒥3 := 𝛾𝜃 :=
𝑀 : 𝐸 dev − 3 𝑀 : 𝐸 dev 𝑁 : 𝐸 dev
= (ln 𝜆)3 cos(6𝜃) ,
8
8

𝒥2 :=

7

(44)

y

x

Figure 1: Orthonormal vectors characterizing the undeformed graphene lattice.
(0)

where 𝐸 dev is the deviatoric part of the logarithmic strain (Appendix A). A similar approach is
used by Sfyris et al. (2014a) to obtain another set of invariants. Here the first and second invariants model the isotropic features of the material, and the third invariant captures anisotropic
^ and 𝑁
^ are defined as
behavior of the material. 𝑀
^ := 𝑥
^ ⊗𝑥
^ −𝑦
^⊗𝑦
^,
𝑀

(45)

^ := 𝑥
^ ⊗𝑦
^+𝑦
^⊗𝑥
^ ,
𝑁

(46)

^ and 𝑦
^ are two orthonormal vectors. The former is aligned along the armchair direction
where 𝑥
(Fig. 1). 𝜃 is the maximum stretch angle relative to the armchair direction and defined as
^) .
𝜃 := arccos (𝑌 1 · 𝑥

(47)

Then, the membrane strain energy density, per unit area of the initial configuration, can be
dil and deviatoric 𝑊 dev parts as
decomposed into pure dilatation 𝑊m
m
dil (𝒥 ) + 𝑊 dev (𝒥 , 𝒥 ; 𝒥 ) .
𝑊m (𝒥1 , 𝒥2 , 𝒥3 ) = 𝑊m
1
2
3
1
m

The dilatation and deviatoric parts of the strain energy density are
[︀
]︀
dil := 𝜀 1 − (1 + 𝛼
𝑊m
^ 𝜖a ) exp(−^
𝛼 𝜖a ) ,
dev := 2 𝜇(𝜖 ) 𝒥 + 𝜂(𝜖 ) 𝒥 ,
𝑊m
a
2
a
3

(48)

(49)

where 𝜇 and 𝜂 are defined as
^

𝜇(𝜖a ) := 𝜇0 − 𝜇1 𝑒𝛽 𝜖a ,
𝜂(𝜖a ) := 𝜂0 − 𝜂1 𝜖2a .

(50)

^ 𝜂0 ,
The membrane model is thus characterized by the seven material constants 𝜀, 𝛼
^ , 𝜇0 , 𝜇1 , 𝛽,
and 𝜂1 . They are given in Tab. 1. In that table, LDA and GGA are abbreviations for local
density approximation and generalized gradient approximation. These are two approximations
used in DFT simulations in order to compute the material constants. Using the relations given
in the previous section, the contra-variant component of the Kirchhoff stress tensor for the cases
of distinct and repeated eigenvalues are
[︀ 2
]︀
1
^
𝛼𝛽
𝜏m
= 𝜀𝛼
^ 𝜖a 𝑒−𝛼^ 𝜖a − 2 (𝜂1 𝜖a 𝒥3 + 𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽 + 2 𝜇 ln 𝜆 𝜒𝛼 𝛽 + 𝜂 𝜇𝛼𝛽 , for 𝜆1 ̸= 𝜆2 ,
4
(51)
8

GGA
LDA

𝛼
^
1.53
1.38

𝜀 [N/m]
93.84
116.43

𝜇0 [N/m]
172.18
164.17

𝜇1 [N/m]
27.03
17.31

𝛽^
5.16
6.22

𝜂0 [N/m]
94.65
86.93

𝜂1 [N/m]
4393.26
3611.53

Table 1: Membrane behavior: Material constants of graphene (Kumar and Parks, 2014).

Exp. Cauchy-Born with FGBP (Arroyo and Belytschko, 2004a)
Exp. Cauchy-Born with SGBP (Arroyo and Belytschko, 2004a)
MM3 (Gupta and Batra, 2010)
Current model with GGA (Kumar and Parks, 2014)
Current model with LDA (Kumar and Parks, 2014)
Ab initio (Kudin et al., 2001)
Experimental (Lee et al., 2008)

Elastic modulus E [N/m]
236
243
340
349
354
345
340

Table 2: Membrane behavior: Comparison of the elastic modulus for various methods. FGBP
= first generation Brenner potential; SGBP = second generation Brenner potential; LDA =
local density approximation; GGA = generalized gradient approximation.
and
[︀ 2
]︀
𝛼𝛽
𝜏m
= 𝜀𝛼
^ 𝜖a 𝑒−𝛼^ 𝜖a 𝑎𝛼𝛽 , for 𝜆1 = 𝜆2 ,

(52)

where 𝜇𝛼𝛽 and 𝜒𝛼 𝛽 are defined as
∑︁

𝜇𝛼𝛽 :=

(0)

𝑓𝑖𝑗 𝑃𝑖𝛼𝛾 𝑇𝛾𝛿 𝑃𝑗𝛿𝛽 ,

(53)

𝑖,𝑗=1,2

𝜒𝛼𝛽 :=

1 𝛼𝛽
1
𝑃
− 2 𝑃2𝛼𝛽 ,
𝜆21 1
𝜆2

(54)

and 𝑓𝑖𝑗 is defined as
𝑓𝑖𝑗

1
, if 𝑖 = 𝑗 ,
2 𝜆2𝑖
:=
ln 𝜆𝑖 − ln 𝜆𝑗
⎪
⎪
, if 𝑖 ̸= 𝑗 .
⎩
𝜆2𝑖 − 𝜆2𝑗
⎧
⎪
⎪
⎨

(55)

(0)

𝑇𝛾𝛿 and the elasticity tensors related to Eqs. (51) and (52) are given in Appendix B.
dil , 𝜇 and 𝜂 with 𝒥 are presented in Fig. 2. The elastic modulus of the
The variation of 𝑊m
1
current model is compared with reference data in Tab. 2. The elastic modulus of the current
model is closer to existing experiments and ab-initio results than the exponential Cauchy-Born
rule.

3.2

Bending energy

In this section, the membrane model of Kumar and Parks (2014) is extended to a shell model
by including a bending strain energy. In the current work, the uniaxial bending model of Lu
et al. (2009) is extended to tensorial form using the tensorial invariants introduced in Sec. 2.
Lu et al. (2009) show that the bending behavior of graphene is linear up to the curvature 1
3

These material parameters were provided by S. Kumar and D. M. Parks on a personal request by the
authors as the material parameters in Kumar and Parks (2014) contain errors.

9

100
50

7 [N/m]

20

W Dil [J/m2 ]

150

GGA
LDA

15
10

100

2 [N/m]

25

GGA
LDA

50

-50
-100

5
0

0

-150
0

0.2

0.4

0.6

0

0

0.1

J1

0.2

0.3

0.4

-200

GGA
LDA

0

0.05

0.1

J1

(a)

(b)

0.15

0.2

0.25

J1

(c)

dil , (b) shear modulus 𝜇
Figure 2: Membrane behavior: Variation of (a) dilatation Energy 𝑊m
and (c) anisotropic part 𝜂 with 𝒥1 for the model of Kumar and Parks (2014).

nm−1 and the material response is isotropic below this limit. They used FGBP and SGBP to
determine the material parameter of their bending model. In the former, only the bond angle
is included and in the latter the dihedral angle is included in the potential. Lu et al. (2009)
report the bending moduli 0.133 and 0.225 nN·nm obtained from FGBP and SGBP (Brenner,
1990; Brenner et al., 2002). The latter is much closer to existing ab-initio energy calcluations,
which determine the bending modulus at 0.238 nN·nm (Kudin et al., 2001). These bending
parameters are summarized in Tab. 3. In all the following simulations the bending modulus is
taken as 0.238 nN·nm unless otherwise mentioned. The bending stiffness of graphene is not due
to its thickness. The main source of the bending stiffness of graphene is the changing of bond
angles and dihedral angles. The model of Canham (1970) is a good candidate for modeling the
bending of graphene. The Canham bending strain energy density per unit current area can be
written as
)︀
𝑐 (︀ 2
𝑤b =
𝜅1 + 𝜅22 ,
(56)
2
where 𝑐 is a material constant. Since 𝜅1 and 𝜅2 depend on the in-plane stretch, model (56)
introduces coupling between membrane and bending deformation. The bending strain energy
per unit of current area can be transformed to the reference area as
)︀
𝑐 (︀ 2
𝑊b = 𝐽
𝜅1 + 𝜅22 .
(57)
2
The final form of the Kirchhoff stress and bending moment tensors due to bending are
]︁
[︁
𝜏b𝛼𝛽 = 𝐽 𝑐 (2𝐻 2 + 𝜅)𝑎𝛼𝛽 − 4𝑐 𝐻 𝑏𝛼𝛽 ,

(58)

and
𝑀0𝛼𝛽 = 𝑐 𝐽 𝑏𝛼𝛽 .

(59)

The elasticity tensors for bending can be found in Sauer and Duong (2017). As seen from (59),
the Canham model leads to a simple linear relationship between curvature and bending moment.
The coefficient of proportionality, 𝑐, corresponds to the bending stiffness. An advantage of the
Canham model is that it does not require the notion of a thickness. An alternative approach
is to use 3D elasticity and apply thickness integration as is considered in the formulation of
Delfani and Shodja (2016).

4

Finite element formulation

The presented graphene material model is implemented within the isogeometric rotation-free
shell formulation of Duong et al. (2017), which is briefly summarized here. This section starts
10

c [nN·nm]

FGBP
0.133

SGBP
0.225

QM
0.238

Table 3: Bending behavior: Bending stiffness according to various atomistic and quantum
models (Lu et al., 2009; Kudin et al., 2001).
with the equilibrium equation and the corresponding discretization approximations. Then,
the weak form is presented. Finally, the stiffness matrices are obtained based on standard
linearization.

4.1

Equilibrium equation

The equilibrium equation for a Kirchhoff-Love shell can be written as
𝑇 𝛼;𝛼 + 𝑓 = 𝜌 𝑣˙

∀𝑥 ∈ ℬ,

(60)

where “;” is the co-variant derivative, 𝑇 𝛼 = 𝜎 T ·𝑎𝛼 , and 𝑓 and 𝑣˙ are body force and acceleration
vectors. 𝑓 and the boundary conditions can be defined as
:= 𝑓 𝛼 𝑎𝛼 + 𝑝 𝑛
¯
:= 𝑢

∀𝑥 ∈ ℬ,
∀ 𝑥 ∈ 𝜕ℬ𝑛 ,

𝑡

¯
:= 𝑛
:= ¯𝑡

𝑚𝜏

:= 𝑚
¯𝜏

∀ 𝑥 ∈ 𝜕ℬ𝜏 ,

𝑓
𝑢
𝑛

∀ 𝑥 ∈ 𝜕ℬ𝑢 ,
(61)

∀ 𝑥 ∈ 𝜕ℬ𝑡 ,

¯ and 𝑛
¯ are prescribed position and normal unit vector at the Dirichlet boundary, and
where 𝑢
¯𝑡 and 𝑚
¯ 𝜏 are tractions and the bending moments at the Neumann boundary.

4.2

Discretization

The continua in current, ℬ , and reference, ℬ0 , configuration can be discretized with NURBS
shape functions and control points as
𝑋 = N X𝑛𝑒 ; 𝑥 = N x𝑛𝑒 ,

(62)

where N and 𝑛𝑒 are the NURBS shape function matrix and the number of control points per
element. Using the same shape functions, the following relations can be obtained
𝛿𝑥
𝑎𝛼
𝛿𝑎𝛼
𝑎𝛼,𝛽
𝑎𝛼;𝛽

=
=
=
=
=

N 𝛿x𝑛𝑒 ,
N,𝛼 x𝑛𝑒 ,
N,𝛼 𝛿x𝑛𝑒 ,
N,𝛼𝛽 x𝑛𝑒 ,
Ñ;𝛼𝛽 x𝑛𝑒 ,

(63)

where “,” is the parametric derivative and N , N,𝛼 , N,𝛼𝛽 , and Ñ;𝛼𝛽 are defined as
N

:= [𝑁1 1, 𝑁2 1, . . . , 𝑁𝑛𝑒 1] ,

N,𝛼 := [𝑁1,𝛼 1, 𝑁2,𝛼 1, . . . , 𝑁𝑛𝑒 ,𝛼 1] ,
N,𝛼𝛽 := [𝑁1,𝛼𝛽 1, 𝑁2,𝛼𝛽 1, . . . , 𝑁𝑛𝑒 ,𝛼𝛽 1] ,
Ñ;𝛼𝛽 := N,𝛼𝛽 − Γ𝛾𝛼𝛽 N,𝛾 ,
11

(64)

where 𝑁𝑖 are NURBS shape functions. Using approximations (Eqs. (62) and (63)), the weak
form can be written as
𝑛el
∑︁

(𝐺𝑒in + 𝐺𝑒int + 𝐺𝑒c − 𝐺𝑒ext ) = 0 , ∀ 𝛿x𝑒 ∈ ℬ ,

(65)

𝑒=1

where 𝑛el is the number of elements. 𝐺𝑒c is related to contact and is discussed in Appendix D.2.
𝐺𝑒in , 𝐺𝑒int , and 𝐺𝑒ext are inertial, internal and external parts and are defined as
∫︁
𝑒
T
𝑒
𝑒
𝐺in := 𝛿x𝑒 fin ; with fin = NT 𝜌 𝑣˙ d𝐴 ,
(66)
ℬ0𝑒

𝐺𝑒int

:=

𝑒
𝛿xT
𝑒 (fint𝜏

+

𝑒
fint𝑀
)

; with

∫︁

𝑒
fint𝜏

𝜏

=

𝛼𝛽

NT
,𝛼 𝑎𝛽

d𝐴 ;

𝑒
fint𝑀

:=

𝛿xT
𝑒

(︀

𝑒
fext0

=

𝑀0𝛼𝛽 ÑT
;𝛼𝛽 𝑛 d𝐴 ,

ℬ0𝑒

ℬ0𝑒

𝐺𝑒ext

∫︁

+

𝑒
fext𝑝

+

𝑒
fext𝑡

+

𝑒
fext𝑚

)︀

.

(67)
(68)

The terms in Eq. (68) are
𝑒
fext0
=

𝑒
fext𝑝
=
𝑒
fext𝑡

∫︁

NT 𝑓 0 d𝐴 ,

ℬ∫︁0𝑒

NT 𝑝 𝑛 d𝑎 ,

𝑒
ℬ∫︁

=

(69)
NT 𝑡 d𝑠 ,

𝜕𝑡 ℬ𝑒 ∫︁

𝛼
NT
,𝛼 𝜈 𝑚𝜏 𝑛 d𝑠 ,

𝑒
fext𝑚
= −

𝜕𝑚 ℬ 𝑒

where along 𝜕𝑚 ℬ bending Neumann BCs are applied. 𝜈 𝛼 are the contra-variant components of
the normal vector on the boundary 𝜕ℬ.

4.3

Stiffness matrix

The stiffness matrix appearing in the linearization of the discretized weak form, can be written
as
k𝑒 = k𝑒𝜏 𝜏 + k𝑒𝜏 𝑀 + k𝑒𝑀 𝜏 + k𝑒𝑀 𝑀 + k𝑒𝜏 + k𝑒𝑀 + k𝑒extp + keextm + k𝑒𝑐 ,
(70)
where the material stiffness matrices are defined by
∫︁
k𝑒𝜏 𝜏 :=
𝑐𝛼𝛽𝛾𝛿 NT
,𝛼 (𝑎𝛽 ⊗ 𝑎𝛾 ) N,𝛿 d𝐴 ,
ℬ0𝑒

k𝑒𝜏 𝑀

∫︁
:=

𝑑𝛼𝛽𝛾𝛿 NT
,𝛼 (𝑎𝛽 ⊗ 𝑛) Ñ;𝛾𝛿 d𝐴 ,

ℬ0𝑒

k𝑒𝑀 𝜏

(71)

∫︁
:=

𝑒

𝛼𝛽𝛾𝛿

ÑT
;𝛼𝛽

(𝑛 ⊗ 𝑎𝛾 ) N,𝛿 d𝐴 ,

ℬ0𝑒

k𝑒𝑀 𝑀

∫︁
:=

𝑓 𝛼𝛽𝛾𝛿 ÑT
;𝛼𝛽 (𝑛 ⊗ 𝑛) Ñ;𝛾𝛿 d𝐴 ,

ℬ0𝑒

12

where the elasticity tensors are given in Appendix B and Sauer and Duong (2017). Furthermore,
the geometrical stiffness matrices are defined as
∫︁
𝛼𝛽
𝑒
NT
N,𝛽 d𝐴 ,
k𝜏 :=
,𝛼 𝜏
(72)
ℬ0𝑒

k𝑒𝑀 := k𝑒𝑀 1 + k𝑒𝑀 2 + (k𝑒𝑀 2 )T ,
with
k𝑒𝑀 1

∫︁
:= −

𝑏𝛼𝛽 𝑀0𝛼𝛽 𝑎𝛾𝛿 NT
,𝛾 (𝑛 ⊗ 𝑛) N,𝛿 d𝐴 ,

ℬ0𝑒

k𝑒𝑀 2 := −

∫︁

(73)

(74)
𝛾
𝑀0𝛼𝛽 NT
,𝛾 (𝑛 ⊗ 𝑎 ) Ñ;𝛼 𝛽 d𝐴 .

ℬ0𝑒

The external tangent matrices k𝑒extp and k𝑒extm can be found in Sauer et al. (2014) and Duong
et al. (2017), and the contact stiffness matrix k𝑒𝑐 is given in Appendix D.2.

5

Elementary model behavior

In this section, the elementary behavior of the model under some simple deformation states is
investigated. Pure dilatation, uniaxial stretch, pure shear and pure bending are considered as
test cases, and the FE implementation is verified by analytical solutions.

5.1

Pure dilatation

The numerical model consists of a square sheet stretched uniformly in two perpendicular directions. The boundary conditions and the lattice orientation are shown in Fig. 3a. The surface
tension4 𝛾 := 12 tr (𝜎) is obtained with the LDA and the GGA parameters and the results are
compared in Fig. 4. The current continuum model is compared with quantum results from the
literature. The difference of LDA and GGA results is small for infinitesimal strains. But, the
response is different for finite strains, where LDA shows a stiffer response. It should be noted
that the material is unstable after the maximum value of the surface tension has been reached.
The shear modulus is zero at this point and become negative afterwards. Displacement control
is needed to go beyond that point.

5.2

Uniaxial stretch

Uniaxial stretching is considered along different directions, i.e. different chiralities. The model
is stretched in one direction and restrained in the perpendicular direction (Fig. 3b). In order to
assess the FE solution, the Cartesian stress components 𝜎11 and 𝜎22 are examined. They are
the normal stress components in the direction of the stretch and perpendicular to it. They are
obtained analytically in Appendix C. The results of stretching the model in the armchair and
the zigzag directions are presented in Figs. 5 and 6. The current continuum model is compared
with quantum results from the literature. Like for dilatational loading, a material instability
occurs after the maximum stress has been reached and the elasticity tensor loses its ellipticity
at this point.
4

The surface tension is equal to the negative in-plane pressure.

13

L

L
uy=0

λ2=λ0

λ2=1

u1=0

2

2

θ

λ1=λ0

L

u1=0

θ

1

(a)

λ1=λ0

1

(b)

u2=0

u2=0

L
λ2=1/λ0

u1=0

2

θ

λ1=λ0

L

1

(c)

(d)
u2=0

Figure 3: Elementary membrane deformations: (a) Pure dilatation; (b) uniaxial stretch; (c)
pure shear; (d) considered FE mesh; linear Lagrangian elements are used and the nodes are
indicated by ∙. Angle 𝜃 denotes the orientation of the graphene lattice with respect to the
direction of larger stretch 𝜆1 , and 𝜆0 > 1 is gradually increased.

14

L

35
30

. [N/m]

25
20
15
10
Wei et al. (2009)
Shodja et al. (2017)
Kumar and Parks (2014) (GGA)
Kumar and Parks (2014) (LDA)

5
0
-5
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

J1

Figure 4: Pure dilatation: Surface tension. The current continuum model is based on the
quantum mechanical results of Kumar and Parks (2014). Other quantum mechanical results
are available in Wei et al. (2009) and Shodja et al. (2017).
40

5

35
4

25

<22 [N/m]

<11 [N/m]

30

20
15
Wei et al. (2009)
Shodja et al. (2017)
Kumar and Parks (2014) (GGA)
Kumar and Parks (2014) (LDA)

10
5

3

2
Wei et al. (2009)
Shodja et al. (2017)
Kumar and Parks (2014) (GGA)
Kumar and Parks (2014) (LDA)

1

0

0
1

1.05

1.1

1.15

1.2

1.25

1

61

1.05

1.1

1.15

1.2

1.25

61

(a)

(b)

Figure 5: Uniaxial stretch in the armchair direction: (a) Stress in the stretched direction, 𝜎11 ;
(b) stress in the perpendicular direction, 𝜎22 . The current continuum model is based on the
quantum mechanical results of Kumar and Parks (2014). Other quantum mechanical results
are available in Wei et al. (2009) and Shodja et al. (2017).
Next, the model is stretched in different directions, relative to the armchair direction, and the
stress variation is examined for a set of stretch ratios in Fig. 7. As can be expected, the material
response is periodic by 𝜋/3. Finally, the standard deviation and average of the stresses in one
period is calculated for the different stretch ratios. The average and standard deviation can be
calculated as
∫︁
3 𝜋/3
(75)
𝜎𝑖𝑗 d𝜃 ,
Avg(𝜎𝑖𝑗 ) = ⟨𝜎𝑖𝑗 ⟩ =
𝜋 0
∫︁ 𝜋/3 √︂
3
(76)
Std(𝜎𝑖𝑗 ) =
(𝜎𝑖𝑗 − ⟨𝜎𝑖𝑗 ⟩)2 d𝜃 .
𝜋
0
The standard deviation is normalized by the average stress and shown versus 𝜆1 in Fig. 7c.
The graph can be interpreted as the error caused by assuming isotropy. This error can reach up
to 80% of the average response. So, an anisotropic model is necessary to capture the material
behavior for finite deformations.

15

50

6
5

40

<22 [N/m]

<11 [N/m]

4
30

20

3
2
1

Wei et al. (2009)
Shodja et al. (2017)
Kumar and Parks (2014) (GGA)
Kumar and Parks (2014) (LDA)

10

Wei et al. (2009)
Shodja et al. (2017)
Kumar and Parks (2014) (GGA)
Kumar and Parks (2014) (LDA)

0

0

-1
1

1.05

1.1

1.15

1.2

1.25

1

1.05

1.1

61

1.15

1.2

1.25

61

(a)

(b)

Figure 6: Uniaxial stretch in the zigzag direction: (a) Stress in the stretched direction, 𝜎11 ;
(b) stress in the perpendicular direction, 𝜎22 . The current continuum model is based on the
quantum mechanical results of Kumar and Parks (2014). Other quantum mechanical results
are available in Wei et al. (2009) and Shodja et al. (2017).

45

61 = 1:05

61 = 1:1

61 = 1:15

61 = 1:2

6

61 = 1:25

61 = 1:15

61 = 1:2

80

61 = 1:25

Std(<)=Avg(<) [%]

<22 [N/m]

<11 [N/m]

61 = 1:1

5

35
30
25

4.5
4
3.5

20
15

61 = 1:05

<11
<22
<11
<22

5.5

40

GGA
0

40

(GGA)
(GGA)
(LDA)
(LDA)

40

20

3

LDA ! " ! "
20

60

60

2.5

GGA
0

LDA ! " ! "
20

40

3 [degree]

3 [degree]

(a)

(b)

60

0

1

1.05

1.1

1.15

1.2

1.25

61

(c)

Figure 7: Uniaxial stretch: Dependence of the stress components on stretch 𝜆1 and angle 𝜃, (a)
𝜎11 ; (b) 𝜎22 . (c) Dependence of the standard deviation of the stress on stretch 𝜆1 .

16

100

0

GGA
LDA

80

<22 [N/m]

<11 [N/m]

-50
60

40

-100
20
GGA
LDA
0

-150
0

0.02

0.04

0.06

0.08

0.1

0

0.12

0.02

0.04

J2

0.06

0.08

0.1

0.12

J2

(a)

(b)

60

0

50

-20

40

-40

<22 [N/m]

<11 [N/m]

Figure 8: Pure shear with pull in the armchair and compression in the zigzag direction (𝜃 = 0
in Fig. 3): (a) Stress in the pull direction, 𝜎11 ; (b) stress in the compression direction, 𝜎22 .

30
20

GGA
LDA

-60
-80

GGA
LDA

10

-100

0

-120
0

0.02

0.04

0.06

0.08

0.1

0.12

0

J2

0.02

0.04

0.06

0.08

0.1

0.12

J2

(a)

(b)

Figure 9: Pure shear with pull in the zigzag and compression in the armchair direction (𝜃 = 𝜋/6
in Fig. 3): (a) Stress in the pull direction, 𝜎11 ; (b) stress in the compression direction, 𝜎22 .

5.3

Pure shear

Pure shear is considered in different directions, i.e. different chiralities. Therefore, the specimen
is pulled in one direction and compressed in the perpendicular direction by 𝜆1 = 𝜆0 > 1 and
𝜆2 = 1/𝜆0 such that 𝐽 = 𝜆1 𝜆2 = 1 (Fig. 3c). The Cartesian stress components 𝜎11 and 𝜎22 are
introduced and analytically derived in Appendix C. The stress variation versus 𝒥2 , in the pull
and compression direction, for the armchair and zigzag direction, are shown in Figs. 8 and 9.
Next, the loading is applied in different directions, relative to the armchair direction, and stress
variation is demonstrated for a set of stretch ratios and GGA and LDA parameters in Fig. 10.
The graphes show the similar periodic behavior as noted in the previous section. Finally, the
standard deviation and average of the stresses in one period are calculated for the different
stretch ratios. The standard deviation is normalized by the average stress and it is shown
versus 𝒥2 in Fig. 10c. Like for uniaxial stretch, the graph can be interpreted as the error caused
by assuming isotropy. As shown, this error can reach up to 17 percent. The anisotropic behavior
is increasing fast even for low values of 𝒥2 .

17

61 = 1:05

61 = 1:1

61 = 1:15

61 = 1:2

0

61 = 1:25

61 = 1:05

61 = 1:1

61 = 1:15

61 = 1:2

<22 [N/m]

<11 [N/m]

80
60
40
GGA

20
0

-50
GGA

20

61 = 1:25

Std(<)=jAvg(<)j [%]

100

LDA ! " ! "

-100

<11
<22
<11
<22

15

(GGA)
(GGA)
(LDA)
(LDA)

10

5

LDA ! " ! "
20

40

60

-150

3 [degree]

0

20

40

60

0

0

0.02

3 [degree]

(a)

0.04

0.06

0.08

0.1

0.12

J2

(b)

(c)

Figure 10: Pure shear: Dependence of the stress components on stretch 𝜆1 and angle 𝜃, (a) 𝜎11 ;
(b) 𝜎22 . (c) Dependence of the standard deviation of the stress on 𝒥2 .

(a)

(b)

Figure 11: Pure bending: (a) boundary conditions (the green, blue and red triangles indicate restriction of displacement dofs in the longitudinal, lateral and perpendicular directions,
respectively; the arrows indicate prescribed rotations); (b) deformed geometry.

5.4

Pure bending

In this section, the bending formulation is validated by considering pure bending. The boundary
conditions are shown in Fig. 11a. The specimen’s width is restrained and it is stress free in
the longitudinal direction. The details of the analytical solution for pure bending can be found
in Sauer and Duong (2017). The resultant force along the longitudinal direction, 𝑁11 (see
Eq. (43)), is zero for pure bending. 𝑁11 = 0 is solved by the standard secant method to obtain
the analytical solution. The final geometry is plotted in Fig. 11b. The convergence of the
total energy is obtained with mesh refinement (Fig. 12a). The bending strain energy per unit
reference area is shown in Fig. 12b. The relation between bending moment and curvature is
linear.

6

Numerical examples

In the previous section, the proposed model is verified by standard tests. In the current section,
several numerical benchmark examples are solved and compared with results from the literature.
Considered are five examples: indentation and peeling of graphene, and torsion, bending and
uniaxial stretch of carbon nanotubes (CNTs).

18

10

0.14

Total energy per unit area [nN/nm]

Relative error of strain energy

100

-2

10-4

10-6

10

-8

100

101

102

103

104

0.12
0.1
0.08
0.06
0.04
0.02
0

105

0

0.2

0.4

0.6

0.8

1

1.2

5 [1/nm]

Number of elements

(a)

(b)

Figure 12: Pure bending: (a) Error of strain energy relative to the finest mesh (384 × 128
quadratic NURBS elements); (b) total energy per unit reference area.

6.1

Indentation of graphene

In order to compare the developed model with experimental results, an indentation problem
is considered in this section. The experimental results are taken from a setup with multiple
cavities. Relative position of these cavities can affect experimental measurements. In current
numerical simulations, only one cavity is modeled and frictional effects are neglected. In the
following subsections, first the specimen details, boundary conditions, and loading are discussed.
Finally, a parameter study is conducted for the adhesion strength and indentor radius. Kumar
and Parks (2015) investigated the same problem using the Morse potential to model substrate
interactions and ABAQUS (Abaqus, 2016) in an explicit dynamic manner. In the current work,
the Lennard-Jones potential and a quasi-static formulation is used.
6.1.1

Specimen details, boundary conditions and substrate adhesion

The problem setup consists of a sheet of graphene on a SiO2 substrate with a micro-cavity.
The boundary conditions and geometry of the substrate are depicted in Fig. 13. For efficiency,
only one quarter of the circular sheet is discretized. Thus, additional symmetry BCs have to
be applied along the symmetry planes. In addition, the outer sheet boundary is fixed in all
directions. The edge of the cavity is smoothed with a fillet radius of 50 nm. The adhesion plays
an important role in the initial stress, stiffness and overall response of the sheet. A relaxation
step is considered before the indentation loading. The adhesion between sheet and substrate
is modeled via van der Waals (vdW) interaction (Appendix D). In the relaxation step, the
adhesion parameter is increased from zero to its final value. At each step, a standard Newton
Raphson iteration is utilized to determine equilibrium. The displacement and stress distribution
of the relaxed geometry are shown in Fig. 14. This stress corresponds to a pre-stress within the
graphene sheet that makes the structure stiffer.
6.1.2

Indentation: Load step

10000 quadratic isogeometric finite elements5 are used over a rectangular domain with the same
number of elements in x and y directions. Before the indentation phase, a relaxation step is
5

The continuum model contains 122,412 nodes while the corresponding atomistic system has about 12 million
atoms, i.e. about 100 times more.

19

RI

RI

R1

R1

R2
(a)

(b)

(c)

Figure 13: Indentation setup: Outer boundary conditions and symmetry boundary conditions
for (a) adhesive substrate and (b) substrate without adhesion; (c) substrate geometry. 𝑅𝐼 is
the indentor radius, 𝑅1 the cavity radius and 𝑅2 the fillet radius of the cavity.

(a)

(b)

Figure 14: Substrate adhesion of graphene: (a) Displacement in the perpendicular direction to
the graphene sheet after relaxation (in nm); (b) tr(𝜎) [N/m]. The adhesion strength Γ and R1
(see Fig. 13a) are set to 0.45 N/m and 500 nm, respectively. The maximum of tr(𝜎) and the
displacement are 4.95 N/m and 11.75 nm, respectively. In the simulation, all displacement dofs
with a radial position greater than 𝑅1 are fixed.

20

(a)

(b)

Figure 15: Graphene indentation: (a) Indented profile (colored by tr(𝜎) [N/m]); (b) zoom of
(a). Adhesion parameter Γ, cavity radius 𝑅1 and indentor radius 𝑅𝐼 are 0.45 N/m, 500 nm and
35 nm, respectively.
conducted. A wide range for adhesion strengths (Γ = 0.1 ∼ 0.45 N/m) can be found in the
literature discussed in Sec. 1. To study the influence of adhesion, the indentation simulation
is conducted for a set of different Γ. The indented geometry is presented in Fig. 15. The
indentor contact force is measured and compared with experimental results in Fig. 16a. A
higher adhesion strength results in a stiffer structure and a higher indentor reaction force.
Next, to study the effect of the indentor size, a set of simulations are conducted for a series
of indentors (Fig. 16b). They have the same force-displacement curve in the beginning of
indentation. But, they have a different response for large indentation. In addition, a larger
indentor has a larger reaction force.
Finally, the problem is solved with zero bending stiffness (pure membrane model) and compared
with the shell model for the substrate with and without adhesion. It is shown that the bending
stiffness does not play an important role for the force-displacement curve (Fig. 17). In addition,
in the cases without adhesion, it is assumed that the edge of the substrate cavity is not filleted
and 𝑅1 is selected to be consistent with the adhesive substrate setup (see Figs. 13a and 13b).

6.2

Peeling of graphene

Peeling between graphene and an adhesive substrate is considered next. A rectangular graphene
strip is considered that is initially positioned in equilibrium distance on a SiO2 substrate. The
strip is restrained on one side and the other side is displaced in the perpendicular direction to
the substrate (𝑢𝑧 ) (Fig. 18a). The first 20 percent of the length of the strip are not absorbed to
the substrate (the adhesion parameter is assumed to be zero in this area). Quadratic NURBS
elements with the same element size in both directions are used in the simulation. The deformed
geometry is plotted in Fig. 18b. The convergence of the peeling force with mesh refinement is
21

2000

2500
!=0.45 N/m
!=0.14 N/m
!=0.10 N/m
Experiment

2000

Force [nN]

Force [nN]

1500

RI =16.5 nm
RI =25 nm
RI =35 nm
Experiment

1000

1500

1000

500
500

0

0
0

20

40

60

80

100

120

0

20

Displacement [nm]

40

60

80

100

120

140

Displacement [nm]

(a)

(b)

Figure 16: Graphene indentation, force-displacement curve of the indentor: (a) Comparison between different adhesion parameters (for fixed 𝑅1 = 500 nm and 𝑅𝐼 = 16.5 nm); (b) comparison
between different indentor radii (for fixed Γ = 0.1 N/m and 𝑅1 = 500 nm). The experimental
results are taken from Kumar and Parks (2015).

2000
Membrane (!=0)
Shell (!=0)
Membrane (!=0.1 N/m)
Shell (!=0.1 N/m)
Experiment

Force [nN]

1500
60
1000 40
20
500

0
0

10

20

0
0

20

40

60

80

100

120

Displacement [nm]

Figure 17: Graphene indentation: Comparison between membrane and shell models considering
various values for Γ (for fixed 𝑅1 = 500 nm and 𝑅𝐼 = 16.5 nm). The experimental results are
taken from Kumar and Parks (2015).

22

shown Fig. 19b. Fig. 19a shows the peeling force against the peeling displacement. The peeling
force exhibits an increasing and a decreasing part. They capture the elastic deformation and
adhesive debonding of the strip, respectively.

(a)

(b)

Figure 18: Graphene peeling: (a) The peeling setup (BC and loading). The size of the strip is
50nm × 10nm. (b) The deformed geometry of the peeled strip colored by tr(𝜎) [N/m].

100

Relative error of peeling force

Reaction force Fz [nN]

0.01

0.008

0.006

0.004

0.002

0

5

10

15

20

10-2

10-4

10-6
101

0
25

16 nm displacement
18 nm displacement
20 nm displacement

uz (nm)

102

103

104

Number of elements

(a)

(b)

Figure 19: Graphene peeling: (a) Peeling force vs. peeling displacement. (b) Error of the peeling
force for various peeling displacements relative to the finest mesh (90 × 60 quadratic NURBS
elements). The LDA and FGBP parameter sets are used together with Γ = 0.14 N/m.

6.3

Deformation of carbon nanotubes (CNTs)

In this section, the proposed shell model is applied to CNTs. First, the relaxation process is
discussed. Then, the CNT behavior under torsional and bending loading is investigated. The
results for torsion are compared and validated with atomistic results from the literature. In
both loading scenarios, the buckling of shell walls is captured and the buckling load is accurately
determined from the ratio of the membrane energy to the total energy. This energy ratio is
determined from the strain energies introduced in Sec. 3.

23

Strain energy [eV/atom]

0.2

0.15

0.1

0.05
Atomistic armchair
Atomistic zigzag
FE shell model with GGA
0
0

1

2

3

4

5

1/R [1/nm]

Figure 20: CNT relaxation: the strain energy per atom obtained after relaxation. 𝑅 denotes the
radius of the CNT according to Eq. (77). For the finite element simulation different chiralities
are used. The atomistic results are taken from Arroyo and Belytschko (2004a).

Figure 21: CNT chirality: CNT(𝑛,0) and CNT(𝑛,𝑛) indicate zigzag and armchair CNTs. 𝑎𝛼
are the primary lattice vectors.
6.3.1

Relaxation of CNTs

The CNT model needs to be relaxed in order to obtain a stress-free initial configuration (Favata
et al., 2016). The internal energy is minimized in this relaxation process. During relaxation
the CNT deforms radially while the length remains almost constant. The minimized strain
energy per atom for CNTs with different chirality is compared with results from the literature
in Fig. 20. The radius for CNT(𝑛,𝑚) can be calculated as
√
3𝑎cc √︀ 2
(77)
𝑅 =
𝑛 + 𝑛 𝑚 + 𝑚2 ,
2𝜋
where the chirality parameters 𝑛 and 𝑚 indicate the number of the unit cells along the primary
lattice vectors (Fig. 21). 𝑎cc = 0.142 nm is the equilibrium length of the carbon-carbon bond.
6.3.2

Torsion of CNT

In this section, the behavior of an initially relaxed CNT under torsional loading is investigated.
The loading can be conducted either by restraining the axial length of CNT (Fig. 22a) or letting
24

the CNT deform in the axial direction (Fig. 22b). A torsion angle is applied to both ends of
CNT, and the simulation is conducted with and without imperfection. The imperfection is
applied as a torque6 in the middle of the CNT (Fig. 22c). The convergence is studied for the
total energy (Fig. 23). The cross section of the model without imperfection remains cylindrical
during torsion, i.e. it does not buckle, and the results converge faster in comparison to the
model with imperfection. The strain energy per atom of CNT is plotted for the perfect and
imperfect structures. The results from the presented model are compared to results from an
atomistic simulation with a perfect crystal in Fig. 24a. In the loading process, the CNT with
the imperfection buckles, but the exact point of the instability is not easily found from the
variation of the total energy. To estimate the instability point more accurately, the ratio of the
membrane energy to the total energy is used. The sharp decreasing of energy ratio determines
the buckling angle around (𝜃b = 14∘ ) (Fig. 24b). The strain energy is measured relative to the
relaxed geometry. Next, a set of simulations is conducted for the different charities and the
results are presented in Fig. 25. The deformed and sliced geometries are compared for different
torsion angles (Fig. 26).

F

F

(a)

(b)

(c)

Figure 22: CNT twisting: BC for (a) axially fixed case; (b) axially free case; (c) imperfection
torque.

Relative error of energy

100

3=15/
3=25/
3=35/

10-1

10-2

10-3

10-4

10-5
102

103

104

Number of elements

Figure 23: CNT twisting: Error of strain energy relative to the finest mesh (100 × 100 quadratic
NURBS elements). CNT(12,6) with the length 6.74 nm is used for the axially fixed case
(Fig. 22a).

6.3.3

Bending of CNT

In this section, the behavior of a CNT under bending loading is investigated. CNT is bent by
applying a bending angle equally at both ends (Fig. 27a). During bending, the end faces of the
6

by applying the force couple (F=1 nN)

25

1
Perfect model
Imperfect model

Atomistic
Perfect model
Imperfect model

0.25

Membrane energy/total energy

Strain energy [eV/per atom]

0.3

0.2
0.15
0.1
0.05
0

0.9
0.8
0.7
0.6
0.5
0.4

0

10

20
30
Torsion angle [degree]

40

50

0

10

20
30
Torsion angle [degree]

(a)

40

50

(b)

Figure 24: CNT twisting: (a) Strain energy per atom; (b) ratio of the membrane energy to the
total energy (measured relative to relaxed geometry). CNT(12,6) with the length 6.74 nm is
used for the axially fixed case (Fig. 22a). The perfect and imperfect models are used for the
simulation. The atomistic results are taken from Sun and Li (2011).

Strain energy [eV/per atom]

0.02
CNT(10,10)
CNT(14,7)
CNT(17,0)

0.015

0.01

0.005

0
0

2

4

6
8
10
Torsion angle [degree]

12

14

16

Figure 25: CNT twisting: Strain energy per atom for different chiralities of CNTs with the
length of 10 nm. CNTs are selected such that their radii are close to each other. The BC is set
to the axially free case (Fig. 22b) and imperfection is applied to the model.
CNT are assumed to be rigid and remain planar, and the CNT can deform in axial direction in
order to attain a state of pure bending loading. The strain energy convergence is studied with
mesh refinement (Fig. 27b). The strain energy per atom is presented in Fig. 28a for perfect
and imperfect structures and the buckling point is determined from the sharp variation in the
ratio of the membrane energy to the total energy (Fig. 28b). It is shown that for the imperfect
model, buckling occurs earlier. tr(𝜎) is compared at different bending angles for the perfect
and imperfect models (Fig. 29). Finally, the models are cut open along the axial direction to
compare the buckled geometry of the perfect and imperfect CNTs (Fig. 30).

6.4

Stretch of chiral CNTs

The deformation of chiral CNTs under axial loading is investigated here. The deformation
includes a uniform extension and twisting along the axial direction of the CNT. The anisotropy
26

(b)

(c)

(d)

75 percent of CNT length 50 percent of CNT length

(a)

Figure 26: CNT twisting: Comparison of cross sections and tr(𝜎) [N/m] at the torsion angles
(a) 𝜃 = 15∘ , (b) 𝜃 = 25∘ , (c) 𝜃 = 35∘ and (d) 𝜃 = 46.5∘ . CNT(12,6) with the length 6.74 nm is
used. The BCs are applied according to the axially fixed case (Fig. 22b).

Relative error of energy

100

10-1

10-2

10-3

10-4
102

3=10/
3=20/
3=30/
103

104

105

Number of elements

(a)

(b)

Figure 27: CNT bending: (a) BC; (b) Error of strain energy relative to the finest mesh (200×200
quadratic NURBS elements). CNT(10,10) with the length 10 nm and imperfection is used.
27

1
0.9

0.05

Membrane energy/total energy

Strain energy [eV/per atom]

0.06

0.04
0.03
0.02
0.01

Perfect model
Imperfect model
5

10

15
20
25
Bending angle [degree]

30

0.7
0.6
0.5
0.4
0.3

0
0

0.8

Perfect model
Imperfect model

0.2

35

0

(a)

5

10

15
20
25
Bending angle [degree]

30

35

(b)

Figure 28: CNT bending: (a) Comparison of strain energy per atom in perfect and imperfect
cases; (b) comparison of the ratio of the membrane energy to the total in perfect and imperfect
cases. CNT(10,10) with the length 10 nm is used.

Perfect model

Imperfect model

(a)

(b)

(c)

(d)

Figure 29: CNT bending: Comparison of tr(𝜎) [N/m] of the perfect and imperfect models for
the bending angles (a) 𝜃 = 10∘ , (b) 𝜃 = 20∘ , (c) 𝜃 = 30∘ and (d) 𝜃 = 34∘ . CNT(10,10) with the
length 10 nm is used.

28

(a)

(b)

Figure 30: CNT bending: Cross sections of (a) the perfect and (b) the imperfect models
(𝜃 = 34∘ ). CNT(10,10) with the length 10 nm is used.
of graphene is the source of this twisting. Three chiralities are selected, and the variation of
the twist angle against the stretch is presented in Fig. 31. The results are compared with the
results of Delfani and Shodja (2013). There is very good agreement for large tube radii, but
for small tube radii differences appear: Delfani and Shodja (2013) find an initial twist at zero
stretch that they argue appears from relaxing the initial energy of the original tube (that is a
rolled graphene sheet with non-zero bending energy). The proposed model, on the other hand,
relaxes the initial bending energy by a change of radius and does not predict an initial twist.

7

Conclusion

A new shell model is developed based on a new anisotropic membrane model and a new isogeometric shell formulation. The model is verified with standard biaxial, uniaxial and bending
tests. It is used for the modeling of indentation and peeling of a graphene sheet, and the torsion, bending and stretching of a CNT. The new model is computationally more efficient than
atomistic methods when the problem size reaches micrometer scale. In addition, the interatomic
potentials presently used within the exponential Cauchy-Born rule (ECB) tend to underestimate the elastic modulus and bending stiffness obtained from quantum mechanics (QM). Also,
molecular dynamics (MD) simulation have the same drawback as the ECB rule. The current
model uses existing QM data to calibrate the material constants. The material response of the
current model is therefore closer to existing ab-initio and experimental results. Thus the current
model improves both accuracy and efficiency. Finally, the current model is fully nonlinear and
can handle large deformations, buckling and postbuckling. The buckling points in torsional and
bending loading are determined from the energy ratio, and the postbuckling behavior of CNTs
is simulated.

Acknowledgment
Financial support from the German Research Foundation (DFG) through grant GSC 111, is
gratefully acknowledged. The authors also thank Maximilian Harmel for checking the manuscript
carefully.

29

angle of twist [Rad/nm]

0.02

0.01

0
Proposed model
-0.01 Delfani and Shodja (2013)

!"!"

-0.02
CNT(8,2)

CNT(12,3)

CNT(16,4)

-0.03
1

1.05

1.1

1.15

1.2

1.25

61

(a)
0.02

0.03

angle of twist [Rad/nm]

angle of twist [Rad/nm]

0.02
0.01
0
Proposed model
-0.01 Delfani and Shodja (2013)

!"!"

-0.02
CNT(6,2)

CNT(9,3)

0.01

0
Proposed model
-0.01 Delfani and Shodja (2013)

-0.02

CNT(12,4)

CNT(6,3)

-0.03

!"!"

CNT(8,4)

CNT(10,5)

-0.03
1

1.05

1.1

1.15

1.2

1.25

1

1.05

1.1

61

1.15

1.2

1.25

61

(b)

(c)

Figure 31: Stretch of chiral CNTs: Twist angle vs. stretch according to the proposed model
and Delfani and Shodja (2013) for various chirality angles. √(a) 𝜙 = 12.82∘ ; (b) 𝜙 = 13.90∘ and
(c) 𝜙 = 19.11∘ . The chirality angle is defined as tan(𝜙) := 3 𝑚/(2𝑛 + 𝑚).

30

A

Tensorial derivative and curvilinear description of logarithmic strain

In this section, the derivative of eigenvalues and a formulation for the tensorial derivative of the
logarithmic strain are given. Then, the logarithmic strain is described in curvilinear coordinates.
The derivative of eigenvalues can be written as (Itskov, 2015)
𝜕Λ𝑖
𝜕Λ𝑖
=
𝐴𝛼 ⊗ 𝐴𝛽 = 𝑃𝑖𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ;
𝜕𝐶
𝜕𝑎𝛼𝛽
where 𝑃𝑖𝛼𝛽 are

with 𝑖 = 1, 2 ,

𝑖 𝐴𝛿𝛽 .
𝑃𝑖𝛼𝛽 = 𝐴𝛼𝛾 𝑃𝛾𝛿

(A.1)

(A.2)

The derivative of the logarithmic strain can be written as (Itskov, 2015)
2
∑︁
𝜕𝐸 (0)
=
𝑓𝑖𝑗 [P𝑖 ⊗ P𝑗 ]s ,
𝜕𝐶

(A.3)

𝑖,𝑗=1

where symmetrization operator (∙)s , transpose operator (∙)t and 𝑓𝑖𝑗 are defined as
𝒜s :=

)︀
1 (︀
𝒜 + 𝒜t ,
2

(𝑎 ⊗ 𝑏 ⊗ 𝑐 ⊗ 𝑑)t := (𝑎 ⊗ 𝑐 ⊗ 𝑏 ⊗ 𝑑) ,
⎧ 1
⎪
⎪
⎨ 2 𝜆2 ; if 𝑖 = 𝑗 ,
𝑖
𝑓𝑖𝑗 :=
ln
𝜆
⎪
𝑖 − ln 𝜆𝑗
⎪
; if 𝑖 ̸= 𝑗 .
⎩
2
𝜆𝑖 − 𝜆2𝑗

(A.4)
(A.5)

(A.6)

It is worth to note that there is another transpose operator (∙)T defined as
(𝑎 ⊗ 𝑏 ⊗ 𝑐 ⊗ 𝑑)T := (𝑏 ⊗ 𝑎 ⊗ 𝑑 ⊗ 𝑐) .

(A.7)

Next, 𝐸 (0) can be expressed in the curvilinear coordinate as
(0)

𝐸 (0) = 𝐸𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ,
with
(0)
𝐸𝛼𝛽

:=

2
∑︁

𝑖
ln 𝜆𝑖 𝑃𝛼𝛽
.

(A.8)

(A.9)

𝑖=1

In addition, the co-variant and contra-variant components of the deviatoric logarithmic strain
can be written as
(︁
)︁
(0)
1 − 𝑃2
𝐸dev 𝛼𝛽 := ln 𝜆 𝑃𝛼𝛽
,
(A.10)
𝛼𝛽
(0) 𝛼𝛽

𝐸dev

(︁
)︁
:= ln 𝜆 𝑃1𝛼𝛽 − 𝑃2𝛼𝛽 .

31

(A.11)

B

Derivation of membrane constitution for curvilinear coordinates

In this section, the Kirchhoff stress and corresponding elastic tensors are derived. They are
related to the membrane energy. The repeated and distinct eigenvalue cases are investigated
and corresponding constitutive laws are given in both cases. Some relations are repeated to
make it easier to follow the derivation.

B.1

Distinct eigenvalues

The case of distinct eigenvalues is considered in this section. The Kirchhoff stress tensor can
be obtained by the chain rule as
(︂
)︂
dil
dev
dev 𝜕𝒥
dev 𝜕𝒥
𝜕𝑊m
𝜕𝑊m
𝜕𝒥1
𝜕𝑊m
𝜕𝑊m
2
3
𝛼𝛽
𝜏m = 2
+2
+2
+2
,
(B.1)
𝜕𝒥1
𝜕𝒥1
𝜕𝑎𝛼𝛽
𝜕𝒥2 𝜕𝑎𝛼𝛽
𝜕𝒥3 𝜕𝑎𝛼𝛽
with
2

dil
𝜕𝑊m
= 2𝜀𝛼
^ 2 𝜖a 𝑒−𝛼^ 𝜖a ,
𝜕𝒥1

2

dev
𝜕𝑊m
^
= 4 𝒥2 𝜇′ (𝜖a ) + 2 𝒥3 𝜂 ′ (𝜖a ) = −4 (𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a + 𝜂1 𝜖a 𝒥3 ) ,
𝜕𝒥1

dev
𝜕𝑊m
= 4 𝜇(𝜖a ) ,
2
𝜕𝒥2

2

(B.2)

dev
𝜕𝑊m
= 2 𝜂(𝜖a ) ,
𝜕𝒥3

𝜕𝒥1
1
= 𝑎𝛼𝛽 ,
𝜕𝑎𝛼𝛽
2
(︂
)︂
𝜕𝒥2
1
1 𝛼𝛽
1 𝛼𝛽
=
ln 𝜆 ,
𝑃 − 2 𝑃2
𝜕𝑎𝛼𝛽
2 𝜆21 1
𝜆2

(B.3)

where 𝜇′ and 𝜂 ′ are defined as
𝜕𝜇
,
𝜕𝜖a
𝜕𝜂
𝜂 ′ :=
.
𝜕𝜖a

𝜇′ :=

(B.4)

Using Eq. (44.3) and the chain rule, the derivative of 𝒥3 can be written as
𝜕𝒥3
1 𝜕(𝛾𝜃3 ) 𝜕𝐸 (0)
:=
:
,
𝜕𝐶
8 𝜕𝐸 (0)
𝜕𝐶
with
𝑇 (0) =

𝜕(𝛾𝜃3 )
𝜕𝐸 (0)

^ +𝑎^ 𝑁
^ .
= 𝑎𝑀^ 𝑀
𝑁

(B.5)

(B.6)

^ and 𝑁
^ can be written in curvilinear coordinates as
𝑀
^ =𝑀
^ 𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ;
𝑀
^ =𝑁
^𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ;
𝑁

^ · 𝐴𝛽 ,
^ 𝛼𝛽 := 𝐴𝛼 · 𝑀
𝑀
^𝛼𝛽
𝑁

32

^ · 𝐴𝛽 .
:= 𝐴𝛼 · 𝑁

(B.7)

Using Eq. (A.9), 𝑎𝑀^ and 𝑎𝑁^ can be written as
(︁
)︁2
(︁
)︁2
^𝛼𝛽 𝐸 (0)𝛼𝛽
^ 𝛼𝛽 𝐸 (0)𝛼𝛽 − 3 𝑁
,
𝑎𝑀^ := 3 𝑀
dev
dev
(︁
)︁ (︁
)︁
^ 𝛼𝛽 𝐸 (0)𝛼𝛽 𝑁
^𝛾𝛿 𝐸 (0)𝛾𝛿 .
𝑎 ^ := −6 𝑀
dev

𝑁

(B.8)

dev

Using Eqs. (B.6) to (B.8), the co-variant components of 𝑇 (0) can be obtained as
(0)

^ 𝛼𝛽 + 𝑎 ^ 𝑁
^
𝑇𝛼𝛽 = 𝑎𝑀^ 𝑀
𝑁 𝛼𝛽 .

(B.9)

[︀
]︀s
Substituting Eqs. (A.3) and (B.6) into Eq. (B.5) and using 𝑋 : [𝑌 ⊗ 𝑍]s = 𝑌 T 𝑋𝑍 T , see
𝜕𝒥3
Itskov (2015),
can be written as
𝜕𝐶
2
[︁
]︁s
𝜕𝒥3
𝜕𝒥3
1 ∑︁
1
(0)
=
𝐴𝛼 ⊗ 𝐴𝛽 =
𝑓𝑖𝑗 𝑃𝑖𝛼𝛾 𝑇𝛾𝛿 𝑃𝑗𝛿𝛽 𝐴𝛼 ⊗ 𝐴𝛽 =: 𝜇𝛼𝛽 𝐴𝛼 ⊗ 𝐴𝛽 ,
(B.10)
𝜕𝐶
𝜕𝑎𝛼𝛽
8
8
𝑖,𝑗=1

where 𝜇𝛼𝛽 is defined as
𝜇𝛼𝛽 :=

2
∑︁

(0)

𝑓𝑖𝑗 𝑃𝑖𝛼𝛾 𝑇𝛾𝛿 𝑃𝑗𝛿𝛽 .

(B.11)

𝑖,𝑗=1

Using the obtained relations, the Kirchhoff stress tensor can be written as
]︁
[︁
1
𝛼𝛽
= 𝜀𝛼
^ 2 𝜖a 𝑒−𝛼^ 𝜖a − 2 (𝜂1 𝜖a 𝒥3 + 𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽 + 2 𝜇 ln 𝜆 𝜒𝛼 𝛽 + 𝜂 𝜇𝛼𝛽 ,
𝜏m
4

(B.12)

where 𝜒𝛼 𝛽 is defined as
(︂

1 𝛼𝛽
1
𝑃
− 2 𝑃2𝛼𝛽
𝜆21 1
𝜆2
In addition, the elastic tensor can be obtained as
𝜒𝛼𝛽 :=

)︂
.

(B.13)

[︁
]︁
^
= 𝜀𝛼
^ 2 (1 − 𝛼
^ 𝜖a ) 𝑒−𝛼^ 𝜖a − 2(𝜂1 𝒥3 + 𝜇1 𝛽^2 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽 𝑎𝛾𝛿
𝑐𝛼𝛽𝛾𝛿
m
[︂
]︂
1
^
−
𝜂1 𝜖a 𝜇𝛾𝛿 + 2 𝜇1 𝛽^ 𝑒𝛽 𝜖a ln 𝜆 𝜒𝛾 𝛿 𝑎𝛼𝛽
2
[︀ 2
]︀
^
+ 2 𝜀𝛼
^ 𝜖a 𝑒−𝛼^ 𝜖a − 2 (𝜂1 𝜖a 𝒥3 + 𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽𝛾𝛿
^
− 2 𝜇1 𝛽^ 𝑒𝛽 𝜖a ln(𝜆) 𝜒𝛼 𝛽 𝑎𝛾𝛿 + 𝜇 𝜒𝛼 𝛽 𝜒𝛾 𝛿

(︃
)︃
)︂
1 𝛼𝛽 𝛾𝛿
1 𝛼𝛽 𝛾𝛿
1 𝜕𝑃1𝛼𝛽
1 𝜕𝑃2𝛼𝛽
𝑃 𝑃2 − 4 𝑃1 𝑃1
− 2
4 𝜇 ln 𝜆
+ 4 𝜇 ln 𝜆
𝜆42 2
𝜆1
𝜆21 𝜕𝑎𝛾𝛿
𝜆2 𝜕𝑎𝛾𝛽
[︁
]︁s
𝜕𝑓
1
1
𝑖𝑗
(0)
𝜂1 𝜖a 𝜇𝛼𝛽 𝑎𝛾𝛿 + 𝜂
𝑃𝑖𝛼𝜉 𝑇𝜉𝜐 𝑃𝑗𝜐𝛽
2
2 𝜕𝑎𝛾𝛿
[︃
]︃s
𝜕𝑃𝑖𝛼𝜉 (0) 𝜐𝛽
1
𝜂 𝑓𝑖𝑗
𝑇 𝑃
2
𝜕𝑎𝛾𝛿 𝜉𝜐 𝑗
[︃
]︃
𝜐𝛽 s
𝜕𝑃
1
𝑗
(0)
𝜂 𝑓𝑖𝑗 𝑃𝑖𝛼𝜉 𝑇𝜉𝜐
2
𝜕𝑎𝛾𝛿
⎡
⎤s
(0)
𝜕𝑇
1
𝜉𝜐
𝜂 𝑓𝑖𝑗 ⎣𝑃𝑖𝛼𝜉
𝑃 𝜐𝛽 ⎦ ,
2
𝜕𝑎𝛾𝛿 𝑗
(︂

+
−
+
+

+

33

(B.14)

where

𝜕𝑃1𝛼𝛽 𝜕𝑃2𝛼𝛽 𝜕𝑓𝑖𝑗
𝜕𝑎𝛾𝛿 , 𝜕𝑎𝛾𝛿 , 𝜕𝑎𝛾𝛿 ,

𝑎𝛼𝛽𝛾𝛿 and 𝐴𝛼𝛽𝛾𝛿 can be written as

𝜕𝑃1𝛼𝛽
𝐴𝛼𝛽𝛾𝛿 + 𝑃2𝛾𝛿 𝐴𝛼𝛽 + (𝑃1𝛾𝛿 − 𝑃2𝛾𝛿 ) 𝑃1𝛼𝛽
,
= −
𝜕𝑎𝛾𝛿
Λ1 − Λ2

(B.15)

𝐴𝛼𝛽𝛾𝛿 + 𝑃1𝛾𝛿 𝐴𝛼𝛽 + (𝑃2𝛾𝛿 − 𝑃1𝛾𝛿 ) 𝑃2𝛼𝛽
𝜕𝑃2𝛼𝛽
,
= −
𝜕𝑎𝛾𝛿
Λ2 − Λ1

(B.16)

⎧
1
⎪
− 4 𝑃𝑖𝛾𝛿 ;
⎪
⎪
⎪
2𝜆𝑖
⎨

𝜕𝑓𝑖𝑗
=
⎪
𝜕𝑎𝛾𝛿
⎪
⎪
⎪
⎩

𝜆2𝑖 −𝜆2𝑗
𝜆2𝑖

if 𝑖 = 𝑗 ,

− 2(ln 𝜆𝑖 − ln 𝜆𝑗 )
2(𝜆2𝑖 − 𝜆2𝑗 )2

𝑃𝑖𝛾𝛿 −

− 2(ln 𝜆𝑖 − ln 𝜆𝑗 )
2(𝜆2𝑖 − 𝜆2𝑗 )2

(B.17)
𝑃𝑗𝛾𝛿 ;

if 𝑖 ̸= 𝑗 ,

)︁
1 (︁ 𝛼𝛾 𝛽𝛿
𝑎 𝑎 + 𝑎𝛼𝛿 𝑎𝛽𝛾 ,
2

(B.18)

)︁
1 (︁ 𝛼𝛾 𝛽𝛿
𝐴 𝐴 + 𝐴𝛼𝛿 𝐴𝛽𝛾 .
2

(B.19)

𝑎𝛼𝛽𝛾𝛿 = −

𝐴𝛼𝛽𝛾𝛿 = −

𝜆2𝑖 −𝜆2𝑗
𝜆2𝑗

(0)

In addition,

𝜕𝑇𝜉𝜐
𝜕𝑎𝛾𝛿

can be obtained as
(0)

𝜕𝑇𝜉𝜐
𝜕C

(0)

=

𝜕𝐸 (0)

(0)

𝜕𝑇𝜉𝜐

𝜕𝑎𝛾𝛿

(0)

𝜕𝑇𝜉𝜐

𝜕𝑇𝜉𝜐
𝜕𝐸 (0)
:
=
𝐴𝛾 ⊗ 𝐴𝛿 ,
𝜕C
𝜕𝑎𝛾𝛿

2
∑︁

=

(B.20)

]︁s
[︁
𝑓𝑘𝑙 𝑃𝑘𝛾𝜌 𝐻𝜉𝜐𝜌𝜔 𝑃𝑙𝜔𝛿 ,

(B.21)

𝑘,𝑙=1
(0)

where 𝐻𝜉𝜐𝜌𝜔 can be related to the derivative of

𝜕𝑇𝜉𝜐

𝜕𝐸 (0)

as

(0)

𝜕𝑇𝜉𝜐
𝜕𝐸

𝐻𝜉𝜐𝜌𝜔

(0)

(B.22)

= 𝐻𝜉𝜐𝜌𝜔 𝐴𝜌 ⊗ 𝐴𝜔 ,

[︁

^ ^
^ ^
^ 𝜉𝜐 𝑀
^ 𝜌𝜔 + 𝑏 ^ ^ 𝑀
^ 𝜉𝜐 𝑁
^𝜌𝜔 + 𝑏 ^ ^ 𝑁
= 6 𝑏𝑀^ 𝑀^ 𝑀
^𝑁
^ 𝑁𝜉𝜐 𝑁𝜌𝜔
𝑀𝑁
𝑁 𝑀 𝜉𝜐 𝑀𝜌𝜔 + 𝑏𝑁

]︁

,

(B.23)

where 𝑏𝑀^ 𝑀^ , 𝑏𝑀^ 𝑁^ , 𝑏𝑁^ 𝑀^ and 𝑏𝑁^ 𝑁^ are defined as
^ 𝜉𝜐 𝐸 (0)𝜉𝜐 ,
𝑏𝑀^ 𝑀^ := 𝑀
dev
(0)𝜉𝜐

,

(0)𝜉𝜐

,

^𝜉𝜐 𝐸
𝑏𝑀^ 𝑁^ := −𝑁
dev

^𝜉𝜐 𝐸
𝑏𝑁^ 𝑀^ := −𝑁
dev

^ 𝜉𝜐 𝐸 (0)𝜉𝜐 .
𝑏𝑁^ 𝑁^ := −𝑀
dev

34

(B.24)

The elastic tensor can then be simplified as
[︁
]︁
^
𝑐𝛼𝛽𝛾𝛿
= 𝜀𝛼
^ 2 (1 − 𝛼
^ 𝜖a ) 𝑒−𝛼^ 𝜖a − 2(𝜂1 𝒥3 + 𝜇1 𝛽^2 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽 𝑎𝛾𝛿
m
−
+
+
+
−

[︀
]︀
[︀
]︀
^
𝜂1 𝜖a 𝑎𝛼𝛽 𝜇𝛾𝛿 + 𝜇𝛼𝛽 𝑎𝛾𝛿 − 2 𝜇1 𝛽^ 𝑒𝛽 𝜖a ln(𝜆) 𝜒𝛼 𝛽 𝑎𝛾𝛿 + 𝑎𝛼𝛽 𝜒𝛾 𝛿
]︁
[︁
^
2 𝜀𝛼
^ 2 𝜖a 𝑒−𝛼^ 𝜖a − 2 (𝜂1 𝜖a 𝒥3 + 𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽𝛾𝛿 + 𝜇 𝜒𝛼 𝛽 𝜒𝛾 𝛿
(︃
)︃
(︂
)︂
1 𝛼𝛽 𝛾𝛿
1 𝛼𝛽 𝛾𝛿
1 𝜕𝑃1𝛼𝛽
1 𝜕𝑃2𝛼𝛽
4 𝜇 ln 𝜆
𝑃 𝑃2 − 4 𝑃1 𝑃1
+ 4 𝜇 ln 𝜆
− 2
𝜆42 2
𝜆1
𝜆21 𝜕𝑎𝛾𝛿
𝜆2 𝜕𝑎𝛾𝛽
{︂
𝜕𝑓𝑖𝑗 𝛼 𝛽
𝑓𝑖𝑗 [︀ 𝛼 𝛽 𝛾𝛿
1
2𝜐𝑖𝑗 (𝑃𝑖 − 𝑃¯𝑖𝛾𝛿 ) + (𝜁𝑗𝛼 𝛽 + 𝜁𝑗𝛽 𝛼 ) 𝑃¯𝑖𝛾 𝛿
𝜂
𝜐𝑖𝑗 −
2
𝜕𝑎𝛾𝛿
Λ𝑖 − Λ¯𝑖
}︂
]︀
𝛿𝛽
𝛾𝛽
𝛾 𝛼 𝛽𝛿
1
𝛼𝛾
𝛼𝛿
𝛿
𝛼
𝛽𝛾
𝛼𝛽𝛾𝛿
¯
𝜁𝑗 + 𝐴 𝜁𝑗 + 𝜁𝑗 𝐴 + 𝜁𝑗 𝐴 ) + 𝐻
,
2 (𝐴
1
2

(B.25)

𝛼𝛽
¯ 𝛼𝛽𝛾𝛿 are defined as
where ¯𝑖 flips with 𝑖 (¯𝑖 is 2 when 𝑖 is 1 and vice versa) and 𝜐𝑖𝑗
, 𝜁𝑗𝛼 𝛽 and 𝐻
𝛼𝛽
𝜐𝑖𝑗
:= 𝑃𝑖𝛼𝜉 𝑇𝜉𝜐 𝑃𝑗𝜐 𝛽 ,

(B.26)

𝜁𝑗𝛼 𝛽 := 𝐴𝛼 𝜉 𝑇𝜉𝜐 𝑃𝑗𝜐 𝛽 ,

(B.27)

¯ 𝛼𝛽𝛾𝛿 := 6[𝑏 ^ ^ 𝜅𝛼𝛽 𝜅𝛾𝛿 + 𝑏 ^ ^ 𝜅𝛼𝛽 𝜆𝛾𝛿 + 𝑏 ^ ^ 𝜆𝛼𝛽 𝜅𝛾𝛿 + 𝑏 ^ ^ 𝜆𝛼𝛽 𝜆𝛾𝛿 ] .
𝐻
𝑀𝑀
𝑀𝑁
𝑁𝑀
𝑁𝑁

(B.28)

𝜅𝛼𝛽 and 𝜆𝛼 𝛽 are defined as
𝜅𝛼𝛽 :=

2
∑︁

^ 𝜉𝜐 𝑃 𝜐 𝛽 ,
𝑓𝑖𝑗 𝑃𝑖𝛼𝜉 𝑀
𝑗

(B.29)

^𝜉𝜐 𝑃 𝜐 𝛽 .
𝑓𝑖𝑗 𝑃𝑖𝛼𝜉 𝑁
𝑗

(B.30)

𝑖,𝑗=1

𝜆𝛼 𝛽 :=

2
∑︁
𝑖,𝑗=1

B.2

Repeated eigenvalues

In addition, the special case of repeated eigenvalues (Λ1 = Λ2 ) should be considered. In this
case, the right Cauchy-Green tensor can be written as

𝜕Λ1
𝜕𝐶

𝐶 = Λ1 𝐼 .

(B.31)

𝜕Λ1 𝜕𝑎𝛼𝛽
𝜕Λ1
𝜕Λ1
=
=
𝐴𝛼 ⊗ 𝐴𝛽 .
𝜕𝐶
𝜕𝑎𝛼𝛽 𝜕𝐶
𝜕𝑎𝛼𝛽

(B.32)

can be written by the chain rule as

In addition,

𝜕Λ1
𝜕𝐶

can be written as

So, the following relation for

𝜕Λ1
= 𝐼 .
𝜕𝐶
𝜕Λ1
𝜕𝑎𝛼𝛽

(B.33)

can be obtained (Simo and Taylor, 1991; Miehe, 1993)
𝜕Λ1
= 𝐴𝛼𝛽 .
𝜕𝑎𝛼𝛽
35

(B.34)

𝜕𝐸 (0)
𝜕𝐶

can be written as
1
𝜕𝐸 (0)
=
[𝐼 ⊗ 𝐼]𝑠 .
𝜕𝐶
2𝜆21

The obtained relation for
𝜕𝒥3
𝜕𝑎𝛼𝛽

𝜕𝒥1
𝜕𝑎𝛼𝛽

(B.35)

is the same for distinct and repeated eigenvalues. But,

should be derived for repeated eigenvalues.

𝜕𝒥2
𝜕𝐶

𝜕𝛾𝑖2
𝜕𝐸 (0)

𝜕𝒥2
𝜕𝐶

(B.36)

can be written as
𝜕𝛾𝑖2
𝜕𝐸

Next,

and

can be expanded as

𝜕𝒥2
𝜕𝒥2 𝜕𝐸 (0)
1 𝜕𝛾𝑖2
𝜕𝐸 (0)
=
:
=
:
,
𝜕𝐶
𝜕𝐶
4 𝜕𝐸 (0)
𝜕𝐶
𝜕𝐸 (0)
where

𝜕𝒥2
𝜕𝑎𝛼𝛽

and

𝜕𝒥2
𝜕𝑎𝛼𝛽

(0)

(0)

(B.37)

= 4𝐸 dev ,

are connected as
𝜕𝒥2
𝜕𝒥2
=
𝐴𝛼 ⊗ 𝐴𝛽 ,
𝜕𝐶
𝜕𝑎𝛼𝛽

(B.38)

Substituting Eqs. (B.37) and (B.35) into Eq. (B.36) and using Eq. (B.38) results in
1
𝜕𝒥2
(0)𝛼𝛽
=
𝐸
= 𝜂 𝛼𝛽 .
𝜕𝑎𝛼𝛽
2𝜆21 dev
The same procedure should be repeated to obtain

𝜕𝒥3
𝜕𝒥3
𝜕𝑎𝛼𝛽 . 𝜕𝐶

(B.39)
can be written as

𝜕𝒥3
𝜕𝒥3 𝜕𝐸 (0)
=
:
.
𝜕𝐶
𝜕𝐶
𝜕𝐸 (0)
Substituting Eqs. (B.6) and (B.35) into Eq. (B.40) results in
]︁𝑠
]︁𝑠
𝜕𝒥3
1 [︁
1 [︁ (0)𝛼𝛽
(0)
=
𝐼
𝑇
𝐼
𝑇
𝐴
⊗
𝐴
,
=
𝛼
𝛽
𝜕𝐶
16𝜆21
16𝜆21

(B.40)

(B.41)

where 𝑇 (0)𝛼𝛽 is defined as
(0)

𝑇 (0)𝛼𝛽 := 𝐴𝛼𝛾 𝑇𝛾𝛿 𝐴𝛽𝛿 .
Furthermore,

So,

𝜕𝒥3
𝜕𝑎𝛼𝛽

𝜕𝒥3
𝜕𝐶

and

𝜕𝒥3
𝜕𝑎𝛼𝛽

(B.42)

are connected as
𝜕𝒥3
𝜕𝒥3
=
𝐴𝛼 ⊗ 𝐴𝛽 .
𝜕𝐶
𝜕𝑎𝛼𝛽

(B.43)

𝜕𝒥3
1
1
𝑇 (0)𝛼𝛽 = 𝜇𝛼𝛽 .
=
2
𝜕𝑎𝛼𝛽
8
16𝜆1

(B.44)

can be obtained as

Using the obtained relations, the stress tensor can be written as
[︁
]︁
1
^
𝛼𝛽
𝜏m
= 𝜀𝛼
^ 2 𝜖a 𝑒−𝛼^ 𝜖a − 2 (𝜂1 𝜖a 𝒥3 + 𝜇1 𝛽^ 𝒥2 𝑒𝛽 𝜖a ) 𝑎𝛼𝛽 + 4 𝜇 𝜂 𝛼𝛽 + 𝜂 𝜇𝛼𝛽 .
4

(B.45)

It should be noted that all 𝒥3 , 𝒥2 , 𝜂 𝛼𝛽 and 𝜇𝛼𝛽 are zero in the case of repeated eigenvalues.
But, they should be kept to take the second derivative and obtain the elastic tensor. It is easy
𝛼𝛽
𝜕𝜂 𝛼𝛽
𝜕𝜂 𝛼𝛽
to show 𝜕𝜇
𝜕𝑎𝛾𝛿 is zero, so only 𝜕𝑎𝛾𝛿 plays a role in the elastic tensor. 𝜕𝑎𝛾𝛿 can be written as
36

(0) 𝛼𝛽

1 𝜕𝐸dev
1
𝜕𝜂 𝛼𝛽
(0) 𝛼𝛽
= − 2 𝐴𝛾𝛿 𝐸dev +
𝜕𝑎𝛾𝛿
2Λ1 𝜕𝑎𝛾𝛿
2Λ1

.

(B.46)

The first term on the right hand side of Eq. (B.46) is zero. The second term can be obtained
from
)︂
(︂
(0)
(0)
𝜕𝐸 dev
𝜕𝐸 dev 𝜕𝐸 (0)
1
𝐼 ⊙𝐼
𝑠
(B.47)
:
[𝐼 ⊗ 𝐼]𝑠 .
=
:
=
[𝐼
⊗
𝐼]
−
(0)
𝜕𝐶
𝜕𝐶
2
2Λ
𝜕𝐸
1
Eq. (B.47) can be further simplified by considering the following relations,
𝐼 ⊗𝐼 :𝑋
𝐼 ⊗𝐼 :𝒜
(𝐼 ⊗ 𝐼)t : 𝑋
𝒜 : (𝐼 ⊗ 𝐼)t
(𝐼 ⊗ 𝐼)t : 𝒜
(𝐼 ⊗ 𝐼)t : (𝐼 ⊗ 𝐼)𝑡

=
=
=
=
=
=

𝑋 :𝐼 ⊗𝐼 =𝑋 ,
𝒜:𝐼 ⊗𝐼 =𝒜 ,
𝑋 : (𝐼 ⊗ 𝐼)t = 𝑋 T ,
𝒜t ,
𝒜TtT ,
𝐼 ⊗𝐼 ,

(B.48)

where 𝑋 and 𝒜 are second and fourth order tensors. The final simplified form of Eq. (B.47) is
(0)

𝜕𝐸 dev
1
=
𝜕𝐶
2Λ1
𝜕𝐸 dev
𝜕𝐶

(B.49)

(0) 𝛼𝛽

(0)

Next,

(︂
)︂
𝐼 ⊙𝐼
[𝐼 ⊗ 𝐼]𝑠 −
.
2

and

𝜕𝐸dev
𝜕𝑎𝛾𝛿

are connected as
(0) 𝛼𝛽

𝜕𝐸dev
𝜕𝐶
(0) 𝛼𝛽

Finally, the basis vectors of

𝜕𝐸dev
𝜕𝑎𝛾𝛿

𝜕𝐸dev
𝜕𝑎𝛾𝛿

(𝐴𝛾 ⊗ 𝐴𝛿 ) .

(B.50)

𝛼𝛽
, 𝜏m
and 𝑐𝛼𝛽𝛾𝛿
can be written as
m

(0) 𝛼𝛽

𝜕𝐸dev
𝜕𝑎𝛾𝛿

𝜕𝐸dev
𝜕𝑎𝛾𝛿

should be rearranged as 𝐴𝛼 ⊗ 𝐴𝛽 ⊙ 𝐴𝛾 ⊗ 𝐴𝛿 to be used

(0) 𝛼𝛽

in the FEM implementation.

(0) 𝛼𝛽

=

= −

]︁
1 [︁ 𝛼𝛽𝛾𝛿
2𝐴
+ 𝐴𝛼𝛽 𝐴𝛾𝛿 ,
4Λ1

[︀ 2
]︀
𝛼𝛽
𝜏m
= 𝜀𝛼
^ 𝜖a 𝑒−𝛼^ 𝜖a 𝑎𝛼𝛽 ,

(B.51)

(B.52)

]︁
[︁
]︁
[︁
]︁
𝜇 [︁ 𝛼𝛽𝛾𝛿
𝛼𝛽 𝛾𝛿
2
−𝛼
^ 𝜖a
𝛼𝛽 𝛾𝛿
2
−𝛼
^ 𝜖a
𝛼𝛽𝛾𝛿
𝑐𝛼𝛽𝛾𝛿
=
𝜀
𝛼
^
(1
−
𝛼
^
𝜖
)
𝑒
𝑎
𝑎
+
2
𝜀
𝛼
^
𝜖
𝑒
2𝐴
+
𝐴
𝐴
.
𝑎
−
a
a
m
Λ21
(B.53)

C

Analytical solution for arbitrary stretch in Cartesian coordinates

In this section, a analytical expression of the Cauchy stress tensor in a Cartesian coordinates
system (Fig. C.1) is obtained for the case of arbitrary stretch along 𝑒1 and 𝑒2 directions. Then,
a relation for pure shear is extracted from the solution of arbitrary stretch as a special case.
This solution is used in the verification of the FE implementation in Sec. 5. It should be noted
37

e3
e2

e1

Figure C.1: Local Cartesian coordinate system
that the general expression for the Cartesian stress components can be directly obtained from
Eq. (40) as
𝜎𝑖𝑗 = 𝑒𝑖 · 𝜎 · 𝑒𝑗 .
(C.1)
This expression fully defines the stress state of the shell. In the following, the particular expressions for 𝜎11 , 𝜎22 and 𝜎12 are derived for the case of a homogeneous membrane stretch, i.e.
the deformation gradient is given by
𝐹 = 𝜆1 𝑒1 ⊗ 𝑒1 + 𝜆2 𝑒2 ⊗ 𝑒2 ;

𝜆1 > 𝜆2 ,

(C.2)

The right Cauchy-Green tensor, right stretch tensor and logarithmic volumetric and deviatoric
strains can be written as
𝐶 = Λ1 𝑒1 ⊗ 𝑒1 + Λ2 𝑒2 ⊗ 𝑒2 ,
(C.3)
𝑈 = 𝜆1 𝑒1 ⊗ 𝑒1 + 𝜆2 𝑒2 ⊗ 𝑒2 ,
𝜖a =

1
ln(𝐽) ,
2

(C.4)
(C.5)

(0)

𝐸 dev = ln 𝜆 ( 𝑒1 ⊗ 𝑒1 − 𝑒2 ⊗ 𝑒2 ) .

(C.6)

The logarithmic stress can be written as
[︁
]︁
𝜂(𝜖a )
(0)
2
−𝛼
^ 𝜖a
′
′
=
𝜀
𝛼
^
𝜖
𝑒
+
2
𝜇
(𝜖
)
𝒥
+
𝜂
(𝜖
)𝒥
𝑆 𝐸 (0) ,
𝐼 + 2𝜇 (𝜖a ) 𝐸 dev +
a
a
2
a
3
(0)
8
dev
𝜕𝐸
(C.7)
where 𝑆 𝐸 (0) is defined as
𝑆 (0) :=

𝜕𝑊

dev

𝑆 𝐸 (0) := 3
dev

{︂[︂(︁

^ : 𝐸 (0)
𝑀
dev

)︁2

(︁

^ : 𝐸 (0)
− 𝑁
dev

)︁2 ]︂

𝑀 −2

[︁(︁

^ : 𝐸 (0)
𝑀
dev

)︁ (︁
)︁]︁ }︂
(0)
^
𝑁 : 𝐸 dev 𝑁 .
(C.8)

^ and 𝑁
^ can be written as
𝑀
^ = 𝑥
^ ⊗𝑥
^ −𝑦
^⊗𝑦
^ = cos(2𝜙) 𝑒1 ⊗ 𝑒1 + sin(2𝜙) (𝑒1 ⊗ 𝑒2 + 𝑒2 ⊗ 𝑒1 ) − cos(2𝜙) 𝑒2 ⊗ 𝑒2 ,
𝑀
(C.9)
^
^ ⊗𝑦
^+𝑦
^⊗𝑥
^ = − sin(2𝜙) 𝑒1 ⊗ 𝑒1 + cos(2𝜙) (𝑒1 ⊗ 𝑒2 + 𝑒2 ⊗ 𝑒1 ) + sin(2𝜙) 𝑒2 ⊗ 𝑒2 ,
𝑁 = 𝑥
(C.10)
^
^
where 𝑥 and 𝑦 are
^ = cos(𝜙) 𝑒1 + sin(𝜙) 𝑒2 ,
𝑥
(C.11)
^ = − sin(𝜙) 𝑒1 + cos(𝜙) 𝑒2 .
𝑦

38

^ . Using the following relations
𝜙 is the angle between 𝑒1 and 𝑥
^ : 𝐸 (0) = cos(2𝜙) ln(𝜆2 ) ,
𝑀
dev
^ : 𝐸 (0) = − sin(2𝜙) ln(𝜆2 ) ,
𝑁
dev

(C.12)

(︁
)︁
^ + sin(4𝜙) 𝒥2 𝑁
^ .
𝑆 𝐸 (0) = 12 cos(4𝜙) 𝒥2 𝑀

(C.13)

Eq. (C.8) can be simplified as

dev

In addition, the logarithmic stress can be simplified as
(0)

𝑆 (0) = 𝑉1 𝐼 + 𝑉2 𝐸 dev + 𝑉3 𝑀 + 𝑉4 𝑁 ,

(C.14)

where 𝑉𝑖 are defined as
𝑉1 := 𝜀 𝛼
^ 2 𝜖a 𝑒−𝛼^ 𝜖a + 2 𝜇′ (𝜖a ) 𝒥2 + 𝜂 ′ (𝜖a ) 𝒥3 ,
𝑉2 := 2𝜇 (𝜖a ) ,
(C.15)

3
𝜂(𝜖a ) cos(4𝜙) 𝒥2 ,
2
3
𝑉4 := 𝜂(𝜖a ) sin(4𝜙) 𝒥2 .
2
𝑉3 :=

𝑆 (0) can be converted to the second Piola-Kirchhoff stress tensor 𝑆 (2) as
𝑆 (2) = 2

𝜕𝐸 (0)
: 𝑆 (0) .
𝜕𝐶

(C.16)

The eigenprojections for the assumed gradient deformation can be written as
𝑃 1 = 𝑒1 ⊗ 𝑒1 ,
𝑃 2 = 𝑒2 ⊗ 𝑒2 .

(C.17)

Considering the following relations,
𝑃1 ⊗ 𝑃1
𝑃1 ⊗ 𝑃2
𝑃1 ⊗ 𝑃1
𝑃1 ⊗ 𝑃2
𝑃1 ⊗ 𝑃1
𝑃1 ⊗ 𝑃2
𝑃1 ⊗ 𝑃1
𝑃1 ⊗ 𝑃2

:𝐼
:𝐼
(0)
: 𝐸 dev
(0)
: 𝐸 dev
^
:𝑀
^
:𝑀
^
:𝑁
^
:𝑁

=
=
=
=
=
=
=
=

𝑃1
0
ln(𝜆) 𝑃 1
0
cos(2𝜙) 𝑃 1
sin(2𝜙) (𝑒1 ⊗ 𝑒2 )
− sin(2𝜙) 𝑃 1
cos(2𝜙) (𝑒1 ⊗ 𝑒2 )

;
;
;
;
;
;
;
;

𝑃2 ⊗ 𝑃2
𝑃2 ⊗ 𝑃1
𝑃2 ⊗ 𝑃2
𝑃2 ⊗ 𝑃1
𝑃2 ⊗ 𝑃2
𝑃2 ⊗ 𝑃1
𝑃2 ⊗ 𝑃2
𝑃2 ⊗ 𝑃1

:𝐼
:𝐼
(0)
: 𝐸 dev
(0)
: 𝐸 dev
^
:𝑀
^
:𝑀
^
:𝑁
^
:𝑁

=
=
=
=
=
=
=
=

𝑃2 ,
0,
− ln(𝜆) 𝑃 2 ,
0,
− cos(2𝜙) 𝑃 2 ,
sin(2𝜙) (𝑒2 ⊗ 𝑒1 ) ,
sin(2𝜙) 𝑃 2 ,
cos(2𝜙) (𝑒2 ⊗ 𝑒1 ) ,

Eq. (C.16) can be simplified as
[︂
]︂
1
1
(2)
𝑆
= 𝑉1
𝑃1 +
𝑃2
Λ1
Λ2
[︂
]︂
ln(𝜆)
ln(𝜆)
+ 𝑉2
𝑃1 −
𝑃2
Λ1
Λ2
[︂
]︂
cos(2𝜙)
cos(2𝜙)
ln(𝜆4 )
+ 𝑉3
𝑃1 −
𝑃2 −
sin(2𝜙) (𝑒1 ⊗ 𝑒2 + 𝑒2 ⊗ 𝑒1 )
Λ1
Λ2
Λ2 − Λ1
[︂
]︂
sin(2𝜙)
ln(𝜆4 )
sin(2𝜙)
+ 𝑉4 −
𝑃1 +
𝑃2 −
cos(2𝜙) (𝑒1 ⊗ 𝑒2 + 𝑒2 ⊗ 𝑒1 ) .
Λ1
Λ2
Λ2 − Λ1
39

(C.18)

(C.19)

𝑆 (2) can be converted to 𝜎 by using the following relation,
𝜎 =

1
𝐹 𝑆 (2) 𝐹 T .
𝐽

(C.20)

Λ1 𝑃 1 ,
Λ2 𝑃 2 ,
𝜆1 𝜆2 𝑒1 ⊗ 𝑒2 = 𝐽 𝑒1 ⊗ 𝑒2 ,
𝜆1 𝜆2 𝑒2 ⊗ 𝑒1 = 𝐽 𝑒2 ⊗ 𝑒1 ,

(C.21)

Considering the following relations,
𝐹
𝐹
𝐹
𝐹

𝑃1 𝐹T
𝑃2 𝐹T
(𝑒1 ⊗ 𝑒2 ) 𝐹 T
(𝑒2 ⊗ 𝑒1 ) 𝐹 T

=
=
=
=

the Cartesian components of 𝜎 then become
1
[𝑉1 + 𝑉2 ln(𝜆) + 𝑉3 cos(2𝜙) − 𝑉4 sin(2𝜙)] ,
𝐽
1
= [𝑉1 − 𝑉2 ln(𝜆) − 𝑉3 cos(2𝜙) + 𝑉4 sin(2𝜙)] ,
𝐽
ln(𝜆4 )
[𝑉3 sin(2𝜙) + 𝑉4 cos(2𝜙)] .
= −
Λ2 − Λ1

𝜎11 =
𝜎22
𝜎12

(C.22)

The analytical solution for pure shear can be easily obtained from Eq. (C.22) by assuming
𝜆2 = 𝜆11 .

D

Coarse grained contact model (CGCM)

This section focuses on the modeling of the substrate. The CGCM is introduced and based
on that, the equivalent contact force and stiffness matrix are obtained for the substrate. Next,
the substrate geometry is formulated based on analytical geometries. The closest point project
(CPP) and principal curvatures are obtained, which are used in the contact formulation.

D.1

Atomic interaction with a half space

The CGCM is used to avoid fully atomistic simulation of interactions between substrate and
graphene. The type of interaction is van-der-Waals (vdW) interaction. It is assumed that
the radius of curvature of any surface point of the substrate is much higher than the cut-off
radius of the used potential. So, a half space approximation is a good estimation. In the
following, an analytical solution for the interaction between graphene and the half space is
given. The Lennard-Jones (L-J) potential is considered to model the atomistic interaction. The
L-J potential between a pair of atoms 𝑖 and 𝑗 can be written as
𝐶1
𝐶2
(D.1)
+ 12 ,
6
𝑟
𝑟
where 𝐶1 and 𝐶2 are constants and 𝑟 is the distance between two atoms. The half space potential
can be written as (Israelachvili, 2011; Sauer and Li, 2007b,a; Aitken and Huang, 2010)
[︃ (︂ )︂
(︂ )︂ ]︃
3 ℎ0 3 1 ℎ0 9
Ψ𝑖𝑗(vdW)h = −Γ
−
,
(D.2)
2 𝑟
2 𝑟
𝜓𝑖𝑗(vdW) = −

where ℎ0 , Γ and 𝑟 are the equilibrium distance, interfacial adhesion energy per unit area and
normal distance of a point of the graphene sheet on the substrate. They are defined as
(︂
)︂ 1
2𝐶2 6
(D.3)
ℎ0 :=
,
5𝐶1
40

Γ :=

𝜋 𝜌s 𝜌g 𝐶1
,
9ℎ30

(D.4)

where 𝜌s and 𝜌g are the number of atoms per unit volume and area of the substrate and
graphene, respectively. Hence, 𝐶1 and 𝐶2 can be written based on ℎ0 and Γ as
𝐶1 =

9ℎ30 Γ
,
𝜋 𝜌s 𝜌g

(D.5)

𝐶2 =

45ℎ90 Γ
.
2𝜋 𝜌s 𝜌g

(D.6)

The various Γ and ℎ0 are provided in the literature (Hossain (2009), Gao et al. (2014), Aitken
and Huang (2010)). It is usual to estimate ℎ0 as the graphite interlayer spacing, which is 0.34
nm.

D.2

Van der Waals contact force and its corresponding contact stiffness matrix

The required contact force and corresponding stiffness matrix for the FE formulation are determined in this section. The contact force can be obtained by taking the derivative from the half
space potential. The contact stiffness matrix can be obtained by the linearization of the contact
force. The details of the linearization of the contact force as well as the contact algorithm can
be found in the literature (Sauer and De Lorenzis, 2013, 2015; Sauer and Wriggers, 2009). Here
the required formulation is shortly reviewed. The formulation is simplified by assuming one of
the bodies to be rigid. Based on this assumption, the contact force per unit current area of
graphene can be written as
𝑓g = −

𝜕Ψ𝑖𝑗(vdW)h
𝜕Ψ𝑖𝑗(vdW)h 𝜕𝑟
=−
= 𝑓g 𝑛p ,
𝜕𝑥g
𝜕𝑟
𝜕𝑥g

(D.7)

where 𝑥g is a point on the graphene surface, and 𝑟, 𝑟, 𝑛p and 𝑓g are defined as
𝑟

:= 𝑥g − 𝑥p ,

𝑟

:= ‖𝑟‖ ,
𝑟
,
𝑟
𝜕Ψ𝑖𝑗(vdW)h
:= −
.
𝜕𝑟

𝑛p :=
𝑓g

(D.8)

𝑥p is the CPP of 𝑥g onto the substrate. The corresponding contact potential can be written as
∫︁
∫︁
𝐺c =
𝛿𝜙g · 𝑓 g d𝑎 =
𝛿𝜙g · 𝑓 0g d𝐴,
(D.9)
𝜕ℬc

𝜕ℬ0c

where 𝛿𝜙g is an admissible variation for the graphene displacement. 𝜕ℬc and 𝜕ℬ0c are the contact boundaries between graphene and the substrate in the current and reference configurations,
respectively. 𝑓 0g and 𝑓 g can be connected by 𝐽 as
𝑓 0g = 𝐽 𝑓 g .

41

(D.10)

𝐺𝑒c can be obtained by changing the integral domain to an element. The nodal contact force
can be obtained as
∫︁
∫︁
T
𝑒
NT 𝑓 0g d𝐴 .
N 𝑓 g d𝑎 = −
fc = −
(D.11)
𝑒
𝜕ℬ0c

𝜕ℬc𝑒

The substrate is assumed to be rigid and 𝜌s is constant. Furthermore, 𝜌0g = 𝐽𝜌g is the density
in the reference configuration and can be considered constant in the linearization. Based on
these assumptions, the contact stiffness matrix can be written as
∫︁
𝜕𝑓 0g
𝑒
kc = −
NT
d𝐴 .
(D.12)
𝜕𝑢𝑒
𝑒
𝜕ℬ0𝑐

𝜕𝑓 0g
can be written as
𝜕𝑢𝑒

where

𝜕𝑓 0g
𝜕𝑥g

𝜕𝑓 0g 𝜕𝑥g
𝜕𝑓 0g
𝜕𝑓 0g
=
=
N,
𝜕𝑢𝑒
𝜕𝑥g 𝜕𝑢𝑒
𝜕𝑥g

(D.13)

(︀ )︀′
𝑓g0
𝜕𝑓 0g
𝑓g0
^ p1 + −1
^ p2 .
^ p1 ⊗ 𝑎
^p ⊗ 𝑎
𝑎
𝑎
= 𝑓g0 𝑛p ⊗ 𝑛p + −1
𝜕𝑥g
𝜅
¯1 + 𝑟
𝜅
¯2 + 𝑟 2

(D.14)

can be written as

^ p𝛼 are the normalized tangent vectors along the principal curvatures of the surface of the
𝑎
substrate. (𝑓g0 )′ is defined as
(︀ 0 )︀′
𝜕𝑓g0
(D.15)
𝑓g :=
.
𝜕𝑟
Convex geometries are assumed to have a positive curvature in the computation of the contact
stiffness. 𝜅
¯ 𝛼 are connected to the principal curvatures of the surface of the substrate 𝜅𝛼 by
𝜅
¯ 𝛼 = −𝜅𝛼 .

D.3

(D.16)

Modeling of the substrate geometry

In this section, an analytical description for the substrate geometry is presented. The substrate
is divided into three area which are a flat, cylindrical and torus part. In each area, the normal
unit vector to the surface, the principal curvatures and directions of the curvature tensor are
given. These quantities are used in computation of the contact force and the contact stiffness
matrix.
D.3.1

Torus part

The torus parametric description can be written as
𝑟(𝜙, 𝜔) = [𝑅 + 𝑟 cos(𝜔)] cos(𝜙) 𝑒1 + [𝑅 + 𝑟 cos(𝜔)] sin(𝜙) 𝑒2 + 𝑅 sin(𝜔) 𝑒3 ,

(D.17)

where 𝑅 and 𝑟 are major and minor radii of the torus. 𝜙 and 𝜔 are parametric coordinates.
The tangent vectors then follow as
𝑎1 = − [𝑅 + 𝑟 cos(𝜔)] sin(𝜙) 𝑒1 + [𝑅 + 𝑟 cos(𝜔)] cos(𝜙) 𝑒2 ,
𝑎2 = −𝑅 cos(𝜙) sin(𝜔) 𝑒1 − 𝑅 sin(𝜙) sin(𝜔) 𝑒2 + 𝑅 cos(𝜔) 𝑒3 .
42

(D.18)

The co-variant components of curvature tensor become
𝑏11 = − [𝑅 + 𝑟 cos(𝜔)] cos(𝜔) ,
𝑏22 = −𝑅 ,
𝑏12 = 𝑏21 = 0 .

(D.19)

The principal curvatures of the torus can be written as
𝑏11
cos(𝜔)
=−
,
𝑎11
𝑅 + 𝑟 cos(𝜔)
1
𝑏22
=− ,
𝜅2 = 𝑏22 =
𝑎22
𝑟
𝜅1 = 𝑏11 =

(D.20)

and 𝜅
¯ 𝛼 can be written as
cos(𝜔)
,
𝑅 + 𝑟 cos(𝜔)
(D.21)
1
𝜅
¯2 = .
𝑟
It should be denoted that the principal directions of the curvature tensor and the tangent vectors
are identical in the torus. The CPP of a point (𝑥0 , 𝑦0 , 𝑧0 ) onto the surface of the torus can be
written as (Li et al., 2016)
𝑥T
(D.22)
p = [𝑥0 𝑤0 , 𝑦0 𝑤0 , 𝑧0 𝑤1 ] ,
𝜅
¯1 =

where 𝑤0 and 𝑤1 are defined as
𝑟
𝑤1 := √︁ √︀
,
( 𝑥20 + 𝑦02 − 𝑅)2 + 𝑧02
(︃
)︃
𝑅
𝑅
𝑤0 := √︀ 2
+ 1 − √︀ 2
𝑤1 .
𝑥0 + 𝑦02
𝑥0 + 𝑦02

(D.23)

The normal unit vector to the surface of the torus can be written as
𝑛T
p = [cos(𝜙) cos(𝜔), sin(𝜙) cos(𝜔), sin(𝜔)] .

(D.24)

A point on the axis and major circle of the torus has infinite number of projections. They are
not considered in the CPP computation.
D.3.2

Cylinder and flat parts

The normal unit vector to the internal area of a cylinder is
𝑛T
p = −[cos(𝜃), sin(𝜃), 0] ,
where 𝜃 is defined as

(︂
𝜃 := arctan

𝑦0
𝑥0

(D.25)

)︂
.

(D.26)

The principal curvatures 𝜅𝛼 of the internal surface of the cylinder are zero and 1/𝑅. The
principal curvatures of a flat area are zero and the computation of its normal unit vector is
trivial.

43

References
Abaqus (2016). Abaqus website. Available at www.3ds.com/products-services/simulia/
products/abaqus/.
Aitken, Z. H. and Huang, R. (2010). Effects of mismatch strain and substrate surface corrugation on morphology of supported monolayer graphene. J. Appl. Phys., 107(12). DOI:
10.1063/1.3437642.
Ansari, R., Ajori, S., and Motevalli, B. (2012). Mechanical properties of defective single-layered
graphene sheets via molecular dynamics simulation. Superlattices Microstruct., 51(2):274–
289.
Arroyo, M. and Belytschko, T. (2002). An atomistic-based finite deformation membrane for
single layer crystalline films. J. Mech. Phys. Solids, 50(9):1941 – 1977.
Arroyo, M. and Belytschko, T. (2004a). Finite crystal elasticity of carbon nanotubes based
on the exponential cauchy-born rule. Phys. Rev. B, 69:115415. DOI: 10.1103/PhysRevB.69.115415.
Arroyo, M. and 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. and Belytschko, T. (2005). Continuum mechanics modeling and simulation of carbon
nanotubes. Meccanica, 40(4):455–469.
Balandin, A. A. (2011). Thermal properties of graphene and nanostructured carbon materials.
Nat. Mater., 10(8):569–581.
Bellido, E. P. and Seminario, J. M. (2010). Molecular dynamics simulations of folding of supported graphene. J. Phys. Chem. C, 114(51):22472–22477.
Boddeti, N. G., Liu, X., Long, R., Xiao, J., Bunch, J. S., and Dunn, M. L. (2013). Graphene
blisters with switchable shapes controlled by pressure and adhesion. Nano Lett., 13(12):6216–
6221. PMID: 24224793.
Boddeti, N. G., Long, R., and Dunn, M. L. (2016). Adhesion mechanics of graphene on textured
substrates. Int. J. Solids Struct., 97:56–74.
Brenner, D. W. (1990). Empirical potential for hydrocarbons for use in simulating the chemical
vapor deposition of diamond films. Phys. Rev. B, 42:9458–9471.
Brenner, D. W., Shenderova, O. A., Harrison, J. A., Stuart, S. J., Ni, B., and Sinnott, S. B.
(2002). A second-generation reactive empirical bond order (rebo) potential energy expression
for hydrocarbons. J. Phys.: Condens. Matter, 14(4):783. DOI: 10.1088/0953-8984/14/4/312.
Bunch, J. and Dunn, M. (2012). Adhesion mechanics of graphene membranes. Solid State
Commun., 152(15):1359 – 1364. Exploring Graphene, Recent Research Advances.
Canham, P. (1970). The minimum energy of bending as a possible explanation of the biconcave
shape of the human red blood cell. J. Theor. Biol., 26(1):61 – 81.
Cao, G. (2014). Atomistic studies of mechanical properties of graphene. Polymers, 6(9):2404.
Cauchy, A.-L. (1851). Note sur l‘équilibre et les mouvements vibratoires des corps solides.
Comptes-Rendus, 32:323–326.
44

Cottrell, J., Hughes, T., and Bazilevs, Y. (2009). Isogeometric Analysis: Toward Integration of
CAD and FEA. Wiley.
Delfani, M. and Shodja, H. (2015). An exact analysis for the hoop elasticity and pressureinduced twist of cnt-nanovessels and cnt-nanopipes. Mech. Mater., 82:47 – 62.
Delfani, M. and Shodja, H. (2016). A large-deformation thin plate theory with application to
one-atom-thick layers. J. Mech. Phys. Solids, 87:65 – 85.
Delfani, M., Shodja, H., and Ojaghnezhad, F. (2013). Mechanics and morphology of singlewalled carbon nanotubes: from graphene to the elastica. Philos. Mag., 93(17):2057–2088.
Delfani, M. R. and Shodja, H. M. (2013). An enhanced continuum modeling of the ideal strength
and the angle of twist in tensile behavior of single-walled carbon nanotubes. J. Appl. Phys.,
114(5):053521.
Duong, T. X., Roohbakhshan, F., and Sauer, R. A. (2017). A new rotation-free isogeometric thin
shell formulation and a corresponding continuity constraint for patch boundaries. Comput.
Methods in Appl. Mech. Eng., 316(Supplement C):43 – 83. Special Issue on Isogeometric
Analysis: Progress and Challenges.
Ericksen, J. L. (1979). On the symmetry of deformable crystals. Arch. Ration. Mech. Anal.,
72(1):1–13.
Fadda, G. and Zanzotto, G. (2000). The arithmetic symmetry of monoatomic 2-nets. Acta
Crystallographica Section A, 56(1):36–48.
Favata, A., Micheletti, A., Podio-Guidugli, P., and Pugno, N. M. (2016). Geometry and selfstress of single-wall carbon nanotubes and graphene via a discrete model based on a 2ndgeneration rebo potential. J. Elast., 125(1):1–37.
Gao, W., Xiao, P., Henkelman, G., Liechti, K. M., and Huang, R. (2014). Interfacial adhesion between graphene and silicon dioxide by density functional theory with van der Waals
corrections. J. Phys. D: Appl. Phys., 47(25):255301. DOI: 10.1088/0022-3727/47/25/255301.
Gao, X., Li, C., Song, Y., and Chou, T.-W. (2016). A continuum mechanics model of multibuckling in graphene - substrate systems with randomly distributed debonding. Int. J. Solids
Struct., 97–98:510 – 519.
Ghaffari, R., Alipour, K., Solgi, S., Irani, S., and Haddadpour, H. (2015). Investigation of
surface stress effect in 3D complex nano parts using FEM and modified boundary cauchyborn method. J. Comput. Sci., 10:1 – 12.
Ghiba, I.-D., Neff, P., and Šilhavỳ, M. (2015). The exponentiated Hencky-logarithmic strain
energy. improvement of planar polyconvexity. Int. J. Non Linear Mech., 71:48–51.
Gupta, S. S. and Batra, R. C. (2010). Elastic Properties and Frequencies of Free Vibrations of
Single-Layer Graphene Sheets. J. Comput. Theor. Nanosci., 7(10):2151–2164.
Han, J., Pugno, N. M., and Ryu, S. (2015). Nanoindentation cannot accurately predict the
tensile strength of graphene or other 2D materials. Nanoscale, 7:15672–15679.
He, Y., Chen, W. F., Yu, W. B., Ouyang, G., and Yang, G. W. (2013). Anomalous interface
adhesion of graphene membranes. Sci. Rep., 3:2660 EP –. DOI: 10.1038/srep02660.
Hossain, M. Z. (2009). Chemistry at the graphene-SiO2 interface. Appl. Phys. Lett., 95(14).
DOI: /10.1063/1.3247964.
45

Ishigami, M., Chen, J. H., Cullen, W. G., Fuhrer, M. S., and Williams, E. D. (2007). Atomic
structure of graphene on SiO2 . Nano Lett., 7(6):1643–1648.
Israelachvili, J. (2011). Intermolecular and Surface Forces: Revised Third Edition. Intermolecular and Surface Forces. Elsevier Science.
Itskov, M. (2015). Tensor Algebra and Tensor Analysis for Engineers: With Applications to
Continuum Mechanics. Mathematical Engineering. Springer International Publishing.
Javili, A., McBride, A., Steinmann, P., and Reddy, B. D. (2014). A unified computational
framework for bulk and surface elasticity theory: a curvilinear-coordinate-based finite element
methodology. Comput. Mech., 54(3):745–762.
Jomehzadeh, E. and Pugno, N. (2015). Bending stiffening of graphene and other 2d materials
via controlled rippling. Compos. Part B: Eng., 83:194 – 202.
Koenig, S. P., Boddeti, N. G., Dunn, M. L., and Bunch, J. S. (2011). Ultrastrong adhesion of
graphene membranes. Nat. Nanotechnol., 6(9):543–546.
Kudin, K. N., Scuseria, G. E., and Yakobson, B. I. (2001). C2 F, BN, and C nanoshell elasticity
from ab initio computations. Phys. Rev. B, 64:235406.
Kumar, S. and Parks, D. M. (2014). On the hyperelastic softening and elastic instabilities
in graphene. Proceedings of the Royal Society of London A: Mathematical, Physical and
Engineering Sciences, 471(2173). DOI: 10.1098/rspa.2014.0567.
Kumar, S. and Parks, D. M. (2015). Strain shielding from mechanically activated covalent
bond formation during nanoindentation of graphene delays the onset of failure. Nano Lett.,
15(3):1503–1510.
Larsson, R. and Samadikhah, K. (2011). Atomistic continuum modeling of graphene membranes.
Comput. Mater. Sci., 50(5):1744–1753.
Lee, C., Wei, X., Kysar, J. W., and Hone, J. (2008). Measurement of the elastic properties and
intrinsic strength of monolayer graphene. Science, 321(5887):385–388.
Li, T. and Zhang, Z. (2010). Snap-through instability of graphene on substrates. Nanoscale
Res. Lett., 5(1):169–173.
Li, X., Wu, Z., Hou, L., Wang, L., and Yue, C. (2016). The accurate method for computing the
minimum distance between a point and an elliptical torus. Computers, 5(1):4.
Lindsay, L. and Broido, D. A. (2010). Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene.
Phys. Rev. B, 81:205441.
Lu, Q., Arroyo, M., and Huang, R. (2009). Elastic bending modulus of monolayer graphene. J.
Phys. D: Appl. Phys., 42(10):102002.
Lu, Z. and Dunn, M. L. (2010). Wan der Waals adhesion of graphene membranes. J. Appl.
Phys., 107(4).
Maboudian, R. and Howe, R. T. (1997). Critical review: Adhesion in surface micromechanical
structures. J. Vac. Sci. Technol., B, 15(1):1–20.

46

Marinho, B., Ghislandi, M., Tkalya, E., Koning, C. E., and de With, G. (2012). Electrical conductivity of compacts of graphene, multi-wall carbon nanotubes, carbon black, and graphite
powder. Powder Technol., 221:351–358.
Miehe, C. (1993). Computation of isotropic tensor functions. Commun. Numer. Methods Eng.,
9(11):889–896.
Montella, G., Govindjee, S., and Neff, P. (2016). The exponentiated Hencky strain energy in
modeling tire derived material for moderately large deformations. J. Eng. Mater. Technol.,
138(3):031008. DOI: 10.1115/1.4032749.
Neff, P., Eidel, B., Osterbrink, F., and Martin, R. (2013). The Hencky strain energy ‖ log 𝑈 ‖2
measures the geodesic distance of the deformation gradient to SO(n) in the canonical leftinvariant riemannian metric on GL (n). PAMM, 13(1):369–370.
Neff, P. and Ghiba, I.-D. (2016). The exponentiated Hencky-logarithmic strain energy: part
III coupling with idealized multiplicative isotropic finite strain plasticity. Continuum Mech.
Thermodyn., 28(1-2):477–487.
Neff, P., Ghiba, I.-D., and Lankeit, J. (2015a). The exponentiated Hencky-logarithmic strain
energy. part I: constitutive issues and rank-one convexity. J. Elast., 121(2):143–234.
Neff, P., Lankeit, J., Ghiba, I.-D., Martin, R., and Steigmann, D. (2015b). The exponentiated
Hencky-logarithmic strain energy. part II: Coercivity, planar polyconvexity and existence of
minimizers. Zeitschrift für angewandte Mathematik und Physik, 66(4):1671–1693.
Park, H. S., Klein, P. A., and Wagner, G. J. (2006). A surface cauchy-born model for nanoscale
materials. Int. J. Numer. Methods Eng., 68(10):1072–1095.
Pitteri, M. (1985). On 𝜈+ 1-lattices. J. Elast., 15(1):3–25.
Qomi, M. J. A., Aghaei, A., and Khoei, A. R. (2011). Multi-scale modeling of surface effect via
the boundary cauchy-born method. Int. J. Numer. Methods Eng., 85(7):827–846.
Roohbakhshan, F., Duong, T. X., and Sauer, R. A. (2016). A projection method to extract
biological membrane models from 3D material models. J. Mech. Behav. Biomed. Mater.,
58:90 – 104. Special issue: Mechanics of biological membranes.
Sauer, R. A. and De Lorenzis, L. (2013). A computational contact formulation based on surface
potentials. Comput. Methods in Appl. Mech. Eng., 253:369 – 395.
Sauer, R. A. and De Lorenzis, L. (2015). An unbiased computational contact formulation for
3D friction. Int. J. Numer. Methods Eng., 101(4):251–280.
Sauer, R. A. and Duong, T. X. (2017). On the theoretical foundations of thin solid and liquid
shells. Math. Mech. Solids, 22(3):343–371.
Sauer, R. A., Duong, T. X., and Corbett, C. J. (2014). A computational formulation for
constrained solid and liquid membranes considering isogeometric finite elements. Comput.
Methods in Appl. Mech. Eng., 271:48 – 68.
Sauer, R. A. and Li, S. (2007a). An atomic interaction-based continuum model for adhesive
contact mechanics. Finite Elem. Anal. Des., 43(5):384 – 396. The Eighteenth Robert J.
Melosh Competition.
Sauer, R. A. and Li, S. (2007b). A contact mechanics model for quasi-continua. Int. J. Numer.
Methods Eng., 71(8):931–962.
47

Sauer, R. A. and Wriggers, P. (2009). Formulation and analysis of a three-dimensional finite
element implementation for adhesive contact at the nanoscale. Comput. Methods in Appl.
Mech. Eng., 198(49 – 52):3871 – 3883.
Sfyris, D., Sfyris, G., and Galiotis, C. (2014a). Curvature dependent surface energy for a free
standing monolayer graphene: Some closed form solutions of the non-linear theory. Int. J.
Non Linear Mech., 67:186 – 197.
Sfyris, D., Sfyris, G., and Galiotis, C. (2014b). Curvature dependent surface energy for free
standing monolayer graphene: Geometrical and material linearization with closed form solutions. Int. J. Eng. Sci., 85:224 – 233.
Sfyris, D., Sfyris, G., and Galiotis, C. (2015). Constitutive modeling of some 2D crystals:
Graphene, hexagonal BN, MoS2 , WSe2 and NbSe2 . Int. J. Solids Struct., 66:98 – 110.
Sfyris, D. and Sfyris, G. I. (2017). Influence of partial blistering on the global and the local
stress and couple stress field for a monolayer graphene resting on substrate. Math. Mech.
Solids, page 1081286516684665. DOI: 10.1177/1081286516684665.
Shodja, H. M., Ojaghnezhad, F., Etehadieh, A., and Tabatabaei, M. (2017). Elastic moduli
tensors, ideal strength, and morphology of stanene based on an enhanced continuum model
and first principles. Mech. Mater., 110:1 – 15.
Simo, J. C. and Taylor, R. L. (1991). Quasi-incompressible finite elasticity in principal stretches.
continuum basis and numerical algorithms. Comput. Methods in Appl. Mech. Eng., 85(3):273
– 310.
Sun, F. and Li, H. (2011). Torsional strain energy evolution of carbon nanotubes and their
stability with encapsulated helical copper nanowires. Carbon, 49(4):1408 – 1415.
Wei, X., Fragneaud, B., Marianetti, C. A., and Kysar, J. W. (2009). Nonlinear elastic behavior
of graphene: Ab initio calculations to continuum description. Phys. Rev. B, 80:205407.
Xu, M., Paci, J. T., Oswald, J., and Belytschko, T. (2012). A constitutive equation for graphene
based on density functional theory. Int. J. Solids Struct., 49(18):2582 – 2589.
Zheng, Q.-S. (1994). Theory of representations for tensor functions: a unified invariant approach
to constitutive equations. Appl. Mech. Rev., 47(11):545–587.

48

