Graphene as a Prototypical Model for Two-Dimensional Continuous Mechanics

: This paper reviews a few problems where continuous-medium theory specialized to two-dimensional media provides a qualitatively correct picture of the mechanical behavior of graphene. A critical analysis of the parameters involved is given. Among other results, a simple mathematical description of a folded graphene sheet is proposed. It is also shown how the graphene–graphene adhesion interaction is related to the cleavage energy of graphite and its C 33 bulk elastic constant.


Introduction
Besides its many potential applications [1], including the improvement of polymer properties and those of other host materials [2,3], graphene is an interesting model material for testing the prediction of two-dimensional (2D) mechanics [4]. In linear elasticity, the simplest formulation of the elastic energy of a 2D medium is the surface integral involving the elements of the 2 × 2 strain tensor and the 2 × 2 curvature tensor κ. Equation (1) relies on four elastic constants: the extensional stiffness A (N/m), the in-plane Poisson coefficient ν, the bending stiffness or flexural rigidity D (J), and the "bending Poisson coefficient" ν (see Section 6). In a previous paper [5], it was shown how experimental measurements together with atomistic calculations could provide reliable values of the elastic constants (A = 365 N/m, ν = 0.20, D = 1.6 eV). These values are in line with recent measurements of in-plane elastic constants of graphene: C ≡ A(1 − ν 2 ) = 342 N/m (2D Young modulus) and ν = 0.19 [6]. The bending modulus D varies between 1.4 eV (ab-initio calculations [7]) to 1.7 eV (mechanical properties of single-wall carbon nanotubes [8], see Section 7). The eV unit used for this modulus refers to its purely atomistic origin: classically, a zero-thickness medium has no rigidity.
With the accumulation of more and more data, it becomes possible to explore the remarkable mechanical properties of graphene, revealed by microscopic manipulations. In most cases, computer modeling based on semi-empirical or density functional theory (DFT) calculations reproduces the experimental observations quite well. In some cases, continuum theory turns out to provide a reliable picture when zooming out the structure from the atomic scale up to a few tens of nanometers. Some examples are illustrated in what follows. In this paper, simple models are favored over detailed calculations. The reader interested in the latter approach may refer to the cited literature and to a recent comprehensive review on the mechanical properties of graphene [9].

Mechanical Properties
Although unstable in free space, graphene is often considered as the strongest material ever found. This assertion needs to be quantified by parameters such as breaking strength, fracture toughness, covalent bond energy, and other mechanical data.
The theoretical tensile strength of graphite in directions parallel to the basal plane is around 100 GPa. The strength of graphene can be estimated by multiplying the bulk value of graphite by the Van der Waals thickness of a single sp 2 carbon plane, to obtain 34 N/m. DFT calculations set up the ultimate tensile strength of graphene along the zigzag and armchair directions at 36 N/m and 40 N/m, respectively [10]. By comparison, the ultimate uniaxial strength of graphene extracted from nano-indentation amounts to 40 N/m [11,12]. This latter value is the result of fitting experimental data with non-linear Hooke's law γ = Cε + C 2 ε 2 , here applied to an isotropic medium subject to uniaxial strain. In that expression, γ is the applied tension (N/m), ε is the strain, C is the first-order 2D Young's modulus introduced in Section 1, and C 2 is a second-order coefficient [11]. With a negative C 2 , the stress-deformation curve has a maximum-by definition the ultimate strength, equal to −C 2 /4C 2 . This represents the largest possible tensile strength graphene can support before breaking. As reviewed in Ref. [13], the introduction of defects (boundaries between large or small grains, vacancies, slits) decreases the strength [12] from 40 down to 1 N/m [4]. By comparison, the tensile strength of a PET (polyethylene terephthalate) film of 9 µm thickness (equivalent to a stack of about 26,500 graphene planes) is 1800 N/m.
Indications exist that graphene undergoes brittle failure [10,14]. Consequently, Griffith's model of crack propagation [15] is more relevant than the theoretical maximum of the γ-ε curve to address the breaking strength of a 2D material [4,16]. Griffith's criterion states that a fracture of initial length l 0 will propagate in the medium when the tension applied perpendicular to the crack length exceeds the critical value γ c = √ 2λC/πl 0 . In that formula, λ is the edge excess energy in the direction of the crack (analog to the surface energy for a 2D material) and C is the 2D Hooke constant. The fracture toughness is by definition K c = √ 2λC. Its value measured on free-standing bilayer graphene is 2.7 × 10 −3 Nm −1/2 [14]. In few-layer graphene, K c should increase linearly with the number of atomic planes, which means that the toughness of monolayer graphene can be estimated at about half the reported value-1.3 × 10 −3 Nm −1/2 (1.2 × 10 −3 Nm −1/2 according to Ref. [4]).
From the energetics point of view, one may estimate the value of the bond breaking energy E b from the cohesive energy of graphene (7.91 eV/atom [17]). The simple formula E c = Z 2 E b with Z = 3 the number of nearest neighbors yields E b = 5.3 eV, which corresponds to a breaking energy G = 3.4 nJ/m to cut the bonds perpendicular to a zigzag direction. We note the excellent agreement with the excess energy reported for the zigzag edge of graphene (half the breaking energy G), λ = 1.7 nJ/m [18]. Using this value of λ together with C = 350 N/m [5] yields K c = 1.1 × 10 −3 Nm −1/2 , in fair agreement with the experimental value reported above for monolayer graphene.
All the above reasoning was based on the implicit hypothesis of flat graphene. Out-of-plane distortions (ripples, wrinkles, etc.) are common features of mechanically transferred graphene [19]. How shape corrugation impacts the mechanical properties and the elastic constants of a 2D material is a new and growing field of research [20][21][22].

