- Email: [email protected]
Sandwich plate model of multilayer graphene sheets for considering interlayer shear effect in vibration analysis via molecular dynamics simulations
Accepted Manuscript
Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations Reza Nazemnezhad , Mojtaba Zare , Shahrokh Hosseini-Hashemi PII: DOI: Reference:
S0307-904X(17)30182-8 10.1016/j.apm.2017.03.033 APM 11674
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
31 July 2016 25 January 2017 14 March 2017
Please cite this article as: Reza Nazemnezhad , Mojtaba Zare , Shahrokh Hosseini-Hashemi , Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.03.033
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
AC
CE
PT
ED
M
AN US
The van der Waals interactions of every two adjacent layers are considered in analysis which results in interlayer shear effect Sandwich model (SM) is used to obtain the governing equations. Molecular Dynamics (MD) simulations are implemented to verify the model. Consistent values for SM parameters are obtained to comply with MD results. SM can predict the vibration behavior of MLGSs well for different values of aspect ratio.
CR IP T
Page 1 of 28
ACCEPTED MANUSCRIPT
Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations
CR IP T
Reza Nazemnezhada,b, Mojtaba Zarea,c*, Shahrokh Hosseini-Hashemia,d
AN US
a- School of Mechanical Engineering, Iran University of Science and Technology, 16842-13114, Narmak, Tehran, Iran b- School of Engineering, Damghan University, 36715-364, Damghan, Semnan, Iran c- Department of Piping, Petrochemical Industries Design and Engineering Company (PIDEC), Eram BLVD, 714 55-618, Shiraz, Iran d- Center of Excellence in Railway Transportation, Iran University of Science and Technology, 16842-13114, Narmak, Tehran, Iran
Abstract
Van der Waals (vdWs) bindings acting as interlayer shear force between graphene layers of multi-layer graphene sheets (MLGSs) are considered in vibration analysis in the present study.
M
To idealize the structure of MLGS incorporating interlayer shear interactions, a sandwich model (SM) is represented which laminates the graphene layers. The layers stick together with vdWs
ED
bonds. The bonds are modeled as core layers between every two adjacent layers. Molecular dynamic (MD) simulation is carried out to validate the results obtained by SM. Afterward, the
PT
values for bending rigidity and layers thicknesses are obtained so as to match SM frequencies with MD results. It is observed the SM can predict the vibration behavior of MLGSs well for
CE
different values of aspect ratio. The present paper deals with a new method of incorporating shear effect and a novel investigation of integrating SM with MD.
AC
Keywords: Multi-layer graphene sheets, Sandwich model, vibration analysis, Interlayer shear effect, Molecular Dynamics. 1. Introduction Graphene, the world’s strongest material, a million times thinner than paper but two hundred times stronger than steel is rapidly rising star on the horizon of materials science and condensed
*
Corresponding author, E-mail: [email protected]
Page 2 of 28
ACCEPTED MANUSCRIPT
matter physics. Flat monolayer of carbon atoms tightly packed into a two-dimensional honeycomb lattice form graphene. This strictly 2D material exhibits exceptionally high crystal and electronic quality and potential applications. Widely used for describing properties of various carbon-based materials, graphene has been studied for sixty years [1-3]. Many possible applications of graphene-based materials, especially as electronic devices, have been reported [4-
CR IP T
6]. Graphenes may exist as separate single-layer or multiple layers held together by the weak vdWs forces. The latter is called multi-layer graphene sheet (MLGS). Stacking the layers of graphene changes the mechanical and electrical properties of MLGSs with respect to monolayer graphene. Hence, analyzing the mechanical behavior of MLGSs is significant. The vdWs binding produces
AN US
two types of forces: shear forces which resist interlayer displacements and tension-compression forces which resist out-of-plane stretches. MLGSs have generally two families of mode shapes: lower classical synchronized modes, and higher modes [7, 8]. Lower classical synchronized modes are those modes at the lower end of the frequency spectrum which exhibit no or very little out-of-plane relative motion between the layers during the vibration and hence are not affected
M
by the interlayer Young's modulus effect of vdWs interactions. The vdWs forces between the layers act as a series of springs connected between the layers’ nodes with extremely large spring
ED
constants. Thereby, all those modes with relative motion between the layers occur at much higher natural frequencies. These findings show, at low mode numbers the interlayer Young's modulus does not have any role in vibrational behavior of MLGSs. On the contrary, researchers
PT
considering the interlayer shear effects on the free vibration of MLGSs have concluded, at low mode numbers graphene layers have interlayer shear deformation or relative sliding and
CE
observed the apparent impacts of this shear deformation on the free vibration of MLGSs. Despite the well investigation of mechanical properties of monolayer graphene, the experimental
AC
results of multi-layer graphene lack due to the limits of experimental technique. Using molecular dynamic (MD) simulations and introducing a new model called multibeam shear model (MBSM), Liu et al. [9] considered the interlayer shear modulus effect on the vibrational behavior of graphene multilayer resonators. They also showed the classical Euler–Bernoulli and Timoshenko theories cannot predict the vibrational behavior of multilayer graphene nanoribbons (MLGNRs). Shen and Wu [10] used the MBSM and MD simulations to investigate the interlayer shear modulus effect on bending property of MLGNRs. Their results exhibited certain difference Page 3 of 28
ACCEPTED MANUSCRIPT
between MD simulations and MBSM in bending response predictions, which becomes more predominant when the layers’ number is smaller than 5. For better investigation of the interlayer shear modulus effect on bending behavior of MLGNRs, Liu et al. [11] employed Newmark’s composite beam theory and extended to multilayer case. Compared to the MD simulations, the Newmark’s beam model predicts much better results than the previous beam model, MBSM .
CR IP T
Nazemnezhad et al. [12] proposed a sandwich beam model to consider the interlayer shear effect of a bilayer graphene nanoribbon for vibration analysis. The model took the in-plane stretch into account. They studied the natural frequencies of a cantilever nanoribbon for a wide range of bending rigidity and graphene layer thickness and used MD simulations to match the best results for the model. Following to this study they investigated the vibration behavior of MLGNRs
AN US
considering vdWs bindings in nonlocal elasticity [13]. They inserted core layers acting as vdWs bonds in between every two graphene layers and modeled the structure according to sandwich beam. They obtained the best values for bending rigidity, core and graphene layers thicknesses as well as nonlocal parameter to match their results with MD. Nazemnezhad [14] and Nazemnezhad et al. [15] also integrated Timoshenko and Reddy beam theories with nonlocal elasticity to study
M
the vibration behavior of MLGNRs. They considered interlayer interactions of vdWs bindings and calibrated the nonlocal parameter for the case of their studies with the aids of MD
thickness in their studies.
ED
simulations. They proposed the most consistent values for bending rigidity and graphene
The literature survey shows the shear force caused by vdWs bindings has been well incorporated
PT
into mechanical behavior of MLGNRs. However, based on the knowledge of authors this impact has not been studied for the case of graphene sheets with multiple layers. The authors believe
CE
this investigation may contribute remarkable results for vibration analysis of MLGSs. The preliminary goal of this paper is to investigate the interlayer shear effect on natural frequencies
AC
of MLGSs. In addition, the main question is: is it possible to extend sandwich beam model to sandwich plate model which predicts the vibration behavior of MLGSs considering vdWs bindings? These are the topics which are discussed here. The MLGS structure is idealized based on the laminated graphene sheets bond together by core layers. The cores represent vdWs bindings. The displacement fields are considered according to classical (Kirchhoff) plate theory. The formulation is derived based on bilayer graphene sheet and the result is extended over MLGSs with n layers. The boundary conditions studied here resemble a cantilever plate, i.e. a Page 4 of 28
ACCEPTED MANUSCRIPT
sheet with three free edges and one clamped edge. Harmonic differential quadrature method is used to solve the governing equations of motion. Since the experimental investigations have not dealt with the case of present study yet, we simulate our model according to MD with the aids of software LAMMPS. The MD results serve as a benchmark to obtain suitable values for effective parameters of our model. The effective parameters are bending rigidity of graphene and
CR IP T
thicknesses of core and graphene layers. You will observe the sandwich model can predict the first two natural frequencies of MLGSs with different layers’ number and aspect ratio. The average error between the present model and MD shows a better perception of low deviation of sandwich model from MD. The present paper embodies a pioneer study in vibration analysis of
AN US
MLGSs incorporating interlayer shear forces caused by vdWs bindings.
2. Problem Formulation
In this section we are about to model a multi-layer graphene sheet depicted in Fig. 1 based on a sandwich plate model depicted in Fig. 2. The sheet consists of n-layers of graphene bonded together by n-1-layers of core. What the core does is indeed representing the vdWs bindings of
M
every two adjacent graphene sheets. The core is assumed to have negligible inertia and
AC
CE
PT
ED
undergoes shear stress alone while its thickness remains unchanged during deformation.
Fig. 1: Atomic structure schematic of a multilayer graphene sheet with arrangement of coordinates
The total thickness of MLGS is nH where H is the thickness of each graphene layer. Now let’s define the parameters appeared in Fig. 2. L is the length of the sheet in x-direction, d is the width in y-direction, hg and hc are the graphene sheet and core thicknesses, respectively. But before we
Page 5 of 28
ACCEPTED MANUSCRIPT
write the equations of motion for MLGS we simply explain the derivation of the equation for a
AN US
CR IP T
bilayer sheet and extend the result to a multi-layer one.
Fig. 2: Sandwich model idealization of a multilayer graphene sheet
A two-dimensional x-z plane view of a bilayer graphene sheet is illustrated in Fig. 3. The sheet
M
consists of a top and a bottom graphene layers with thicknesses hg and a core layer with thickness hc. The two graphene layers are assumed to be homogeneous and modeled according to
ED
Kirchhoff plate theory. The vdWs bindings are so perfect that no delamination occurs at the interfaces of core and graphene layers. Consequently, there is no relative flexural displacement
PT
between graphene layers. z
Top layer
AC
CE
hg hc hg
Core x L
Bottom layer
Fig. 3: Two dimensional representation of a bilayer graphene sheet
The displacement fields for the top and bottom layers of the sheet, according to the classical plate theory are as follows: (1)
Page 6 of 28
ACCEPTED MANUSCRIPT
*
+
(2)
*
+
(3)
where the subscript t and b stand for top and bottom layers, respectively and the asterisk sign
and bottom layers, respectively.
and plus sign ,
and
for the displacement equations of top are the displacement vector components of
CR IP T
may be interchanged to minus sign
top and bottom neutral plane. Furthermore, the variations of coordinate along the thickness of sheet are
and
for the top and bottom layers,
respectively. With the assumption of no delamination, one can deduce the lateral displacements
AN US
of all layers are equal and independent of z. Therefore, we may write:
(4)
in which c denotes the core. Considering linear distribution of the top and bottom layers’ inplane displacement, the core longitudinal and transverse displacement fields may be obtained as
ED
M
follows:
(
)
(5)
(
)
(6)
PT
The governing equations of motion can be derived by the Hamilton’s principle. The terms in the principle are referred to APPENDIX A. After some mathematical manipulations and
AC
:
CE
simplifications, we arrive to the following governing equations of motion:
Page 7 of 28
ACCEPTED MANUSCRIPT
(
)
(
(
(
) )
)
(7) )
(
(
)
(
(
(
(
) )
)
)
(
(
)
M
(
)
PT ⁄
CE
(9)
)
(
(10) )
(
In which
(8)
)
ED
(
)
)
AN US
(
(
)
CR IP T
(
(
)(
))
(11)
is the bending rigidity,
∫
⁄ ⁄
and
⁄
is the mass density of graphene. Furthermore, the boundary conditions studied here are specified
AC
as follows:
Free edge at
( (
:
) )
( (
) (12)
)
Page 8 of 28
ACCEPTED MANUSCRIPT
) (
Clamped edge at
)
(
(
)
)
:
)
(
(
:
Free edge at (
)
CR IP T
(
)
(
)
(
)
AN US
(
) (
)
(
)
(
(
)
)
(13)
(14)
Furthermore, if a plate has free corners, an additional boundary condition should be satisfied as
Free corners at
ED
M
[16]:
(15)
PT
The same procedure may be carried out to achieve the governing equations of motion for an n-
CE
layers graphene sheet. The final results are expressed as follows:
AC
(
) (
( (
(
(
( (
) )
) (
(16)
)
)
)
)
(17)
( (
) )
Page 9 of 28
)
(18)
ACCEPTED MANUSCRIPT
)
(
(
(
(
)
)
(
)
)
)
)
AN US
(
is the two
and free for all other edges are written as follows:
(
)
ED
)
(
) )
(
(
)
)
CE
(23)
PT
(
At
(22)
. Again the boundary conditions
At (
(21)
M
of MLGs for clamped edge at
(24)
AC
At (
)
(
dimensional Laplace operator and
(20)
is the shear modulus for the ith core layer and in the last equation
where
)
(
( )
(19)
)
)
(
)
(
( (
)
CR IP T
(
)
(
) (
(
) )
(
At corners
Page 10 of 28
(
)
)
(25)
ACCEPTED MANUSCRIPT
(26)
3. Solution procedure We now proceed with solving Eqs. (16)-(22) via the harmonic differential quadrature method (HDQM) [17]. The weighting functions of the HDQ are the trigonometric Sine and Cosine
CR IP T
functions. In this method, a linear summation of the values of weighted functions at all of the grid points in the entire domain of a spatial variable is used to approximate the derivative of a function with respect to the spatial variable at a given discrete point. For a two dimensional analysis we consider a function [
defined in a rectangular domain
[
]
]. The grid is obtained by taking N and M points along the x and y coordinates, respectively. are expressed in APPENDIX A.
AN US
The derivatives of the function
For generalizing the problem, the non-dimensional parameters are defined as below:
M
(27)
where
√
the cores, i.e.
{
and
}
(28)
is the natural frequency. We also assume the same shear modulus for all
PT
{ ̅ ̅ ̅}
ED
Also we consider the harmonic motion with respect to the time for displacements as:
. Therefore, applying Eqs. (A.5)-(A.10) into Eqs. (16)-(22) we obtain the
AC
CE
non-dimensional discretized form of the equations of motion as follows:
Page 11 of 28
ACCEPTED MANUSCRIPT
(
∑ ̅
)
̅
∑∑
(
(
∑
(
)
(
)
∑ ̅
)
∑∑ (
AC
CE
PT
(
[
(
)
(29)
(
))
(30)
∑ ̅
)
̅
∑∑
( (
)
̅
ED
)
)
(
M
(
(
∑
∑
))
AN US
(
CR IP T
∑
)
(
)]
∑
(
)
) (31)
Page 12 of 28
ACCEPTED MANUSCRIPT
∑ ̅
∑
(
(
))
(
)
)
∑ ̅
(
AN US
∑
∑ ̅
CR IP T
)
)
̅
∑∑
( (
(
)
(
)
(
))
M
(
)
̅
∑∑
(
(32)
AC
CE
PT
∑ ̅
)
ED
(
( (
(33)
∑
)
)
̅
∑∑
(
(
(
))
∑ ̅
) (34)
Page 13 of 28
ACCEPTED MANUSCRIPT
∑
(
∑∑
)
̅
(
)
∑
( ∑
(
(
∑ ̅
)
)
)
∑ ̅
(∑
(
)
∑ ̅
AN US
(
(
CR IP T
∑ ̅
)
)
)
(35)
[ [ [
] [ ] [
] { ]{ ] {
ED
M
The assembly of the Eqs. (29) through (35) and the quadrature analog equations of the boundary conditions (Eqs. (23)-(25) and (15)) mentioned in APPENDIX A yields the following matrix equation: } } }
[
[ ]
[
] [
[ ]
{ ]{ ] {
} } }
PT
Where the subscripts B and D stand for boundary and domain and {
(36)
} and {
} represent the
follows: {
}
{
} {
}
{
}
}
AC
{
CE
boundary and domain degrees of freedom of the displacement fields, respectively defined as
{
}
(37) {
Page 14 of 28
{
}
{
}
{
}
Note the term
CR IP T
ACCEPTED MANUSCRIPT
(
is the abbreviation for
(36) as a system of two linear equations with the variables {
). Considering Eq.
} and {
}, after doing some
mathematical manipulations and collecting terms, the characteristic equation for the natural
[
]
[
][
] [
]
[
]
AN US
frequencies of the multi-layer graphene sheet is expressed as: [
][
] [
] {
}
Eq. (38) is an eigenvalue problem with the characteristic value eigenvector {
{
}
(38)
associated with the
}. The eigenvalues and subsequently the natural frequencies of the multi-layer
M
graphene sheet can be obtained with the aids of a self-developed Mathematica program. 4. Simulation procedure
A wide set of properties and applications of graphene is related to interaction between graphene
ED
layers. Thermal conductivity, super lubricity, fast diffusion and self-retracting motion of graphene are among disparate properties and phenomena induced by interlayer van der Waals
PT
interaction. For felicitous consideration of the phenomena, nanoelectromechanical systems and realistic simulation of graphene-based systems, exhaustive information on the interaction
CE
between graphene layers is necessary The existing continuum models for multilayer graphene sheets decline to reflect the shear energy
AC
between the graphene layers since the discrete graphene lattice registry is not incorporated in the continuum models. Interlayer interaction is significant to the bending properties of multi-layer material. Since the very small thickness of the graphene results in a very small bending rigidity, the interlayer shear energy may play a very pivotal role in bending and vibration of multilayer graphene sheets and must be considered in the related continuum models. Interlayer shear modulus is three orders lower than in-plane Young’s modulus, thus the bending behavior is dominated by the interlayer shear. In multilayer graphene sheets considered in this model, the Page 15 of 28
ACCEPTED MANUSCRIPT
interlayer van der Waals binding is strong enough that negates the flexibility for interlayer shear deformation or relative sliding. Since, neither classical continuum models nor experimental techniques are able to account for the interlayer interactions, we perform the MD simulation to obtain a valid theoretical benchmark for our analysis. This section of the study is devoted to the simulation of MLGSs via molecular dynamics. The
CR IP T
major reason of this activity is to validate the natural frequencies results obtained by sandwich plate theory, since no valid experimental results have been published for the present problem. In this regard, LAMMPS software is used to perform the simulation. But before we go on, we ought to select an appropriate potential for mechanical properties of graphene which can consider the bonding atomic interaction in MLGSs well. The potential chosen here is the familiar AIREBO
AN US
(Adaptive Intermolecular Reactive Empirical Bond Order) which has broadly been used to study the mechanical properties of carbon nanotubes and graphene [18-20] and is evaluated to describe both the intralayer and interlayer interactions between carbon atoms. The parameter of E_LJ term in the potential is epsilon_CC in which its value indicates the strength of the vdWs interactions in the potential. Nazemnezhad et al. [12] developed their simulation based on two
M
values of the parameter for vibration analysis of multi-layer graphene nanoribbons. In this section we want to perform the simulations based on these two values to see whether the
ED
sandwich model can generally predict the vibration behavior of MLGSs. In this regard, this value sets to 2.84 meV first and the simulation is carried out. Afterwards this parameter is changed from 2.84 meV to 45.44 meV as done by Shen and Wu [21] for investigating the effect of
PT
interlayer shear on bending property of multi-layer graphene nanoribbons. In the cantilever MLGSs studied here the armchair and zigzag directions are along the length and width,
CE
respectively. The cross-section has different heights concerning to the number of layers n (n=2,…, 5). The graphene monolayer thickness and width are 0.335 nm and 5 nm respectively,
AC
while the aspect ratio is varied. After structural relaxation we use the following method for the MD simulation. At the beginning of the simulation a displacement of 1 nm is applied to the free end of each cantilever GS to initialize a transverse vibration. The deformed shape of a four-layer GS is represented in Fig. 4 as a sample.
Page 16 of 28
ACCEPTED MANUSCRIPT
Fig. 4. Initial deflection of a four-layer cantilever graphene sheet with an amplitude of 1 nm at its free end (
The system then evolves in a NVE ensemble. The equations of motion are integrated numerically using the Velocity-Verlet algorithm with a time step of 1 fs. The free end transversal motion of
CR IP T
the GSs is tracked. Fig 5 brings to you the graph of end transversal motion versus time for a four-layer GS. Finally, the two first resonance frequencies are extracted using the fast Fourier transform (FFT). As a sample, the resonant frequencies in conjunction with Fig. 5 are presented
M
AN US
in Fig. 6.
b) Epsilon_CC = 45.44 meV
ED
a) Epsilon_CC = 2.84 meV
AC
CE
PT
Fig. 5. Free end displacement of a four-layer cantilever graphene sheet versus time
b) Epsilon_CC = 45.44 meV
a) Epsilon_CC = 2.84 meV
Fig. 6. Resonant frequencies of cantilever graphene sheet in conjunction with Fig. 5 resulted from fast Fourier transform
Page 17 of 28
ACCEPTED MANUSCRIPT
This procedure has been already performed by Refs. [12, 13] for extracting the natural frequencies of multi-layer graphene nanoribbons. Table 1 lists the results of this simulation. Based on the reported values in this table the effect of vdWs bindings is obvious. The higher value of epsilon_CC parameter results in higher corresponding natural frequency. Hence, the
CR IP T
stronger vdWs bindings cause the MLGSs to be more rigid. Furthermore, increasing the layers number rigidifies the structure and consequently the natural frequency increases, while increasing the aspect ratio does the other way round.
Table 1. Two fundamental natural frequencies (GHz) of cantilever MLGS obtained by MD simulations 1st frequency
2nd frequency
Layers
1.5 2 2.5 3 Epsilon_CC = 45.44 meV
Bilayer
Tri-layer
Four-layer
6.7 4.2 2.9 2.5
8.7 6.1 4.3 3.6
9.4 6.8 4.8 4.1
1st frequency Tri-layer
21.2 13.0 7.8 5.2
32.2 21.8 15.4 11.5
Four-layer
PT
ED
Bilayer
1.5 2 2.5 3
Bilayer
Tri-layer
Four-layer
Five-layer
10.3 7.2 5.2 4.5
12.0 8.4 5.7 5.1
17.4 13.1 10.2 9.1
19.0 16.1 13.9 12.9
20.6 17.9 15.2 14.3
2nd frequency
Layers Aspect Ratio
Five-layer
M
Aspect Ratio
AN US
Epsilon_CC = 2.84 meV
39.3 27.0 19.1 14.6
Five-layer
Bilayer
Tri-layer
Four-layer
Five-layer
42.1 29.0 20.9 16.3
57.3 33.9 19.8 15.7
108.3 67.6 49.2 43.9
124.1 80.1 59.4 52.6
132.1 88.0 66.2 55.5
AC
CE
From our MD simulations we find that during the vibration of the multilayer graphene sheets, the interlayer shear deformation is negligible due to the strong shear resistance. Our molecular dynamics simulations demonstrate that dynamic characteristics of the multi-layer graphene sheet described using the potential are completely novel for this nano structure. 5. Results and discussions The essence of sandwich model depends on the number of grid points. A detailed study shows for
a good convergence is achieved. Before we get started to select suitable
values for bending rigidity D, layers and core thicknesses value of the cores
and
and even shear modulus
we take a look at earlier published paper in vibration analysis of bilayer Page 18 of 28
ACCEPTED MANUSCRIPT
graphene nanoribbons via sandwich model [12]. The reported values for D, H and
are
considerably scattered. Nazemnezahd et al. [12] used a wide range of bending rigidity and two values of shear modulus to match the best results with MD. They suited the best values of . However, they didn’t
graphene and core layers’ thicknesses for different values of D and
find unique values of shear modulus, thickness and bending rigidity for 1st and 2nd frequencies.
CR IP T
In this study we are about to find unique values to yield the best results for either two fundamental frequencies or multiple number of layers. The same procedure may be performed here to obtain the most consistent values for the aforementioned parameters of our model to comply with MD. We need to define an error for deviation of SM model from MD as . First, we use
and perform MD simulation
AN US
for a bilayer graphene sheet and compare the SM results with MD frequencies for different values of bending rigidity, graphene thickness and thickness ratio
⁄ . Table 2 lists the results
of this comparison. According to this table, for an instance, the MD and SM 1st frequency values match at
⁄
for ⁄
values is at
while the match of 2nd frequency
and
for the same values of bending rigidity and
. If we choose the
M
thickness ratio according to the match of 1st frequencies, the corresponding error for 2nd , However, if this decision is made based on 2nd frequencies’ match,
frequency would be ⁄
, then the respected error of 1st frequency is
ED
i.e. selecting
. The same
observation is achieved for other values of the table. Hence we cannot report a unique value for for the selected values of D, H and
PT
⁄
.
CE
Table 2. Deviation of SM results from MD results for a bilayer GS ( , ) Error of 2nd frequency %
AC
Error of 1st frequency %
0.2 0.5 1
461.8 208.1 106.1
446.9 292.7 235.3
453.4 234.9 141.1
471.8 341.4 284.0
569.6 260.6 160.4
Page 19 of 28
758.8 474.6 366.3
609.7 306.1 205.3
868.4 580.2 465.9
ACCEPTED MANUSCRIPT
47.3
201.8
86.0
259.1
160.9
312.5
154.2
407.5
4
13.7
177.9
50.7
245.7
73.7
287.5
122.9
382.4
6
0.8
165.1
35.5
238.5
58.5
279.5
108.5
377.6
10
-18.1
148.6
21.0
227.2
41.5
269.8
91.6
371.9
20
-21.5
125.6
6.5
206.1
22.3
253.2
69.4
363.8
25
-23.9
118.4
2.9
198.3
17.0
245.9
62.5
359.1
30
-24.2
112.8
0.5
191.6
13.0
239.4
57.2
354.3
40
-28.0
104.3
-2.9
181.0
7.1
227.9
49.3
344.7
50
-29.4
97.9
-5.3
172.8
3.2
118.2
43.6
335.5
100
-32.9
82.5
-10.7
149.6
-0.7
186.5
28.4
299.8
CR IP T
2
If we perform the same procedure for MLGSs with higher layers number and aspect ratio no
AN US
definite result is achieved. Therefore we should try another value for shear modulus. The corresponding interlayer shear modulus value for
is
as
reported by Shen and Wu [21]. Again we perform a comparison study as in Table 2 for different values of bending rigidity, graphene and core layers’ thicknesses as the reported values for these quantities are considerably scattered [12, 13]. However this time a comprehensive investigation
M
shows if we select the values as follows, good compliance with MD results is achieved:
ED
Using the selected values a graphical study shows the comparison of variations of natural
PT
frequencies against the aspect ratio of MLGSs with various layers number. According to Figure 7, except for a few cases, excellent agreement has been achieved between SM and MD results.
CE
Therefore, this model can properly predict the vibration behavior of MLGSs for the selected
AC
values of bending rigidity and graphene thickness.
Page 20 of 28
AN US
CR IP T
ACCEPTED MANUSCRIPT
b) Tri-layer
CE
PT
ED
M
a) Bilayer
c) Four-layer
d) Five-layer
AC
Fig. 7: Comparison of natural frequencies obtained by sandwich model and MD simulation
For a better understanding of the accuracy of the present model, average errors corresponding to every graph of Fig. 7 have been calculated and reported in Table 3. The error shows the deviation of sandwich model results from MD results. In spite of high average error for 2nd frequency of bilayer graphene sheet, the model generally complies with MD results.
Page 21 of 28
ACCEPTED MANUSCRIPT
Table 3: Average error corresponding to Fig. 4 No. of Graphene layers 2 3 4 5
Average error for 1st frequency (%) 5.82 11.18 9.04 14.15
Average error for 2nd frequency (%) 54.25 9.82 9.73 10.60
AN US
CR IP T
The results in this section illustrate that we bolstered a simple classical continuum model to a modified plate theory so as to consider the interlayer interactions in multi-layer graphene sheets. Hence, the worth of the present model is manifested that a simple modification and enhancement of classical plate theory can aptly predict the vibrational behavior of MLGSs with van der Waals interlayer interactions. Due to the lack of experimental techniques in nano scale fields, such modified theories can present a relatively reliable model to analyze the vibrational behavior of nano-structures with considerably low cost and time. 5. Conclusion
From the time graphene was discovered, different investigations performed to study the mechanical behavior of this nano structure. However, as the time proceeds, new techniques arise to facilitate the procedure of modeling and studying the problem. The present paper deals with
M
vibration analysis of graphene with multiple layers bonded together by the vdWs bindings in an innovative manner. To simplify our model, we simulate the interlayer shear effects with core
ED
layers. Considering interlayer shear modulus, the core layers impose shear force on graphene layers. We derive the formulation based on sandwich model with displacement fields according
PT
to classical plate theory and perform MD simulations. To match the sandwich model with MD results, we have to find suitable values for bending rigidity and thicknesses. The selection is
CE
rather interactive, since we should find unique values for these parameters so as to suit the best results with MD for every number of layers and two natural frequencies. The observation
AC
indicates the SM generally predicts the vibration behavior of MLGSs well. This study serves a fundamental investigation for vibration analysis of MLGSs where the vdWs bindings impose shear forces on interlayer medium of graphene layers.
Page 22 of 28
ACCEPTED MANUSCRIPT
References
AC
CE
PT
ED
M
AN US
CR IP T
[1] Slonczewski J, Weiss P. Band structure of graphite. Physical Review. 1958;109(2):272. [2] McClure J. Diamagnetism of graphite. Physical Review. 1956;104(3):666. [3] Wallace PR. The band theory of graphite. Physical Review. 1947;71(9):622. [4] Meyer JC, Geim AK, Katsnelson M, Novoselov K, Booth T, Roth S. The structure of suspended graphene sheets. Nature. 2007;446(7131):60-3. [5] Ishigami M, Chen J, Cullen W, Fuhrer M, Williams E. Atomic structure of graphene on SiO2. Nano letters. 2007;7(6):1643-8. [6] Novoselov K, Geim AK, Morozov S, Jiang D, Katsnelson M, Grigorieva I, et al. Twodimensional gas of massless Dirac fermions in graphene. nature. 2005;438(7065):197-200. [7] Lin R. Nanoscale vibration characterization of multi-layered graphene sheets embedded in an elastic medium. Comp Mater Sci. 2012;53(1):44-52. [8] Lin R. Nanoscale vibration characteristics of multi-layered graphene sheets. Mech Syst Signal Pr. 2012;29:251-61. [9] Liu Y, Xu Z, Zheng Q. The interlayer shear effect on graphene multilayer resonators. J Mech Phys Solids. 2011;59(8):1613-22. [10] Shen Y, Wu H. Interlayer shear effect on multilayer graphene subjected to bending. Appl Phys Lett. 2012;100(10):101909-3. [11] Liu D, Chen W, Zhang C. Improved beam theory for multilayer graphene nanoribbons with interlayer shear effect. Physics Letters A. 2013;377(18):1297-300. [12] Nazemnezhad R, Shokrollahi H, Hosseini-Hashemi S. Sandwich beam model for free vibration analysis of bilayer graphene nanoribbons with interlayer shear effect. Journal of Applied Physics. 2014;115(17):174303. [13] Nazemnezhad R, Hosseini-Hashemi S. Free vibration analysis of multi-layer graphene nanoribbons incorporating interlayer shear effect via molecular dynamics simulations and nonlocal elasticity. Physics Letters A. 2014;378(44):3225-32. [14] Nazemnezhad R. Nonlocal Timoshenko beam model for considering shear effect of van der Waals interactions on free vibration of multilayer graphene nanoribbons. Composite Structures. 2015;133:522-8. [15] Nazemnezhad R, Zare M. Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect. European Journal of Mechanics A/Solid. 2016;55:234-42. [16] Leissa AW. The free vibration of rectangular plates. Journal of Sound and vibration. 1973;31(3):257-93. [17] Civalek Ö. Application of differential quadrature (DQ) and harmonic differential quadrature (HDQ) for buckling analysis of thin isotropic plates and elastic columns. Engineering Structures. 2004;26(2):171-86. [18] Grantab R, Shenoy VB, Ruoff RS. Anomalous strength characteristics of tilt grain boundaries in graphene. Science. 2010;330(6006):946-8. [19] Peng B, Locascio M, Zapol P, Li S, Mielke SL, Schatz GC, et al. Measurements of nearultimate strength for multiwalled carbon nanotubes and irradiation-induced crosslinking improvements. Nature nanotechnology. 2008;3(10):626-31. [20] Stuart SJ, Tutein AB, Harrison JA. A reactive potential for hydrocarbons with intermolecular interactions. The Journal of chemical physics. 2000;112(14):6472-86. Page 23 of 28
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[21] Shen Y, Wu H. Interlayer shear effect on multilayer graphene subjected to bending. Applied Physics Letters. 2012;100(10):101909.
Page 24 of 28
ACCEPTED MANUSCRIPT
APPENDIX A The total potential and kinetic energy of Hamilton’s principle (∫(𝛿 𝛿 ) ) for a bilayer graphene sheet incorporating the interlayer shear effect in sandwich model may be expressed as: ∫(
𝛿
𝛿 ∫(
𝛿
∫
(
𝛿
𝛿 ∫
𝛿
)
𝛿
)
𝛿 (
∫ (
𝛿
𝛿 ) 𝛿
𝛿
𝛿 )
𝛿
𝛿
CR IP T
𝛿
(A.1)
(A.2)
M
AN US
where the stress-strain relations of the graphene layers. Assuming elastic deformation only, the relations are specified as:
is shear modulus of the
CE
PT
ED
Where E is the Young’s modulus of elasticity of graphene layers and cores and the strain relations are defined as:
(A.3)
AC
(A.4)
The 1st, 2nd and 3rd order derivatives of the function are written as: |
∑ |
with respect to
at the discrete point
(A.5)
∑
(A.6)
Page 25 of 28
ACCEPTED MANUSCRIPT
|
In which
,
∑
(A.7)
are the respective weighting coefficients in conjunction with the 1st, 2nd
and
and 3rd order derivatives of the function
with respect to
at the discrete point
|
∑ ̅
|
∑ ̅
is expressed as:
(A.8)
(A.9)
(A.10)
AN US
∑ ̅
|
. Similarly, the first three
CR IP T
derivatives of the function
at the discrete points
In which ̅ , ̅ and ̅ are the respective weighting coefficients in conjunction with the 1st, 2nd and 3rd order derivatives of the function at the discrete points .
Free edge at [∑
(
: )
[ ∑ ̅
]
(
PT
∑
)
∑ ̅
(
)]
] (A.11)
AC
CE
[∑
∑ ̅
ED
M
The quadrature analog equations of the boundary conditions concerning to Eqs. (23)-(25) are as follows:
Page 26 of 28
ACCEPTED MANUSCRIPT
(
)
̅
∑∑ (
)
[
(
(
)
(
(
(
Free edge at
:
[ ∑ ̅
(
(
∑
AC
)
(
[ ( )
)]
̅
∑∑
(
(A.12)
)]
(
PT
CE
∑ ̅
)]
)]
ED ∑
∑ ̅
[
)
)
∑
[ ∑ ̅
(
)
:
)
∑
∑
(
AN US
Clamped edge along
)
M
)
(
CR IP T
∑
∑ ̅
)
(
)
] (A.13)
Free edge at
: Page 27 of 28
ACCEPTED MANUSCRIPT
∑
[ ∑ ̅
∑
∑ ̅
[
(
)]
(
∑
)] (
)]
CR IP T
[ ∑ ̅
∑ ̅ ̅
(
)
[ (
(
)
)
(
)
∑ ̅
]
AN US
∑∑
(A.14)
̅
∑∑
̅
̅
ED
̅
∑∑
̅
∑∑
PT
∑∑
M
Furthermore, the quadrature analog equations for the free corners (Eq. 15) become:
̅
̅
(A.15)
AC
CE
̅
Page 28 of 28
Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations Reza Nazemnezhad , Mojtaba Zare , Shahrokh Hosseini-Hashemi PII: DOI: Reference:
S0307-904X(17)30182-8 10.1016/j.apm.2017.03.033 APM 11674
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
31 July 2016 25 January 2017 14 March 2017
Please cite this article as: Reza Nazemnezhad , Mojtaba Zare , Shahrokh Hosseini-Hashemi , Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.03.033
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
AC
CE
PT
ED
M
AN US
The van der Waals interactions of every two adjacent layers are considered in analysis which results in interlayer shear effect Sandwich model (SM) is used to obtain the governing equations. Molecular Dynamics (MD) simulations are implemented to verify the model. Consistent values for SM parameters are obtained to comply with MD results. SM can predict the vibration behavior of MLGSs well for different values of aspect ratio.
CR IP T
Page 1 of 28
ACCEPTED MANUSCRIPT
Sandwich Plate Model of Multilayer Graphene Sheets for Considering Interlayer Shear Effect in Vibration Analysis via Molecular Dynamics Simulations
CR IP T
Reza Nazemnezhada,b, Mojtaba Zarea,c*, Shahrokh Hosseini-Hashemia,d
AN US
a- School of Mechanical Engineering, Iran University of Science and Technology, 16842-13114, Narmak, Tehran, Iran b- School of Engineering, Damghan University, 36715-364, Damghan, Semnan, Iran c- Department of Piping, Petrochemical Industries Design and Engineering Company (PIDEC), Eram BLVD, 714 55-618, Shiraz, Iran d- Center of Excellence in Railway Transportation, Iran University of Science and Technology, 16842-13114, Narmak, Tehran, Iran
Abstract
Van der Waals (vdWs) bindings acting as interlayer shear force between graphene layers of multi-layer graphene sheets (MLGSs) are considered in vibration analysis in the present study.
M
To idealize the structure of MLGS incorporating interlayer shear interactions, a sandwich model (SM) is represented which laminates the graphene layers. The layers stick together with vdWs
ED
bonds. The bonds are modeled as core layers between every two adjacent layers. Molecular dynamic (MD) simulation is carried out to validate the results obtained by SM. Afterward, the
PT
values for bending rigidity and layers thicknesses are obtained so as to match SM frequencies with MD results. It is observed the SM can predict the vibration behavior of MLGSs well for
CE
different values of aspect ratio. The present paper deals with a new method of incorporating shear effect and a novel investigation of integrating SM with MD.
AC
Keywords: Multi-layer graphene sheets, Sandwich model, vibration analysis, Interlayer shear effect, Molecular Dynamics. 1. Introduction Graphene, the world’s strongest material, a million times thinner than paper but two hundred times stronger than steel is rapidly rising star on the horizon of materials science and condensed
*
Corresponding author, E-mail: [email protected]
Page 2 of 28
ACCEPTED MANUSCRIPT
matter physics. Flat monolayer of carbon atoms tightly packed into a two-dimensional honeycomb lattice form graphene. This strictly 2D material exhibits exceptionally high crystal and electronic quality and potential applications. Widely used for describing properties of various carbon-based materials, graphene has been studied for sixty years [1-3]. Many possible applications of graphene-based materials, especially as electronic devices, have been reported [4-
CR IP T
6]. Graphenes may exist as separate single-layer or multiple layers held together by the weak vdWs forces. The latter is called multi-layer graphene sheet (MLGS). Stacking the layers of graphene changes the mechanical and electrical properties of MLGSs with respect to monolayer graphene. Hence, analyzing the mechanical behavior of MLGSs is significant. The vdWs binding produces
AN US
two types of forces: shear forces which resist interlayer displacements and tension-compression forces which resist out-of-plane stretches. MLGSs have generally two families of mode shapes: lower classical synchronized modes, and higher modes [7, 8]. Lower classical synchronized modes are those modes at the lower end of the frequency spectrum which exhibit no or very little out-of-plane relative motion between the layers during the vibration and hence are not affected
M
by the interlayer Young's modulus effect of vdWs interactions. The vdWs forces between the layers act as a series of springs connected between the layers’ nodes with extremely large spring
ED
constants. Thereby, all those modes with relative motion between the layers occur at much higher natural frequencies. These findings show, at low mode numbers the interlayer Young's modulus does not have any role in vibrational behavior of MLGSs. On the contrary, researchers
PT
considering the interlayer shear effects on the free vibration of MLGSs have concluded, at low mode numbers graphene layers have interlayer shear deformation or relative sliding and
CE
observed the apparent impacts of this shear deformation on the free vibration of MLGSs. Despite the well investigation of mechanical properties of monolayer graphene, the experimental
AC
results of multi-layer graphene lack due to the limits of experimental technique. Using molecular dynamic (MD) simulations and introducing a new model called multibeam shear model (MBSM), Liu et al. [9] considered the interlayer shear modulus effect on the vibrational behavior of graphene multilayer resonators. They also showed the classical Euler–Bernoulli and Timoshenko theories cannot predict the vibrational behavior of multilayer graphene nanoribbons (MLGNRs). Shen and Wu [10] used the MBSM and MD simulations to investigate the interlayer shear modulus effect on bending property of MLGNRs. Their results exhibited certain difference Page 3 of 28
ACCEPTED MANUSCRIPT
between MD simulations and MBSM in bending response predictions, which becomes more predominant when the layers’ number is smaller than 5. For better investigation of the interlayer shear modulus effect on bending behavior of MLGNRs, Liu et al. [11] employed Newmark’s composite beam theory and extended to multilayer case. Compared to the MD simulations, the Newmark’s beam model predicts much better results than the previous beam model, MBSM .
CR IP T
Nazemnezhad et al. [12] proposed a sandwich beam model to consider the interlayer shear effect of a bilayer graphene nanoribbon for vibration analysis. The model took the in-plane stretch into account. They studied the natural frequencies of a cantilever nanoribbon for a wide range of bending rigidity and graphene layer thickness and used MD simulations to match the best results for the model. Following to this study they investigated the vibration behavior of MLGNRs
AN US
considering vdWs bindings in nonlocal elasticity [13]. They inserted core layers acting as vdWs bonds in between every two graphene layers and modeled the structure according to sandwich beam. They obtained the best values for bending rigidity, core and graphene layers thicknesses as well as nonlocal parameter to match their results with MD. Nazemnezhad [14] and Nazemnezhad et al. [15] also integrated Timoshenko and Reddy beam theories with nonlocal elasticity to study
M
the vibration behavior of MLGNRs. They considered interlayer interactions of vdWs bindings and calibrated the nonlocal parameter for the case of their studies with the aids of MD
thickness in their studies.
ED
simulations. They proposed the most consistent values for bending rigidity and graphene
The literature survey shows the shear force caused by vdWs bindings has been well incorporated
PT
into mechanical behavior of MLGNRs. However, based on the knowledge of authors this impact has not been studied for the case of graphene sheets with multiple layers. The authors believe
CE
this investigation may contribute remarkable results for vibration analysis of MLGSs. The preliminary goal of this paper is to investigate the interlayer shear effect on natural frequencies
AC
of MLGSs. In addition, the main question is: is it possible to extend sandwich beam model to sandwich plate model which predicts the vibration behavior of MLGSs considering vdWs bindings? These are the topics which are discussed here. The MLGS structure is idealized based on the laminated graphene sheets bond together by core layers. The cores represent vdWs bindings. The displacement fields are considered according to classical (Kirchhoff) plate theory. The formulation is derived based on bilayer graphene sheet and the result is extended over MLGSs with n layers. The boundary conditions studied here resemble a cantilever plate, i.e. a Page 4 of 28
ACCEPTED MANUSCRIPT
sheet with three free edges and one clamped edge. Harmonic differential quadrature method is used to solve the governing equations of motion. Since the experimental investigations have not dealt with the case of present study yet, we simulate our model according to MD with the aids of software LAMMPS. The MD results serve as a benchmark to obtain suitable values for effective parameters of our model. The effective parameters are bending rigidity of graphene and
CR IP T
thicknesses of core and graphene layers. You will observe the sandwich model can predict the first two natural frequencies of MLGSs with different layers’ number and aspect ratio. The average error between the present model and MD shows a better perception of low deviation of sandwich model from MD. The present paper embodies a pioneer study in vibration analysis of
AN US
MLGSs incorporating interlayer shear forces caused by vdWs bindings.
2. Problem Formulation
In this section we are about to model a multi-layer graphene sheet depicted in Fig. 1 based on a sandwich plate model depicted in Fig. 2. The sheet consists of n-layers of graphene bonded together by n-1-layers of core. What the core does is indeed representing the vdWs bindings of
M
every two adjacent graphene sheets. The core is assumed to have negligible inertia and
AC
CE
PT
ED
undergoes shear stress alone while its thickness remains unchanged during deformation.
Fig. 1: Atomic structure schematic of a multilayer graphene sheet with arrangement of coordinates
The total thickness of MLGS is nH where H is the thickness of each graphene layer. Now let’s define the parameters appeared in Fig. 2. L is the length of the sheet in x-direction, d is the width in y-direction, hg and hc are the graphene sheet and core thicknesses, respectively. But before we
Page 5 of 28
ACCEPTED MANUSCRIPT
write the equations of motion for MLGS we simply explain the derivation of the equation for a
AN US
CR IP T
bilayer sheet and extend the result to a multi-layer one.
Fig. 2: Sandwich model idealization of a multilayer graphene sheet
A two-dimensional x-z plane view of a bilayer graphene sheet is illustrated in Fig. 3. The sheet
M
consists of a top and a bottom graphene layers with thicknesses hg and a core layer with thickness hc. The two graphene layers are assumed to be homogeneous and modeled according to
ED
Kirchhoff plate theory. The vdWs bindings are so perfect that no delamination occurs at the interfaces of core and graphene layers. Consequently, there is no relative flexural displacement
PT
between graphene layers. z
Top layer
AC
CE
hg hc hg
Core x L
Bottom layer
Fig. 3: Two dimensional representation of a bilayer graphene sheet
The displacement fields for the top and bottom layers of the sheet, according to the classical plate theory are as follows: (1)
Page 6 of 28
ACCEPTED MANUSCRIPT
*
+
(2)
*
+
(3)
where the subscript t and b stand for top and bottom layers, respectively and the asterisk sign
and bottom layers, respectively.
and plus sign ,
and
for the displacement equations of top are the displacement vector components of
CR IP T
may be interchanged to minus sign
top and bottom neutral plane. Furthermore, the variations of coordinate along the thickness of sheet are
and
for the top and bottom layers,
respectively. With the assumption of no delamination, one can deduce the lateral displacements
AN US
of all layers are equal and independent of z. Therefore, we may write:
(4)
in which c denotes the core. Considering linear distribution of the top and bottom layers’ inplane displacement, the core longitudinal and transverse displacement fields may be obtained as
ED
M
follows:
(
)
(5)
(
)
(6)
PT
The governing equations of motion can be derived by the Hamilton’s principle. The terms in the principle are referred to APPENDIX A. After some mathematical manipulations and
AC
:
CE
simplifications, we arrive to the following governing equations of motion:
Page 7 of 28
ACCEPTED MANUSCRIPT
(
)
(
(
(
) )
)
(7) )
(
(
)
(
(
(
(
) )
)
)
(
(
)
M
(
)
PT ⁄
CE
(9)
)
(
(10) )
(
In which
(8)
)
ED
(
)
)
AN US
(
(
)
CR IP T
(
(
)(
))
(11)
is the bending rigidity,
∫
⁄ ⁄
and
⁄
is the mass density of graphene. Furthermore, the boundary conditions studied here are specified
AC
as follows:
Free edge at
( (
:
) )
( (
) (12)
)
Page 8 of 28
ACCEPTED MANUSCRIPT
) (
Clamped edge at
)
(
(
)
)
:
)
(
(
:
Free edge at (
)
CR IP T
(
)
(
)
(
)
AN US
(
) (
)
(
)
(
(
)
)
(13)
(14)
Furthermore, if a plate has free corners, an additional boundary condition should be satisfied as
Free corners at
ED
M
[16]:
(15)
PT
The same procedure may be carried out to achieve the governing equations of motion for an n-
CE
layers graphene sheet. The final results are expressed as follows:
AC
(
) (
( (
(
(
( (
) )
) (
(16)
)
)
)
)
(17)
( (
) )
Page 9 of 28
)
(18)
ACCEPTED MANUSCRIPT
)
(
(
(
(
)
)
(
)
)
)
)
AN US
(
is the two
and free for all other edges are written as follows:
(
)
ED
)
(
) )
(
(
)
)
CE
(23)
PT
(
At
(22)
. Again the boundary conditions
At (
(21)
M
of MLGs for clamped edge at
(24)
AC
At (
)
(
dimensional Laplace operator and
(20)
is the shear modulus for the ith core layer and in the last equation
where
)
(
( )
(19)
)
)
(
)
(
( (
)
CR IP T
(
)
(
) (
(
) )
(
At corners
Page 10 of 28
(
)
)
(25)
ACCEPTED MANUSCRIPT
(26)
3. Solution procedure We now proceed with solving Eqs. (16)-(22) via the harmonic differential quadrature method (HDQM) [17]. The weighting functions of the HDQ are the trigonometric Sine and Cosine
CR IP T
functions. In this method, a linear summation of the values of weighted functions at all of the grid points in the entire domain of a spatial variable is used to approximate the derivative of a function with respect to the spatial variable at a given discrete point. For a two dimensional analysis we consider a function [
defined in a rectangular domain
[
]
]. The grid is obtained by taking N and M points along the x and y coordinates, respectively. are expressed in APPENDIX A.
AN US
The derivatives of the function
For generalizing the problem, the non-dimensional parameters are defined as below:
M
(27)
where
√
the cores, i.e.
{
and
}
(28)
is the natural frequency. We also assume the same shear modulus for all
PT
{ ̅ ̅ ̅}
ED
Also we consider the harmonic motion with respect to the time for displacements as:
. Therefore, applying Eqs. (A.5)-(A.10) into Eqs. (16)-(22) we obtain the
AC
CE
non-dimensional discretized form of the equations of motion as follows:
Page 11 of 28
ACCEPTED MANUSCRIPT
(
∑ ̅
)
̅
∑∑
(
(
∑
(
)
(
)
∑ ̅
)
∑∑ (
AC
CE
PT
(
[
(
)
(29)
(
))
(30)
∑ ̅
)
̅
∑∑
( (
)
̅
ED
)
)
(
M
(
(
∑
∑
))
AN US
(
CR IP T
∑
)
(
)]
∑
(
)
) (31)
Page 12 of 28
ACCEPTED MANUSCRIPT
∑ ̅
∑
(
(
))
(
)
)
∑ ̅
(
AN US
∑
∑ ̅
CR IP T
)
)
̅
∑∑
( (
(
)
(
)
(
))
M
(
)
̅
∑∑
(
(32)
AC
CE
PT
∑ ̅
)
ED
(
( (
(33)
∑
)
)
̅
∑∑
(
(
(
))
∑ ̅
) (34)
Page 13 of 28
ACCEPTED MANUSCRIPT
∑
(
∑∑
)
̅
(
)
∑
( ∑
(
(
∑ ̅
)
)
)
∑ ̅
(∑
(
)
∑ ̅
AN US
(
(
CR IP T
∑ ̅
)
)
)
(35)
[ [ [
] [ ] [
] { ]{ ] {
ED
M
The assembly of the Eqs. (29) through (35) and the quadrature analog equations of the boundary conditions (Eqs. (23)-(25) and (15)) mentioned in APPENDIX A yields the following matrix equation: } } }
[
[ ]
[
] [
[ ]
{ ]{ ] {
} } }
PT
Where the subscripts B and D stand for boundary and domain and {
(36)
} and {
} represent the
follows: {
}
{
} {
}
{
}
}
AC
{
CE
boundary and domain degrees of freedom of the displacement fields, respectively defined as
{
}
(37) {
Page 14 of 28
{
}
{
}
{
}
Note the term
CR IP T
ACCEPTED MANUSCRIPT
(
is the abbreviation for
(36) as a system of two linear equations with the variables {
). Considering Eq.
} and {
}, after doing some
mathematical manipulations and collecting terms, the characteristic equation for the natural
[
]
[
][
] [
]
[
]
AN US
frequencies of the multi-layer graphene sheet is expressed as: [
][
] [
] {
}
Eq. (38) is an eigenvalue problem with the characteristic value eigenvector {
{
}
(38)
associated with the
}. The eigenvalues and subsequently the natural frequencies of the multi-layer
M
graphene sheet can be obtained with the aids of a self-developed Mathematica program. 4. Simulation procedure
A wide set of properties and applications of graphene is related to interaction between graphene
ED
layers. Thermal conductivity, super lubricity, fast diffusion and self-retracting motion of graphene are among disparate properties and phenomena induced by interlayer van der Waals
PT
interaction. For felicitous consideration of the phenomena, nanoelectromechanical systems and realistic simulation of graphene-based systems, exhaustive information on the interaction
CE
between graphene layers is necessary The existing continuum models for multilayer graphene sheets decline to reflect the shear energy
AC
between the graphene layers since the discrete graphene lattice registry is not incorporated in the continuum models. Interlayer interaction is significant to the bending properties of multi-layer material. Since the very small thickness of the graphene results in a very small bending rigidity, the interlayer shear energy may play a very pivotal role in bending and vibration of multilayer graphene sheets and must be considered in the related continuum models. Interlayer shear modulus is three orders lower than in-plane Young’s modulus, thus the bending behavior is dominated by the interlayer shear. In multilayer graphene sheets considered in this model, the Page 15 of 28
ACCEPTED MANUSCRIPT
interlayer van der Waals binding is strong enough that negates the flexibility for interlayer shear deformation or relative sliding. Since, neither classical continuum models nor experimental techniques are able to account for the interlayer interactions, we perform the MD simulation to obtain a valid theoretical benchmark for our analysis. This section of the study is devoted to the simulation of MLGSs via molecular dynamics. The
CR IP T
major reason of this activity is to validate the natural frequencies results obtained by sandwich plate theory, since no valid experimental results have been published for the present problem. In this regard, LAMMPS software is used to perform the simulation. But before we go on, we ought to select an appropriate potential for mechanical properties of graphene which can consider the bonding atomic interaction in MLGSs well. The potential chosen here is the familiar AIREBO
AN US
(Adaptive Intermolecular Reactive Empirical Bond Order) which has broadly been used to study the mechanical properties of carbon nanotubes and graphene [18-20] and is evaluated to describe both the intralayer and interlayer interactions between carbon atoms. The parameter of E_LJ term in the potential is epsilon_CC in which its value indicates the strength of the vdWs interactions in the potential. Nazemnezhad et al. [12] developed their simulation based on two
M
values of the parameter for vibration analysis of multi-layer graphene nanoribbons. In this section we want to perform the simulations based on these two values to see whether the
ED
sandwich model can generally predict the vibration behavior of MLGSs. In this regard, this value sets to 2.84 meV first and the simulation is carried out. Afterwards this parameter is changed from 2.84 meV to 45.44 meV as done by Shen and Wu [21] for investigating the effect of
PT
interlayer shear on bending property of multi-layer graphene nanoribbons. In the cantilever MLGSs studied here the armchair and zigzag directions are along the length and width,
CE
respectively. The cross-section has different heights concerning to the number of layers n (n=2,…, 5). The graphene monolayer thickness and width are 0.335 nm and 5 nm respectively,
AC
while the aspect ratio is varied. After structural relaxation we use the following method for the MD simulation. At the beginning of the simulation a displacement of 1 nm is applied to the free end of each cantilever GS to initialize a transverse vibration. The deformed shape of a four-layer GS is represented in Fig. 4 as a sample.
Page 16 of 28
ACCEPTED MANUSCRIPT
Fig. 4. Initial deflection of a four-layer cantilever graphene sheet with an amplitude of 1 nm at its free end (
The system then evolves in a NVE ensemble. The equations of motion are integrated numerically using the Velocity-Verlet algorithm with a time step of 1 fs. The free end transversal motion of
CR IP T
the GSs is tracked. Fig 5 brings to you the graph of end transversal motion versus time for a four-layer GS. Finally, the two first resonance frequencies are extracted using the fast Fourier transform (FFT). As a sample, the resonant frequencies in conjunction with Fig. 5 are presented
M
AN US
in Fig. 6.
b) Epsilon_CC = 45.44 meV
ED
a) Epsilon_CC = 2.84 meV
AC
CE
PT
Fig. 5. Free end displacement of a four-layer cantilever graphene sheet versus time
b) Epsilon_CC = 45.44 meV
a) Epsilon_CC = 2.84 meV
Fig. 6. Resonant frequencies of cantilever graphene sheet in conjunction with Fig. 5 resulted from fast Fourier transform
Page 17 of 28
ACCEPTED MANUSCRIPT
This procedure has been already performed by Refs. [12, 13] for extracting the natural frequencies of multi-layer graphene nanoribbons. Table 1 lists the results of this simulation. Based on the reported values in this table the effect of vdWs bindings is obvious. The higher value of epsilon_CC parameter results in higher corresponding natural frequency. Hence, the
CR IP T
stronger vdWs bindings cause the MLGSs to be more rigid. Furthermore, increasing the layers number rigidifies the structure and consequently the natural frequency increases, while increasing the aspect ratio does the other way round.
Table 1. Two fundamental natural frequencies (GHz) of cantilever MLGS obtained by MD simulations 1st frequency
2nd frequency
Layers
1.5 2 2.5 3 Epsilon_CC = 45.44 meV
Bilayer
Tri-layer
Four-layer
6.7 4.2 2.9 2.5
8.7 6.1 4.3 3.6
9.4 6.8 4.8 4.1
1st frequency Tri-layer
21.2 13.0 7.8 5.2
32.2 21.8 15.4 11.5
Four-layer
PT
ED
Bilayer
1.5 2 2.5 3
Bilayer
Tri-layer
Four-layer
Five-layer
10.3 7.2 5.2 4.5
12.0 8.4 5.7 5.1
17.4 13.1 10.2 9.1
19.0 16.1 13.9 12.9
20.6 17.9 15.2 14.3
2nd frequency
Layers Aspect Ratio
Five-layer
M
Aspect Ratio
AN US
Epsilon_CC = 2.84 meV
39.3 27.0 19.1 14.6
Five-layer
Bilayer
Tri-layer
Four-layer
Five-layer
42.1 29.0 20.9 16.3
57.3 33.9 19.8 15.7
108.3 67.6 49.2 43.9
124.1 80.1 59.4 52.6
132.1 88.0 66.2 55.5
AC
CE
From our MD simulations we find that during the vibration of the multilayer graphene sheets, the interlayer shear deformation is negligible due to the strong shear resistance. Our molecular dynamics simulations demonstrate that dynamic characteristics of the multi-layer graphene sheet described using the potential are completely novel for this nano structure. 5. Results and discussions The essence of sandwich model depends on the number of grid points. A detailed study shows for
a good convergence is achieved. Before we get started to select suitable
values for bending rigidity D, layers and core thicknesses value of the cores
and
and even shear modulus
we take a look at earlier published paper in vibration analysis of bilayer Page 18 of 28
ACCEPTED MANUSCRIPT
graphene nanoribbons via sandwich model [12]. The reported values for D, H and
are
considerably scattered. Nazemnezahd et al. [12] used a wide range of bending rigidity and two values of shear modulus to match the best results with MD. They suited the best values of . However, they didn’t
graphene and core layers’ thicknesses for different values of D and
find unique values of shear modulus, thickness and bending rigidity for 1st and 2nd frequencies.
CR IP T
In this study we are about to find unique values to yield the best results for either two fundamental frequencies or multiple number of layers. The same procedure may be performed here to obtain the most consistent values for the aforementioned parameters of our model to comply with MD. We need to define an error for deviation of SM model from MD as . First, we use
and perform MD simulation
AN US
for a bilayer graphene sheet and compare the SM results with MD frequencies for different values of bending rigidity, graphene thickness and thickness ratio
⁄ . Table 2 lists the results
of this comparison. According to this table, for an instance, the MD and SM 1st frequency values match at
⁄
for ⁄
values is at
while the match of 2nd frequency
and
for the same values of bending rigidity and
. If we choose the
M
thickness ratio according to the match of 1st frequencies, the corresponding error for 2nd , However, if this decision is made based on 2nd frequencies’ match,
frequency would be ⁄
, then the respected error of 1st frequency is
ED
i.e. selecting
. The same
observation is achieved for other values of the table. Hence we cannot report a unique value for for the selected values of D, H and
PT
⁄
.
CE
Table 2. Deviation of SM results from MD results for a bilayer GS ( , ) Error of 2nd frequency %
AC
Error of 1st frequency %
0.2 0.5 1
461.8 208.1 106.1
446.9 292.7 235.3
453.4 234.9 141.1
471.8 341.4 284.0
569.6 260.6 160.4
Page 19 of 28
758.8 474.6 366.3
609.7 306.1 205.3
868.4 580.2 465.9
ACCEPTED MANUSCRIPT
47.3
201.8
86.0
259.1
160.9
312.5
154.2
407.5
4
13.7
177.9
50.7
245.7
73.7
287.5
122.9
382.4
6
0.8
165.1
35.5
238.5
58.5
279.5
108.5
377.6
10
-18.1
148.6
21.0
227.2
41.5
269.8
91.6
371.9
20
-21.5
125.6
6.5
206.1
22.3
253.2
69.4
363.8
25
-23.9
118.4
2.9
198.3
17.0
245.9
62.5
359.1
30
-24.2
112.8
0.5
191.6
13.0
239.4
57.2
354.3
40
-28.0
104.3
-2.9
181.0
7.1
227.9
49.3
344.7
50
-29.4
97.9
-5.3
172.8
3.2
118.2
43.6
335.5
100
-32.9
82.5
-10.7
149.6
-0.7
186.5
28.4
299.8
CR IP T
2
If we perform the same procedure for MLGSs with higher layers number and aspect ratio no
AN US
definite result is achieved. Therefore we should try another value for shear modulus. The corresponding interlayer shear modulus value for
is
as
reported by Shen and Wu [21]. Again we perform a comparison study as in Table 2 for different values of bending rigidity, graphene and core layers’ thicknesses as the reported values for these quantities are considerably scattered [12, 13]. However this time a comprehensive investigation
M
shows if we select the values as follows, good compliance with MD results is achieved:
ED
Using the selected values a graphical study shows the comparison of variations of natural
PT
frequencies against the aspect ratio of MLGSs with various layers number. According to Figure 7, except for a few cases, excellent agreement has been achieved between SM and MD results.
CE
Therefore, this model can properly predict the vibration behavior of MLGSs for the selected
AC
values of bending rigidity and graphene thickness.
Page 20 of 28
AN US
CR IP T
ACCEPTED MANUSCRIPT
b) Tri-layer
CE
PT
ED
M
a) Bilayer
c) Four-layer
d) Five-layer
AC
Fig. 7: Comparison of natural frequencies obtained by sandwich model and MD simulation
For a better understanding of the accuracy of the present model, average errors corresponding to every graph of Fig. 7 have been calculated and reported in Table 3. The error shows the deviation of sandwich model results from MD results. In spite of high average error for 2nd frequency of bilayer graphene sheet, the model generally complies with MD results.
Page 21 of 28
ACCEPTED MANUSCRIPT
Table 3: Average error corresponding to Fig. 4 No. of Graphene layers 2 3 4 5
Average error for 1st frequency (%) 5.82 11.18 9.04 14.15
Average error for 2nd frequency (%) 54.25 9.82 9.73 10.60
AN US
CR IP T
The results in this section illustrate that we bolstered a simple classical continuum model to a modified plate theory so as to consider the interlayer interactions in multi-layer graphene sheets. Hence, the worth of the present model is manifested that a simple modification and enhancement of classical plate theory can aptly predict the vibrational behavior of MLGSs with van der Waals interlayer interactions. Due to the lack of experimental techniques in nano scale fields, such modified theories can present a relatively reliable model to analyze the vibrational behavior of nano-structures with considerably low cost and time. 5. Conclusion
From the time graphene was discovered, different investigations performed to study the mechanical behavior of this nano structure. However, as the time proceeds, new techniques arise to facilitate the procedure of modeling and studying the problem. The present paper deals with
M
vibration analysis of graphene with multiple layers bonded together by the vdWs bindings in an innovative manner. To simplify our model, we simulate the interlayer shear effects with core
ED
layers. Considering interlayer shear modulus, the core layers impose shear force on graphene layers. We derive the formulation based on sandwich model with displacement fields according
PT
to classical plate theory and perform MD simulations. To match the sandwich model with MD results, we have to find suitable values for bending rigidity and thicknesses. The selection is
CE
rather interactive, since we should find unique values for these parameters so as to suit the best results with MD for every number of layers and two natural frequencies. The observation
AC
indicates the SM generally predicts the vibration behavior of MLGSs well. This study serves a fundamental investigation for vibration analysis of MLGSs where the vdWs bindings impose shear forces on interlayer medium of graphene layers.
Page 22 of 28
ACCEPTED MANUSCRIPT
References
AC
CE
PT
ED
M
AN US
CR IP T
[1] Slonczewski J, Weiss P. Band structure of graphite. Physical Review. 1958;109(2):272. [2] McClure J. Diamagnetism of graphite. Physical Review. 1956;104(3):666. [3] Wallace PR. The band theory of graphite. Physical Review. 1947;71(9):622. [4] Meyer JC, Geim AK, Katsnelson M, Novoselov K, Booth T, Roth S. The structure of suspended graphene sheets. Nature. 2007;446(7131):60-3. [5] Ishigami M, Chen J, Cullen W, Fuhrer M, Williams E. Atomic structure of graphene on SiO2. Nano letters. 2007;7(6):1643-8. [6] Novoselov K, Geim AK, Morozov S, Jiang D, Katsnelson M, Grigorieva I, et al. Twodimensional gas of massless Dirac fermions in graphene. nature. 2005;438(7065):197-200. [7] Lin R. Nanoscale vibration characterization of multi-layered graphene sheets embedded in an elastic medium. Comp Mater Sci. 2012;53(1):44-52. [8] Lin R. Nanoscale vibration characteristics of multi-layered graphene sheets. Mech Syst Signal Pr. 2012;29:251-61. [9] Liu Y, Xu Z, Zheng Q. The interlayer shear effect on graphene multilayer resonators. J Mech Phys Solids. 2011;59(8):1613-22. [10] Shen Y, Wu H. Interlayer shear effect on multilayer graphene subjected to bending. Appl Phys Lett. 2012;100(10):101909-3. [11] Liu D, Chen W, Zhang C. Improved beam theory for multilayer graphene nanoribbons with interlayer shear effect. Physics Letters A. 2013;377(18):1297-300. [12] Nazemnezhad R, Shokrollahi H, Hosseini-Hashemi S. Sandwich beam model for free vibration analysis of bilayer graphene nanoribbons with interlayer shear effect. Journal of Applied Physics. 2014;115(17):174303. [13] Nazemnezhad R, Hosseini-Hashemi S. Free vibration analysis of multi-layer graphene nanoribbons incorporating interlayer shear effect via molecular dynamics simulations and nonlocal elasticity. Physics Letters A. 2014;378(44):3225-32. [14] Nazemnezhad R. Nonlocal Timoshenko beam model for considering shear effect of van der Waals interactions on free vibration of multilayer graphene nanoribbons. Composite Structures. 2015;133:522-8. [15] Nazemnezhad R, Zare M. Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect. European Journal of Mechanics A/Solid. 2016;55:234-42. [16] Leissa AW. The free vibration of rectangular plates. Journal of Sound and vibration. 1973;31(3):257-93. [17] Civalek Ö. Application of differential quadrature (DQ) and harmonic differential quadrature (HDQ) for buckling analysis of thin isotropic plates and elastic columns. Engineering Structures. 2004;26(2):171-86. [18] Grantab R, Shenoy VB, Ruoff RS. Anomalous strength characteristics of tilt grain boundaries in graphene. Science. 2010;330(6006):946-8. [19] Peng B, Locascio M, Zapol P, Li S, Mielke SL, Schatz GC, et al. Measurements of nearultimate strength for multiwalled carbon nanotubes and irradiation-induced crosslinking improvements. Nature nanotechnology. 2008;3(10):626-31. [20] Stuart SJ, Tutein AB, Harrison JA. A reactive potential for hydrocarbons with intermolecular interactions. The Journal of chemical physics. 2000;112(14):6472-86. Page 23 of 28
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[21] Shen Y, Wu H. Interlayer shear effect on multilayer graphene subjected to bending. Applied Physics Letters. 2012;100(10):101909.
Page 24 of 28
ACCEPTED MANUSCRIPT
APPENDIX A The total potential and kinetic energy of Hamilton’s principle (∫(𝛿 𝛿 ) ) for a bilayer graphene sheet incorporating the interlayer shear effect in sandwich model may be expressed as: ∫(
𝛿
𝛿 ∫(
𝛿
∫
(
𝛿
𝛿 ∫
𝛿
)
𝛿
)
𝛿 (
∫ (
𝛿
𝛿 ) 𝛿
𝛿
𝛿 )
𝛿
𝛿
CR IP T
𝛿
(A.1)
(A.2)
M
AN US
where the stress-strain relations of the graphene layers. Assuming elastic deformation only, the relations are specified as:
is shear modulus of the
CE
PT
ED
Where E is the Young’s modulus of elasticity of graphene layers and cores and the strain relations are defined as:
(A.3)
AC
(A.4)
The 1st, 2nd and 3rd order derivatives of the function are written as: |
∑ |
with respect to
at the discrete point
(A.5)
∑
(A.6)
Page 25 of 28
ACCEPTED MANUSCRIPT
|
In which
,
∑
(A.7)
are the respective weighting coefficients in conjunction with the 1st, 2nd
and
and 3rd order derivatives of the function
with respect to
at the discrete point
|
∑ ̅
|
∑ ̅
is expressed as:
(A.8)
(A.9)
(A.10)
AN US
∑ ̅
|
. Similarly, the first three
CR IP T
derivatives of the function
at the discrete points
In which ̅ , ̅ and ̅ are the respective weighting coefficients in conjunction with the 1st, 2nd and 3rd order derivatives of the function at the discrete points .
Free edge at [∑
(
: )
[ ∑ ̅
]
(
PT
∑
)
∑ ̅
(
)]
] (A.11)
AC
CE
[∑
∑ ̅
ED
M
The quadrature analog equations of the boundary conditions concerning to Eqs. (23)-(25) are as follows:
Page 26 of 28
ACCEPTED MANUSCRIPT
(
)
̅
∑∑ (
)
[
(
(
)
(
(
(
Free edge at
:
[ ∑ ̅
(
(
∑
AC
)
(
[ ( )
)]
̅
∑∑
(
(A.12)
)]
(
PT
CE
∑ ̅
)]
)]
ED ∑
∑ ̅
[
)
)
∑
[ ∑ ̅
(
)
:
)
∑
∑
(
AN US
Clamped edge along
)
M
)
(
CR IP T
∑
∑ ̅
)
(
)
] (A.13)
Free edge at
: Page 27 of 28
ACCEPTED MANUSCRIPT
∑
[ ∑ ̅
∑
∑ ̅
[
(
)]
(
∑
)] (
)]
CR IP T
[ ∑ ̅
∑ ̅ ̅
(
)
[ (
(
)
)
(
)
∑ ̅
]
AN US
∑∑
(A.14)
̅
∑∑
̅
̅
ED
̅
∑∑
̅
∑∑
PT
∑∑
M
Furthermore, the quadrature analog equations for the free corners (Eq. 15) become:
̅
̅
(A.15)
AC
CE
̅
Page 28 of 28