Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect

Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect

Accepted Manuscript Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect Reza Nazemne...

2MB Sizes 0 Downloads 44 Views

Accepted Manuscript Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect Reza Nazemnezhad, Mojtaba Zare PII:

S0997-7538(15)00127-8

DOI:

10.1016/j.euromechsol.2015.09.006

Reference:

EJMSOL 3236

To appear in:

European Journal of Mechanics / A Solids

Received Date: 26 April 2015 Revised Date:

17 July 2015

Accepted Date: 18 September 2015

Please cite this article as: Nazemnezhad, R., Zare, M., Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect, European Journal of Mechanics / A Solids (2015), doi: 10.1016/j.euromechsol.2015.09.006. 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

Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect

a- School of Engineering, Damghan University, 36716-41167, Damghan, Iran

RI PT

Reza Nazemnezhada*, Mojtaba Zareb,c

b- School of Mechanical Engineering, Iran University of Science and Technology, 16842-13114, Narmak, Tehran, Iran

c- Department of Piping, Petrochemical Industries Design and Engineering Company (PIDEC), Eram BLVD, 714 55-618, Shiraz, Iran

SC

Abstract

The aim of this study is to investigate the interlayer shear effect on nonlocal free vibration of multilayer

M AN U

graphene nanoribbons (MLGNRs) based on Reddy beam theory. The major novelty of the study is using the higher order shear deformation theory of Reddy to take the van der Waals’ interaction of the layers into account and calibrating the nonlocal parameter with the aids of molecular dynamics (MD) simulations to match the best results. Furthermore, the shear correction factor of nonlocal Timoshenko theory is calculated for modeling the MLGNRs. MLGNRs have broad applications in mechanical devices,

TE D

such as resonators, sensors and actuators. Accordingly the present results can be used as a reference for future works in which the interlayer shear effects has substantial influences on the mechanical behavior of MLGNRs. The present study outstands for its modeling simplicity and low computational time and cost as well as excellent accuracy.

AC C

EP

Keywords: Free vibration, Nonlocal Reddy model, Graphene nanoribbon, Interlayer shear effect.

1. Introduction

Graphene has attracted remarkable attention in the field of nanoscience and nanotechnology due to its extraordinary mechanical [1], electrical [2], and thermal [3] properties. Graphenes may exist as single layer or multilayer structures. Graphenes have attracted tremendous attention in both their twodimensional and one-dimensional forms where the latter is obtained by patterning the layer into strips or ribbons [4].

*

Corresponding author. Email address: [email protected]. Tel./Fax: +98 233 522 0414.

Page 1 of 17

ACCEPTED MANUSCRIPT Graphene nanoribbons (GNRs) are finite-width graphene sheets. Multilayer graphene nanoribbons (MLGNRs) are made of single layer nanoribbon held together by the weak van der Waals (vdWs) forces. These weak vdWs forces can produce interlayer Young’s and shear moduli for MLGNRs. A literature survey shows that the interlayer Young’s modulus effects on the mechanical behavior of MLGNRs have

investigating the interlayer shear modulus effects can be found [16-21].

RI PT

been studied using a wide range of experimental and theoretical approaches [5-15], but a few works

Results of the works investigating the interlayer Young’s modulus effects on the free vibration of MLGs have shown that MLG layers have generally two families of mode shapes: lower classical synchronized

SC

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

M AN U

during the vibration and hence are not affected by the interlayer Young’s modulus effect of vdWs interactions. Higher modes occur at much higher frequency due to the physical fact that the vdWs forces between the layers act as a series of springs connected between the layers of nodes and these springs have very large spring constants, thereby enhancing all those modes with relative motions between the layers to much higher natural frequencies. These findings show that at low mode numbers the interlayer

TE D

Young’s modulus does not have any role in vibrational behavior of MLGNRs. On the contrary, researchers considering the interlayer shear effects on free vibration of MLGNRs have concluded that at low mode numbers GNR layers have interlayer shear deformation or relative sliding [16, 17] and

EP

observed the apparent impacts of this shear deformation on the free vibration of MLGNRs. Although the works considering the interlayer shear modulus effects on mechanical behavior of MLGNRs have verified their theoretical results with the aid of molecular dynamics (MD) simulations, there are

AC C

some criticisms on them: doing verification for limited data [18, 19, 21] and needing high computational volume [16, 18, 19]. Furthermore, research findings on measurements of the elastic constants show relatively considerable scatter in the values of geometrical and mechanical properties of MLGNRs. This indicates that the precision of theoretical methods is dependent on the selected values of geometrical and mechanical properties. The mentioned points highlight requirement of considering the interlayer shear effect with more details. The aim of the present study is to analyze the free vibration of MLGNRs incorporating the interlayer shear modulus based on the nonlocal Reddy beam theory. Contrary to theoretical techniques with large

Page 2 of 17

ACCEPTED MANUSCRIPT formulations and computations implemented to study the interlayer shear effect like sandwich beam model, the present work represents a simple model with low computational cost and time. The nonlocal Reddy beam theory is a well-known theory, does not have high computational volume contrary to the sandwich formulations [16, 18] or Newmark’s theory [19].

Moreover it does not require a shear

RI PT

coefficient to account for the variations of the shear stress across the cross section contrary to the Timoshenko theory. In order to analyze the free vibration of MLGNRs, the governing equations and boundary conditions are derived by implementing the Hamilton’s principle, and they are solved by the harmonic differential quadrature method (HDQM). Then, free vibration of MLGNRs is simulated by the

SC

molecular dynamics simulator software, LAMMPS, and the first two natural frequencies are extracted. In the next step, the first two frequencies obtained by the nonlocal Reddy model are compared to the MD

M AN U

