Finite element buckling analysis of multi-layered graphene sheets on elastic substrate based on nonlocal elasticity theory

Accepted Manuscript Finite Element Buckling Analysis of Multi-Layered Graphene Sheets on Elastic Substrate Based on Nonlocal Elasticity Theory A. Anjomshoa, A.R. Shahidi, B. Hassani, E. Jomehzadeh PII: DOI: Reference:

S0307-904X(14)00139-5 http://dx.doi.org/10.1016/j.apm.2014.03.036 APM 9924

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

19 March 2013 3 December 2013 16 March 2014

Please cite this article as: A. Anjomshoa, A.R. Shahidi, B. Hassani, E. Jomehzadeh, Finite Element Buckling Analysis of Multi-Layered Graphene Sheets on Elastic Substrate Based on Nonlocal Elasticity Theory, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm.2014.03.036

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.

Finite Element Buckling Analysis of Multi-Layered Graphene Sheets on Elastic Substrate Based on Nonlocal Elasticity Theory A. Anjomshoaa*, A. R. Shahidib, B. Hassania, E. Jomehzadehc a

Department of Mechanical Engineering, Ferdowsi University of Mashhad, 91775-1111 Mashhad, Iran

b

Department of Mechanical Engineering, Isfahan University of Technology, 84156-83111 Isfahan, Iran

c

Department of Mechanical Engineering, Graduate University of Advanced Technology, 7631133131 Kerman, Iran

Abstract

Graphene-polymer nano-composites are one of the most applicable engineering nanostructures with superior mechanical properties. In the present study, a finite element (FE) approach based on the size dependent nonlocal elasticity theory is developed for buckling analysis of nanoscaled multi-layered graphene sheets (MLGSs) embedded in polymer matrix. The van der Waals (vdW) interactions between the graphene layers and graphene-polymer are simulated as a set of linear springs using the Lennard–Jones potential model. The governing stability equations for nonlocal classical orthotropic plates together with the weighted residual formulation are employed to explicitly obtain stiffness and buckling matrices for a multi-layered super element of MLGS. The accuracy of the current finite element analysis (FEA) is approved through a comparison with molecular dynamics (MD) and analytical solutions available in the literature. Effects of nonlocal parameter, dimensions, vdW interactions, elastic foundation, mode numbers and boundary conditions on critical in-plane loads are investigated for different types of MLGS. It is found that buckling loads of MLGS are generally of two types namely In-Phase (INPH) and Out-of-Phase (OPH) loads. The INPH loads are independent of interlayer vdW interactions while the OPH loads depend on vdW interactions. It is seen that the decreasing effect of nonlocal parameter on the OPH buckling loads dwindles as the interlayer vdW interactions become

*

Corresponding author. E-mail address: [email protected] or [email protected] (A. Anjomshoa).

1

stronger. Also, it is found that the small scale and polymer substrate have noticeable effects on the buckling loads of embedded MLGS. Keywords: Buckling analysis; Multi-layered graphene sheet; Nonlocal elasticity theory; Finite element method; Elastic substrate. 1. Introduction Plate like elements in micro and nano scales form a novel class of structural components which have been integrated into different micro-electro-mechanical systems (MEMS) and nanoelectro-mechanical systems (NEMS). Graphene is a typical carbon plate like nanostructure with lots of applications due to its unique mechanical, chemical and electrical properties. The covalent bond of carbon atoms induces a Young modulus in range of 1 Tera-Pascal for graphene sheet (GS) reported by experiment [1], and atomic or molecular simulations [2,3]. These intrinsic properties lead to functional graphene based nanostructures utilized in different NEMS applications as sensors [4] or reinforcements in composites [5,6]. It is seen that graphene appears in two types of structure, namely single-layered graphene sheet (SLGS) and multi-layered graphene sheet (MLGS). Due to specific properties of nanostructures (graphene based and non-graphene based), a wide range of experimental, numerical and theoretical studies has been performed on their mechanical properties by different research communities. Among these methods, experimental measurements are hard to reproduce and depend on the development of equipments for manipulation of nanosized objects. In the case of computational nano-mechanics, different numerical techniques based on semi-empirical approaches such as quantum chemical ab-initio methods [2,7], density functional theory (DFT) [7], molecular dynamics (MD) [8], multi-scale modeling [9], etc. have been developed for analysis of nanostructures with various molecular or atomic structures. Although these semiempirical numerical methods produce results which are in good agreement with experiment, they are computationally restricted by the number of total atoms and bonds included in the system. Hence, modeling of nanostructures as an equivalent continuum structure is being the focus of

2

interest since it provides a balance between accuracy and efficiency. The continuum models benefit from the computational efficiency of continuum theories and at the same time produce desired results comparable to atomistic models. Due to their computational efficiency, continuum models can be effectively used to simulate very small to very large systems. However, these continuous equivalent models of nanostructures must satisfy all the conditions affecting the physical properties of structures at very small scales such as “nano” scale. In fact, when the dimensions of a system decrease to nanometer, they become comparable to the inter-atomic or inter-molecular spacing of that system, and the system cannot be considered as a continuous medium any more. Also, at small size the influence of inter-atomic and inter-molecular cohesive forces on static and dynamic responses tends to be significant. The collection of these effects is known as the “Quantum”, “Small scale” or “Size” effect [10]. Since the classical continuum models are unable to account for quantum effects, different modified size dependent theories such as strain gradient theory [11], couple stress theory [12], surface elasticity theory [13] and nonlocal elasticity theory [14,15] have been developed. Among these, it is seen that the nonlocal elasticity theory can produce well-matched results with lattice dynamics [15]. Using the nonlocal elasticity theory, static and dynamic analysis of microtubules (MTs), single-walled carbon nanotubes (SWCNTs) and multi-walled carbon nanotubes (MWCNTs) have been investigated by different researchers [16-23]; In the case of graphene sheets and nanoplates, buckling of simply supported orthotropic single-layered graphene sheets was reported by Murmu and Pradhan [24] using Navier’s approach. Aksencer and Aydogdu [25] investigated the vibration and buckling of nanoplates with different boundary conditions by Levy method. Pradhan [26] reported the buckling analysis and small scale effect of biaxially compressed graphene sheets using non-local elasticity theory and Levy’s approach. Jomehzadeh and Saidi [27,28] considered large amplitude vibration of multi-layered graphene sheets with different edge supports and in nonlinear polymer matrix. They used the Hamilton’s principle to obtain the coupled nonlinear governing vibration equations of MLGS according to von-Karman geometrical model. Further, they decoupled the nonlocal governing equations associated with three-dimensional vibration of nano3

plates [29]. Using the well-known Galerkin method, Babaie and Shahidi [30] studied the vibrations of quadrilateral MLGSs. Further, they investigated the small scale effect on the buckling of quadrilateral nanoplates using the same solution procedure [31]. Buckling of embedded simply supported MLGSs was investigated by Pradhan and Phadikar [32] using the nonlocal classical plate model and Navier’s approach. Murmu et al. [33,34] also studied the buckling of double nanoplate systems and bonded double nanoplate systems via nonlocal theory.

In the present study, buckling analysis of orthotropic MLGS rested on elastic substrate is carried out employing the size dependent nonlocal elasticity theory and widely applicable finite element method (FEM). Due to the numerous abilities, the numerical finite element method is broadly employed for static and dynamic analysis of plate structures with different geometries, boundary conditions, variations in thickness, and stiffness [35]. In the case of nonlocal theory, finite element formulation was derived by Phadikar and Pradhan [36] for vibration and buckling of nanobeams and nanoplates. Adali [37] reported the variational principle and natural boundary conditions for vibration and buckling of embedded orthotropic MLGSs via nonlocal classical plate theory. Also, Ansari et al. [38] developed a nonlocal finite element model for vibration analysis of isotropic embedded multi-layered graphene sheets using the Mindlin plate theory. Here, by considering Lennard-Jones interactions between graphene layers and graphene-polymer, governing stability equations for MLGSs under biaxial compression are developed based on the size dependent nonlocal continuum theory. From the weak forms of governing equations, stiffness and buckling matrices for a multilayered element of MLGS are explicitly obtained using the Galerkin method. Then, buckling loads are calculated through solving the global finite element equation. The validation and accuracy of the current FEA are verified by comparing its results with existing analytical and molecular dynamics solutions in literature. To the best of author’s knowledge, finite element buckling analysis of MLGSs rested on elastic substrate and based on nonlocal elasticity theory has not been reported in available literature. It worth noting that the results from current finite element analysis of MLGS

4

can be considered as a useful source since other carbon nanostructures such as single- and multiwalled carbon nanotubes and fullerene are assumed to be special configurations of multi-layered graphene sheets.

