Buckling analysis of non-uniform thickness nanoplates in an elastic medium using the isogeometric analysis

Accepted Manuscript Buckling analysis of non-uniform thickness nanoplates in an elastic medium using the Isogeometric Analysis T. Banh-Thien, H. Dang-Trung, L. Le-Anh, V. Ho-Huu, T. Nguyen-Thoi PII: DOI: Reference:

S0263-8223(16)31460-X http://dx.doi.org/10.1016/j.compstruct.2016.11.092 COST 8053

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

6 August 2016 25 November 2016 30 November 2016

Please cite this article as: Banh-Thien, T., Dang-Trung, H., Le-Anh, L., Ho-Huu, V., Nguyen-Thoi, T., Buckling analysis of non-uniform thickness nanoplates in an elastic medium using the Isogeometric Analysis, Composite Structures (2016), doi: http://dx.doi.org/10.1016/j.compstruct.2016.11.092

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Buckling analysis of non-uniform thickness nanoplates in an elastic medium using the Isogeometric Analysis T. Banh-Thien, H. Dang-Trung, L. Le-Anh, V. Ho-Huu, T. Nguyen-Thoi* Division of Computational Mathematics and Engineering, Institute for Computational Science, Ton Duc Thang University, Vietnam Faculty of Civil Engineering, Ton Duc Thang University, Vietnam

Abstract The paper presents a new numerical approach for buckling analysis of non-uniform thickness nanoplates in an elastic medium using the isogeometric analysis (IGA). By ignoring the van der Waals interaction between two adjacent plates, non-uniform thickness nanoplates are described as a single-layered graphene sheet. The governing differential equation of the nanoplates is derived by the nonlocal theory in which the nonlocal stress-strain relation is used to capture the nonlocal mechanics caused by small size effects. The governing equation is then discretized into algebraic equations and solved by using IGA procedure to determine the critical buckling load. By using the non-uniform rational B-splines, IGA easily satisfies the required continuity of the partial differential equations in buckling analysis. Several numerical examples are solved and compared with those of previous publications to illustrate the performance of IGA for buckling analysis of nanoplates. Keywords: nanoplate; non-local theory; buckling load; isogeometric analysis.

1. Introduction Since carbon nanotubes (CNTs) were discovered by Iijima [1] in 1991, nanomaterials have attracted a huge amount of attention from the scientific community. Due to their extraordinary mechanical, chemical, thermal, electrical and electronic properties, structures made from nanomaterials such as nanobeams [2], nanorods [3], nanoribbons [4], nanoplates [5], nanoshells [6], etc. have been studied and applied to many engineering disciplines. Their applications as sensors [7], [8], atomic dust detectors [9], enhancer of surface image resolution [10], and so on, have been widely developed in the past few years. Among *

Corresponding author. Tel.: +84 933 666 226. E-mail address: [email protected] (T. Nguyen-Thoi)

nanostructures, nanoplate is considered as the fundamental structure. Graphene sheets consist of carbon atoms connected to each other via covalent bonds in the honeycomb lattice. Because many of the carbon-based nanostructures are considered as deformed graphene sheets, the mechanical behavior of nanoplates is very important in the field of nanomaterials. Owing to the fact that controlled experiments at nanoscale are difficult to execute, theoretical analysis and numerical simulations become favorable among the scientific community to simulate the behaviors of nanostructures. In terms of numerical simulation, two approaches are commonly used for mechanical analysis of nanostructures namely molecular dynamic modelling and continuum modelling. As molecular dynamic simulations are highly computationally expensive, continuum model is rapidly developed in the past few years. Various size-dependent continuum models capturing small-scale effects were proposed for nanostructures analysis such as couple stress theory [11], strain gradient theory [12], modified couple stress theory [13], etc. Among them, nonlocal elasticity theory proposed by Eringen [14], [15] has emerged as one of the most promising theories for analyzing smallscale structures. In this theory, the stress at a point is assumed as a function of the strain at all points in the problem domain. This helps capture the small scale effects when analyzing structures. On the other front of numerical simulations, Hughes et al. [16] recently proposed a novel numerical method based on weak-form Galerkin formulation named isogeometric analysis (IGA) for numerical simulation. This method integrates both of the geometric description and finite element approximation by the same basic function space of B-Spline or Non-Uniform Rational B-Splines (NURBS). One of advantages of IGA compared to other existing numerical methods includes flexible control of basic functions continuity. This helps fulfill the high-order continuity requirement of partial differential equations (PDEs) naturally. In addition, the IGA produces the seamless integration of computer aid design and finite element analysis tools. This simplifies the costly model generation procedure which is one of the major bottlenecks in engineering analysis-design. The IGA then has been widely developed and applied to many engineering disciplines such as fluid-structure interaction [17], [18], [19], fluids and turbulence [20], [21], shape optimization [22], [23], electromagnetics [24], composite material [25] and fracture mechanics [26]. IGA has also been widely applied to the plate structures analysis in association with various plate models such as classical plate theory [27], [28], first-order shear deformation theory [29], [30] highorder shear deformation theory [31], [32], [33], layer-wise [34], [35], etc.

Regarding to the numerical simulation of nanopates, by using the non-local continuum plate model, Pradhan and Murmu et al. [36], [37] studied the buckling of SLGS using differential quadrature method (DQM). Ravari et al. [38] investigated buckling of nanoplates by using the finite difference method (FDM). Ansari et al. [39] studied the vibration of SLGSs by using the generalized differential quadrature (GDQ) method. In the abovementioned studies, the solutions are obtained from strong-form differential equation. This is intuitive and easy to be implemented for simple problems. However, mathematical formulas might become unwieldy when it comes to solving more complicated problem such as ones having moving boundaries or an unstructured grid. The paper thus proposes a new simulation of buckling for nanoplates by using the IGA. With the flexiblity in controling of basic functions continuity and the seamless integration of computer aid design, the IGA is expected to overcome the above-mentioned drawbacks of existing methods. By ignoring the van der Waals interactions between any two adjacent layers, a non-uniform thickness nanoplate is described as a SLGS [40]. Both Winkler-type [41] and Pasternak-type foundations [42] are used to simulate the interaction between SLGS and elastic medium. In this work, the governing equations of nanoplates subjected to uniaxial and biaxial compressions are derived by combining the classical plate and nonlocal theory. IGA is then employed to determine the critical buckling load for the plates. The influences of the scale coefficient, plates’ length, the aspect ratio on buckling behavior of nanoplates are also investigated. The obtained results are compared with those by previous publications to illustrate the accuracy and reliability of the proposed method.

