Postbuckling analysis of multi-layered graphene sheets under non-uniform biaxial compression

Physica E 47 (2013) 197–206

Contents lists available at SciVerse ScienceDirect

Physica E journal homepage: www.elsevier.com/locate/physe

Postbuckling analysis of multi-layered graphene sheets under non-uniform biaxial compression Ali Farajpour a,n, Alireza Arab Solghar b, Alireza Shahidi a a b

Department of Mechanical Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran Department of Mechanical Engineering, School of Engineering, Vali-e-Asr University of Rafsanjan, Iran

H I G H L I G H T S c c c c c

Nonlinear stability of multi-layered graphene sheets is studied. Graphene sheet is subjected to non-uniform in-plane loading through its thickness. Small scale and nonlinearity effects are taken into consideration. Six types of armchair and zigzag graphene sheets are investigated. Postbuckling loads are intensely sensitive to non-uniform and nonlocal parameters.

a r t i c l e i n f o

abstract

Article history: Received 22 July 2012 Received in revised form 16 September 2012 Accepted 26 October 2012 This article is lovingly dedicated to the memory of Mehdi Farajpour that his life was based on good thoughts, good words and good deeds. Ouderj’s mountains will never forget his amazing hunting stories. Available online 7 November 2012

In this article, the nonlinear buckling characteristics of multi-layered graphene sheets are investigated. The graphene sheet is modeled as an orthotropic nanoplate with size-dependent material properties. The graphene film is subjected by non-uniformly distributed in-plane load through its thickness. To include the small scale and the geometrical nonlinearity effects, the governing differential equations are derived based on the nonlocal elasticity theory in conjunction with the von Karman geometrical model. Explicit expressions for the postbuckling loads of single- and double-layered graphene sheets with simply supported edges under biaxial compression are obtained. For numerical results, six types of armchair and zigzag graphene sheets with different aspect ratio are considered. The present formulation and method of solution are validated by comparing the results, in the limit cases, with those available in the open literature. Excellent agreement between the obtained and available results is observed. Finally, the effects of nonlocal parameter, buckling mode number, compression ratio and non-uniform parameter on the postbuckling behavior of multi-layered graphene sheets are studied. & 2012 Elsevier B.V. All rights reserved.

Keywords: Nonlinear stability Graphene sheet Biaxial compression Nonlocal plate theory

1. Introduction In recent years, nanostructural elements such as carbon nanotubes (CNTs) [1], gold nanorings (GNRs) [2] and graphene sheets (GSs) [3] are commonly used as components in the micro electro-mechanical systems (MEMS) and nano electro-mechanical systems (NEMS) due to their superior mechanical, chemical and electronic properties. There are many works using fully atomistic approaches that investigated the static and dynamic behaviors of graphene membranes [4–8]. In contrast to statistical considerations, some researchers used the classical plate theory to investigate the vibration and buckling of carbon nanosheets [9,10].

n

Corresponding author. Tel.: þ9809138633720 E-mail addresses: [email protected], [email protected] (A. Farajpour). 1386-9477/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physe.2012.10.028

However, there are recent research works that pointing out to the failure of classical continuum based methods to study the mechanical properties of nanostructures [11,12]. The traditional continuum theory is a scale free theory and thus cannot predict the static and dynamic characteristics of nanomaterials properly. Therefore, the use of classical continuum models may be uncertain in the analysis of structural elements in nanoscale such as graphene sheets. Recently, modified continuum mechanics approach has been widely used to study the mechanical behavior of nanostructure [13–19]. The main reason for this is that controlled experiments on nanoscale are difficult to perform. On the other hand, molecular dynamic (MD) simulations are highly computationally expensive especially in the analysis of nanomaterials with large numbers of atoms inside them. It is well-known that the size effect plays a prominent role in the mechanical characteristics of structural elements at small scales. As the dimensions of these structures are reduced, the

198

A. Farajpour et al. / Physica E 47 (2013) 197–206

effects of inter-molecular and inter-atomic interactions must be considered to predict mechanical behavior properly. Thus, it is necessary to modify the classical continuum theory to include small scale effects. In Eringen’s nonlocal elasticity theory [20,21], the small scale effects are taken into consideration by assuming the stress tensor at a given point as a function not only of the strain tensor at that point but also a function of the strain tensors at all other points of entire domain occupied by the material. The classical elasticity theory is a typical case of the nonlocal elasticity theory in which stress state at an arbitrary point depends only on the strain state at that point. Chen et al. [22] investigated microcontinuum field theories such as couple stress theory, micromorphic theory, nonlocal theory, Cosserat theory, etc., employing lattice dynamics and molecular dynamics (MD) simulations. They reported that the nonlocal continuum models are reasonable from the physical and atomic points of view. Further, it has been shown that the nonlocal elasticity theory is quite accurate and reliable for the mechanical analysis of carbon nanotubes [23] and graphene sheets [24,25] by comparing the obtained results with those of molecular dynamics (MD) modeling. Understanding the vibration and buckling behaviors of GSs is essential for proper design of many MEMS and NEMS devices which utilize graphene sheets like resonators [3], mass sensors and atomistic dust detectors [26]. Behfar and Naghdabadi [9] used the classical elasticity theory to determine the natural frequencies and associated mode shapes of orthotropic multi-layered graphene sheets (MLGS) embedded in an elastic medium. They obtained van der Waals (vdW) interaction coefficient between two adjacent layers of MLGS using the Lennard–Jones model and the principle of virtual work. In addition, He et al. [10] investigated the nanoscale vibrational analysis of MLGSs used as nanoscale resonators. They presented explicit formula for the linear vdW interaction between all layers of GS. Duan and Wang [27] obtained an exact closed-form solution for the axisymmetric bending of circular single-layered graphene sheet (SLGS) via the nonlocal continuum plate model. In another work, molecular dynamics method and classical plate theory were employed to develop an equivalent continuum model for large deflection of circular SLGS with a concentrated force exerted at the centre [28]. Aghababaei and Reddy [29] developed a nonlocal third-order shear deformation plate theory for the linear vibration and bending of rectangular nanoplates. Further, Eringen’s nonlocal elasticity theory has been used to study the in-plane vibration [30] and thermal buckling [31] of nanoplates. Farajpour et al. [32] investigated the buckling behavior of nanoscale circular plates under uniform radial compression. They also studied the linear buckling response of orthotropic SLGS under linearly varying inplane load via the nonlocal elasticity theory [33]. An analytical method for determining the natural frequencies of nonlocal double-nanoplate-system (NDNPS) was introduced by Murmu and Adhikari [34]. Based on the nonlocal elasticity theory, elastic wave propagation in a SLGS supported by an elastic foundation was also investigated [35]. Carbon nanotubes and graphene sheets are always found in the form of single-layered or multi-layered nanostructures. A double or multi-layered graphene sheets (DLGSs/MLGSs) is composed of two or more layers of SLGS which are bonded by van der Waals (vdW) interaction. Nonlocal effects on the linear vibration and buckling of MLGSs were discussed by some researchers [36–38]. Recently, a nonlinear elastic plate model without including the small scale effects has been developed to analyze the vibration behavior of MLGSs [39]. Shen et al. [40] studied the nonlinear vibration response of SLGSs in thermal environments employing the nonlocal plate model and von Karman geometrical nonlinearity. They used MD simulation to

