Curvature-dependent surface energy for free-standing monolayer graphene

Abstract

The purpose of the present work is to give a continuum model that can capture bending effects for free-standing graphene monolayers taking material’s symmetry properly into account. Starting from the discrete picture of graphene modelled as a hexagonal 2-lattice, we give the arithmetic symmetries. Confined to weak transformation neighbourhoods one is able to work with the geometric symmetry group. Use of the Cauchy–Born rule allows the transition from the discrete case to the continuum case. At the continuum level we use a surface energy that depends on an in-plane strain measure, the curvature tensor and the shift vector. Dependence of the energy on the curvature tensor allows for incorporating bending effects into the model. Dependence on the shift vector is motivated by the fact that discretely graphene is a 2-lattice. We lay down the complete and irreducible set of invariants for this surface energy amenable to available representation theory. This way we obtain the expression for the surface stress as well as the surface couple stress tensor, the first being responsible for the in-plane deformations and the second for the out-of-plane motions. Forms for the elasticities of the material are given accompanied by the field equations. The model, in its simplest form, predicts 13 independent scalar variables in the constitutive relations to be observed in experiments. The framework presented is valid for both materially and geometrically nonlinear theories. We also present the case where symmetry changes at the continuum level, without taking into account how energy behaves at the transition regime.

Keywords

Monolayer graphene anisotropy hyperelasticity curvature tensor shift vector

1. Introduction

Graphene is a two-dimensional sheet that constitutes the building unit of all graphitic forms of matter, such as graphite, carbon nanotubes and carbon fibres. For modelling graphene many different approaches at different scales can be found in the literature, including first-principle calculations [1, 2], atomistic calculations [3, 4] and continuum mechanics [5 –7]. Furthermore, mixed atomistic formulations with finite elements are being reported for graphene [8 –10] based on the earlier notion of a quasi-continuum [11, 12].

In recent literature there is a debate related to the use of continuum mechanics for describing nanomaterials in general and graphene in particular. In [13] it is argued that the plate idealization fails for monolayer graphene due to the pure bending character of the model that does not take into account in-plane stretching. With an error of approximately 6% for N = 3 (N is the number of layers) the plate phenomenology works well. Another work [14] supports the failure of continuum mechanics when rippling of graphene takes place. The arguments are the failure of the classical continuum framework in taking into account plates’ thickness and also the inability of the model presented in [15 –17] to determine the amplitude of a wave-like solution that represents the out-of-plane displacement in buckling.

It is worth mentioning in this respect that these approaches [13 –17] are based on Foppl–von Karman plate theory where the potential energy scales to h⁴ (h is plate thickness) and thus appears not to be able to support significant stretching [18]. It also assumes an isotropic response under in-plane loading and a transversely isotropic response for the out-of-plane [19] loading, in line with the fundamental approach of Timoshenko and Gere [20]. For the purpose of modelling significant stretching more refined plate theories have been reported [18] that take into account bending effects as well as stretching in the nonlinear regime. The exact assumptions of the approaches in [13 –17] under the perspective of nonlinear incremental elasticity can be found in the recent work of Steigmann and Ogden [21].

From the mathematical point of view a promising tool for bridging the length scales appears to be the technique of Gamma convergence [22]. Using this theory Friesecke et al. [23] proposed a hierarchy of plate models under pure bending derived from nonlinear elasticity highlighting the fact that plate models are able to take into account the thickness of the layer they model. More refined continuum models (Cosserat surfaces) are proposed by the same authors [24, 25] in order to model thin films at the nano scale.

The mathematical theory of surface elasticity is established by Gurtin and Murdoch [26]. This pure membrane approach is incapable of taking into account out-of-plane deformations. Generalization of this framework to take into account bending effects is given by Steigmann and Ogden [27]. These authors propose a surface energy which depends, apart from a surface measure of the deformation, on the curvature tensor as well, in trends similar to previous works [28, 29]. The curvature tensor is a measure of the out-of-plane deformations the surface suffers and this way bending effects are introduced into the framework. Steigmann and Ogden, in the same work [27], also describe a rigorous way of defining the notion of material symmetry for curvature-dependent surface energies in line with Noll’s fundamental work [30]. Implications of such energies for nanostructures are studied in Chhapadia et al. [31].

We here adopt the framework of Steigmann and Ogden [27] and utilize a surface energy function depending on three arguments for a free-standing monolayer graphene. The first one is an in-surface strain measure describing changes happening on the surface. The second argument is the curvature tensor which describes the out-of-surface motions and introduces bending effects into the model. The third argument is the shift vector. The motivation for assuming the shift vector as an independent variable comes from the work of Pitteri and Zanzotto ([32] and references therein). These authors utilize an energy function depending on the shift vector when modelling a multilattice. We note that for graphene a similar assumption is made by E and Ming [33]. Section 4 presents the basic prerequisites from surface elasticity modelling as well as curvature-dependent surface energies. Section 5 gives material symmetry and field equations for bodies described by curvature-dependent surface energies.

Using the surface energy, calculation of the surface stress and the surface couple stress tensor at the continuum level is possible. This way the number of independent relations to be observed in experiments becomes available. The surface stress tensor is responsible for in-plane motions while out-of-plane motions are due to the surface couple stress tensor. The elasticities of the material can then be calculated and we also give the field equations characterizing the problem; these are the momentum and the moment of momentum equations as well as an equation for the evaluation of the shift vector. These are the contents of Section 6.