2. Governing equation for nanoplates using the nonlocal theory To formulate the governing equations, a graphene sheet resting on the elastic medium as shown in Figure 1 is examined. The plate’s dimensions in Cartesian coordinates with three coordinates (x,y,z) are taken along the length, width, thickness of the plate, respectively. According to the nonlocal theory [14], the stress field σ ij at a reference point x is produced by local strain εij at the source point x’ by: σ ijn l ( x ) =

∫

ψ

V (x ′ )

( x − x ′ ,α )C

ijk l

ε kl ( x ′ ) d V ( x ′ )

(1)

where

σijnl , εkl and

C ijkl

are, respectively, the stress, strain and fourth-order elasticity tensors;

ψ ( x − x′ ,α ) is the nonlocal modulus or attenuation function which incorporates into constitutive equations; x − x′ represents the Euclidean distance between the reference point and the source point; α = e 0 a l x is the material constant which depends on the internal length a (lattice parameter, granular size, distance between C-C bonds), external characteristics lengths lx and e0 is a constant which could be calibrated from experimental results. When the internal characteristic length (C–C bonds length) is negligible compared to the external characteristic length, α approaches zero, then the nonlocal elasticity theory reduces to the local classical elasticity theory. Because the integral in Eq. (1) is difficult to implement into the calculation, a simplified equation is used to replace Eq. (1) as a basic of all nonlocal constitutive formulation. This simplified equation is given by [43]: 1 − ( e 0 a ) 2 ∇ 2  σ ij = C ijkl ε kl  

(2)

where e0a is the nonlocal scale coefficient denoting the small scale effect on the response of

(

)

2 2 2 2 2 structures in nanosize; ∇ 2 is the Laplacian operator given by ∇ = ∂ / ∂x + ∂ / ∂y for

two-dimensional space. The basic stress-strain relationship of an orthotropic graphene sheet in Cartesian coordinates can be expressed in the following form:

σ xxnl  σ xxnl, xx + σ xxnl, yy   E11 E12 0  ε xx  2  nl  nl      nl σ yy  − ( e0 a ) σ yy, xx + σ yy, yy  =  E21 E22 0  ε yy  nl  nl   nl  0 E33  γ xy  σ xy  σ xy, xx + σ xy, yy   0

(3)

where

E11 =

E1 1 −ν12ν 21

,

E22 =

E2 1 −ν12ν 21

,

E12 = E21 =

ν12 E2 , 1 −ν12ν 21

E33 = G12

(4)

in which E1 and E2 are Young’s modulus; G12 = E1 / 2(1+ν12 ) is the shear modulus and ν 1 2 , ν 2 1 are Poisson’s ratio.

The displacement field at any point of the plate is expressed as:

u ( x, y, z) = u0 ( x, y) − z w0,x ( x, y)

(5)

v ( x, y, z ) = v0 ( x, y) − zw0, y ( x, y)

(6)

w( x, y, z ) = w0 ( x, y)

(7)

where ( u0 , v0 , w0 ) are the displacement of an arbitrary point on the plate’s middle surface. The small strains and moderate rotations in the strain-displacement relations take the form:

εxx = u, x + w,2x / 2,

ε yy = v, y + w,2y / 2,

γ xy = u, y + v,x + w,xw, y

(8)

Substituting Eqs. (5)- (7) into Eq. (8), one obtains:

(

)

1 ε xx = ε xx0 + z ε xx

(9)

ε yy = ε 0yy + zε1yy

(10)

γ xy = γ xy0 + zγ xy1

(11)

(

)

0 0 0 1 1 1 where ε xx , ε yy , γ xy and ε xx , ε yy , γ xy are the membrane and bending strains defined as

ε xx0 = u0,x + w0,2 x / 2,

ε 0yy = u0, y + w0,2 y / 2,

ε 1xx = − w0, xx ,

γ xx0 = u0, y + v0, x + w0, x w0, y

ε 1yy = − w0, yy ,

γ 1xy = −2 w0, xy

(

(12)

)

The non-local resultant forces N xx , N yy , N xy and the non-local resultant moments

(M

xx

, M yy , M xy ) are defined by  ( N xx , N yy , N xy ) , ( M 

xx

,M

yy

,M

xy

) = ∫

h/2 − h /2

(σ

nl xx

,σ

nl yy

nl , σ xy ) (1, z ) d z

(13)

Substituting Eqs. (3) and (9)-(11) into Eq. (13), one obtains:

 N xx   N xx, xx + N xx, yy   B11 B12     2  N yy  − (e0 a)  N yy , xx + N yy, yy  =  B21 B22     0  N xy   N xy , xx + N xy, yy   0

0  ε xx0    0  ε yy0   B33  γ xy0 

(14)

 M xx   M xx, xx + M xx, yy   D11 D12     2  M yy  − (e0 a)  M yy, xx + M yy, yy  =  D21 D22     0  M xy   M xy, xx + M xy , yy   0

(15)

0  ε 1xx    0  ε 1yy  D33  γ 1xy 

where B nm = hE nm , D n m = h 3 E n m / 1 2 are the extensional and flexural rigidities, respectively. The virtual strain energy of the plate is given by: h /2

δU = ∫

∫ (σ

δε xx + σ yynl δε yy + σ yynl δγ xy ) dxdydz

nl xx

Ω − h /2

(16)

h /2

=

∫ (N

0 xx

0 yy

0 xy

1 xx

1 yy

δε + N yyδε + N xyδγ + M xxδε + M yyδε + M xyδγ

xx

1 xy

) dxdy

− h /2

The virtual word done by external forces is written as

δV = −∫ qδ w0 dxdy − ∫ pδ w0 dxdy Ω

(17)

Ω

where q is the distributed transverse load at the surface; p = k w w 0 − k s ∇ 2 w 0 represents the pressure–deflection relationship created by plate-foundation interaction by Pasternak model; kw is Winkler foundation modulus and ks is the foundation stiffness of Pasternak model.

It’s noted that the principle of virtual work is independent of the constitutive relations [37]. As a result, this principle can be used to derive the equilibrium equations of the nonlocal plates. By substituting Eqs. (16) and (17) into the principle of virtual displacements for static, the equilibrium equations can be obtained as: δ u0 :

Nxx,x + Nxy, y = 0

(18)

δ v0 :

Nxy,x + Nyy, y = 0

(19)

δ w0 :

M xx , xx + 2 M xy , xy + M yy , yy + ( N xx w0, x + N xy w0, y )

,x

+ ( N xy w0, x + N yy w0, y ) + q − k w w0 + k s ( w0, xx + w0, yy ) = 0

(20)

,y

Substituting Eqs. (15) into Eq. (20), the governing differential equation for the stability of nanoplate is given as

(D

11

w0, xx + D12 w0, yy )

, xx

