Abstract

We consider how the bending stiffness of a multilayer graphene sheet relies on its bending geometry, including the in-plane length L and the curvature κ. We use an interlayer shear model to characterize the periodic interlayer tractions due to the lattice structure. The bending stiffness for the sheet bent along a cylindrical surface is extracted via an energetic consideration. Our discussion mainly focuses on trilayer sheets, particularly the complex geometry-dependency of their interlayer stress transfer behavior and the overall bending stiffness. We find that L and κ dominate the bending stiffness, respectively, in different stable regions. These results show good quantitative agreement with recent experiments where the stiffness was found to be a non-monotonic function of the bending angle (i.e., Lκ). Besides, for a given in-plane length, the trilayer graphene in the flat state (κ → 0) is found to have the maximum bending stiffness. According to our analytical solution to the flat state, the bending stiffness of trilayer graphene sheet can vary by two orders of magnitude. Furthermore, once multilayer graphene sheets are bent along a cylindrical surface with small curvature, the sheets perform similar characteristics. Though the discussion mainly focuses on the trilayer graphene, the theoretical framework presented here can be readily extended for various van der Waals materials beyond graphene of arbitrary layer numbers.

1 Introduction

Two-dimensional (2D) materials such as graphene are a class of atomically thin sheets that have promised extraordinary electrical, thermal, and mechanical properties [13]. Multilayer 2D materials, also known as van der Waals (vdW) materials, have been recently developed by vertically stacking individual 2D materials, which have provided a useful routine for the unique properties of these layers being well-exploited in electronic devices [46]. The vdW materials are thin and very susceptible to bending deformations, just like a few pieces of paper [7,8]. As a result, they have shown various out-of-plane deformation forms in applications, such as ripples, buckles, scrolls, folds, bubbles, tents, and so on [914]. Though sometimes undesired for device performances, these forms have found exciting applications for mechanics metrology, novel strained electronics, and in situ transmission electron imaging [15,16] as well as enabled a number of anomalous phenomena, including flexoelectricity, pseudomagnetic fields, and photon emission [1719]. In these scenarios, bending stiffness of the vdW materials plays a crucial role in controlling out-of-plane geometry as well as the associated functionalities.

As for multilayer sheets, on the assumption that the layers are perfectly bonded without interlayer slip (i.e., plane sections remain plane upon bending), bending stiffness would have a cubic relationship with thickness and be linearly proportional to Young’s modulus [20,21]. However, multilayer vdW materials feature strong in-plane covalent bonding and weak interlayer vdW interactions [21,22], so that layers are readily subjected to relative slips, as evidenced by recent observations on the structural super-lubrication, self-tearing, and self-rotation [2325]. Such slip would make the bending stiffness of vdW materials significantly overestimated by traditional plate theory [21,2629]. More intriguingly, recent reports have revealed that the bending stiffness of multilayer graphene sheets is a variable depending on their bending geometry, rather than a constant material property [27,29]. On the one hand, for multilayer graphene sheets with prescribed layer numbers, the theoretical prediction showed that the bending stiffness increases with the increasing in-plane length [29]. On the other hand, experiments and density functional theory (DFT) calculations indicated that the bending angle has a decisive effect on the bending stiffness [27]. These findings not only support the invalid cubic relationship between bending stiffness and thickness claimed in recent reports [21,2729] but also may explain the discrepancies among different measurements on even the same type of vdW materials in these reports [8,27,3034].

These advances have posed a question about how the bending stiffness of a vdW material depends on its two geometric descriptors—the in-plane length L and the bending angle θ. Alternatively, to characterize the bending geometry, we may use the in-plane and out-of-plane length scale: the length L and the radius of bending curvature 1/κ (since θ = ). The essence of this question is how much the bending geometry of a vdW material modifies the classical “plane sections remain plane” assumption via interlayer shear. Apparently, besides the bending geometry, modifications are expected due to the intrinsic material properties, including the bending stiffness of a single layer, the lattice constant, the in-plane stiffness, the interlayer spacing, and the interlayer shear traction-separation law. This has been implied in a recent experimental work where among graphene, hexagonal boron nitride (hBN), and molybdenum disulfide (MoS2) multilayers with comparable thickness, the graphene system shows the largest Young’s modulus but the smallest bending stiffness because of its weak interlayer shear resistance [21]. In this work, we adopt the material properties of graphene sheets (which have been relatively well-studied) and investigate the collective effect of L and κ on their bending stiffness that has been unclear previously. However, one may find that the proposed theoretical framework applies to various vdW materials, where the application can be achieved simply by the substitution of specific material parameters.

The rest of the paper is organized as follows. In Sec. 2, we analyze the free body diagram of each layer in a multilayer, which is bent along a cylindrical surface. We assume fixed interlayer distance and constant amplitude for the interlayer shear traction as a sinusoidal function of the interlayer slip. The equilibrium equations for these layers are derived, and the effective bending stiffness is expressed based on the stored energy. In Sec. 3, we discuss the stability of interlayer slip for trilayer graphene as an exemplary case. In Sec. 4, we predict specific values of bending stiffness for trilayer graphene sheets and compare these with experimental measurements. A more detailed bending stiffness map is also provided based on the two geometric factors. We further explore the multilayer graphene of three or more layers in Sec. 5. Finally, the concluding remarks are presented in Sec. 6.

2 A Bending Model for Multilayer Graphene Sheets

When a multilayer sample is bent along a cylindrical surface, and the interlayer is perfectly bonded without slip, the deflection curve will only occur in one plane, just like a pure bending of a beam with a rectangular cross section. There exist three characteristics or hypotheses [35]: (1) The neutral surface exists, in which longitudinal fibers do not change in length, and the neutral axis (longitudinal) passes through the centroid of the cross-sectional area. As for the straight sample with a rectangular cross section, the neutral surface is the middle plane. (2) All cross sections of a straight beam remain plane and perpendicular to the longitudinal axis during deformation. (3) Any deformation of the cross section within its own plane can be neglected. However, if the interlayer action is limited, the second hypothesis above will no longer hold [21]. For example, when we bend a book along a cylindrical surface, whose interlayer may be deemed as ideal lubrication with no resistance to slip, the cross section at the free end which is a plane before deformation becomes a surface and not perpendicular to the longitudinal axis (as seen in more detailed descriptions as the  Appendix). As for a multilayer graphene sheet whose interlayer is neither perfectly bonded nor ideally lubricated, the interlayer shear and slip would coexist, as shown in Fig. 1(a). Based on the observation of the multilayer sample bending along a cylindrical surface, we have the following four hypotheses in order to simplify the problem:

  1. The material is homogeneous and linear elastic.

  2. The middle plane of the whole multilayer sample is a neutral surface.

  3. Any deformation of the cross section within its own plane, including the change of interlayer distance, can be neglected.

  4. The interlayer shear response is isotropic, represented by the response in the zigzag direction.

Fig. 1
(a) Bright-field STEM image of 12-layer graphene bent to 12 deg, where the number of atomic columns in the arc is higher for the outer graphene layers than for the inner layers owing to interlayer slip, implying that the cross section is no longer a plane. (a) is adapted with permission from the literature [27]. (b) Schematic diagram of the multilayer sample bent along a cylindrical surface without interlayer sliding, where the part of the longitudinal section of the sample before and after bending deformation are shown respectively.
Fig. 1
(a) Bright-field STEM image of 12-layer graphene bent to 12 deg, where the number of atomic columns in the arc is higher for the outer graphene layers than for the inner layers owing to interlayer slip, implying that the cross section is no longer a plane. (a) is adapted with permission from the literature [27]. (b) Schematic diagram of the multilayer sample bent along a cylindrical surface without interlayer sliding, where the part of the longitudinal section of the sample before and after bending deformation are shown respectively.
Close modal

On account of our third hypothesis, we only take interlayer shearing, one of the interlayer failure modes, into consideration but ignore two other failure modes, rippling and kink buckling [29].

Besides, we only discuss the condition of odd layers, though the scenarios of even layers may be analyzed in a similar way. Suppose that the total layer number is N = 2n + 1, n ∈ N+. We denote the middle layer by i = 0, and z0 denotes the corresponding thickness direction coordinate value of the middle plane. From the middle layer to the curved convex side, i gradually increases, and to the concave side, i decreases. Thus, i belongs to {− n, …, −1, 0, 1, …, n}, and the corresponding thickness direction coordinate value is zi = id, where d denotes the interlayer distance and equals 0.34 nm for multilayer graphene.

