Application of modified couple-stress theory to stability and free vibration analysis of single and multi-layered graphene sheets

Journal Pre-proof Application of modified couple-stress theory to stability and free vibration analysis of single and multi-layered graphene sheets

Zahra Shafiei, Saeid Sarrami-Foroushani, Fatemeh Azhari, Mojtaba Azhari

PII:

S1270-9638(19)31133-2

DOI:

https://doi.org/10.1016/j.ast.2019.105652

Reference:

AESCTE 105652

To appear in:

Aerospace Science and Technology

Received date:

27 April 2019

Revised date:

9 October 2019

Accepted date:

15 December 2019

Please cite this article as: Z. Shafiei et al., Application of modified couple-stress theory to stability and free vibration analysis of single and multi-layered graphene sheets, Aerosp. Sci. Technol. (2019), 0, 105652, doi: https://doi.org/10.1016/j.ast.2019.105652.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier.

Application of modified couple-stress theory to stability and free vibration analysis of single and multi-layered graphene sheets

Zahra Shafiei1, Saeid Sarrami-Foroushani1, Fatemeh Azhari2, Mojtaba Azhari1*

1

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

2

Department of Civil Engineering, Monash University, Melbourne, VIC3800, Australia

Submitted to:

Aerospace Science and Technology

* Corresponding author: Department of Civil Engineering, Isfahan University of Technology, Isfahan, 84156-83111, Iran. Tel.: +98 3133913804. Fax: +98 3133912700. E-mail address: [email protected] (M. Azhari)

Abstract With the development of nanotechnology, research activities on carbon nanostructures have increased rapidly. In recent years, due to the extraordinary mechanical properties of graphene sheets, there has been a growing interest in investigating the mechanical response of these carbon nanostructures. In this article, the modified couple-stress theory (MCST) is first employed to study the free vibration and mechanical buckling of single-layered graphene sheets (SLGSs). To this end, SLGS is modeled as a nanoplate and the two-variable refined plate theory (TVRPT) is adopted to extend the finite strip method (FSM) formulation. The natural free vibration frequency and mechanical buckling loads of the sheet are then obtained by solving the proper eigenvalue problems. Mechanical buckling and free vibration of multilayered graphene sheets (MLGSs) are also investigated considering the effects of van der Waals (vdW) bonds between the layers. Modified couple-stress theory is applied to consider the small-scale effects of the graphene sheets. The results obtained by the proposed method are validated against those available in the literature. Finally, a comprehensive parametric study is performed to investigate the effects of different parameters such as loading schemes, nanoplate dimensions and boundary conditions.

Keywords: Graphene sheets; Modified couple-stress theory; Refined plate theory; Mechanical buckling; Free vibration; Finite strip method.

2

1. Introduction Graphene is a carbon material with nanoscale structure which has superior mechanical, chemical, electrical, and thermal properties. The abundant application of single and multilayered graphene sheets (SLGSs and MLGSs, respectively) necessitates a deep understanding of the behavior of these nanostructures. Experimental investigations on nanostructures are limited because of their small size and the expensive nature of conducting such investigations through atomic simulations [1-7]. These limitations have led to development of different modeling approaches for studying the behavior of these nanostructures based on the mechanics of continuous media [8-17]. In some of these proposed models, such as the classical elasticitybased models, the SLGS is modelled as a plate, and the governing differential equations are solved to evaluate the response of the sheet. However, since the classical continuum mechanics does not consider the effects of intermolecular forces, a discrepancy is found between the results of these models and the experimental results. To solve this problem, several theories have been developed to take into account the size effects. In Ref. [8], Cosserat brothers developed the theory of deformable bodies, in which a classical continuous environment consisting of a series of substructures with the ability to transform independently was considered. In Refs. [12, 13], the couple-stress theory was proposed assuming that the rotation of substructures are the same as their surroundings. In this theory, small-scale effects are taken into account by introducing higher-order gradients in the strain energy. To this end, two length scale parameters (LSP) are applied to higher order gradients, which are usually difficult to find. In Ref. [15], the modified couple-stress theory (MCST) which assumes only one LSP was proposed by Yang et al. Subsequently, several robust models were proposed based on the MCST, such as the Bernoulli-Euler beam model [18, 19], the Timoshenko beam model [20, 21], the Reddy-Levinson beam model [22], the doubly-curved shell model [23], the Kirchhoff plate model [24], and the Mindlin plate model [25]. Using these models, many studies have

3

been carried out to investigate the bending, buckling and free vibration of different microbeams [18-22] and microplates [23-31]. In case of modelling the plate behavior under different loading conditions, the classical plate theory (CPT) and the shear deformation theories are the most widely used theories. In Ref [32], Shimpi combined these two theories and proposed the two-variable refined plate theory (TVRPT), which has been employed in several studies [3343]. Huu-Tai Thai and Seung-Eock Kim [41-43] presented a closed form solution by using TVRPT in Levy-type form to analyzs free vibration and buckling of orthotropic plates. Narendar and Gopalakrishnan [33, 34] used TVRPT and considered the nonlocal effects to examine the nanoscale buckling of both isotropic and anisotropic SLGSs. They considered biaxial loading, simply supported boundary conditions and the Navier's method in their study. Satish et al. [35] analyzed the thermal vibration of orthotropic nanoplates using TVRPT, based on the nonlocal continuum model. Huu-Tai Thai and Seung-Eock Kim [44, 45] used TVRPT to examine free vibration, bending and buckling of laminate composite plates. The TVRPT was rewritten by introducing displacements in terms of trigonometric functions, and the bending, buckling, and vibration of FGM plates were studied [36-38]. Sarrami-Foroushani and Azhari [39] studied the nonlocal vibration and buckling of thick rectangular nanoplates using the finite strip method (FSM) based on TVRPT. In Ref [40], Sobhy used the TVRPT to study the free vibration, mechanical buckling and thermal buckling of MLGSs. Hygrothermal vibration of orthotropic double-layered graphene sheets embedded in an elastic medium was studied by Sobhy [46] considering the two-variable plate theory. In this study, mechanical buckling and free vibration of SLGSs and MLGSs are investigated based on the modified couple-stress theory (MCST) developed by Yang et al. [15]. TVRPT is used to consider the shear deformation effects, and finite strip formulation is developed and employed to study nanostructures. Numerical results are presented by evaluating the influence of length scale parameters, plate dimensions, higher modes, different

4

boundary conditions, different load patterns and their interactions on the behavior of nanostructures. This paper is organized into the following sections. In Section 2, derivation of the governing equations based on the MCST and TVRPT are presented. In Section 3, mathematical formulation of the finite strip method is presented and the method is developed to study the buckling and vibration of graphene sheets. Numerical results and discussions are demonstrated in Section 4 for two types of SLGSs and MLGSs. Finally, concluding remarks are presented in Section 5.

2. The governing equations In this section, derivation of the governing equations using the modified couple-stress theory (MCST) and the two-variable refined plate theory (TVRPT) is explained. 2.1. Modified Couple-Stress Theory (MCST) According to the MCST developed for the linear isotropic materials, density of the strain energy, U , is expressed as [15]:

U



1 2 O >trε @  P εTε  l 2 χ T χ 2



(1)

Where ε , χ are known as strain tensor and symmetric part of curvature tensor, respectively, and tr shows the sum of the elements on the main diagonal of the strain tensor. Also, l is the length scale parameter (LSP), and O and P are the Lame's constants. Considering the higher order gradients, the strain energy, U, stored in a deformed elastic body is defined as:

U

1 V ij H ij  mij Fij dv 2 ³v





(2)

in which V ij and H ij represent the components of the stress and strain tensors, respectively.

5

The symmetric part of the couple-stress and curvature tensors are also shown by mij and Fij , respectively. The relationships between V ij and H ij , as well as mij and Fij are written as:

V ij

OH kkGij  2PH ij

(3a)

mij

2l 2 PFij

(3b)

where G ij is the Kroncker delta function. The components of strain and curvature tensors in terms of the displacement field are expressed as:

H ij

1 (ui , j  u j ,i ) 2

1 wui wu j ( )  2 wx j wxi

(4a)

F ij

1 (Ti , j  T j ,i ) 2

1 wTi wT j ( )  2 wx j wxi

(4b)

in which ui and Ti represent the displacement and rotation vectors, respectively. The rotation vector is also defined as:

Ti

1 eijk uk , j 2

wu 1 eijk k wx j 2

(5)

in which eijk is the anti-symmetric permutation symbol. 2.2 Two-Variable Refined Plate Theory (TVRPT) Based on the TVRPT introduced by Shimpi [32], the displacement field of a plate in the Cartesian coordinate system (x, y, z), in which the z axis is perpendicular to the surface of the plate, is assumed as: u( x, y, z) u ( x, y)  zwb, x  F ( z)ws, x

(6a)

v( x, y, z) v ( x, y)  zwb, y  F ( z)ws, y

(6b)

w( x, y, z ) wb ( x, y )  ws ( x, y )

(6c)

where u, v and w, are displacements in the x, y, and z directions respectively. u and v are the mid-plane displacements which are ignored in this study as the studied graphene sheets are

6

isotropic and homogeneous. wb and ws are the bending and shear components of the lateral displacement w, respectively. Therefore, according to the assumed displacement field, the unknown variables to be found are wb and ws . In Section 3.1, these variables will be defined in terms of x and y coordinates as well as a set of new unknown variables using the finite strip method. In this theory, F(z) is defined in a way that the transverse shear stresses, W xz and W yz , take parabolic variations through the plate thickness, h, and are zero at the free surfaces of the plate (z=-h/2 and z=h/2). Thus, it is assumed that: F ( z)

(7)

ª1 z 5 z º h « ( )  ( )3 » ¬4 h 3 h ¼

Using Eqs. (4) and (6), the non-zero components of strain and curvature tensors are obtained as:

Hx

2 w 2 wb ª 1 z 5 z 3 º w ws z 2  h « ( )  ( ) » 2 wx ¬ 4 h 3 h ¼ wx

(8a)

Hy

z

2 w 2 wb ª 1 z 5 z 3 º w ws  h  ( ) ( ) «¬ 4 h 3 h »¼ wy 2 wy 2

(8b)

J xy

2H xy

2 z

J yz

2H yz

z 2 º wws ª5 «¬ 4  5( h ) »¼ wy

(8d)

J xz

2H xz

z 2 º wws ª5 «¬ 4  5( h ) »¼ wy

(8e)

2 w 2 wb ª 1 z 5 z º w ws  2h « ( )  ( )3 » wywx ¬ 4 h 3 h ¼ wywx

(8c)

and

F xx

