Nonlocal plate model for nonlinear bending of bilayer graphene sheets subjected to transverse loads in thermal environments

Composite Structures 98 (2013) 294–302

Contents lists available at SciVerse ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Nonlocal plate model for nonlinear bending of bilayer graphene sheets subjected to transverse loads in thermal environments Yu-Mou Xu a, Hui-Shen Shen a,b,⇑, Chen-Li Zhang a a b

School of Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200030, People’s Republic of China State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200030, People’s Republic of China

a r t i c l e

i n f o

Article history: Available online 23 November 2012 Keywords: Graphene sheet Nonlinear bending Nonlocal plate model

a b s t r a c t This paper investigates the nonlinear bending behavior of a bilayer rectangular graphene sheet subjected to a transverse uniform load in thermal environments. The bilayer graphene sheet (BLGS) is modeled as a nonlocal double-layered plate which contains small scale effect and van der Waals interaction forces. The geometric nonlinearity in the von Kármán sense is adopted. The thermal effects are included and the material properties are assumed to be size-dependent and temperature-dependent, and are obtained from molecular dynamics (MD) simulations. The small scale parameter e0a is estimated by matching the deflections of graphene sheets observed from the MD simulation results with the numerical results obtained from the nonlocal plate model. The numerical results show that the stacking sequence has a moderate effect, whereas the temperature change as well as the aspect ratio has a significant effect on the nonlinear bending behavior of BLGSs. The results reveal that the small scale parameter reduces the static large deflections of BLGSs, and the small scale effect also plays an important role in the nonlinear bending of BLGSs. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Recently, a member of carbon family known as graphene sheet (GS) has drawn considerable attention [1]. Graphene is a twodimensional atomic crystal which consists of carbon atoms arranged in a hexagonal lattice. Owing to their exceptional thermal, electrical and mechanical properties [2–5], GSs can be used in a variety of ways such as flexible cancer sensors [6], mechanical resonators [7], conductive electrodes for solar cells [8] and novel composite materials [9,10]. As a consequence, the analysis of GSs is a fundamental issue in the study of nanocomposites. Much attention has been devoted to the linear and nonlinear bending behaviors of both suspended and supported graphene sheets. Poot et al. [11] measured the nano-mechanical properties of graphene sheets that were suspended over circular holes by using an atomic force microscope (AFM). They found that graphene sheets can sustain very large bending and stretching prior to the occurrence of fracture. This means the large deflection analysis of GSs needs to be performed. Duan and Wang [12] investigated the deformation of a single layer circular graphene sheet under a central point load by molecular mechanics simulations. They concluded that the nonlinear bending and stretching behaviors of ⇑ Corresponding author at: School of Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200030, People’s Republic of China. E-mail address: [email protected] (H.-S. Shen). 0263-8223/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2012.10.041

single layer graphene sheets (SLGSs) can be well described by the von Kármán plate theory when the material properties are properly chosen. Hemmasizadeh et al. [13] also studied the large deflection of circular SLGSs loaded at the centre based on an equivalent continuum model. Scarpa et al. [14] studied both circular and rectangular graphene sheets subjected to point loading based on a special equivalent atomistic-continuum model. They found that the rectangular SLGSs show a different distribution of the geometrical and mechanical properties, compared to the circular configuration. On the other hand, Wang [15] re-examined the bending rigidity of rectangular graphene sheets subjected to a central point load by using molecular mechanics simulations. Shen and Wu [16] investigated the effects of interlayer shear and ripple on the bending properties of multi-layer graphene sheets (MLGSs) through molecular dynamics (MD) simulations. It is worth noting that most studies cited above are for SLGSs except for [16], and the material properties of GSs are usually assumed to be independent of temperature. Although the classical continuum models are efficient in GSs bending analysis through their relatively simple formulas [12– 14], the mismatch in deformation profiles between MD calculations and continuum mechanics methods is about 89% at a central deflection of SLGSs [17]. It has been reported that the small scale effect must play an important role in the nanoscale structures, but this small scale effect has been ignored when classical local continuum theory was adopted. A number of elasticity theories have been proposed to account for the size effect in the micro/nano-scale

295

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

structures, e.g. nonlocal elasticity [18], modified strain gradient elasticity [19], and modified couple stress theory [20]. Among those, the nonlocal elasticity contains two additional material length scale parameters e0 and a, the modified strain gradient elasticity contains three additional material length scale parameters li (i=0,1,2) and the modified couple stress theory contains only one additional material length scale parameter l, which is the special case of the modified strain gradient elasticity when l0=l1=0 and l2=l. These theories allow for the integration of small scale effects into classical continuum models. However, Sun et al. [21] found that there exists a noticeable inconsistency between atomistic simulation and the strain gradient elasticity solution for the bending of micro/nano-scale structures. In contrast, the nonlocal elasticity of Eringen [18] is simple and convenient to apply for predicting the nonlinear response of the micro/nano-scale structures. The nonlinear bending behaviors of SLGSs subjected to a transverse uniform load in thermal environments were recently studied by Shen et al. [22] based on a nonlocal orthotropic plate model. In the present work, we focus our attention on the nonlinear bending behavior of bilayer graphene sheets (BLGSs) subjected to a transverse uniform load in thermal environments. The BLGS is modeled as a nonlocal double-layered plate which contains small scale effect and van der Waals interaction forces. The nonlinear bending analysis is based on thin plate theory with a von Kármán-type of kinematic nonlinearity. The thermal effects are also included and the material properties are assumed to be anisotropic, sizedependent and temperature-dependent, and are obtained from MD simulations. The small scale parameter e0a is estimated by matching the deflections of graphene sheets observed from the MD simulation results with the numerical results obtained from the nonlocal plate model. The numerical illustrations show the full nonlinear bending response of BLGSs subjected to a transverse uniform load under different sets of environmental conditions. 2. Nonlocal plate model for BLGSs Consider a BLGS modeled as a double-layered plate which is exposed to elevated temperature and subjected to a transverse uniform load q. The top and bottom sheets are assumed to have the same length Lx, width Ly, but may have different thickness hI and hII. Depending upon the chirality of GSs, the top and bottom sheets are assumed to be either armchair or zigzag type. The definitions of armchair and zigzag are the same as those used in carbon nanotubes [23]. As usual, the coordinate system has its origin at the corner of the top sheet, as shown in Fig. 1. Let U; V and W be the plate displacements parallel to a right-hand set of axes (X, Y, Z), where X is longitudinal and Z is perpendicular to the plate. Let FðX; YÞ be the stress function for the stress resultants defined by N x ¼ F;YY ; N y ¼ F;XX and N xy ¼ F;XY , where a comma denotes partial differentiation with respect to the corresponding coordinates. It has been reported that the material properties at nanoscales are size dependent [23]. In order to incorporate the small scale effect, continuum plate models need to be refined. This may be accomplished by using the nonlocal continuum theory. In the theory of nonlocal elasticity, the constitutive relations of nonlocal elasticity for 3D problems are expressed as [18]

  1  s2 L2x r2 rij ¼ C ijkl ekl