determine the nonlocal scaling effect parameter and anisotropic size-dependent material properties. In their work, numerical results were obtained for six types of armchair and zigzag SLGSs with three different values of aspect ratio. Jomehzadeh et al. [41] investigated the large amplitude vibration of DLGSs embedded in a nonlinear polymer matrix. To the best of authors’ knowledge, up to now the nonlinear stability of MLGSs using the nonlocal elasticity theory is not studied in the open literature. In addition, the graphene’s layers may be subjected to different in-plane edge loading. This motivates us to investigate the postbuckling behavior of multi-layered graphene films subjected to non-uniform edge loading through its thickness. Using nonlocal elasticity theory and von Karman’s strain-displacement relations, nonlinear governing equations of MLGSs are derived. Nonlinear vdW interactions are taken into consideration employing Lennard–Jones pair potential. Galerkin method is used to obtain the postbuckling loads of MLGS with all edges simply supported. The effect of small size and nonlinearity on the postbucking loads of nanoplates through considering various parameters such as nonlocal parameter, non-uniform parameter, mode number and compression ratio are examined and discussed. The results show that the nonlocal and nonuniform parameters have prominent effects on the postbuckling behavior of graphene sheets.

2. Mathematical modeling 2.1. Nonlocal nonlinear plate model In the following section, a nonlocal continuum model is developed for the nonlinear stability analysis of graphene sheets (GSs). A rectangular single-layered graphene sheet with uniform thickness (h) is shown in Fig. 1. Cartesian coordinate frame with axes x, y and z used for the graphene sheet is also shown in the figure. The length and width of the nanoplate are denoted by ‘x and ‘y , respectively. In the present formulation, the nonlocal stress resultants are defined as follows Z h=2 Z h=2 Z h=2 snl N yy ¼ snl Nxy ¼ tnl Nxx ¼ xx dz, yy dz, xy dz, h=2

Mxx ¼

Z

h=2 h=2

h=2

snl M yy ¼ xx zdz,

Z

h=2 h=2

h=2

snl M xy ¼ yy zdz,

Z

h=2 h=2

tnl xy zdz

ð1Þ

nl nl where snl xx and syy are the nonlocal normal stresses and txy is the nonlocal shear stress. It is well-known that unlike the large scale structures, mechanical behavior of nanostructures such as carbon nanotubes [16], nanorods [18] and graphene sheets [19] depends considerably on their sizes. The traditional local elasticity theory is a scale free theory and thus cannot predict the mechanical characteristics of nanomaterials properly. In recent years, many researchers have employed the nonlocal elasticity theory of Eringen [20,21] to overcome the shortcomings of the classical elasticity theory. This theory is based on a simple physical

Fig. 1. Discrete model of a rectangular single-layered graphene sheet (SLGS).

A. Farajpour et al. / Physica E 47 (2013) 197–206

concept that the components of stress tensor at a given point are a function not only of strain tensor at that point but also are a function of strain tensors at all other points in the domain. According to the nonlocal elasticity theory, the basic stress– strain equation for a Hookean solid neglecting the body force is expressed by the following relationship ZZZ   snl K 9xx0 9, l scl ð2Þ ij ¼ ij dV: V

  cl are the nonlocal and The parameters snl ij and sij ¼ C ijkl ekl classical stress tensors, respectively. Also, Cijkl and ekl are the fourth-order elasticity and strain tensors, respectively. The integration extends over the whole volume of the nanostructure (V). The term K is the nonlocal modulus or attenuation function which has the dimension of (length)  3. 9x x0 9 denotes the distance between x and x0 in the Euclidean form. l( ¼e0lc/le) is a constant that depends on the characteristic length ratio lc/le, where lc is an internal characteristic length (e.g., lattice parameter, granular distance, distance between C–C bonds) and le is an external characteristic length (e.g., crack length, wave length). The value of parameter e0 is crucial to calibrate the nonlocal model with experimental results or molecular dynamics (MD) simulation results. It is difficult to apply Eq. (2) for solving the nonlocal elasticity problems. Therefore, the following differential relation proposed by Eringen [21] is often used

snl ðe0 lc Þ2 r2 snl ¼ scl

where E1 , 1n12 n21

E12 ¼ E21 ¼

E22 ¼

n12 E2 , 1n12 n21

and twisting curvature kxy can be written as

kxx ¼ 

@2 w , @x2

kyy ¼ 

@2 w , @y2

@2 w @x@y

kxy ¼ 2

ð8Þ

Eliminating the in-plane displacement components u and v from Eq. (7), the following compatibility equation is obtained for the deformation of middle surface !2 2 2 @2 e0xx @ e0yy @ g0xy @2 w @2 w @2 w ¼ þ   2 ð9Þ @x@y @x@y @x @y2 @y2 @x2 The strain-displacement relations (kinematic equations) are independent of constitutive equations and therefore can be used to derive the nonlocal governing equations. In view of Eqs. (1), (4)–(6), the nonlocal stress resultants can be written in terms of mid-plane strains and curvatures as follows 9 9 2 8 8 38 0 9 N > N > 0 > > > B B12 > > > > exx > = = 6 11 < xx > < xx > 7< e0 = Nyy ðe0 lc Þ2 r2 Nyy 0 7 B21 B22 yy ¼6 ð10aÞ 5> > > 4 > > > > > ; ; ; : N xy > : N xy > : g0 > B33 > 0 0 xy

8 M > > < xx Myy > > :M xy

9 > > =

8 M > > < xx 2 2 M yy ðe0 lc Þ r > > > >M ; : xy

9 > > =

2

D11 6 D21 ¼6 4 > > ; 0

D12 D22 0

38 > kxx 7< 0 7 kyy 5> D33 : kxy 0

9 > = > ; ð10bÞ

ð3Þ