results and the appropriate values for MLGNRs geometrical and mechanical properties as well as the calibrated nonlocal parameters are obtained by using the least square fitting procedure. Lastly, since the value of shear coefficient is not yet reported for nanoscale structures modeled based on the Timoshenko theory, the geometrical and mechanical properties as well as the nonlocal parameter values obtained in the previous steps are used to calibrate the shear correction factor for the MLGNRs. Despite the

TE D

simplicity of the present model, it can predict accurately the vibration behavior of MLGNRs.

2. Problem Formulation

In this section, equations of motion of multilayer graphene nanoribbons modeled based on the nonlocal

EP

Reddy and Timoshenko beam theories are obtained.

AC C

2.1 Nonlocal Reddy formulation

According to the third order shear deformation theory of Reddy [22] the displacement fields for a multilayer graphene nanoribbon depicted in Fig. 1 are expressed as: , ,  =  ,  +  , ; −0.5ℎ ≤  ≤ 0.5ℎ , ,  = 0

, ,  = , 

(1)

Page 3 of 17

ACCEPTED MANUSCRIPT in which u, v and w are the displacement variables along the x, y and z axes, respectively, ,  is the

rotation of the normal to mid-plane about the y axis and ,  is the cubic variation of , ,  across the thickness [22]. The equations of motion can be derived based on Hamilton’s Principle. For the problem here the principle

01

/

   02

(

(.),

+(.),

RI PT

is stated as: 1 # % # %   + 2   − !" $ + " $ &' --.- - = 0 # # +(.)* 2



(.)*

(2)

where A=bh is the MLGNR cross section area with thickness h and width b,  is the stress tensor,  is the

SC

strain tensor and ρ=2260 kg.m-3 is the mass density of monolayer GNR. Considering the strain-

displacement relations and substituting the displacement fields into Eq. (2) and collecting the coefficients

#3 #%

#% − 4 − 56% ! % & + 67 ! % &8 = 0 # # #

#9 #%

#% − 3; − 567 ! % & + 6< ! % &8 = 0 # # # #4 #% − =! %& = 0 # #

M AN U

of  ,  and , the following equations of motion are obtained:

(3)

(4)

(5)

(.)*

+(.)*

9=

(.)*

+(.)*

4=

(.)*

 -=

  -=

;=

(.)*

+(.)*

(6)

(7)

(8)

AC C

+(.)*

 -=

EP

3=

TE D

where 6> = ?  > -=, @ = 2,4,6, and:

  % -=

(9)

in which Q and M are the force and moment stress resultants, respectively and P and R are the higherorder stress resultants and appear only in the higher-order shear deformation theories. To consider the small-scale effect, the nonlocal elasticity theory of Eringen [23] is used. In the nonlocal theory, the one dimensional nonlocal operator ℒ is defined as: ℒ = 1−D

#% # %

(10)

Page 4 of 17

ACCEPTED MANUSCRIPT in which µ is the nonlocal parameter. Furthermore, we assume harmonic motion with respect to the time as: E , , , , , F = EG, H, IFJ KLM 0

(11)

where N> is the natural frequency. To generalize the solution, the following non-dimensional terms are introduced:

I  O = ℎ% G ; HQ = H; R = ; G P P

RI PT

O = I

(12)

in which L represents the nanoribbon length along the x axis. Now applying the nonlocal operator, Eq.

become:

O D 6% % % # % HQ D 6% 6) % # % G 6% % P% % U=6% P% N> & % + "6) − N> $ % + ! N − & HQ ST, #R ST, #R ST, > ST,

"67 −

X1 −

O 6% 6) P% % 6% 6 UP% U=6% P% #I O− N> − 3 &G =0 ST, ST, ST, #R

M AN U

+!

(13)

O D 6% 67 % # % HQ D 6% 6V % # % G 6% 67 P% % 6% % UP% N> $ % + "6V − N> $ % + ! N> − 3 & HQ ST, #R ST, #R ST, ST, O 6% 6V P% % 6% 6) UP% U6 P% #I O− % +! N> − 9 &G =0 ST, ST, ST, #R

O #HQ 36 #G O D % #%I P% % O =0 N> Y + + + N I U #R % #R = #R U >

TE D

!6% −

SC

(10), on Eqs. (3)-(5) and using Eqs.(6)-(9), (11) and (12), the dimensionless nonlocal equations of motion

(14) (15)

where Db is the bending rigidity, G interlayer shear modulus of GNRs or equivalent shear modulus of nonlocal Reddy model, and 6% , 67 , 6<  ℎ%

EP

6 , 6) , 6V  =

Furthermore, the boundary conditions of a cantilever MLGNR can be obtained as

AC C

O = HQ = I O =0 Clamped end at X=0 ⇒ G

(16)

O = 9Q = 4Q = 0 Free end at X=1 ⇒ 3

(17)

O , 9Qand 4Q are the stress resultants in terms of the non-dimensional displacement fields where 3

and are obtained by implementing Eq. (12) in Eqs. (6)-(8) as follows: O 6 #HQ 6 #G O = [ \ % ! & + 7 ! &] 3 % P #R Pℎ #R

(18)

O 67 #HQ 6< #G 9Q = [ \ ! & + % ! &] P #R Pℎ #R

(19)

Page 5 of 17

ACCEPTED MANUSCRIPT O O #I 36% G 4Q = U ^= 5HQ + ! &8 + % _ #R ℎ

(20)

2.2 Nonlocal Timoshenko formulation

ρA

∂ 2w ∂ 4w  ∂ψ ∂ 2w  − µρA − kAG + 2  =0  ∂t 2 ∂x 2 ∂t 2  ∂x ∂x 

∂w  ∂ 2ψ ∂ 4ψ ∂ 2ψ  kA G ψ + − µρI − Db =0  + ρI 2 2 2 ∂x  ∂t ∂x ∂t ∂x 2 

(21)

SC

be given as:

RI PT

Based on the nonlocal Timoshenko theory [24, 25], the governing equations of motions of MLGNRs can