ð1Þ

where s = e0a/Lx, rij and eij are the stress and strain tensors, and Cijkl is the elastic modulus tensor of classical isotropic elasticity, e0 is a material constant, and a and Lx are the internal and external characteristic lengths, respectively. The distinct difference between the classical and nonlocal elasticity theories lies in the presence of small scale parameter e0a in the nonlocal theory.

X Y

Z

(a) q

X

Z

(b)

Fig. 1. A bilayer graphene sheet under a transverse uniform pressure.

For carbon nanotubes the characteristic length a may be taken as the length of the C–C bond, i.e. a = 0.142 nm, or taken as other material properties. It is a better way to use e0a as a single scale coefficient that captures the small scale effect on the response of structures in nanosize [22]. Applying Eq. (1) to classical plate theory, the Kármán-type nonlinear differential equations for the BLGS, including small scale effects, van der Waals interaction forces and thermal effects, have readily been derived and can be expressed in terms of a stress function F, and a transverse displacement W. They are

  ðIÞ ðIÞ e L 11 ðW I Þ  e L 14 ððMT ÞI Þ ¼ 1  s2 L2x r2 ½e LðW I ; F I Þ þ C 1 ðW II  W I Þ þ C 3 ðW II  W I Þ3 þ q 1 ðIÞ ðIÞ e L 21 ðF I Þ  e LðW I ; W I Þ L 23 ððNT ÞI Þ ¼  e 2   ðIIÞ ðIIÞ e L 11 ðW II Þ  e LðW II ; F II Þ  C 1 ðW II L 14 ððM T ÞII Þ ¼ 1  s2 L2x r2 ½e  W I Þ  C 3 ðW II  W I Þ3  1 ðIIÞ e LðW II ; W II Þ L 23 ððN T ÞII Þ ¼  e L 21 ðF II Þ  e 2 where (with J = I, II)

h    @4   i @4   @4 ðJÞ e þ D22 J 4 L 11 ðÞ ¼ D11 J 4 þ 2 D12 J þ 2 D66 J 2 2 @X @X @Y @Y 4 h     i @4    @4 @ ðJÞ    e L 21 ðÞ ¼ ðA22 ÞJ 4 þ 2 A12 J þ A66 J þ A11 J 4 @X @X 2 @Y 2 @Y 2  2  2     @ @ @ ðJÞ e M Txy þ 2 MTy L 14 ððMT ÞJ Þ ¼ 2 M Tx J þ 2 J J @X@Y @X @Y

ð2Þ ð3Þ

ð4Þ ð5Þ

296

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

        @2 @2 1 ðJÞ e L 23 ððNT ÞJ Þ ¼ 2 A12 NTx J þ A22 NTy þ2  A66 NTxy J J 2 @X@Y @X  2      @ þ 2 A11 NTx J þ A12 NTy J @Y @2 @2 @2 @2 @2 @2 e LðÞ ¼ 2 2 þ 2 2 @X@Y @X@Y @X @Y @Y @X 2

r2 ðÞ ¼

@2 @X

2

þ

pffiffiffi!2 4 3 C1 ¼ 9d " # 24e0 13p r0 14 1 7p r0 8 1  2  3 d 3 d ðz1  z2 Þ6 r0 ðz1  z2 Þ12 pffiffiffi!2 4 3 9d " # r 16 r 10 4e0 1 1 0 0  4 780p  126p d d r0 ðz1  z2 Þ14 ðz1  z2 Þ8

2

Q 11

6 7 6 6 Ay ðTÞ 7 ¼ 6 6 Q 12 5 4 4 Axy ðTÞ

J

c4 þ s4 2c2 s2

c2 s2 c4

4c2 s2 4c2 s2

72 3 7 Q 11 7 76 Q 7 76 12 7 7 6 3 3 3 2 2 7 cs  c s cs 2csðc  s Þ 74 Q 22 5 7 c3 s  cs3 c3 s 2csðc2  s2 Þ 7 5 Q 66 2c2 s2 c2 s2 ðc2  s2 Þ2

ð10Þ

Q 16

Q 12 Q 22 Q 26

Q 16

32

c2

s2

E22 ðTÞ ; ð1  m12 m21 Þ

Q 22 ¼

Q 12

Q 66 ¼ G12 ðTÞ

ð11aÞ

and

c ¼ cos h;

s ¼ sin h

ð11bÞ

where h is the skew angle with respect to the plate X axis. All four edges are assumed to be simply supported with no inplane displacements, so that the boundary conditions are X = 0, Lx;

ð12aÞ

Ly

Z

Lx

@2W J @X

2

 D12

@U dXdY ¼ @X

Z

Ly

@2W J @Y

Z

Lx

2

8" <

  þ MTx J ¼ 0 ðJ ¼ I; IIÞ

A11

@2 F J

þ A12

@2 FJ

# 

1 @W J 2 @X

ð12bÞ !2

@Y 2 @X 2      dXdY ¼ 0 ðJ ¼ I; IIÞ  A11 NTx J þ A12 N Ty

0

0

:

0

J

ð8Þ

ð12cÞ

Y = 0, Ly;

W J ¼ 0 ðJ ¼ I; IIÞ My ¼ D12 Z 0

Lx

Z 0

Ly

@2W J @X

2

ð12dÞ

 D22

@2W J @Y

2

  þ M Ty ¼ 0 ðJ ¼ I; IIÞ J

8" # !2 < @2 FJ @2 FJ 1 @W J  A22 þ A12 2 2 2 @Y @X @Y 0 0 :        A12 N Tx J þ A22 N Ty dYdX ðJ ¼ I; IIÞ

@V dYdX ¼ @Y

Z

Lx

Z

ð12eÞ

Ly

J

ð12fÞ