here (e0lc)2 is the nonlocal parameter which incorporates the small scale effect into the formulation.r2 denotes the Laplacian operator and is given by r2(*) ¼q2(*)/qx2 þq2(*)/qy2 for a twodimensional (2D) space. Using Eq. (3), the nonlocal constitutive relationship for orthotropic nanoplates in Cartesian coordinates can be written in the following matrix form 8 nl 9 8 nl 9 2 9 38 s > s > e > > E11 E12 0 > > > = = < xx > < xx > < xx > = nl nl 7 syy ðe0 lc Þ2 r2 syy ¼ 6 0 5 eyy ð4Þ 4 E21 E22 > > > > > > > :g > ; ; ; : tnl > : tnl > 0 0 E33 xy xy xy

E11 ¼

199

E2 , 1n12 n21 E33 ¼ G12

ð5Þ

here E1 and E2 are the Young’s modulus in the x and y directions, respectively; G12 denotes the shear modulus and n12, n21 are Poisson’s ratios. The postbuckling problem is treated as a geometrically nonlinear problem in the context of von Karman’s assumptions. In this way and according to the classical plate theory (CPT), the strain-displacement relations at an arbitrary point of the nanoplate can be written as

exx ¼ e0xx þzkxx eyy ¼ e0yy þ zkyy gxy ¼ g0xy þ zkxy

ð6Þ 0 ij

In the above equation e and kij (i,j¼x,y) are the mid-plane strains and curvatures, respectively. The components of strain tensor in the middle plane of the plate are expressed as     @u 1 @w 2 @v 1 @w 2 @u @v @w @w þ þ þ e0xx ¼ þ , e0yy ¼ , g0xy ¼ @x 2 @x @y 2 @y @y @x @x @y ð7Þ where u, v and w represent the displacement of point (x,y,0) along x, y and z directions, respectively. The bending curvatures kxx, kyy

where B11 ¼

E1 h E2 h n12 E2 h , B22 ¼ , B12 ¼ B21 ¼ , ð1n12 n21 Þ ð1n12 n21 Þ ð1n12 n21 Þ

B33 ¼ G12 h

ð11aÞ 3

D11 ¼

3

E1 h E2 h , D22 ¼ , 12ð1n12 n21 Þ 12ð1n12 n21 Þ 3

D12 ¼ D21 ¼

3

n12 E2 h G12 h , D33 ¼ 12ð1n12 n21 Þ 12

ð11bÞ

where Bij(i,j ¼1,2) and Dij(i,j¼1,2) are called the extensional and bending stiffnesses of the graphene sheet, respectively. Also, B33 and D33 are the shear and torsional rigidities of GS, respectively. It should be noted that when the nonlocal parameter is set to zero, e0lc ¼0, the Eq. (10) reduces to that of the local plate theory. Using relation n12E2 ¼ n21E1 [9], the mid-plane strains can be determined from Eq. (10a) in the following form 9 8 0 9 2 38 exx > 0 N xx > > > v12 =E1 1=E1 > > > = 1 = < 0 > < 6 7 2 2 eyy ¼ 0 7 N yy v12 =E1 1=E2 1ðe0 lc Þ r 6 4 5 > h > > > > 1=G12 > ; ; : g0 > : N xy > 0 0 xy

ð12Þ Substituting Eq. (12) into Eq. (9) yields " !  1 @2 N 1 1 @2 Nyy n12 @2 Nxx @2 Nyy xx 2 2 1ðe0 lc Þ r þ  þ h E1 @y2 E2 @x2 E1 @x2 @y2 !2 # 1 @2 Nxy @2 w @2 w @2 w   2 ¼ G12 @x@y @x@y @x @y2

ð13Þ

Eq. (13) is the nonlocal compatibility equation for the nonlinear analysis of micro/nanoplates made of orthotropic materials. Using the principle of virtual work, one obtains @N xx @Nxy þ ¼ 0, @x @y

ð14aÞ

@N xy @N yy þ ¼ 0, @x @y

ð14bÞ

200

A. Farajpour et al. / Physica E 47 (2013) 197–206

  @2 M xy @2 M yy @2 M xx @ @w @w N xx þN xy þ þ2 þqþ 2 2 @x @x @y @x@y @x @y   @ @w @w Nyy þ Nxy ¼0 þ @y @y @x

ð14cÞ

where q is the transverse force per unit area. It should be noted that, unlike the small-deflection theory, the stress resultants (Nxx, Nyy and Nxy) in nanoplate with large deflections depend not only on in-plane external edge loading but also on the stretching of mid-plane due to its bending. The Airy’s stress function may be introduced as follows Nxx ¼

@2 C , @y2

N yy ¼

@2 C , @x2

Nxy ¼ 

@2 C @x@y

ð15Þ

in which C ¼ ch. It is obvious that the above expressions for the stress resultants satisfy exactly Eqs. (14a) and (14b). Substituting Eq. (15) into Eqs. (13) and (14c) and taking into account the relations Eqs. (10b) and (8), one can obtain the following nonlocal governing equations for the nonlinear stability of SLGSs under general in-plane loading " #     1 @4 c 1 n12 @4 c 1 @4 c 2 2 1ðe0 lc Þ r þ 2  þ E2 @x4 2G12 E1 @x2 @y2 E1 @y4 !2 @2 w @2 w @2 w  2 ð16Þ ¼ @x@y @x @y2 @4 w @4 w @4 w D11 4 þ 2ðD12 þ 2D33 Þ 2 2 þ D22 4 @x @x @y @y !   q @2 c @2 w @2 c @2 w @2 c @2 w 2 2 þ þ h 1ðe0 lc Þ r 2 ¼0 h @y2 @x2 @x@y @x@y @x2 @y2 ð17Þ It is observed that four governing differential Eqs. (13) and (14a–c) are reduced to two governing differential Eqs. (16) and (17) by introducing the stress function C for in-plane stress resultants. Eqs. (16) and (17) can be interpreted as the nonlocal compatibility and equilibrium equations, respectively. These partial differential equations in terms of transverse displacement and stress function are coupled and nonlinear. It should be noted that traditional local governing differential equations for the large deflection of plates can be obtained by setting the nonlocal parameter equal to zero in the above equations. 2.2. Modeling of MLGSs considering nonlinear vdW interactions Fig. 2 illustrates an orthotropic multi-layered graphene sheet (MLGS) with N layers. The interaction between each two layers of MLGS is described by van der Waals (vdW) interaction. Every graphene’s layer consists of several hexagons which each of them contains six carbon atoms (Fig. 1). The Lennard–Jones pair potential can be used to calculate the vdW force between two carbon atoms: "   6 #   ^ 12 s s^ PELJ d ¼ 4e^  ð18Þ d d where PELJ denotes the Lennard–Jones potential energy that depends upon the distance between interacting atoms (d). The ^ are constants which are determined by the depth of terms e^ and s potential and equilibrium distance, respectively. The vdW force can be obtained from the first derivative of the potential energy (18) with respect to distance d: "   7 #     24e^ s ^ 13 d  s^ F LJ d ¼  2  PELJ d ¼ ð19Þ s^ d d dd

