Continuum Plate Theory and Atomistic Modeling to Find the Flexural Rigidity of a Graphene Sheet Interacting with a Substrate

Using a combination of continuum modeling, atomistic simulations, and numerical optimization, we estimate the flexural rigidity of a graphene sheet. We consider a rectangular sheet that is initially parallel to a rigid substrate. The sheet interacts with the substrate by van der Waals forces and deflects in response to loading on a pair of opposite edges. To estimate the flexural rigidity, we model the graphene sheet as a continuum and numerically solve an appropriate differential equation for the transverse deflection. This solution depends on the flexural rigidity. We then use an optimization procedure to find the value of the flexural rigidity that minimizes the difference between the numerical solutions and the deflections predicted by atomistic simulations. This procedure predicts a flexural rigidity of 0.26 nN nm = 1.62 eV.


Introduction
A graphene sheet is a one-atom thick lattice of carbon atoms arranged in hexagonal rings; see Figure 1. Graphite, a naturally occurring form of carbon, has the structure of stacks of graphene sheets, and consequently the properties of graphene have long been of interest [1]. In recent years, interest in graphene has intensified, driven by the discovery of novel forms of carbon like buckyballs and singlewalled and multiwalled carbon nanotubes. Graphene is the fundamental structure for each of these new forms of carbon. For example, each wall of a nanotube is a graphene sheet wrapped into a cylinder.
Although graphene is a fundamental structure for other forms of carbon, for many years attempts to synthesize isolated graphene sheets failed [2,3]. However, recently both mechanical and chemical methods have been developed for isolating individual graphene sheets [4][5][6][7]. Isolated graphene is especially exciting for its novel electronic transport properties [8][9][10][11], and significant experimental work is currently underway to develop nanoelectromechanical (NEM) systems, such as resonators, switches, and valves, based on isolated graphene [12,13]. Electronic transport in graphene is coupled to its state of deformation. Hence to design NEM devices, it is essential to understand the basic mechanics of graphene. [14][15][16][17].
In this paper, we estimate the flexural rigidity of graphene by comparing atomistic and continuum models of a rectangular graphene sheet initially parallel to and supported by a rigid substrate.
The geometry of our problem seems fundamental for the study of the mechanics of graphene because mechanical exfoliation, a technique for isolating individual graphene from bulk graphite, yields graphene sheets supported by a rigid substrate [18]. Also, several proposed electronic devices are based on graphene in a geometry similar to that of our problem [19]. The sheet deflects in response to two loads: a transverse load from its interaction with the substrate by van der Waals forces, and a lateral load acting parallel to the substrate on a pair of opposite edges, which we shall refer to as the loaded edges. The van der Waals forces are derived from a 6-12 Lennard-Jones potential. We assume that the two loaded edges stay a fixed distance from the substrate. We consider sheets of various lengths.
A continuum description of this problem is to model the graphene sheet as a thin, isotropic, linearly elastic plate. The  We also model the deflection of a graphene sheet over a rigid substrate atomistically. The deflection of the sheet is determined by a combination of molecular dynamics and energy minimization techniques. We then use an optimization procedure to find the value of the flexural rigidity that minimizes the difference between the continuum and atomistic solutions. The procedure has three steps: an atomistic simulation to determine the deflection, at a chosen edge load, of a graphene sheet of a prescribed length whose loaded edges are a prescribed distance from the substrate; the construction of an objective function that computes, for a given flexural rigidity, the difference between the deflection determined in the first step and the numerical solution to the continuum equation at the same edge load; and optimization to find the flexural rigidity that minimizes the objective function.
This procedure predicts a flexural rigidity of about 0.26 nN nm = 1.62 eV. This estimate for the flexural rigidity is in reasonable agreement with the value of 0.238 nN nm derived in [20]. However, our reported value is significantly higher than the value of 0.136 nN nm derived in [21] and somewhat higher than the value of 0.192 nN nm reported in [22]. See also the discussion and estimates of flexural rigidity in [23].
In the next section, we recall the continuum model and corresponding governing equation for the transverse deflection of a graphene sheet. In Section 3, after first discussing the atomistic simulations, we present our main results using optimization to determine the flexural rigidity of graphene under boundary conditions for which the governing equation is a plate equation. Section 4 contains  analogous results under periodic boundary conditions on the unloaded edges. A final section briefly summarizes this work.