+ ( D12 w0, xx + D22 w0, yy )

, yy

+ ( 4 D33 w0, xy )

, xy

− ( N xx w0, x + N xy w0, y ) + ( N xy w0, x + N yy w0, y ) + q − k w w0 + k s ( w0, xx + w0, yy )  ,x ,y  

(21)

2 + ( e0 a ) ∇ 2  ( N xx w0, x + N xy w0, y ) + ( N xy w0, x + N yy w0, y ) + q − k w w0 + k s ( w0, xx + w0, yy )  = 0 ,x ,y  

Note that the governing differential equation returns to classical plate equation when the scale coefficient e0a is set to be zero. In this study, the load N

xx

and N y y are considered to

have the relationships as: N xx = N 0 ,

N yy = µ N 0 ,

N xy = q = 0

(22)

It’s assumed that the thickness of SLGS varies in x-axis and y-axis as shown in Figure 2 following:

h = h0λ ( x, y) where h0 is constant and

(23)

λ = λ ( x, y ) ≥ 1 is a function in the x,y-axis. Due to the fact that the

rigidities of the nanoplate depend on the thickness, from Eq. (23) we have: D11 = D110 λ 3 ( x , y ) ,

0 D 22 = D 22 λ 3 ( x, y )

D12 = D120 λ 3 ( x , y ) ,

0 D 33 = D33 λ 3 ( x, y )

(24)

with D n0m = h03 E n m / 1 2 . By substituting Eqs. (22) and (24) into Eq. (21), one obtains A1w0,xxxx + A2 w0, yyyy − A3 w0, xxx − A4 w0, yyy − A5 w0, xx − A6 w0, yy + A7 ( w0, xxyy + w0, yyxx ) − A8 w0, xxy − A9 w0, xyy − A10 ( w0, xy + w0, yx ) + A11 w0

(

(25)

)

+ N 0 ( e0 a ) w0,xxxx + µ w0, yyyy + (1 + µ ) ( w0, xxyy + w0, yyxx ) / 2 − ( w0, xx + µ w0, yy ) = 0   2

where A1 , A2 , A3 , A4 , A5 , A6 , A7 , A8 , A9 , A10 and A11 are defined in Appendix 1 . By using Galerkin method, multiplying with the test function integrating over the domain, the Eq. (25) becomes

φ ∈Η20 ( Ω) and

∫

Ω

{A w 1

0, xxxx

+ A2 w0, yyyy − A3 w0, xxx − A4 w0, yyy − A5 w0, xx − A6 w0, yy

+ A7 ( w0, xxyy + w0, yyxx ) − A8 w0, xxy − A9 w0, xyy − A10 ( w0, xy + w0, yx ) + A11 w0

(

(26)

}

)

2 + N 0 ( e0 a ) w0, xxxx + µ w0, yyyy + (1 + µ ) ( w0, xxyy + w0, yyxx ) / 2 − ( w0, xx + µ w0, yy )  φ d Ω = 0  

Applying integration by parts theorem, the weak form of Eq. (25) can be obtained as:

∫ A dΩ+ ∫ B d∂Ω = 0 Ω

(27)

∂Ω

where A = A1w0, xxφ, xx + A2 w0, yy φ, yy + A3 w0, xxφ, x + A4 w0, yyφ, y + A5 w0, xφ x + A6 w0, yφ, y + A7 ( w0, yyφ, xx + w0, xxφ, yy ) + A8 w0, xxφ, y + A9 w0, yyφ, x + A10 ( w0, yφ, x + w0, xφ y )

(

+ A11 w0φ + N 0 ( e0 a ) w0, xxφ, xx + µ w0, yyφ, yy + (1 + µ ) ( w0, yyφ, xx + w0, xxφ, yy ) / 2 2

(28)

)

+ N 0 ( w0, xφ, x + µ w0, yφ, y )

and

( + ( A (w φ − w φ ) + A (w φ − w φ ) − ( A w − N (( w φ ) n + µ ( w φ ) n ) + N (e a ) (( w φ − w

) )φ ) n

B = A1 ( w0, xxxφ − w0, xxφ, x ) + A7 ( w0, xyyφ − w0, yyφ, x ) − ( A3 w0, xx + A5 w0, x + A9 w0, yy + A10 w0, y ) φ nx 2

0, yyy

0, yy , y

7

0, xxy

0, xx , y

4

2

0

0, x

x

0, y

y

0

0

0, xxx

0, yy

+ A6 w0, y + A8 w0, xx + A10 w0, x

0, xxφ, x ) nx + µ ( w0, yyyφ − w0, yyφ, y ) n y

+ (1 + µ ) ( w0, xyyφ − w0, yyφ, x ) nx / 2 + (1 + µ ) ( w0, yxxφ − w0, xxφ, y ) n y / 2

y

(29)

)

in which nx , ny are direction cosines of the unit normal on the boundary.

3. Isogeometric nanoplate formulation for buckling analysis

In this section, fundamental components of IGA (i.e. knot vector and basis functions) and its implementation for nanoplates formulations are presented. More details could be found in [14] for interested readers. 3.1. A brief on knot vector and basis functions

{

}

In IGA, a knot vector [Ξ] = ξ1 , ξ2 ,..., ξn+ p +1 is defined as a non-decreasing sequence of parameter values ξi , ξ i ∈ R where ξi is called ith knot in the parametric space; p is the order

of the B-spline and n is the number of basis functions. The B-spline basis function is C∞ continuous inside a knot span and Cp-1 continuous at each unique knot. Each repetition of a knot decreases one degree of continuity at the repeated knot. Degree p of basis functions have up to p − 1 continuous derivatives.

[Ξ] = {0,0,0,0.5,1,1,1} . The univariate B-spline basis functions Ni, p ( ξ ) are defined by the Cox-de Boor recursive formula and is given by: For p = 0

1 Ni,0 (ξ ) =  0

if ξi ≤ ξ ≤ ξi +1

otherwise.

(30)

For p ≥ 1 the definition is

N i, p (ξ ) =

ξ −ξ ξ − ξi Ni , p−1 (ξ ) + i + p+1 N (ξ ) ξi + p − ξ i ξi+ p +1 − ξi+1 i+1, p −1

(31)

The B-spline surfaces are defined by the tensor product of univariate B-spline basis functions in two parametric dimensions ξ and η associated with two-knot vectors

[Ξ] = {ξ1 , ξ2 ,..., ξn+ p+1} and [Θ] = {η1,η2 ,...,ηm+q+1} , expressed by n×m

n ×m

I =1

I =1

S (ξ ,η ) = ∑ N i , p (ξ ) M j ,q (η ) Pi , j = ∑ N IB (ξ ,η ) PI