(22)

where ρ, A, µ and w are the same variables defined earlier, ψ denotes the rotation of the cross-section, k is

M AN U

the shear correction factor, G interlayer shear modulus of GNRs or equivalent shear modulus of nonlocal Timoshenko model, and I is the second moment of area.

For a cantilever MLGNR, the boundary conditions are defined as at x=0 ⇒w (0 ) = ψ (0 ) = 0 at x=L ⇒ M ( L ) = Q ( L ) = 0

(23) (24)

respectively, and given by M = Db

TE D

in which, M and Q are the nonlocal bending moment resultant and the nonlocal force resultant, ∂ψ ∂ 3ψ ∂ 2w + µρI + µρA ∂x ∂x ∂t 2 ∂t 2

EP

∂w  ∂ 3w  Q = kA G ψ +  + µρA ∂x  ∂x ∂t 2 

Again we assume harmonic motion with respect to the time and consider the non-dimensional terms

AC C

similar to the previous section. Hence the non-dimensional form of Eqs. (21) and (22) will be written as follows: −! −!

O O `=U #Ψ D =N2@ `=U #2 I O =0 + & 2 − + =N2@ PI P P P #R #R

O O D 6 N2@ T, #2 Ψ #Ψ O =0 + % & 2 + `=U + `=U + 6 N2@ Ψ % P P #R #R

Page 6 of 17

(25)

(26)

ACCEPTED MANUSCRIPT 3. Solution Procedure In order to avoid repetitive representations, we solve only the governing equations of nonlocal Reddy model, Eqs. (13)-(15), via the harmonic differential quadrature method (HDQM). The nature of the HDQ method is that the derivative of a function with respect to a spatial variable at a given discrete point is approximated as a linear summation of the values of weighted harmonic Sine and Cosine functions at all

RI PT

of the grid points in the entire domain of that variable. To this end, we consider a function b defined

in the domain 0 ≤  ≤ P discretized by N points along the x coordinate. Then the nth order derivatives of

the function b with respect to x at the point xi can be expressed as: h

SC

-> bK  > = c dKe bfe g ; k = 1,2, … , m; @ = 1,2, … , m − 1 - > eij

(27)

in which dKe are the respective weighting coefficients in conjunction to the nth order derivatives of the

M AN U

>

function b, at the discrete points xi. For more details about the HDQM and the respective weighting O , HQ coefficients one can refer to Ref. [26]. Here, the function b can represents any of the functions G

O . To find the numerical approximation for Eqs. (13)-(15) we apply the HDQ formula (Eq. (27)) to and I

these equations and obtain their discretized forms at any point Xi as follows: h

h

eij

eij

= N>% n

h

h

h

eij

= N>% n

h

h

eij

eij

(28)

h