Being aware of the molecular theories of elasticity, where the energy depends on the lattice vectors, we stress that our analysis is confined to weak transformation neighbourhoods [34]. This way the classical theory of invariants for continuum mechanics can be utilized. This is compatible with molecular theories when the Cauchy–Born rule is enforced, and also compatible with the global theory of Ericksen [35, 36]. In this respect, we give a short account of the arithmetic symmetries of graphene as a two-dimensional lattice in Section 2 based on the work of Fadda and Zanzotto [37]. Section 3 describes the basic prerequisites for application of the Cauchy–Born rule. Confinement to weak transformation neighbourhoods and the Cauchy–Born rule are the basic assumptions for the passage from the discrete picture to the continuum case.

Section 7 treats changes of symmetry at the continuum level without taking into account what happens at the transition regime. Finally, the paper ends in Section 8 with some concluding remarks. Einstein summation convention is enforced for repeated indices. The dot and tensor product are defined in line with common practise as a ⋅ b = a_ib_i, a ⊗ b = a_ib_j. Contraction in two indices is denoted by the symbol ‘:’, for example, for two tensors of second order it is A: B = A_ijB_ij. Contraction in more than two indices is denoted by •.

2. Arithmetic symmetries of graphene

The importance of the difference between arithmetic and geometric symmetries for crystals stems from the fundamental work of Ericksen [35]. We refer to the book of Pitteri and Zanzotto ([32], and references therein) for a nice exposition of the topic.

A three-dimensional simple lattice in $ℛ^{3}$ is defined as (see e.g. [32])

ℒ (e_{a}) = {x \in ℛ^{3} : x = M^{a} e_{a}, a = 1, 2, 3, M^{a} \in Z},

(1)

e_a being the lattice vectors and $Z$ the space of positive integers. The geometric symmetry group of $ℒ$ is [32, 35]

\begin{array}{l} P (e_{a}) & = & {Q \in O : ℒ (Q e_{a}) = ℒ (e_{a}) \\ = & {Q \in O : Q e_{a} = m_{a}^{b} e_{b}, m \in G L (3, Z)}, \end{array}

(2)

with $G L (3, Z)$ and $O$ being the general linear and the orthogonal groups, respectively. Essentially, this group gives all orthogonal transformations Q that map $ℒ$ to itself. Here one finds, for a three-dimensional lattice, the seven crystal systems and the 28 crystallographic point groups continuum mechanics utilizes.

Due to the fact that conjugacy in $G L (3, Z)$ is more stringent than conjugacy in $O$ , a finer description of symmetry is given by the arithmetic symmetry which is defined as [32, 35]

L (e_{a}) = {m \in G L (3, Z) : m_{a}^{b} e_{b} = Q e_{a}, Q \in P (e_{a})} .

(3)

This group gives all distinct types of lattices that are compatible with a given geometric group. Essentially, this is a finer description of symmetry that can differentiate between Bravais lattice types within the same crystal system.

Multilattice is a generalization of a simple lattice in the sense that it is the finite union of translates of some suitable simple lattice

ℳ (p_{i}, e_{a}) = \cup_{i = 0}^{n - 1} ℒ (p_{i}, e_{a}) .

(4)

The particular case of a 2-lattice is the union of two simple lattices

ℳ (e_{a}, p) = ℒ (e_{a}) \cup {p + ℒ (e_{a})},

(5)

$ℒ (e_{a})$ being the simple lattice generated by the basis e_a, while the vector

p = p^{1} e_{1} + p^{2} e_{2}

(6)

is called the shift vector and gives the separation of the two simple lattices constituting $ℳ$ . The geometric and arithmetic symmetry groups of multilattices are defined in a similar fashion to the corresponding definitions of simple lattices and we refer to [32] for more information.

Following the classification of 2-lattices by Fadda and Zanzotto [37], we treat monolayer graphene as a hexagonal monoatomic 2-lattice with unit cell of the form in Figure 1. The lattice and shift vectors are depicted in Figure 2 and defined as

e_{1} = (\sqrt{3} l, 0), e_{2} = (\frac{\sqrt{3}}{2} l, \frac{3}{2} l), p = (\frac{\sqrt{3}}{2} l, \frac{1}{2} l),

(7)

Figure 1.

The unit cell of a hexagonal 2-lattice [37].

Figure 2.

The lattice and shift vectors of graphene.

l being the lattice size, namely the interatomic distance at ease which is approximately 1.42 Å. The two simple hexagonal lattices are

\begin{array}{l} L_{1} (l) = {x \in ℛ^{2} : x = n^{1} e_{1} + n^{2} e_{2}, (n^{1}, n^{2}) \in Z^{2}}, \\ L_{2} (l) = p + L_{1} (l) . \end{array}

(8)

Searching for the arithmetic symmetry of 2-lattices equal seeking matrices μ such that [37]

μ^{T} K μ = K,

(9)

K being the lattice metric

K = ε_{a} \cdot ε_{b}, a, b = 1, 2, 3, ε_{1} = e_{1}, ε_{2} = e_{2}, ε_{3} = p .

(10)

In matrix form the lattice metric for graphene is

K = (\begin{matrix} 3 l^{2} & \frac{3}{2} l^{2} & \frac{3}{2} l^{2} \\ \frac{3}{2} l^{2} & \frac{3}{2} l^{2} & \frac{3}{2} l^{2} \\ \frac{3}{2} l^{2} & \frac{3}{2} l^{2} & l^{2} \end{matrix}) .

(11)

Following [37] the arithmetic symmetry group obtained from the proper fixed set in such a case is described by the matrices

(\begin{matrix} - 1 & - 1 & - 1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}), (\begin{matrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}), (\begin{matrix} - 1 & - 1 & - 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}),

(12)

(\begin{matrix} 1 & 0 & 0 \\ - 1 & - 1 & - 1 \\ 0 & 0 & 1 \end{matrix}), (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}), (\begin{matrix} 0 & 1 & 0 \\ - 1 & - 1 & - 1 \\ 0 & 0 & 1 \end{matrix}),