Fig. 2. Continuum model of a N-layered graphene sheet.

In order to study the postbuckling behavior of MLGSs, a nonlinear model may be developed for the vdW interaction between any two layers of graphene sheet by expressing the third-order Taylor expansion of FLJ around the equilibrium position and integrating the results over the entire plate. In this way, the van der Waals force acting on the ith layer is obtained [41]. qðiÞ ¼ vdW

N X

N   X  3 cij wi wj þ eij wi wj

j¼1

ð20Þ

j¼1

where cij ¼ 

" # pffiffiffi!2  14  8 ^ ^ 4 3 24e^ 13p s 1 7p s 1   12  6 , 9lc 3 lc s^ 2 3 lc zi zj zi zj ð21aÞ

" # pffiffiffi!2  16  10 ^ 4 3 336e^ 65p s 1 s^ 1 eij ¼   14 2p  8 9lc 6 lc lc s^ 4 zi zj zi zj ð21bÞ in which zj ¼ zj =lc (j¼1,2,y,N). The internal characteristic length is taken to be the carbon–carbon (C–C) bond length (lc ¼0.142 nm). Also, the equilibrium distance between any two adjacent layers of MLGS is 0.34 nm [39]. It is assumed that all graphene’s layers are flat and have the same thickness, width, length and material properties. Using Eqs. (16), (17) and (20), the nonlocal nonlinear governing equations of rectangular MLGS can be expressed as " #     1 @4 c 1 n12 @4 c1 1 @4 c1 2 2 1 1ðe0 lc Þ r þ2  þ E2 @x4 2G12 E1 @x2 @y2 E1 @y4 !2 @2 w1 @2 w1 @2 w1 ¼  , ð22aÞ @x@y @x2 @y2 @4 w1 @4 w1 @4 w1 þ2ðD12 þ 2D33 Þ 2 2 þ D22 4 @x @x @y @y4 ( N N   X X 2 2  1ðe0 lc Þ r c1k ðw1 wk Þ þ e1k ðw1 wk Þ3

D11

k¼1

k¼1

@2 c1 @2 w1 @2 c1 @2 w1 @2 c1 @2 w1 þ 2 þh 2 2 @x@y @x@y @y @x @x2 @y2

!) ¼0

" #     1 @4 c 1 n12 @4 c2 1 @4 c2 2 2 þ 2  þ 1ðe0 lc Þ r2 E2 @x4 2G12 E1 @x2 @y2 E1 @y4 !2 @2 w2 @2 w2 @2 w2 ¼  , @x@y @x2 @y2

ð22bÞ

ð23aÞ

A. Farajpour et al. / Physica E 47 (2013) 197–206

 @4 w2 @4 w2 @4 w2  2 2 þ 2ðD12 þ 2D33 Þ 2 2 þ D22  1ðe0 lc Þ r @x4 @x @y @y4 ( N N X X @2 c2 @2 w2 c2k ðw2 wk Þ þ e2k ðw2 wk Þ3 þ h  @y2 @x2 k¼1 k¼1 !) @2 c2 @2 w2 @2 c2 @2 w2 ¼0 þ 2 @x@y @x@y @x2 @y2

here PN is the in-plane load in the x direction exerted on the middle surface of Nth layer (i.e., at z=_ ¼ N1).

D11

3. Solution procedure ð23bÞ

^ 

" #    1 @4 c 1 n12 @4 cN 1 @4 cN N þ2  þ E2 @x4 2G12 E1 @x2 @y2 E1 @y4 !2 @2 wN @2 wN @2 wN  , @x@y @x2 @y2 2

1ðe0 lc Þ r

¼

2

 @4 wN @4 wN @4 wN  2 þ2ðD12 þ 2D33 Þ 2 2 þ D22  1ðe0 lc Þ r2 4 4 @x @x @y @y ( N N X X @2 cN @2 wN cNk ðwN wk Þ þ eNk ðwN wk Þ3 þh  @y2 @x2 k¼1 k¼1 !) @2 cN @2 wN @2 cN @2 wN þ 2 ¼0 @x@y @x@y @x2 @y2

An approximate expression for the transverse deflection of each nanoplate’s layer with all edges simply supported and movable can be written in the following form [41]     mpx npy sin ð27Þ wi ðx,yÞ ¼ W i sin

‘x

ð24aÞ

cos ð24bÞ

PN P 1

  2mpx

‘x

þ

 2   E 1 m‘ y W i 2npy cos   2 2 ‘y 32ðn‘x Þ ½1 þ4p2 n2 e0 lc =‘y  ð28Þ

ð25Þ

where _ denotes the equilibrium interlayer space between two adjacent layers of GS. P1 is the magnitude of compressive load per unit length acting in the plane of first layer (i.e., at z ¼0). Also, b is the numerical loading factor. By changing the value of b, one obtains various cases of non-uniform in-plane loading as shown in Fig. 3. For example, if b¼0, then the case of uniform compressive force is obtained. By taking b¼1, the compressive load varies linearly along the z direction. Furthermore, when the load factor is equal to two (b¼2), the external load is parabolically distributed across the thickness of the nanoplate. The parameter a can be easily determined as follows

P 1 ðN1Þb

i  P1 h E2 ðn‘x W i Þ 1 þ aði1Þb y2 þ mi x2 þ  2  2 2h 32 m‘y ½1 þ 4p2 m2 e0 lc =‘x  2

ci ¼ 

As seen from the above relations, the nonlinear buckling of a graphene sheet made of N layers which are bonded by vdW interaction is described by a system of 2N coupled nonlinear differential equations. In this work an attempt is made to analyze the postbuckling behavior of MLGSs under non-uniform in-plane external loading. The in-plane compressive force acting at the edges of the nanoplate, is considered to be constant in the x and y directions but varies arbitrary in the z direction (see Fig. 3). Thus, the external load is distributed across the thickness of graphene according to the following law

a¼

‘y