where M x and M y are, respectively, the bending moments per unit width and per unit length of the sheet. It is worth noting that the thermal forces N T and the thermal moments MT are also included in the boundary condition of Eq. (12). h i In the above equations the reduced stiffness matrices Aij and h i Dij ði; j ¼ 1; 2; 6Þ are functions of temperature, in the present case, determined through relationship

J

3

4c2 s2



3

2 3 Ax ðTÞ Z þhJ =2 7 6 7 6 Ay ðTÞ 7 ð1; ZÞDTdZ M Ty 7 7 ¼ 4 5 5 hJ =2 T Axy ðTÞ J M xy

Ax ðTÞ

s4

E11 ðTÞ ; ð1  m12 m21 Þ m21 E11 ðTÞ ¼ ; ð1  m12 m21 Þ

0

where DT = T  T0 is temperature rise from some reference temperature T0 at which there are no thermal strains, and

2

3

2c2 s2

Q 11 ¼

Z

where zi ¼ Z i =d; d ¼ 0:142 nm is the C–C bond length, e0 is the depth of the potential, and r0 is a parameter that is determined by the equilibrium distance. For each sheet, four independent material constants (Young’s moduli E11, E22, shear modulus G12, and Poisson’s ratio m12) are needed to describe their mechanical properties. It is assumed that the effective Young’s moduli E11 and E22, shear modulus G12 and thermal expansion coefficients a11 and a22 of GSs are temperature-dependent, whereas Poisson’s ratio m12 depends weakly on temperature change and is assumed to be a constant. The temperature field considered is assumed to be a uniform distribution over the sheet surface and varied in the thickness direction. The thermal forces N T and moments M T caused by elevated temperature are defined by

N Txy

7 6 6 6 Q 12 7 6 c2 s2 7 6 6 6 Q 7 6 s4 6 22 7 6 7¼6 6 6 Q 16 7 6 c3 s 7 6 6 7 6 3 6 4 Q 26 5 4 cs Q 66 c 2 s2

Mx ¼ D11

M Tx

c4

W J ¼ 0 ðJ ¼ I; IIÞ

ð7bÞ

NTx 6 6 NT 6 y 4

2

ð7aÞ

and

2

3

ð6Þ

2

The geometric nonlinearity in the von Kármán sense is given in terms of e L() in Eqs. (2)–(5). Small scale parameter e0a is included in Eqs. (2) and (4). Eqs. (2)–(5) are coupled and should be solved simultaneously. The nonlocal plate model described by Eqs. (2)– (5) reduces to the local plate model when the small scale parameter e0a vanishes. In Eqs. (2) and (4), C 1 and C 3 are, respectively, the linear and nonlinear van der Waals interaction coefficients and can be defined by [24]

C3 ¼

Q 11

where

@2 @Y

2

3

A ¼ A1 ;

D ¼ D

ð13Þ

where Aij and Dij are the stiffnesses of each sheet, defined by

"

76 7 a11 ðTÞ 7 s2 c2 7 Q 26 7 6 5 4 5 a22 ðTÞ Q 66 J 2cs 2cs

# ð9Þ J

in which a11 and a22 are the thermal expansion coefficients in the longitudinal and transverse directions, and Q ij are the transformed elastic constants, defined by

ðAij ; Dij ÞJ ¼

Z

þhJ =2

hJ =2

ðQ ij ðTÞÞJ ð1; Z 2 ÞdZ ði; j ¼ 1; 2; 6; and J ¼ I; IIÞ

ð14Þ

Note that for the armchair graphene sheet h = 0°, and for the zigzag graphene sheet h = 90°, and therefore, A16 ; A26 ; D16 and D26 are zero-valued, and for both of them, no extension-flexural coupling exists, and Bij are all zero-valued.

297

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

Having developed the theory, we are in a position to solve Eqs. (2)–(5) with boundary conditions (12). Before proceeding, it is convenient first to define the following dimensionless quantities (with J = I, II) WJ        i1=4 A11 J A22 J       "   #1=2 D12 J þ 2 D66 J D22 J FJ ðJÞ ðJÞ   F J ¼       1=2 ; c12 ¼ ; c14 ¼    D11 J D11 J ½ D11 J D22 J        "   #1=2    A12 J þ A66 J =2 A11 J A12 J ðJÞ ðJÞ    cðJÞ ; c24 ¼    ; c5 ¼     22 ¼ A22 J A22 J A22 J      T T Ax ; Ay   L2 qL4 J J ðJÞ x cðJÞ ; kq ¼         x        1=4 T1 ; cT2 ¼ 2 h 4 D p D  D  i1=2 p D D 11 I 11 I 22 I A11 I A22 I 11 J 22 J x¼p

X Y Lx ; y ¼ p ; b ¼ ; W J ¼ h Lx Ly Ly

ðM x ; M y ÞJ ¼

L2



D11 J

D22 J

ðM x ; M y ÞJ

x p2 D  hD  D  A  A  i1=4 11 J 22 J 11 J 11 J 22 J

4

h i h       i1=2 1 Lx ðJÞ ðJÞ C1 ; C3 ¼    C 1 ; C 3 D11 J D22 J A11 J A22 J ; D11 J p

ð15Þ

in which ATx and ATy are defined by

"

ATx

#

DT ¼ 

ATy

Z

hJ =2

hJ =2

J



Ax Ay



DTdZ

ð16Þ

#



2 @2F J @2FJ 1 @W J   c c 5 24 2 @y2 @x2 @x 0 0    2 ðJÞ ðJÞ ðJÞ ðJÞ þ c24 cT1  c5 cT2 DT dxdy ¼ 0

Z pZ p

3. Theoretical solutions