(13)

since for the fixed set intersecting the fundamental domain we have K₁₃ = K₂₃ = K₁₁/2. Evaluating the eigenvalues of these matrices one will find the values $1, - 1, e^{i \frac{π}{3}}, e^{- i \frac{π}{3}}$ . Thus, these matrices describe identity transformation, reflection transformation and rotations by 60 and −60 degrees.

Following [32], the fixed set mentioned above is the collection of all those metrics for which equation (9) holds. The fundamental domain is a region in the space of all metrics such that the orbit of the arithmetic symmetry group has only one element in the space of {∊_a}. Symmetry will first break to the rhombic arithmetic type and then to the oblique arithmetic type [37].

For the analysis to follow it is instructive to delineate the conditions under which arithmetic symmetry may reduce to the geometric symmetry [34, 38, 39]. This reconciliation can be accomplished when one is confined to weak transformation neighbourhoods. In short, for simple lattices say that $ℬ$ is the nine-dimensional space of all bases e_a. Then one can prove [34, 38, 39] that there exists a neighbourhood $N_{1} \supset ℬ$ such that the action of L(e_a) coincides with that of P(e_a). Namely, when one is confined to such neighbourhoods there is no need to distinguish between geometric and arithmetic symmetry groups. Similar arguments hold for multilattices as well.

To sum up, when graphene is modelled as a hexagonal 2-lattice, its arithmetic symmetries are given by the matrices of equations (12) and (13). These matrices describe identity transformation, reflection transformation and rotations by 60 and −60 degrees.

3. Passage to the continuum: The Cauchy–Born rule

For passing the discrete picture of a lattice to its continuum analog, the Cauchy–Born rule should be used (see e.g. [36] and references therein). Limitations of this rule can be found in [40] (see also [33]). In the present work validity of this rule is enforced.

According to this rule, for multilattices, atomic motion of the lattice agrees with the gross deformation, while the shift vector is free to adjust so as to reach equilibrium. Thus, if we assume the existence of a stored energy function φ for the multilattice we have

φ = φ (e_{a}, p) .

(14)

The Cauchy–Born rule then states that the reference and the current lattice vectors $e_{a}^{0}$ and e_a, respectively, are related according to the formula [32]

e_{a} = F e_{a}^{0},

(15)

F = ∂x/∂X being the well-known deformation gradient of continuum mechanics. The shift vector is adjusted so as to reach equilibrium. Using minimization arguments one may show that [32]

\frac{\partial φ}{\partial p} = 0 .

(16)

This equation plays the role of the field equation for the evaluation of p. So, using the Cauchy–Born rule at the level of the energy one is able to write

φ = φ (F, p) .

(17)

Classical invariance then would require

φ (F, p) = φ (Q F H^{T}, p H^{T}), H \in L (e_{a}), Q \in Q, F \in N_{e_{a}},

(18)

but when confined to a weak transformation neighbourhood,

φ (F, p) = φ (Q F H^{T}, p H^{T}), H \in P (e_{a}), Q \in Q, F \in N_{1} .

(19)

The action of H is due to material symmetry, while the action of Q is due to frame indifference. So, confined to this neighbourhood, material symmetry uses the geometric symmetry group, namely the crystal systems continuum mechanics is akin to. Frame indifference then leads to

φ = φ (C, p),

(20)

where C = F^TF is the right Cauchy–Green deformation tensor. So, the Cauchy–Born rule permits the transition from a lattice to its continuum counterpart (see Section 6.2, [32]). What is new for multilattices is the dependence on the shift vector as well [32, 33]. The motivation for the transformation rule for the shift vector under the action of the symmetry group in equations (18) and (19) is given in Section 5.

To sum up, there are two crucial assumptions for the validity of the framework to follow: the Cauchy–Born rule is enforced and we confine ourselves to weak transformation neighbourhoods. This way one is able to work with an energy of the form

φ = φ (C, p),

(21)

augmented properly to take into account bending effects, while symmetry is the one continuum mechanics uses.

4. Curvature-dependent surface energy

The arguments of Section 3 pertain to classical three-dimensional (bulk) elasticity. In this section we lay down the kinematics of surface elasticity [27, 41]. We assume in the reference configuration a smooth surface A₀, described by Y = Y (θ¹, θ²), where Y is the position vector for the point of the surface from the origin and the parameters θ^α α = 1, 2 are curvilinear coordinates on the surface. Covariant and contravariant vectors are defined by

A_{α} = Y_{, α}, A^{α} \cdot A_{β} = δ_{β}^{α},

(22)

where $δ_{β}^{α}$ is the Kronecker delta in two dimensions. Deformation of the surface brings point Y to point y(θ¹, θ²) on surface A, in the current configuration. The covariant and contravariant vectors related to the surface A are then defined:

a_{α} = y_{, α}, a^{α} \cdot a_{β} = δ_{β}^{α} .

(23)

The linear mapping that maps vectors in the tangent plane of A₀ to those of A is the surface deformation gradient defined by

F_{s} = a_{α} \otimes A^{α} .

(24)

The surface right Cauchy–Green deformation tensor is then defined as

C_{s} = F_{s}^{T} \cdot F_{s} .

(25)

The surface curvature tensors in the reference and the current configuration can be expressed as

b_{0} = b_{0 α β} A^{α} \otimes A^{β},

(26)

b = b_{α β} a^{α} \otimes a^{β} .

(27)

Surface divergence is defined by

\nabla_{s} () = \nabla () - n (n \cdot \nabla ()),

(28)

where ∇ is the common divergence operator in three dimensions, while n is the outward unit normal for the surface A. The outward unit normal for the surface A₀ is denoted by N.