In most experiments, the thickness of multilayer graphene sheets is small relative to its in-plane size. Thus, when the sample is bent along a cylindrical surface, we treat it as a multilayer beam. Besides, for the symmetry, we only take the half into consideration, and the symmetrical cross section without relative interlayer displacement can be treated as one end clamped.

We then consider the case of the perfectly bonded interlayer, where the second hypothesis in traditional theory still holds, as shown in Fig. 1(b). The length of the longitudinal axis in ith layer changes from L to L(1 + zi/R), scilicet L(1 + ziκ), where R is the curvature radius of the middle surface after bending, and κ = 1/R is the corresponding curvature.

Having considered the two limiting cases, we move on to the real bending stiffness of a multilayer sheet. To do so, we consider the energy stored in the multilayer sheet when it is bent along a cylindrical surface. There are various ways to enforce this specific displacement field for each layer in the sheet, including sandwiching the sheet between two rigid concentric circles, applying distributed transverse pressure, and so on. Here, however, we use an alternative way that the middle layer is mapped to a cylindrical surface without change in the in-plane length [29], by which the problem becomes to solve ordinary differential equations (ODEs) for merely the in-plane displacement field of each layer. The out-of-plane displacement fields would be prescribed by the curvature of the cylindrical surface and our hypothesis (3).

We first apply a linearly distributed load and necessary moments to the end of 2n + 1 individual, separated layers. Such loads guarantee that the assembly of these layers into a multilayer behaves exactly the same as the perfectly bonded sheet with a curvature of κ for the middle layer. We then confine the assembled multilayer between two rigid, frictionless concentric circles with the inner radius of (1/κnd), and the outer radius of (1/κ + nd). Finally, we convert the interlayer from the perfectly bonded case into realistic cases simply by applying a relaxation load Pi at the end of each layer
Pi=E2DLL(1+κzi)L(1+κzi)=E2Dκzi1+κzi
(1)
where E2D is 2D in-plane stiffness of graphene, which is taken as 340 N · m−1 [3]. The use of the relaxation load cancels out the mechanical response (stretching or compression) originating from the initially applied linearly distributed load. Note that when the problem is proceeded using the initial undeformed configuration, the magnitude of relaxation load just equals to that of the initially applied load. However, we will use the deformed configuration assuming no interlayer slip as the reference, so that the denominator in Eq. (1) takes (1 + κzi). This would stretch/compress the sheet into a length of L if the interlayer is ideally lubricated. This process ensures that the strain energy stored in the multilayer sheet purely results from the bending of the sheet onto a cylinder. Note that we do not apply relaxation moment at the free end or relaxation pressure on the surface because these terms would not contribute to the stored energy in the deformed multilayer sheet under the assumption of rigid confinement and our hypothesis (3). In other words, the out-of-plane displacements (work conjugate of pressure), including the bending angle at the edge (work conjugate of edge moment), have been prescribed.
Referring to the results of the interlayer shear response based on molecular dynamics (MD) simulations, we choose a sinusoidal function to describe the relationship between interlayer shear traction τ and relative displacement [29,36]. The sinusoidal function reflects the periodicity of the lattice, which is similar to the classic Peierls–Nabarro dislocation model [37]. For A–B stacking graphene sheets, according to our fourth assumption, we consider a case of one layer sliding on the other static layer and describe the interlayer shear traction as
τ=τmaxsin(2πl0u)
(2)
where u is the displacement of the sliding layer, l0=3lcc is the periodical length where lc−c is the C–C bond length and taken as 0.14 nm, and the amplitude τmax can be deduced from molecular dynamics (MD) simulations on interlayer shear response [29]. τmax=Fmaxsydc, where Fmaxsy is the maximum average shear force per atom, and dc=4/(33lcc2)39.3Atom/nm2 [29]. According to Fig. 11(c) in the literature [29], we can approximatively obtain Fmaxsy=0.015nN/atom. Thus, τmax ≅ 590 MPa. Besides, we ignore the edge effect, which should be trivial for typical experimental samples whose bending length is much larger than l0 [29].

Equation (2), along with the line of our fourth hypothesis, actually dictates that samples are bent along the zigzag direction. We shall show shortly that this simple form allows for a number of explicit understandings on the important role played by the bending geometry in the bending stiffness of a multilayer sheet. The consideration of an arbitrary bending direction, however, requires a more complicated interlayer shear traction law, which should present preciser predictions but may preclude any analytical efforts. We thus leave this in a future study.

In the following modeling process, when considering the state after bending with imaginary perfectly bonded interlayer as the reference configuration, we add a bar above the symbol of all quantities, such as u¯, except for those that do not change with the coordinate system transformation. In particular, the in-plane displacement of the middle plane of the whole sample is always zero. However, for the sake of unity, it is still denoted by u¯0.