("

c224 b2

y = 0, p;

WJ ¼ 0

ð22dÞ

ðM y ÞJ ¼ 0

ð22eÞ

Z p Z p (" 2 @ F 0

J @x2

0

1 ðIÞ L21 ðF I Þ ¼  c24 b2 LðW I ; W I Þ 2 h ðIIÞ ðIIÞ L11 ðW II Þ ¼ð1  s2 p2 r2 Þ c14 b2 LðW II ; F II Þ  C 1 ðW II  W I Þ i ðIIÞ C 3 ðW II  W I Þ3 1 ðIIÞ L21 ðF II Þ ¼  c24 b2 LðW II ; W II Þ 2

X

Wðx; y; eÞ ¼

ðJÞ

 2 @ 4 @4 @4 ðJÞ ðJÞ þ 2c12 b2 2 2 þ c14 b4 4 4 @x @x @y @y

ðJÞ

@4 @4 @4 ðJÞ ðJÞ þ 2c22 b2 2 2 þ ðc24 Þ2 b4 4 4 @x @x @y @y

L21 ðÞ ¼ LðÞ ¼

@2 @2 þ b2 2 2 @x @y

X

X

ej f ðjÞ ðx; yÞ;

j¼0

ej kj

ð23Þ

j¼1

where e is a small perturbation parameter and the first term of w(j)(x, y) is assumed to have the form ð1Þ

ð17Þ

ð1Þ

ð1Þ wII ðx; yÞ

¼

ð1Þ a11

ð24aÞ

sin mx sin ny

ð24bÞ ð1Þ ð1Þ a11 =A11

ð18Þ

ð19Þ ð20Þ

where (m, n) is the deflected mode, ¼ v1 the details of which may be found in Appendix A. Substituting Eq. (23) into Eqs. (17)–(20) and collecting the terms of the same order of e, a set of perturbation equations is obtained. By using Eq. (24) to solve these perturbation equations of each order, the amplitudes of the terms w(j)(x, y) and f(j)(x, y) are determined step by step, and kj can be determined by the Galerkin procedure. As a result, up to 3rd-order asymptotic solutions can be obtained as

h i ð1Þ W I ¼ e A11 sinmxsinny h i ð3Þ ð3Þ ð3Þ þ e3 A13 sinmxsin3nyþA31 sin3mxsinnyþA33 sin3mxsin3ny þOðe4 Þ ð25Þ

ð21Þ

2 2 ð0Þ x ð0Þ y  b00 F I ¼ B00 2 2   2 2 ð2Þ x ð2Þ y ð2Þ ð2Þ þ e2 B00  b00 þ B20 cos 2mx þ B02 cos 2ny 2 2

þ Oðe4 Þ

The boundary conditions of Eq. (12) become

x ¼ 0;

ð22fÞ

wI ðx; yÞ ¼ A11 sin mx sin ny

@2 @2 @2 @2 @2 @2 2 þ 2 2 2 2 @x @y @x@y @x@y @y @x

r2 ðÞ ¼

ej wðjÞ ðx; yÞ; Fðx; y; eÞ ¼

j¼1

where (with J = I, II)

L11 ðÞ ¼

#

2 @2FJ 1 2 @W J c b  2 24 @y2 @y

By virtue of the fact that DT is assumed to be uniform, the thermal coupling in Eqs. (17)–(20) vanishes, but terms in DT intervene in Eqs. (22c) and (22f). Applying Eqs. (17)-(22), the nonlinear bending response of BLGSs subjected to a transverse uniform load in thermal environments is now determined by means of a two step perturbation technique, where the small perturbation parameter has no physical meaning at the first step, and is then replaced by a dimensionless deflection at the second step. The essence of this procedure, in the present case, is to assume that

kq ¼

h ðIÞ ðIÞ L11 ðW I Þ ¼ð1  s2 p2 r2 Þ c14 b2 LðW I ; F I Þ þ C 1 ðW II  W I Þ i ðIÞ þ C 3 ðW II  W I Þ3 þ kq

 c5 b2

h i o ðJÞ ðJÞ ðJÞ þ cT2  c5 cT1 DT dydx ¼ 0

J

where Ax and Ay are given in detail in Eq. (9). The nonlinear Eqs. (2)–(5) may then be written in dimensionless form as

ð22cÞ

p;

WJ ¼ 0

ð22aÞ

ðMx ÞJ ¼ 0

ð22bÞ

ð26Þ

h i ð1Þ W II ¼e a11 sin mx sin ny h ð3Þ ð3Þ ð3Þ þ e3 a11 sin mx sin ny þ a13 sin mx sin 3ny þ a31 sin 3mx sin ny i ð3Þ ð27Þ þa33 sin 3mx sin 3ny þ Oðe4 Þ

298

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

2 2 ~ð0Þ x e ð0Þ y  b F II ¼  B 00 00 2 2   2 2 ð2Þ y 2 ~ð2Þ x þ bð2Þ cos 2mx þ bð2Þ cos 2ny e þ e  B 00 b 20 02 00 2 2

þ Oðe4 Þ

ð28Þ

and

kq ¼ ek1 þ e3 k3 þ Oðe4 Þ

ð29Þ

All coefficients in Eqs. (25)–(29) are related and can be written ð1Þ as functions of A11 , so that Eqs. (25) and (29) can be rewritten as

   3 ð1Þ ð1Þ W I ¼ W ð1Þ ðx; yÞ A11 e þ W ð3Þ ðx; yÞ A11 e þ   

ð30Þ

and

   3 ð1Þ ð1Þ kq ¼ kqð1Þ A11 e þ kqð3Þ A11 e þ   

ð31Þ

From Eqs. (30) and (31) the load–central deflection relationship can be written as

q L4 ð1Þ W m  0 x ¼ AW hI D11 I hI

! þ

ð3Þ AW

Wm hI

!3 þ 

ð32Þ

or

Wm qL4 ð1Þ ¼ BW   x hI D11 I hI ðiÞ

! ð3Þ

 BW

qL4   x D11 I hI

!3 þ 

ð33Þ

ðiÞ

It is noted that AW and BW ði ¼ 1; 3; . . .Þ are related to the material properties and are all functions of temperature with details being given in Appendix A. 4. Numerical results and discussion Numerical results are presented in this section for BLGSs subjected to a transverse uniform load in thermal environments. As mentioned previously [22] common mechanical properties of materials such as Young’s modulus, shear modulus, Poisson’s ratio, as well as the effective thickness of the sheet, are all concepts of continuum mechanics, and should be determined firstly. However, a large variation of Young’s modulus E and Poisson’s ratio m, as well as effective thickness h was obtained and reported in the open literature [25]. The wide dispersion of the mechanical properties of SLGSs can be attributed principally to the uncertainty associated to the thickness of these nanostructures. For the majority of models used, the assumed thickness of the graphene layer is 0.34 nm. The 0.34 nm value provides in-plane Young’s modulus of the order of 1 TPa [4]. In contrast, Hemmasizadeh et al. [13] used a mixed MD–continuum mechanics model to obtain E = 0.939 TPa and m = 0.19 with h = 0.1317 nm. Huang et al. [23] used the second generation Brenner potentials to calculate the in-plane Young’s modulus, Poisson’s ratio and thickness of SLGSs and to obtain E = 2.99 TPa and m = 0.397 with h = 0.0811 nm. Their value of Young’s modulus is found to be about 3 times as large as that of Hemmasizadeh et al. [13]. It has been reported that the material properties of SLGSs are anisotropic, chirality- and size-dependent and temperature-dependent [4,25]. Therefore, all material properties and effective thickness of a SLGS need to be carefully determined, otherwise the results may be incorrect. The MD simulations are first carried out for SLGSs, in which Newtonian equations of motion governed by interatomic interactions are solved numerically to determine the trajectories of a large number of atoms. A Velocity-Verlet algorithm is used to integrate the equations of motion and a basic time step of 0.5 fs is employed to guarantee good conservation of energy. In our MD simulations

we use the LAMMPS code, in which the adaptive intermolecular reactive empirical bond order potential (AIREBO) is adopted to describe the short-range covalent C–C interactions, the long-range van der Waals interaction (LJ terms) and torsion interactions [26]. System temperature conversion is carried out by the Nose– Hoover feedback thermostat [27]. To begin the MD simulation, the SLGS is initially optimized and freely relaxed to reach the minimum energy configuration. The deformation of the SLGS is carried out in a quasi-static way by gradually increasing the applied load in a small increment and allowing the SLGS to relax fully until the next equilibrium configuration is reached. The perfect SLGSs subjected to uniaxial tensile force and tangential force are simulated under temperature varying from 300 to 700 K. The bending tests are then carried out and the effective thickness of SLGSs can be determined uniquely, as reported in [25]. Eight types of armchair and zigzag SLGSs with three different values of aspect ratio are considered. From our MD simulation results the material properties of SLGSs under thermal environmental conditions T = 300, 500, 700 K are obtained numerically. Typical results are listed in Table 1. It is noted that the effective thickness for the armchair graphene sheets is 0.0744–0.1132 nm, while for the zigzag graphene sheets is 0.0801–0.121 nm. After getting the material properties and effective thickness of SLGSs, the bending tests are then carried out for BLGSs. Four types of BLGSs are configurated. For Type AA, the top and bottom sheets are both armchair. Similarly, for type ZZ, the top and bottom sheets are both zigzag. For Type AZ, the top sheet is armchair while the bottom sheet is zigzag. For Type ZA, the stacking sequence is inversed, i.e. the top sheet is zigzag while the bottom sheet is armchair. From our MD simulation results the large deflections of these four types of BLGSs under thermal environmental conditions T = 300, 500, 700 K are obtained numerically. The large deflections of these four types of BLGSs with or without nonlinear van der Waals interaction coefficient C 3 are listed and compared in Table 2. By taking e0 = 2.84 meV and r0 = 0.34 nm [26], from Eq. (7) we have the van der Waals interaction constants C 1 ¼ 99:9943 GPa=nm and C 3 ¼ 47142:62 GPa=nm3 . It can be seen that the nonlinear van der Waals interaction coefficient C 3 has almost no effect on the nonlinear bending behavior of BLGSs, and it may be reasonably neglected, as previously reported in [28–30]. Fig. 2 presents load–deflection curves of four types of BLGSs subjected to a transverse uniform pressure at 300 K. The graphene sheet aspect ratio is taken to be about b = 4.5, which consists of Armchair sheet IV and Zigzag sheet IV, as defined in Table 1. It can be seen that the load–deflection curves of AA and ZZ BLGSs are close. It can also be seen that the ZA BLGS will have the highest, while the AZ BLGS will have the lowest load–deflection curves among the four. Hence, only ZA and AZ BLGSs will be discussed in the following examples. Fig. 3 presents the load–deflection curves of ZA and AZ (Lx = 21.9777 nm, Ly = 4.9126 nm) BLGSs subjected to a transverse uniform load under thermal environmental condition T = 300 K. The small scale parameter e0a = 0.0, 0.4, 0.6 and 0.8 nm. Note that here e0a = 0 represents the local plate model. It is found that the load–deflection curves are decreased with increase in the small scale parameter e0a. In other words, the effect of the small scale on the nonlinear bending response becomes more significant as the parameter e0a becomes larger. However, there are no experiments conducted to determine the value of e0a for GSs. In present study, we give the estimation of parameter e0a by matching the deflections of graphene sheets observed from the MD simulation results with the numerical results obtained from the nonlocal plate model, as shown in the next examples. Fig. 4 presents the load–deflection curves for the same two BLGSs subjected to a transverse uniform load under thermal

299

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302 Table 1 Material properties of SLGSs in thermal environments. G12 (TPa)

a11 (106/K)

a22 (106/K)

Armchair sheet I: Lx = 12.1248 nm, Ly = 4.8933 nm, h = 0.1132 nm, m12 = 0.188 300 2.7833 2.854 500 2.765 2.8269 700 2.7385 2.7827

1.0868 1.1219 1.0866

1.705 1.914 2.004

1.977 2.172 2.376

Armchair sheet II: Lx = 15.1039 nm, Ly = 4.893 nm, h = 0.1023 nm, m12 = 0.209 300 3.011 3.0697 500 2.9912 3.0205 700 2.9619 3.001

1.2122 1.2414 1.3001

1.767 1.922 2.159

2.288 2.294 2.884

Armchair sheet III: Lx = 18.0764 nm, Ly = 4.8909 nm, h = 0.084 nm, m12 = 0.229 300 3.5694 3.6289 500 3.5119 3.5833 700 3.4762 3.5595

1.4397 1.5476 1.5119

2.287 2.423 2.45

2.354 2.496 2.69

Armchair sheet IV: Lx = 21.9131 nm, Ly = 4.895 nm, h = 0.0744 nm, m12 = 0.192 300 3.8868 3.9674 500 3.7903 3.9113 700 3.7634 3.8306

1.7753 1.6532 1.6801

1.583 1.735 2.022

2.164 2.259 2.34

Zigzag sheet I: Lx = 11.7698 nm, Ly = 4.8993 nm, h = 0.121 nm, m12 = 0.222 300 2.62 2.562 500 2.562 2.5207 700 2.5207 2.5041

1.0827 1.0909 1.0248

1.985 2.28 2.382

1.568 1.992 2.191

Zigzag sheet II: Lx = 14.9609 nm, Ly = 4.8974 nm, h = 0.1058 nm, m12 = 0.233 300 2.9294 2.8821 500 2.8639 2.7883 700 2.8261 2.7694

1.1906 1.2571 1.3516

2.118 2.529 2.771

1.873 2.007 2.175

Zigzag sheet III: Lx = 17.9063 nm, Ly = 4.8957 nm, h = 0.0919 nm, m12 = 0.207 300 3.2439 3.2004 500 3.2318 3.1012 700 3.1556 3.0577

1.328 1.4037 1.4037

2.013 2.438 2.484

1.813 2.22 2.386

Zigzag sheet IV: Lx = 21.8391 nm, Ly = 4.8976 nm, h = 0.0801 nm, m12 = 0.205 300 3.5592 3.5217 500 3.5081 3.4832 700 3.4956 3.3708

1.5236 1.5231 1.6854

1.833 2.016 2.252

1.316 1.819 2.111

T (K)

E11 (TPa)

E22 (TPa)

Table 2 Large deflections of BLGSs with or without nonlinear van der Waals interaction coefficient. C 3 ðGPa=nm3 Þ

AA-I AA-II AA-III AA-IV ZZ-I ZZ-II ZZ-III ZZ-IV AZ-I AZ-II AZ-III AZ-IV ZA-I ZA-II ZA-III ZA-IV

W (nm) q0 = 30 MPa

q0 = 45 MPa

q0 = 60 MPa

47142.62 0 47142.62 0 47142.62 0 47142.62 0

0.10892345294447 0.10892345296955 0.11415142509019 0.11415142509833 0.11921297644729 0.11921297644868 0.12273948099999 0.12273948100039

0.12939282194457 0.12939282197580 0.13418973698576 0.13418973699564 0.13869299465245 0.13869299465410 0.14213602455788 0.14213602455835

0.14534878214759 0.14534878218356 0.14988073158111 0.14988073159234 0.15403226900995 0.15403226901180 0.15745313708439 0.15745313708491

47142.62 0 47142.62 0 47142.62 0 47142.62 0

0.10824687500033 0.10824687504457 0.11503988366702 0.11503988367752 0.12056851267216 0.12056851267466 0.12427185339487 0.12427185339550

0.12938468701688 0.12938468707287 0.13545482930581 0.13545482931860 0.14063709970603 0.14063709970900 0.14411524919932 0.14411524920007

0.14582616279746 0.14582616286247 0.15142830726817 0.15142830728273 0.15641598712935 0.15641598713270 0.15977116974643 0.15977116974727

47142.62 0 47142.62 0 47142.62 0 47142.62 0

0.10703307132398 0.10703306748271 0.11389660706352 0.11389660631228 0.11752312023844 0.11752311965567 0.12087125414680 0.12087125394444

0.12731472846772 0.12731472366747 0.13390913219055 0.13390913127857 0.13680264679656 0.13680264610529 0.14001621786293 0.14001621762547

0.14311657914158 0.14311657360471 0.14957885929202 0.14957885825581 0.15197910697568 0.15197910619983 0.15513165709199 0.15513165682702

47142.62 0 47142.62 0 47142.62 0 47142.62 0

0.11254950983790 0.11254951259647 0.11727474547038 0.11727474613772 0.12442410351849 0.12442410393211 0.12730098885619 0.12730098901218

0.13407094890075 0.13407095236095 0.13793975544464 0.13793975625561 0.14495927949797 0.14495927998955 0.14754557354429 0.14754557372754

0.15082996914830 0.15082997314719 0.15411704331367 0.15411704423561 0.16111615280040 0.16111615335269 0.16352362264570 0.16352362285031

300

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

0.20

0.20

ZA BLGSs Lx = 21.9777 nm, Ly = 4.9126 nm β= 4.47

BLGSs T = 300 K 0.15

0.15

W (nm)

W (nm)

3 0.10

AA (β= 4.48) ZZ (β= 4.47) AZ (β= 4.47) ZA (β= 4.46)

0.05

2 1

0.10

1:T=300 K, e0 a= 0.0 nm 2:T=500 K, e0 a= 0.2 nm 3:T=700 K, e0 a= 0.8 nm

0.05 Nonlocal plate model MD simulations

0.00

0

20

40

60

0.00

80

0

20

q 0 (MPa)

40

60

80

q0 (MPa)

(a)

Fig. 2. Load–deflection curves of four types of BLGSs subjected to uniform pressure at 300 K.

0.20

AZ BLGSs L x = 21.9777 nm, L y= 4.9126 nm β= 4.47

0.20

ZA BLGSs L x = 21.9777 nm, L y = 4.9126 nm T= 300 K

3

W (nm)

W (nm)

0.15

0.15

2 1

0.10

1:T=300 K, e0 a= 0.0 nm 2:T=500 K, e0 a= 0.1 nm 3:T=700 K, e0 a= 0.4 nm

0.10 0.05

1: e0 a = 0.0 nm 2: e0 a = 0.4 nm 3: e0 a = 0.6 nm 4: e0 a = 0.8 nm

0.05

Nonlocal plate model MD simulations 0.00 0

0.00

0

20

40

60

80

20

40

60

80

q0 (MPa)

(b)

q0 (MPa)

(a)

Fig. 4. Load–deflection curves of BLGSs under different thermal environmental conditions: (a) ZA type; (b) AZ type.

0.20

AZ BLGSs L x = 21.9777 nm, Ly = 4.9126 nm T= 300 K

W (nm)

0.15

0.10

1: e0a = 0.0 nm 2: e0a = 0.4 nm 3: e0a = 0.6 nm 4: e0a = 0.8 nm

0.05

0.00

0

20

40

60

80

q0 (MPa)

(b) Fig. 3. Load–deflection curves of BLGSs under different values of small scale parameter e0a subjected to uniform pressure at 300 K: (a) ZA type; (b) AZ type.

environmental conditions T = 300, 500 and 700 K. The temperature-dependent material properties are adopted, as given in Table 1. It can be seen that the deflections are increased with increase in temperature. This is due to the fact that the variation of temperature reduces the elastic moduli of the SLGS and degrades the

strength of the BLGS, and causes the larger deflections. The MD simulation results are also included for direct comparison. Through comparison, we find that the load–deflection curves obtained from the nonlocal plate model and MD simulations can match very well if the small scale parameters are properly chosen. For example, e0a = 0.0 nm at T = 300 K for both ZA and AZ BLGSs. This means the local thin plate model can provide an accurate prediction of the nonlinear bending behavior of BLGSs at a reference temperature T = 300 K. By nonlocal plate model, for AZ BLGSs, the small scale parameters need to be taken about e0a = 0.1 nm at T = 500 K and about e0a = 0.4 nm at T = 700 K. This means at higher temperature environmental conditions, i.e. beyond the reference temperature, only nonlocal plate model can provide an accurate prediction of the nonlinear bending behavior of BLGSs. Fig. 5 shows the effect of aspect ratio on the load–deflection curves of ZA and AZ BLGSs subjected to a transverse uniform load at T = 700 K. Like in the case of SLGSs [22], the load–deflection curves are increased with increase in the sheet aspect ratio. It is found that the difference between MD simulation results and the predictions from local plate model is decreased with increase in the sheet aspect ratio. By nonlocal plate model, the small scale parameters should be chosen about e0a = 1.4, 1.3 and 0.8 nm for ZA BLGSs and about e0a = 1.1, 0.9 and 0.4 nm for AZ BLGSs with b = 3.08, 3.69 and 4.77, respectively. This is due to the fact that

301

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302

Acknowledgments

0.20

3

ZA BLGSs T = 700 K

2

This work is supported in part by the National Natural Science Foundation of China under Grant 51279103. The authors are grateful for this financial support.

1

W (nm)

0.15

0.10

Appendix A

1:β= 3.08, e0 a= 1.4 nm 2:β= 3.69, e0 a= 1.3 nm 3:β= 4.47, e0 a= 0.8 nm

In Eqs. (32) and (33)

0.05 Nonlocal plate model

ð1Þ

AW ¼

MD simulations 0.00

p6 16

4096

ð3Þ

0

20

40

60

80

q0 (MPa)

BW ¼

ð3Þ Aw ¼

mnQ 11 ;

p

18 m3 n3 Q 4 11

p6 16

ð1Þ

ð1Þ

BW ¼ 1=AW ;

mnQ 13 ;

Q 13

ðA:1Þ

and (with J = I, II)

(a)

   2 ðIÞ ðIÞ ðIÞ Q 11 ¼ m4 þ 2c12 m2 n2 b2 þ c14 n4 b4 þ R11 C 1 ð1  v1 Þ   ðIÞ ðIÞ ðIÞ  c14 R11 cT1 m2 þ cT2 n2 b2 DT

0.20

AZ BLGSs T = 700 K

3 1

2

W (nm)

0.15

Q 13 ¼ ½ða313 þ a331  a333 ÞQ 11 þ H2    0.10

1:β= 3.08, e0 a= 1.1 nm 2:β= 3.69, e0 a= 0.9 nm 3:β= 4.47, e0 a= 0.4 nm Nonlocal plate model MD simulations 0

20

40

60

I

ðhI Þ2      1=2 D22 I A11 I A22 I

      2   2   2 ðIÞ ðIÞ ðIÞ ðIÞ ðIÞ 3 c24  c5 m4 þ c24 n4 b4 þ 4c5 c24 m2 n2 b2 c H2 ¼R11  2  2 > 16 ðIÞ ðIÞ > : c c24  c5 9 ðIÞ ðIÞ C 1 a311 þ C 3 ð1  v1 Þ3 16 8 > > <1

0.05

0.00

D11



80

q0 (MPa)

(b) Fig. 5. Effect of aspect ratio on the load–deflection curves of BLGSs at 700 K: (a) ZA type; (b) AZ type.

the material properties of SLGSs are size-dependent and temperature-dependent. It is clear that a large range of values for the small scale parameter e0a is possible due to different BLGSs under different thermal environmental conditions, and these values will be successfully used in solving nonlinear bending problems of BLGSs. 5. Concluding remarks Nonlinear bending response of BLGSs subjected to a transverse uniform load in thermal environments has been presented on the basis of a nonlocal plate model. The nonlocal stress incorporating the small scale parameter e0a is introduced into the governing equations. An approximate solution is obtained in an explicit form that can be served as benchmark for verifying numerically obtained solutions based on other continuum mechanics methods. The results confirm that the nonlinear van der Waals interaction coefficient has almost no effect on the nonlinear bending behavior of BLGSs, and may be neglected. The results show that the nonlinear bending behavior of BLGSs is sensitive to the small scale parameter e0a. The results reveal that the stacking sequence has a moderate effect, whereas the temperature change as well as the aspect ratio has a significant effect on the nonlinear bending response of BLGSs.

ðIÞ 14 ðIÞ 24

   2 ðJÞ ðJÞ ðJÞ g 11 ¼ m4 þ 2c12 m2 n2 b2 þ c14 n4 b4   ðJÞ ðJÞ ðJÞ ðJÞ  c14 R11 cT1 m2 þ cT2 n2 b2 DT þ R11 C 1    2 ðJÞ ðJÞ ðJÞ g 13 ¼ m4 þ 18c12 m2 n2 b2 þ 81 c14 n4 b4   ðJÞ ðJÞ ðJÞ ðJÞ  c14 R13 cT1 m2 þ 9cT2 n2 b2 DT þ R13 C 1    2 ðJÞ ðJÞ ðJÞ g 31 ¼ 81m4 þ 18c12 m2 n2 b2 þ c14 n4 b4   ðJÞ ðJÞ ðJÞ ðJÞ  c14 R31 9cT1 m2 þ cT2 n2 b2 DT þ R31 C 1 h   i ðJÞ ðJÞ ðJÞ g 33 ¼ 81 m4 þ 2c12 m2 n2 b2 þ c14 n4 b4   ðJÞ ðJÞ ðJÞ ðJÞ  9c14 R33 cT1 m2 þ cT2 n2 b2 DT þ R33 C 1 R11 ¼ 1 þ s2 p2 ðm2 þ n2 b2 Þ; R31 ¼ 1 þ s

2

2

2

R13 ¼ 1 þ s2 p2 ðm2 þ 9n2 b2 Þ;

2 2

p ð9m þ n b Þ

R33 ¼ 1 þ 9s2 p2 ðm2 þ n2 b2 Þ

1 a313 ¼ 16

ðIÞ

  ðIIÞ c c14 ðIIÞ ðIÞ ðIÞ ðIIÞ ðIIÞ ðIÞ 3 3ð1  v1 Þ3 g 13 C 3  C 1 C 3 R13 þ m4 14 ðIÞ g 13 þ ðIIÞ C 1 R13 v1

1 3ð1  v1 Þ a331 ¼ 16

c24

ðIÞ

ðIIÞ

ðIÞ

ðIIÞ

g 13 g 13  C 1 C 1 R213 3



c24

R13

   ðIIÞ ðIÞ ðIÞ ðIIÞ ðIÞ ðIÞ ðIIÞ ðIIÞ ðIIÞ ðIÞ g 31 C 3  C 1 C 3 R31 þ n4 b4 c14 c24 g 31 þ c14 c24 C 1 R31 v31 g 31 g 31  C 1 C 1 R231 ðIÞ

ðIIÞ

ðIÞ

ðIIÞ

R31

302

Y.-M. Xu et al. / Composite Structures 98 (2013) 294–302 ðIÞ ðIIÞ

a333 ¼ 

a311 ¼

ðIÞ

ðIIÞ

1 C 3 g 33  C 1 C 3 R33 ð1  v1 Þ3 ðIÞ R33 ; ðIIÞ ðIÞ ðIIÞ 16 g 33 g 33  C 1 C 1 R233

ðIIÞ

v1 ¼

C 1 R11 ðIIÞ g 11

R11 ðIIÞ

16g 11    8 9 2  2   2  2 ðIIÞ ðIIÞ ðIIÞ ðIIÞ ðIIÞ > > > 3 c24  c5 m4 þ c24 n4 b4 þ 4c5 c24 m2 n2 b2 > < = cðIIÞ 3 3 14  9C ðIIÞ ð1  v Þ  v ðA:2Þ     1 1 3 2 2 ðIIÞ > > ðIIÞ ðIIÞ c > > 24 : ; c24  c5

References [1] Novoselov KS. Nobel lecture: graphene: materials in the flatland. Rev Modern Phys 2011;83:837–49. [2] Reddy CD, Rajendran S, Liew KM. Equilibrium configuration and continuum elastic properties of finite sized graphene. Nanotechnology 2006;17: 864–70. [3] Scarpa F, Adhikari S, Phani AS. Effective elastic mechanical properties of single layer graphene sheets. Nanotechnology 2009;20:065709. [4] Ni Z, Bu H, Zou M, Yi H, Bi K, Chen Y. Anisotropic mechanical properties of graphene sheets from molecular dynamics. Physica B 2010;405:1301–6. [5] Zhang YY, Wang CM, Cheng Y, Xiang Y. Mechanical properties of bilayer graphene sheets coupled by sp3 bonding. Carbon 2011;49:4511–7. [6] Zhang B, Cui T. An ultrasensitive and low-cost graphene sensor based on layerby-layer nano self-assembly. Appl Phys Lett 2011;98:073116. [7] Bunch JS, van der Zande AM, Verbridge SS, Frank IW, Tanenbaum DM, Parpia JM, et al. Electromechanical resonators from graphene sheets. Science 2007;315:490–3. [8] Wang X, Zhi L, Müllen K. Transparent, conductive graphene electrodes for dyesensitized solar cells. Nano Lett 2008;8:323–7. [9] Stankovich S, Dikin DA, Dommett GHB, Kohlhaas KM, Zimney EJ, Stach EA, et al. Graphene-based composite materials. Nature 2006;442:282–6. [10] Hsiao MC, Liao SH, Yen MY, Teng CC, Lee SH, Pu NW, et al. Preparation and properties of a graphene reinforced nanocomposite conducting plate. J Mater Chem 2010;20:8496–505. [11] Poot M, van der Zant HSJ. Nanomechanical properties of few-layer graphene membranes. Appl Phys Lett 2008;92:063111. [12] Duan WH, Wang CM. Nonlinear bending and stretching of a circular graphene sheet under a central point load. Nanotechnology 2009;20:075702.

[13] Hemmasizadeh A, Mahzoon M, Hadi E, Khandan R. A method for developing the equivalent continuum model of a single layer graphene sheet. Thin Solid Films 2008;516:7636–40. [14] Scarpa F, Adhikari S, Gil AJ, Remillat C. The bending of single layer graphene sheets: the lattice versus continuum approach. Nanotechnology 2010;21:125702. [15] Wang Q. Simulations of the bending rigidity of graphene. Phys Lett A 2010;374:1180–3. [16] Shen YK, Wu HA. Interlayer shear effect on multilayer graphene subjected to bending. Appl Phys Lett 2012;100:101909. [17] Xu X, Liao K. Molecular and continuum mechanics modeling of graphene deformation. Mater Phys Mech 2001;4:148–51. [18] Eringen AC. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. J Appl Phys 1983;54:4703–10. [19] Lam DCC, Yang F, Chong ACM, Wang J, Tong P. Experiments and theory in strain gradient elasticity. J Mech Phys Solids 2003;51:1477–508. [20] Yang F, Chong ACM, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39:2731–43. [21] Sun ZH, Wang XX, Soh AK, Wu HA, Wang Y. Bending of nanoscale structures: inconsistency between atomistic simulation and strain gradient elasticity solution. Comput Mater Sci 2007;40:108–13. [22] Shen H-S, Shen L, Zhang C-L. Nonlocal plate model for nonlinear bending of single-layer graphene sheets subjected to transverse loads in thermal environments. Appl Phys A 2011;103:103–12. [23] Huang Y, Wu J, Hwang KC. Thickness of graphene and single-wall carbon nanotubes. Phys Rev B 2006;74:245413. [24] He XQ, Wang JB, Liu B, Liew KM. Analysis of nonlinear forced vibration of multi-layered graphene sheets. Comput Mater Sci 2012;61:194–9. [25] Shen L, Shen H-S, Zhang C-L. Temperature-dependent elastic properties of single layer graphene sheets. Mater Des 2010;31:4445–9. [26] Stuart SJ, Tutein AB, Harrison JA. A reactive potential for hydrocarbons with intermolecular interactions. J Chem Phys 2000;112:6472–86. [27] Hoover WG. Canonical dynamics: equilibrium phase-space distributions. Phys Rev A 1985;31:1695–7. [28] Shen H-S, Zhang C-L. Torsional buckling and postbuckling of double-walled carbon nanotubes by nonlocal shear deformable shell model. Compos Struct 2010;92:1073–84. [29] Shen H-S, Zhang C-L. Nonlocal shear deformable shell model for post-buckling of axially compressed double-walled carbon nanotubes embedded in an elastic matrix. J Appl Mech ASME 2010;77:041006. [30] Shen H-S, Zhang C-L, Xiang Y. Nonlocal shear deformable shell model for thermal postbuckling of axially compressed double-walled carbon nanotubes. Philos Mag 2010;90:3189–214.