w 2 wb 3 5 z 2 w 2 ws (  ( ) ) wxwy 8 2 h wxwy

(9a)

7

w 2 wb 3 5 z 2 w 2 ws (  ( ) ) wxwy 8 2 h wxwy

F yy



F xy

º w2 w2 w2 wF ( z ) 1 ª w2 ws ) » «( 2  2 )( wb  ws )  ( 2  2 )( wb  wz 4 ¬ wy wx wy wx ¼

(9c)

F xz

1 z wws (10 2 ) 4 h wy

(9d)

F xz

z wws 1 (10 2 ) 4 h wx

(9e)

(9b)

By introducing Eqs. (8) and (9) into Eqs. (3a) and (3b) derived using the MCST, the nonzero components of the stress and couple-stress tensors are obtained as:

Vx



w 2 wb w 2 ws Ez w 2 wb E 1 z 5 z 3 w 2 ws   h   [ ] [ ( ) ( ) ][ ] Q Q wy 2 wx 2 wy 2 1 Q 2 wx 2 1 Q 2 4 h 3 h

(10a)

Vx



w 2 wb w 2 ws Ez w 2 wb E 1 z 5 z 3 w 2 ws   h   [ ] [ ( ) ( ) ][ ] Q Q wx 2 wy 2 wx 2 1 Q 2 wy 2 1 Q 2 4 h 3 h

(10b)

W xy



2 Ez w 2 wb ª 1 z 5 z º w ws  h « ( )  ( )3 » ( ) 1 Q wxwy ¬ 4 h 3 h ¼ wxwy

(10c)

W yz

E ª5 z º ww  5( ) 2 » s « h ¼ wy 2(1 Q ) ¬ 4

(10d)

W xz

E ª5 z º ww  5( ) 2 » s « h ¼ wx 2(1 Q ) ¬ 4

(10e)

and

w 2 wb 3 5 z 2 w 2 ws (  ( ) ) ) wxwy 8 2 h wxwy

mxx