(32)

where Ni, p ( ξ ) and M j ,q (η ) are univariate B-spline basis function defined as Eqs. (30)-(31);

NIB (ξ ,η) is the B-spline shape function associated with the node I = i ( p +1) + j ; PI is the bidirectional control net. An example of 1D and 2D basis functions formed from knot vector

[Ξ] = {0,0,0,0.5,1,1,1} is illustrated in Figure 3. To present exactly the geometry of a curve, non-uniform rational B-splines (NURBS) is used. Different from B-spline, each control point of NURBS has an additional value called g

individual weight wi, j . Then NURBS surface can be expressed as

n×m

S (ξ ,η ) = ∑ N I ( ξ ,η ) PI

with N I ( ξ ,η ) =

I =1

N IB wIg n ×m

∑N

B J

(33)

wJg

J =1

3.2 NURBS formulation for buckling analysis of nanoplates The displacement w is approximated by NURBS basis functions as n× m

w = ∑ N I ( ξ ,η ) wI

(34)

I =1

where wI is the displacement at node I. Substituting Eq. (34) into Eq. (27), the discretized system of equations for buckling analysis of nanoplates can be expressed as

∑{ K

E IJ

+ N0 K IJG }wI = 0

(35)

I

where

K IJE

 A1 N I , xx N J , xx + A2 N I , yy N J , yy + A3 N I , xx N J , x + A4 N I , yy N J , y + A5 N I , x N J , x    = ∫  + A6 N I , y N J , y + A7 ( N I , xx N J , yy + N I , yy N J , xx ) + A8 N I , xx N J , y + A9 N I , yy N J , x  dΩ  Ω   + A N N + N N + A N N ( ) 10 I,y J ,x I ,x J,y 11 I J  

(

( e0 a )2 N I , xx N J , xx + µ N I , yy N J , yy + (1 + µ ) ( N I , yy N J , xx + N I , xx N J , yy ) / 2 K = ∫  + ( N I , x N J ,x + µ N I , y N J , y ) Ω  G IJ

) dΩ

(36)

(37)

 

Eq. (35) can be reduced to the following matrix form:

{K E

+ N 0K G }w = 0

(38)

nel

in which Kα = ∑Kαe (α = E, G ) is defined by e=1

K eE = ∫ e B T DB dΩ; K eG = ∫ e B TG D G B G dΩ Ω

Ω

(39)

where B , B G , D and DG are derivative matrices of basis function, defined in Appendix 2. The jacobian matrix expressing the relationship between physical coordinate system and NURBS parameter space is obtained by

 x ,ξ J=  y ,ξ

x ,η  y ,η 

(40)

where x,ξ = ∑ I RI ,ξ (ξ ,η )Px ;

x,η = ∑ I RI ,η ( ξ ,η )Px

y,ξ = ∑ I RI ,ξ (ξ ,η )Py ;

y,η = ∑ I RI ,η ( ξ ,η )Py

n×m

n×m

n ×m

n×m

(41)

in which Px , Py are the compositions of control net. The derivative components ∂ , x , ∂ , y , ∂ , x x , ∂ , y y , ∂ , xy of derivative matrices in physical space are calculated from derivative components in parametric space by using the inverse matrix T (Appendix 2) as

{∂

T

,x

T

, ∂ , y , ∂ , xx , ∂ , yy , ∂ , xy } = T −1 {∂ ,ξ , ∂ ,η , ∂ ,ξξ , ∂ ,ηη , ∂ ,ξη }

(42)

4. Numerical results

In this section, numerical results obtained by the present method are presented to illustrate the robustness and effectiveness of the IGA for buckling analysis of non-uniform thickness orthotropic graphene sheet embedded in the elastic medium. The properties of the orthotropic graphene sheet in this paper are similar to those in the reference [44]. The fourth-order basis functions (p = q = 4) along with 12 control points in each direction are used to approximate the displacement in the examples 4.1 and 4.2. In the examples 4.3, 4.4 and 4.5, the displacement is approximated by third-order basis functions and 11 control points in each direction. In all examples, the meshing 8×8 and the weight values { wi } =1 are used. To investigate the effect of nonlocal theory on the behavior of the nanoplates, buckling load ratio (or critical load ratio) is calculated as:

λ=

Buckling load from nonlocal theory Buckling load from local theory

(43)

In order to verify the accuracy of the proposed method, the numerical results obtained by the paper are compared with those by Pradhan et al. [45], [37], Ravari et al. [38], M. Eisenberger [46] and Farajpour et al. [40]. The effect of the length L, thickness t of the plate,

the scale coefficient e0a , the elastic medium parameters kw and ks on the buckling load of uniaxial and biaxial compressed plates are respectively investigated in the following sections to provide a more comprehensive view of nanoplates buckling behavior. 4.1 The effect of length, small scale and boundary conditions on the buckling load of SLSG First of all, the effect of the plates’ length and scale coefficient on buckling behavior is investigated. For this purpose, a simply supported square nanoplate subjected to an uniaxial load N 0 , as shown in Figure 4 (a) is considered. The material properties of the plate are given by: Young’s modulus E = 1.06 TPa and Poisson’s ratio v = 0.25. The plate has no defect and is considered to have a constant temperature during loading. The Winkler and shear factors are ignored. The thickness of the plate is h = 0.34 nm and the length of the plate is varied from

5.0

e 0 a = {0.0

to 0.5

1.0

45.2896 1.5

nm.

The

scaling

coefficient

is

taken

as

2.0} nm .

The variations of buckling load ratio corresponding to different values of plates’ length and scale coefficients are shown in Table 1. In addition, a comparison between the results by the present method and those by Pradhan et al. [37] is shown in Figure 5. It can be seen that the results by the IGA in association with the nonlocal theory agree very well with those by the DQM used in research of Pradhan et al. [37]. The scaling coefficient e0 a has a significant effect on the buckling behavior of nanoplates. In specific, the buckling load ratio in the case e 0 a = 2 and L = 5 n m is six times less than that in the case e 0 a = 0 . Hence, it could be

concluded that small-scale effects crucially influence the buckling characteristics of nanoplates. In addition, when the plate’s length increases, the buckling load by nonlocal theory converges with classical theory ( e 0 a = 0 ). This phenomenon is also obtained by molecular mechanics modelling. It’s also seen that the buckling load of the nanoplates based on nonlocal elasticity theory is always smaller than that of local theory. The effect of the scale coefficient on the buckling behavior of the plate with various aspect ratios l x l y is also considered. The buckling load and scale coefficient are nondimensionalized by

λ1 = N0lx2 D;

ψ = e0a ly

(44)

