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

Authors: M. W. Roberts; C. B. Clemons; J. P. Wilber; G. W. Young; A. Buldum; D. D. Quinn
Journal: Journal of Nanotechnology
 (2010)

---

## Abstract

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 nNnm=1.62 eV.


---

## Body

## 1. Introduction

A graphene sheet is a one-atom thick lattice of carbon atoms arranged in hexagonal rings; see Figure1. 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 single-walled 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.Figure 1
Schematic of the structure of a graphene sheet.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–7]. Isolated graphene is especially exciting for its novel electronic transport properties [8–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–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 basic material constant is the flexural rigidity of the plate. For this model, we consider two sets of boundary conditions on the pair of unloaded edges of the plate. In one case, each of the two unloaded edges is prescribed to remain straight and at the same distance from the substrate as the loaded edges. The governing equation is then a partial differential equation. In the second case, we enforce periodic boundary conditions, for which the plate equation can be simplified to the classical Euler-Bernoulli beam equation with a nonlinear term arising from the van der Waals interaction. For both sets of boundary conditions, we numerically solve the governing equation of the continuum model. The solution depends on our choice of the flexural rigidity.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 Section3, 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.


## 2. 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(1)α(∂4u∂x4+2∂4u∂x2∂y2+∂4u∂y4)-T∂2u∂x2=q.
(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.Table 1
List of variables and constants.
SymbolDescriptionUnitsuDeflection of the free graphene sheetnmlLength of the free graphene sheet in thex directionnmwWidth of the free graphene sheet in they directionnmdInitial separation between the two graphene sheetsnmTLateral force on the free graphene sheetnN/nmf(u+d)Transverse force on the free graphene sheetnN/nm2αFlexural Rigidity of a graphene sheetnN nmc12, c6Constants for the Lennard-Jones (12-6) potentialnNnm13,   nN nm7σDensity of a graphene sheetatoms/nm2Figure 2
Schematic of the system of graphene sheets.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 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(2)W(s)=∫02π∫0∞σV(s2+r2)rdrdθ,
where V(ξ):=c12ξ-12-c6ξ-6 is the van der Waals energy between two nonbonded carbon atoms a distance ξ apart and σ=38.177nm-2. The constants in V are(3)c12=3.859×10-9nNnm13,c6=2.43×10-6nNnm7
(numerical values for c6, c12, and σ are taken from [25]). It follows that(4)f(s)=-σ∂W(s)∂s=2πσ2(c12s11-c6s5).Figure 3
Schematic of a single carbon atom which is distances     from an infinite graphene sheet.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(5)f(u+d)=2πσ2(c12(u+d)-11-c6(u+d)-5),
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.Figure 4
Plot of the transverse force, (4): σ = 38.177 atoms/nm2, c12=3.859×10-9 nN nm13, and c6=2.43×10-6 nN nm7.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–28]. Our choice for f¯, which describes the interaction of the plate with the foundation, is specific to graphene.


## 3. 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 edgesy=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.


### 3.1. 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 Figure5. 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 edges—the 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-4K). 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.Figure 5
Top view of the top graphene layer. The positions of atoms on the loaded edges of the layer are prescribed during the simulation. Atoms in the middle are free to relax after each displacement of the right edge.


### 3.2. 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(6)u(0,y)=u(l,y)=0,y∈[0,w],u(x,0)=u(x,w)=0,x∈[0,l].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 tox is zero at the boundaries x=0, l. Hence neither of the familiar boundary conditions ∂u/∂x=0 or ∂2u/∂x2=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(7)∂2u∂x2(0,y)=∂2u∂x2(l,y)=β,y∈[0,w],∂2u∂y2(x,0)=∂2u∂y2(x,w)=γ,x∈[0,l],
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 loadT*. 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 settingT 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,Eabs, at a particular grid point by(8)Eabs=|ucont-uatom|,
where ucont is the solution of (1) and uatom is the deflection predicted by the atomistic simulation. Likewise, the relative error, Erel, is defined by(9)Erel=Eabs|uatom+eps|,
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 distanced=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.Figure6 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).Interpolated atomistic simulation results (a) and continuum model solutions (b) with optimizedα for a 6.4 nm by 6.8 nm graphene sheet with d=0.31 nm and T=0 nN/nm.
(a)(b)Table2 summarizes our results for the boundary conditions considered in this section. The absolute error for all of these cases was less than 2×10-3nm. The relative error for the cases with a separation of d=0.37nm was less than 0.02, and the relative error for the cases with a separation of d=0.31nm 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.26nNnm=1.62eV.Table 2
Calculated flexural rigidities for 2D simulations.
Length (nm)Width (nm)Separation (nm)Flexural rigidity (nN nm)3.16.80.310.26013.16.80.370.26536.46.80.310.26126.46.80.370.2527