A curvature-dependent surface energy is an energy of the form [27 –29]

W = W (C_{s}, b_{0}) .

(29)

So, aside from the in-surface strain measure C_s, dependence on b₀ is assumed. This dependence allows the modelling of bending effects since it takes into account out-of-plane deformations. For monolayer graphene this energy should be augmented to

W = W (C_{s}, b_{0}, p),

(30)

where p is the shift vector. The need for the dependence on the shift vector stems from the fact that graphene is a multilattice. A similar dependence for graphene can be found in the work of E and Ming [33].

5. Material symmetry and field equations

Material symmetry for curvature-dependent surface energies is a subject tackled elegantly by Steigmann and Ogden [27, Section 6] based on the earlier work of Murdoch and Cohen [29]. For such energies, elements of the symmetry group are pairs which leave the response of the surface invariant to superposed deformations. Steigmann and Ogden [27] concluded that for surfaces with constant non-negative curvature (namely, for planes and spheres) amenability to available representation theory is possible when the symmetry group reads {R, 0}, while for the energy it then holds that

W (C_{s}, b_{0}) = W (R C_{s} R^{T}, R b_{0} R^{T}), R R^{T} = I, \det R = 1 .

(31)

The form of the energy should be augmented for graphene by taking into account dependence on the shift vector. By assuming the shift vector to be a tangent vector of the referential surface we see that action of the symmetry group would change p as follows:

p^{*} = p R^{T} .

(32)

Justification of this formula stems from the way the surface deformation gradient changes under the action of R. F_s, being a two-point tensor, is affected up to its ‘leg’ on the reference surface, $F_{s}^{*} = F_{s} R^{T}$ . So, only a tangent vector on the reference surface is affected by the action of the symmetry group. The assumption that p is a referential tangent vector results then in equation (32). So, collectively, the action of the symmetry group for a curvature-dependent surface energy for monolayer graphene reads

W (C_{s}, b_{0}, p) = W (R C_{s} R^{T}, R b_{0} R^{T}, p R^{T}) .

(33)

Field equations for materials described by curvature-dependent surface energies are studied in [31]. According to these authors, the momentum equation for the static case and in the absence of body forces reads

σ^{bulk} \cdot n + \nabla_{S} T_{S} = 0,

(34)

where T_S is the surface first Piola–Kirchhoff stress tensor defined by [27]:

T_{S} = \frac{\partial \bar{W}}{\partial F_{S}}

(35)

where $W = \bar{W} (F_{S}, b_{0}, p)$ , while σ ^bulk is the Cauchy stress tensor for the bulk material surrounding the surface. Here, since we speak about a free-standing surface, there is no bulk material surrounding the graphene, so the bulk stress tensor should be set equal to zero, σ ^bulk = 0. The surface first Piola–Kirchhoff stress tensor is related to the second Piola–Kirchhoff surface stress tensor, S_S, according to the formula

S_{S} = F_{S}^{- 1} \cdot T_{S} .

(36)

The second Piola–Kirchhoff stress tensor can also be written as

S_{S} = \frac{\partial W}{\partial C_{S}} .

(37)

The momentum equation in the absence of body forces and inertia can be expressed using the second Piola–Kirchhoff stress tensor as

{\bar{\nabla}}_{S} S_{S} = 0,

(38)

where ${\bar{\nabla}}_{S} () = \nabla^{C} () - N (N \cdot \nabla^{C} ())$ is the surface gradient when the bulk gradient ∇^C is taken with respect to the metric C of the bulk. In this respect, the bulk divergence in equation (28) is taken with respect to the metric G of the bulk reference configuration. Essentially, equation (38) is the surface analog of the momentum equation written using the second Piola–Kirchhoff stress tensor in classical elasticity [42 –44].

The moment of momentum balance, in the absence of body couples and inertia, when setting σ ^bulk = 0, reads [31]

\nabla_{S} F_{S} m_{S} - \nabla_{S} (F_{S} \cdot S_{S} \times y) = 0,

(39)

where the surface couple stress tensor is defined as [27]

m_{S} = \frac{\partial W}{\partial b_{θ}} .

(40)

The symbol × in equation (39) denoted the cross-product of the three-dimensional space.

For the shift vector the field equation reads [32, 33]

\frac{\partial W}{\partial p} = 0 .

(41)

6. Application to graphene

The case of graphene considered as a three-dimensional material has as symmetry group class number 27 denoted by D_6h or P6/mmm [45, 46]. When viewed as a two-dimensional continuum we assume that it belongs to class number 10 of the classification adopted in [46, 47]; this group is hereafter denoted $D_{6 h}$ . The generators of this group are R(2π/6) and R_j, R(θ) representing rotation, θ being the angle of rotation and R_j being a two-dimensional reflection transformation. These generators describe the same transformation as the matrices of equations (12) and (13). In this section we give the complete and irreducible representation of an energy of the form of equation (33) when the symmetry group is $D_{6 h}$ .

The structure tensor for $D_{6 h}$ is denoted by P₆ and defined by [48]

P_{6} = R e {(i + i j)}^{6},

(42)

or equivalently as

P_{6} = U \otimes U \otimes U - (L \otimes U \otimes L + L \otimes L \otimes U),

(43)

where U = a₁ ⊗ a₁ − a₂ ⊗ a₂, L = a₁ ⊗ a₂ − a₂ ⊗ a₁, for {a₁, a₂} an orthonormal basis vector. It can also be written as

P_{6} = R e [e^{i 6 θ} {(c_{1} + i c_{2})}^{6}],

(44)