where m and n represent the buckling half-wave number in the x and y directions, respectively. Wi is the transverse displacement of center of ith layer. Substituting Eq. (27) into the compatibility Eq. (16) leads to

D11

 b P x ðzÞ ¼ P1 ½1 þ a z=_ 

201

ð26Þ

here (mi ¼Py/Px) is the compression ratio of ith layer. From the above relation, it can be easily seen that the nonlocal parameter (e0lc) appears in the coefficients of particular solution of Airy’s stress function. Furthermore, the in-plane force resultants can be obtained by substituting Eq. (28) into Eq. (15), i.e.,   E1 hðmpW i Þ2 2npy cos ,   2 2 ‘y 8‘x ½1 þ 4p2 n2 e0 lc =‘y  h i E2 hðnpW i Þ2 b h NðiÞ  2 i yy ¼ mi P 1 1 þ aði1Þ  2 8‘y 1 þ 4p2 m2 e0 lc =‘x   2mpx ð29Þ ,NðiÞ cos xy ¼ 0

h i b NðiÞ xx ¼ P 1 1 þaði1Þ 

‘x

It should be noted that prior to buckling, every nanoplate’s layer is loaded by a uniform external force in its plane. However, the force distribution over the edges of any layer will be nonuniform after buckling. Applying the general procedure of Galerkin’s method yields the following Z ‘ y Z ‘x 0

0

Gi ½fðx,yÞfðx,yÞdxdy ¼ 0

ð30Þ

in which the operator Gi is the right hand side of nonlinear equilibrium equation, i.e., " # @4 @4 @4 GðÞ ¼ W i D11 4 ðÞ þ 2ðD12 þ 2D33 Þ 2 2 ðÞ þ D22 4 ðÞ @x @x @y @y

Fig. 3. Schematic representation of a MLGS under different types of non-uniform in-plane load.

202

A. Farajpour et al. / Physica E 47 (2013) 197–206

" 

N X

#

" # N   X 2 2 3 cik ðW i W k Þ 1ðe0 lc Þ r ðÞ eik ðW i W k Þ

k¼1

k¼1

"     @2 2 2 2 2 ðÞ  1ðe0 lc Þ r ðÞ3 W i 1ðe0 lc Þ r NðiÞ xx ðx,yÞ @x2 # @2 @2 ðiÞ þ 2N ðiÞ ð  Þ þ N ð x,y Þ ð x,y Þ ð  Þ ð31Þ xy yy @x@y @y2 ð Þ is the The term  f(x,y)   basic function and is given by f x,y ¼ sin mpx=‘x sin npy=‘y . Substituting for nonlinear operator from Eq. (31) into Eq. (30), as well as using relations (29) and integrating over the entire domain of GS, an eigenvalue equation is obtained as follows   ½HP 1 ½F fW g ¼ f0g ð32Þ where 2

Hði,iÞ ¼

Hði,jÞ ¼

2 2

p m n   w2 n2  m2 þ mi w2 n2

½1 þ p  2  n 2

m þ 2w2 Q 12 þ w4 Q 22  n m " # 3p2 ð1n12 n21 Þ m4 w4 Q 22 n4 2  þ W max þ  2 4 m þ mi w2 n2 1 þ4w2 ðnplÞ2 1 þ4ðmplÞ2 4p

2



1 m2 þ

mi w

2 n2

2 4cij þ9W max eij  , 4 2 m 2 þ i 2 n2 b

p

N  X



F ði,iÞ ¼ 1 þ aði1Þ , F ði,jÞ ¼ 0,

 2 4cik þ9W max eik ,

ð33Þ

here i,j ¼1,2,y,N and [H] and [F] are the N  N square matrix. Also, {W} indicates the column matrix of central displacements Wi. In the above relations, the non-dimensional parameters are defined as follows 4

2 4

2 cij ‘x eij h ‘x P1 ‘x e0 lc ‘x , cij ¼ , eij ¼ , l¼ , w¼ , D11 D11 D11 ‘x ‘y W D12 þ2D33 2G12 W i ¼ i , Q 12 ¼ ¼ n21 þ ð1n12 n21 Þ, h D11 E1 D22 E2 Q 22 ¼ ¼ D11 E1

P1 ¼



n

P2

ð36Þ

 SLGS

, g1 and g2 are defined by 2

SLGS

4c þ9W e ; gi ¼ 212 2 max2 12 4p m þ mi w n2

¼ P SLGS 9m ¼ m gi , i

i ¼ 1,2

ð37Þ

SLGS

It is observed that there are two sets of nonlocal solution for the nonlinear stability of DLGSs. In general, both two nondimensional buckling loads depend on vdW interactions. However, when two layers of DLGS are subjected to the same in-plane force, nonlinear buckling loads can be expressed as   p2 m4 þ 2w2 Q 12 m2 n2 þ w4 Q 22 n4 PDLGS ¼ P SLGS ¼    2 m2 þ mw2 n2 ½1 þ p2 l m2 þ w2 n2  " # 3p2 ð1n12 n21 Þ m4 w4 Q 22 n4 2  þ W max þ  2 4 m þ mw2 n2 1 þ4w2 ðnplÞ2 1 þ4ðmplÞ2 ð39Þ

ð34Þ

The Galerkin’s method transforms the stability problem into a standard eigenvalue problem. The nonlinear buckling loads of MLGS (P 1 ) are the eigenvalues of Eq. (32) that can be found using standard eigenvalue extraction techniques. Note that the matrix [F] reduces to identity matrix (Iij ¼ dij) in the case of uniformly distributed loading (a¼0). Explicit expressions can be obtained for the buckling load of single- and double-layered grapheme sheets. In the case of SLGS, the system of algebraic Eq. (32) is reduced to one equation with one unknown coefficient W1. From this equation, one can easily determine the following equation for non-dimensional postbuckling load

p2 m2 n2 P SLGS ¼   2 2 2 ½1 þ p l m þ w2 n2  m2 þ mw2 n2

 2  n 2

m þ2w2 Q 12 þ w4 Q 22  n m " # 3p2 ð1n12 n21 Þ m4 w4 Q 22 n4 2  þ  2 þ W max 4 m þ mw2 n2 1 þ 4w2 ðnplÞ2 1 þ 4ðmplÞ2

SLGS

,

SLGS

for ia j for ia j

 n in which the term P 1

where m1 and m2 are the compression ratios of lower and upper layers of DLGS, respectively. To determine the nontrivial solution of two homogeneous Eq. (36), its determinant must be equal to zero. Thus, the following relation is obtained for the nonlinear buckling loads of DLGSs n  n  n   1 ð1 þaÞ P 1 P 1 DLGS ¼ þ P2 SLGS SLGS 2ð1 þ aÞ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  n  n ð38Þ  P2 2 þ4ð1 þ aÞg1 g2 7 ½ð1 þaÞ P 1