2. Governing Equations 2.1 Nonlocal Elasticity Theory The theory of nonlocal elasticity was originally used by Eringen to study surface waves and screw dislocation in solids [14,15]. Unlike the classical elasticity theory, the nonlocal elasticity theory assumes that stress at a reference point in a body depends not only on strains at that point, but also on strains at all other points of the body. The nonlocal elasticity theory is in accordance with the atomic theory of lattice dynamics [15]. In the limit when the effect of strains at points other than reference point is neglected, the nonlocal elasticity theory conforms to the classical (local) theory of elasticity. According to Eringen [14,15], the differential form of nonlocal constitutive equation is written as

(1 − τ

2 2 2 le ∇

∇ 2 (•) =

)σ

( nl ) ij

= σ ij(l ) = Sijkl ε kl

(1a)

∂2 ∂2 ∂2 ( • ) + ( • ) + (•) ∂x 2 ∂y 2 ∂z 2

(1b)

Here, σ ij(nl ) , σ ij(l ) , Sijkl and ε k l are nonlocal stress tensor, local stress tensor, elasticity tensor and strain tensor, respectively. Considering orthotropic behavior for the graphene layers, we will have

(1 − τ

2 2 2 le ∇

)σ

( nl ) xx

=

1 1 −ν x ν y

(E

x

ε xx + ν x E y ε yy

)

(2a)

5

1

(1 − τ

2 2 2 le ∇

)σ

( nl ) yy

=

(1 − τ

2 2 2 le ∇

)σ

( nl ) xy

= 2G xy ε xy

1 −ν x ν y

(ν

y

E x ε xx + E y ε yy

(2b)

)

(2c)

In which, Ex and Ey are Young’s moduli, ν x and ν y are Poisson’s ratios and Gxy is the shear modulus of the graphene layer, respectively. In Eqs. (1a) and (2) τ = e 0 l i / l e is the nonlocal modulus which represents the intensity of small scale effect. It depends on a characteristic length ratio li / le where l i is an internal characteristic length (lattice parameter, size of grain, granular distance or length of atomic bond) and le is an external characteristic length of the system which can be the wavelength or one of the sample dimensions. For example for graphene sheets, the internal characteristic length can be the length of covalent bond between two adjacent carbon atoms i.e. lc-c shown in Fig. 1(a). The parameter e0 is a constant which needs to be determined for each material independently from experiments or by matching dispersion curves of plane waves with those obtained from lattice dynamics [15]. It is obvious that by putting τ = 0 in Eq. (1a), the well-known classical constitutive equation for Hookean solids can be obtained.

2.2 Nonlocal Plate Model for MLGS Schematics of an n-layered graphene sheet under biaxial compression together with the graphengraphene and graphene-polymer interactions are shown in Fig. 1. The Cartesian coordinate system is chosen for deriving the equilibrium equations. The displacement fields for jth GS at time t according to classical plate theory (CLPT) are

u( j ) = u0( j ) ( x, y, t) − zw,(xj ) , v( j ) = v0( j ) ( x, y, t ) − zw,(yj ) , w( j ) = w( j ) ( x, y, t ),

j = 1,2,..., n.

(3)

Here u0(j)(x,y,t), v0(j)(x,y,t) and w(j)(x,y,t) denote displacement of the point (x,y) placed in the midsurface of jth graphene layer along x, y and z directions, respectively. Also, j=1 and j=n denote the

6

bottommost layer and topmost layer, respectively. The corresponding strain components based on displacement fields in Eq. (3) are obtained as

ε xx( j ) =

∂u0( j ) ∂v0( j ) ∂ 2 w( j ) ∂ 2 w( j ) 1  ∂u0( j ) ∂v0( j ) ∂ 2 w( j ) ( j) ( j) , , 2 −z = − = + − ε z ε z  yy xy ∂x ∂y 2  ∂y ∂x ∂x∂y ∂x 2 ∂y 2

ε zz( j ) = ε xz( j ) = ε (yzj ) = 0

  , 

(4)

Considering the small scale effect through nonlocal theory, the nonlocal stress and nonlocal moment resultants are defined as

h 2

( nl ) N xx =

∫σ

( nl ) xx dz ,

( nl ) ( nl ) i.e. N xx − µ∇ 2 N xx =

h − 2

h 2

( nl ) N yy =

∫σ −

( nl ) yy dz ,

h 2

h 2

( nl ) N xy =

∫σ −

( nl ) xy dz ,

( nl ) i.e. N yy − µ∇ 2 N (yynl ) =

h 1 − ν xν y

(E u

( j) x 0, x

+ ν x E y v0,( jy)

(5a)

h Exν y u0,( jx) + E y v0,( jy) 1 − ν xν y

(

(

)

( nl ) ( nl ) i.e. N xy − µ∇ 2 N xy = Gxy h u0,( jy) + v0,( jx)

) (5b)

)

h 2

(5c)

h 2

( nl ) M xx =

∫ zσ −

( nl ) xx dz ,

( nl ) ( nl ) i.e. M xx − µ∇ 2 M xx = − D11 w,(xxj ) − D12 w,(yyj )

h 2

(5d)

h 2

( nl ) M yy =

∫ zσ −

( nl ) yy dz ,

( nl ) ( nl ) i.e. M yy − µ∇ 2 M yy = − D12 w,(xxj ) − D22 w,(yyj )

h 2

(5e)

7

h 2

∫ zσ

( nl ) M xy =

−

( nl ) xy dz ,

( nl ) ( nl ) i.e. M xy − µ∇ 2 M xy = −2 D66 w,(xyj )

h 2

(5f)

Where, h represents the thickness of graphene layer and µ= τ 2 le2 = (e0li)2 is the nonlocal parameter. Also, Dij are the bending rigidities of graphene layer defined as

D11 =

Ex h3

(

12 1 − ν xν y

)

, D12 =

ν y E x h3

(

12 1 − ν xν y

)

=

ν x Ey h3

(

12 1 −ν xν y

)

, D22 =

E y h3

(

12 1 − ν xν y

)

, D66 =

Gxy h3 12

(6)

Using the Hamilton’s principle, the governing equations of motion for each graphene layer can be obtained as

(N ) + (N ) ( nl ) xx

( nl ) xy

,x

− ρ hu ,(ttj ) = 0

,y

(N ) + (N )

,y

(M )

)

( nl ) + M yy