where c₁ = cos(θ)a₁ + sin(θ)a₂, c₂ = −sin(θ)a₁ + cos(θ)a₂. This tensor is an irreducible tensor since $D_{6 h}$ is compact. In two dimensions it has only two independent components [46],

P_{111111} = \cos (6 θ), P_{211111} = \sin (6 θ) .

(45)

These two components introduce the anisotropy and for graphene, they model the zigzag and armchair directions. Since θ = 2π/6 we have P₁₁₁₁₁₁ = cos(2π) = 1, P₂₁₁₁₁₁ = sin(2π) = 0.

Thus, for a graphene monolayer modelled using a curvature-dependent surface energy we have

W = W_{a n i s o t r o p i c} (C_{S}, b_{0}, p),

(46)

where the anisotropy stems from the fact that the symmetry group is $D_{6 h}$ and not the full rotation group. Reduction to an isotropic function is done through the use of the principle of isotropy of space [46], that gives

W = W_{a n i s o t r o p i c} (C_{S}, b_{0}, p) = W_{i s o t r o p i c} (C_{S}, b_{0}, P_{6}, p) .

(47)

Thus, we take an isotropic function at the expense of using the structure tensor (that describes the anisotropy) as an additional argument.

The complete and irreducible representation of such a scalar function under the group $D_{6 h}$ consists of the following quantities [46, 48]:

\begin{array}{l} I_{1} = tr C_{S}, I_{2} = \det C_{S}, I_{3} = tr (Π_{6}^{C_{S}} C_{S}), I_{4} = tr (C_{S} b_{0}), \\ I_{5} = tr (Π_{6}^{b_{0}} b_{0}), I_{6} = tr b_{0}, I_{7} = \det b_{0}, I_{8} = p \cdot C_{S} p, \\ I_{9} = p \cdot b_{0} p, I_{10} = p \cdot π_{6}^{p}, I_{11} = p \cdot p, I_{12} = tr (Π_{6}^{p} C_{S}), I_{13} = tr (Π_{6}^{p} b_{0}) . \end{array}

(48)

Thus, in general for such a model of graphene we have the following expression for the energy:

W = \tilde{W} (I_{i}), i = 1, 2, \dots, 13 .

(49)

The term $Π_{6}^{A}$ for a symmetric tensor of second order A is defined as ([46], indices ranging from 1 to 2)

Π_{6}^{A} = P_{i j k l m n} A_{k l} A_{m n} c_{i} \otimes c_{j},

(50)

and renders a second-order tensor. The term $π_{6}^{z}$ with respect to the vector z is defined as

π_{6}^{z} = P_{i j k l m n} z_{j} z_{k} z_{l} z_{m} z_{n} c_{i},

(51)

while for $Π_{6}^{z}$ we have

Π_{6}^{z} = P_{i j k l m n} z_{k} z_{l} z_{m} z_{n} c_{i} \otimes c_{j} .

(52)

The material parameters related to I₆, I₇ describe pure bending effects since detb₀, trb₀ are the mean and the Gaussian curvature of the surface, respectively. The term denoted by I₅ describes the effect of the armchair and the zigzag direction of graphene at bending. The parameters related to I₁, I₂ are related to pure stretching, while those related to the term I₃ describe the anisotropy (zigzag, armchair directions) to stretching. The parameter related to I₄ describes coupling between bending and stretching response. The terms I₈ and I₉ describe the effect of the in-plane and the out-of-plane deformations, respectively, on the shift vector. The term I₁₀ describes the way anisotropy affects the shift vector, while I₁₁ describes changes related to the shift vector solely. Terms I₁₂ and I₁₃ are related to coupling of anisotropy with the shift vector for the in-plane and the out-of-plane deformations, respectively. The zigzag and armchair directions of a graphene sheet are depicted in Figures 3 and 4.

Figure 3.

The armchair direction of graphene is introduced into the mathematical framework through the tensors $Π_{6}^{C_{S}}, Π_{6}^{b_{0}}, π_{6}^{p}, Π_{6}^{p}$ .

Figure 4.

The zigzag direction of graphene is introduced into the mathematical framework through the tensors $Π_{6}^{C_{S}}, Π_{6}^{b_{0}}, π_{6}^{p}, Π_{6}^{p}$ .

For evaluating the surface stress tensor and the surface couple stress tensor one has to determine the derivatives

S_{S} = \frac{\partial \tilde{W}}{\partial C_{S}}, m_{S} = \frac{\partial \tilde{W}}{\partial b_{0}} .

(53)

The term $\partial \tilde{W} / \partial p$ is also important due to equation (41). Using the expressions of equation (48) in equation (53), after some calculations we obtain

\begin{array}{l} S_{S} = \frac{\partial \tilde{W}}{\partial I_{1}} I + \frac{\partial \tilde{W}}{\partial I_{2}} [tr (C_{S}) 1 - C_{S}] + 3 \frac{\partial \tilde{W}}{\partial I_{3}} P_{6} : (C_{S} \otimes C_{S}) + \frac{\partial \tilde{W}}{\partial I_{4}} b_{0} \\ + \frac{\partial \tilde{W}}{\partial I_{8}} p \otimes p + \frac{\partial \tilde{W}}{\partial I_{12}} Π_{6}^{p}, \end{array}

(54)

\begin{array}{l} m_{S} = \frac{\partial \tilde{W}}{\partial I_{6}} I + \frac{\partial \tilde{W}}{\partial I_{7}} [tr (b_{0}) 1 - b_{0}] + 3 \frac{\partial \tilde{W}}{\partial I_{5}} P_{6} : (b_{0} \otimes b_{0}) + \frac{\partial \tilde{W}}{\partial I_{4}} C_{S} \\ + \frac{\partial \tilde{W}}{\partial I_{9}} p \otimes p + \frac{\partial \tilde{W}}{\partial I_{13}} Π_{6}^{p}, \end{array}