Governing Equations for the Transverse Deflection of a Graphene Sheet
To derive a continuum equation for a graphene sheet over a rigid substrate, we model the sheet as a thin, isotropic, linearly elastic rectangular plate. An appropriate theory for such plates is described in [24]. To apply this theory, we refer to Figure 2, which shows a schematic of our system. In this figure, l and w are the length and width of the sheet. We let u(x, y) denote the transverse deflection of the plate at the point (x, y). The plate deflects in response to a compressive load T applied to the pair of opposite edges at x = 0 and x = l. Also, there is a transverse load q on the plate, but no shearing force. The equilibrium equation can be written as (See [24]. We note that this model tracks only the transverse deflection u and ignores any in-plane shortening or extension of the sheet.) In (1), α is the flexural rigidity of the plate. For convenience, we include Table 1, which lists the variables and constants we use.
In (1), we must specify the form of the transverse load q, which in our problem describes the interaction of the graphene sheet with the rigid substrate. This interaction arises because, when the sheet and the substrate are sufficiently close, each atom on the substrate exerts a force on Density of a graphene sheet atoms/nm 2 each carbon atom on the sheet. We let f (s) denote the total van der Waals force per unit area exerted at a point on the sheet by all the atoms on the substrate when the point and the substrate are a distance s apart, as shown in Figure 3. To find f , we assume the atoms are continuously distributed on the substrate and the sheet with an atomic density σ, and we assume the substrate is an infinite plane. We get the total energy between a point on the sheet and the substrate by evaluating the improper integral where V (ξ) := c 12 ξ −12 − c 6 ξ −6 is the van der Waals energy between two nonbonded carbon atoms a distance ξ apart and σ = 38.177 nm −2 . The constants in V are c 12 = 3.859 × 10 −9 nN nm 13 , c 6 = 2.43 × 10 −6 nN nm 7 (numerical values for c 6 , c 12 , and σ are taken from [25]). It follows that From (4), if the loaded edges of a sheet are prescribed to be a distance d from the substrate and u is the deflection of the sheet at a point, then the transverse force per unit area at that point is which is used for q in (1). This transverse force may be attractive or repulsive. Figure 4 shows the graph of f . As the separation between the sheets approaches zero, the transverse force is repulsive and grows without bound. When the separation is at the equilibrium spacing of 0.34 nm, the transverse force is zero. As the separation between the sheets increases beyond this equilibrium spacing, the transverse force is attractive but quickly decays toward 0. Note that the problem we have formulated is a version of the classical problem of a plate or beam on a nonlinearly elastic foundation [24,[26][27][28]. Our choice for f , which describes the interaction of the plate with the foundation, is specific to graphene.

Optimization for Plate Deflections
In this section, we study the problem assuming that the sheet does not deflect at the unloaded edges, that is, the edges y = 0 and y = w (see Figure 2). We solve the continuum model, which is the plate equation (1), numerically. To determine the flexural rigidity α, we compare this solution to the results of atomistic simulations, the details of which we present next.

Solutions from Atomistic
Simulation. In our atomistic simulation, we model the interaction between the bonded atoms within the graphene sheet by the Tersoff-Brenner empirical potential [29]. Also, to describe the van der Waals interaction between the atoms in the sheet and the the atoms in the rigid substrate, we use Girifalco's Lennard-Jones potential V , which is given above, following (2).
We simulated the lateral loading of a graphene sheet by prescribing the positions of the atoms in the five rows closest to each of the loaded edges of the sheet; see Figure 5. At the start of a given simulation, these edge atoms were placed a prescribed distance d from the substrate. During a simulation, the edge atoms were moved inward toward the center of the sheet in increments of 0.01 nm while their distance from the substrate was kept unchanged. After each increment, the remaining atoms in the sheet were relaxed. During this relaxation, corresponding atoms along the unloaded edgesthe top and bottom edges in Figure 5-were identified to enforce the boundary conditions along these edges. The Nordsieck predictor-corrector integration algorithm was used to perform molecular dynamics (MD) annealing and  quenching simulations. Vibrational relaxations remained at even very low temperatures (1 × 10 −4 K). Therefore, the MD simulations were continued with steepest descent and conjugate gradient energy minimization algorithms. The steepest descent algorithm resulted in the best minimization. As the position of the prescribed atoms was incremented, the lateral load was determined numerically by computing the total energy of the system as a function of the edge displacement. A given simulation was continued until a desired lateral load was attained. The displacement at the edges y = 0 and y = w is fixed at zero. We simulate sheets 6.8 nm wide.

Boundary Conditions and Numerical Solutions for the Continuum Model.
To make a quantitative comparison of the deflections from the atomistic simulations and the numerical solutions of (1), we must specify the boundary conditions for (1) correctly, that is, in a way that reflects the boundary conditions in the simulations. Because the simulations are performed by prescribing the positions of strips of atoms near the edges of the sheet, we do not consider these atoms when determining the boundary conditions for (1). Rather, we define the left boundary of the sheet in the simulation, that is, the position in the sheet that corresponds to x = 0 in the continuum model, as the innermost row of atoms on the left whose position is prescribed. The right boundary in the simulation is defined analogously. It is clear that at the edges of the sheet we should enforce u 0, y = u l, y = 0, y ∈ [0, w], Determining the other boundary conditions is less straightforward. Computing derivatives numerically from the atomistic simulation data indicates that neither the first nor the second partial derivative of deflection with respect to x is zero at the boundaries x = 0, l. Hence neither of the familiar boundary conditions ∂u/∂x = 0 or ∂ 2 u/∂x 2 = 0 is appropriate. An analogous situation holds for the partials with respect to y at the boundaries y = 0, w. To handle this, we incorporate finding the second derivatives of the numerical solution at x = 0, l and y = 0, w into the optimization procedure. Specifically, at the two pairs of opposite edges of the sheet, we assume that the appropriate second derivatives of the deflection are constants and we determine these constants by the optimization. Hence we assume where the values of β and γ are to be determined as part of the optimization procedure that finds the value of the flexural rigidity α.
As part of the optimization procedure, (1) subject to the boundary conditions (6), (7) is solved numerically. We solve (1) using second-order central-difference discretizations. We use the method of successive approximations to handle the nonlinear transverse force, f . To implement this method, the system is solved iteratively. The solution to a given iteration, u, serves as an approximate solution, u, for the next iteration. For the successive approximations, f is approximated by a two-term Taylor expansion around u. At each step in the iteration, the result of the Taylor expansion is a system of linear equations, which is solved exactly to find u. Iterations are continued until the difference between u and u becomes sufficiently small. For the numerical solutions in this paper, the summation of the difference over the entire plate had to be less than 1 × 10 −5 to stop the iteration; see [30] for details.
We describe the optimization procedure next. The first step in this procedure is to perform an atomistic simulation for a graphene sheet of a chosen length and distance from the substrate. The data from this simulation are the positions of all the atoms on the sheet over some range of lateral loads. From the data we take the atomic positions corresponding to a chosen lateral load T * . The simulation data are on a relatively coarse grid, and hence to compare the data to the numerical solutions we first interpolate the position of the sheet between the atomic positions. We use cubic splines for this interpolation. We note that only the atoms in the relaxed region and on the two lines that form the boundary between the prescibed region and the relaxed region in Figure 5 were included in the cubic spline interpolation of the simulation results. After setting T in (1) equal to T * above, we define an objective function that, for its inputs, takes the flexural rigidity α and the two constants β and γ from the boundary condition (7). The objective function solves (1) numerically and computes the total absolute error between this numerical solution and the interpolated data. The MATLAB function fminsearch is used to find the combination of α, β, and γ that minimizes the objective function, that is, that minimizes the error between the numerical solution of (1) and the interpolated deflection data from the simulation. Note that in addition to a combination of α, β, and γ that minimizes the objective function, the procedure yields the numerical solution to (1) for this α, β, and γ.
To perform the optimization, we define the absolute error, E abs , at a particular grid point by where u cont is the solution of (1) and u atom is the deflection predicted by the atomistic simulation. Likewise, the relative error, E rel , is defined by where eps is machine epsilon (2.2204 × 10 −16 in MATLAB). We note that the maximum of the absolute error between the simulation data and the numerical solutions over the entire sheet was generally less than 1 × 10 −3 nm. The maximum relative error over the entire sheet was generally less than 0.1. We performed atomistic simulations for graphene sheets of two lengths, 3.2 nm and 6.4 nm. In each case the sheets had width 6.8 nm. For each length, we ran simulations keeping the prescribed atoms a distance d = 0.31 nm from the substrate as well as simulations with d = 0.37 nm. From the various simulations, we took configurations of the sheet corresponding to zero lateral load along the loaded edge. These configurations served as the basis for the optimization procedure. Figure 6 compares a configuration from one of the simulations with a corresponding numerical solution. Figure 6(a) shows a configuration constructed by interpolating MD simulation data; Figure 6(b) shows the numerical solution of (1) in which the value for α was determined using the data from Figure 6(a). Table 2 summarizes our results for the boundary conditions considered in this section. The absolute error for all of these cases was less than 2 × 10 −3 nm. The relative error for the cases with a separation of d = 0.37 nm was less than 0.02, and the relative error for the cases with a separation of d = 0.31 nm was less than 0.04. We note the slight variation in the predicted flexural rigidities in Table 2.
We believe these deviations occur because near the minima, the graph of the objective function is very "flat", that is, the function is nearly constant. Hence even with the tolerance set at machine epsilon, the algorithm cannot locate the minima precisely. Using the results in Table 2, we conclude that α is approximately 0.26 nN nm = 1.62 eV.

Optimization for Beam-Like Deflections
In this section we present results analogous to those of the previous section but under periodic boundary conditions on the unloaded edges. Specifically, we assume in both the atomistic simulations and the continuum problem that the two unloaded edges undergo the same deflection. In this case, the atomistic simulations yield deflections of the sheet that are independent of y. Hence we seek solutions to (1) that are independent of y, so that the plate equation reduces to Although (10) is formally the same as the equation for an Euler-Bernoulli beam on an elastic foundation, we continue to interpret α as the flexural rigidity of plate rather than the bending stiffness of a beam. We note that, because this problem is a version of the classical problem of a beam on an elastic foundation, one expects the continuum model to predict the buckling of the sheet at some critical load. Buckling is not observed here because, in the simulations we report, the edge load does not exceed the critical load; see [31], in which a continuum model is used to study the buckling of a graphene sheet supported by a rigid substrate in a geometry similar to that studied here.
Equation (10) is an ordinary differential equation, and its solution u describes the deflection of a cross-section of the sheet. Based upon Table 2, we assume that the flexural rigidity α = 0.26 nN nm, and we compare the corresponding solution to the results of atomistic simulations.
We ran atomistic simulations for graphene sheets of various lengths and of various distances d from the substrate. Figure 7 shows the data obtained from a typical simulation. The boundary conditions yield deflections of the sheet that are independent of y. Hence all the simulations discussed in this section were run on a narrow 1 nm wide strip of graphene.
Here also we face the issue of correctly specifying conditions on the solution u of (10) at x = 0 and l that match the conditions at the boundary in the atomistic simulations.
We require that the transverse deflection of the graphene sheet is zero on the boundary, that is, Issues similar to those discussed in Section 3 arise when determining the other boundary condition at x = 0, l. Just as above, we incorporate finding the second derivative of the numerical solution at x = 0 and l into the optimization procedure by taking the remaining boundary conditions for (10) as where β is the unknown constant. The value of β is then found by the optimization procedure, which is similar to the procedure outlined in Section 3. In this case, to solve the continuum equation, we use MATLAB's built-in 1D boundary value-problem solver (bvpsolve). The solver uses a Labatto IIIc quadrature algorithm to find numerical solutions on an arbitrary grid [32]. We now solve (10) subject to (11), (12) for sheets of various lengths and at various distances d from the substrate. In particular, we consider sheets of several lengths between 1.7 nm and 21 nm with d = 0.31 nm. We also consider sheets of the same lengths but with d = 0.37 nm. For all the cases just mentioned, we used simulation data corresponding to a lateral load, T, of zero. Figures 8 and 9 depict numerical solutions determined as part of the optimization method. In Figures 8(a) and 9(a), the distance d = 0.3 nm is less than the equilibrium distance of approximately 0.34 nm, and hence the transverse force is repulsive and the deflection is positive. In Figures 8(b) and 9(b), the distance d = 0.37 nm is greater than the equilibrium distance, and hence the transverse force is attractive and the deflection is negative. In both cases, the deflection in the middle of the sheet is close to the equilibrium spacing of 0.34 nm. Note that the scale on the x-axis is much larger than the scale on the u-axis, and hence the deflection of the sheet in both cases is in fact small. These figures illustrate that for each case a value of β can be found so that the atomistic and continuum modeling are in good agreement. Hence these results confirm the value of α = .26 nN nm for the flexural rigidity found in the previous section.
According to the continuum model, α is a material property and hence is independent of the lateral load. To test this, we also considered data for deflections corresponding to nonzero lateral loads. We considered both tensile and compressive lateral loads for sheets 3.4 nm and 6.8 nm in length at separations of d = 0.31 nm and d = 0.37 nm. We discovered that for any given length and separation, the value α = 0.26 nN nm resulted in deflections that matched closely with those calculated by the atomistic simulations.

Summary and Conclusions
Using an optimization procedure based upon a comparison between atomistic simulations and continuum plate theory, we estimate the flexural rigidity for a graphene sheet over a rigid substrate. This procedure yields flexural rigidities of about 0.26 nN nm = 1.62 eV.