Since the effect the normal cohesion is neglected, the equilibrium of in-plane stress can be simplified to an one-dimensional (1D) problem, as shown in Fig. 2. It is convenient for notations to introduce the sign function
sign(i)={1i>00i=01i<0
(3)
Fig. 2
Schematic diagrams of the force equilibrium in the tangential direction of a layer on the curved convex (i = 1) (a) and concave (i = −1) (b) side in a multilayer sample bent along a cylindrical surface
Fig. 2
Schematic diagrams of the force equilibrium in the tangential direction of a layer on the curved convex (i = 1) (a) and concave (i = −1) (b) side in a multilayer sample bent along a cylindrical surface
Close modal
The interlayer shear traction can be expressed as
τ¯i={τmax{sin[2πl0(u¯i1+κziu¯isign(i)1+κzi1)]}i{n,n}τmax{sin[2πl0(u¯i1+κziu¯i11+κzi1)]+sin[2πl0(u¯i1+κziu¯i+11+κzi+1)]}i{n,0,n}
(4)
where the change of periodical length l0 has been considered to match the geometry of deformed lattice [29].
Accordingly, the equilibrium equation can be written as
dσ¯idx¯iτ¯i=0,
(5)
where the symbol σ refers to the stress resultant rather than stress. According to the geometrical relationship on the assumption of small deformation, the strain is
ε¯i=du¯idx¯i
(6)
and Hooke's law in one dimension is
σ¯i=E2Dε¯i
(7)
We then obtain the equilibrium equation in terms of displacement
d2u¯idx¯i2τ¯iE2D=0
(8)
subject to the boundary conditions
{u¯i(0)=0E2Ddu¯idx¯i|x¯=L(1+κzi)=Pi
(9)
In order to make the range of independent quantities in differential equations consistent, we map the coordinate system by introducing
x~=x¯i1+κzi[0,L]
(10)

For the unity of symbols, we also introduce u~i(x~)=u¯i(x¯i) and τ~i=τ¯i.

By this, the governing equation, Eq. (8), becomes
d2u~idx~2(1+κzi)2τ~iE2D=0
(11)
and boundary conditions, Eq. (9), become
{u~i(0)=0E2D1+κzidu~idx~|x~=L=Pi
(12)
The displacement solution
u~i=u~i(x~)
(13)
can be obtained either analytically or numerically.
We then take the state before any bending deformation as the reference configuration, and remap the coordinate system, introducing
x=x~=x¯i1+κzi[0,L]
(14)
The displacement solution becomes
ui(x)=u~i(x~)+κzix~
(15)
The corresponding strain solution is
εi(x)=duidx
(16)
We now can compute the in-plane strain energy of an individual layer by
Vεi=12E2D0Lεi2dx
(17)
and the bending energy of the layer by
Vκi=12B(κ1+κzi)2L
(18)
where B is the bending stiffness of monolayer graphene and taken as 1.44 eV in our theoretical analyses [38]. In particular, for the middle layer
{Vε0=0Vκ0=12Bκ2L
(19)
For multilayer graphene, the deformation energy is the sum of in-plane strain energy and bending energy of each layer. The total energy includes the deformation energy and the cohesion energy that results from the change in interlayer distance. On account of our third assumption, the cohesive energy can be neglected; that is, the total energy is equal to the total deformation energy
V=i=nn(Vεi+Vκi)
(20)
Therefore, the effective bending stiffness of a multilayer graphene bent onto a cylinder can be extracted
Deff=2Vκ2L
(21)
Here, we note that the bending stiffness defined by Eq. (21) is a nominal quantity, which is consistent with the definition of bending energy in experimental measurement and DFT calculations [27]. In Secs. 4 and 5.1.2, we will compare our results with the experimental measurement and DFT calculations.

3 Stability Analysis of Interlayer Slip Within a Trilayer Graphene Sheet

The trilayer sheet is a good starting point because of not only its simple formulation but also that the associated discussion could provide a basis for sheets of more than three layers. In this section, we, therefore, analyze the interfacial shear stress transfer (in particular, the stability of interlayer slip) within a trilayer graphene sheet.

3.1 Qualitative Theory.

Stability analysis is performed with the first layer on the curved concave side as a representative. From Eqs. (2) and (8), we can obtain the governing equation
d2u¯1dx¯12τmaxE2Dsin(2π(1κd)l0u¯1)=0
(22)
Introduce ls = (1 − κd)l0/(2π), and dimensionless quantities, ϑ=u¯1/ls, η=x¯1/ls, and ω1=(lsτmax)/E2D. Then, the governing equation can be written as
d2ϑdη2ω12sinϑ=0
(23)
We introduce T = ω1η so that the governing equation becomes
d2ϑdT2sinϑ=0
(24)
or
{dϑdT=αdαdT=sinϑ
(25)
With the mapping, Tp = T, ϑp=ϑ+π, and αp=dϑp/dTp=α, Eq. (23) is able to transformed into the equation of nonlinear pendulum [39], which may help us to understand our differential system. In the classical pendulum system, Tp, ϑp, and αp correspond to the time, the angle of pendulum, and the angular velocity, and ω12 corresponds to the ratio of the acceleration of gravity to the length of the frictionless pivoted massless rod. A particular feature of Eq. (24) is the critical points at (kπ, 0), kZ, which are located at the center when k is odd and the saddle otherwise.
Integrating Eq. (24) once, we get
12(dϑdT)2+cosϑ=C0,aconstantforT[0,Tb]
(26)
where Tb = ω1L(1 − κd)/ls, and the subscript ‘b’ refers to the free end. C0 is a constant for a specific trajectory in the phase diagram, which corresponds to the total energy of the pendulum, including the kinetic and potential energy.

As shown in Fig. 2, where the reference configuration is temporarily supposed as a stress-free state, the relaxation load is applied as an active external load depending on the bending curvature in the model. Thus, the sign of the stress at the clamped end should be the same as that of the relaxation load despite the presence of interlayer shear force. In addition, the overall deformation response of the studied layer should be consistent with the relaxation load. These two judgments can also be obtained from the phase diagram of the differential equation, as shown in Fig. 3(a). For example, for the case of the curved concave side we consider, where the state after bending without interlayer slip is used as the reference configuration and the relaxation load is a tensile load, the overall displacement and strain should be positive. This is equivalent to that any place on the concave side of the trilayer sample except the clamped end correspondingly locates in the first quadrant of the phase diagram, where ϑ and α are corresponding to the displacement and strain, respectively. Thus, for the case of the curved concave side, C0 > 02/2 + cos 0 = 1. In the pendulum system, this critical C0 = 1 corresponds to the maximum possible potential energy, which is obtained from the highest point of the system. When the energy is larger than this critical value, the pendulum can swing over; otherwise, it can only swing back and forth [39]. In our system, the clamped end of the sample exists. In the phase diagram, as illustrated in Figs. 2(b) and 2(c), the boundary condition at the clamped end corresponds to the axis, ϑa=ϑ|T=0=0, where “a” refers to the clamped end. This corresponds to that the pendulum is able to swing over, and the total energy must over the maximum possible potential energy. In this way, we can also obtain the result C0 > 1.

Fig. 3
(a) The phase field (arrows) and development path (colored solid lines) derived from several points (colored solid points) with increasing T. Phase diagrams with (b) Tb = 1.5 and (c) Tb = 2.5 respectively, where the colored solid lines represent the development path with different specific ϑ at the free end. The black dashed line in (b) and (c) connects the free end of each colored line and shows the relationship between α and ϑ at the free end. From the bottom to the top, the colored solid lines in (b) and (c) are based on ϑb from 1, 2, 3, 4, 5, 6, and 7, respectively.
Fig. 3
(a) The phase field (arrows) and development path (colored solid lines) derived from several points (colored solid points) with increasing T. Phase diagrams with (b) Tb = 1.5 and (c) Tb = 2.5 respectively, where the colored solid lines represent the development path with different specific ϑ at the free end. The black dashed line in (b) and (c) connects the free end of each colored line and shows the relationship between α and ϑ at the free end. From the bottom to the top, the colored solid lines in (b) and (c) are based on ϑb from 1, 2, 3, 4, 5, 6, and 7, respectively.
Close modal
The discussed curved concave side corresponds to the first quadrant in the phase diagram in Fig. 3(a) (as illustrated by Figs. 3(b) and 3(c)), which implies that dϑ/dT=α>0. Thus, ϑ monotonically increases as T increases. The effect of curvature may first appear in the relaxation stress and then the boundary condition at the free end. The second formula of Eq. (9) further determines
αb=α(Tb)=dϑdT|T=Tb=κdω1(1κd)
(27)

The combination of the boundary condition at the clamped end, ϑa=0, and Eq. (27) implies that different from the pendulum system with initial boundary conditions, our system has moving boundary conditions.

In the phase diagram, two different trajectories cannot intersect. We take two arbitrary trajectories with Tb1 = Tb2 but ϑb1>ϑb2, where the subscripts “1” and “2” refer to these two trajectories. As shown in Figs. 3(b) and 3(c), αa increases with increasing ϑb, i.e., αa1 > αa2. Because these two trajectories cannot intersect, α1|ϑ=ϑb2>α2|ϑ=ϑb2=αb2. According to Eq. (25), dα/dϑ=sinϑ/α0 when ϑπ, while dα/dϑ=0 when ϑ=π. In the process of mechanical deformation, κ is the direct variable. Correspondingly, we use αb for analyses. ϑb=π corresponds to αb ≥ 2. Thus, when αb ≤ 2 (and ϑπ), dα/dϑ>0 (i.e., αb1>α1|ϑ=ϑb2>αb2) can be ensured. This implies that αb increases monotonically with the increasing ϑb. We suppose that the curvature increases from zero. Thus, when αb increases from zero but lower than 2 and ϑb(0,π), the system is stable without the change of monotonicity. Introduce another set of dimensional quantities, x^=x¯1/l0, L^=L/l0, u^=u¯1/l0, κ^=κd, then αb = 2 corresponds to κ^=κ^c0.015976. Besides, when Tb is less than a certain critical value, αb becomes a monotonically increasing function of ϑb in the first quadrant of the entire phase diagram, as shown in Fig. 3(b). Finding this critical value of Tb is an inverse problem, and we will obtain this value analytically in Sec. 3.2. It is the possible transformation of monotonicity that requires us to use ϑb rather than αb in boundary conditions to obtain the results in Figs. 3(b) and 3(c). When the function of αb and ϑb is monotonic, it is reasonable to choose either αb or ϑb as a boundary condition, while the choice of αb may be more straightforward.

Furthermore, we can obtain the distribution of the relative slip displacement between layers through the phase diagram, as shown by the colored solid line in Figs. 3(b) and 3(c). Although the displacement of the interlayer slip first increases monotonically from the clamped end, we may observe the phenomenon of oscillation when the length and curvature are large enough to ensure that ϑ stably reaches π, which we attribute to the sinusoidal relationship between interlayer shear traction and interlayer relative slip displacement. This oscillation phenomenon does not exist when the interlayer shear constitutive relationship is elastic or elastoplastic [40].

Note that two crucial dimensionless parameters behind the dimensionless system, Eqs. (24) or (25) are ls and ω1, which consist of the factors on geometry (κ and L) and material properties (l0, d, E2D, and τmax). This work focuses on the effect of geometry on the bending stiffness.

3.2 Critical Bending Length.

For i = −1, the boundary condition, corresponding to the governing equation, Eq. (22), is
{u¯1(0)=0du¯1dx¯1|x¯1=L(1κd)=κd1κd
(28)
We add the pending boundary condition, u¯1b=u¯1(L(1κd)) and make the boundary condition dimensionless (which corresponds to Eq. (23)) to obtain
{ϑ(0)=0ϑ(Lls(1κd))=ϑbdϑdη|η=Lls(1κd)=κd1κd
(29)
Equation (26) implies
C0=(12(dϑdT)2+cosϑ)|T=ω1Lls(1κd)=12(κdω1(1κd))2+cosϑb
(30)
From the analyses in Sec. 3.1, we know that α > 0 (i.e., dϑ/dη>0). Combining Eqs. (26) and (30), we have
dϑdη=dϑd(T/ω1)=(κd1κd)2+2ω12(cosϑbcosϑ)
(31)
or
dη=dϑ(κd1κd)2+2ω12(cosϑbcosϑ)
(32)
where η changes from 0 to L(1 − κd)/ls while accordingly ϑ changes from 0 to ϑb. We take η = 0 as the initial state and integrate Eq. (32), obtaining
η=2F(ϑ2|4ω12(κd1κd)24ω12sin2(ϑb2))(κd1κd)24ω12sin2(ϑb2)
(33)
where F is the first kind of incomplete elliptic integral, which is defined as F(φ|m)=0φ(1msin2ζ)12dζ [41]. At the clamped end
dϑdη|η=0=(κd1κd)2+2ω12(cosϑb1)=(κd1κd)24ω12sin2(ϑb2)
(34)
We introduce β=dϑ/dη and βa=(dϑ/dη)|η=0 so that Eq. (33) can be converted into
η=2F(ϑ2|4ω12βa2)βa
(35)
At the free end, Eq. (35) becomes
Lls(1κd)=2F(ϑb2|4ω12βa2)βa
(36)
Here, we introduce the implicit function about ϑb and κ
G(ϑb,κ)=F(ϑb2|4ω12βa2)L(1κd)βa2ls=0
(37)
We can obtain the inflection point, (ϑb*,αb*), by solving
d2ϑbdαb2=0
(38)
The critical value of Tb discussed at the end of Sec. 3.1 corresponds to the case that the inflection point discussed previously coincides with the stationary point, so let
(dϑbdαb)*=0
(39)
to obtain the corresponding sample length L, from which we may further deduce the corresponding dimensionless quantities introduced before, ϑb*4.188, αb*2.320, Tb*1.994, and L^*39.14. Besides, we may consider the physical implication of these parameters by taking L^* as 39.

According to the analyses on stability in Sec. 3.1, when Tb<Tb*, that is, L^L^*, ϑb and αb are mapped one on one, while when a non-monotonic function relationship between ϑb and αb appears, a bifurcation of saddle point appears, which corresponds to a snap-through sliding in terms of mechanical response [29].

Now, it is clear that when L^L^c=39 or κ^κ^c=0.015976, the interlayer slip behavior is stable. For the convenience of discussion, the phase diagram is roughly divided into four regions ([0,κ^c]×[0,L^c], [0,κ^c]×(L^c,+), (κ^c,+)×[0,L^c], (κ^c,+)×(L^c,+)), and the bifurcation phenomenon can only arise in the fourth region ((κ^c,+)×(L^c,+)), as shown in Fig. 4. Because the typical experimental measurement results are in stable regions, which will be shown in Sec. 4, we focus on the stable regions. It should be pointed out that there is a small stable part in the fourth region. However, in order to briefly discuss the continuous shear behavior, i.e., the stable part, we do not discuss this stable but small region.

Fig. 4
Region division of phase diagram based on stability. Here, we have used the dimensionless curvature κ^=κd and bending length L^=L/l0.
Fig. 4
Region division of phase diagram based on stability. Here, we have used the dimensionless curvature κ^=κd and bending length L^=L/l0.
Close modal

In the stable regions, ϑ=ϑ(η,κ) can be obtained from the implicit function, Eq. (37), which is equivalent to that the displacement solution u = u(x; κ, L) is obtained. Then, following the process of Sec. 2.2, we can deduce the effective bending stiffness Deff = Deff(κ, L), which implies that the bending stiffness is dominated by both curvature and bending length.

4 Theoretical Predictions and Comparison With Experimental Results

The theoretical analysis in Sec. 3.2 has suggested that the bending stiffness is dominated by both curvature and bending length. As for the case that the sample is bent along a cylindrical surface, the product of the curvature and the bending length is the bending angle. A question arises about whether it is possible that the bending stiffness is controlled solely by the bending angle, as previously concluded from the general trend of DFT and experimental results [27].

We try to answer the question by using our model to obtain the effective bending stiffness of trilayer graphene. This model characterizes the bending behavior of a multilayer with a fixed curvature while the experiment is slightly different as it presents the bending near a step, as shown in Fig. 5(a) [27]. In the literature [27], the model to deduce the bending stiffness assumes that there is a detached part of length lp, which includes a straight part and two curved parts with curvature radius Rp and bending angle θ (Fig. 5(b)). Besides, Hp denotes the height of the step and Lp denotes the length of the projection of the detached part onto the hBN substrate. Hp, Rp, and θ have been viewed as three independent quantities. As a result, lp and Lp can be, respectively, described as
lp=2Rpθ+Hp2Rp+2Rpcosθsinθ
(40)
and
Lp=2Rpsinθ+Hp2Rp+2Rpcosθsinθcosθ
(41)
Fig. 5
(a) One scanning transmission electron microscope (STEM) image of trilayer graphene (N = 3) over a bilayer hBN step (H = 2). (b) Schematic diagram of multilayer graphene over the hBN step. (a) and (b) are adapted with permission from the literature [27]. In order to consider the effect of the straight part on the bending stiffness, we use the adaptation model of (b) as shown in (c), which combines one curved part (LI = Rpθ) with half straight detached part (LII), and the further considered straight part on hBN (3LII). Thus, the average curvature in (c) is (LI + 4LII)/θ.
Fig. 5
(a) One scanning transmission electron microscope (STEM) image of trilayer graphene (N = 3) over a bilayer hBN step (H = 2). (b) Schematic diagram of multilayer graphene over the hBN step. (a) and (b) are adapted with permission from the literature [27]. In order to consider the effect of the straight part on the bending stiffness, we use the adaptation model of (b) as shown in (c), which combines one curved part (LI = Rpθ) with half straight detached part (LII), and the further considered straight part on hBN (3LII). Thus, the average curvature in (c) is (LI + 4LII)/θ.
Close modal

Obviously, the curved part in the experiment does not have two free ends (or the case that one end is clamped, and the other is free), but have adjacent straight parts. We first denote the length of one curved part as LI, which is equal to Rpθ, and the half-length of the straight detached part as LII. Then, further taking effect from the straight part on hBN into consideration, we consider an additional three times LII and average the curvature, as shown in Fig. 5(c). Due to the symmetry, the bending length in experiments is adopted as LI + 4LII in our theoretical analyses, and the corresponding curvature is (LI + 4LII)/θ.

We compare the theoretical prediction with the experimental measurement in Fig. 6(a). Our model agrees with the experimental results well and reveals the phenomenon that it is not necessary that a smaller bending angle corresponds to a smaller stiffness, which may conflict with the whole trend implied by other experimental data [27]. For example, as shown in Fig. 6(a), the sample with the smallest bending angle (12deg) does not have the largest bending stiffness. We attribute this novel phenomenon to the effect of bending length. To further elucidate the mechanism behind this phenomenon, we consider both curvature and bending length to deduce the bending stiffness based on our model, as shown in Fig. 6(b). We also mark the corresponding experimental points with white cross-shaped symbols and plot contour of bending angle with white dashed curves. For a specific bending length, bending stiffness indeed decreases initially with the increase of curvature (or bending angle), but may also experience a slight increase when the bending angle is about 30deg40deg. Along the contour curve of 20deg, we observe that for a specific bending angle, the bending stiffness increases solely with the increase of bending length. Moreover, as for the sample with the smallest bending angle (12deg), we reveal that it is the smallest bending length that makes the corresponding bending stiffness smaller than that of the sample, which owns a larger bending angle (16deg). Further theoretical analyses on the bending stiffness of trilayer graphene based on our model will be shown in Sec. 5.1.

Fig. 6
(a) Experimental measurement results and corresponding theoretical prediction results, whose L^ and θ have been adapted to our model. The dashed lines in (a) connect the corresponding points, implying the trend. (b) Contour plot for the bending stiffness of trilayer graphene as a function of two governing parameters: κ^ and L^. The points which correspond to the samples in (a) are marked with white crosses in (b). Contour curves of 14 deg, 20 deg, 30 deg, and 40 deg are also plotted in (b). The experimental data are from the literature [27].
Fig. 6
(a) Experimental measurement results and corresponding theoretical prediction results, whose L^ and θ have been adapted to our model. The dashed lines in (a) connect the corresponding points, implying the trend. (b) Contour plot for the bending stiffness of trilayer graphene as a function of two governing parameters: κ^ and L^. The points which correspond to the samples in (a) are marked with white crosses in (b). Contour curves of 14 deg, 20 deg, 30 deg, and 40 deg are also plotted in (b). The experimental data are from the literature [27].
Close modal

In our analyses, an additional 3LII is included to consider the effect from the straight part on hBN substrate. If we choose other additional length other than 3LII, such as 2LII or 4LII, the value of the predicted bending stiffness shall be different because of the different bending length. However, the whole trends based on different choices are all similar in a robust way, and agree well with the whole trend of experimental data. Furthermore, to match the experimental data on values, we choose 3LII.

We attribute the difference between the theoretical prediction and the experimental measurement, especially when the bending angle of the sample is relatively small, to the following two factors.

First of all, in the experiment, the graphene may not be bent along the zigzag direction; that is, the interlayer shear action may differ from our fourth hypothesis to some extent. This also contributes to the difference between the theoretical prediction results and the experimental measurement results.

Second, the model to deduce the bending stiffness from experimental measurement is clever and simple, which, however, may ignore some secondary factors [27]. For example, when the bending angle is small, the part that is not obviously detached may be counted as a detached part in the model shown in Fig. 5(b). The bending stiffness is obtained from the balance between the bending energy and interface energy in experimental measurement [27]
Deff=RpHp2Rp(1cosθ)sin2θΓ
(42)
where Γ is the interfacial adhesion energy between graphene and hBN [42]. Thus, the overestimated detached length may lead to a systematic error of the experimental measurement results. Besides, the ruck may form near the upper step edge when the bending angle is relatively large, which is hard to be precisely described in the model [27].

Despite some differences, our theoretical predictions are in good agreement with the experimental measurement results. One of our key findings is presented: the bending stiffness is not a monotonic function of the bending angle due to the fact that the bending length also plays a vital role.

5 Further Results and Discussion

5.1 Trilayer Graphene

5.1.1 Theoretical Results.

In stable zones shown in Fig. 4, that is, κ^<κ^c or L^<L^c, we solve our model, revealing the mechanism clearly that the bending stiffness is determined by both the bending length and curvature, as shown in Fig. 7.

Fig. 7
Contour plots for the bending stiffness of trilayer graphene as a function of two governing parameters: L^ and κ^. (a) κ^<κ^c and L^<104. (b) L^<L^c and κ^<0.05.
Fig. 7
Contour plots for the bending stiffness of trilayer graphene as a function of two governing parameters: L^ and κ^. (a) κ^<κ^c and L^<104. (b) L^<L^c and κ^<0.05.
Close modal

We find that when the bending length is relatively large (κ^<κ^c and L^L^c), or the curvature is relatively small (L^<L^c and κ^κ^c), the bending stiffness is dominated by the bending length, increasing monotonically as bending length increases and only decreases insensitively with the increase of curvature. This find is consistent with the previous conclusion that the bending stiffness increases with the increase of length to some extent but extends the conclusion to a wider region, where the bending length can exceed 9.7 nm (L^40) [29].

Meanwhile, we can observe from Figs. 7(a) and 7(b) that when the curvature gradually decreases, regardless of the bending length, the bending stiffness always converges to an initial value, and this initial value is a monotonically increasing function of the bending length, which will be further discussed in Sec. 5.1.2.

The effect from curvature on the bending stiffness gradually becomes apparent as we further consider the condition of relatively large curvature (but L^<L^c) in the stable region. As shown in Fig. 7(b), when the curvature is close to the critical curvature (κ^κ^c), it is interesting that the increase of the bending length has a small effect on the increase of the bending stiffness, and even there exists an area where as the bending length increases, the bending stiffness decreases under the same bending curvature. In these regions with relatively large curvature, the curvature dominates the change of the bending stiffness.

Here, we take microscale and nanoscale bending behaviors as two typical examples to illustrate the implication of the above-mentioned results. On the one hand, when we consider the case of a large bending length, such as L^=104, i.e., 2.4 μm, which corresponds to a microscale bending behavior, for κ^=κ^c, the bending angle is up to 6528 deg. In other words, even if this trilayer graphene sample is wound for more than 18 turns around the cylindrical surface, the sample can still amazingly maintain its bending stiffness approximately the same as its initial value. On the other hand, when the bending length is a few nanometers, even if the bending angle is relatively small, the corresponding curvature can be considerable, and there may be a great difference between the bending stiffness and the initial value. For example, for a sample with bending length of 5 nm, when κ^ reaches κ^c, the bending angle is only about 14deg. It should be noted that the bending angle mentioned earlier is considered for the half sample when the sample has two free ends, or the whole sample when the sample has one clamped end and one free end. These discussions imply that the possible changes of bending stiffness ought to be considered when we design microscale or nanoscale out-of-plane deformation of vdW materials in order to precisely induce the anomalous physical phenomena or meet the functional requirements of the devices.

5.1.2 A Solution to the Flat State.

As discussed in Sec. 5.1.1, for a specific sample (a specific bending length), when the bending curvature tends to zero, the bending stiffness converges to a constant initial value. We call the state that the bending curvature is zero the flat state. For the flat state, κ → 0, the governing equation, Eq. (22), can degenerate into
d2u¯1dx¯122πτmaxE2Dl0u¯1=0
(43)
on account of 2πu¯1/[(1κd)l0]1 and 1 − κd ≈ 1. Use dimensionless quantities, u^=u¯1/l0, x^=x¯/l0, L^=L/l0, and κ^=κd, and introduce ω2=2πτmaxl0/E2D. Equation (43) then becomes
d2u^dx^2ω22u^=0
(44)
and the corresponding boundary condition becomes
{u^(0)=0du^dx^|x^=L^=κ^
(45)
Solving this differential equation and obtaining the displacement solution
u^=κ^sinh(ω2x^)ω2cosh(ω2L^)
(46)
Combining Eqs. (15) and (16), we obtain the strain solution
ε=κ^+κ^cosh(ω2x^)cosh(ω2L^)
(47)
The solution of the convex side can be deduced similarly. Combining Eqs. (17)(21), we obtain
Dκ^0eff=3B+E2Dd2(2+1cosh2(ω2L^)3sinh(ω2L^)ω2L^cosh(ω2L^))
(48)

Here, we take L^[1,104] as an example to illustrate the bending stiffness of trilayer graphene in the flat state, as shown in Fig. 8(a). On the one hand, the bending stiffness increases monotonically with the bending length in the flat state. One the other hand, the bending stiffness seems to have an upper bound and a lower bound. The upper bound can be analytically deduced from Eq. (48) as (3B + 2E2Dd2), and the lower bound is easily known as 3B. In other words, the bending stiffness of trilayer graphene predicted by our model can range from 4.32 eV to 495 eV, which we suppose that can partially explain the wide divergence on the measurements of bending stiffness [8,27,3034]. Revisiting the single-plate continuum model which is modified for the discrete nature [27,31], Dupeff=NB+E2Dd2(N3N)/12, we find the upper bound we predict for trilayer graphene is just the case of N = 3. Similarly, the lower bound we predict for trilayer graphene can also be seen as in the case of N = 3 for Dlowereff=NB under the ideal lubricated interlayer condition. What is more, on account of BE2Dd2/12 [38], the ratio between the upper bound and the lower bound can significantly exceed N3.

Fig. 8
(a) Plot of the effective bending stiffness versus bending length in the flat state for L^∈[1,104]. (b) Plot of the effective bending stiffness versus bending length in the flat state for L^<L^c, and schematic illustration of the modulation to bending stiffness from curvature.
Fig. 8
(a) Plot of the effective bending stiffness versus bending length in the flat state for L^∈[1,104]. (b) Plot of the effective bending stiffness versus bending length in the flat state for L^<L^c, and schematic illustration of the modulation to bending stiffness from curvature.
Close modal

For the sample with a specific bending length, as the curvature increases from the flat state, the bending stiffness generally decreases, and tends to the lower bound, as illustrated in Fig. 8(b). We may combine Figs. 7 and 8 to obtain a comprehensive understanding of how the bending length and curvature jointly determine the bending stiffness.

It should be pointed out that for a specific bending length, the bending stiffness is not a monotonic decreasing function of the bending curvature in the strict sense. Here, we use the sample length L^ equal to 20 as an example to show the transition of monotonicity. As shown in Fig. 9(a), the effective stiffness Deff decreases with the increase of the bending curvature κ^ in general, but when the bending curvature κ^ is about 0.04, there exists a small area where the monotonicity changes, which is magnified in the inset figure of Fig. 9(a).

Fig. 9
(a) Plot of the effective bending stiffness Deff versus curvature κ^ for L^=20. The inset figure of (a) magnifies the transition of monotonicity when κ^ is about 0.04. (b) Plot of the effective bending stiffness Deff versus bending angle θ for L^=13.5. The inset figure of (b) magnifies the transition of monotonicity when θ is about 35deg.
Fig. 9
(a) Plot of the effective bending stiffness Deff versus curvature κ^ for L^=20. The inset figure of (a) magnifies the transition of monotonicity when κ^ is about 0.04. (b) Plot of the effective bending stiffness Deff versus bending angle θ for L^=13.5. The inset figure of (b) magnifies the transition of monotonicity when θ is about 35deg.
Close modal

If we multiply the curvature by the bending length, we can also obtain the figure that shows the relationship between the bending stiffness and the bending angle. The format of the figure remains unchanged, but the abscissa is linearly mapped. We additionally take the sample length L^ equal to 13.5 as an example and consider the change of bending stiffness with respect to the bending angle, as shown in Fig. 9(b). Because of the reduction of bending length, the bending stiffness of the sample with bending length L^=13.5 is generally smaller than that of the sample with bending length L^=20. Furthermore, in DFT calculation, the relationship between the bending stiffness and the bending angle can also be obtained by compressing the sample to form a ruck [27], and the bending curve must experience the change of concavity, which implies that the corresponding bending length in our model is smaller than the real length of sample along the bending direction in DFT calculation. By comparison, it can be found that the DFT calculation result with length L^=20 is similar to the result in Fig. 9(b) [27].

5.2 Multilayer Graphene With More Than Three Layers.

According to the discussion in Sec. 2.2, we numerically solve our model for multilayer graphene with more than three layers under the condition of small curvature. Figure 10(a) shows the relationship between the bending stiffness and bending length for graphene with five layers, revealing that the bending stiffness is insensitive to curvature when the curvature is small, which is similar to the case of trilayer—considering that the larger layer number, the greater the overall slip under the same curvature, we use κ^=0.01/N to represent the small curvature bending state of different layers. We observe that the bending stiffness always has an upper bound as the bending length increases. Thus, for samples with more than three layers, we approximate the bending stiffness of L^=104 as the upper bound and normalize the bending stiffness, as shown in Fig. 10(b). It can be observed that for the cases with different layers under small curvature bending, the increasing trends of the bending stiffness with the bending length are generally similar to that of trilayer in the flat state, and merely the bending length which is corresponding to a rapid increase of the bending stiffness increases with the increase of the layer number.

Fig. 10
(a) Bending stiffness versus bending length for five layers under small curvature bending. (b) Normalized bending stiffness versus bending length for 5, 15, and 25 layers, respectively, under small curvature bending, and the comparison with the solution to trilayer in the flat state. In (a) and (b), L^∈[1,104].
Fig. 10
(a) Bending stiffness versus bending length for five layers under small curvature bending. (b) Normalized bending stiffness versus bending length for 5, 15, and 25 layers, respectively, under small curvature bending, and the comparison with the solution to trilayer in the flat state. In (a) and (b), L^∈[1,104].
Close modal

6 Conclusion

In this work, the mechanism of how bending stiffness of multilayer graphene with a specific layer number is jointly determined by the bending length and curvature is investigated theoretically. We establish a model to describe the multilayer graphene sheet bent along the cylindrical surface and obtain the energy expression form of bending stiffness. For trilayer graphene, our model theoretically predicts the experimental measurement results well and reveals that the bending stiffness is dominated by the bending length and curvature, respectively, in different stable regions, rather than solely by the bending angle.

As for the implications, we call for consideration of the possible change of bending stiffness, when we design the various microscale or nanoscale out-of-plane deformation patterns of vdW materials. On the one hand, when the bending length is relatively large (L^L^c and κ^<κ^c), or the curvature is relatively small (κ^κ^c and L^<L^c), the bending stiffness is dominated by the bending length. In this region, the bending stiffness increases monotonically as bending length increases, and only decreases insensitively with the increase of curvature. On the other hand, the effect of curvature on the bending stiffness becomes apparent when the curvature is relatively large (κ^κ^c and L^<L^c). These two results imply that while the microscale trilayer graphene is able to maintain its initial bending stiffness well even if it is rolled a dozen times, the bending stiffness of nanoscale trilayer graphene has to experience a significant change when the bending angle reaches only a dozen degrees.

The flat state, in which the bending stiffness reached its maximum for a specific bending length, is analytically explored for the trilayer graphene sheet, which indicates that the bending stiffness of trilayer graphene can range from 3B to (3B + 2E2Dd2), i.e., from 4.32 eV to 495 eV. We suppose that this large range of two magnitude orders is able to partially explain the wide divergence on the measurements of bending stiffness before. Furthermore, when multilayer graphene samples are bent along a cylindrical surface with small curvature, samples with more than three layers share similar characteristics with trilayer graphene. Although we have focused on multilayer graphene, our results should be extended to other vdW materials, giving fundamental guidance on the rational design of functional vdW structures with various out-of-plane deformation patterns.

Acknowledgment

This work was supported by the NSF of China through Grants Nos. 11890681, 11672301, 11521202, 11890682, and 11832010, and the Strategic Priority Research Program of Chinese Academy of Sciences through Grant No. XDB36000000. X. M. thanks Fei Pan, Pinshane Y. Huang, Guorui Wang, Yuan Hou, Hao Long, and Hongtian Qiu for their helpful discussions. We thank Zhaohe Dai for their useful comments and edits on the manuscript. We also thank Pinshane Y. Huang for data from the literature [27].

Appendix. Bending a Book Along a Cylindrical Surface

In order to explore the basic geometric relationship of the behavior that the sample is bent along a cylindrical surface, we bend a book, as shown in Fig. 11. In the experiment, we try our best to ensure that there are no obvious gaps between the sample layers; that is, the interlayer distance (the distance between the middle planes of adjacent layers) is equal to the thickness of each layer, so as to meet the third hypothesis of our model. If we choose some points in each layer and draw straight lines perpendicular to the deflection curve of the corresponding layer through these points, then these straight lines will converge to the same point, the center of curvature, which we denote as the point O′.

Fig. 11
Bending a book along a cylindrical surface. The middle dashed arc (green) corresponds to the middle plane of the book. The dashed curve in the northwest part (red) online) corresponds to the free end of the book. Other dashed curves (blue) correspond to the profile of the book except for the free end. All dashed curves are obtained from Eq. (A1). (Color version online.)
Fig. 11
Bending a book along a cylindrical surface. The middle dashed arc (green) corresponds to the middle plane of the book. The dashed curve in the northwest part (red) online) corresponds to the free end of the book. Other dashed curves (blue) correspond to the profile of the book except for the free end. All dashed curves are obtained from Eq. (A1). (Color version online.)
Close modal
In Fig. 11, we take the centroid of the top surface of the cylinder as the polar coordinate pole and let the pole diameter point to the clamped end of the book from the pole, establishing the polar coordinate system. L and t, respectively, denote the width and thickness of the book, and R denotes the curvature radius of the middle plane of book, whose value can be obtained from the relation that the radius of cylindrical top surface equals Rt/2. The coordinates of any point in the middle plane of the ith layer after bending is
{r=R+zizi[t2,t2]θ=lrl[0,Li]
(A1)
where l denotes the length of the arc along the path of the corresponding middle plane from the point to the clamped end, and Li denotes the length of the corresponding middle plane after bending. For the case that the interlayer is ideally lubricated, the middle plane of each layer remains neutral during the bending deformation, that is, Li = L.

As for Eq. (A1), by setting l = L, we can theoretically obtain the parametric equation of the cross section at the free end. Similarly, the parametric equations of other cross sections can be obtained, and some representative curves of profiles are drawn in Fig. 11 with dashed lines. These profiles, which are theoretically predicted, are in good agreement with those of experiments. While the cross section at the clamped end remains perpendicular to each layer, the cross section at the free end, whose profile corresponds to the dashed curve in the northwest part of Fig. 11, is no longer a plane and not perpendicular to each layer.

We derive Eq. (A1) solely from the geometry, and the properties of the material, including interlayer properties, do not enter into the discussion. Therefore, Eq. (A1) is valid irrespective of the properties of material, which we can use to consider the bending behavior with finite interlayer shear effect.

References

1.
Castro Neto
,
A. H.
,
Guinea
,
F.
,
Peres
,
N. M. R.
,
Novoselov
,
K. S.
, and
Geim
,
A. K.
,
2009
, “
The Electronic Properties of Graphene
,”
Rev. Mod. Phys.
,
81
(
1
), pp.
109
162
. 10.1103/RevModPhys.81.109
2.
Balandin
,
A. A.
,
Ghosh
,
S.
,
Bao
,
W.
,
Calizo
,
I.
,
Teweldebrhan
,
D.
,
Miao
,
F.
, and
Lau
,
C. N.
,
2008
, “
Superior Thermal Conductivity of Single-Layer Graphene
,”
Nano Lett.
,
8
(
3
), pp.
902
907
. 10.1021/nl0731872
3.
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
), pp.
385
388
. 10.1126/science.1157996
4.
Geim
,
A. K.
, and
Grigorieva
,
I. V.
,
2013
, “
Van der Waals Heterostructures
,”
Nature
,
499
(
7459
), pp.
419
425
. 10.1038/nature12385
5.
Novoselov
,
K. S.
,
Mishchenko
,
A.
,
Carvalho
,
A.
, and
Castro Neto
,
A. H.
,
2016
, “
2D Materials and Van der Waals Heterostructures
,”
Science
,
353
(
6298
), p.
aac9439
. 10.1126/science.aac9439
6.
Liu
,
Y.
,
Huang
,
Y.
, and
Duan
,
X.
,
2019
, “
Van der Waals Integration Before and Beyond Two-Dimensional Materials
,”
Nature
,
567
(
7748
), pp.
323
333
. 10.1038/s41586-019-1013-x
7.
Reis
,
P. M.
,
2015
, “
A Perspective on the Revival of Structural (in)Stability With Novel Opportunities for Function: From Buckliphobia to Buckliphilia
,”
ASME J. Appl. Mech.
,
82
(
11
), p.
111001
. 10.1115/1.4031456
8.
Huang
,
R.
,
2020
, “
Bending With Slip
,”
Nat. Mater.
,
19
(
3
), pp.
257
258
. 10.1038/s41563-020-0604-0
9.
Bao
,
W.
,
Miao
,
F.
,
Chen
,
Z.
,
Zhang
,
H.
,
Jang
,
W.
,
Dames
,
C.
, and
Lau
,
C. N.
,
2009
, “
Controlled Ripple Texturing of Suspended Graphene and Ultrathin Graphite Membranes
,”
Nat. Nanotechnol.
,
4
(
9
), pp.
562
566
. 10.1038/nnano.2009.191
10.
Castellanos-Gomez
,
A.
,
Roldan
,
R.
,
Cappelluti
,
E.
,
Buscema
,
M.
,
Guinea
,
F.
,
van der Zant
,
H. S. J.
, and
Steele
,
G. A.
,
2013
, “
Local Strain Engineering in Atomically Thin MoS2
,”
Nano Lett.
,
13
(
11
), pp.
5361
5366
. 10.1021/nl402875m
11.
Zheng
,
J.
,
Liu
,
H.
,
Wu
,
B.
,
Guo
,
Y.
,
Wu
,
T.
,
Yu
,
G.
,
Liu
,
Y.
, and
Zhu
,
D.
,
2011
, “
Production of High-Quality Carbon Nanoscrolls With Microwave Spark Assistance in Liquid Nitrogen
,”
Adv. Mater.
,
23
(
21
), pp.
2460
2463
. 10.1002/adma.201004759
12.
Lopez-Bezanilla
,
A.
,
Campos-Delgado
,
J.
,
Sumpter
,
B. G.
,
Baptista
,
D. L.
,
Hayashi
,
T.
,
Kim
,
Y. A.
,
Muramatsu
,
H.
,
Endo
,
M.
,
Achete
,
C. A.
,
Terrones
,
M.
, and
Meunier
,
V.
,
2012
, “
Geometric and Electronic Structure of Closed Graphene Edges
,”
J. Phys. Chem. Lett.
,
3
(
15
), pp.
2097
2102
. 10.1021/jz300695h
13.
Koenig
,
S. P.
,
Boddeti
,
N. G.
,
Dunn
,
M. L.
, and
Bunch
,
J. S.
,
2011
, “
Ultrastrong Adhesion of Graphene Membranes
,”
Nat. Nanotechnol.
,
6
(
9
), pp.
543
546
. 10.1038/nnano.2011.123
14.
Dai
,
Z.
,
Hou
,
Y.
,
Sanchez
,
D. A.
,
Wang
,
G.
,
Brennan
,
C. J.
,
Zhang
,
Z.
,
Liu
,
L.
, and
Lu
,
N.
,
2018
, “
Interface-Governed Deformation of Nanobubbles and Nanotents Formed by Two-Dimensional Materials
,”
Phys. Rev. Lett.
,
121
(
26
), p.
266101
. 10.1103/PhysRevLett.121.266101
15.
Rogers
,
J. A.
,
Someya
,
T.
, and
Huang
,
Y.
,
2010
, “
Materials and Mechanics for Stretchable Electronics
,”
Science
,
327
(
5973
), pp.
1603
1607
. 10.1126/science.1182383
16.
Yuk
,
J. M.
,
Park
,
J.
,
Ercius
,
P.
,
Kim
,
K.
,
Hellebusch
,
D. J.
,
Crommie
,
M. F.
,
Lee
,
J. Y.
,
Zettl
,
A.
, and
Alivisatos
,
A. P.
,
2012
, “
High-Resolution EM of Colloidal Nanocrystal Growth Using Graphene Liquid Cells
,”
Science
,
336
(
6077
), pp.
61
64
. 10.1126/science.1217654
17.
Chandratre
,
S.
, and
Sharma
,
P.
,
2012
, “
Coaxing Graphene to be Piezoelectric
,”
Appl. Phys. Lett.
,
100
(
2
), pp.
023114
. 10.1063/1.3676084
18.
Klimov
,
N. N.
,
Jung
,
S.
,
Zhu
,
S.
,
Li
,
T.
,
Wright
,
C. A.
,
Solares
,
S. D.
,
Newell
,
D. B.
,
Zhitenev
,
N. B.
, and
Stroscio
,
J. A.
,
2012
, “
Electromechanical Properties of Graphene Drumheads
,”
Science
,
336
(
6088
), pp.
1557
1561
. 10.1126/science.1220335
19.
Branny
,
A.
,
Kumar
,
S.
,
Proux
,
R.
, and
Gerardot
,
B. D.
,
2017
, “
Deterministic Strain-Induced Arrays of Quantum Emitters in a Two-Dimensional Semiconductor
,”
Nat. Commun.
,
8
(
1
), p.
15053
. 10.1038/ncomms15053
20.
Timoshenko
,
S.
, and
Woinowsky-Krieger
,
S.
,
1959
,
Theory of Plates and Shells
,
McGraw-Hill
,
New York
.
21.
Wang
,
G.
,
Dai
,
Z.
,
Xiao
,
J.
,
Feng
,
S.
,
Weng
,
C.
,
Liu
,
L.
,
Xu
,
Z.
,
Huang
,
R.
, and
Zhang
,
Z.
,
2019
, “
Bending of Multilayer Van der Waals Materials
,”
Phys. Rev. Lett.
,
123
(
11
), p.
116101
.
22.
Wang
,
G.
,
Dai
,
Z.
,
Wang
,
Y.
,
Tan
,
P.
,
Liu
,
L.
,
Xu
,
Z.
,
Wei
,
Y.
,
Huang
,
R.
, and
Zhang
,
Z.
,
2017
, “
Measuring Interlayer Shear Stress in Bilayer Graphene
,”
Phys. Rev. Lett.
,
119
(
3
), p.
036101
. 10.1103/PhysRevLett.119.036101
23.
Hod
,
O.
,
Meyer
,
E.
,
Zheng
,
Q.
, and
Urbakh
,
M.
,
2018
, “
Structural Superlubricity and Ultralow Friction Across the Length Scales
,”
Nature
,
563
(
7732
), pp.
485
492
. 10.1038/s41586-018-0704-z
24.
Annett
,
J.
, and
Cross
,
G. L.
,
2016
, “
Self-Assembly of Graphene Ribbons by Spontaneous Self-Tearing and Peeling From a Substrate
,”
Nature
,
535
(
7611
), pp.
271
275
. 10.1038/nature18304
25.
Wang
,
D.
,
Chen
,
G.
,
Li
,
C.
,
Cheng
,
M.
,
Yang
,
W.
,
Wu
,
S.
,
Xie
,
G.
,
Zhang
,
J.
,
Zhao
,
J.
,
Lu
,
X.
,
Chen
,
P.
,
Wang
,
G.
,
Meng
,
J.
,
Tang
,
J.
,
Yang
,
R.
,
He
,
C.
,
Liu
,
D.
,
Shi
,
D.
,
Watanabe
,
K.
,
Taniguchi
,
T.
,
Feng
,
J.
,
Zhang
,
Y.
, and
Zhang
,
G.
,
2016
, “
Thermally Induced Graphene Rotation on Hexagonal Boron Nitride
,”
Phys. Rev. Lett.
,
116
(
12
), p.
126101
. 10.1103/PhysRevLett.116.126101
26.
Shen
,
Y.
, and
Wu
,
H.
,
2012
, “
Interlayer Shear Effect on Multilayer Graphene Subjected to Bending
,”
Appl. Phys. Lett.
,
100
(
10
), pp.
101909
. 10.1063/1.3693390
27.
Han
,
E.
,
Yu
,
J.
,
Annevelink
,
E.
,
Son
,
J.
,
Kang
,
D. A.
,
Watanabe
,
K.
,
Taniguchi
,
T.
,
Ertekin
,
E.
,
Huang
,
P. Y.
, and
van der Zande
,
A. M.
,
2019
, “
Ultrasoft Slip-Mediated Bending in Few-Layer Graphene
,”
Nat. Mater.
,
19
(
3
), pp.
305
309
. 10.1038/s41563-019-0529-7
28.
Qu
,
W.
,
Bagchi
,
S.
,
Chen
,
X.
,
Chew
,
H. B.
, and
Ke
,
C.
,
2019
, “
Bending and Interlayer Shear Moduli of Ultrathin Boron Nitride Nanosheet
,”
J. Phys. D: Appl. Phys.
,
52
(
46
), p.
465301
.
29.
Pan
,
F.
,
Wang
,
G.
,
Liu
,
L.
,
Chen
,
Y.
,
Zhang
,
Z.
, and
Shi
,
X.
,
2019
, “
Bending Induced Interlayer Shearing, Rippling and Kink Buckling of Multilayered Graphene Sheets
,”
J. Mech. Phys. Solids
,
122
, pp.
340
363
. 10.1016/j.jmps.2018.09.019
30.
Chen
,
X.
,
Yi
,
C.
, and
Ke
,
C.
,
2015
, “
Bending Stiffness and Interlayer Shear Modulus of Few-Layer Graphene
,”
Appl. Phys. Lett.
,
106
(
10
), p.
101907
.
31.
Koskinen
,
P.
, and
Kit
,
O. O.
,
2010
, “
Approximate Modeling of Spherical Membranes
,”
Phys. Rev. B
,
82
(
23
), p.
235420
. 10.1103/PhysRevB.82.235420
32.
Lindahl
,
N.
,
Midtvedt
,
D.
,
Svensson
,
J.
,
Nerushev
,
O. A.
,
Lindvall
,
N.
,
Isacsson
,
A.
, and
Campbell
,
E. E.
,
2012
, “
Determination of the Bending Rigidity of Graphene via Electrostatic Actuation of Buckled Membranes
,”
Nano Lett.
,
12
(
7
), pp.
3526
3531
. 10.1021/nl301080v
33.
Zhang
,
D.-B.
,
Akatyeva
,
E.
, and
Dumitrica
,
T.
,
2011
, “
Bending Ultrathin Graphene at the Margins of Continuum Mechanics
,”
Phys. Rev. Lett.
,
106
(
25
), p.
255503
. 10.1103/PhysRevLett.106.255503
34.
Zhao
,
J.
,
Deng
,
Q.
,
Ly
,
T. H.
,
Han
,
G. H.
,
Sandeep
,
G.
, and
Rummeli
,
M. H.
,
2015
, “
Two-Dimensional Membrane as Elastic Shell With Proof on the Folds Revealed by Three-Dimensional Atomic Mapping
,”
Nat. Commun.
,
6
(
1
), p.
8935
. 10.1038/ncomms9935
35.
Hibbeler
,
R. C.
,
1994
,
Mechanics of Materials
,
Prentice Hall
,
Englewood Cliffs, NJ
.
36.
Wang
,
S.
,
Chen
,
Y.
,
Ma
,
Y.
,
Wang
,
Z.
, and
Zhang
,
J.
,
2017
, “
Size Effect on Interlayer Shear Between Graphene Sheets
,”
J. Appl. Phys.
,
122
(
7
), pp.
074301
. 10.1063/1.4997607
37.
Xu
,
G.
, and
Zhang
,
C.
,
2003
, “
Analysis of Dislocation Nucleation From a Crystal Surface Based on the Peierls–Nabarro Dislocation Model
,”
J. Mech. Phys. Solids
,
51
(
8
), pp.
1371
1394
. 10.1016/S0022-5096(03)00067-X
38.
Wei
,
Y.
,
Wang
,
B.
,
Wu
,
J.
,
Yang
,
R.
, and
Dunn
,
M. L.
,
2013
, “
Bending Rigidity and Gaussian Bending Stiffness of Single-Layered Graphene
,”
Nano Lett.
,
13
(
1
), pp.
26
30
. 10.1021/nl303168w
39.
Ochs
,
K.
,
2011
, “
A Comprehensive Analytical Solution of the Nonlinear Pendulum
,”
Eur. J. Phys.
,
32
(
2
), pp.
479
490
. 10.1088/0143-0807/32/2/019
40.
Peng
,
S.
, and
Wei
,
Y.
,
2016
, “
On the Influence of Interfacial Properties to the Bending Rigidity of Layered Structures
,”
J. Mech. Phys. Solids
,
92
, pp.
278
296
. 10.1016/j.jmps.2016.04.005
41.
Abramowitz
,
M.
, and
Stegun
,
I. A.
,
1965
,
Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables
,
Dover Publications
,
New York
.
42.
Sanchez
,
D. A.
,
Dai
,
Z.
,
Wang
,
P.
,
Cantu-Chavez
,
A.
,
Brennan
,
C. J.
,
Huang
,
R.
, and
Lu
,
N.
,
2018
, “
Mechanics of Spontaneously Formed Nanoblisters Trapped by Transferred 2D Crystals
,”
Proc. Natl. Acad. Sci. U. S. A.
,
115
(
31
), pp.
7884
7889
. 10.1073/pnas.1801551115