The results obtained by the present method for both cases of fully simply supported plate and fully clamped plate are compared with those by Ravari et al. [38] and shown in Figure 6. It can be seen that the results by IGA agree well with those by FDM from Ravari et al. [38] for both cases of boundary conditions with different values of aspect ratio. It’s also observed that as the non-dimensional scale coefficient

ψ

increases, the non-dimensional

buckling load λ1 decreases for both boundary conditions. In addition, when the nondimensional aspect ratio l x l y increases, the non-dimensional buckling load λ1 decreases in both cases. These trends coincide with conclusions in previous sections where buckling load ratio and scale coefficient are not non-dimensionalized. This again demonstrates the stability of IGA for buckling analysis of nanoplates. Figure 7 shows the non-dimensional buckling load values for different boundary conditions, in which the plate size is fixed with length l x = l y = 4 5.2 9 nm. It can be seen that the non-dimensional buckling load λ1 obtained the

highest value for the fully clamped nanoplate and the lowest one for the fully simply supported nanoplate as expected. 4.2 Effects of aspect ratio on buckling of uniaxial and biaxial compressed SLGS The nonlocal effects on the buckling behavior of nanoplates under one and two directional compressions are investigated in this section. For this purpose, a nanoplate having the same properties as in section 4.1 is considered. The uniaxial and biaxial load cases as shown in Figure 4 are solved with three different plate’s side lengths: L = 5, 10, 15 nm . The compression ratio µ is assumed to be unity for all cases. The results of buckling load ratio for both uniaxial and biaxial compressions with different side lengths are presented in Table 2 and Figure 8. The figure shows that the results by IGA agree well with those by DQM used in the study of Pradhan et al. [45]. This again illustrates that IGA is reliable and stable for both uniaxial and biaxial buckling analyses of nanoplates. It can be seen that for all cases, the nonlocal solutions are always smaller than that of the local counterparts. This coincides with the phenomenon observed in section 4.1. It’s also seen that the small-scale effects in the case of biaxial compression are relatively less than those in the case of uniaxial compression, especially for large scale coefficient ( e0 a = 2 ). However, this gap gradually diminishes when the plate’s length increases. At the plate’s length L = 25 nm , the buckling loads for both uniaxial and biaxial compressions have the

same values. This is understandable because the size-effect reduces with the increase of plate’s length. The effects of aspect ratio on buckling of uniaxial and biaxial compression graphene sheets are next conducted. The length of the plate is taken as 10 nm. Aspect ratio lx l y is varied from 0.2 to 2. In this example, the scale coefficient e0a is chosen to be 2, and the Winkler and shear factors are ignored. The results of load ratio with respect to the aspect ratio are presented in Figure 9 in comparison with those by Pradhan et al. [45]. It can be observed that the load ratio decreases when the aspect ratio increases from 0.2 to 2.0. In addition, the effects of small scale in the case of biaxial compression are less than those in the case of uniaxial one for all aspect ratios. The difference between the load ratio of uniaxial compression and biaxial compression becomes larger when the aspect ratio increases. It is also seen that the present results are closely equal to those by Pradhan et al. This again reinforces the accuracy of IGA for buckling analysis of nanoplates. Next, to further demonstrate the advantage of using the IGA, an investigation on approximation order and mesh size for the use of the IGA is carried out. The geometry of the considered plate are L = 45 nm, H = L / 3 and t = 0.34 . The scaling coefficient is e0 a = 0.5 while other material properties of the plate are similar to those given in the above examples. The bulking loads corresponding to different values of mesh size and order of displacement field approximation are plotted in Figure 10. It can be seen in Figure 10 that when the displacement field approximation order is larger than 5, the effect of mesh size to the results becomes less significant. This is one of the most useful aspects of the IGA. It does not only control the mesh-size but also control the order of the field approximation by using prefinement. And because this p-refinement is very easy to be carried out in the IGA, it helps the IGA achieved the high accuracy with coarse meshes. This would helps significantly reduce computational cost for analyzing structures. 4.3 The effect of elastic medium parameters on the buckling load of SLGS In this section, the effect of elastic medium on bucking load is investigated. For this purpose, the non-dimensional Winkler parameter KW and the non-dimensional shear parameter K S are defined as: KW = k wlx4 / D ,

K S = k s lx2 / D

(45)

The results presented in Figure 11 show the variation of critical load ratio with respect to different values of Winkler stiffness in comparison with those by Pradhan’s [37]. In this case, the non-dimensional shear parameter is kept as constant K S = 0 . It can be seen that for higher values of e0a , the critical load ratios are lower. It is also observed that the critical load ratio increases and decreases unequally. In addition, the effect of the approximation (polynomial order) and the mesh size on critical load ratio for the case e0a = 2 is shown in Figure 12. It is seen from this figure that for small values of Winkler modulus parameter, the mesh size and approximation order have slight effect on the bulking load. However, for higher values of modulus parameter, they have more significant effect on the critical bulking load. Figure 13 shows the effect of shear modulus parameter on buckling load ratio, corresponding to the case K W = 200 . Similar to the Figure 11, when the parameter e0a increases, the critical load ratio decreases. Figure 14 depicts the influence of aspect ratio on the critical buckling load for SLGS embedded in the elastic medium. The elastic medium parameter is fixed as K W = 10, K S = 2 while the aspect ratio lx l y increases from 0.2 to 1.0. From the figure, it is seen that for higher aspect ratio, the effect of e0a on the critical load ratio values increases. In other words, the effect of nonlocal parameters for the square SLGS is stronger than that for the rectangular one. 4.4 Nanoplates with linear thickness variation In this example, the effect of the plate thickness on buckling load is studied. A simply

(

supported isotropic classic plates, e 0 a = 0 , with linearly variable thickness h = h0 1 + ay ly

)

is considered. In comparison with those by Eisenberger [46] and Farajpour [40], the nondimensional buckling load is calculated as:

λ 2 = N 0l y2 / π 2 D110

(46)

The results are listed in Table 3 compared to those by Eisenberger [46]. In this Table, the Winkler and shear factors are ignored. From the table, it is seen that the results by the present method are relatively equal to those by Eisenberger [46]. It is also observed that the effect of thickness on λ2 becomes stronger when the coefficient a increases.