( j) y w, y

) + (N

( nl ) xy

( nl ) xx

( nl ) yy

,x

, xx

(

+ N xx w,(xj )

(

( nl ) + 2 M xy

) + (N ,x

(7a)

− ρ hv ,(ttj ) = 0

, xy

(

,y

(7b)

)

, yy

( j) xy w, y

) + (N ,x

( j) yx w, x

)

,y

+ qe( j ) − ρ hw,(ttj ) = 0

(7c)

Here, ρ denotes the density of each GS which is assumed to be the same for all layers and qe(j) represents the external load induced from the vdW interactions and/or elastic foundation on jth layer defined as below

n

q (e j ) =

∑ c ∆w ij

ij

− k (f j ) w( j ) ,

j = 1, 2,..., n.

(8)

i =1

8

Here, the first term in the right hand side of the equation represents the vdW pressure exerted on the layer j by other graphene layers in which Δwij= w(i)-w(j) and cij is the graphene-graphene vdW coefficient. The second term denotes the vdW interaction between the carbon atoms of graphene sheet and molecules of polymer substrate in which kf(j) is the graphene-polymer vdW coefficient. The graphene-graphene and graphene-polymer vdW interactions coefficients can be obtained using the interaction potential energy, as a function of the inter-atomic spacing using the Lennard–Jones model [39,40] as

U LJ

 σ 12  σ 6  = 4ε   −     d    d  (9)

where, ULJ denotes the Lennard-Jones potential, σ and ε are parameters changing for different types of material and d is distance between interacting atoms that have the Lennard-Jones interaction. The van der Waals forces can be determined by taking the derivative of the potential energy. Differentiating the Eq. (9) with respect to distance d and expressing the Taylor expansion of the results around the equilibrium position, the vdW graphene-graphene and graphene-polymer coefficients can be obtained as [28]

 24ε cij = cgraphene (i )− graphene ( j ) = ρC2  2C −C  σ C −C

   13πσ C14−C 7πσ C8 −C −   12   3
j −
i 3
j −
i 

(

)

(

)

 , 6 

i ≠ j.

(10)

and

 8εC −CH 2 k ( j ) = k graphene ( j )− polymer = ρC ρCH 2  2  σ C −CH  2

  13πσ C14−CH 7πσ C8 −CH 2 2  − 5   11
11 5
j   j

( )

9

( )

   

(11)

(

Here, ρC = 4 / 3 3l cc2

) is the area density of the carbon atoms in which l

c-c

is the length of carbon–

carbon bond shown in Fig. 1a. Also, ρCH 2 is the volume density of polymer molecules (number of polymer molecules per unit volume) which is equal to 31.00nm-3 for the –CH2– unit in polyethylene [28,40]. In Eq. (10) and (11), ħj is the distance between jth graphene layer and elastic foundation in thickness direction when all the graphene layers and polymer substrate are in their equilibrium positions. For graphene-polyethylene in equilibrium we have:

εC-C = 0.002390eV, σC-C = 0.3415nm,

εC-CH2= 0.004656eV and σC-CH2 = 0.3825nm [28,40]. The equilibrium distances for embedded multilayered graphene sheet are reported by Lu et al. [40] from minimizing the total cohesive energy between the layers. For example, for a specific model of double layered graphene sheet (DLGS) in a polyethylene medium the equilibrium distances i.e. ħ2-ħ1 and ħ1 are found to be 0.997σC-C and 0.851σC-C, respectively [28,40]. Putting the appropriate values in Eq. (10) and (11), the graphenegraphene and graphene-polymer vdW coefficients will be obtained as c12= 84.16GPa.nm-1, k(1)= 28.49GPa.nm-1 and k(2)= 0.77GPa.nm-1.

Assuming identical material properties and geometry for all graphene layers and using Eq. (5) and (7) together with adjacent equilibrium criterion [41], the decoupled form of stability equation for each graphene layer will be obtained as below

 n j) j) j) L(w( j ) ) = D11w,(xxxx + 2( D12 + 2 D66 )w,(xyxy + D22 w,(yyyy − (1 − µ∇2 )  cij ∆wij − k (f j ) w( j )  i =1 + N xx w,(xj ) + N yy w,(yj ) + Nxy w,(yj ) + N yx w,(xj ) − ρ hw,(ttj )  = 0, j = 1,2,..., n.  ,x ,y ,x ,y

∑

(

) (

) (

) (

)

(12)

In the above governing stability equation the small scale effect is taken into account through the nonlocal parameter µ. It is obvious that by putting this parameter equal to zero, the governing stability equations for embedded MLGSs according to the classical continuum theory will be achieved. 10

3. Finite Element Solution A four nodes rectangular bicubic Hermitian plate element [35] with total sixteen degrees of freedom shown in Fig. 2 is adopted for the current finite element analysis. The associated nodal displacement vector for each element in jth layer is defined as

{ } { dˆ ( j )

e

= [ χ ]1

T

[ χ ]2 [ χ ]3 [ χ ]4 }



[ χ ]i = Wi( j )

,



∂Wi ( j ) ∂x

∂Wi ( j ) ∂y

∂ 2Wi ( j )   , i = 1, 2,3, 4. ∂x∂y  (13)

Here, i represents the node number. Using this nodal displacement vector, the corresponding transverse displacement of mid-surface of the jth GS inside the element domain is defined as 16

W ( j ) ( x, y) =

∑dˆ

( j) ( j) i Ni ( x, y)

(14)

i =1

In which, Ne(j) is the shape function vector defined as

Ne( j ) ( x, y) = [ Φ][ J ]−1

(15)

where, [Ф] and [J] are, respectively, polynomial vector and the coordinate matrix defined as

[Φ] =[1, x , y , x 2 , xy , y 2 , x 3, x 2 y , xy 2 , y 3 , x 3 y , x 2 y 2 , xy 3, x 3 y 2 , x 2 y 3 , x 3 y 3 ]

(16)

,

T

[ J ] = {λ}1 {λ}2 {λ}3 {λ}4 

{

T

, {λ}i = [ Φ]

[Φ]T,x [Φ]T, y [Φ]T,xy }

, i = 1,2,3,4. ( x =x i , y = y i )

(17)

Assuming pure biaxial compression i.e. Nxy = Nyx = 0 in Eq. (12) and employing the Galerkin method, yields 11

∫ L (W

( j)

)

( x, y ) δ W ( j ) dAe = 0,

j = 1,2,..., n.

(18)

Ae

here, δW(j) denotes the weight function. After transferring the integration in Eq. (18) to region Ae and boundary Γe of the element using the divergence theorem, following weak form will be obtained for each element of the jth graphene layer

∫Π

( j)

∫

( x, y)dA + Λ( j ) (s )ds = 0,

j = 1, 2,..., n.

(19)

Γe

Ae

in which,

( j) Π( j) (x, y) = D11W,xx δW,xx( j) + D12 W,(xxj)δW,(yyj) +W,(yyj)δW,(xxj)  + D22W,(yyj)δW,(yyj) + 4D66W,(xyj)δW,(xyj) n

∑

− cij (∆Wij )δW( j) + µ ∆Wij,xδW,x( j) +∆Wij,yδW,(yj)  + k(f j) W( j)δW( j) + µ W,(xj)δW,(xj) +W,(yj)δW,(yj)    i=1

(

)

)

(

( j) +Nxx W,(xj)δW,(xj) + µ W,xx δW,(xxj) +W,xy( j)δW,(xyj)  + Nyy W,(yj)δW,(yj) +µ W,(xyj)δW,xy( j) +W,(yyj)δW,(yyj)   

(

)

(

)

(20)

and

{

(

)

( j) ( j) Λ( j) (s) = δW( j) D11W,xxx nx + D12 W,(xyyj) nx +W, xxy ny + D22W, (yyyj) ny + 2D66W,(xyyj) nx + 2D66W, (xxyj) ny n

) ∑c ( ∆W

(

− µk(f j) W,x( j) nx +W,(yj)ny + µ

ij , xnx

ij

i =1

)

(

)

( j) ( j) + ∆Wij, yny − Nxx W, x( j )nx − µ W,xxx nx +W, xyy nx  

∂ − Nyy W, (yj )ny − µ W, (xxyj ) ny +W, (yyyj ) ny  − δW( j ) D11W,(xxj)nx + D12W,(yyj)ny + 2D66W,(xyj)ny  ∂x

)}

(

(

( j) ( j) nx +W,xy ny +µNxx W,xx

) } − ∂∂y (δW ){D W ( j)

(

( j) 12 ,xx nx

){

(

+ D22W,(yyj)ny + 2D66W,(xyj) nx + µNyy W,(xyj) nx +W,(yyj) ny

)} (21)

12

here, nx and ny are components of unit normal vector on the boundary of element. It is found that the boundary term Λ(j)(s) vanishes for clamped edges. Further, if we take the same assumption as the previous researchers [24], i.e. M(NL) = M(nl) - µ∇2M(nl) = 0 along the simply edges, the boundary term can also be neglected for GS with simply supported edges.

Substituting Eq. (14) into Eq. (19), a set of simultaneous linear equations will be yielded for each multilayered element of MLGS as below

( K

stiff

 + [CvdW ] − Nxx [ Bbuck ]

) {dˆ} = {0} e

(22)

e

Where,

 K stiff

[ K (1) ] [0] . .  (2)  [0] [ K ] . .  . . .   =  . . . e  .  .  [0] [0] . .   [0] [0] . .

n (1,i)  (1 − δ1i )[C ]  i=1   −[C(1,2) ]   .   C = . [ vdW ]e  .    −[C(n−1,1) ]    −[C( n,1) ]  

∑

[0]   . [0] [0]  . .   . .   . . .  . [ K ( n −1) ] [0]   . [0] [ K ( n ) ] .

[0]

(23a)

   n   (1 − δ2i )[C(2,i ) ] . . . −[C(2,n−1) ] −[C(2,n) ]  i =1  . . . .   . . . .  . . . .  n  . . . (1 − δn−1i )[C( n−1,i) ] −[C(n−1,2) ] −[C( n−1,n) ]   i =1  n (n,n −1) ( n,i)  ( n,2) . . . ] (1 − δni )[C ] −[C ] −[C  i =1  −[C(1,2) ]

−[C(1,n−1) ]

. . .

−[C(1,n) ]

∑

∑

∑

13

(23b)

[ Bbuck ]e

K ij(J ) =

[ B (1) ] [0] . .  ( 2)  [0] [ B ] . .  . . .  = . . .  . .   [0] [0] . .  [0] . .  [0]

∫∫ D

11N i , xx N j , xx

[0]   . [0] [0]  . .   . .   . . .  . [ B ( n −1) ] [0]   . [0] [ B ( n ) ] .

[0]

(

(23c)

)

+ D12 N i ,xx N j , yy + N i , yy N j ,xx + D22N i , yy N j , yy + 4D66 N i ,xy N j ,xy

Ae

(

))

+ k f(J ) N i N j + µ N i ,x N j ,x + N i , y N j , y  dAe 

Cij( I , J ) = cIJ

(

∫∫  N

i

(

(23d)

)

N j + µ N i , x N j , x + N i , y N j , y  dAe 

(23e)

Ae

BijJ =

∫∫{N

i,x N j,x

(

)

(

)

}

+ µ Ni, xx N j , xx + Ni , xy N j , xy + β  Ni , y N j , y + µ Ni, xy N j , xy + Ni, yy N j , yy  dAe 

Ae

{dˆ} = { {dˆ } {dˆ } ... {dˆ } {dˆ } (1)

e

T e

(2)

T e

( n−1)

T

( n)

e

T e

(23f)

T

}

(23g)

Here, β =Nyy/Nxx is the load ratio and δij are the components of identity matrix.

After the assembling process, the global finite element equation will be written as

([ K ] + [C] − Nxx [ B])ML {d}ML = {0}

(24)

14

where, [K]ML, [C]ML, [B]ML are, respectively, global stiffness, inter-layer and buckling matrices of MLGS and {d}ML is the global nodal displacement vector. The eigen-values obtained from Eq. (24) are the buckling loads of biaxially compressed MLGS on polymer substrate.

4. Numerical Results and Discussions

4.1. Convergence and Validation Studies In order to verify the accuracy of the present finite element solution, buckling loads of monolayer and double-layered graphene sheets under in-plane compression are compared with existing molecular dynamics and analytical solutions in the literature. Since the reported elastic properties for graphene sheets are scattered [42,43], six types of graphene sheet are taken here from the literature [44-48]. The mechanical properties of these graphene sheets are presented in Table 1. In fact, determination of the exact elastic properties of GSs has been a challenging issue between the research communities during the recent decades. This leads to the range 0.5-1.67 Tpa for GS’s Young modulus and 0.149-1.144 for its Poisson ratio obtained from experimental measurements, atomistic simulations and theoretical studies based on inter-atomic potential fields [1,3,42,43,49]. For example, the differences between the amounts of elastic constants reported by theoretical studies come from the different types of potential fields taken by different researchers to obtain the lattice deformation under special loading conditions [49]. In this field, the so-called “Yakobson’s paradox” refers to the contradicting values of the elastic properties reported in the literature for CNTs and GSs [50]. By the way, in the current finite element analysis, we take the graphene layers in Table 1 since it is seen that some of them produce results which are in good agreement with atomistic simulations [46,47,48]. The finite element buckling solutions for different rectangular SLGSs and DLGSs made of these materials under various loading and boundary conditions are then compared with the reported analytical and MD solutions in Table 2 and 3, respectively. It is found from Table 2 that the 15

results calculated via Eq. (24) for a grid with 20×20 elements are in a good agreement with those obtained through the Navier’s approach [24,32] for both single and double layered graphene sheets. Hence, this grid is adopted to generate all the other results herein. Further, from Table 3 it is found that the current FE solutions are very close to the MD solutions for a biaxially compressed simply supported square monolayer armchair graphene sheet made of material (II) if the nonlocal parameter µ =1.85nm2 (reported by Ansari et al. [46]) is taken. The FE solutions obtained here are also very close to those obtained by Ansari et al. [46] through the Galerkin method. In addition, we found that the buckling loads of bridge GS under uniaxial compression will be in line with those atomistic solutions obtained by Sakhaee-Pour [51] if the material Armchair-I, -II and Zigzag-I, -II are taken. We obtained the nonlocal parameters 0.47nm2, 0.54nm2 and 0.65nm2 from matching the atomistic solutions with nonlocal solutions for Armchair-II, Zigzag-I, and Zigzag-II, respectively. These values are in line with the range 0.22-0.67nm reported by Shen et al. [48]. They showed that the range of values for nonlocal parameter depends obviously on the material properties taken for graphene sheets. Wang and Wang [52] also pointed out that the nonlocal parameter should be smaller than 2.0nm for carbon nano-tubes. Also, Ansari et al. [47] showed that the nonlocal parameter depends on the chirality of the GS and the edge conditions. They reported the value of nonlocal parameter in range of 0.87-1.41nm2 for zigzag and 0.71-1.34nm2 for armchair GSs. Since all nonlocal parameters in Table 3 are less than 2.00nm2, the range 0-2.00nm2 is adopted here to consider the small scale effect on the buckling behavior of all MLGSs.

It is found from Table 2 that the buckling loads of a DLGS, in general, are categorized into two distinct groups (N1, N2). One includes loads which are independent of interlayer interactions known as In-Phase (INPH) loads (N1), and the other includes loads which depend on the interlayer vdW interactions known as Out-of-Phase (OPH) loads (N2). For better understanding, buckling mode shapes of SLGS and DLGS under uniform biaxial compression are depicted in Fig. 3 for wave numbers m=n=1, m=2,n=1 and m=n=2. By comparing the INPH mode shapes of DLGS with the corresponding 16

ones of SLGS in this figure, the agreement of the INPH results for DLGS in Table 2 with those obtained for single layered graphene sheet will be better understood. Also, from Table 2 and Fig. 3 it is revealed that the OPH buckling loads are much higher in amount than the INPH loads due to the presence of strong interlayer vdW interactions.

Another point in Table 2 which should be noted is that the classical continuum theory ( µ = 0.00nm2) overestimates the critical buckling loads of GSs which can lead to deficient design of nanocomposites and NEMS utilizing MLGSs as load-bearing components. Consequently, it is established that nonlocal assumptions should be considered for accurately estimating the critical buckling loads of graphene sheets at very small dimensions. For better understanding the influence of non-locality, a comprehensive study is performed in the following sections to examine all the effects of material properties, size, buckling modes, interlayer interactions, boundary conditions and polymer substrate on the buckling loads of embedded MLGS. The focus of the paper is on the buckling behavior of DLGSs. However, some results are also presented for other multilayered graphene sheets.

4.2. Buckling Loads of Nonlocal Multilayered Graphene Sheets: Effect of nonlocal parameter

Taking the material properties in Table 1 and interlayer vdW coefficients from Eq. (10), buckling loads for six types of double layered graphene sheet under uniform biaxial compression are calculated from Eq. (24) and the results are presented in Table 4 for a=b=10nm and nonlocal parameters 0.00 and 2.00nm2. To show the intensity of nonlocal effect, the buckling load ratios (BLRs) are also presented in this table in the parentheses below the nonlocal results. Here, the BLR is the ratio of the results obtained by nonlocal continuum theory (µ=2.00 nm2) to those obtained by classical continuum theory (µ=0.00 nm2). As it is seen in this table the BLR is less than unity for nonlocal results. This means that the nonlocal results are always smaller than the classical ones. In fact, as the nonlocal parameter increases the buckling resistance of the equivalent nonlocal 17

continuum plate model of MLGS decreases. This is due to the small scale effect which is imposed to the system by assumptions in the nonlocal continuum theory. As mentioned before the classical continuum mechanics fails to take the account of such effects. It is also found from Table 4 that the largest and smallest values of buckling loads belong to Material-I and Armchair-I, respectively. Moreover, it is found from this table that as the mode number increases the INPH and OPH buckling load ratios decrease. This implies the importance of the small scale effect in higher modes. The nonlocal INPH buckling loads of Material-II, Armchair-II, Zigzag-I and Zigzag-II based on the introduced nonlocal parameters in Table 3 are respectively found to be 0.4860, 0.1003, 0.1042, 0.1064 nN.nm-1for m=n=1 and 1.0786, 0.3197, 0.3232, 0.3173 nN.nm-1 for m=n=2. Similarly, the nonlocal OPH buckling loads are found to be 853.1980, 852.8122, 852.8161, 852.8183 nN.nm-1 for m=n=1 and 214.2465, 213.4976, 213.5012, 213.4953 nN.nm-1 for m=n=2.

Buckling loads of triple-layered, four-layered and five-layered graphene sheets (TLGS, FLGS and FILGS) are tabulated in Table 5 and 6 for wave numbers m=n=1 and m=n=2, respectively. The associated mode shapes are depicted in Fig. 4a and 4b. It is seen in these figures that there are three distinct types of buckling mode shape for a TLGS, four types for FLGS and five types for FILGS. Consequently, it can be concluded that, generally, for an n-layered graphene sheet there are n distinct mode types among which only one is independent of vdW interactions (INPH mode). It is again seen here that the OPH buckling loads are much higher than the INPH loads due to the interlayer interactions and the magnitude of OPH buckling loads depends on the associated mode shape to which MLGS buckles. The increase in the share of vdW interactions between the layers of MLGS results in the higher OPH buckling loads. This is the reason by which one can understand why the amount of OPH load N5 is higher than N4 and N4 is higher than N3 and so on. From Tables 5 and 6, it is found again that similar to DLGS the buckling loads of TLGS, FLGS and FILGS decrease by increasing in nonlocal parameter. But this reduction is more noticeable for INPH loads with higher 18

mode numbers. Also as the share of interlayer interactions increases (higher OPH modes e.g. N3, N4 and N5) the small scale effect decreases. Moreover, it is seen in these tables that the amounts of OPH buckling load for all materials are close to each other in comparison to the INPH ones. Such trend is more noticeable when the wave number is low enough and the share of vdW interactions in the buckled shape of MLGS is high enough (higher OPH modes). This means that for such case it doesn’t matter what type of material is chosen for graphene layers because the effect of vdW interactions on the buckling loads is more prominent than the flexural stiffness of each graphene layer.

4.3. Buckling Loads of Nonlocal Double Layered Graphene Sheets: Effect of size and interlayer vdw interactions

To show the intensity of the small scale effect in different sizes, critical buckling loads of different square DLGSs are plotted versus the length in Fig. 5a and 5b. It is seen in these figures that the buckling loads and small scale effect decrease by increasing in the length of the sample. The reduction in the small scale effect in larger graphene sheets comes from the size dependent essence of the nonlocal elasticity theory. In fact, by increasing the external characteristic length of the square GS (here the length a), the small scale effect decreases while the internal characteristic length is assumed to be unchanged. By further increasing in the side length, the results converge to those obtained through classical theory of elasticity (µ=0.0nm2) for all materials. This is logical since for large enough systems the classical continuum theory is valid.

To show the influence of interlayer vdW interaction and nonlocal parameter on the results, buckling load ratios associated with OPH modes for a simply supported square DLGS with length 10nm are plotted against vdW coefficient c12 in Fig. 6a and 6b for nonlocal parameter 2.00nm2. The 19

figure shows that the nonlocal effect is more prominent for smaller values of interlayer vdW coefficient c12. This means that for a DLGS with strong enough vdW interlayer interaction, the classical continuum mechanics can be accurately used instead of the more complicated nonlocal elasticity theory. For the value c12=84.16GPa.nm-1 calculated from Eq. (10), the relative error percent due neglecting the nonlocal assumptions for Material-I, Material-II, Armchair-I, Armchair-II, Zigzag-I and Zigzag-II are respectively found to be %0.0396, %0.0220, %0.0030, %0.0036, %0.0038, %0.0039 for m=n=1 and %1.3440, %0.7529, %0.1041, %0.1256, %0.1321, %0.1376 for m=n=2. From Fig. 6 it is also found that the nonlocal effect is more noticeable for Material-I and Material-II.

It is found from the buckling analysis of MLGSs in Table 4, 5 and 6 that as the wave number increases, the buckling loads associated with OPH modes decrease. In order to show the influence of interlayer vdW interaction and nonlocal parameter on this trend, OPH buckling loads of a simply supported square DLGS made of Material-II are plotted in Fig. 7a and 7b for nonlocal parameter µ=0.00nm2 and µ=2.00nm2, respectively. It is seen in this figure that as the inter-layer interactions increase the buckling loads increase for all mode numbers. In addition, for smaller values of vdW coefficient, approximately for Log (c12) < 8.2168Gpa.m-1 in the local case (µ=0.00nm2), the mode shape with m=n=1 is the critical mode shape but for greater values of c12 the critical mode shape is different from m=n=1. By further increasing in the amount of c12, the order of critical mode shapes is completely reversed. For example, it is found from this figure that in the local case (µ=0.00nm2) the order of three mode shapes m=n=1, m=2,n=1 and m=n=2 changes in three points with vdW coefficients:

Log(C(11,21))=8.2168Gpa.m-1,

Log(C(11,22))=8.4188Gpa.m-1

and

Log(C(21,22))=8.8159Gpa.m-1. It is found that in the nonlocal case with µ=2.00nm2 the previous values

change

to

Log(C(11,21))=7.7708Gpa.m-1,

Log(C(11,22))=7.8645Gpa.m-1

and

Log(C(21,22))=8.1064Gpa.m-1. This means that by increasing in the nonlocal parameter the change in mode order happens for a DLGS with weaker interlayer interaction. The amounts of vdW coefficient in the change points for other graphene materials are tabulated in Table. 7.

20

4.4. Buckling Loads of Nonlocal DLGS: Effect of boundary conditions

To show the effect of different boundary conditions (BCs), critical buckling load ratios of a nonlocal DLGS are presented in Fig. 8 for different edge conditions. Here, symbolic notations are used for boundary conditions. For instance, the notation SCSC denotes a rectangular DLGS with all layers simply supported along y=0 and y=b and clamped along x=0 and x=a. The figure reveals that the nonlocal effect is more considerable for DLGS with more clamped edges. For example, the error percent through using the classical continuum mechanics for µ=2.00nm2 is respectively found to be % 28.30 and % 51.15 for SSSS and CCCC boundaries.

4.5. Buckling Loads of Nonlocal DLGS: Effect of polymer substrate and mode number

Considering the vdW interaction between carbon atoms of each graphene sheet of DLGS and – CH2– units of polymer matrix, critical buckling loads of six types of DLGS embedded in polymer matrix are plotted against the nonlocal parameter in Fig. 9a and 9b. It is seen in these figures that small scale effect reduces the buckling loads of embedded DLGSs. Also, it is seen in Fig. 9b that the largest amounts of buckling load belong to Zigzag-II.

To show the effect of polymer matrix on OPH buckling loads, local OPH buckling loads of freestanding square DLGS (without elastic foundation) and embedded DLGS are plotted in Fig. 10 for different sizes. It is seen in this figure that the buckling loads of embedded DLGS are greater than those in free case especially at larger sizes. It is also found from this figure that as the side length of DLGS increases the OPH buckling loads change from a descending trend to an ascending one. This is because of the effect of interlayer vdW interactions which increases during the increase in DLGS side length. In fact, for small lengths the vdW interlayer interaction has inconsiderable effect on the

21

results and similar to INPH loads the OPH buckling loads decrease by increasing in the side length. But after a specific value for the side length (approximately a > 5.00nm), the buckling loads increase with increasing in the length due to the great influence of interlayer vdW interaction. Similar trend is seen for embedded DLGS in Fig. 10. The increase in the difference between the two cases, i.e. freestanding and embedded, by length comes also from the greater influence of graphene-polymer vdW interaction in the larger DLGSs.

Effect of small scale in higher mode numbers is illustrated in Figs. 11a and 11b, respectively, for free-standing and embedded DLGS with the properties of Material-II. It is seen in these figures that the buckling load ratios in all mode numbers decrease by increasing the nonlocal parameter. From these figures it is also found that the small scale effect is more pronounced in higher buckling modes. So, it is concluded that the use of classical continuum mechanics in higher buckling modes can lead to unreliable results. Comparing the amounts of BLR associated with each mode number in Figs. 11a and 11b, it is found that the graphene-polymer interaction reduces the small scale effect on the buckling loads in these modes. To show the effect of graphene-polymer vdW interaction coefficient kf(2) on the results, critical buckling loads of DLGS are plotted against nonlocal parameter µ in Fig. 12 considering zero and nonzero values for kf(2). It is found from this figure that the effect of vdW coefficient kf(2) on the results is negligible. Based on this conclusion, the variations of critical buckling loads are plotted in Fig. 13 against vdW coefficient kf(1) for different nonlocal parameters to show the influence of graphene-polymer vdW interaction on the stability of embedded DLGSs. It is seen in this figure that the critical buckling loads of DLGS increase with increase in vdW coefficient kf(1). However, it is found that this ascending trend is more noticeable for smaller values of nonlocal parameter.

To clearly show the effect of polymer substrate and small scale, critical buckling load ratios of nonlocal embedded DLGS are plotted in Fig. 14 against graphene-polymer vdW coefficient kf(1). As it is seen in this figure, the intensity of small scale effect depends obviously on the amount of 22

graphene-polymer vdW interaction coefficient. In fact, the small scale effect on buckling loads of DLGS can decrease or increase with changing in the amount of polymer coefficient. Such behavior was also seen in the buckling analysis of embedded single layer graphene sheet [53]. The reason for this behavior is the change in the configuration of the buckled graphene sheet with increasing the value of graphene-polymer coefficient kf(1). In fact, as the graphene-polymer coefficient increases from zero an increase in the buckling load ratio is observed which means that the small scale effect decreases. By further increase in kf(1) to a specific value, the DLGS buckles into a higher mode shape. This change in the mode shape to a higher mode, as it was shown previously in Figs. 11a and 11b, leads to an increase in the small scale effect and reduces the buckling load ratios.

5. Conclusions In the present study, finite element buckling analysis of embedded MLGS has been carried out using the size dependent Eringen’s nonlocal elasticity theory. The graphene-graphene and graphenepolymer vdW interactions are modeled as linear springs which their stiffnesses are obtained from the Lennard–Jones interatomic potential. A bicubic Hermitian plate element with four nodes and total sixteen degrees of freedom has been adopted due to its swift convergence. The accuracy of the current FE model has been verified through a comparison with molecular dynamics and analytical solutions available in the literature. From the results it is concluded that the small scale has a decreasing effect on both INPH and OPH buckling loads of MLGS however this reduction is more pronounced for INPH buckling loads with higher mode numbers. As the intensity of interlayer interaction increases the small scale effect decreases. Also, it is seen that for low enough wave numbers, strong enough interlayer interactions and higher OPH modes the change in the material properties, chosen for graphene layers, has inconsiderable effect on the amount of buckling loads because of the much stronger influence of interlayer interactions in comparison to the flexural stiffness of graphene sheets. Further, it is found that the bottom graphene-polymer vdW interaction 23

together with the nonlocal parameter has a remarkable effect on the critical buckling loads of embedded DLGS.

References

[1] [2] [3] [4]

[5] [6]

[7] [8]

[9]

C. Lee, X. Wei, J.W. Kysar, J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science 321 (2008) 385–388. K.N. Kudin, G.E. Scuseria, B.I. Yakobson, C2F, BN and C nanoshell elasticity from ab initio computations, Phys. Rev. B 64 (2001) 235406. A. Sakhaee-Pour, Elastic properties of single-layered graphene sheet, Solid State Commun. 149 (2009) 91-95. R. Arsat, M. Breedon, M. Shafiei, P.G. Spizziri, S. Gilje, R.B. Kaner, K. Kalantar-Zadeh, W. Wlodarski, Graphene-like nanosheets for surface acoustic wave gas sensor applications, Chem. Phys. Lett. 467(4/6) (2009) 344-347. T. Kuilla, S. Bhadra, D. Yao, N.H. Kim, S. Bose, J.H. Lee, Recent advances in graphene based polymer composites, Prog. Polym. Sci. 35 (2010) 1350-1375. P. Zhua, Z.X. Leia, K.M. Liew, Static and free vibration analyses of carbon nanotube-reinforced composite plates using finite element method with first order shear deformation plate theory, Compos. Struct. 94(4) (2012) 1450–1460. T.J. Chuang, P.M. Anderson, M.K. Wu, S. Hsieh, Nanomechanics of materials and structures, Springer, Dordrecht the Netherlands, 2006. A. Sears, R.C. Batra, Macroscopic properties of carbon nanotubes from molecular-mechanics simulations, Phys. Rev. B 69(23) (2004) 235406. M.R. Ayatollahi, S. Shadlou, M.M. Shokrieh, Multiscale modeling for mechanical properties of carbon nanotube reinforced nanocomposites subjected to different types of loading, Compos. Struct. 93(9) (2011) 2250–2259.

[10] T. Chang, H. Gao, Size-dependent elastic properties of a single-walled carbon nanotube via a molecular mechanics model, J. Mech. Phys. Solids 51 (2003) 1059-1074. [11] N.A. Fleck, J.W. Hutchinson, Strain gradient plasticity, Adv. Appl. Mech. 33 (1997) 296–358. [12] F. Yang, A.C.M. Chong, D.C.C. Lam, P. Tong, Couple stress based strain gradient theory for elasticity, Int. J. Solids Struct. 39 (2002) 2731–2743. [13] P. Lu, L.H. He, H.P. Lee, C. Lu, Thin plate theory including surface effects, Int. J. Solids Struct. 43 (2006) 4631-4647. [14] A.C. Eringen, On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves, J. App. Phys. 54 (1983) 4703-4710. [15] A.C. Eringen, Nonlocal Continuum Field Theories, Springer, NewYork, 2002. [16] Ö. Civalek, C. Demir, Bending analysis of microtubules using nonlocal Euler–Bernoulli beam theory, Appl. Math. Model. 35 (2011) 2053–2067.

24

[17] Ö. Civalek, B. Akgöz, Free vibration analysis of microtubules as cytoskeleton components: nonlocal euler-bernoulli beam modeling, Sci. Iran 17(5) (2010) 367-375. [18] Ö. Civalek, C. Demir, B. Akgöz, Free vibration and bending analyses of cantilever microtubules based on nonlocal continuum model, Math. Comput. Appl. 15 (2010) 289-298. [19] M. Simsek, Vibration analysis of a single-walled carbon nanotube under action of a moving harmonic load based on nonlocal elasticity theory, Physica E 43(1) (2010) 182-191. [20] M. Simsek, Forced vibration of an embedded single-walled carbon nanotube traversed by a moving load using nonlocal Timoshenko beam theory, Steel Compos. Struct. 11(1) (2011) 5976. [21] A. Benzair, A. Tounsi, A. Besseghier, H. Heireche, N. Moulay, L. Boumia, The thermal effect on vibration of single-walled carbon nanotubes using nonlocal Timoshenko beam theory, J. Phys. D: Appl. Phys. 41 (2008) 225404. [22] K. Amara, A. Tounsi, I. Mechab, E.A. Adda-Bedia, Nonlocal elasticity effect on column buckling of multiwalled carbon nanotubes under temperature field, Appl. Math. Model. 34 (2010) 3933–3942. [23] L.J. Sudak, Column buckling of multiwalled carbon nanotubes using nonlocal continuum mechanics, J. Appl. Phys. 94 (2003) 7281-7287. [24] T. Murmu, S.C. Pradhan, Buckling of biaxially compressed orthotropic plates at small scales, Mech. Res. Commun. 36 (2009) 933–938. [25] T. Aksencer, M. Aydogdu, Levy type solution method for vibration and buckling of nanoplates using nonlocal elasticity theory, Physica E 43 (2011) 954–959. [26] S.C. Pradhan, Buckling analysis and small scale effect of biaxially compressed graphene sheets using non-local elasticity theory, Sadhana 37 (2012) 461-480. [27] E. Jomehzadeh, A.R. Saidi, A study on large amplitude vibration of multilayered graphene sheets, Comput. Mater. Sci. 50 (2011) 1043–1051. [28] E. Jomehzadeh, A.R. Saidi, N.M. Pugno, Large amplitude vibration of a double layered graphene sheet embedded in a nonlinear polymer matrix, Physica E 44(10) (2012) 1973-1982. [29] E. Jomehzadeh, A.R. Saidi, Decoupling the nonlocal elasticity equations for three dimensional vibration analysis of nano-plates, Compos. Struct. 93 (2011) 1015–1020. [30] H. Babaie, A.R. Shahidi, Vibration of quadrilateral embedded multilayered graphene sheets based on nonlocal continuum models using the Galerkin method, Acta Mech. Sinica 27(6) (2011) 967-976. [31] H. Babaie, A.R. Shahidi, Small-scale effects on the buckling of quadrilateral nanoplates based on nonlocal elasticity theory using the Galerkin method, Arch. Appl. Mech. 81(8) (2011) 10511062. [32] S.C. Pradhan, J.K. Phadikar, Scale effect and buckling analysis of multilayered graphene sheets based on nonlocal continuum mechanics, J. Comput. Theor. Nanos. 7 (2010) 1948-1954. [33] T. Murmu, J. Sienz, S. Adhikari, C. Arnold, Nonlocal buckling of double-nanoplate-systems under biaxial compression, Compos. Part B: Eng. 44(1) (2013) 84–94. [34] T. Murmu, J. Sienz, S. Adhikari, C. Arnold, Nonlocal buckling behavior of bonded doublenanoplate-systems, J. Appl. Phys. 110(8) (2011) 084316-9. [35] R. Szilard, Theories and applications of plate analysis: classical, numerical and engineering methods, John Wiley & Sons, inc., Hoboken, New Jersey, 2004. [36] J.K. Phadikar, S.C. Pradhan, Variational formulation and finite element analysis for nonlocal elastic nanobeams and nanoplates, Comput. Mater. Sci. 49 (2010) 492-499. 25

[37] S. Adali, Variational principles for nonlocal continuum model of orthotropic graphene sheets embedded in an elastic medium, Acta Math. Sci. 32 (1) (2012) 325–338. [38] R. Ansari, R. Rajabiehfard, B. Arash, Nonlocal finite element model for vibrations of embedded multi-layered graphene sheets, Comput. Mater. Sci. 49 (2010) 831–838. [39] J.E. Lennard-Jones, The determination of molecular fields: from the variation of the viscosity of a gas with temperature, Proc. R. Soc. London 106 (1924) 441-462. [40] W.B. Lu, J. Wu, J. Song, K.C. Hwang, L.Y. Jiang, Y. Huang, A cohesive law for interfaces between multi-wall carbon nanotubes and polymers due to the van der Waals interactions, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3261–3267. [41] D.O. Brush, B.O. Almroth, Buckling of bars, plates, and shells, McGraw-Hill, New-York, 1975. [42] F. Scarpa, S. Adhikari, A. Srikantha Phani, Effective elastic mechanical properties of single layer graphene sheets, Nanotechnology 20 (2009) 065709. [43] S.K. Georgantzinos, G.I. Giannopoulos, N.K. Anifantis, Numerical investigation of elastic mechanical properties of graphene structures, Mater. Des. 31 (2010) 4646–4654. [44] K. Behfar, R. Naghdabadi, Nanoscale vibrational analysis of a multi-layered graphene sheet embedded in an elastic medium, Compos. Sci. Technol. 65 (2005) 1159–1164. [45] S. Kitipornchai, X.Q. He, K.M. Liew, Continuum model for the vibration of multilayered graphene sheets, Phys. Rev. B 72 (2005) 075443. [46] R. Ansari, H. Rouhi, Explicit analytical expressions for the critical buckling stresses in a monolayer graphene sheet based on nonlocal elasticity, Solid State Commun. 152 (2012) 56– 59. [47] R. Ansari, S. Sahmani, B. Arash, Nonlocal plate model for free vibrations of single-layered graphene sheets, Phys. Lett. A 375 (2010) 53–62. [48] L. Shen, H.S. Shen, C.L. Zhang, Nonlocal plate model for nonlinear vibration of single layer graphene sheets in thermal environments, Comput. Mater. Sci. 48 (2010) 680–685. [49] Q. Lu, R. Huang, Nonlinear mechanics of single-atomic-layer graphene sheets, Int. J. Appl. Mech. 1(3) (2009) 443–467. [50] R. Chowdhury, S. Adhikari, F. Scarpa, M.I. Friswell, Transverse vibration of single-layer graphene sheets, J. Phys. D: Appl. Phys. 44 (2011) 205401. [51] A. Sakhaee-Pour, Elastic buckling of single-layered graphene sheet, Comput. Mater. Sci. 45 (2009) 266–270. [52] Q. Wang, C.M. Wang, The constitutive relation and small scale parameter of nonlocal continuum mechanics for modelling carbon nanotubes, Nanotechnology 18 (2007) 075702. [53] S.C. Pradhan, T. Murmu, Small scale effect on the buckling analysis of single-layered graphene sheet embedded in an elastic medium based on nonlocal plate theory, Physica E 42 (2010) 1293–1301.

26

Figure Captions

Fig. 1: Schematics of an embedded MLGS: (a) MLGS under biaxial compression; (b) MLGS stacked up from single graphene layers; (c) Continuum plate models of graphene layers and the graphene-graphene and graphene-polymer interactions in MLGS (LN: Layer Number). Fig. 2: Bicubic Hermitian rectangular plate element with total sixteen degrees of freedom. Fig. 3: Buckling mode shapes for SLGS and DLGS: (a) m=n=1; (b) m=2, n=1; (c) m=n=2. Fig. 4: Buckling mode shapes for TLGS, FLGS and FILGS: (a) m=n=1; (b) m=n=2. Fig. 5: Nonlocal effect on buckling loads of DLGS in different sizes: (a) Material-I and –II; (b) Armchair-II and Zigzag-II. Fig. 6: Effect of vdW interaction coefficient c12 on nonlocal buckling loads of DLGS; (a) m=n=1; (b) m=n=2. Fig. 7: Effect of graphene-graphene vdW interaction coefficient c12 on the order of OPH modes of DLGS (Material-II); (a) µ=0.00nm2; (b) µ=2.00nm2. Fig. 8: Effect of different boundary conditions on critical buckling loads of nonlocal DLGS.

27

Fig. 9: Small scale effect on the critical buckling loads of DLGS in the presence of polymer matrix: (a) Material-I, -II; (b) Armchair-I, -II and Zigzag-I, -II. Fig. 10: The influence of polymer matrix on the OPH buckling loads of DLGS for different sizes (Material-II). Fig. 11: Small scale effect on the buckling load ratios of free-standing and embedded DLGS in higher mode numbers: (a) Free-standing ; (b) Embedded. Fig. 12: The influence of vdW interaction between the top graphene layer of DLGS and the polymer matrix on the amount of critical buckling loads (Material-II). Fig. 13: Effect of bottom graphene-polymer vdW interaction coefficient kf (1) on the critical buckling loads of DLGS (Material-II). Fig. 14: Effect of graphene-polymer coefficient and small scale on the critical buckling load ratios of embedded DLGS (Material-II).

Table Captions

Table 1: Material properties of graphene sheets

Table 2: Convergence and comparison of the buckling loads (nN.nm-1) of simply supported square double layered orthotropic graphene sheet under uniform biaxial compression (a=10.2nm, c12=45Gpa.nm-1) 28

Table 3: Critical buckling loads of monolayer armchair and zigzag graphene sheets obtained from the current finite element analysis and atomistic simulations (results from a grid with 20×20 elements)

Table 4: Local and nonlocal buckling loads for six types of simply supported square double layered graphene sheets (a=10nm, c12=84.16 Gpa.nm-1)

Table 5: Local and nonlocal buckling loads for six types of simply supported square triple-layered, four-layered and five-layered graphene sheets (a=10nm, m=n=1)

Table 6: Local and nonlocal buckling loads for six types of simply supported square triple-layered, four-layered and five-layered graphene sheets (a=10nm, m=n=2)

Table 7: The values of interlayer vdW coefficient in three change points for each graphene material

Table 1

Material

E1(Gpa)

E2(Gpa)

ν12

G12(Gpa)

ρ (g.cm-3)

h (nm)

References

I

1765

1588

0.300

678.85

2.300

0.340

[24], [32], [44]

II

1000

1000

0.160

431.03

2.250

0.340

[45], [46], [47]

29

Armchair-I

2434

2473

0.197

1039

6.316

0.129

[27], [28], [48]

Armchair-II

2154

2168

0.202

923

5.727

0.143

[27], [28], [48]

Zigzag-I

2145

2097

0.223

938

5.624

0.145

[27], [28], [48]

Zigzag-II

2067

2054

0.204

913

5.482

0.149

[27], [28], [48]

Table 2

m n

µ 2

(nm )

Number of Elements

Navier’s Solution

5×5

10×10

15×15

18×18

20×20

1 1

1.1498

1.1497

1.1497

1.1497

1.1497

1.1497 a,b

1 2

2.7744

2.7704

2.7701

2.7701

2.7701

2.7701 a,b

2 1

2.9549

2.9499

2.9497

2.9496

2.9496

2.9496 a,b

2 2

4.6015

4.5992

4.5990

4.5989

4.5989

4.5989 a,b

1 1

0.6537

0.6537

0.6537

0.6537

0.6537

0.6537 a,b

1 2

0.9566

0.9561

0.9561

0.9561

0.9561

0.9561 a,b

2 1

1.0188

1.0181

1.0181

1.0181

1.0181

1.0181 a,b

2 2

1.1401

1.1398

1.1396

1.1396

1.1396

1.1396 a,b

1 1

475.9100

475.5249

475.5166

475.5158

475.5155

475.5153b

1 2

192.6092

192.5228

192.5192

192.5179

192.5170

192.5163 b

2 1

192.8131

192.7054

192.6991

192.6977

192.6964

192.6958 b

2 2

123.7073

123.2406

123.2292

123.1911

123.1906

123.1903 b

N1 0.00

4.00

N2

0.00

30

4.00

1 1

475.3190

475.0278

475.0202

475.0195

475.0194

475.0192 b

1 2

190.6870

190.6935

190.7000

190.7011

190.7019

190.7023 b

2 1

190.7493

190.7555

190.7620

190.7631

190.7640

190.7643 b

2 2

120.0114

119.7355

119.7320

119.7317

119.7312

119.7309 b

a Values for single layered graphene sheet of Ref. [24]. b Values for double layered graphene sheet of Ref. [32].

Table 3

Materials

a×b

Atomistic Simulations

(nm×nm)

Material-II

Material-II

2

µ=0.00nm

2

µ=1.85nm

4.99×4.99

2.6647

1.0803

(1.0803) *

1.0837a

8.08×8.08

1.0163

0.6518

(0.6517)

0.6536 a

10.77×10.77

0.5720

0.4351

(0.4350)

0.4331 a

14.65×14.65

0.3092

0.2642

(0.2642)

0.2609 a

18.51×18.51

0.1937

0.1750

(0.1750)

0.1714 a

22.35×22.35

0.1328

0.1238

(0.1238)

0.1191 a

26.22×26.22

0.0965

0.0916

(0.0916)

0.0889 a

30.04×30.04

0.0735

0.0707

(0.0707)

0.0691 a

31

33.85×33.85

0.0579

0.0561

(0.0561)

0.0554 a

37.81×37.81

0.0464

0.0453

(0.0452)

0.0449 a

41.78×41.78

0.0380

0.0372

(0.0372)

0.0372 a

45.66×45.66

0.0318

0.0313

(0.0313)

0.0315 a

Armchair-I 2

Armchair-II 2

µ=0.00nm

µ=0.47nm

10×10

0.1779

0.1811

0.1854 b

10×20

0.1785

0.1817

0.1892 b

20×10

0.0463

0.0510

0.0488 b

20×20

0.0479

0.0513

0.0499 b

Zigzag-I

Zigzag-II 2

2

µ=0.54nm

µ=0.65nm

10×10

0.1850

0.1856

0.1917 c

10×20

0.1857

0.1862

0.1957 c

20×10

0.0530

0.0534

0.0503 c

20×20

0.0533

0.0538

0.0514 c

* the values in parentheses are those reported in [46] through the Galerkin method. a Molecular dynamics solutions for simply supported armchair monolayer graphene sheets under biaxial compression taken from [46]. b Values for bridge monolayer armchair graphene sheet under uniaxial compression taken from [51]. c Values for bridge monolayer zigzag graphene sheet under uniaxial compression taken from [51].

32

Table 4

Materials

µ m n

2

(nm )

I

II

Armchair-I

Armchair-II

Zigzag-I

Zigzag-II

0.00

1.1962

0.6635

0.0907

0.1096

0.1153

0.1200

2.00

0.8576

0.4757

0.0651

0.0786

0.0826

0.0861

(0.7169)*

(0.7170)

(0.7178)

(0.7172)

(0.7164)

(0.7175)

0.00

853.9082

853.3754

852.8027

852.8215

852.8272

852.8320

2.00

853.5696

853.1876

852.7771

852.7905

852.7946

852.7980

(0.9996)

(0.9997)

(1.0000)

(0.9999)

(0.9999)

(0.9999)

0.00

4.7847

2.6540

0.3630

0.4383

0.4611

0.4802

2.00

1.8552

1.0290

0.1407

0.1699

0.1788

0.1862

(0.3877)

(0.3877)

(0.3876)

(0.3876)

(0.3878)

(0.3878)

0.00

217.9627

215.8320

213.5410

213.6163

213.6390

213.6581

2.00

215.0331

214.2070

213.3187

213.3479

213.3567

213.3641

(0.9866)

(0.9925)

(0.9990)

(0.9987)

(0.9987)

(0.9986)

INPH (N1) 1

1

OPH (N2) 1

1

INPH (N1) 2

2

OPH (N2)

2

2

* the values in parentheses are the buckling load ratios.

33

Table 5

Materials

µ Buckling Load Type

2

(nm )

I

II

Armchair-I

Armchair-II

Zigzag-I

Zigzag-II

0.00

1.1962

0.6635

0.0907

0.1096

0.1153

0.1200

2.00

0.8576

0.4757

0.0651

0.0786

0.0826

0.0861

0.00

427.5521

427.0194

426.4467

426.4655

426.4712

426.4760

2.00

427.2134

426.8315

426.4209

426.4344

426.4385

426.4419

0.00

1280.2641

1279.7314

1279.1586

1279.1775

1279.183 2

1279.1879

2.00

1279.9252

1279.5433

1279.1326

1279.1461

1279.150 2

1279.1537

0.00

1.1962

0.6635

0.0907

0.1096

0.1153

0.1200

2.00

0.8576

0.4757

0.0651

0.0786

0.0826

0.0861

0.00

250.9497

250.4170

249.8442

249.8631

249.8688

249.8736

TLGS N1

N2

N3

FLGS N1

N2

34

N3

N4

2.00

250.6111

250.2291

249.8185

249.8320

249.8361

249.8396

0.00

853.9081

853.3754

852.8027

852.8215

852.8272

852.8320

2.00

853.5693

853.1874

852.7768

852.7903

852.7943

852.7978

0.00

1456.8665

1456.3338

1455.7611

1455.7799

1455.785 6

1455.7904

2.00

1456.5276

1456.1457

1455.7350

1455.7485

1455.752 6

1455.7560

0.00

1.1962

0.6635

0.0907

0.1096

0.1153

0.1200

2.00

0.8576

0.4757

0.0651

0.0786

0.0826

0.0861

0.00

164.0496

163.5170

162.9442

162.9631

162.9687

162.9735

2.00

163.7110

163.3291

162.9185

162.9320

162.9360

162.9395

0.00

590.4056

589.8729

589.3002

589.3190

589.3247

589.3295

2.00

590.0669

589.6850

589.2743

589.2878

589.2919

589.2954

0.00

1117.4106

1116.8779

1116.3051

1116.3240

1116.329 7

1116.3344

2.00

1117.0717

1116.6898

1116.2792

1116.2927

1116.296 8

1116.3002

0.00

1543.7666

1543.2339

1542.6611

1542.6800

1542.685 6

1542.6904

2.00

1543.4276

1543.0457

1542.6351

1542.6486

1542.652 6

1542.6561

FILGS

N1

N2

N3

N4

N5

35

Table 6

Materials

µ Buckling Load Type

2

(nm )

I

II

Armchair-I

Armchair-II

Zigzag-I

Zigzag-II

0.00

4.7847

2.6540

0.3630

0.4383

0.4611

0.4802

2.00

1.8552

1.0290

0.1407

0.1699

0.1788

0.1862

0.00

111.2778

109.1153

106.7928

106.8692

106.8921

106.9115

2.00

107.7685

106.9396

106.0481

106.0774

106.083

106.0937

0.00

324.1235

321.9645

319.6432

319.7196

319.7425

319.7619

2.00

319.5824

318.7535

317.8620

317.8913

317.9002

317.9076

0.00

4.7847

2.6540

0.3630

0.4383

0.4611

0.4802

2.00

1.8552

1.0290

0.1407

0.1699

0.1788

0.1862

0.00

67.2034

65.0342

62.7101

62.7864

62.8094

62.8288

2.00

63.9004

63.0715

62.1800

62.2093

62.2182

62.2256

0.00

217.6993

215.5396

213.2180

213.2944

213.3173

213.3367

2.00

213.6755

212.8465

211.9551

211.9844

211.9932

212.0007

0.00

368.2061

366.0473

363.7260

363.8023

363.8253

363.8447

2.00

363.4505

362.6216

361.7301

361.7594

361.7683

361.7757

TLGS N1

N2

N3

FLGS N1

N2

N3

N4

FILGS

36

N1

N2

N3

N4

N5

0.00

4.7847

2.6540

0.3630

0.4383

0.4611

0.4802

2.00

1.8552

1.0290

0.1407

0.1699

0.1788

0.1862

0.00

45.5378

43.3456

41.0185

41.0949

41.1178

41.1372

2.00

42.3144

41.4855

40.5940

40.6233

40.6322

40.6396

0.00

151.9264

149.7656

147.4436

147.5200

147.5429

147.5623

2.00

148.2214

147.3925

146.5010

146.5303

146.5391

146.5466

0.00

283.4730

281.3138

278.9924

279.0688

279.0917

279.1111

2.00

279.1296

278.3006

277.4092

277.4385

277.4473

277.4548

0.00

389.8976

387.7388

385.4176

385.4940

385.5169

385.5363

2.00

385.0365

384.2076

383.3161

383.3454

383.3543

383.3617

Table 7

Materials

I

µ

Log (C(11,21))

Log (C(11,22))

Log (C(21,22))

(nm2)

(Gpa.m-1)

(Gpa.m-1)

(Gpa.m-1)

0.00

8.4911

8.6753

9.0512

2.00

8.0544

8.1179

8.3083

37

II

Armchair-I

Armchair-II

Zigzag-I

Zigzag-II

0.00

8.2168

8.4188

8.8159

1.85

7.8012

7.8942

8.1444

2.00

7.7708

7.8645

8.1064

0.00

7.3390

7.5553

7.9739

2.00

6.9025

7.0028

7.2632

0.00

7.4318

7.6371

8.0334

0.47

7.2911

7.4520

7.8022

2.00

6.9827

7.0793

7.3365

0.00

7.4541

7.6613

8.0577

0.54

7.3025

7.4626

7.8054

2.00

7.0082

7.1038

7.3519

0.00

7.4641

7.6778

8.0746

0.65

7.2930

7.4439

7.7908

2.00

7.0204

7.1211

7.3764

38

Fig. 1 (a)

Fig. 1 (b) 39

Fig. 1 (c)

40

Fig. 2

Fig. 3 (a)

41

Fig. 3 (b)

Fig. 3 (c)

42

Fig. 4 (a)

43

Fig. 4 (b)

Fig. 5 (a)

44

Fig. 5 (b)

45

Fig. 6 (a)

46

Fig. 6 (b)

47

Fig. 7 (a)

48

Fig. 7 (b)

49

Fig. 8

50

Fig. 9 (a)

51

Fig. 9 (b)

52

Fig. 10

53

Fig. 11 (a)

Fig. 11 (b) 54

Fig. 12

55

Fig. 13

56

Fig. 14

58