k¼1

mw

SLGS

 n Pi

 2 l2 m2 þ



smaller than their local counterparts. Again, the classical solution for the postbuckling of local plate under in-plane compression can be obtained by setting the nonlocal parameter equal to zero (l ¼0) in Eq. (35). In other words, the scale coefficient (nonlocal parameter) transforms the traditional equations of classical continuum mechanics into the corresponding governing equations of nonlocal continuum mechanics. In a similar way, the eigenvalue equation of double-layered graphene sheets (DLGSs) can be written in the following matrix form 3 02  n  1 )  " # ( P g 1 0 0 7   B6 1 SLGS  1 C W1  P ¼ 5 @4 A 1 DLGS n W2 0 1 þa 0 g2 P2

and 2

4c12 þ 9W e12 PDLGS ¼ P SLGS  2  2 max2 2  2p m þ mw n   p2 m4 þ2w2 Q 12 m2 n2 þ w4 Q 22 n4 ¼   2 m2 þ mw2 n2 ½1 þ p2 l m2 þ w2 n2  " # 3p2 ð1n12 n21 Þ m4 w4 Q 22 n4 2  þ W max þ  2 4 m þ mw2 n2 1 þ 4w2 ðnplÞ2 1 þ 4ðmplÞ2 2



4c12 þ 9W max e12   2p2 m2 þ mw2 n2

ð40Þ

From Eqs. (39) and (40), it is found that the first set of nonlinear buckling loads is independent of vdW interactions. These postbuckling loads are exactly the same as those of single-layered graphene sheets. However, the second set of nonlinear buckling loads depends on the vdW interactions between the two layers of the plate.

ð35Þ

It should be noted that the first term on the right-hand side of Eq. (35) indicates the critical load for the linear stability of SLGSs. Consequently, it is obvious that the nonlinear buckling load of the nanoplate is greater than the linear buckling load. Further, it can be concluded that the nonlocal postbuckling loads are always

4. Results and discussion In order to validate the numerical results, the linear nonlocal buckling loads of a square SLGS under biaxial compression (m ¼1) obtained by Eq. (35) are compared with those of molecular dynamic (MD) simulations as reported by Ansari and Rouhi [25].

A. Farajpour et al. / Physica E 47 (2013) 197–206

The material properties of GS are taken as E ¼1 TPa and n ¼0.19 [25]. The thickness of nanoplate is assumed to be constant (h¼0.34 nm). Fig. 4 shows the critical buckling loads of SLGS versus the variation of nanoplate’s length for various values of scale coefficient (e0lc)2. The buckling results of MD simulation [25] are also shown in the figure. It is observed that the present results are in good agreement with MD results for (e0lc)2 E2 nm. Thus the present formulation and solution approach can be trusted. As another example, the small deflection of multi-layered graphene sheet (MLGS) under in-plane edge loading is considered. All layers of MLGS have the same geometrical and material properties as follows [9] E1 ¼ 1765 Gpa,

E2 ¼ 1588 Gpa,

n12 ¼ 0:3, n21 ¼ 0:27,

‘x ¼ ‘y ¼ 10:2 nm, h ¼ 0:34 nm In addition, it is assumed that all nanoplate’s layers are subjected to the same biaxial loading (m ¼1). The linear buckling solutions of MLGS are compared with the available results reported by Pradhan and Phadikar [36] in Fig. 5. As mentioned in the foregoing, they investigated the linear stability of MLGS under uniform biaxial compression by means of nonlocal elasticity theory. The linear vdW interaction coefficient between any two layers of GS is taken as ci,i þ 1 ¼  45 Gpa nm  1[9]. Note that in the work of Pradhan and Phadikar [36], the vdW interaction

203

between two adjacent layers of graphene sheet is only taken into consideration (cij ¼0 for 9i  j9 41). From Fig. 5, it is clearly seen that present results for buckling load ratio (nonlocal buckling load/local buckling load) exactly match with those of exact ones in all cases. In recent years, a large variation of mechanical properties has been reported for the graphene sheets. Furthermore, different researchers employed different values of thickness to study the mechanical behavior of GSs [40]. In addition, it is well-known that the material properties of nanostructures are size-dependent. Thus, in the present work, the Young’s modulus, Poisson’s ratio and effective thickness of both armchair and zigzag graphene sheets with different sizes are considered (see Table 1). These elastic properties and effective thicknesses were obtained by MD simulations as reported in the paper of Shen et al. [40]. Here, in order to calibrate the nonlinear nonlocal plate model, molecular dynamics (MD) simulator ‘‘LAMMPS’’ is employed to obtain postbuckling loads of armchair monolayer GSs. A VelocityVerlet algorithm with a 0.5 fm s time step is used to integrate the Newton’s equations of motion to calculate trajectories of particles. MD simulation is carried out based on the adaptive intermolecular reactive empirical bond order potential (AIREBO) [42]. AIREBO is a potential energy which is widely used to describe intermolecular interactions in condensed-phased hydrocarbon structures such as graphene. The in-plane biaxial load is applied in a quasi-static way by proportional increase in its value. In order to simulate simply supported boundary conditions, it is assumed that one layer of carbon atoms at the edges of GS have no motion in the z direction (see Fig. 6). However, these carbon atoms have rigid-body motion in the x–y plane. At the beginning of MD simulation, energy minimization has been carried out for the square armchair SLGS III (Table 1) to reach the equilibrium configuration. The postbuckling loads obtained by MD simulations are listed in Table 2 for different magnitudes of maximum lateral deflection. To make reasonable comparison, analytical results of both classical and nonlocal plate theory are also listed in the table. Various values of nonlocal parameter in the range of 0.25–1 nm are considered. This range is taken because the calibrated values of nonlocal parameter for vibration analysis of Table 1 Geometrical and material properties of six types of GSs [40].

Fig. 4. Comparison of present results with those of MD simulation [25] for isotropic SLGS under biaxial compression.

Fig. 5. Comparison of present results with those calculated by Pradhan and Phadikar [36] for linear buckling of MLGS under uniform edge loading.

Type of GS

‘x ðnmÞ ‘y ðnmÞ h (nm) E1 (GPa) E2 (GPa) G12 (GPa) n12

AGS I AGS II AGS III ZGS I ZGS II ZGS III

9.519 6.995 4.888 9.496 7.065 4.855

4.844 4.847 4.855 4.877 4.887 4.888