Next, the non-dimensional buckling load value λ2 of the first four modes is investigated for a square nanoplate. The thickness varies in the y-axis with λ ( x, y ) = 1 + y 4l y . The plate is loaded by uniaxial compressive force in the x-axis, µ = 0 . The length of the square nanoplate is assumed to be 5, 7.5, 10, 12.5, 15, 17.5 and 20 nm and the scale coefficient is equal to 2 nm. The results of the present method are shown in Figure 15 in comparison with those by Farajpour [40]. It is seen that when the value of plate’s length increases, the nondimensional buckling load also increases. The non-dimensional buckling load λ2 is higher at higher modes. It can also be observed from Figure 16 that when the non-dimensional scale coefficient increases, the non-dimensional buckling load λ2 decreases for all modes. In addition, when the values of scale coefficient are greater than or equal to 2 ( e0a ≥ 2) , the nondimensional buckling load depends on the mode numbers. The effects of the aspect ratio in conjunction with the percentage change (100a%) in the non-dimensional buckling load λ2 are presented in Figure 17. The results show that the non-dimensional buckling load is proportional to the increase of the percentage change. This effect is stronger in the case of lower aspect ratio. It also can be seen from the figure that the results by the IGA are closely equal to those by Farajpour et al. [40]. This again strengthens the robustness of present method in solving buckling problem for nanoplates. In addition, it’s worth noting that the continuity order of the IGA could be controlled by knot vectors. As a result, the IGA has high potential to be extended to problems related to high-ordered differential equations. 4.5 Nanoplates with linear thickness variation is embedded in elastic medium In order to further illustrate the capability of the present method, a thin orthotropic graphene plate with inconstant thickness is next considered. The thickness variation is given by the rule

h = 1+ ay ly , ( h0 = 1) in the x-axis. Figure 18 depicts the critical buckling load ratio of a biaxial compressed square graphene sheet with length = 15 nm, scale coefficient e0a = 2 and the shear modulus parameter K S = 0. From the table, it is seen that the effect of KW on buckling ratio is relatively negligible for the clamped plate but it becomes considerable in the simply supported case. In general, the buckling ratio of simply supported plates is bigger than that of clamped plates regardless of Winkler modulus parameters.

In Figure 19, the effects of the non-dimensional shear modulus parameter K S on fully simply supported square graphene sheet is presented. The graphene sheet is embedded in an elastic medium with the Winkler modulus parameter K W = 100 . The variation in thickness

h

is given by h = h0 (1+ 0.25y / ly) with h0 = 0.34 nm. The plate is subjected to biaxial compressive forces in the x, y-directions. The length is taken as L=15 nm. It can be seen that when the scale coefficient increases, the impact of the shear parameter K S on critical load ratio also increases. Figure 20 illustrates the influence of KW on the buckling load with respect to different values of scale coefficient. In this case, the value of non-dimensional K S is fixed at 10 while the non-dimensional Winkler modulus parameter is taken from 0 to 500. Similarly, when the scale coefficient increases, the effect of the shear parameter K w on critical load ratio increases. When the scale coefficient is close to 0, the impact of K w becomes unconsiderable. The effect of the parameter a on the critical load ratio is next demonstrated for simply 2

supported SLGS under biaxial compression. The scale coefficient is taken as ( e0a ) = 0.0, 0.5, 0.1, 1.5 and 2.0 nm. The stiffness is characterized by K S = 0 for the Pasternak foundation model and by K W = 100 for the Winkler foundation model. Minimum thickness

( h0 ) and dimension of SLGS are taken as 0.34 and 15, respectively. The results are presented in Figure 21. It is seen that the higher the scale coefficient value is, the smaller the buckling load ratio becomes.

5. Conclusion

The paper contributes a novel numerical approach for buckling analysis of nanoplates subjected to in-plane loadings using IGA in association with the nonlocal theory. The graphene sheets in the paper are modeled as nanoplates which comply with classical plate theory. The small scale effects (i.e influence of small distance, connective forces between atoms) are taken into account and presented in the relationship of the plates’ stress and strain. Numerical performance is conducted to investigate the influences of the scale coefficient, the Winkler modulus, the shear modulus of the elastic medium, plate dimensions and aspect ratio on buckling behavior of nanoplates. The numerical results by the present

method agree well with those by previous publications for different cases. Correlative relations between small scale effects (i.e scale coefficients, elastic medium parameters, plates dimensions, aspect ratio, etc.) and buckling load found in the present study are also similar to those by previous publications. This shows that IGA is suitable and robust for modelling buckling behavior of the nanoplates. In addition, because the construction of high-ordered continuous basis function could be easily controlled through knot vector, IGA has high potential to solve more complex problems related to high-ordered partial differential equations with high accuracy. This opens a promising future for IGA in analyzing mechanical behaviors of structures at micro and nano scales.

Acknowledgements

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.99-2014.11.

References

[1]

Iijima S. Helical microtubules of graphitic carbon. Nature 1991;354:56–8.

[2]

Lu P, Lee HP, Lu C, Zhang PQ. Dynamic properties of flexural beams using a nonlocal elasticity model. J Appl Phys 2006;99.

[3]

Aydogdu M. A general nonlocal beam theory: Its application to nanobeam bending, buckling and vibration. Phys E Low-Dimensional Syst Nanostructures 2009;41:1651– 5.

[4]

Liu J, Ma Z, Wright AR, Zhang C. Orbital magnetization of graphene and graphene nanoribbons. J Appl Phys 2008;103.

[5]

Duan WH, Wang CM. Exact solutions for axisymmetric bending of micro/nanoscale circular plates based on nonlocal plate theory. Nanotechnology 2007;18:385704.

[6]

Wang Q, Wang CM. The constitutive relation and small scale parameter of nonlocal continuum mechanics for modelling carbon nanotubes. Nanotechnology 2007;18:75702.

[7]

Sakhaee-Pour A, Ahmadian MT, Naghdabadi R. Vibrational analysis of single-layered graphene sheets. Nanotechnology 2008;19:85702.

[8]

Li C, Chou T-W. Mass detection using carbon nanotube-based nanomechanical resonators. Appl Phys Lett 2004;84.

[9]

Sakhaee-Pour A, Ahmadian MT, Vafai A. Applications of single-layered graphene sheets as mass sensors and atomistic dust detectors. Solid State Commun 2008;145:168–72.

[10] Dai H, Hafner JH, Rinzler AG, Colbert DT, Smalley RE. Nanotubes as nanoprobes in