Graphene Folding
Folding of graphene has often been observed in high-resolution electron microscopy images [23][24][25]. It has been shown by electron nano-diffraction that graphene folds preferentially along either the armchair or the zigzag direction [26]. The folded sheet forms a bilayer with continuous curved connection. The elastic energy in the curved fold is a penalty that can be more then balanced by the Van der Waals adhesion of the overlapping flat layers. When this happens, the adhesion energy acts as a driving force to increase the area of the overlap. In the absence of pinning defect, the system will evolve spontaneously until the overlapping area reaches a maximum. With a simple dimensional reasoning, one may estimate the average curvature radius of the folded edge to be R ≈ √ D/σ, where D is the flexural rigidity of graphene (1.6 eV) and σ is the interlayer binding energy. In the case of graphite, the Van der Waals cohesion can be measured accurately from the cleavage energy (the energy needed to separate two parts of graphite along a basal plane): σ Gr = 0.37 J/m 2 [27]. Experimentation is much more difficult to conduct for bilayer graphene. Electronic-structure calculations yield dispersed values that may depend strongly on the particular AB, AA, or Moiré stacking [28]. To circumvent the problem, we have used a relationship demonstrated in the Appendix within the Girifalco model [29]: σ + σ Gr = hC 33 /20, where C 33 is the elastic constant of graphite along the c axis and h is the interlayer distance. Using the experimentally-measured value of σ Gr [27] and C 33 = 38.7 GPa [30], we obtain σ = 0.28 J/m 2 for the graphene-graphene interaction (see also Table 2). The folding mechanism can be examined with a simple continuous-mechanical model [24,26,[31][32][33]. Assuming a graphene ribbon of width L that is folded back along a line perpendicular to its length, neglecting the in-plane strain energy, we simply have to minimize the difference between the bending energy and the Van der Waals energy where D is the flexural rigidity, σ is the interlayer binding energy, C is the curve drawn by the folded edge (see Figure 1), l 0 is the initial length of this curve, and κ is the local curvature. During the evolution of the folding, the curve changes in such a way that its length l = C ds decreases, thereby increasing the overlapping area and then the Van der Waals adhesion. At the same time, the bending energy increases until an equilibrium is reached. We have a variational problem to solve [25]. This problem can be further simplified by assuming an a-priori mathematical form of the curve C , evaluating the right-hand side of Equation (2) and minimizing the result with respect to the parameter(s) on which the curve depends. We have chosen a one-parameter racket-like curve whose expression in polar coordinates (ρ, ϕ) is where h is the Van der Walls spacing between two graphene planes, α is a free parameter in the interval (0, 1). The shape of the curve is shown in Figure 1b. ϕ = 0 corresponds to the apex of the folded edge where the local curvature radius is R = h/[cos(απ/2)(α 2 + 3)]. The initial length l 0 plays no role in the minimization of Equation (2), because it is a constant term. For mathematical reasons, l 0 was chosen as the length of the curve defined by Equation (3) for α = 0. Equation (3) transforms Equation (2) in a function of the curvature parameter α: where f α = cos(αϕ/2) 4 cos(απ/2) cos(ϕ/2) , the prime marks denoting derivatives with respect to ϕ. The first integral involves the expression of the curvature κ in polar coordinates. The second integral is the variation of length l 0 − l of the curve. Both integrals require a numerical evaluation. The function (Equation 4) depends on the dimensionless parameter m = D/σh 2 . The minimization of ∆E/L with respect to α leads to the results listed in Table 1 (see Appendix). Table 1. Variation of the dimensionless parameter m = D/σh 2 , folding parameter α (Equation 3), curvature radius R of the folded edge at its apex, and the bending energy per unit length w b (first term in Equation 2) versus graphene parameters (flexural rigidity D, interlayer binding energy σ, inter-sheet distance h). It is interesting to investigate how the geometry of the folded edge changes with the number of layers. Due to the inter-layer interactions, the bending rigidity of bilayer graphene is expected to be much greater than two times the bending modulus of an isolated graphene plane [34]. According to first-principle molecular dynamics calculations [35], the bending modulus D of graphene is multiplied by a factor of 60 on going from single-layer to double-layer and then scales with the cube of the number of layers, in agreement with the continuum elasticity of plates that gives a flexural rigidity proportional to the cube of the thickness. Whereas the flexural rigidity of graphite flakes with thickness above 2.4 nm has been shown to follow the expected cubic law [36], a controversy exists. In case the atomic planes slip over each other, continuum theory predicts that D for multilayer graphene is proportional to the number of planes, which has been confirmed by some atomistic simulations [37,38]. On the experimental side, atomic force microscopy (AFM) measurement of the deformation of convex-buckled suspended graphene ribbons yielded D = 35.5 eV for bilayer graphene, with relative uncertainty of 50% [39]. The last row of Table 1 corresponds to this parameter; h has been doubled because the inter-sheet distance must be measured from the plane at mid-thickness. The obtained curvature radius R at the folded edge apex also refers to the surface at mid thickness. The value of σ used for the bilayer-bilayer interaction was taken from Table 2. Table 2. Binding energy per unit area of an n-layer graphite sheet deposited on another n-layer graphite sheet for different values of n.

Self -Folding of Graphene
The previous section was devoted to graphene in free space-a situation that is not of direct interest. However, most of what has been predicted remains qualitatively correct for graphene deposited on a substrate. Of course, the shape of the folded sheet changes from that of Figure 1b to that of Figure 1a due to the underlying material [38]. The case of SiO 2 is very frequent: the driving force for folding is the difference σ gr-gr − σ gr-SiO 2 . This difference can be small and even negative according to recent accurate measurements: σ gr-SiO 2 = 0.45 J/m 2 for monolayer graphene and 0.31 J/m 2 for few-layer (2-5) graphene [40]. Accordingly, folding should be favorable for few-layer graphene, but not for monolayer graphene (see Table 2). Nanoindentation of few-layer graphene on SiO 2 by a triangular AFM tip gives rise to interesting observations. When the tip has small-amplitude horizontal movements during the punch, three graphene ribbons are peeling and growing away from the edges of the indentation [41]. The driving force is the difference between the graphene-graphene and graphene-silica adhesion energies. While the fold moves away from its origin, the ribbon gets narrower so as to reduce the bending energy of the fold. The tapering angle θ is given by the equation [41] sin where w b is the bending energy per unit width of the folded region and G is the bond breaking energy introduced in Section 2. The value of θ measured experimentally for bi-to tetra-layer graphene samples was 12 • [41]. Assuming that the bond breaking energy scales linearly with the number of layers, one obtains G ≈ 42 eV/nm for a bilayer. One therefore needs w b ≈ 4.5 eV/nm to obtain a tapering angle of 12 • for a self-folded bilayer. This value is a factor of 8.3 smaller than that given in the last row of Table 1. However, the bending constant D of the self-folded bilayer could be around 3 eV-not 35 eV-if the atomic layers were bending without coupling, as predicted by molecular dynamics [38]. This would bring w b close to the value estimated from the tapering angle. Instead of making a single fold, a monolayer graphene can scroll, but this transformation works only in solution with alcohol, which helps the graphene layer to detach from the substrate [42]. The gr-gr interlayer interaction is about a factor of two larger in a scroll than in a single fold, because there is an interaction on both faces of graphene in most part of the scroll. Then, Van der Waals interaction favors the scroll over the fold, the more so that the bending energy in a scroll is not a strong penalty compared to a folded sheet [5,43].

Mechanical Tearing
Graphene flakes obtained by mechanical exfoliation have a complex shape [44] that sometimes present triangular indentation of their edges. Such indentation can be obtained by applying a tearing load parallel to the graphene layers so as to separate the planes. The tearing force propagates initial cracks, with the consequence that ribbons are folded back during the peeling process and removed from what remains on the substrate. At the same time, the width of the peeled-off ribbons gets narrower until they reduce to nothing. The tapering angle predicted by continuum elasticity [45]-extrapolated to graphene in the case of low adhesion-is governed by the equation [35] sin with D the flexural rigidity, σ the adhesion energy of the graphene layer on its substrate, and G the bond breaking energy. The expression is different from Equation (5), chiefly because the torn ribbon is stretched by the applied force. The latter must overcome the adhesion strength σ. Both D and G depend on the number N of atomic planes, with a more or less linear dependence of G on N.
Experiments [35] reveal that θ is about 7 • for a monolayer graphene and about 20 • for a bilayer, both on SiO 2 . Using the parameter values given above for monolayer graphene on SiO 2 (σ gr-SiO 2 = 0.45 J/m 2 , D = 1.6 eV, G = 3.4 nJ/m), one obtains θ = 9 • . For bilayer graphene, G can be multiplied by a factor of two, σ gr-SiO 2 can be reduced to 0.3 J/m 2 , as reported in Section 4. Then, the observed value θ = 20 • requires D = 30 eV-quite close to the value of 35 eV measured for the bending rigidity of bilayer graphene [39]. In view of this result, mechanical tearing is likely to bend the atomic layers by keeping them coupled, by opposition with self-folding where the apparent bending modulus was found to be much smaller (Section 4).

Gaussian Curvature
The bending contribution to the elastic energy of a 2D body (Equation 1) involves a "bending Poisson coefficient" ν that enters the equation in front of the Gaussian curvature κ 11 κ 22 − κ 2 12 . The coefficient ν plays no role in nanotubes or in folded sheets, where the Gaussian curvature is zero. Alternatively, a Gaussian bending modulus D ≡ −(1 − ν )D can be introduced. It is often used in membrane theory [46], but remains difficult to measure or to evaluate numerically [47]. It has been demonstrated theoretically [48] that the coefficient ν of graphene differs from its in-plane Poisson coefficient ν. Using the parameters of second-generation Brenner potential, Davini et al. [48] arrived at ν = 0.42 and ν = 0.40, still quite similar. With the same set of parameters, the flexural rigidity obtained is D = 1.40 eV, which according the the above definition yields D = −0.8 eV [49]. By comparison, DFT calculations for carbon nanotubes and fullerene structures give D = −1.52 eV and D = 1.44 eV (ν = 0.47) [7] or D = −0.7 eV and D = 1.6 eV within the tight-binding formalism [46]. Recently, finite-temperature molecular dynamics has been performed for a graphene sheet containing a free edge in order to have access to the Gaussian bending modulus D [47]. The numerical data obtained could only be understood by attributing a small positive value to this modulus, D ≈ 0.03 eV (ν > 1). This surprising result emphasizes the dominant role of entropic effects. At zero temperature, both integrands in Equation (1) must be positive definite expressions. This condition imposes A > 0, D > 0, and both ν and ν to be the interval (−1, +1).

Collapse of Nanotubes
Under hydrostatic pressure, an isolated single-wall carbon nanotube deforms first radially, then flattens and eventually collapses [50]. For large diameter, the collapse is realized when two opposite sides of the deformed nanotube come close to the Van de Waals distance of 0.34 nm. Combining Raman spectroscopy of carbon nanotubes under hydrostatic pressure and atomistic calculations, Torres et al. [8] were able to determine the deformation and the collapse of the cylindrical structure. For each nanotube, two critical pressures were observed: the onset of ovalisation, P 1 , which for most nanotubes corresponds to a bifurcation point in the radius versus pressure curve, and the collapse pressure P 2 . The authors of this study have shown that the two pressures are both related to a universal function where d is the nanotube diameter at 1 atm, D is the flexural rigidity of graphene, and β is an empirical parameter. The relations between P 1 , P 2 , and P c (d) are P 1 = P c (d) and P 2 ≈ 1.5P c (d). By fitting all the available data, Torres et al. arrived at D = 1.7 ± 0.2 eV and β = 0.44 nm. According to Equation (7), the critical pressure P c (d) increases with increasing d, reaches a maximum at d = √ 5/3β ≈ 0.57 nm, and then decreases. Asymptotically, P c (d) ∼ 24D/d 3 , which is what continuous-mechanics predicts for a thin-wall elastic ring (Lévy-Carrier law) [51].

Conclusions and Perspectives
Continuous mechanics specialized to two-dimensional bodies seems to work qualitatively well for graphene, and most likely for related materials like boron nitride.. If the picture that emerges from this approach is reasonably coherent with the observation, experimental measurements and computed values of the relevant parameters may vary substantially. The quality of graphene [52], the composition and corrugation of the substrate [53], the temperature [54], etc. all play a role in the mechanical behavior of graphene. The presence of wrinkles or ripples [19] introduces an additional complexity that may affect the picture [55]. In particular, thermal fluctuations in suspended graphene, in addition to requiring a statistical approach, may amplify non-linear effects [56].
Among the linear-elastic constants, the flexural rigidity of graphene and especially the Gaussian modulus remains the most difficult to measure mechanically. The situation is worse for few-layer graphene, where it is difficult to appreciate to what degree the bent atomic layers are coupled. In the domain of mechanics, a full understanding of how graphene fails under tensile or bending strain is another important topic [14,57]. Works devoted to these problems are on the increase, and there is no doubt that important progress will be made in the coming months.
Acknowledgments: Helpful discussions with Dr. Evgeni Ivanov and Dr. Fatemeh Ahmadpoor during the writing of the paper are greatly acknowledged. This work has been partly supported by the H2020 MCA RISE project "Graphene-3d", Grant No. 734164, funded by the European Union.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Van der Waals Binding Energy between Few-Layer Graphite Sheets
Consistent with the continuous-medium theory used to describe the elastic properties of graphene, the Girifalco model [29] is an interesting approach for addressing the Van der Waals cohesion of layered materials. The central idea of the model is to treat each atomic layer as a continuous plane containing N atoms per unit area with which an external atom interacts through a Lennard-Jones 6-12 potential. If d is the distance, the interaction energy between the plane and the atom is V(d) = b/10d 10 − a/4d 4 , where a and b are two positive constants. Considering the external atom as belonging to a graphene plane located at a distance h above another graphene plane, it follows from the expression of V(d) that the binding energy per unit area between the two parallel planes is where the index 1, 1 means the interaction between a monolayer of graphene and another monolayer of graphene. The same procedure can be developed for the interaction energy between an n-layer graphite sheet and another n-layer graphite, assuming a regular spacing h between the successive n + n atomic planes. Summing up the interactions over the 2n − 1 inter-plane distances yields the following expression of the interaction energy At the limit n → ∞, the Van der Waals binding energy of graphite is obtained: with ζ k the Riemann zeta function of power k. Using the parameters a, b, h, and N given in Ref. [5], the values listed in Table 2 are obtained. The result for graphite agrees remarkably well with available experimental data [27,58]. Equations (8) and (10) readily lead to the relation with C 33 the elastic constant of graphite along the c axis, where A = aζ 4 , B = bζ 10 , knowing that B/A = h 6 and C 33 = 6N A/h 5 [5].