(55)

\begin{array}{l} \frac{\partial W}{\partial p} = 2 \frac{\partial \tilde{W}}{\partial I_{8}} C_{S} p + 2 \frac{\partial \tilde{W}}{\partial I_{9}} b_{0} p + 6 \frac{\partial \tilde{W}}{\partial I_{10}} P_{6} • (p \otimes p \otimes p \otimes p \otimes p) + \frac{\partial \tilde{W}}{\partial I_{11}} p \\ + 4 \frac{\partial \tilde{W}}{\partial I_{12}} [P_{6} : (p \otimes p \otimes p)] : C_{S} + 4 \frac{\partial \tilde{W}}{\partial I_{13}} [P_{6} : (p \otimes p \otimes p)] : b_{0} . \end{array}

(56)

By making the simplest possible assumption that $\tilde{W}$ is linear with respect to the invariants I_i, i = 1, 2, 3,…, 13, we take

\begin{array}{l} S_{S} = α I + β [tr (C_{S}) 1 - C_{S}] + 3 γ P_{6} : (C_{S} \otimes C_{S}) \\ + δ b_{0} + θ p \otimes p + ρ Π_{6}^{p}, \end{array}

(57)

m_{S} = ε I + ζ [tr (b_{0}) 1 - b_{0}] + 3 η P_{6} : (b_{0} \otimes b_{0}) + δ C_{S} + ι p \otimes p + τ Π_{6}^{p},

(58)

\begin{array}{l} \frac{\partial W}{\partial p} = & θ C_{S} p + ι b_{0} p + 6 λ P_{6} • (p \otimes p \otimes p \otimes p \otimes p) + ξ p \\ + 4 ρ [P_{6} : (p \otimes p \otimes p)] : C_{S} + 4 τ [P_{6} : (p \otimes p \otimes p)] : b_{0} . \end{array}

(59)

The Greek letters α, β, γ, δ, θ, ρ, ε, ζ, η, ι, τ, λ, ξ are material parameters to be determined by experiments.

Written with respect to indices ranging from 1 to 2 these formulae are

\begin{array}{l} S_{S_{A B}} = & α δ_{A B} + β [t r (C_{S}) δ_{A B} - C_{S_{A B}}] + 3 γ P_{A B C D E F} C_{S_{E F}} C_{S_{C D}} + δ d_{0_{A B}} \\ + θ p_{A} p_{B} + ρ P_{A B C D E F} p_{C} p_{D} p_{E} p_{F}, \end{array}

(60)

\begin{array}{l} m_{S_{A B}} = & ε δ_{A B} + ζ [t r (b_{0}) δ_{A B} - b_{0_{A B}}] + 3 η P_{A B C D E F} b_{0_{E F}} b_{0_{C D}} + δ C_{S_{A B}} \\ + ι p_{A} p_{B} + τ P_{A B C D E F} p_{C} p_{D} p_{E} p_{F}, \end{array}

(61)

\begin{array}{l} \frac{\partial W}{\partial p_{A}} = & θ C_{S_{A B}} p_{A} + ι b_{0_{A B}} p_{B} + 6 λ P_{A B C D E F} p_{B} p_{C} p_{D} p_{E} p_{F} + ξ p_{A} \\ + 4 ρ P_{A B C D E F} p_{D} p_{E} p_{F} C_{S_{B C}} + 4 τ P_{A B C D E F} p_{D} p_{E} p_{F} b_{0_{B C}} . \end{array}

(62)

The elasticities of this model are given by the following fourth-order tensors:

A = \frac{\partial^{2} W}{\partial C_{S}^{2}}, ℬ = \frac{\partial^{2} W}{\partial b_{0}^{2}}, C = \frac{\partial^{2} W}{\partial C_{S} \partial b_{0}} .

(63)

Quantities of the first term are related to the in-plane motion and those of the second to the out-of-plane motion, while those of the third relate to the coupling between in-plane and out-of-plane motions.

7. Changes of symmetry

Recent ab initio calculations [49] reveal changes of symmetry for graphene when subjected to strain. More specifically, for strain at arbitrary directions the point group reduces from D_6h to C_2h, and for strain direction from 0⁰ to 30⁰, to D_2h. Similar reductions of symmetry for the general case of carbon nanotubes are reported in [6]. In this section we present the constitutive relations when such symmetry changes are present without taking into account what happens to the transition regime between the symmetry changes.

Groups D_2h, C_2h as three-dimensional point groups occupy group numbers 8 and 5, respectively, in common tables of compact point groups in three dimensions [46]. Since we here model graphene as a two-dimensional material surface, these groups should be related to their counterparts in two dimensions. So, viewed as two-dimensional groups, D_2h and C_2h should correspond to the group C_2v with generators −1 = R(π), R_e = 1 − 2r⊗r ; R_e is a reflection transformation with respect to the unit tensor r, following the two-dimensional classification of point symmetry groups by Zheng [47, 48]. This case corresponds to orthotropy in two dimensions. The structure tensor in this case is the tensor M = i ⊗ I − j ⊗ j, with {i,j} an orthonormal base in two dimensions.

Thus, in this case we have

W_{a n i s o t r o p i c} (C_{S}, b_{0}, p) = W_{i s o t r o p i c} (C_{S}, b_{0}, p, M) .

(64)

The complete and irreducible representation of such a scalar function under the group C_2v consists of the following quantities [48]:

\begin{array}{l} I_{1} = tr C_{S}, I_{2} = \det C_{S}, I_{3} = tr (M C_{S}), I_{4} = tr (C_{S} b_{0}), \\ I_{5} = tr (M b_{0}), I_{6} = tr b_{0}, I_{7} = \det b_{0}, I_{8} = p \cdot p, \\ I_{9} = p \cdot M p, I_{10} = p \cdot C_{S} p, I_{11} = p \cdot b_{0} p . \end{array}

(65)

Thus, for the energy we obtain

W = \tilde{W} (I_{i}), i = 1, 2, \dots, 11 .

(66)

The surface stress and surface couple stress that correspond to this energy are calculated as

\begin{array}{l} S_{S} = & \frac{\partial \tilde{W}}{\partial I_{1}} I + \frac{\partial \tilde{W}}{\partial I_{2}} [tr (C_{S}) 1 - C_{S}] \\ + \frac{\partial \tilde{W}}{\partial I_{3}} M : (C_{S} \otimes C_{S}) + \frac{\partial \tilde{W}}{\partial I_{4}} b_{0} + \frac{\partial \tilde{W}}{\partial I_{10}} p \otimes p, \end{array}

(67)

\begin{array}{l} m_{S} = & \frac{\partial \tilde{W}}{\partial I_{4}} I + \frac{\partial \tilde{W}}{\partial I_{7}} [tr (b_{0}) 1 - b_{0}] + \frac{\partial \tilde{W}}{\partial I_{5}} M : (b_{0} \otimes b_{0}) \\ + \frac{\partial \tilde{W}}{\partial I_{4}} C_{S} + \frac{\partial \tilde{W}}{\partial I_{11}} p \otimes p . \end{array}

(68)

For the term ∂W/∂p we have

\frac{\partial W}{\partial p} = \frac{\partial \tilde{W}}{\partial I_{8}} p + \frac{\partial \tilde{W}}{\partial I_{9}} M p + \frac{\partial \tilde{W}}{\partial I_{10}} C_{S} p + \frac{\partial \tilde{W}}{\partial I_{11}} b_{0} p .

(69)

By making the simplest possible assumption that $\tilde{W}$ is linear with respect to the invariants I_i, i = 1, 2, 3, …, 11, we take

S_{S} = α G + β [tr (C_{S}) 1 - C_{S}] + γ M : (C_{S} \otimes C_{S}) + δ H + θ p \otimes p,

(70)

m_{S} = ε G + ζ [tr (H_{S}) 1 - H_{S}] + η M : (b_{0} \otimes b_{0}) + δ C_{S} + ι p \otimes p

(71)

and

\frac{\partial W}{\partial p} = λ p + μ M p + θ C_{S} p + ι b_{0} p .

(72)

Now the effect of anisotropy at the level of the constitutive law is introduced through the terms where the structure tensor, M, is present. Namely the material parameters γ, η and μ measure the effect of anisotropy at the in-plane motion, the out-of-plane motion and the shift vector, respectively. For these cases the number of independent material constants to be observed in experiments is 11: the Greek letters α, β, γ, δ, θ, ε, ζ, η, ι, λ, μ.

8. Conclusion

The purpose of the present paper is to give a generic framework for modelling graphene monolayers at the continuum level. The starting point is the discrete nature of graphene: a hexagonal 2-lattice in two dimensions. Passage to the continuum rests upon two fundamental assumptions: the Cauchy–Born rule and restriction to weak transformation neighbourhoods. When passing to the continuum, we use curvature-dependent surface energies to capture bending effects. The model is general enough in the sense that it is valid for both material and geometrical nonlinearities, properly taking into account graphene’s symmetry. We evaluate the stress and the couple stress tensor and lay down the field equations describing such a model.

Footnotes

Acknowledgements

Valuable discussions with G Dassios, C Papaggelis, G Kalosakas (Patras, Greece) and GI Sfyris (Paris, France) are really appreciated. Special thanks go to E Koukaras (Patras, Greece) for drawing the figures as well as for his comments regarding the manuscript.

Funding

This research has been co-financed by the European Union (European Social Fund – ESF) and Greek national funds through the operational program ‘Education and Lifelong Learning’ of the National Strategic Reference Framework (NSRF) – Research Funding Program: ERC-10 ‘Deformation, Yield and Failure of Graphene and Graphene-Based Nanocomposites’.

References

Kudin

Scuseria

Yakobson

. C₂ F, BN and C nanoshell elasticity from ab-initio calculations. Phys Rev B 2001; 64: 234406.

Liu

Ming

. Ab initio calculations of ideal strength and phonon instability of graphene under tension. Phys Rev B 2007; 76: 064120.

Zakharchenko

Katselnson

Fasolino

. Finite temperature lattice properties of graphene beyond the quasiharmonic approximation. Phys Rev Lett 2009; 102: 046808.

Zhao

Min

Alum

. Size and dependent elastic properties of graphene nanoribbons under uniaxial tension. Nano Lett 2009; 9: 3012–3015.

Cadelano

Palla

Giordano

. Nonlinear elasticity of monolayer graphene. Phys Rev Lett 2009; 102: 235502.

Yakobson

Brabec

Berhholc

. Nanomechanics of carbon nanotubes: Instabilities beyond linear response. Phys Rev Lett 1996; 76: 2511–2514.

Huang

. Nonlinear mechanics of single-atomic-layer graphene sheets. Int J Appl Mech 2009; 1: 443–467.

Arroyo

Belytschko

. Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule. Phys Rev B 2004; 69: 115415.

Arroyo

Belytschko

. An atomistic based finite deformation membrane for single layer crystalline films. J Mech Phys Sol 2002; 50: 1941–1977.

10.

Theodosiou

Saravanos

. Molecular mechanics of graphene using finite elements. Eur J Comput Mech 2013; 22: 59–78.

11.

Tadmor

Ortiz

Phillips