3U6% =P% 3U6% =P% 9U6) P% j O O RK  HQRK  − c dKe I fRe g − G [ [ [ h

eij

D 67 D 6V P% 67 P% 6V % % O O RK o c dKe HQfRe g + c dKe G fRe g − HQRK  − G [ [ [ [

j % O c dKe HQfRe g + c dKe I fRe g + eij

eij

eij

EP eij

AC C

h

h

D 6% D 6) P% 6% P% 6) % % O O RK o c dKe HQfRe g + c dKe G fRe g − HQRK  − G [ [ [ [

% % O 67 c dKe HQfRe g + 6V c dKe G fRe g − eij

h

U=P% U=P% 3U6 P% j O O RK  HQRK  − c dKe I fRe g − G [ [ [

TE D

% % O 6% c dKe HQfRe g + 6) c dKe G fRe g −

h

eij

h

36 P% D j O % O O RK  + c dKe G fRe g = N>% n− I c dKe I fRe go = U U eij

eij

(29)

(30)

Furthermore, applying the HDQM on the boundary conditions, Eqs. (16) and (17), results in their discretized form as:

O Rj  = I O Rj  = 0 HQRj  = G h

h

eij

eij

(31)

O = 6% c Xdhe HQfRe gY + 6) c Xdhe G O fRe gY = 0 3

(32)

Page 7 of 17

ACCEPTED MANUSCRIPT h

h

eij

eij

O fRe gY = 0 9Q = 67 c Xdhe HQfRe gY + 6V c Xdhe G h

O fRe gY + 4Q = HQRh  + c Xdhe I eij

(33)

36 O Rh  = 0 G =

(34)

RI PT

It should be noted that for the edge X=0 the discrete point is X1 and for X=1 the discrete point is XN. Developing Eqs. (13)-(15) and (31)-(34), the following assembled matrix equations are obtained: q` r q`,s r ET, F ET F 0 0 p ,, tu v = N>% pq3 r q3 rt u , v q`s, r q`ss r ETs F ETs F s, ss

(35)

and domain degrees of freedom (DOF), respectively as follows:

O Rj , HQRj , G O Rj , I O Rh , HQRh , G O Rh  F ET, F = EI

SC

where the subscripts b and d stand for boundary and domain and ET, F and ETs F represent the boundary

O R% , … , I O Rh+j , HQR% , … , HQRh+j , G O R% , … , G O Rh+j  F ETs F = EI

(36) (37)

M AN U

Let’s consider Eq. (35) as a system of two linear equations with the variables ET, F and ETs F. After doing some mathematical simplifications and collecting terms, the characteristic equation for the natural frequencies of the MLGNR is obtained as follow:

q3ss r − q3s, rq`,, r+j q`,s r+j q`ss r − q`s, rq`,, r+j q`,s rETs F = N>% ETs F

(38)

TE D

where N>% is the eigenvalue of Eq. (38) and can be obtained with the aids of a self-developed MATLAB program. The square roots of these eigenvalues give the natural frequencies of the cantilever MLGNRs.

4. MD simulations

EP

In the MD simulations for free transverse vibration of cantilever MLGNRs, the bonding atomic interaction in MLGNRs is described by the Adaptive Intermolecular Reactive Empirical Bond Order

AC C

(AIREBO) potential [27]. The potential has been widely used to study the mechanical properties of carbon nanotubes, graphene and diamond [27-29]. LAMMPS software is used to carry out the MD simulations. The cantilever MLGNR studied here is armchair-type with the armchair and zigzag direction along the length and width, respectively and has rectangular cross-section with different heights concerning to the number of layers n (n=2,…, 5). The thickness and width of the monolayer GNR are 0.34 nm and 2 nm, respectively. The MLGNR length ranges from 8 to 16 nm in the simulations. The MD simulation procedure carried out by Liu et al. [17] is used in this study. After structural relaxation, a displacement of 1 nm is applied to the free end of each cantilever GNR at the beginning of the simulation

Page 8 of 17

ACCEPTED MANUSCRIPT to initialize a transverse vibration. 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 the GNRs is tracked. Finally, the first two resonance frequencies are extracted using the fast Fourier transform (FFT).

RI PT

5. Results and discussion

Before we pass on the results, we should point out that the values of the bending rigidity, nanoribbon thickness and shear modulus are the important parameters which affect the natural frequencies of

SC

MLGNR. The bending rigidity and interlayer shear modulus of GNR are chiefly sensitive to chirality, size and shape of GNR, temperature, stacking pattern, sp3 bindings, deflection, and ripples [30-36]. But the

M AN U

values reported in the literature are considerably scattered (see Tables 1 and 2).

Table 1: Interlayer shear modulus reported in the literature. Stacking Chirality Method pattern Armchair MD simulation (Dreiding force field) MD simulation Bilayer graphene [21] AB Armchair (AIREBO with epsilon_CC=45.44 meV) MD simulation AB Armchair (AIREBO with epsilon_CC=2.84 meV) MD simulation AA Armchair (AIREBO with epsilon_CC=45.44 meV) Graphite [37] AB Experiment Graphite [38] Random Experiment Few-layer graphene [33] AB DFT (LCAO-S2+vdW) Random DFT (LCAO-S2+vdW) AA DFT (LCAO-S2+vdW) Molecular Mechanics Bilayer graphene [1] AB Armchair (Morse and Lenard-Jones potentials) Molecular Mechanics AB Zigzag (Morse and Lenard-Jones potentials) Bilayer graphene [32] Armchair MD simulation (AIREBO)

EP

TE D

Type of graphitic material Bilayer graphene [17]

G (GPa) 0.25 4.6 0.288 -3.8 5±3 0.18-0.35 4.8 0.19-0.34 -3.8 0.482 0.609 3.01

AC C

Since the accuracy of the analytical results depends on the selected values for bending rigidity, interlayer shear modulus and the thickness of nanoribbon, we should browse through appropriate methods to select these values. The present paper deals with a cantilever MLGNRs. According to the literature, the nonlocal parameter shows a decreasing effect on all natural frequencies of nanobeams with clamped and/or simplysupported conditions at all ends. For the stiffer boundaries, this effect is more predominant leading to the higher reduction on all resonant frequencies values. However, for nanobeams with at least one free end, it is observed that the nonlocal parameter causes negligible increasing effect on the fundamental resonant frequency while decreasing higher-order natural frequencies. This means that the fundamental resonant

Page 9 of 17

ACCEPTED MANUSCRIPT frequency of a cantilever GNR is not basically dependent on the nonlocal parameter. Therefore, the bending rigidity, nanobeam thickness and interlayer shear modulus should be chosen such a way that the difference between the fundamental frequency of MLGNRs obtained by the MD simulations and classical Reddy model (i.e. the nonlocal parameter µ is set to zero) is to be the least. To achieve this goal, the least-

RI PT

square fitting procedure is used and the values of bending rigidity, nanoribbon thickness and interlayer shear modulus are obtained as D=1.2 eV [39]; h=0.34 nm [40]; G=3.01 GPa [32], respectively. Table 2. Thickness and bending rigidity of Graphene

0.3400 0.3400 0.0620 0.0570 0.1317 0.0840 0.6900

SC

0.3350 0.3400 0.3400 0.3400

bending rigidity (eV) Method Atomistic-continuum modeling [1] MD simulation (DRIEDING force field) [17] DFT [28] Atomistic-continuum modeling (Morse force field) [34] Molecular mechanics [31] DFT [46] Atomistic-continuum modeling [48] Atomistic-continuum modeling (AMBER force field) [34] Spring based structural mechanics approach [51] DFT [52] DFT [54] DFT [56] Experiment [39] Analytical method based on secondgeneration Brenner potential without the effect of dihedral angles [59] AIREBO force field [62]

M AN U

Value 0.3400 0.3350 0.3400 0.3350

TE D

Thickness (nm) Method Structural mechanics [41] Molecular mechanics [42] Equivalent-continuum modeling [43] Ultrasonic, sonic resonance, and static test methods [38] Lattice-dynamical model [44] Continuum mechanics approach [45] Molecular Dynamics [47] Tight binding molecular dynamics [49] Structural mechanics: stiffness matrix method [50] Molecular dynamics [40] Empirical potential [53] Atomistic simulation [55] Equivalent continuum model [57] Structural mechanical models [58] Equivalent-continuum modeling [60] Structural mechanics: FE method [61] Molecular dynamics [63]

0.1470

Ab inito computations [65] Local density approximation model [66] Ring theory continuum mechanics [67] Continuum shell modeling [69]

0.0890 0.0750

AC C

EP

0.0660

Continuum model for longwavelength phonons [70] Ab initio calculation [71]

0.0617

0.0750 0.0870 0.0665

Value 2.18 2.14 1.44 1.13

0.36-2.65 1.48 1.62 1.21 2.18 1.61 1.40-1.46 1.52 1.20 0.69

1.12

Atomistic-continuum modeling (COMPASS force field) [64] AIREBO force field [36] AIREBO force field [66]

1.78

Continuum elasticity and the tight-binding atomistic simulation [68] Analytical method based on secondgeneration Brenner potential [59] AIREBO force field [65]

1.40

Analytical method based on first-generation Brenner potential [59]

0.84

0.83 1.17

1.40 1.02

Fig. 2 shows the first two fundamental frequencies of cantilever armchair MLGNRs obtained by the classical and nonlocal Reddy model as well as the MD simulations. According to this figure, by selecting the appropriate values for bending rigidity, nanoribbon thickness and interlayer shear modulus, the first classical frequencies (the nonlocal Reddy model with µ=0) are close to those obtained by the MD simulations. Therefore, it is deduced that the classical Reddy model can accurately predict the

Page 10 of 17

ACCEPTED MANUSCRIPT fundamental resonant frequencies of the cantilever MLGNRs with any number of layers if the bending rigidity, nanoribbon thickness and the interlayer shear modulus are chosen as 1.2 eV, 0.34 nm and 3.01 GPa, respectively. However, there is a considerable difference between the second frequencies obtained

by the MD simulations and those of the classical Reddy model. Furthermore, the difference is not equal

RI PT

for different number of layers of MLGNRs. Therefore, a calibrated nonlocal parameter should be found for the second frequencies of nonlocal Reddy model for each of the number of layers. This procedure can be accomplished by using the least square method and as depicted in Fig. 2 the nonlocal Reddy results have matched very well with those of the MD simulations. The calibrated value of the nonlocal parameter

SC

for each of the number of layers is also valid for all values of the MLGNR length. This observation is like the one reported by Ansari et al. [72] and Arash and Wang [13] for the nonlocal parameter calibration of

M AN U

single- and double-layered graphene sheets with simply supported and clamped edges. For better investigation of the accuracy of the obtained results compared to the MD simulations, an average error has been reported in Table 3. This table shows the average error between the classic and nonlocal second frequencies and those obtained by the MD simulations for different lengths of the MLGNR. It is seen from Table 3 that there is a considerable error between the classic results and those of

TE D

MD simulations whereas the errors significantly decrease by using the nonlocal theory. Table 3: Average error between the second frequencies of the classic and nonlocal Reddy models and the MD simulations for different lengths of the cantilever MLGNR. Calibrated nonlocal parameter (nm2) 6.0 8.8 11.0 12.3

Average Error (%) Classic Nonlocal 30.20 4.56 41.89 7.34 52.89 8.08 59.80 8.63

AC C

EP

Number of layers 2 3 4 5

Another study in Fig. 3 shows the variations of the first two frequencies of the cantilever MLGNR versus the number of layers for L=12 nm. According to this figure, by increasing the number of MLGNR layers the natural frequency increases, as expected. Furthermore, by increasing the number of layers, the difference between the second frequency obtained by the MD simulations and those obtained by the classical Reddy model increases. To reduce the difference, a higher value of the nonlocal parameter is required to match the nonlocal Reddy results with the MD results as reported in Table 3. Since finding the calibrated nonlocal parameter is laborious for different number of layers, we have suggested a curve in Figure 4 which can predict the variations of the nonlocal parameter versus the

Page 11 of 17

ACCEPTED MANUSCRIPT number of layers. As depicted in this figure and mentioned earlier, by increasing the number of layers the calibrated nonlocal parameter increases to predict the best results comparing to the MD simulations. The last investigation performed in this study deals with obtaining the shear correction factor of the nonlocal Timoshenko theory for MLGNRs as done for macro-scale beams [73]. To this end, the values of

RI PT

nonlocal parameter, geometrical properties and the related natural frequencies obtained by the nonlocal Reddy theory is used as the known variables for the nonlocal Timoshenko theory while the shear correction factor is the unknown one of the problem. After solving the problem, the values of shear correction factors for the first and second frequencies of a 12 nm length MLGNR with different number

SC

of layers are listed in Table 4. According to this table the values of shear correction factor is not dependent on the thickness of MLGNRs for the first frequency as the behavior has been reported in [73]

M AN U

while it is the other way round for the shear correction factor of the second frequency. Furthermore, the results tabulated in Table 4 are close to the values of shear correction factor for micro-beams reported in the literature [74-76]. This indicates that one can use the results of shear correction factor of micro-beams with good precision for nano-beams. This study is a preliminary investigation for obtaining the shear correction factor of nonlocal Timoshenko theory in nanoscale field and can be considered as a prelude to

TE D

present other methods and broader values of shear correction factor for the future studies. Table 4: Shear correction factor of nonlocal Timoshenko theory for the first two natural frequencies of MLGNRs

Number of layers

First Frequency

Second Frequency

0.833 0.833 0.833 0.833

0.879 0.876 0.873 0.874

AC C

EP

2 3 4 5

Shear correction factor

Conclusion

The interlayer shear effect of vdWs interactions on nonlocal free vibration of MLGNRs is investigated based on Reddy beam theory. To this end, the natural frequencies of MLGNRs with various layers number and lengths are calculated by using the nonlocal Reddy theory and MD simulations. The results show that the nonlocal Reddy beam model will properly predict the free vibration of MLGNRs if three considerations are taken into account. First, the values of bending rigidity, shear modulus, and thickness of every GNR layer should be taken as 1.2 eV, 3.01 GPa, and 0.34 nm, respectively. Second, the calibrated nonlocal parameter (D) for the first frequencies should be zero, and third, the calibrated

Page 12 of 17

ACCEPTED MANUSCRIPT nonlocal parameter for the second frequencies should be chosen from D ≈ 6.544 Pxyq1.279 @r, where n is the number of layers of the MLGNR. Finally, the shear correction factor of the nonlocal Timoshenko theory is obtained for the values of nonlocal parameter, geometrical properties and the related natural frequencies of MLGNRs calculated by the nonlocal Reddy theory. The present study features the low

RI PT

computational costs, low time-consuming process and the simplicity of modeling the interlayer shear effect via the simple nonlocal Reddy beam theory compared with other theoretical and numerical methods like sandwich beam model with large computations.

SC

Appendix A

Effect of nonlocal parameter on natural frequencies of cantilever nanobeams

M AN U

The nonlocal governing equation of motion of Euler-Bernoulli nanobeams is given as [77]

∂ 2 M nl ∂ 2w ρA + =0 ∂x 2 ∂t 2 where

∂ 2w ∂ 2w − μρA ∂x 2 ∂t 2

By assuming

(A.2)

w ( x,t ) = X ( x ) eiωt and by substituting this definition in Eq. (A.1), the solution of

TE D

M nl = EI

(A.1)

governing equation is obtained as

X( x ) = A1 Sin( λ1 x ) + A2Cos( λ1 x ) + A3Sinh( λ2 x ) + A4Cosh( λ2 x );

EP

( μρAω )

2

2

λ2 =

+ 4EIρAω2

2EI

( μρAω )

AC C

λ1 =

μρAω2 +

− μρAω2 +

2

2EI

2

;

+ 4EIρAω2

(A.3)

;

Boundary conditions for a cantilever nanobeam are given as Clamped end:

X (0) = 0;

dX ( 0) dx

=0

(A.4)

Free end:

Page 13 of 17

ACCEPTED MANUSCRIPT  nl  d2 X 2  = 0 ⇒ M L ( )   EI 2 + μρAω X  = 0   dx x =L  nl  d3 X  ∂M ( L ) 2 dX  = 0 ⇒  EI 3 + μρAω  =0  ∂x dx  x = L  dx 

(A.5)

Bernoulli nanobeams

)

(

) (

)

( EIλ

2 1

)(

)

− μρAω2 EIλ22 + μρAω2 × ( 2λ1 λ2Cos ( λ1 L ) Cosh ( λ2 L ) −



2 2

)

)

− λ21 Sin ( λ1 L ) Sinh ( λ2 L ) = 0

SC

(

2 2 λ1 λ2  ( EI ) λ41 + λ42 + 2μρAEIω2 λ22 − λ21 + 2 μρAω2  +  

RI PT

Applying boundary conditions results in the following characteristic equation for cantilever Euler-

(A.6)

Based on the above outlined formulations and by aids of the MATHEMATICA program solver a self-

M AN U

developed computer program is written by which the natural frequencies of the nanobeam can be obtained.

In Table A.1 the first two natural frequencies are given for various values of the nanobeam length and nonlocal parameter. It can be concluded from Table A.1 that the nonlocal parameter shows negligible increasing effects on the fundamental frequencies but significant decreasing effects on the second

TE D

frequencies of cantilever nanobeams.

EP

Table A.1: The first two natural frequencies of cantilever Euler-Bernoulli nanobeams for various nonlocal parameters and nanobeam lengths. Length Nonlocal parameter 1st Frequency 2nd Frequency (nm2) (GHz) (GHz) 0 2 4

160.61 162.85 165.27

1006.54 838.93 724.27

10

0 2 4

102.79 103.70 104.65

644.18 570.37 513.84

20

0 2 4

25.70 25.75 25.81

161.05 155.90 151.14

30

0 2 4

11.42 11.43 11.44

71.58 70.54 69.54

AC C

8

Page 14 of 17

ACCEPTED MANUSCRIPT References

AC C

EP

TE D

M AN U

SC

RI PT

[1] Kordkheili SH, Moshrefzadeh-Sani H. Mechanical properties of double-layered graphene sheets. Computational Materials Science. 2013;69:335-43. [2] Mattevi C, Eda G, Agnoli S, Miller S, Mkhoyan KA, Celik O, et al. Evolution of electrical, chemical, and structural properties of transparent and conducting chemically derived graphene thin films. Advanced Functional Materials. 2009;19(16):2577-83. [3] Balandin AA, Ghosh S, Bao W, Calizo I, Teweldebrhan D, Miao F, et al. Superior thermal conductivity of single-layer graphene. Nano letters. 2008;8(3):902-7. [4] Jiao L, Zhang L, Wang X, Diankov G, Dai H. Narrow graphene nanoribbons from carbon nanotubes. Nature. 2009;458(7240):877-80. [5] Alibeigloo A. Three-dimensional free vibration analysis of multi-layered graphene sheets embedded in elastic matrix. Journal of Vibration and Control. 2013;19(16):2357-71. [6] Wang J, Tian M, He X. Analysis of free vibration of embedded multi-layered graphene sheets. Procedia Engineering. 2012;31:641-6. [7] Lin R. Nanoscale vibration characteristics of multi-layered graphene sheets. Mechanical Systems and Signal Processing. 2012;29:251-61. [8] Lin R. Nanoscale vibration characterization of multi-layered graphene sheets embedded in an elastic medium. Computational Materials Science. 2012;53(1):44-52. [9] He X, Wang J, Liu B, Liew KM. Analysis of nonlinear forced vibration of multi-layered graphene sheets. Computational Materials Science. 2012;61:194-9. [10] Arghavan S, Singh A. Effects of van der Waals interactions on the nonlinear vibration of multilayered graphene sheets. Journal of Physics D: Applied Physics. 2012;45(45):455305. [11] Wang J, He X, Kitipornchai S, Zhang H. Geometrical nonlinear free vibration of multi-layered graphene sheets. Journal of Physics D: Applied Physics. 2011;44(13):135401. [12] Chandra Y, Chowdhury R, Scarpa F, Adhikaricor S. Vibrational characteristics of bilayer graphene sheets. Thin solid films. 2011;519(18):6026-32. [13] Arash B, Wang Q. Vibration of single-and double-layered graphene sheets. Journal of Nanotechnology in Engineering and Medicine. 2011;2(1):011012. [14] Ansari R, Arash B, Rouhi H. Nanoscale vibration analysis of embedded multi-layered graphene sheets under various boundary conditions. Computational Materials Science. 2011; 50(11):3091100. [15] Ansari R, Arash B, Rouhi H. Vibration characteristics of embedded multi-layered graphene sheets with different boundary conditions via nonlocal elasticity. Composite Structures. 2011;93(9):2419-29. [16] 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. [17] Liu Y, Xu Z, Zheng Q. The interlayer shear effect on graphene multilayer resonators. Journal of the Mechanics and Physics of Solids. 2011;59(8):1613-22. [18] 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. [19] 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. [20] Rokni H, Lu W. A continuum model for the static pull-in behavior of graphene nanoribbon electrostatic actuators with interlayer shear and surface energy effects. Journal of Applied Physics. 2013;113(15):153512. [21] Shen Y, Wu H. Interlayer shear effect on multilayer graphene subjected to bending. Applied Physics Letters. 2012;100(10):101909. [22] Leung A. An improved third order beam theory. Journal of sound and vibration. 1990;142(3):527-8. [23] Eringen AC, Edelen D. On nonlocal elasticity. International Journal of Engineering Science. 1972;10(3):233-48. [24] Hosseini-Hashemi S, Fakher M, Nazemnezhad R. Surface effects on free vibration analysis of nanobeams using nonlocal elasticity: a comparison between Euler–Bernoulli and Timoshenko. J Solid Mech. 2013;5(3):290-304.

Page 15 of 17

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[25] Wang C, Zhang Y, He X. Vibration of nonlocal Timoshenko beams. Nanotechnology. 2007;18(10):105401. [26] 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. [27] Stuart SJ, Tutein AB, Harrison JA. A reactive potential for hydrocarbons with intermolecular interactions. The Journal of chemical physics. 2000;112(14):6472-86. [28] Grantab R, Shenoy VB, Ruoff RS. Anomalous strength characteristics of tilt grain boundaries in graphene. Science. 2010;330(6006):946-8. [29] Peng B, Locascio M, Zapol P, Li S, Mielke SL, Schatz GC, et al. Measurements of near-ultimate strength for multiwalled carbon nanotubes and irradiation-induced crosslinking improvements. Nature nanotechnology. 2008;3(10):626-31. [30] Kang JW, Lee S. Molecular dynamics study on the bending rigidity of graphene nanoribbons. Computational Materials Science. 2013;74:107-13. [31] Giannopoulos GI. Elastic buckling and flexural rigidity of graphene nanoribbons by using a unique translational spring element per interatomic interaction. Computational Materials Science. 2012;53(1):388-95. [32] Zhang Y, Wang C, Cheng Y, Xiang Y. Mechanical properties of bilayer graphene sheets coupled by sp3 bonding. Carbon. 2011;49(13):4511-7. [33] Savini G, Dappe Y, Öberg S, Charlier J-C, Katsnelson M, Fasolino A. Bending modes, elastic constants and mechanical stability of graphitic systems. Carbon. 2011;49(1):62-9. [34] Wang Q. Simulations of the bending rigidity of graphene. Physics Letters A. 2010;374(9):11803. [35] Zhao H, Min K, Aluru N. Size and chirality dependent elastic properties of graphene nanoribbons under uniaxial tension. Nano letters. 2009;9(8):3012-5. [36] Liu P, Zhang Y. Temperature-dependent bending rigidity of graphene. Applied Physics Letters. 2009;94(23):231912. [37] Bosak A, Krisch M, Mohr M, Maultzsch J, Thomsen C. Elasticity of single-crystalline graphite: inelastic X-ray scattering study. Phys Rev B. 2007;75(15):153408. [38] Blakslee O, Proctor D, Seldin E, Spence G, Weng T. Elastic constants of compression-annealed pyrolytic graphite. Journal of Applied Physics. 1970;41(8):3373-82. [39] Nicklow R, Wakabayashi N, Smith H. Lattice dynamics of pyrolytic graphite. Physical Review B. 1972;5(12):4951. [40] Jin Y, Yuan F. Simulation of elastic properties of single-walled carbon nanotubes. Composites Science and Technology. 2003;63(11):1507-15. [41] Chang T, Gao H. Size-dependent elastic properties of a single-walled carbon nanotube via a molecular mechanics model. Journal of the Mechanics and Physics of Solids. 2003;51(6):1059-74. [42] Cho J, Luo J, Daniel I. Mechanical characterization of graphite/epoxy nanocomposites by multiscale analysis. Composites science and technology. 2007;67(11):2399-407. [43] Sakhaee-Pour A. Elastic properties of single-layered graphene sheet. Solid State Communications. 2009;149(1):91-5. [44] Al-Jishi R, Dresselhaus G. Lattice-dynamical model for graphite. Physical Review B. 1972;5(12):4951. [45] Reddy C, Rajendran S, Liew K. Equilibrium configuration and continuum elastic properties of finite sized graphene. Nanotechnology. 2006;17(3):864. [46] Roberts M, Clemons C, Wilber J, Young G, Buldum A, Quinn D. Continuum Plate Theory and Atomistic Modeling to Find the Flexural Rigidity of a Graphene Sheet Interacting with a Substrate. J Nanotechnol. 2010;2010:1-8. [47] Lu JP. Elastic properties of carbon nanotubes and nanoropes. Physical Review Letters. 1997;79(7):1297. [48] Scarpa F, Adhikari S, Gil A, Remillat C. The bending of single layer graphene sheets: the lattice versus continuum approach. Nanotechnology. 2010;21(12):125702. [49] Hernandez E, Goze C, Bernier P, Rubio A. Elastic properties of C and B x C y N z composite nanotubes. Physical Review Letters. 1998;80(20):4502. [50] Li C, Chou T-W. A structural mechanics approach for the analysis of carbon nanotubes. International Journal of Solids and Structures. 2003;40(10):2487-99.

Page 16 of 17

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[51] Wei Y, Wang B, Wu J, Yang R, Dunn ML. Bending rigidity and Gaussian bending stiffness of single-layered graphene. Nano lett. 2012;13(1):26-30. [52] Muñoz E, Singh AK, Ribas MA, Penev ES, Yakobson BI. The ultimate diamond slab: GraphAne versus graphEne. Diam Relat Mater. 2010;19(5):368-73. [53] Brenner DW. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. Physical Review B. 1990;42(15):9458. [54] Arroyo M, Belytschko T. Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule. Phys Rev B. 2004;69(11):115415. [55] Huang Y, Wu J, Hwang K. Thickness of graphene and single-wall carbon nanotubes. Physical review B. 2006;74(24):245413. [56] Sánchez-Portal D, Artacho E, Soler JM, Rubio A, Ordejón P. Ab initio structural, elastic, and vibrational properties of carbon nanotubes. Phys Rev B. 1999;59(19):12678. [57] Hemmasizadeh A, Mahzoon M, Hadi E, Khandan R. A method for developing the equivalent continuum model of a single layer graphene sheet. Thin Solid Films. 2008; 516(21):7636-40. [58] Scarpa F, Adhikari S, Phani AS. Effective elastic mechanical properties of single layer graphene sheets. Nanotechnology. 2009;20(6):065709. [59] Cadelano E, Giordano S, Colombo L. Interplay between bending and stretching in carbon nanoribbons. Phys Rev B. 2010;81(14):144105. [60] Odegard GM, Gates TS, Nicholson LM, Wise KE. Equivalent-continuum modeling of nanostructured materials. Composites Science and Technology. 2002;62(14):1869-80. [61] Tserpes K, Papanikos P. Finite element modeling of single-walled carbon nanotubes. Composites Part B: Engineering. 2005;36(5):468-77. [62] Lu Q, Arroyo M, Huang R. Elastic bending modulus of monolayer graphene. J Phys D Appl Phys. 2009;42(10):102002. [63] Yakobson BI, Brabec C, Bernholc J. Nanomechanics of carbon tubes: instabilities beyond linear response. Physical review letters. 1996;76(14):2511. [64] Wang Q, Liew K. Molecular mechanics modeling for properties of carbon nanotubes. J Appl Phys. 2008;103(4):046103. [65] Kudin KN, Scuseria GE, Yakobson BI. C2F, BN, and C nanoshell elasticity from ab initio computations. Physical Review B. 2001;64(23):235406. [66] Tu Z-c, Ou-Yang Z-c. Single-walled and multiwalled carbon nanotubes viewed as elastic tubes with the effective Young’s moduli dependent on layer number. Physical Review B. 2002;65(23):233407. [67] Vodenitcharova T, Zhang L. Effective wall thickness of a single-walled carbon nanotube. Physical Review B. 2003;68(16):165401. [68] Koskinen P, Kit OO. Approximate modeling of spherical membranes. Phys Rev B. 2010;82(23):235420. [69] Pantano A, Parks DM, Boyce MC. Mechanics of deformation of single-and multi-wall carbon nanotubes. Journal of the Mechanics and Physics of Solids. 2004;52(4):789-821. [70] Goupalov S. Continuum model for long-wavelength phonons in two-dimensional graphite and carbon nanotubes. Physical Review B. 2005;71(8):085420. [71] Wang L, Zheng Q, Liu JZ, Jiang Q. Size dependence of the thin-shell model for carbon nanotubes. Physical review letters. 2005;95(10):105501. [72] Ansari R, Sahmani S, Arash B. Nonlocal plate model for free vibrations of single-layered graphene sheets. Physics Letters A. 2010;375(1):53-62. [73] Heyliger P, Reddy J. A higher order beam finite element for bending and vibration problems. Journal of sound and vibration. 1988;126(2):309-26. [74] Hutchinson J. Shear coefficients for Timoshenko beam theory. Journal of Applied Mechanics. 2001;68(1):87-92. [75] Jensen JJ. On the shear coefficient in Timoshenko's beam theory. Journal of Sound and Vibration. 1983;87(4):621-35. [76] Cowper G. The shear coefficient in Timoshenko’s beam theory. Journal of applied mechanics. 1966;33(2):335-40. [77] Hosseini-Hashemi S, Fakher M, Nazemnezhad R. Surface effects on free vibration analysis of nanobeams using nonlocal elasticity: a comparison between Euler-Bernoulli and Timoshenko. J Solid Mech. 2013;5(3):290-304.

Page 17 of 17

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Fig. 1: Schematic of a multilayer graphene nanoribbon

b) Trilayer

EP

c) Fourlayer

TE D

M AN U

a) Bilayer

SC

RI PT

ACCEPTED MANUSCRIPT

d) Fivelayer

AC C

Figure 2: Comparison of the first two natural frequencies between MD simulations and the classical and nonlocal Reddy model for MLGNR.

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

Figure 3: Variations of the frequencies of MLGNRs versus the number of layers for L=12 nm.

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

Figure 4: Variations of calibrated nonlocal parameter versus number of layers

ACCEPTED MANUSCRIPT Multilayer graphene nanoribbons are modeled as a single Reddy nanobeam.



Interlayer shear modulus of MLGNRs is considered as Reddy’s shear modulus.



Appropriate geometrical and mechanical properties of MLGNRs are selected.



Nonlocal parameter is calibrated for the first and second frequencies of MLGNRs.



Shear correction factor of the nonlocal Timoshenko model is obtained for MLGNRs.

AC C

EP

TE D

M AN U

SC

RI PT