## 3.1. 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 Figure5. 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 edges—the 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-4K). 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.Figure 5
Top view of the top graphene layer. The positions of atoms on the loaded edges of the layer are prescribed during the simulation. Atoms in the middle are free to relax after each displacement of the right edge.


## 3.2. 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(6)u(0,y)=u(l,y)=0,y∈[0,w],u(x,0)=u(x,w)=0,x∈[0,l].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 tox is zero at the boundaries x=0, l. Hence neither of the familiar boundary conditions ∂u/∂x=0 or ∂2u/∂x2=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(7)∂2u∂x2(0,y)=∂2u∂x2(l,y)=β,y∈[0,w],∂2u∂y2(x,0)=∂2u∂y2(x,w)=γ,x∈[0,l],
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 loadT*. 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 settingT 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,Eabs, at a particular grid point by(8)Eabs=|ucont-uatom|,
where ucont is the solution of (1) and uatom is the deflection predicted by the atomistic simulation. Likewise, the relative error, Erel, is defined by(9)Erel=Eabs|uatom+eps|,
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 distanced=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.Figure6 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).Interpolated atomistic simulation results (a) and continuum model solutions (b) with optimizedα for a 6.4 nm by 6.8 nm graphene sheet with d=0.31 nm and T=0 nN/nm.
(a)(b)Table2 summarizes our results for the boundary conditions considered in this section. The absolute error for all of these cases was less than 2×10-3nm. The relative error for the cases with a separation of d=0.37nm was less than 0.02, and the relative error for the cases with a separation of d=0.31nm 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.26nNnm=1.62eV.Table 2
Calculated flexural rigidities for 2D simulations.
Length (nm)Width (nm)Separation (nm)Flexural rigidity (nN nm)3.16.80.310.26013.16.80.370.26536.46.80.310.26126.46.80.370.2527


## 4. 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 ofy. Hence we seek solutions to (1) that are independent of y, so that the plate equation reduces to(10)αd4udx4-Td2udx2=f(u+d).
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 distancesd 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.Figure 7
Atomistic simulation of the deflection of a 6 nm graphene sheet withd=0.31 nm. Note that the continuum model tries to match only the deflected part of the sheet, that is, from x=0 to about x=5.6 in this example.Here also we face the issue of correctly specifying conditions on the solutionu 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,(11)u(0)=u(l)=0.
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(12)d2udx2(0)=d2udx2(l)=β,
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.Figures8 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.Results for an optimization withα=.26 nN nm, l=1.7 nm, and T=0 nN/nm. In (a), d=0.31 nm. In (b), d=0.37 nm. Note that the atoms on the edges of the layer, whose positions are prescribed, are not plotted because the continuum model does not describe the deflection of these atoms.
(a)(b)Results for an optimization withα=.26 nN nm, l=6.8 nm, and T=0 nN/nm. In (a), d=0.31 nm. In (b), d=0.37 nm. Note that the atoms on the edges of the layer, whose positions are prescribed, are not plotted because the continuum model does not describe the deflection of these atoms.
(a)(b)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.


## 5. 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 about0.26nNnm=1.62eV.



---

*Source: 868492-2010-07-20.xml*