scanning probe microscopy. Nature 1996;384:147–50. [11] Zhou S, Li Z. Length scales in the static and dynamic torsion of a circular cylindrical micro - bar. Jounral Shandong Univ Technol 2001;31:401–7. [12] Fleck N, Hutchinson J. Strain gradient plasticity. Adv Appl Mech 1997;33:295–361. [13] Yang F, Chong ACM, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39:2731–43. [14] Eringen AC. Nonlocal polar elastic continua. Int J Eng Sci 1972;10:1–16. [15] Eringen AC. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. J Appl Phys 1983;54. [16] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194:4135–95. [17] Bazilevs Y, Calo VM, Hughes TJR, Zhang Y. Isogeometric fluid-structure interaction: theory, algorithms, and computations. Comput Mech 2008;43:3–37. [18] Bazilevs Y, Gohean JR, Hughes TJR, Moser RD, Zhang Y. Patient-specific isogeometric fluid–structure interaction analysis of thoracic aortic blood flow due to implantation of the Jarvik 2000 left ventricular assist device. Comput Methods Appl Mech Eng 2009;198:3534–50. [19] Gomez H, Hughes TJR, Nogueira X, Calo VM. Isogeometric analysis of the isothermal Navier–Stokes–Korteweg equations. Comput Methods Appl Mech Eng 2010;199:1828–40. [20] Bazilevs Y, Akkerman I. Large eddy simulation of turbulent Taylor–Couette flow using isogeometric analysis and the residual-based variational multiscale method. J Comput Phys 2010;229:3402–14. [21] Bazilevs Y, Michler C, Calo VM, Hughes TJR. Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes. Comput Methods Appl Mech Eng 2010;199:780–90. [22] Wall WA, Frenzel MA, Cyron C. Isogeometric structural shape optimization. Comput Methods Appl Mech Eng 2008;197:2976–88. [23] Qian X. Full analytical sensitivities in NURBS based isogeometric shape optimization. Comput Methods Appl Mech Eng 2010;199:2059–71. [24] Buffa A, Sangalli G, Vázquez R. Isogeometric analysis in electromagnetics: B-splines approximation. Comput Methods Appl Mech Eng 2010;199:1143–52. [25] Manalo AC. Behaviour of fibre composite sandwich structures under short and asymmetrical beam shear tests. Compos Struct 2013;99:339–49. [26] Nguyen-Thanh N, Valizadeh N, Nguyen MN, Nguyen-Xuan H, Zhuang X, Areias P, et al. An extended isogeometric thin shell analysis based on Kirchhoff–Love theory. Comput Methods Appl Mech Eng 2015;284:265–91. [27] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger KU, Bazilevs Y, et al. Rotation free isogeometric thin shell analysis using PHT-splines. Comput Methods Appl Mech Eng 2011;200:3410–24. [28] Shojaee S, Izadpanah E, Valizadeh N, Kiendl J. Free vibration analysis of thin plates by using a NURBS-based isogeometric approach. Finite Elem Anal Des 2012;61:23– 34.

[29] Beirão da Veiga L, Buffa A, Lovadina C, Martinelli M, Sangalli G. An isogeometric method for the Reissner–Mindlin plate bending problem. Comput Methods Appl Mech Eng 2012;209–212:45–53. [30] Valizadeh N, Natarajan S, Gonzalez-Estrada OA, Rabczuk T, Bui TQ, Bordas SPA. NURBS-based finite element analysis of functionally graded plates: Static bending, vibration, buckling and flutter. Compos Struct 2013;99:309–26. [31] Tran L V, Thai CH, Nguyen-Xuan H. An isogeometric finite element formulation for thermal buckling analysis of functionally graded plates. Finite Elem Anal Des 2013;73:65–76. [32] Thai CH, Ferreira AJM, Bordas SPA, Rabczuk T, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. Eur J Mech - A/Solids 2014;43:89–108. [33] Thai CH, Nguyen-Xuan H, Bordas SPA, Nguyen-Thanh N, Rabczuk T. Isogeometric Analysis of Laminated Composite Plates Using the Higher-Order Shear Deformation Theory. Mech Adv Mater Struct 2015;22:451–69. [34] Thai CH, Ferreira AJM, Carrera E, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a layerwise deformation theory. Compos Struct 2013;104:196–214. [35] Guo Y, Nagy AP, Gürdal Z. A layerwise theory for laminated composites in the framework of isogeometric analysis. Compos Struct 2014;107:447–57. [36] Pradhan SC, Phadikar JK. Small scale effect on vibration of embedded multilayered graphene sheets based on nonlocal continuum models. Phys Lett A 2009;373:1062–9. [37] Pradhan SC, Murmu T. Small scale effect on the buckling analysis of single-layered graphene sheet embedded in an elastic medium based on nonlocal plate theory. Phys E Low-Dimensional Syst Nanostructures 2010;42:1293–301. [38] Karamooz Ravari MR, Talebi S, Shahidi AR. Analysis of the buckling of rectangular nanoplates by use of finite-difference method. Meccanica 2014;49:1443–55. [39] Ansari R, Sahmani S, Arash B. Nonlocal plate model for free vibrations of singlelayered graphene sheets. Phys Lett A 2010;375:53–62. [40] Farajpour A, Danesh M, Mohammadi M. Buckling analysis of variable thickness nanoplates using nonlocal continuum mechanics. Phys E Low-Dimensional Syst Nanostructures 2011;44:719–27. [41] Winkler E. Theory of elasticity and strength. Dominicus Prague, Czechoslov 1867. [42] Pasternak PL. On a new method of analysis of an elastic foundation by means of two foundation constants. Gos Izd Lit Po Stroit I Arkhitekture, Moscow 1954. [43] Lu P, Zhang PQ, Lee HP, Wang CM, Reddy JN. Non-local elastic plate theories. Proc R Soc London A Math Phys Eng Sci 2007;463:3225–40. [44] Liew KM, He XQ, Kitipornchai S. Predicting nanovibration of multi-layered graphene sheets embedded in an elastic matrix. Acta Mater 2006;54:4229–36. [45] Pradhan SC, Murmu T. Small scale effect on the buckling of single-layered graphene sheets under biaxial compression via nonlocal continuum mechanics. Comput Mater Sci 2009;47:268–74. [46] Eisenberger M, Alexandrov A. Buckling loads of variable thickness thin isotropic plates. Thin-Walled Struct 2003;41:871–89.

Appendix 1. The parameters in nanoplate equation A1 =

( e0a )

2

ks + λ 3D110

(1.1)

A2 =

( e0a )

2

ks + λ 3D220

(1.2)

A3 =

−6D110 λ 2λ,x

(1.3)

A4 =

−6D220 λ2λ, y

(1.4)

A5 =

( e0 a )

2

k w + k s − 3D110 ( 2 λλ,2x + λ 2λ, xx ) − 3 D120 ( 2 λλ,2y + λ 2λ, yy )

(1.5)

A6 =

( e0 a )

2

0 k w + k s − 3D120 ( 2 λλ,2x + λ 2λ, xx ) − 3D22 ( 2λλ,2y + λ 2λ, yy )

(1.6)

A7 =

( e0 a )

2

0 k s + λ 3 ( D120 + 2 D 33 )

(1.7)

A8 =

−6λ 2λ, y ( D120 + 2 D330 )

(1.8)

A9 =