2Pl 2 (

m yy

2Pl 2 (

mxy

º w2 w2 w2 wF ( z ) 1 2 ª w2 P l «( 2  2 )( wb  ws )  ( 2  2 )( wb  ws ) » wz 2 wy wx ¬ wy wx ¼

(11c)

mxz

1 2 z ww P l (10 2 s ) 2 h wy

(11d)

(11a)

w 2 wb 3 5 z 2 w 2 ws (  ( ) ) ) wxwy 8 2 h wxwy

(11b)

8

m yz

z ww 1 2 P l (10 2 s ) 2 h wx

(11e)

Therefore, the system's strain energy can be extracted according to the TVRPT by substituting Eqs. (8) to (11) in Eq. (2). Then, by determining the potential and kinetic energy, the governing equations of motion can be derived based on the Hamilton's principle as:

G3

t1

³t

G (U  V p  T ) dt

0

(12)

In this function, U is the strain energy defined by Eq. (2), V p is the strain energy of the external in-plane forces acting on the edges of the nanoplate, which is defined as: Vp

1 ª N x ( w, x )2  N y ( w, y )2  2 N xy ( w, x )( w, y ) ºdA ¼ 2³¬

(13)

where Nx, Ny and Nxy are the in-plane forces per unit length.

T is the kinetic energy corresponding to the free vibration of plates, which can be written in terms of the time of vibration (t) as [47]: T

1 m ª(u ,t )2  (v ,t )2  ( w,t )2 º¼ dA 2 ³A ¬

(14)

in which m is the mass of the plate per unit area, u and v are the mid-plane displacements which are ignored in this section, and w wb  ws .

3. The solution method In this study, for the first time, TVRPT and MCST are incorporated into the Finite Strip Method to investigate the buckling and free vibration behaviors of rectangular SLGSs and MLGSs with different boundary conditions. The procedure of this implementation is provided in details in the following sections. 3.1. Discretization Fig. 1 shows a single strip of length L and width b in the rectangular coordinate system (x,

9

y, z) with two nodal lines of i and j. The strip nodal degrees of freedom are also shown in this figure.

Fig. 1. The strip and degrees of freedom

Considering the displacement field defined in Eq. (6) based on the TVRPT, the bending and shear deflections of the strip are respectively defined as: r

wb ( x, y )

¦ X nb ( x).Ynb ( y )

(15a)

n 1 r

ws ( x, y )

¦ X ns ( x).Yns ( y)

(15b)

n 1

in which the subscripts and superscripts b and s represent the deformations due to the bending and shear, respectively, and r is the number of harmonic modes. X n ( x) is the suitable polynomial shape function in transverse direction and Yn ( y ) is the appropriate trigonometric shape function in longitudinal direction, which are chosen such that the boundary conditions are satisfied. Table 1 presents the trigonometric shape functions used in this study for different boundary conditions in the longitudinal direction (y direction). In this table, S, C and F respectively represent the simply supported, clamped and free boundary conditions at two ends of the plate.

10

Table 1. The trigonometric shape functions used in finite strip method Boundary

Trigonometric shape function ( Yn ( y ) )

conditions SS

sin

mS y L

CC

sin

mS y Sy sin L L

SC

sin

(m  1)S y m  1 mS y ( )sin L m L

CF

1  cos

FF

Y1 ( y ) 1

(m  1)S y 2L

Y2 ( y) 1 

P1 0

Ym ( y) sin(

Dm

P2 1

2y L

Pm y L

)  sinh(

Pm y L

)  D m[cos(

sin Pm  sinh Pm , Pm cos Pm  cosh Pm

2m  3 S 2

Pm y L m

)  cosh(

Pm y L

)]

3, 4,...

Moreover, using the Hermitian shape functions as polynomial shape functions in the transverse direction ( X n ( x) ), Eqs. (15a) and (15b) may be rewritten as:

ª 3x 2 2 x3 2 x 2 x3 b b       2 )(Ti ) n  (1 )( w ) ( x « ¦ i n b b2 b3 b n 1¬ r

wb ( x, y )

º 3x 2 2 x3 x 2 x3 ( 2  3 )( wib ) n  (  2 )(Tib ) n » Ynb (y) b b b b ¼

11

(16a)

ª 3x 2 2 x3 2 x 2 x3 s s       2 )(Ti ) n  (1 )( w ) ( x « ¦ i n b b2 b3 b n 1¬ r

ws ( x, y )

º 3x 2 2 x3 x 2 x3 ( 2  3 )( wis ) n  (  2 )(Tis ) n » Yns (y) b b b b ¼

(16b)

where ( wi )n , (Ti )n , (w j )n and (T j )n are the deflections and rotations of each nodal line corresponding to the nth mode. Eqs. (16a) and (16b) can be rewritten in the vector form as: r

wb ( x, y )

¦ Lbnδbn

(17a)

n 1 r

ws ( x, y )

¦ Lsnδns

(17b)

n 1

in which

Lbn

Lsn

Ln

ª 3x 2 2 x3 2 x 2 x3 3 x 2 2 x3 x 2 x3 º      ),(  3 ), (  2 ) » Yn (1 ), ( x « b b2 b2 b b ¼ b2 b3 b ¬

(18)

δbn and δ ns , which are the displacement vectors related to the nth mode, are also given by

T

δbn

­ b ww b b ww b ½ ® wi , ( )i , w j , ( ) j ¾ wx wx ¿n ¯

T

δ ns

­ s ww s s ww s ½ ® wi , ( )i , w j , ( ) j ¾ wx wx ¿n ¯

Where Ti

(w, x )i and T j

­ winb ½ ° b° °° Tin °° ® b ¾ ° w jn ° ° b ° °¯T jn °¿

b °­ δin °½ ® b ¾ °¯δ jn °¿

(19a)

­ wins ½ ° s ° °° Tin °° ® s ¾ ° w jn ° ° s ° ¯°T jn ¿°

s °­ δin °½ ® s ¾ °¯δ jn °¿

(19b)

(w, x ) j The displacement vector including all degrees of freedom

for one strip can be finally written as:

12

­δbi ½ ° b° °°δ j °° ® s¾ ° δi ° ° s° ¯°δ j ¿°

δn

r

¦ ^wib

Tib

wbj T bj

wis Tis

wsj T js

n 1

`

T n

(20)

3.2. Mechanical Buckling of Single-Layered Nanoplates In this section, finite strip relations for the mechanical buckling of SLGSs are derived. To this end, Eq. (17) is used to rewrite the displacement field proposed in Eq. (6) as:

­u ½ ° ° ® v ¾ N nδn ° w° ¯ ¿

(21)

in which N n is the shape function matrix and is defined as:

Nn

ª  z (Ln, x ) F ( z )(Ln, x ) º « » «  z (Ln, y ) F ( z )(L n, y ) » « L » Ln n ¬ ¼

(22)

By inserting Eq. (21) into Eq. (4), the strain and curvature matrices are derived as:

ε B nδ n

χ

(23a)

JBcnδ n

(23b)

in which Bn, J and B'n matrices are defined as:

Bn

F(z)(L n, xx ) º ª  z (L n, xx ) «  z (L ) F(z)(L n, yy ) »» n , yy « « 2 z (L n, xy ) 2 F(z)(L n, xy ) » « » 0 (F(z),z  1)(L n, y ) » « « » 0 (F(z),z  1)(L n, x ) ¼» ¬«

(24a)

13

J

Bcn

ª2 «0 1« «0 4« «0 «¬ 0

0 0 0 0º 2 0 0 0 »» 0 1 0 0» » 0 0 1 0» 0 0 0 1 »¼

(24b)

2(L n, xy ) (1  F(z), z )(L n, xy ) ª º « » (F(z), z  1)(L n, xy ) « 2(L n, xy ) » « 2(L » n , yy  L n , xx ) (1  F(z), z )(L n , yy  L n , xx ) » « « » 0 (F(z), zz )(L n, x ) « » 0 ( F(z), zz )(L n, y ) «¬ »¼

(24c)

According to Eq. (3) and using Eq. (23), σ and m can be written as:

σ Dε DB nδ n m

2Pl 2 χ

(25a)

2Pl 2 JBcnδn

(25b)

where D is the stiffness matrix of the material. In case of isotropic materials, D is defined as:

ª1 «Q E « «0 D 1 Q 2 « «0 «¬ 0

Q 1 0 0 0

º » » » (1 Q ) / 2 » 0 (1 Q ) / 2 0 » 0 0 (1 Q ) / 2 »¼ 0 0

0 0 0

0 0 0

(26)

in which E and Q are the modulus of elasticity and Poisson's ratio, respectively. Considering Eqs. (23) and (25), the strain energy defined in Eq. (2) may be expressed as: U

1 (ε Tσ  χ Tm)dv ³ 2 1 (δTn BTn DB mδm  δTn BcnT J T 2Pl 2 JBcmδ m )dv 2³

1 T ªδ (K  K c)δ º ¼ 2¬

(27)

Thus, according to Eq. (27), the total stiffness matrix of a single strip, Kt ,strip , can be defined as: K t ,strip

K  Kc

³ (Bn DBm ) dv  ³ (Bcn J T

T T

2Pl 2 JBcm ) dv

(28)

By subtituting Eq. (17) into the potential energy defined by Eq. (13), the following strip

14

geometric stiffness matrix, KG,strip , is derived as: K G ,strip

³  BGn NBGn dA T

(29)

in which BGn and N are defined as: BGn

N

ª L n, x «L ¬ n, y ª Nx « ¬« N xy

L n, x º L n, y »¼

(30a)

N xy º » N y ¼»

(30b)

The stiffness and geometric matrices derived for each strip correspond to the degrees of freedom shown in Fig. 1. These degrees of freedom lie on the nodal lines, which are the common borders of the strips (Fig. 2). Therefore, the compatibility equations must be satisfied along the nodal lines for these degrees of freedom. Z,w

X

Nodal lines

j i

N

L

Y

Fig. 2. Nodal line in the border of the strips

During this process, the total stiffness corresponding to each degree of freedom is determind by sumation of the corresponding stiffness components of each strip. Thus, the stiffness and geometric matrices of the entire nanoplate, called K tt and K Gt , respectively, are determined and as a result, the strain and potential energy of plate are obtained. By inserting these energies in Eq. (12) , the following eigenvalue equation is obtained in the absence of the kinetic energy:

15

K

tt



 K Gt δ

0

(31)

Eq. (31) is a typical eigenvalue problem which is solved to find the critical buckling load by vanishing the determinant of K tt  K Gt . In this case, δ is a vector in the null-space of

K tt  K Gt matrix, representing the eigenmodes or the buckling modes of the nanoplate. 3.3. Vibration of Single-Layered Nanoplates In order to consider the effect of time (t) and the natural free vibration frequency (ω) in the displacement field, the following relations are assumed for the bending and shear deflections of the strip [48]:

wb ( x, y, t )

wb ( x, y )eiZt

(32a)

ws ( x, y, t )

ws ( x, y )eiZt

(32b)

According to Eq. (14), the mass matrix of the strip, M strip , is defined as:

M strip

³A m Z LnLmdA 2 T

(33)

where, L n and L m can be obtained from Eq. (18), and m and n represent the harmonic function mode. The mass matrix of the entire nanoplate, M t , is then obtained by assembling the mass matrices of all strips as stated in the extraction of stiffness matrix process. Similar to the previous section, using Eq. (12) in the absence of potential energy, natural free vibration frequencies of the plate are obtained by solving the following eigenvalue problem:

K

tt



 Mt δ

0

(34)

in which, K tt is the stiffness matrix of the plate predefined in Eq. (31). 3.4. Multi-Layered Nanoplates Placement of single-layer graphene sheets on each other generates graphite, one of the most stable structures of carbon. The main reason for the placement of these plates on each other is

16

the four capacity of carbon atoms. In fact, each carbon atom is bonded to three other atoms on graphene sheet and forms a regular hexagon. Thus, at each atom, a bond capacity remains and as graphene sheets approach each other, the remaining capacity in each atom is bonded to the front plate and forms a multi-layered graphene sheet (MLGS). Unlike the carbon-carbon (CC) bonds inside the plates, this bond is weak and known as van der Waals (vdW) type. When analyzing the behavior of MLGSs, they can be considered as a stack of SLGSs; however, it is important to simulate the chemical vdW bonds between the sheets. In a general case, each layer can be closely modeled as a plate and the chemical bonding forces between the layers could be modeled as vdW links. Deformation of layers causes an increase in the pressure between the layers and this pressure, Ri, can be expressed as: NL

Ri

¦ cij (wi  w j )

(35)

j 1

where wi and w j are displacements of the ith and jth layers, respectively, and NL represents the number of the layers. Also, the coefficient cij , which could be assumed as the stiffness of vdW links, is derived by differentiating the potential function given in Ref. [49] as: cij

(

4 3 2 24H V 8 ­° 3003S ) ( ) ® 9acc V 2 acc °¯ 256

( 1) k § 5 · V 6 1 35S ¦ ¨ ¸ ( ) 12  8 hij k 0 2k  1 © k ¹ acc 5

( 1) k § 2 · 1 ½° ¦ ¨ ¸ 6¾ k 0 2 k  1 © k ¹ hij ° ¿ 2

(36)

in which acc = 1.42 Å is the C-C bond length, ε and σ are parameters associated with the physical properties of the material and are assumed to be ε = 2.968 meV and σ = 3.407Å, and hij is the normalized distance between the ith and jth layers and is defined as:

hij

zi  z j

(37)

acc

where z denotes the coordinate of the layers in direction of the thickness. The added potential energy due to the pressure between the layers, U van , is then obtained by:

17

U van

NL ½° 1 °­ NL 2 c w dA cij w2j dA¾  ®¦ ij ³A i ¦ ³ A 2 ¯° j 1 j 1 ¿°

(38)

Using Eq. (38), the stiffness of the system corresponding to vdW bonds is derived as:

K van

K van1  K van 2

NL









(¦ cij ) ³ LTn L m dA  ³ cij LTn L m dA j 1

A

A

(39)

in which Kvan1 and Kvan2 stiffness matrices are related to the degrees of freedom of the ith and jth layers, respectively. Considering the vdW bonds, the total stiffness matrix of MLGSs, i.e. Ktml, may finally be written as:

K tml

(K van 2 )2 ª(K tt )1  (K van1 )1 « (K tt )2  (K van1 ) 2 « (K van 2 )1 « « « (K van 2 )1 (K van 2 )2 ¬

(K van 2 ) NL

º » (K van 2 ) NL » » » (K tt ) NNL  (K van1 ) NL »¼

(40)

As a result, geometric stiffness matrix of a MLGS, i.e. KGml, is obtained as:

K Gml

0 ª(K Gt )1 « 0 (K Gt ) 2 « « « 0 ¬ 0

º 0 »» » » (K Gt ) NL ¼ 0

(41)

Similarly, the assembled mass matrix of a MLGS, i.e. M tml , could be written as:

M Gml

0 ª(M Gt )1 « 0 (M Gt ) 2 « « « 0 ¬ 0

º » 0 » » » (M Gt ) NL ¼ 0

(42)

Finally, the critical buckling load and the natural free vibration frequency of the MLGS are obtained by solving the following eigenvalue problems:

 K tml  K Gml δ

0

(43a)

 K tml  MGml δ

0

(43b)

18

4. Results and Discussion In the previous sections, the finite strip formulations were extended based on the MCST and TVRPT to evaluate the vibration and buckling of SLGSs and MLGSs. This section presents a number of numerical experiments solved based on these pre-defined formulations by computer programming in MATLAB software. The thickness and material properties of the graphene sheets used in all examples except for the validation section, are presented in Table 2.

Table 2. material properties of the graphene sheets Quantity

Value

Unit

Thickness (h)

0.34

nm

Young's modulus (E)

1.06

TPa

Poisson's ratio (ν)

0.25

-

Mass density (ρ)

2250

kg/m3

4.1. Convergence study Since the accuracy of the finite strip method depends on the number of strips, a convergence study is first carried out for a 10nm square isotropic SLGS with all edges simply supported (SSSS). For evaluating the convergence of mechanical buckling loads, a SSSS SLGS is subjected to uniaxial uniform compression loading as shown in Fig. 3a.

(a)

(b) SLGS

19

(c)

σy

σy

σx

(d)

(e) MLGS

Fig. 3. Nanoplate load patterns for SLGSs under: (a) uniaxial and (b) biaxial compression, and also (c) shear loadings; and MLGSs under: (d) uniaxial and (e) biaxial compression loadings

Table 3 shows the buckling loads per unit length of the nanoplate for different values of the LSP over thickness ratio (l/h) and numbers of strips (1 to 10 strips). Comparing the values of the buckling loads obtained using different numbers of strips, the convergence trend of the method can be observed with increasing the number of strips. Moreover, it can be seen that for l/h ratios of 0, 0.5 and 1, the results are converged after using, 4, 6 and 8 strips, respectively. Thus, when using higher length scale parameters (l), more strips are required to converge.

Table 3. Convergence of mechanical buckling load (N/m) for square (a = 10nm) all edges simply supported SLGS Number of strips l/h

1

2

4

6

8

10

0

1.5465

1.4563

1.4534

1.4532

1.4532

1.4532

0.5

3.2494

3.0974

3.0925

3.0923

3.0922

3.0922

1

8.3575

8.0199

8.0092

8.0087

8.0086

8.0086

Another convergence study is also performed for the free vibration analysis of the same SLGS. Non-dimensional natural free vibration frequencies of the SLGS defined in Eq. (44), are presented in Table 4 for different numbers of strips and different values of LSP to thickness ratio (l/h).

20

Z Z

a2 h

U

(44)

E

It can be also observed from this table that for plates with different values of l/h, by increasing the number of strips to six, the results are quite converged.

Table 4. Convergence of non-dimensional natural frequency ( Z ) of square (a = 10nm) all edges simply supported SLGS Number of strips l/h

1

2

4

6

8

10

0

6.0527

5.8735

5.8677

5.8673

5.8673

5.8673

0.5

8.7732

8.5660

8.5593

8.5588

8.5588

8.5588

1

14.0696

13.7837

13.7742

13.7742

13.7742

13.7742

4.2. Validation of the proposed method The accuracy of the method presented in this paper is evaluated by comparing the results with those available in relevant references. The number of harmonic shape functions used in the finite strip method is different depending on the boundary conditions, applied load type and the aspect ratio of the nanoplate. The study is first carried out for the mechanical buckling of a square SLGS with two kinds of boundary conditions, i.e. all edges simply supported (SSSS) and all edges clamped (CCCC). The results are compared with the corresponding ones obtained in Ref. [50], in which the behavior of the functionally graded materials (FGM) micro-plates were evaluated under different loading patterns and diverse boundary conditions using a refined quasi-3D isogeometric analysis. In this section, the results chosen for the validation study are those given for a special case of FGM micro-plates, with 100% homogenous ceramic and no metal. The geometrical and material properties of this micro-plate are presented in Table 5.

21

Table 5. The geometrical and material properties of micro-plate for validation of mechanical buckling of the proposed method Quantity Value Unit Thickness (h)

17.6

μm

Young's modulus (E)

380

GPa

Poisson's ratio (ν)

0.3

-

dimensions of plate (a, L)

100h

-

It should be noted that we used one harmonic shape function in the finite strip method for SSSS and eight harmonic shape functions for CCCC boundary conditions. Table 6 presents the values of non-dimensional critical buckling load defined in Eq. (45) for the square SSSS and CCCC SLGS micro-plates with different l/h ratios.

N cr

12 N cr a 2 (1 Q 2 ) Em h3

(45)

In Eq. (45), N cr is the critical buckling load and Em is Young's modulus of aluminum which is assumed to be 70GPa. As can be seen from the difference percentage values shown in Table 6, i.e. less than 1.1%, the present results are in an excellent agreement with those obtained in Ref. [50].

Table 6. Dimensionless critical buckling load for square micro-plate under uniaxial loading l/h

Boundary conditions SSSS Present

CCCC Ref. [50]

Difference (%)

Present method

Ref. [50]

method

Difference (%)

0

107.0993

107.0958

0.0033

285.0190

283.7646

0.4421

0.4

179.0898

179.0825

0.0041

427.9718

425.0550

0.6862

0.6

269.0734

269.0653

0.0030

606.4589

601.2918

0.8593

1

557.0264

557.0082

0.0033

1177.1795

1164.4998

1.0889

As the second validation problem, the free vibration of a square SSSS SLGS is studied using one harmonic shape function. The non-dimensional frequencies obtained for this problem

22

using the present method are compared with the corresponding ones obtained in Ref. [50], in which the FG micro-plates with different sizes are studied. Here, similar to the previous validation example, FG plates with 100% ceramic material (no metal) are used. The geometrical and material properties of this micro-plate are presented in Table 7. The values of non-dimensional natural free vibration frequency defined in Eq. (44) are presented in Table 8 for the square SSSS SLGS micro-plates with the width of a=10nm and different l/h ratios. Considering the low values (less than 0.01%) of the difference percentage between the results of the present method and those of Ref. [50], it can be concluded that an excellent agreement exists between them.

Table 7. The geometrical and material properties of micro-plate for validation of free vibration of the proposed method Quantity Value Unit Thickness (h)

0.1

μm

Young's modulus (E)

380

GPa

Poisson's ratio (ν)

0.3

-

Dimensions of plate (a, L)

100h

-

mass density (ρ)

2250

kg/m3

Table 8. Dimensionless natural frequencies for square SSSS micro-plate l/h

0

0.2

0.4

0.6

0.8

1.0

Present method

5.9717

6.4540

7.7221

9.4654

11.4690

13.6190

Ref. [50]

5.9712

6.4534

7.7215

9.4646

11.4682

13.6178

Difference (%)

0.00837

0.00930

0.00777

0.00845

0.00698

0.00881

Similar to the validation process performed for the SLGSs, the accuracy of the proposed method for the MLGSs is evaluated by analyzing their mechanical buckling and free vibration responses. The geometrical and material properties of the MLGS examples used for validation are presented in Table 2. Depending on the number of layers, different types of buckling and

23

vibration are observed due to the vdWs bonds between the layers, such that for a n-layered graphene sheets, there are n types of buckling and vibration. The first three modes for all types of buckling and vibration for two-layered, three-layered and four-layered graphene sheets are schematically shown in Fig. 4. As shown in this figure, the first buckling mode of each type is the same as that in the SLGSs, as no relative displacement occurs between the layers and vdW bonds forces do not have any contribution in this mode.

First mode of first type

second mode of first type

third mode of first type

First mode of second type

second mode of second type

third mode of second type

(a)

First mode of first type

second mode of first type

third mode of first type

First mode of second type

second mode of second type

third mode of second type

First mode of third type

second mode of third type

third mode of third type

(b)

First mode of first type

second mode of first type

24

third mode of first type

First mode of second type

second mode of second type

third mode of second type

First mode of third type

second mode of third type

third mode of third type

First mode of forth type

second mode of forth type

third mode of forth type

(c) Fig. 4. Different buckling and vibration shapes for: (a) two-layered, (b) three-layered, and (c) four-layered graphene sheets

For validation of the mechanical buckling analysis of MLGSs by the proposed method, it is assumed that a three-layered square graphene sheets (3LGSs) with the width of a=10nm and all edges simply supported (SSSS) is subjected to uniaxial compression (Fig. 3). In the analysis of a single-layered graphene sheet with SSSS boundary conditions, it was sufficient to use only one harmonic shape function. However, for the analysis of MLGSs, in order to capture higher modes of each type of buckling, it is necessary to use more shape functions. Thus, twelve harmonic shape functions are used in this study. Table 9 presents the critical buckling load values of this problem obtained using the present method for three types of buckling (Fig. 4b). The results are also compared with the corresponding ones given in Ref. [51], in which multilayered graphene sheets are analyzed based on the nonlocal Eringen theory using the finite strip method. The low values of difference percentages between the results of the present method and those of Ref. [51] represent the high accuracy of the proposed method for the mechanical buckling analysis of MLGSs.

Table 9. First mode critical buckling load (N/m) for square SSSS threelayered GS (a = 10nm) under uniaxial loading Type of buckling

First

Second

Third

Present method

1.4532

117.1160

355.9261

25

Ref. [51]

1.4620

117.2359

356.0642

Difference (%)

0.602

0.102

0.039

In order to validate the accuracy of the present method for the vibration analysis of MLGSs, two-layered and three-layered square graphene sheets (2LGSs and 3LGSs) with the width of a=10nm and all edges simply supported (SSSS), are considered. The values of natural free vibration frequency obtained by the present method are reported in Table 10 and compared with the corresponding ones given in Ref. [51]. It can be observed that an excellent agreement exists between the results, showing the high accuracy of the present method in the vibration analysis of MLGSs.

Table 10. Natural frequencies (THz) for square SSSS MLGSs (a = 10nm) Natural frequency

Two-layered graphene sheet Present method Ref. [51]

Three-layered graphene sheet Present method Ref. [51]

Z1

0.06891

0.06912

0.06891

0.06912

Z2

2.62801

2.68328

1.82576

1.86414

Z3

----

-----

3.21827

3.28597

4.3. Mechanical Buckling of SLGSs In this section, mechanical buckling of isotropic SLGSs under different load patterns illustrated in Fig. 3, is studied. In all assessments, the width of the plate, i.e. a, is assumed to be 10 nm and different boundary conditions shown in Fig. 5 are considered.

S F

F

S

S

S

S C

C

C

S

S

S

C

S

S

SSFF

SCSC

SSCC

SSSS

Fig. 5. Boundary conditions labeling pattern

26

In order to obtain accurate results, the number of strips is selected in accordance with the convergence study section and different numbers of harmonic shape functions are assumed for different boundary conditions. Considering various LSP to thickness ratios (l/h), the nondimensional critical buckling load described in Eq. (46) is computed for SLGSs with different values of aspect ratios (L/a) and different boundary conditions.

N cr

N cr a 2 Eh3

(46)

Tables 11 to 13 show the results of Ncr for SLGSs with various boundary conditions under three types of loading pattern: uniaxial uniform compression (Fig. 3a), equal biaxial uniform compression (Fig. 3b) and uniform shear loading (Fig. 3c), respectively. For all loading patterns and boundary conditions, it is observed that the critical buckling load values of SLGSs are steadily increased by increasing l/h values. Particularly, applying LSP leads to the greatest increase in the Ncr values for SLGSs with SSSS boundary conditions subjected to uniaxial compression load (Table 11), as well as those with SSCF boundary conditions under the other two loading patterns, i.e. biaxial uniform compression and uniform shear loading (Tables 12 and 13).

Table 11. Non-dimensional buckling forces for SLGSs under uniform uniaxial loading Boundary conditions

L/a

CCCC

CCCS

SSCC

l/h 0

0.2

0.4

0.6

0.8

1

1

8.6762

9.8374

13.3145

19.1000

27.1925

37.5921

1.5

7.2012

8.2262

11.2930

16.3930

23.5244

32.6876

2

6.7874

7.7737

10.7248

15.6317

22.4918

31.3063

1

6.9751

7.9583

10.9022

15.8006

22.6519

31.4568

1.5

6.0594

6.9562

9.6404

14.1053

20.3496

28.3739

2

5.4011

6.2414

8.7571

12.9419

18.7938

26.3135

1

6.6264

7.5937

10.4872

15.2968

22.0209

30.6601

27

SCCS

SSCS

SSSS

CCFF

SSCF

SSFF

CFFF

1.5

6.1605

7.0952

9.8927

14.5446

21.0490

29.4060

2

6.0304

6.9507

9.7042

14.2820

20.6816

28.9039

1

5.4167

6.2584

8.7777

12.9683

18.8288

26.3591

1.5

5.1303

5.9440

8.3798

12.4302

18.0941

25.3717

2

4.9239

5.7212

8.1078

12.0766

17.6258

24.7562

1

4.9937

5.8014

8.2194

12.2423

17.8677

25.0966

1.5

4.7160

5.4937

7.8210

11.6904

17.1011

24.0532

2

4.8617

5.6502

8.0099

11.9338

17.4201

24.4699

1

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

1.5

3.7758

4.4314

6.3976

9.6745

14.2618

20.1593

2

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

1

3.4237

3.7677

4.7182

6.2272

8.2979

10.9368

1.5

0.0941

1.6823

2.1074

2.7653

3.6568

4.7863

2

0.8559

0.9474

1.1874

1.5503

2.0354

2.6460

1

1.4838

1.7966

2.6141

3.8570

5.5335

7.6578

1.5

1.1720

1.4798

2.3033

3.5572

5.2412

7.3690

2

1.4838

1.7966

2.6141

3.8570

5.5335

7.6578

1

0.8468

0.9443

1.1862

1.5374

1.9992

2.5767

1.5

0.3730

0.4188

0.5273

0.6745

0.8593

1.0854

2

0.2086

0.2352

0.2964

0.3749

0.4683

0.5794

1

0.2139

0.2371

0.2972

0.3860

0.5014

0.6442

1.5

0.0943

0.1051

0.1320

0.1697

0.2163

0.2722

2

0.0528

0.0590

0.0742

0.0946

0.1191

0.1471

Table 12. Non-dimensional buckling forces for SLGSs under equal biaxial compression Boundary conditions

L/a

CCCC

CCCS

SSCC

l/h 0

0.2

0.4

0.6

0.8

1

1

4.5804

5.1970

7.0440

10.1173

14.4157

19.9396

1.5

3.5716

4.0207

5.3629

7.5939

10.7128

14.7200

2

3.4021

3.7946

4.9637

6.9034

9.6130

13.0941

1

3.7343

4.2667

5.8605

8.5111

12.2173

16.9794

1.5

2.3820

2.7168

3.7187

5.3845

7.7137

10.7061

2

2.0335

2.2944

3.0747

4.3713

6.1840

8.5130

1

3.3215

3.7984

5.2280

7.6095

10.9423

15.2265

28

SCCS

SSCS

SSSS

CCFF

SSCF

SSFF

CFFF

1.5

3.2682

3.6683

4.8689

6.8698

9.6706

13.2717

2

3.3172

3.6839

4.7842

6.6180

9.1851

12.4854

1

2.9079

3.3625

4.7232

6.9862

10.1504

14.2160

1.5

2.1485

2.4607

3.3949

4.9479

7.1194

9.9097

2

1.9476

2.2001

2.9545

4.2088

5.9630

8.2173

1

2.3172

2.6967

3.8342

5.7282

8.3781

11.7841

1.5

1.9756

2.2560

3.0975

4.5000

6.4629

8.9866

2

1.8784

2.1134

2.8181

3.9928

5.6375

7.7519

1

1.7440

2.0587

3.0030

4.5766

6.7793

9.6113

1.5

1.2616

1.4723

2.1045

3.1583

4.6332

6.5294

2

1.0926

1.2541

1.7387

2.5464

3.6774

5.1313

1

2.4809

3.0963

4.7119

6.1924

8.2314

10.8340

1.5

0.9762

1.2681

2.1053

2.7461

3.6155

4.7199

2

0.4990

0.6697

1.1675

1.5398

2.0107

2.6045

1

1.0446

1.3204

2.0033

2.9910

4.3017

5.9550

1.5

0.5818

0.7688

1.2829

2.0621

3.0932

4.3850

2

0.4215

0.5559

0.9445

1.5650

2.4072

3.4715

1

0.8353

0.9414

1.1857

1.5318

1.9865

2.5556

1.5

0.3689

0.4176

0.5271

0.6711

0.8514

1.0722

2

0.2069

0.2345

0.2964

0.3730

0.4637

0.5715

1

0.2131

0.2369

0.2972

0.3852

0.4995

0.6406

1.5

0.0941

0.1051

0.1320

0.1692

0.2155

0.2705

2

0.0528

0.0590

0.0742

0.0946

0.1186

0.1464

Table 13. Non-dimensional buckling forces for SLGSs under shear loading Boundary conditions

L/a

CCCC

CCCS

SSCC

l/h 0

0.2

0.4

0.6

0.8

1

1

12.3978

14.2695

19.8583

29.1357

42.0971

58.7453

1.5

9.7868

11.2003

15.4154

22.4037

32.1608

44.6891

2

8.7703

10.0573

13.8929

20.2450

29.1053

40.4762

1

11.3767

13.1337

18.3828

27.1006

39.2819

54.9291

1.5

8.0879

9.3106

12.9652

19.0359

27.5196

38.4178

2

7.3407

8.4602

11.7041

17.0236

24.4399

33.9588

1

10.7143

12.3836

17.3762

25.6736

37.2726

52.1748

1.5

9.2187

10.6264

14.8201

21.7646

31.2403

43.3988

29

SSCS

SSSS

CCFF

SSCF

SSFF

2

8.5737

9.8230

13.5492

19.7241

28.3403

39.3999

1

9.1728

10.6939

15.2421

22.7995

33.3619

46.9310

1.5

7.7034

8.8644

12.3371

18.1109

26.1844

36.5585

2

6.9646

8.0406

11.2459

16.5394

23.8746

33.2330

1

8.0176

9.4277

13.6428

20.6439

30.4256

42.9895

1.5

6.1132

7.1170

10.1231

15.1245

22.1193

31.1078

2

5.6699

6.5311

9.1070

13.3901

19.3789

27.0746

1

6.5285

7.7732

11.2222

16.6554

24.0683

33.1331

1.5

2.9406

3.5663

4.9416

7.0644

9.9243

13.5355

2

1.5726

1.8583

2.6208

3.7770

5.3173

7.2502

1

4.4093

5.4476

8.2418

12.5613

18.4227

25.8647

1.5

2.5668

3.2099

4.9829

7.7595

11.5399

16.3429

2

2.0426

2.5224

3.8582

5.9673

8.8481

12.5130

1

3.7475

4.4117

6.1226

8.6709

12.0457

16.2701

1.5

1.4325

1.7251

2.4672

3.5351

4.9119

6.6084

2

0.7350

0.9020

1.3223

1.9120

2.6530

3.5493

In Figs Fig. 6 and Fig. 7, the effect of boundary conditions on variation of the Ncr values with length to width ratio (L/a) is demonstrated for square SLGSs with two different LSP to thickness ratios (i.e. l/h=0 and 1). Thus, in each figure, three different boundary conditions including SSSS, SSCC and CCCC are considered for three different loading patterns including uniaxial, biaxial and shear loadings. The general trend of all curves illustrated in Figs Fig. 6 and Fig. 7 indicates that the critical buckling load decreases by increasing L/a ratios and eventually becomes constant after a certain value of L/a. Moreover, as expected, by reducing the degrees of freedom at the edges of the nanoplates, i.e. CCCC boundary conditions as opposed to the SSSS ones, the critical buckling load is considerably increased (up to ~3 times). Moreover, as the L/a ratio increases, the effect of the boundary conditions of the loaded edges on the buckling force becomes less pronounced compared to other edges. Therefore, for plates with high L/a aspect ratios and SSCC and CCCC boundary conditions (the unloaded edges are CC), the buckling forces are very close to each

30

other. 7 ssss l/h=0 sscc cccc

6

Ncr (N/m)

5 4 3 2 1 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a

(a) 6 ssss l/h=0 sscc cccc

5

Ncr (N/m)

4 3 2 1 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a

(b) 14 ssss l/h=0 sscc cccc

12

Ncr (N/m)

10 8 6 4 2 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a (c) Fig. 6. Buckling force vs. aspect ratio for SLGSs with l/h=0 under: (a) uniaxial and (b) biaxial compression, and (c) shear loadings

31

30 ssss l/h=1 sscc cccc

Ncr (N/m)

25 20 15 10 5 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a

(a) 25

ssss l/h=1 sscc cccc

Ncr (N/m)

20 15 10 5 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a

(b) 80 ssss l/h=1 sscc cccc

Ncr (N/m)

60 40 20 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

L/a

(c) Fig. 7. Buckling force vs. aspect ratio for SLGSs with l/h=1 under: (a) uniaxial and (b) biaxial compression, and (c) shear loadings

Fig. 8 demonstrates the variation of Ncr values with thickness to LSP ratio (h/l) for square

32

SLGSs with L=a=10nm and SSSS, SSCC and CCCC boundary conditions under the three afore-mentioned load patterns. From the curves shown in Fig. 8, it can be concluded that by increasing the h/l ratio, the critical buckling load value is dramatically decreased. However, it can be observed that after certain h/l values, the critical buckling load values start to converge and eventually become constant.

20 17.5 15 12.5 10 7.5 5

6 4 2

2.5 0

ssss sscc cccc

8 Ncr (N/m)

Ncr (N/m)

10

ssss sscc cccc

0

5

10

15

0

20

h/l

0

5

10

15

20

h/l

(a)

(b) 25 ssss sscc cccc

22.5

Ncr (N/m)

20 17.5 15

12.5 10 7.5 5 2.5 0

0

5

10

15

20

h/l

(c) Fig. 8. Buckling force vs. thickness to LSP ratio (h/l) for SLGSs with L=a=10nm under: (a) uniaxial and (b) biaxial compression, and (c) shear loadings

Another important topic of discussion in this section is the interaction of axial and shear loads. In Fig. 9, the interaction curves of the uniaxial loads in two orthogonal directions

33

normalized with respect to their corresponding critical buckling loads, i.e. Nx/Ncrx and Ny/Ncry, as well as those of normalized shear and uniaxial loads, i.e. Nxy/Ncrxy and Ny/Ncry, are illustrated for SLGSs with different l/h and L/a ratios and SSSS boundary conditions. It can be observed that, while the effects of l/h and L/a ratios on the interaction curves of the Nx/Ncrx and Ny/Ncry are quite considerable, they have negligible effect on the interaction curves of the Nxy/Ncrxy and Ny/Ncry. This is due to the phase difference between the axial and shear forces that does not exist in axial forces. In fact, considering the interaction of the axial forces, since the wavelengths are perpendicular, the corresponding axial forces have strong effects on each other [52].

2

l/h=1 l/h=0 l/h=1 l/h=0 l/h=1 l/h=0

Nx/Ncrx

1.5

L/a=1 L/a=1 L/a=2 L/a=2 L/a=4 L/a=4

1

0.5

0 -1

-0.5

0 Ny/Ncry (a)

34

0.5

1

1.5

l/h=1 L/a=1 l/h=0 L/a=1 l/h=1 L/a=2 l/h=0 L/a=2 l/h=1 L/a=4 l/h=0 L/a=4

Nxy/Ncrxy

1

0.5

0 -1

-0.5

0

0.5

1

Ny/Ncry (b) Fig. 9. Interaction curves for buckling of SLGS nano-plate with SSSS boundary conditions under: (a) biaxial compression loading and (b) combined compression and shear loading

4.4. Free Vibration of SLGSs In this section, applying the equations derived in section 3.3, the free vibration of SLGSs with different boundary conditions is evaluated. In Table 14, the variations of non-dimensional natural free vibration frequency ( Z defined in Eq. (44)) with l/h ratio are presented for SLGSs with various boundary conditions and L/a ratios of 1 and 2. The width of the SLGS is assumed to be a=10nm and other material and geometrical properties are as given in Table 2. According to Table 14, by increasing the l/h ratio, the natural free vibration frequencies of the modeled SLGSs are increased up to ~2.6 times. Changing the dimensions of the sheet from square to rectangular is also accompanied by a decrease in the values of natural frequency. This decrease becomes more significant when increasing the LSP to thickness ratio (l/h).

Table 14. Non-dimensional natural frequencies of SLGSs with different boundary conditions, L/a and l/h ratios l/h L/a Boundary conditions 0 0.2 0.4 0.6 0.8 1

35

1

2

SSFF

2.8908

3.0528

3.4219

3.8959

4.4427

5.0435

SSCF

3.8270

4.2109

5.0793

6.1697

7.3902

8.6937

SSSS

5.8674

6.3742

7.6995

9.5044

11.5682

13.7742

CCFF

6.7227

7.0447

7.8823

7.8823

10.4794

12.0456

SSCS

7.0206

7.5667

9.0071

10.9923

13.2796

15.7377

SSCC

8.5789

9.167

10.7375

12.9368

15.4965

18.2676

CCCS

9.4705

10.1074

11.8121

14.1983

16.983

19.998

CCCC

10.6806

11.3608

13.1929

15.7729

18.7974

22.0835

SSFF

0.7173

0.7617

0.8555

0.9619

1.0752

1.1959

SSCF

1.7306

1.9507

2.4581

3.0833

3.7653

4.4813

SSSS

3.6713

3.9333

4.6315

5.6051

6.7356

7.9565

CCFF

1.6834

1.7669

1.9773

2.2646

2.6025

2.9756

SSCS

5.149

5.4594

6.2994

7.4924

8.8975

10.4292

SSCC

7.0545

7.4325

8.4651

9.9516

11.7186

13.659

CCCS

5.4696

5.8059

6.7133

7.9983

9.5089

11.1539

CCCC

7.2943

7.694

8.7835

10.3459

12.2023

14.2363

In order to show the effect of length to width ratio of the SLGSs (L/a) on their natural free vibration frequency, variations of Z with L/a for SLGSs with SSSS, SSCC and CCCC boundary conditions are plotted in Fig. 10. It can be observed that the natural free vibration frequencies decrease by increasing the L/a aspect ratio for all three boundary conditions. However, after certain values of L/a, which varies for different l/h values and boundary conditions, the Z values start to converge and are eventually stabilized at constant values. It is obvious from these curves that the natural free vibration frequency converges faster for the CCCC boundary conditions than for the SSSS. The effect of the length scale parameter (l) on the natural free vibration frequency is also depicted in Fig. 11 for square SLGS nanoplates with the three afore-mentioned boundary conditions, i.e. SSSS, SSCC and CCCC. According to this figure, by increasing the h/l ratios, the natural free vibration frequency decreases in all three types of boundary conditions. Therefore, it can be concluded that after a significant increase in the thickness of the nanoplate (large values of h/l), the natural free vibration frequency is no longer affected by the length 36

scale parameter.

(a)

(b) Fig. 10. Non-dimensional natural frequency for SLGSs with different aspect ratios: (a) l/h=0 and (b) l/h=1

37

Fig. 11. Non-dimensional natural frequency vs. h/l for different boundary conditions

4.5. Mechanical Buckling of MLGSs Based on the MCST and the finite strip formulation developed in section 3.4, the stability of multi-layered graphene sheets is evaluated in this section, considering the effects of vdW bonds between the layers. In this section, three types of MLGSs including two, three and fourlayered graphene sheets, i.e. 2LGSs, 3LGSs and 4LGSs, respectively, are considered. The MLGSs are subjected to the loading patterns shown in Fig. 3 and their boundary conditions are assumed to be SSSS and SSCC. In order to obtain converged results, each layer is divided into 10 strips, and 15 and 25 harmonic shape functions are used for the square and rectangular MLGS nanoplates, respectively. The non-dimensional buckling loads (defined in Eq. (46)) of the afore-mentioned MLGSs are presented in Tables 15 and 16. According to these tables, buckling loads for the higher buckling types are significantly different from those of the first type in which the vdW bonds do not exist. Meanwhile, by increasing the number of layers, the difference between critical buckling loads of different types decreases. In MLGSs, middle layers are affected by vdW bonds from both sides and as a result, the layers have neutralizing effect on each other. In addition, by increasing the l/h ratio, critical buckling loads are increased in all types of buckling. However, the rate of increasing significantly reduces for the higher types in which 38

the effects of vdW bonds are prominent. Moreover, increasing the L/a aspect ratio generally leads to a decrease in the buckling loads for all types of buckling; however, this reduction is not significant in the second, third and fourth buckling types. In fact, the increase in the dimensions affects the critical force. However, the magnitudes of the critical forces are very large when considering the van der Waals force which diminish the effect of dimensions’ ratio.

Table 15. Non-dimensional buckling forces of different MLGSs under uniaxial loading

BCs

SSSS

Nu mbe r of laye rs

2

L/a

1

1.5

2

3

1

1.5

2

4

1

1.5

Buckli ng type

l/h 0

0.2

0.4

0.6

0.8

1

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

119.2070

126.5518

146.3603

174.3437

208.4982

245.7177

First

3.7758

4.4314

6.3976

9.6745

14.2618

20.1593

Second

119.2070

126.6679

146.3123

174.3437

207.8616

245.1071

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

119.2070

126.5518

146.3603

174.3437

208.4982

245.7177

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

86.5724

91.7701

105.3049

125.9891

149.0443

177.8054

Third

120.0339

151.4633

175.0729

209.4667

250.3087

297.0534

First

3.7758

4.4314

6.3976

9.6745

14.2618

20.1593

Second

86.5724

91.7701

105.3441

125.1200

149.0443

176.4922

Third

141.4151

150.7216

175.0729

209.8692

250.3087

295.2722

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

86.5724

91.6660

105.3049

125.0605

149.0467

176.0234

Third

141.4187

150.6198

175.0729

209.4667

250.3087

295.0105

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

67.4534

71.5914

84.0044

97.6798

116.1735

137.6275

Third

118.3098

125.9407

145.2526

172.9417

206.6672

243.8866

Forth

149.0186

159.1789

185.4293

222.2523

266.4905

313.2352

First

3.7758

4.4314

6.3976

9.6745

14.2618

20.1593

Second

67.3456

71.1565

81.6977

97.5627

116.0415

139.7985

Third

118.2858

125.6380

145.1177

172.9417

206.1915

243.0882

39

2

SSCC

2

1

1.5

2

3

1

1.5

2

4

1

1.5

2

Forth

149.0186

158.9302

185.3921

221.8524

265.3864

313.2352

First

3.4881

4.1174

6.0059

9.1531

13.5586

19.2227

Second

67.4534

71.1452

81.6977

97.1474

116.1735

137.6275

Third

118.3098

125.5576

145.2526

172.9417

206.5023

243.4557

Forth

149.0186

159.4324

185.2006

221.9883

265.1408

313.2352

First

6.6264

7.5937

10.4872

15.2968

22.0209

30.6601

Second

119.3251

126.9848

146.6263

174.8262

209.3346

246.9066

First

6.1605

7.0952

9.8927

14.5446

21.0490

29.4060

Second

119.3251

126.8304

146.5913

174.8262

208.6573

246.3579

First

6.0304

6.9507

9.7042

14.2820

20.6816

28.9039

Second

119.3251

126.7102

146.6263

174.8262

208.8983

246.8638

First

6.6264

7.5937

10.4872

15.2968

22.0209

30.6601

Second

86.7392

91.9727

105.6693

126.5508

150.0337

179.2033

Third

120.1352

150.7475

175.3060

209.8889

251.0340

298.0903

First

6.1605

7.0952

9.8927

14.5446

21.0490

29.4060

Second

86.7392

91.9727

105.6878

125.7141

150.0337

178.0797

Third

141.5826

150.8553

175.3055

210.3103

251.0340

296.4053

First

6.0304

6.9507

9.7042

14.2820

20.6816

28.9039

Second

86.7392

91.8878

105.6695

125.6718

150.0337

177.5584

Third

141.5126

150.7475

175.3060

209.8891

251.0340

296.1178

First

6.6264

7.5937

10.4872

15.2968

22.0209

30.6601

Second

67.6567

71.8358

82.1394

98.3507

121.0413

139.3314

Third

118.4276

126.0876

145.5186

173.4242

207.5035

245.0754

Forth

149.1062

159.3066

185.6624

222.6748

267.2158

314.2721

First

6.1605

7.0952

9.8927

14.5446

21.0490

29.4060

Second

67.5638

71.4177

82.1394

98.3327

117.2419

139.3026

Third

118.4276

125.8005

145.3966

173.4242

206.9889

244.3390

Forth

149.1062

159.9576

185.6624

222.2581

266.0529

314.2721

First

6.0304

6.9507

9.7042

14.2820

20.6816

28.9039

Second

67.6567

71.4157

82.1394

97.8901

117.3871

139.3314

Third

118.4276

125.7160

145.5186

173.4242

207.2791

244.7401

Forth

149.1062

159.0584

185.4199

222.3860

265.8208

314.2721

Table 16. Non-dimensional buckling forces for MLGSs under bi-axial loading

40

BCs

SSSS

Number of layers

L/a

2

1

1.5

2

3

1

1.5

2

4

1

1.5

2

SSCC

2

1

1.5

2

Buckling type

l/h 0

0.2

0.4

0.6

0.8

1

First

1.7440

2.0587

3.0030

4.5766

6.7793

9.6113

Second

118.0268

125.5821

144.5755

171.6614

204.3282

240.8034

First

1.2616

1.4723

2.1045

3.1583

4.6332

6.5294

Second

118.0268

125.2304

144.3900

171.6614

204.0670

239.7135

First

1.0926

1.2541

1.7387

2.5464

3.6774

5.1313

Second

118.0268

125.1649

144.5362

171.6614

204.3282

239.9014

First

1.7440

2.0587

3.0030

4.5766

6.7793

9.6113

Second

85.2405

90.3583

103.1990

123.2488

145.0160

172.4778

Third

140.4754

149.3854

173.3397

206.9121

246.4577

291.4946

First

1.2616

1.4723

2.1045

3.1583

4.6332

6.5294

Second

85.2405

90.2321

103.4207

122.3669

147.5480

170.4980

Third

140.3227

149.4084

173.3397

207.1118

246.4577

289.8817

First

1.0926

1.2541

1.7387

2.5464

3.6774

5.1313

Second

85.2405

90.0650

103.1990

122.1689

145.0160

170.3907

Third

140.3573

149.3852

173.3397

206.9121

246.4577

289.8574

First

1.7440

2.0587

3.0030

4.5766

6.7793

9.6113

Second

65.8637

70.1580

79.4895

95.0398

111.7052

132.3342

Third

117.1384

124.6208

143.4812

170.2811

202.5338

239.0088

Forth

147.9911

157.8741

183.5936

219.5420

262.3906

308.4162

First

1.2616

1.4723

2.1045

3.1583

4.6332

6.5294

Second

65.8637

69.5908

79.4895

94.2494

112.1005

131.4279

Third

117.1384

124.2120

143.2110

170.2811

202.4291

237.7390

Forth

147.9911

157.7025

183.5936

219.3346

261.8995

308.4162

First

1.0926

1.2541

1.7387

2.5464

3.6774

5.1313

Second

65.8637

70.1597

79.4895

94.0386

111.7052

131.7738

Third

117.1382

124.1815

143.3111

170.2811

202.5338

237.8266

Forth

147.9911

157.8744

183.5360

219.5420

261.5207

307.9472

First

3.3215

3.7984

5.2280

7.6095

10.9423

15.2265

Second

118.0472

125.6466

144.7311

171.9771

204.8922

241.6774

First

3.2682

3.6683

4.8689

6.8698

9.6706

13.2717

Second

118.0472

125.2895

144.5462

171.9771

204.6222

240.6083

First

3.3172

3.6839

4.7842

6.6180

9.1851

12.4854

41

3

1

1.5

2

4

1

1.5

2

Second

118.0472

125.2263

144.6922

171.9771

204.8922

240.8046

First

3.3215

3.7984

5.2280

7.6095

10.9423

15.2265

Second

85.2801

90.4526

103.4083

123.6566

145.7135

173.6131

Third

140.5104

149.4351

173.4707

207.1847

246.9500

292.3016

First

3.2682

3.6683

4.8689

6.8698

9.6706

13.2717

Second

85.2801

90.3120

103.6296

122.7690

145.7135

171.6187

Third

140.3328

149.4562

173.4707

207.3871

246.9500

290.6788

First

3.3172

3.6839

4.7842

6.6180

9.1851

12.4854

Second

85.2801

90.1509

103.4085

122.5745

145.7135

171.4980

Third

140.3717

149.4351

173.4707

207.1847

246.9500

290.6476

First

3.3215

3.7984

5.2280

7.6095

10.9423

15.2265

Second

66.1726

70.2384

79.7475

95.5119

112.5412

133.5994

Third

117.1610

124.6745

143.6380

170.5979

203.1002

239.8853

Forth

148.0031

157.9176

183.7177

219.8048

262.8670

309.1660

First

3.2682

3.6683

4.8689

6.8698

9.6706

13.2717

Second

65.9132

69.7113

79.7475

94.7405

112.9140

132.7372

Third

117.1610

124.2725

143.3687

170.5979

202.9862

231.4362

Forth

148.0031

157.7479

183.7177

219.5950

262.3632

309.1660

First

3.3172

3.6839

4.7842

6.6180

9.1851

12.4854

Second

65.9614

69.6165

79.7475

94.5280

112.5412

133.0989

Third

117.7822

124.2442

143.4685

170.5979

203.1002

238.7334

Forth

148.0031

157.7925

183.6596

219.8048

261.9886

308.7172

In Figs. Fig. 12 and Fig. 13, the effect of L/a aspect ratio on the critical buckling load of MLGSs is demonstrated for a three-layered SSCC graphene sheet subjected to uniaxial and biaxial compressive loadings, respectively. In these examples, the width of the sheets in the x direction is assumed to be a=10nm and LSP is assumed to be equal to h. According to these figures, it can be concluded that while for both loading conditions the critical buckling load decreases by increasing L/a ratios, the rate of this decrease for biaxial compression loading is less than that for uniaxial compression.

42

79 78 Ncr(N/m)

77 76 75 74 73 0.5

1

1.5

2

2.5

3

3.5

4

2.5

3

3.5

4

L/a

(a) 126 125.5

Ncr(N/m)

125 124.5 124 123.5 123 0.5

1

1.5

2 L/a

(b) Fig. 12. Variation of the critical buckling load vs. aspect ratio for three-layered SSCC graphene sheets under uniaxial compression loading with l/h=1: (a) second type, (b) third type

73

Ncr(N/m)

72.5 72 71.5 71 0.5

1

1.5

2

2.5 L/a (a)

43

3

3.5

4

123.5

Ncr(N/m)

123 122.5 122 121.5 121 0.5

1

1.5

2

2.5

3

3.5

4

L/a (b) Fig. 13. Variation of the critical buckling load vs. aspect ratio for three-layered SSCC graphene sheets under biaxial compression loading with l/h=1: (a) second type, (b) third type

The effect of the length scale parameter is also depicted in Fig. 14 for a three-layered square 10nm graphene sheet with SSCC boundary conditions subjected to uniaxial and biaxial compression loadings. It is observed that for small values of h/l, the critical buckling load significantly varies, and these variations are more significant in the higher types of buckling. However, by increasing the thickness of the sheet above specific values, the change in the LSP no longer affects the critical buckling load results.

140 first form 120

second form third form

Ncr(N/m)

100 80 60 40 20 0

0

5

10 h/l

(a)

44

15

20

140 first form second forn third form

120

Ncr(N/m)

100 80 60 40 20 0

0

5

10

15

20

h/I

(b) Fig. 14. Variation of the critical buckling load vs. h/l ratio for three-layered SSCC graphene sheets under: (a) uniaxial and (b) biaxial compression loadings

4-6. Vibration of MLGSs To investigate the free vibration of MLGSs, a two-layered graphene sheet with free boundary conditions (FFFF) is considered. The trigonometric shape function used for this boundary condition is given in Table 1. Each layer is divided into 10 strips and 6 harmonic shape functions are used. In this investigation, the width of the sheet is assumed to be a=10nm and the aspect ratios of 1, 1.5 and 2 are considered. The non-dimensional natural free vibration frequency ( Z ) of FFFF two-layered grapheme sheets are presented in Table 17. Similar to the buckling problem, two types of vibrational deformation are available for a two-layered graphene sheet and thus, two Z values are reported in this table. According to Table 17, the vdW bonds increase the second type of natural free vibration frequency more significantly than the first type. By considering the small-scale effects, the first-type natural free vibration frequencies increase significantly; however, there are very small changes in the results of the second type which is attributed to overcoming the effects of vdW bonds.

Table 17. Non-dimensional natural frequencies for FFFF two-layered graphene sheets

45

L/a l/h

1

1.5

2

Z1

Z2

Z1

Z2

Z1

Z2

0

6.0118

223.7551

2.7855

223.6928

1.6206

223.6802

0.2

6.5736

223.7709

3.0619

223.6953

1.7239

223.6810

0.4

7.5377

223.8013

3.4816

223.7008

1.9370

223.6872

0.6

7.8369

223.8116

3.8097

223.7068

2.1499

223.6847

0.8

8.1353

223.8222

4.1147

223.7122

2.3330

223.6865

1

8.4366

223.8330

4.3687

223.7170

2.4391

223.6743

The effect of L/a aspect ratio on the natural free vibration frequency of FFFF two-layered grapheme sheets is also evaluated for two different values of LSP. The results are illustrated in Fig. 15-a and Fig. 15-b for the first and the second types of vibration, respectively. According to these figures, the value of natural free vibration frequency decreases when the aspect ratio increases. Moreover, the natural free vibration frequency is negligible for L/a≤1, while for 1≤L/a≤2, the free vibration frequency of the sheet dramatically decreases. For aspect ratios above 2, the slope of the curves reduces until the values of natural free vibration frequency converges. While for the first type of vibration the convergence is observed for L/a≥4, in the second type, convergence occurs at smaller values of L/a ratio, i.e. L/a≥2. It is also interesting that by increasing the aspect ratio, the effects of LSP on the natural free vibration frequency vanishes and the curves approach each other, especially in the second type of vibration. In order to better investigate the effect of the LSP on the natural free vibration frequency, the frequency variation is provided in terms of h/l values for a 10nm square two-layered graphene sheet nanoplate. The results are shown in Fig. 16-a and Fig. 16-b for two types of vibration. According to this figure, it can be concluded that the natural free vibration frequency decreases by increasing h/l ratio. Although the variation patterns in both frequency types are the same, the variation in the first type of vibration, in the absence of the vdW bonds between the layers, is greater than the second type.

46

(a)

(b) Fig. 15. Variation of the non-dimensional natural frequency vs. aspect ratio for two-layered graphene sheets for the: (a) first and, (b) second types of vibration

47

(a)

(b) Fig. 16. Variation of the non-dimensional natural frequency vs. h/l ratio for two-layered graphene sheets for the: (a) first and (b) second types of variation

5. Conclusions In this study, the mechanical buckling and free vibration behavior of SLGSs and MLGSs were evaluated and discussed in details. The critical mechanical buckling loads of the nanoplates mainly depend on the loading pattern and boundary conditions. Thus, SLGSs and MLGSs with different boundary conditions were evaluated under uniaxial compression, biaxial compression and shear loads and the corresponding results were discussed. The small-scale effect was considered by applying the length scale parameter (LSP) to the formulations. According to the results, taking the small-scale effects into account leads to increased mechanical buckling loads and natural free vibration frequencies. Moreover, the effect of length-to-width aspect ratio of the sheets were studied. It was shown that by increasing the aspect ratio above certain values, the critical buckling loads as well as the natural free vibration frequency start to converge and eventually become constant. The buckling load interactions were also evaluated for SLGSs. It was observed that while in the interaction of axial in-plane compressive loads, the effects of the aspect ratio and LSP are significant, the effects of these

48

two parameters were sharply reduced for the interaction between axial and shear forces. In the study of MLGSs, it was shown that the number of buckling and vibration types is equal to the number of the layers, such that their first types are the same as those of SLGSs. Furthermore, in the analysis of MLGSs, a significant change occurs in the values of critical buckling load and natural free vibration frequency when the vdW bonds between the layers appear and affect the total stiffness matrix. In fact, the effects of vdW bonds in MLGSs is so dominant over all other parameters, that by changing the rest of parameters, no significant change is resulted in the values of critical buckling load and the natural free vibration frequency.

49

References: [1]

S. Iijima, C. Brabec, A. Maiti, and J. Bernholc, "Structural flexibility of carbon nanotubes," The Journal of chemical physics, vol. 104, no. 5, pp. 2089-2092, 1996.

[2]

B. Yakobson, M. Campbell, C. Brabec, and J. Bernholc, "High strain rate fracture and C-chain unraveling in carbon nanotubes," Computational Materials Science, vol. 8, no. 4, pp. 341-348, 1997.

[3]

D. Sánchez-Portal, E. Artacho, J. M. Soler, A. Rubio, and P. Ordejón, "Ab initio structural, elastic, and vibrational properties of carbon nanotubes," Physical Review B, vol. 59, no. 19, p. 12678, 1999.

[4]

Y. Li, S. Wang, Q. Wang, and M. Xing, "A comparison study on mechanical properties of polymer composites reinforced by carbon nanotubes and graphene sheet," Composites Part B: Engineering, vol. 133, pp. 35-41, 2018.

[5]

D. Savvas and G. Stefanou, "Determination of random material properties of graphene sheets with different types of defects," Composites Part B: Engineering, vol. 143, pp. 47-54, 2018.

[6]

S. Singh and B. Patel, "Nonlinear elastic properties of graphene sheet using MM3 potential under finite deformation," Composites Part B: Engineering, vol. 136, pp. 8191, 2018.

[7]

X.-H. Hou, X.-J. Xu, J.-M. Meng, Y.-B. Ma, and Z.-C. Deng, "Elastic constants and phonon dispersion relation analysis of graphene sheet with varied Poisson's ratio," Composites Part B: Engineering, vol. 162, pp. 411-424, 2019.

[8]

E. Cosserat and F. Cosserat, "Theory of deformable bodies, No. 11561 in NASA-TTF," NASA Center for Aerospace Information (CASI), 1968.

[9]

A. C. Eringen, "Theory of micropolar plates," Zeitschrift für angewandte Mathematik und Physik ZAMP, vol. 18, no. 1, pp. 12-30, 1967.

50

[10]

A. C. Eringen and E. Suhubi, "Nonlinear theory of simple micro-elastic solids—I," International Journal of Engineering Science, vol. 2, no. 2, pp. 189-203, 1964.

[11]

A. C. Eringen, "Nonlocal polar elastic continua," International journal of engineering science, vol. 10, no. 1, pp. 1-16, 1972.

[12]

R. Mindlin and H. Tiersten, "Effects of couple-stresses in linear elasticity," Archive for Rational Mechanics and analysis, vol. 11, no. 1, pp. 415-448, 1962.

[13]

R. A. Toupin, "Elastic materials with couple-stresses," Archive for Rational Mechanics and Analysis, vol. 11, no. 1, pp. 385-414, 1962.

[14]

R. D. Mindlin, "Second gradient of strain and surface-tension in linear elasticity," International Journal of Solids and Structures, vol. 1, no. 4, pp. 417-438, 1965.

[15]

F. Yang, A. Chong, D. C. C. Lam, and P. Tong, "Couple stress based strain gradient theory for elasticity," International Journal of Solids and Structures, vol. 39, no. 10, pp. 2731-2743, 2002.

[16]

M. H. Hajmohammad, M. S. Zarei, M. Sepehr, and N. Abtahi, "Bending and buckling analysis of functionally graded annular microplate integrated with piezoelectric layers based on layerwise theory using DQM," Aerospace Science and Technology, vol. 79, pp. 679-688, 2018.

[17]

R. M. R. Reddy, W. Karunasena, and W. Lokuge, "Free vibration of functionally graded-GPL reinforced composite plates with different boundary conditions," Aerospace Science and Technology, vol. 78, pp. 147-156, 2018.

[18]

N. Shafiei and M. Kazemi, "Buckling analysis on the bi-dimensional functionally graded porous tapered nano-/micro-scale beams," Aerospace Science and Technology, vol. 66, pp. 1-11, 2017.

51

[19]

S. Kong, S. Zhou, Z. Nie, and K. Wang, "The size-dependent natural frequency of Bernoulli–Euler micro-beams," International Journal of Engineering Science, vol. 46, no. 5, pp. 427-437, 2008.

[20]

H. Ma, X.-L. Gao, and J. Reddy, "A microstructure-dependent Timoshenko beam model based on a modified couple stress theory," Journal of the Mechanics and Physics of Solids, vol. 56, no. 12, pp. 3379-3391, 2008.

[21]

Y. Fu and J. Zhang, "Modeling and analysis of microtubules based on a modified couple stress theory," Physica E: Low-dimensional Systems and Nanostructures, vol. 42, no. 5, pp. 1741-1745, 2010.

[22]

H. Ma, X.-L. Gao, and J. Reddy, "A nonclassical Reddy-Levinson beam model based on a modified couple stress theory," International Journal for Multiscale Computational Engineering, vol. 8, no. 2, 2010.

[23]

F. Z. Jouneghani, P. M. Dashtaki, R. Dimitri, M. Bacciocchi, and F. Tornabene, "Firstorder shear deformation theory for orthotropic doubly-curved shells based on a modified couple stress elasticity," Aerospace Science and Technology, vol. 73, pp. 129147, 2018.

[24]

G. C. Tsiatas, "A new Kirchhoff plate model based on a modified couple stress theory," International Journal of Solids and Structures, vol. 46, no. 13, pp. 2757-2764, 2009.

[25]

L.-L. Ke, Y.-S. Wang, J. Yang, and S. Kitipornchai, "Free vibration of size-dependent Mindlin microplates based on the modified couple stress theory," Journal of Sound and Vibration, vol. 331, no. 1, pp. 94-106, 2012.

[26]

L.-L. Ke and Y.-S. Wang, "Flow-induced vibration and instability of embedded doublewalled carbon nanotubes based on a modified couple stress theory," Physica E: Lowdimensional Systems and Nanostructures, vol. 43, no. 5, pp. 1031-1039, 2011.

52

[27]

E. Jomehzadeh, H. Noori, and A. Saidi, "The size-dependent vibration analysis of micro-plates based on a modified couple stress theory," Physica E: Low-dimensional Systems and Nanostructures, vol. 43, no. 4, pp. 877-883, 2011.

[28]

J. Kim and J. Reddy, "Analytical solutions for bending, vibration, and buckling of FGM plates using a couple stress-based third-order theory," Composite Structures, vol. 103, pp. 86-98, 2013.

[29]

H.-T. Thai and S.-E. Kim, "A size-dependent functionally graded Reddy plate model based on a modified couple stress theory," Composites Part B: Engineering, vol. 45, no. 1, pp. 1636-1645, 2013.

[30]

L. Yin, Q. Qian, L. Wang, and W. Xia, "Vibration analysis of microscale plates based on modified couple stress theory," Acta Mechanica Solida Sinica, vol. 23, no. 5, pp. 386-393, 2010.

[31]

M. Mirsalehi, M. Azhari, and H. Amoushahi, "Stability of thin FGM microplate subjected to mechanical and thermal loading based on the modified couple stress theory and spline finite strip method," Aerospace Science and Technology, vol. 47, pp. 356366, 2015.

[32]

R. P. Shimpi, "Refined plate theory and its variants," AIAA journal, vol. 40, no. 1, pp. 137-146, 2002.

[33]

S. Narendar, "Buckling analysis of micro-/nano-scale plates based on two-variable refined plate theory incorporating nonlocal scale effects," Composite Structures, vol. 93, no. 12, pp. 3093-3103, 2011.

[34]

S. Narendar and S. Gopalakrishnan, "Scale effects on buckling analysis of orthotropic nanoplates based on nonlocal two-variable refined plate theory," Acta Mechanica, vol. 223, no. 2, pp. 395-413, 2012.

53

[35]

N. Satish, S. Narendar, and S. Gopalakrishnan, "Thermal vibration analysis of orthotropic nanoplates based on nonlocal continuum mechanics," Physica E: Lowdimensional Systems and Nanostructures, vol. 44, no. 9, pp. 1950-1962, 2012.

[36]

Y. Wang and D. Wu, "Free vibration of functionally graded porous cylindrical shell using a sinusoidal shear deformation theory," Aerospace Science and Technology, vol. 66, pp. 83-91, 2017.

[37]

A. Tounsi, M. S. A. Houari, and S. Benyoucef, "A refined trigonometric shear deformation theory for thermoelastic bending of functionally graded sandwich plates," Aerospace science and technology, vol. 24, no. 1, pp. 209-220, 2013.

[38]

M. Ameur, A. Tounsi, I. Mechab, and A. A. El Bedia, "A new trigonometric shear deformation theory for bending analysis of functionally graded plates resting on elastic foundations," KSCE Journal of Civil Engineering, vol. 15, no. 8, pp. 1405-1414, 2011.

[39]

S. Sarrami-Foroushani and M. Azhari, "Nonlocal buckling and vibration analysis of thick rectangular nanoplates using finite strip method based on refined plate theory," Acta Mechanica, vol. 227, no. 3, pp. 721-742, 2016.

[40]

M. Sobhy, "Generalized two-variable plate theory for multi-layered graphene sheets with arbitrary boundary conditions," Acta Mechanica, vol. 225, no. 9, pp. 2521-2538, 2014.

[41]

H.-T. Thai and S.-E. Kim, "Levy-type solution for buckling analysis of orthotropic plates based on two variable refined plate theory," Composite Structures, vol. 93, no. 7, pp. 1738-1746, 2011.

[42]

H.-T. Thai and S.-E. Kim, "Analytical solution of a two variable refined plate theory for bending analysis of orthotropic Levy-type plates," International Journal of Mechanical Sciences, vol. 54, no. 1, pp. 269-276, 2012.

54

[43]

H.-T. Thai and S.-E. Kim, "Levy-type solution for free vibration analysis of orthotropic plates based on two variable refined plate theory," Applied Mathematical Modelling, vol. 36, no. 8, pp. 3870-3882, 2012.

[44]

H.-T. Thai and S.-E. Kim, "Free vibration of laminated composite plates using two variable refined plate theory," International Journal of Mechanical Sciences, vol. 52, no. 4, pp. 626-633, 2010.

[45]

S.-E. Kim, H.-T. Thai, and J. Lee, "A two variable refined plate theory for laminated composite plates," Composite Structures, vol. 89, no. 2, pp. 197-205, 2009.

[46]

M. Sobhy, "Hygrothermal vibration of orthotropic double-layered graphene sheets embedded in an elastic medium using the two-variable plate theory," Applied Mathematical Modelling, vol. 40, no. 1, pp. 85-99, 2016.

[47]

A. C. Ugural, Stresses in beams, plates, and shells. CRC press, 2009.

[48]

C. Y. Wang and C. Wang, Structural vibration: exact solutions for strings, membranes, beams, and plates. CRC Press, 2016.

[49]

X. He, S. Kitipornchai, and K. Liew, "Resonance analysis of multi-layered graphene sheets used as nanoscale resonators," Nanotechnology, vol. 16, no. 10, p. 2086, 2005.

[50]

H. X. Nguyen, T. N. Nguyen, M. Abdel-Wahab, S. Bordas, H. Nguyen-Xuan, and T. P. Vo, "A refined quasi-3D isogeometric analysis for functionally graded microplates based on the modified couple stress theory," Computer Methods in Applied Mechanics and Engineering, vol. 313, pp. 904-940, 2017.

[51]

S. Sarrami-Foroushani and M. Azhari, "Nonlocal vibration and buckling analysis of single and multi-layered graphene sheets using finite strip method including van der Waals effects," Physica E: Low-dimensional Systems and Nanostructures, vol. 57, pp. 83-95, 2014.

55

[52]

R. Plank and W. Wittrick, "Buckling under combined loading of thin, flat‐walled structures by a complex finite strip method," International Journal for Numerical Methods in Engineering, vol. 8, no. 2, pp. 323-339, 1974.

56