0.129 0.143 0.156 0.145 0.149 0.154

2434 2154 1949 2145 2067 1987

2473 2168 1962 2097 2054 1974

1039 923 846 938 913 857

Fig. 6. Buckling mode shape of monolayer AGS III.

0.197 0.202 0.201 0.223 0.204 0.205

204

A. Farajpour et al. / Physica E 47 (2013) 197–206

Table 2 Comparison of nonlinear buckling load (nN nm  1) obtained via Eq. (35) with those of MD simulation for orthotropic monolayer AGS III under biaxial compression (m ¼ 1). Loading step

Wmax (nm)

Classical continuum model

Nonlocal continuum model

MD simulation

e0 ‘c (nm)

700 800 900 1000 1080

0.2839 0.3669 0.3752 0.4162 0.4510

1.1841 1.6125 1.6613 1.9186 2.1580

0.25

0.5

0.75

0.85

1

1.0970 1.4850 1.5293 1.7624 1.9791

0.9026 1.2051 1.2396 1.4214 1.5904

0.7015 0.9228 0.9480 1.0809 1.2046

0.6307 0.8252 0.8475 0.9643 1.0730

0.5375 0.6983 0.7167 0.8133 0.9031

0.6785 0.7754 0.8723 0.9693 1.0468

Fig. 8. Change of non-dimensional buckling load with non-dimensional pffiffiffi lateral   deflection W max =h for different values of non-uniform parameter ðe0 lc ¼ 2 nmÞ.

Fig. 7. Change of load ratio (N ¼ m ¼ 1) with non-dimensional  lateral deflection pffiffiffi (Wmax/h) for (a) armchair sheets e0 lc ¼ 2 nm and (b) zigzag  pffiffiffi ghraphene graphene sheets e0 lc ¼ 2=2 nm .

armchair SLGS at the room temperature (300 K) have been obtained in the range 0.27–0.67 nm [40]. As seen from Table 2, the classical continuum model fails to obtain the nonlinear buckling loads of the GS, while, the nonlocal continuum model with appropriate value of nonlocal parameter ðe0 ‘c  0:85 nmÞ provides more accurate results. The variation of buckling load ratio with non-dimensional transverse deflection (Wmax/h) is plotted for three types of armchair and zigzag graphene sheets (AGSs and ZGSs) in (Fig. 7a and b), respectively. The results are calculated for the first and second mode numbers (m¼n ¼1,2). It is found that as the nondimensional lateral deflection increases from 0 to 3, the buckling load ratio increases slowly for AGS I and ZGS I (w E2). In contrast, the buckling load ratio reduces by increasing non-dimensional lateral deflection for AGS II, AGS III, ZGS II and ZGS III (1r w r1.5). It is also found that the nonlocal effect is more noticeable in the second buckling mode in comparison with the first buckling mode. This phenomenon is because of small wavelength effect for higher modes. At smaller wavelengths (higher mode

Fig. 9. Change of non-dimensional buckling load with nonlocal parameter (e0lc) for different values of non-uniform parameter ðW max ¼ hÞ:

numbers), the interaction between atoms increases and it causes an increase in the small scale effects. To illustrate the influence of non-uniform parameter (b) on the postbuckling load, a five-layered armchair graphene sheet III under biaxial compression (a¼ m ¼ 1) is considered. Figs. 8 and 9 depict the  variation of non-dimensional postbuckling load of first 2 layer, P 1 ¼ P 1 ‘x =D11 , with non-dimensional transverse displacement at the center of nanoplate and nonlocal parameter, respectively. Five different values of non-uniform parameter (b¼ 0,0.5,1,1.5,2) are taken into consideration. The value of

A. Farajpour et al. / Physica E 47 (2013) 197–206

205

load of GS is presented. The numerical results are calculated for thep five-layered armchair grapheme sheet III (see Table 1). A value ffiffiffi of 2 nm is taken for the small scale parameter (e0lc). The nanoplate is loaded by in-plane edge force which varies according to the law P ðiÞ x ¼ iP 1 . Four different values in the range of 0–1.5 are chosen for the compression ratio. From Fig. 11, it can be seen that the non-dimensional postbuckling load increases with increasing the non-dimensional transverse displacement (Wmax/h) from 0 to 3. Also, it can be seen that as the compression ratio increases from 0 to 1.5, the rate of increase in critical buckling load with respect to lateral deflection decreases. Furthermore, it is obvious that increase in the value of compression ratio leads to a decrease in non-dimensional postbuckling loads of MLGSs.

Fig. 10. Change of non-dimensional buckling load with nonlocal parameter (e0lc) for various mode numbers.

5. Conclusions In this paper, a nonlinear nonlocal plate model has been developed to study the postbuckling behavior of multi-layered graphene sheets (MLGSs) loaded by non-uniformly distributed inplane force across the thickness. The geometrical nonlinearity was taken into account with the use of von Karman’s assumptions. Nonlinear van der Waals (vdW) interactions between any two layers of a MLGS were modeled based on the Lennard–Jones pair potential. Numerical results have been presented for six types of orthotropic graphene sheet with all edges simply supported. From the results following conclusions are noticeable:

 The nonlocal parameter has a decreasing effect on the critical postbuckling loads of MLGSs.

 The critical postbuckling load of MLGSs decreases with increasing non-uniform parameter (b) from 0 to 2.

 Small scale effects decrease with increasing the non-dimensional lateral deflection at the plate centre for AGS I and ZGS I (w E2).

 Small scale effects increase with increasing the nonFig. 11. Change of non-dimensional buckling load with non-dimensional lateral deflection (Wmax/h) for different values of compression ratio.

nonlocal parameter is taken in the range of 0–2 nm. One logical reason for this is that Wang and Wang [43] reported that the scale factor (e0lc) of a single-walled carbon nanotube (SWCNT) must be smaller than 2.0 nm. In addition, Shen et al. [40] obtained small scale parameter for vibration analysis of various types of SLGS with the use of MD simulations. In their work, it is reported that the value of nonlocal parameter changes from 0.22 nm to 0.67 nm depending upon geometrical properties of GS. From (Figs. 8 and 9), it is clearly seen that critical buckling load decreases with increasing non-uniform parameter from 0 to 2. Furthermore, it is observed that nonlocal parameter has a decreasing effect on the critical postbuckling loads of MLGSs. In order to study the effect of small scale on higher buckling modes, the non-dimensional postbuckling load of five-layered AGS III versus the variation of nonlocal parameter is plotted for first four mode numbers (m ¼n¼ 1,2,3,4) in Fig. 10. The maximum value of lateral deflection of nanoplate is assumed to be W max ¼ h. The five-layered AGS III is subjected to biaxial compressive force (m ¼1) that increases linearly across the thickness of the plate (a¼b ¼1). It can be seen that nonlocal effects are more significant at higher buckling modes. It is also found that the gap between the curves gradually decreases with increasing nonlocal parameter from 0 nm to 2 nm. This means that the influence of mode number on the nonlinear stability of graphene sheets decreases with increasing nonlocal effects. In Fig. 11, the effect of compression ratio in conjunction with non-dimensional lateral deflection on the critical postbuckling

 