. Quasicontinuum analysis of defects in solids. Phil Mag A 1996; 73: 1529–1563.

12.

Tadmor

Smith

Bernstein

. Mixed finite element and atomistic formulation for complex systems. Phys Rev B 1999; 59: 235–245.

13.

Zhang

D-B

Anatyeva

Dumitrica

. Bending ultrathin graphene at the margins of continuum mechanics. Phys Rev Lett 2011; 106: 255503.

14.

Tapastzo

Dumitrica

Kim

. Breakdown of continuum mechanics for nanometre-wavelength rippling of graphene. Nat Phys 2012; 8: 739–742.

15.

Cedra

Ravi-Chandar

Mahadevan

. Buckling cascades in free sheets. Nat 2002; 419: 579–580.

16.

Cedra

Mahadevan

. Geometry and physics of wrinkling. Phys Rev Lett 2003; 90: 074302-1.

17.

Puntel

Deseri

Fried

. Wrinkling of a stretched thin sheet. J Elast 2011; 105: 137–170.

18.

Steigmann

. A well posed finite-strain model for thin elastic sheets with bending stiffness. Math Mech Sol 2013; 18: 103–112.

19.

Ortiz

Gioia

. The morphology and folding patterns of buckling-driven thin film blisters. J Mech Phys Sol 1994; 42: 531–559.

20.

Timoshenko

Gere

. Theory of elastic stability. Mineola, NY: Dover, 2009.

21.

Steigmann

Ogden

. Classical plate buckling theory as the small-thickness limit of three dimensional incremental elasticity. ZAMM 2014; 94: 7–20.

22.

Braides

. Gamma convergence for beginners. Oxford: Oxford University Press, 2002.

23.

Friesecke

James

Muller

. A hierarchy of plate models derived from nonlinear elasticity by gamma convergence. Arch Rat Mech Anal 2006; 180: 183–236.

24.

Bhattacharya

James

. A theory of thin films of martensitic materials with application to microactuators. J Mech Phys Sol 1999; 47: 531–576.

25.

Friesecke

James

. A scheme for the passage from atomic to continuum theory for thin films, nanotubes and nanorods. J Mech Phys Sol 2000; 48: 1519–1540.

26.

Gurtin

Murdoch

. A continuum theory of elastic material surfaces. Arch Rat Mech Anal 1975; 57: 291–323.

27.

Steigmann

Ogden

. Elastic surface-substrate interactions. Proc R Soc Lond A 1999; 455: 437–474.

28.

Cohen

DeSilva

. Nonlinear theory of elastic surfaces. J Math Phys 1966; 7: 246–253.

29.

Murdoch

Cohen

. Symmetry consideration for material surfaces. Arch Rat Mech Anal 1979; 72: 61–97.

30.

Noll

. A mathematical theory of the mechanical behavior of continuous media. Arch Rat Mech Anal 1858; 2: 197–226.

31.

Chhapadia

Mohammadi

Sharma

. Curvature-dependent surface energy and implications for nanostructures. J Mech Phys Sol 2011; 59: 2103–2115.

32.

Pitteri

Zanzotto

. Continuum models for phase transition and twinning in crystals. Boca Raton, FL: Chapman and Hall, 2003.

33.

Ming

. Cauchy–Born rule and the stability of crystalline solids: Static problem. Arch Rat Mech Anal 2007; 183: 241–297.

34.

Pitteri

. Reconciliation of local and global symmetries of crystals. J Elast 1984; 14: 175–190.

35.

Ericksen

. On the symmetry of deformable crystals. Arch Rat Mech Anal 1979; 72: 1–13.

36.

Ericksen

. On the Cauchy-Born rule. Math Mech Sol 2008; 13: 199–220.

37.

Fadda

Zanzotto

. The arithmetic symmetry of monoatomic 2-nets. Acta Cryst 2000; A56: 36–48.

38.

Ericksen

. Some phase transitions in crystals. Arch Rat Mech Anal 1980; 73: 99–124.

39.

Pitteri

. On ν +1 lattices. J Elast 1985; 15: 3–25.

40.

Friesecke

Theil

. Validity and failure of the Cauchy-Born hypothesis in a two dimensional mass spring lattice. J Nonlin Sci 2002; 12: 445–478.

41.

Wang

Z-Q

Zhao

Y-P

Huang

Z-P

. The effects of surface tension on the elastic properties of nano structures. Int J Eng Sci 2010; 48: 140–150.

42.

Marsden

Hughes

TJR

. Mathematical foundations of elasticity. Mineola, NY: Dover Publications, 1994.

43.

Simo

Marsden

. On the rotated stress tensor and the material version of the Doyle-Ericksen formula. Arch Rat Mech Anal 1984; 86: 213–231.

44.

Simo

Marsden

Krishnaprasad

. The Hamiltonian structure of nonlinear elasticity: The material and convective representation of solids, rods and plates. Arch Rat Mech Anal 1987; 104: 125–183.

45.

Dresselhaus

Jorio

. Group theory: Application to the physics of condensed matter. New York, NY: Springer, 2007.

46.

Zheng

Q-S

. Theory of representation for tensor functions: A unified invariant approach to constitutive theories. Appl Mech Rev 1994; 47: 545.

47.

Zheng

Q-S

Boehler

. The description, classification, and reality of material and physical symmetries. Acta Mech 1994; 102: 73–89.

48.

Zheng

Q-S

. Two dimensional tensor function representation for all kinds of material symmetry. Proc R Soc Lond A 1997; 433: 127–138.

49.

Mohr

Papagellis

Maultzsch

. Two dimensional electronic and vibrational band structure of uniaxially strained graphene from ab initio calculations. Phys Rev B 2009; 80: 205410.