−6λ 2λ, x ( D120 + 2D330 )

(1.9)

A10 =

−6D330 ( 2λλ,xλ, y + λ 2λ,xy )

(1.10)

A11 =

kw

(1.11)

Appendix 2. The matrix components  N I , xx   N I , xx  N    I , yy    N I , yy  B =  N I , x  , BG =   N   NI,x   I ,y   N I , y   N I 

 A1 A  7 D=0  0  0

A7

A3

A8

A2

A9

A4

0 0

A5 A10

A10 A6

0

0

0

 1 2  ( e0 a )  (1 + µ ) / 2 DG =   [0]2×2 

 x,ξ x  ,η  T =  x,ξξ   x,ηη x  ,ξη

y,ξ y,η

0 0

0 0 2

2

y,ηη

,η

,η

y,ξη

x,ξ x,η

y,ξ y,η

,ξ

,ξ

2

 

(2.2)

    1 0   0 µ    

[0]2×2

   2 x,ξ y,ξ    2 x,η y,η  y,ξ x,η + y,η x,ξ 

(2.3)

0 0

(x ) ( y ) (x ) ( y )

y,ξξ

0 0  0   0  A11 

(1 + µ ) / 2  µ

(2.1)

2

(2.4)

TABLES Table 1. Buckling load ratio of simply-supported square nanoplate under uniaxial compression Length (nm)

e0 a

(nm)

5

10

15

20

25

30

35

40

45

0.0

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

0.5

0.8351

0.9530

0.9785

0.9878

0.9922

0.9945

0.9960

0.9969

0.9976

1.0

0.5254

0.8351

0.9193

0.9530

0.9694

0.9785

0.9841

0.9878

0.9903

1.5

0.2805

0.6925

0.8351

0.9001

0.9337

0.9530

0.9650

0.9730

0.9785

2.0

0.1587

0.5254

0.7402

0.8351

0.8878

0.9193

0.9395

0.9530

0.9625

Table 2. Buckling load ratio of simply-supported square nanoplate under biaxial load e0 a (nm)

Length (nm) 5

10

25

0.0

1.00

1.00

1.00

0.5

0.8352

0.9530

0.9922

1.0

0.5588

0.8352

0.9694

1.5

0.3602

0.6925

0.9337

2.0

0.2405

0.5588

0.8878

Table 3. The non-dimensional buckling load values of uniformly compressed isotropic plate with variable thickness in y-axis: lx l y

0.5

0.7

0.9

1.1

1.3 1.5

1.6

2.0

2.4

2.8

Method

a

0.125

0.25

0.5

0.75

1.0

Present solution

7.4621

8.7531

11.5687

14.6953

18.1368

Eisenberger [46]

7.4645

8.7633

11.6112

14.7942

18.3175

Present solution

5.4194

6.3869

8.5627

11.0657

13.9017

Eisenberger [46]

5.4199

6.3891

8.5741

11.0979

13.9730

Present solution

4.8428

5.7224

7.7327

10.0858

12.7877

Eisenberger [46]

4.8413

5.7165

7.7111

10.0460

12.7381

Present solution

4.8363

5.7242

7.7729

10.1928

12.9874

Eisenberger [46]

4.8327

5.7097

7.7163

10.0765

12.8101

Present solution

5.1324

6.0814

8.2833

10.8964

13.9201

Eisenberger [46]

5.1266

6.0574

8.1885

10.6973

13.6053

Present solution

5.1931

6.1249

8.2311

10.6672

13.4394

Eisenberger [46]

5.1931

6.1252

8.2344

10.6817

13.4811

Present solution

5.0296

5.9363

7.9941

10.3853

13.1163

Eisenberger [46]

5.0291

5.9345

7.9891

10.3820

13.1280

Present solution

4.7912

5.6667

7.6782

10.0452

12.7726

Eisenberger [46]

4.7887

5.6565

7.6394

9.9678

12.6602

Present solution

4.9549

5.8682

7.9820

10.3858

13.1169

Eisenberger [46]

4.9502

5.8490

7.9067

10.3289

13.1362

Present solution

4.8131

5.6891

7.6953

10.0481

12.7533

Eisenberger [46]

4.8864

5.7686

7.7771

10.1247

12.8279

FIGURES

L

b

(a)

(b)

Figure 1: (a) Graphene sheet modeled as a continuum nanoplates. (b) Graphene sheet embedded in elastic medium is represented by Pasternak model.

Figure 2: Rectangular SLGS with non-uniform thickness subjected to biaxial compression.

1

1

0.8 0.6

0.5 0.4

0 0

0.2

0 0.5

0 0

0.2

0.4

0.6

0.8

1

0.5 1 1

Figure 3: B-spline basis functions in 1D and 2D corresponding knot vectors.

(a)

(b)

Figure 4: Nanoplates subjected to uniaxial (a) and biaxial (b) compressions.

Figure 5: Effect of length on the critical buckling load ratio of a simly supported SLGS for different small scale coefficients.

(a) Fully simply supported plate

(b) Fully clamped plate

Figure 6: The change of non-dimensional buckling load λ1 with non-dimensional nonlocal variable ψ .

Figure 7: The change of non-dimensional buckling load λ1 with non-dimensional nonlocal variable ψ for the different boundary condition cases.

(a) L = 5 nm

(b) L = 10 nm

(c)

L = 25 nm

Figure 8: The effects of scale coefficient on axial and biaxial compressions.

Figure 9: The effects of aspect ratio on load ratio λ for uniaxial and biaxial compressions.

Figure 10: The effects of the polynomial order and the number of the elements on critical load ratio.

Figure 11: Effect of Winkler modulus

Figure 12 The effect of the polynomial order and

parameter on critical buckling load for simply

the numbering of the element with different

supported square nanoplate.

Winkler modulus parameters.

Figure 13: Effect of shear modulus parameter on critical buckling load for simply supported square nanoplate.

Figure 14: Effect of aspect ratio on critical buckling load for simply supported square nanoplate.

Figure 15: Effect of length on critical buckling load λ2 for the first four mode numbers.

Figure 16: Effect of scale coefficient on critical buckling load λ2 for the first four mode number.

Figure 17: Effect of aspect ratio on critical buckling load for simply supported edge boundary conditions.

Figure 18: The change of critical load ratio λ with respect to the parameter a in different boundary conditions.

Figure 19: Effect of the shear parameter on Figure 20: Effect of Winkler modulus parameter critical buckling load with Winkler modulus on critical buckling load with the shear parameter parameter K W = 100 .

K S = 10 .

Figure 21: Effect of parameter a on critical load ratio for simply supported edge boundary 2

conditions with square scale effect ( e0 a ) increased from 0 to 2.