dimensional lateral deflection at the plate centre for AGS II, AGS III, ZGS II and ZGS III (1r w r1.5). The effect of small length scale is higher at higher buckling modes. The influence of mode number on the nonlinear stability of graphene sheets decreases with increasing nonlocal effects.

Acknowledgment The authors gratefully acknowledge the valuable comments made on an earlier version of the manuscript by an anonymous reviewer which significantly improved the paper. References [1] M. Li, H.X. Tang, M.L. Roukes, Nature Nanotechnology 2 (2007) 114. [2] C.Y. Tsai, S.P. Lu, J.W. Lin, P.T. Lee, Applied Physics Letters 98 (2011) 153108. [3] J.S. Bunch, A.M. van der Zande, S.S. Verbridge, I.W. Frank, D.M. Tanenbsum, J.M. Parpia, H.G. Craighead, P.L. McEuen, Science 315 (2007) 490. [4] B.I. Yakobson, C.J. Brabec, J. Bernholc, Physical Review Letters 76 (1996) 2511. [5] H.J. Hwang, K.R. Byun, J.W. Kang, Physica E: Low-dimensional Systems and Nanostructures 23 (2004) 208. [6] H.J. Hwang, J.W. Kang, Physica E: Low-dimensional Systems and Nanostructures 27 (2005) 163. [7] O. Leenaerts, B. Partoens, F.M. Peeters, Applied Physics Letters 93 (2008) 193107. [8] W.H. Duan, C.M. Wang, Nanotechnology 20 (2009) 075702. [9] K. Behfar, R. Naghdabadi, Composites Science and Technology 65 (2005) 1159. [10] X.Q. He, S. Kitipornchai, K.M. Liew, Nanotechnology 16 (2005) 2086. [11] L. Tapaszto´, T. Dumitrica, S.J. Kim, P. Nemes-Incze, C. Hwang, L.P. Biro´, Nature Physics (2012), http://dx.doi.org/10.1038/NPHYS2389. [12] D.B. Zhang, E. Akatyeva, T. Dumitrica, Physical Review Letters 106 (2011) 255503.

206

A. Farajpour et al. / Physica E 47 (2013) 197–206

[13] J. Peddieson, R. Buchanan, R.P. McNitt, International Journal of Engineering Science 41 (2003) 305. [14] L.J. Sudak, Journal of Applied Physics 94 (2003) 7281. [15] J.N. Reddy, S.D. Pang, Journal of Applied Physics 103 (2008) 023511. [16] H. Heireche, A. Tounsi, A. Benzair, M. Maachou, E.A.Adda Bedia, Physica E: Low-dimensional Systems and Nanostructures 40 (2008) 2791. [17] T. Murmu, S.C. Pradhan, Mechanics Research Communications 36 (2009) 933. [18] M. Danesh, A. Farajpour, M. Mohammadi, Mechanics Research Communications 39 (2012) 23. [19] T. Aksencer, M. Aydogdu, Physica E: Low-dimensional Systems and Nanostructures 43 (2011) 954. [20] A.C. Eringen, D.G.B. Edelen, International Journal of Engineering Science 10 (1972) 233. [21] A.C. Eringen, Journal of Applied Physics 54 (1983) 4703. [22] Y. Chen, J.D. Lee, A. Eskandarian, International Journal of Solids and Structures 41 (2004) 2085. [23] W.H. Duan, C.M. Wang, Y.Y. Zhang, Journal of Applied Physics 101 (2007) 024305. [24] R. Ansari, S. Sahmani, B. Arash, Physics Letters A 375 (2010) 53. [25] R. Ansari, H. Rouhi, Solid State Communications 152 (2012) 56. [26] A. Sakhaee-Pour, M.T. Ahmadian, A. Vafai, Solid State Communications 145 (2008) 168. [27] W.H. Duan, C.M. Wang, Nanotechnology 18 (2007) 385704. [28] A. Hemmasizadeh, M. Mahzoon, E. Hadi, R. Khandan, Thin Solid Films 516 (2008) 7636.

[29] R. Aghababaei, J.N. Reddy, Journal of Sound and Vibration 326 (2009) 277. [30] T. Murmu, S.C. Pradhan, Physica E: Low-dimensional Systems and Nanostructures 41 (2009) 1628. [31] P. Malekzadeh, A.R. Setoodeh, A.Alibeygi Beni, Composite Structures 93 (2011) 2083. [32] A. Farajpour, M. Mohammadi, A.R. Shahidi, M. Mahzoon, Physica E: Lowdimensional Systems and Nanostructures 43 (2011) 1820. [33] A. Farajpour, A.R. Shahidi, M. Mohammadi, M. Mahzoon, Composite Structures 94 (2012) 1605. [34] T. Murmu, S. Adhikari, Composites Part B Engineering 42 (7) (2011) 1901. [35] H. Liu, J.L. Yang, Physica E: Low-dimensional Systems and Nanostructures 44 (2012) 1236. [36] S.C. Pradhan, J.K. Phadikar, Journal of Computational and Theoretical Nanoscience 7 (2010) 1948. [37] S. Adali, Journal of Theoretical and Applied Mechanics 49 (3) (2011) 621. [38] H. Babaei, A.R. Shahidi, Acta Mechanica Sinica 27 (6) (2011) 967. [39] J. Wang, X. He, S. Kitipornchai, H. Zhang, Journal of Physics D 44 (2011) 135401. [40] L. Shen, H.S. Shen, C.L. Zhang, Computational Materials Science 48 (2010) 680. [41] E. Jomehzadeh, A.R. Saidi, N.M. Pugno, Physica E: Low-dimensional Systems and Nanostructures (2012), http://dx.doi.org/10.1016/j.physe.2012.05.015. [42] S.J. Stuart, A.B. Tutein, J.A. Harrison, Journal of Chemical Physics 112 (2000) 6472. [43] Q. Wang, C.M. Wang, Nanotechnology 18 (2007) 075702.