Naturally stabilized nodal integration meshfree formulations for analysis of laminated composite and sandwich plates

Accepted Manuscript Naturally stabilized nodal integration meshfree formulations for analysis of laminated composite and sandwich plates Chien H. Thai, A.J.M. Ferreira, H. Nguyen-Xuan PII: DOI: Reference:

S0263-8223(17)31203-5 http://dx.doi.org/10.1016/j.compstruct.2017.06.049 COST 8636

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

15 April 2017 13 June 2017 20 June 2017

Please cite this article as: Thai, C.H., Ferreira, A.J.M., Nguyen-Xuan, H., Naturally stabilized nodal integration meshfree formulations for analysis of laminated composite and sandwich plates, Composite Structures (2017), doi: http://dx.doi.org/10.1016/j.compstruct.2017.06.049

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.

Naturally stabilized nodal integration meshfree formulations for analysis of laminated composite and sandwich plates Chien H. Thai1,2*, A.J.M. Ferreira3, H. Nguyen-Xuan4,5* 1

Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Vietnam 2

3

Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam

Departamento de Engenharia Mecanica, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal 4

5

Institute of Research and Development, Duy Tan University, Da Nang, Vietnam

Department of Architectural Engineering, Sejong University, 98 Gunja-dong, Gwangjin-gu, Seoul 143-747,

South Korea

ABSTRACT This article addresses a naturally stabilized nodal integration (NSNI) meshfree formulation for static, free vibration and buckling analyses of laminated composite and sandwich plates based on the higher order shear deformation theory (HSDT). The crucial idea is to directly compute integrations at nodes similar to the direct nodal integration (DNI) but the instability existing in DNI is avoided and the computation cost is significantly reduced when compared to traditional high-order Gauss quadrature schemes. Due to the displacements, strains and stresses directly computed at nodes, computation time of the present approach is equivalent with the DNI. Being different from the stabilized conforming nodal integration technique using a divergence theorem to evaluate gradient strains by the surface-to-boundary integration, NSNI is a naturally implicit gradient expansion. Variational consistency in Galerkin weak form is explicitly described. Thanks to the satisfied Kronecker delta function property, the moving Kriging integration (MKI) shape function is used to impose directly essential boundary conditions as same as the finite element method (FEM). A variety of numerical examples with various complex geometries, aspect ratios, stiffness ratios and boundary conditions are studied. Numerical results show the effectiveness of the present method in comparison with some existing ones in the literature. *

Corresponding author at: Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam

(Chien H. Thai, Email: [email protected]) and Institute of Research and Development, Duy Tan University, Da Nang, Vietnam (H. Nguyen-Xuan, Email: [email protected]).

1

Keywords: Naturally stabilized nodal integration; moving Kriging interpolation; meshfree method; laminated composite and sandwich plates. 1. Introduction Composite materials are widely used in various engineering applications due to their excellent features in the high stiffness and strength-to-weight ratios, long fatigue life, wear resistance, lightweight, etc… In engineering applications, these materials are used in beam, plate and shell structures. Therefore, analysis of its behaviors on displacements, stresses, free vibration frequency and buckling loading factors are always required. It has been found in the literature that, computational theories addressing to composite and sandwich plate structures can be generally divided into three groups according to the continuum-based three-dimensional (3D) elasticity theory, the discrete layer theory based on 2D computational model and the equivalent single layer (ESL) theory based on 2D computational model. The 3D elasticity theory was assumed that each lamina is taken as a complete 3D solid. It can provide the actual behavior of plates but takes the high computational cost. In addition, the finding of an exact solution for problems with complex geometries and arbitrary boundary conditions is not a feasible task. The 3D elasticity solutions for several simple cases of composite plates were investigated in [1-4]. The second one considering as the layer-wise/zigzag theory is that each layer is taken into account in the analysis. Thus, the number of approximation variables depends on the number of layers and the computational cost is very significant when increased the number of layers. Some contributions of the discrete layer theories have been published in [5-7]. For practical applications of the 3D elasticity theory and the discrete layer theory, some literature reviews [8-11] have been clearly presented. The last one was assumed that all laminated layers are unified as a homogeneous single layer. Hence, the number of unknowns is independent of the number of layers. The ESL theories include the classical laminated plate theory (CLPT), the first order shear deformation theory (FSDT) [12-14] and the higher order shear deformation theory (HSDT). The CLPT leads to inaccurate results for laminated composite plates because it is not taken into account for effects of transverse shear strains. The FSDT needs a shear correction factor to match the strain energy derived from FSDT solution with those from 3D elasticity solution because of assumptions of constantly transverse shear stresses through the plate thickness. In computation, it is difficult to determine an optimized value for shear correction factor because of depending on the material properties, geometries and boundary conditions of problems [15]. The HSDT models involve the third order shear deformation theory (TSDT) [16-19], the fifth

2

order shear deformation theory (FiSDT) [20,21], the seventh order shear deformation theory [22], the trigonometric shear deformation theory [23-26], the inverse trigonometric shear deformation theory [27], the exponential shear deformation theory (ESDT) [28,29], and so on. It was clearly observed that almost these higher order shear deformation theories request the C1-continuity and involve slope components related to the derivation of the transverse displacement. Numerical methods as the standard finite element method, the meshfree method are therefore difficult to impose directly boundary condition for these slope components. In this paper, we address a C0-type higher order shear deformation model. Meshfree methods based on the Galerkin weak form have been applied and developed for analysis of laminated composite and sandwich plates. A main problem of meshfree methods with the Galerkin weak form is the choice of efficient quadrature techniques. It is very difficult to integrate exactly for rational shape functions in meshfree methods. In particular, two types of integration based on the Galerkin weak are the Gauss quadrature and nodal integrations. The first one requires background cells to compute the numerical quadrature integrations. Generally, the traditional high-order Gauss quadrature scheme is employed to achieve stable and accurate solutions but it takes very much computation time. Using the lower-order quadrature retains less CPU times but the numerical solutions can be not converged and stabilized [30,31]. The last one related to the nodal integration is found to be an appropriate choice. In the directly nodal integration technique (DNI), gradients are evaluated directly at the nodes. It was found in this DNI that the underlying solution is unstable due to leading to rank deficiency of stiffness matrices. An improved version of the DNI is known as the stabilized conforming nodal integration (SCNI) technique proposed by Chen et al. [30]. This technique used the smoothing operation to convert to boundary integrations according to the divergence theorem. Therefore, gradients are not directly evaluated at the nodes but through vertices/nodes of conforming representative nodal domains or Voronoi diagrams, which are not defined as approximate variables. A new version of the SCNI is called a modified SCNI (MSCNI) by an additional stabilization term. This technique has been applied to various methods and problems [32, 33]. In other way, the least-squares method was presented to alleviate the instability in DNI, as shown in [34,35]. However, this method contains second order derivatives of approximation functions and requires an additional stabilization parameter to guarantee convergent and stable solutions. In addition, Taylor series expansions of the displacement fields [36,37] are one of ways to reduce the instability in DNI, but high order derivatives are included in underlying formulations resulting in increasing the computational cost. Furthermore, other form of stabilization consisting of the Taylor expansion and a displacement smoothing was

3

introduced, as described in [38, 39]. The nodal integration technique is used for this method yielding stable solutions. Recently, Hillman and Chen [33] proposed a naturally stabilized nodal integration (NSNI) technique which is a naturally implicit gradient expansion and does not involve high order derivatives in approximate formulations. More details for developments and applications of meshfree methods in over 20 years can be found in [40]. In this study, we exploit the advantages of the NSNI for meshfree method based on moving Kriging integration (MKI) shape functions for laminated composite and sandwich plate analysis. The outline of the paper is organized as follows. In section 2, basic equations of laminated composite plate are introduced. In section 3, the NSNI for the laminated plates is presented. Section 4 defines geometries and essential boundary conditions. Numerical results are given in Section 5. Section 6 concludes some remarks. 2. Basic equations 2.1 Kinematics of laminated composite plates

 h h We consider a plate carrying a domain V = Ω ×  − ,  , where Ω∈ R 2 is the mid-plane of the  2 2 plate and h is the plate thickness. According to the generally higher order shear deformation theory [29], the displacement field of any points in the plate can be described by u ( x, y , z ) = u 0 ( x, y ) + zu1 ( x, y ) + f ( z )u 2 ( x, y )

(1)

where

u   u0   w0, x  θ x    0     2   1 u =  v ; u =  v0  ; u = − w0, y ; u = θ y   w w   0  0    0    

(2)

in which u0 , v0 , w0 , θ x and θ y are the in-plane, transverse displacements and the rotation components in the y-z, x-z planes, respectively, ‘,x’ or ‘,y’ defines the derivation of any functions following x or y direction and f ( z ) is any function depending to the z-coordinate.

Compatible strain fields which require C1 continuity of the transverse displacement are obtained from Eq. (1) by using the displacement-strain relationship. To avoid the usual difficulties associated with C1 continuity requirement together with boundary conditions for the slope components related to

4

the derivation of displacements ( w0, x and w0, y ), which are not defined as the approximation variables, the following additional assumptions are enforced w0, x = β x and w0, y = β y

(3)

Substituting Eq. (3) into Eq. (2), the displacement field based on C0-type higher order shear deformation theory is defined as

 u0  βx  θ x      2   1 u =  v0 ; u = − β y  ; u = θ y  w  0 0  0     0

(4)

From Eq. (1), the transverse displacement does not depend on the coordinates of the thickness direction so that the transverse normal stress and strain are assumed to be equal zero. Bending and shear strains according to the C0-type higher order shear deformation theory are calculated as follows ε = {ε xx

T

ε yy γ xy } = ε 0 + z ε1 + f ( z )ε 2 and γ = {γ xz

T

γ yz } = ε s 0 + f ′( z )ε s1

(5)

where  θ x, x   u0, x   β x, x     1   2  ε =  v0, y  ; ε = −  β y , y  ; ε =  θ y , y  θ + θ  u + v  β + β  y, x  y, x  0, x   0, y  x, y  x, y 0

(6)

 w0, x − β x  s1 θ x  εs 0 =  ; ε =   w0, y − β y  θ y 

in which f ′( z ) is the derivation of the function f ( z ) . In this study, two functions across the plate thickness including to TSDT and ITSDT are investigated, as shown in Table 1. Table 1: Two functions and derivations. Model TSDT [17]

ITSDT [26]

f ′( z )

f ( z) f (z) = z −

4 z3 3h 2

f ′( z ) = 1 −

4z2 h2

 4z2  f ′ ( z ) = 2 /  2 + 1 − 1  h  Using the Hooke law in the case of plane stress, constitutive equations of each orthotropic layer in  2z  f ( z ) = h tan −1   − z  h 

the local coordinate system are defined by

5

σ 1( k )  Q  ( k )   11 Q12 σ 2  Q21 Q22  (k )   0 τ 12  =  0  (k )   0 0 τ 13   0  ( k )   0 τ 23 

0 0

0 0

Q66

0

0 Q55

0

0

0  0  0   0  Q44 

(k )

ε 1( k )   (k )  ε 2   (k )  γ 12   (k )  γ 13   (k )  γ 23 

(7)

where, subscripts 1 and 2 are the directions of the fiber and in-plane normal to fiber, respectively while subscript 3 indicates the direction normal to the plate. Qij( k ) is given as follows Q11( k ) =

E1( k )

, Q12( k ) = (k )

1 −ν 12( k )ν 21

ν 12( k ) E2( k ) E2( k ) (k ) (k ) (k ) (k) Q , = , Q66 = G12( k ) , Q55( k ) = G13( k ) , Q44 = G23 22 (k ) (k ) (k ) (k ) 1 −ν 12 ν 21 1 − ν 12 ν 21

(8)

(k ) where E1( k ) , E2(k ) are the Young modulus in the 1 and 2 directions, respectively, and G12( k ) , G23 , G13( k ) (k ) are the shear modulus in the 1-2, 2-3 and 1-3 planes, respectively, and ν12( k ) and ν 21 are Poisson’s

ratios. In practical application, several orthotropic layers are bonded together with a laminate and each layer must be transformed to the laminate coordinate system. The stress-strain relationship is defined by σ xx( k )  Q11( k )  (k )   (k ) σ yy  Q21  (k )   (k ) τ xy  = Q61 τ ( k )   0  xz( k )   τ yz   0

Q12( k )

Q16( k )

0

(k ) Q22

Q26( k )

0

(k ) 62

(k ) 66

0

0 0

(k ) 55 (k ) 45

Q

Q

0 0

Q Q

(k ) 0  ε xx    (k )  0  ε yy   (k )  0  γ xy   (k ) Q54( k )  γ xz   (k )  Q44( k )  γ yz 

(9)

where Qij( k ) is transformed material constant. A detailed presentation can be found in [41]. The essential equations for bending analysis of the laminated composite plates under a transverse loading q0 are obtained by the Galerkin weak form as follows T

ε 0    δ  ε1  Ω  2 ε 

∫

 Ab  b B E 

E  ε 0     Db F1  ε1  dΩ + F1 H  ε 2    Bb

T

s0 s ε   A δ  s1   s Ω ε   B

∫

Bs  ε s 0     dΩ = Ds  ε s1 

∫ δ w q dΩ

(10)

( z ) ) Qijk dz where (i,j =1,2,6)

(11)

Ω

0 0

where

( A , B , D , E , F1 , H ) b ij

b ij

b ij

ij

ij

ij

(k )

=∫

h/ 2

− h /2

(1, z, z , f ( z), z f ( z), f 2

6

2

s (k ) ij

(A ,B ,D ) s ij

s ij

=∫

h/ 2

(1, f ′( z ), f ′ ( z ) ) Q dz where (i,j =4,5) 2

− h/ 2

k ij

Similarly, the dynamic version of the Galerkin weak form for the laminated composite plates can also be expressed by T

ε 0    δ  ε1  Ω  2 ε 

∫

u 0    δ  u1  Ω  2 u 

T

∫

 Ab  b B E 

Bb

E  ε 0     Db F1  ε1  dΩ + F1 H  ε 2   

 I1 I  2  I 4

I2 I3 I5

T

s0 s ε   A δ  s1   s Ω  ε   B

∫

B s  ε s 0     dΩ = ... Ds   ε s1  (12)

u0  I 4  &&   1  &&  dΩ I 5   u 2 I 6  u && 

where h/ 2

( I1 , I 2 , I3 , I 4 , I5 , I6 ) = ∫ ρ ( z ) (1, z, z 2 , f ( z), zf ( z ), f 2 ( z ))dz

(13)

− h/ 2

For linear buckling analysis, the Galerkin weak form for the laminated composite plates under an in-plane loading can be described as T

ε 0    δ  ε1  Ω  2 ε 

∫

 Ab  b B E  T

 w0, x  h δ  Ω  w0, y   

∫

E  ε 0     Db F1  ε1  dΩ + F1 H  ε 2    Bb

T

s0 ε  δ  s1  Ω  ε 

∫

 As  s  B

Bs  ε s 0     dΩ + ... Ds   ε s1 

(14)

0  σ x0 σ xy   w0, x   0   dΩ = 0 0 w σ xy σ y   0, y 

0 where σ x0 , σ 0y and σ xy are the pre-buckling loadings corresponding with x, y and x-y directions.

2.2 Moving Kriging interpolation shape functions The moving Kriging interpolation shape functions are chosen to interpolate kinematic variables in laminated composite plates. The domain Ω is discretized by a set of nodes xI (I = 1, . . . , NP), in which NP is the quantity of the nodes. The moving Kriging interpolation of a variable u ( x) , defined by u h ( x) , can be written as follows NP

u h ( x ) = ∑ φI ( x)uI I =1

where φI and u I are the shape function and the unknown coefficient associated with node I, and

7

(15)

m

n

j =1

k =1

φI ( x) = ∑ p j ( x)AjI + ∑ rk (x)BkI or φI ( x) = p T ( x) A + r T ( x)B

(16)

in which n and m are a number of nodes in a support domain Ω x ∈ Ω (see Figure 1) and the order of polynomial basis function, respectively. Four terms of A , B , p( x) and r( x) are defined as follows

A = ( P T R -1P )-1 P T R -1 , B = R -1 (I − PA) p ( x ) = {1 x

x2

y

xy

(17)

T

y 2 ...} and r(x) =[ R(x1,x) R(x2, x) L R(xn, x)]T

where I is a unit matrix with the size n × n . For two-dimensional problem, the polynomial basis function of quadratic form is chosen as p ( x ) = {1 x

y

x2

T

y 2 } (m=6)

xy

(18)

In addition, two matrices P and R are expressed by  p1 ( x1 ) ... ... P =  ...  p1 (x n ) ...

pm ( x1 )   R( x1 , x1 ) ... R( x1 , x n )   ...  and R =  ... ... ...   R (x n , x1 ) ... R (x n , x n )  pm ( x n ) 

(19)

where R (x I , x J ) is a correlation function and is defined by R (x I , x J ) =

2 1  h E ( u (x I ) − u h (x J ) )   2 

(20)

in which E defines an expected value of a random function. The correlation function R (x I , x J ) as multi-quadrics, thin plate splines, Gaussian, etc… can be chosen to build MKI shape function. In this paper, we use the Gaussian function described by R (x I , x J ) = e

 βr −  IJ  a0

  

2

(21)

where rIJ = x I - x J , the correlation parameter β is related to the variance σ 2 of the normal distribution function by β 2 = 1 / 2σ 2 and a0 is the scale factor used to normalize the distance. The maximum distance between a pair of nodes in the support domain is chosen for this value a0 . Following previous studies [42-44], the correlation function with a normalized distance is stable and is not influenced by the correlation parameter β . For this reason, the correlation parameter is fixed at 1 ( β = 1 ) for this study.

8

Figure 1. Domain representation and support domain of 2D model. The first-order derivative of the MKI shape functions is written as follows m

n

m

n

j =1

k =1

j =1

k =1

φI , x ( x) = ∑ p j , x ( x)AjI + ∑ rk , x (x)BkI and φI , y (x) = ∑ p j , y (x) AjI + ∑ rk , y ( x)BkI

(22)

For meshfree methods, there is a support domain (influence domain) to choose a set of nodes to build shape functions. In this paper, a circular influence domain with a support size is defined as dm = α dc

(23)

where dc denotes an average distance between nodes and α is a scale factor. In numerical computation, the radius of influence domain is enough large to cover a sufficient number of nodes for constructing shape functions. In this study, a circular support domain with a radius fixed at 2.4 times nodal spacing ( α = 2.4 ) is used in the numerical computation. However, this factor can be chosen differently between 2.0 and 4.0 [45]. 2.3 Nodal integration formulations for laminated composite plates based on the C0-type HSDT The weak forms in Eqs. (10), (12) and (14) require numerical quadrature, e.g. Gauss quadrature. Generally speaking, the shape functions in meshfree methods are rational functions of the spatial coordinates. Hence, the high-order Gauss quadrature scheme is used to compute stiffness, mass, geometry matrices and force vector. As known, such a quadrature scheme consumes much computational time. An appropriate choice related to the nodal integration has been therefore recommended. It is so-called direct nodal integration (DNI) in the literature because shape functions are directly evaluated at the nodes. For nodal integration technique, the domain Ω is divided into a set of non-overlapping sub-domains Ω I (I=1, 2,..., NP), where each of which includes a node NP

( Ω = ∑ Ω I ), as shown in Figure 2. For any distributed nodes, a Voronoi diagram can be generated I =1

9

automatically. The weak forms for static, free vibration and buckling analyses can be directly computed at the nodes by T

ε 0    δ  ε1  I =1 ΩI  2  ε  NP

∑∫

 Ab  b B E 

ε 0    δ  ε1  I =1 ΩI  2 ε 

E  ε 0     Db F1  ε1 dΩ I + 2 F1 H    ε  Bb

T

NP

∑∫

u 0  NP   δ  u1  I =1 Ω  2 I u 

T

∑∫

T

ε 0  NP   δ  ε1  I =1 Ω  2 I ε 

∑∫

 Ab  b B E   I1 I  2  I 4  Ab  b B E  T

 w0 I , x  h δ  I =1 Ω I   w0 I , y  NP

∑∫

∑∫

E  ε 0     Db F1  ε1 dΩ I + F1 H  ε 2   

Bb

I2 I3 I5

T

s0 ε  δ  s1  I =1 Ω I  ε  NP

As  s  B

B s  ε s 0    dΩ I = Ds   ε s1 

T

s0 ε  δ  s1   ε  I =1 Ω I NP

∑∫

As  s  B

NP

∑ ∫ δw

0 I q0 dΩ I

B s  ε s 0    dΩ I = ... Ds   ε s1  (25)

&&0  I 4  u  1  I 5   && u dΩ I  2 I 6  u && 

E  ε 0     Db F1  ε1 dΩ I + F1 H  ε 2    Bb

(24)

I =1 ΩI

T

s0 s ε   A δ  s1   s I =1 Ω  ε   B I NP

∑∫

Bs  ε s 0    dΩ I + ... Ds  ε s1 

0  σ x0 σ xy   w0 I , x   0  dΩ I = 0 0  σ xy σ y  w0 I , y 

MKI node

Voronoi domain

Figure 2. Geometry of a representative node domain and its Voronoi domain.

10

(26)

3. Naturally stabilized nodal integration Particularly, the direct nodal integration leads to unstable solutions due to rank deficiency of approximated matrices. To overcome this shortcoming, Hillman and Chen [33] proposed a naturally stabilized nodal integration (NSNI) technique which is a naturally implicit gradient expansion of the displacements field in Eq. (15) as follows NP

NP

I

I

u h ( x ) ≈ ∑ φI ( x )q I + ∑ ∇φI ( x ). ( x − x ) q I

(27)

where ‘.’ defines an inner product and ∇φI ( x ) is the derivation of shape function φI ( x ) . To reduce the order of high-order derivations in the approximation formulation, a new shape function φˆI is constructed such that the following condition is satisfied NP

NP

∑ φˆ ( x )q = ∑ ∇φ ( x )q I

I

I

I

I

= uˆ h ( x )

(28)

I

The condition in Eq. (28) is easily obtained in the MKI by m

n

m

n

j =1

k =1

j =1

k =1

φˆIx ( x) = ∑ p j , x ( x)AjI + ∑ rk , x ( x)BkI ; φˆIy (x) = ∑ p j , y ( x) AjI + ∑ rk , y (x)BkI

(29)

The implicit gradient expansion by the first order Taylor series expansions brings naturally to a stabilized integration technique when evaluating the displacements and strains at each node. This expression is similar to the finite element context presented by Liu [46], but the present displacement and strain are directly evaluated at nodes in the Cartesian coordinates instead of in the isoparametric ones. Substituting Eq. (27) into Eq. (6), bending and shear strains with nodal integration zones at x = x I are expressed by ε 0I ( u h ( x ) ) ≈ εI0 ( u h ( x I ) ) + ( x − xI ) ε 0xI ( uˆ hx ( x I ) ) + ( y − yI ) ε 0yI ( uˆ hy ( x I ) ) ; ε1I ( u h ( x ) ) ≈ εI1 ( u h ( x I ) ) + ( x − xI ) ε1xI ( uˆ hx ( x I ) ) + ( y − y I ) ε1yI ( uˆ hy ( x I ) ) ; ε 2I ( u h ( x ) ) ≈ εI2 ( u h ( x I ) ) + ( x − xI ) ε 2xI ( uˆ hx ( x I ) ) + ( y − yI ) ε 2yI ( uˆ hy ( x I ) ) ;

(30)

ε sI 0 ( u h ( x ) ) ≈ εIs 0 ( u h ( x I ) ) + ( x − xI ) ε sxI0 ( uˆ hx ( x I ) ) + ( y − yI ) ε syI0 ( uˆ hy ( x I ) ) ; ε sI1 ( u h ( x ) ) ≈ εIs1 ( u h ( x I ) ) + ( x − xI ) ε sxI1 ( uˆ hx ( x I ) ) + ( y − y I ) ε syI1 ( uˆ hy ( x I ) )

where NP

NP

NP

I

I

I

u h ( x I ) = ∑ φI ( x )q I ; uˆ hx ( x I ) = ∑ φˆIx ( x I )q I ; uˆ hy ( x I ) = ∑ φˆIy ( x I )q I

A compact form of the bending and shear strains can be rewritten as

11

(31)

ε 0I ≈ B 0I q I + ( x − xI ) B 0Ix q I + ( y − y I ) B 0Iy q I ; ε1I ≈ B1I q I + ( x − xI ) B1Ix q I + ( y − y I ) B1Iy q I ; ε 2I ≈ B 2I q I + ( x − xI ) B 2Ix q I + ( y − yI ) B 2Iy q I ; ε sI 0 ≈ B sI 0 q I + ( x − xI ) B sIx0 q I + ( y − yI ) B sIy0 q I ;

(32)

ε sI1 ≈ B Is1q I + ( x − xI ) B sIx1q I + ( y − y I ) B sIy1q I

where φI , x  0 B I = 0 φI , y  φˆIy , x  B 0Iy =  0 ˆ φIy , y

0

φI , y φI , x 0 φˆ

Iy , y

φˆ

Iy , x

0 0 0 0 0  0 0 0 0 0 ; 0 0 0 0 0 

φˆIx , x  B 0Ix =  0 ˆ φIx , y

0 0 0 0 0  0 0 0 0 0 ;  0 0 0 0 0

 0 0 0 0 0 φI , x 0    φI , y  ; B = − 0 0 0 0 0 0  0 0 0 0 0 φI , y φ I , x   

 0 0 0 0 0 φˆIx , x  B1Ix = −  0 0 0 0 0 0  ˆ  0 0 0 0 0 φIx , y 0 0 0 φI , x  2 B I = − 0 0 0 0 0 0 0 φI , y   0 0 0 φˆIy , x  B 2Iy =  0 0 0 0  ˆ  0 0 0 φIy , y

0 0 φˆIx, x B sIx0 =  ˆ 0 0 φIx, y

0

φI , y φI , x 0 φˆ

Iy , y

φˆIy , x

0 0 −φˆIx 0 0 0

 0 0 0 φI , x 0 B sI1 =  φI , y 0 0 0 0

0 φˆ

Ix , y

φˆ

Ix , x

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

1 I

0 φˆ

  Ix , y  ;  φˆIx , x 

 0 0 0 0 0 φˆIy , x  B1Iy = −  0 0 0 0 0 0  ˆ  0 0 0 0 0 φIy , y

0 0  0 0 ; 0 0

 0 0 0 φˆIx , x 0  B 2Ix =  0 0 0 0 φˆIx, y  ˆ ˆ  0 0 0 φIx , y φIx, x

0 0  0 0 ;  0 0

 ; −φˆIx 

0

0 0 ; 0 0 

0 0 φI , x B sI 0 =  0 0 φI , y

0 0 0 φˆIx , x B sIx1 =  0 0 0 0

0 0 0 φˆIy , x B sIy1 =  0 0 0 0

0 φˆ

Iy , y

0 0  0 0 ;  0 0 0  ; −φI 

0 0 −φI 0 0 0

0 0 φˆIy , x B sIy0 =  ˆ 0 0 φIy , y

0 0 −φˆIy 0 0 0

0 φˆ

Ix , y

  Iy , y  ;  φˆIy , x 

0 φˆ

 ; −φˆIy  0

0 0 ; 0 0 

0 0 ; 0 0 

Inserting Eq. (27) into Eq. (1), the displacement fields u 0 , u1 and u 2 can be written by u 0I = N 0I q I + ( x − xI ) N 0Ixq I + ( y − y I ) N 0Iy q I ;

12

(33)

u1I = N1I q I + ( x − xI ) N1Ix q I + ( y − y I ) N1Iy q I ;

(34)

u 2I = N 2I q I + ( x − xI ) N 2Ix q I + ( y − y I ) N 2Iy q I

where φI 0 N I =  0  0

0

0

φI

0

0 0 0 0 0 0 0 0  ; 0 0 0 0 

φI

φˆIy 0  N 0Iy =  0 φˆIy   0 0

0 0 φˆ

Iy

0 0 0 0  0 0 0 0 ;  0 0 0 0

 0 0 0 0 0 φˆIx  N1Ix = −  0 0 0 0 0 0   0 0 0 0 0 0 0 0 0 φ I 2 N I = 0 0 0 0 0 0 0 0

0

φI 0

0   φˆIx  ;  0  

0 0 0 0  ; 0 0 

φˆIx  N 0Ix =  0   0

0 φˆ

Ix

0

0 0 0 0  0 0 0 0 ;  0 0 0 0 

0 0 φˆ

Ix

 0 0 0 0 0 φI N = −  0 0 0 0 0 0  0 0 0 0 0 0 1 I

(35)

0 φI  ; 0 

0 0 0 0 0 φˆIy 0    N1Iy = − 0 0 0 0 0 0 φˆIy  ;   0 0 0 0 0 0 0   0 0 0 φˆIx 0  N 2Ix =  0 0 0 0 φˆIx   0 0 0 0 0

 0 0 0 φˆIy 0  N 2Iy =  0 0 0 0 φˆIy   0 0 0 0 0

0 0  0 0 ;  0 0 

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

The derivation of the transverse displacements is also described by  w0 I , x  g g g   = B I q I + ( x − xI ) B Ixq I + ( y − y I ) B Iy q I  w0 I , y 

(36)

where  0 0 φI,x B Ig =   0 0 φI,y

0 0 0 0 ; 0 0 0 0 

0 0 φˆIy,x B Iyg =  0 0 φˆIy,y

0 0 φˆIx,x B Ixg =  0 0 φˆIx,y 0 0 0 0  0 0 0 0 

13

0 0 0 0 ; 0 0 0 0 

(37)

Substituting the expansions in Eqs. (32), (34) and (36) into Eqs. (24), (25) and (26), respectively, and ignoring the zero terms (

∫ ( x − x ) dΩ = ∫ ( y − y ) dΩ I

I

ΩI

I

= 0 ), the weak forms for static, free

I

ΩI

vibration and buckling analyses with naturally stabilized nodal integration technique are rewritten in the compact form of equations:

Kq = F

(38)

(K - ω M )q = 0

(39)

(K - λ

(40)

2

cr

Kg )q = 0

where K , M , K g and F are the global stiffness matrix, mass matrix, geometry stiffness matrix and force vector, respectively, in which K = K 0 + K x + K y ; M = M 0 + M x + M y and K g = K 0g + K gx + K gy

(41)

and  B 0 T  A b I NP    1   b  K0 = B   B  I I =1  2    B I   E 

E   B 0I  T   1  B sI 0  b D F1  B I  +  s1  B F1 H  B 2I   I   

A s  s  B

 s0     B B I    WI ; D s   B sI1    

 B 0 T  A b Ix NP    1   b  Kx = B   B  Ix I =1  2    B Ix   E 

E  B 0Ix  T   1  B sIx0  b D F1 B Ix  +  s1  B F1 H  B 2Ix   Ix   

A  s  B

 s0    B  B Ix    M Ix ; D s   B sIx1    

 B 0 T  A b   Iy     K y =  B1Iy   B b  I =1  2    B Iy   E 

E  B 0Iy  T   1   B sIy0  b D F1 B Iy  +  s1  B F1 H  B 2Iy   Iy   

A s  s  B

 s0 B s  B Iy     s1  M Iy ; D s   B Iy    

Bb

∑

∑

NP

∑

Bb

Bb

NP

F=

∑ q {0 0

0 w0 I

s

s

s

T

0 0 0 0} WI

I =1 T

 N 0I  NP   1  M0 = N I  I =1  2   N I 

∑

 I1 I  2  I 4

I2

I3 I5

0 I 4  N I    I 5   N1I WI ; I 6   N 2I   

T

 N0Ix  NP   1  Mx = N Ix  I =1  2   N Ix 

∑

14

 I1 I  2  I 4

I2

I3 I5

0 I 4   N Ix    I 5   N1Ix M Ix ; I 6   N 2Ix   

(42)

T

 N 0Iy  NP   1  My =  N Iy  I =1  2   N Iy 

∑ NP

K gx = h

∑(B ) g Ix

I =1

in which, WI =

∫

T

 I1 I  2  I 4

 σ x0  0 σ xy

0 I 4   N Iy    I 5   N1Iy M Iy ; I 6   N 2Iy   

I2

I3 I5

NP

K 0g

NP

 B gIx M Ix ; σ y0 

ΩI

∑( )

2

∫ ( x − xI ) dΩ I ,

T

I =1

σ xy0 

dΩ I , M Ix =

=h

B gI

K gy = h

∑(B ) g Iy

T

I =1

M Iy =

ΩI

2

∫ ( y − yI ) dΩ I

 σ x0  0 σ xy

σ xy0 

 σ x0  0 σ xy

σ xy0 

 B gI WI ;

σ y0 

 B Iyg M Iy ;

σ 0y 

in Eq. (42) are the area and the

ΩI

second moments of inertia of each nodal domain, respectively. In addition, ω and λcr are the natural frequency and the critical buckling load, respectively. 4. Geometries and imposing essential boundary conditions A Matlab mesh generator program with unstructured three-node triangular (T3) mesh is available by Koko [47]. The 2D geometries and T3 meshes are easy generated by this program. To apply the NSNI technique, a Voronoi mesh is used as shown in Figure 2. The area and second moments of inertia of each nodal domain Ω I are automatically calculated through the vertices of this representative nodal domain. The essential boundary conditions are imposed easily and directly as in FEM because the MKI shape function is satisfied the Kronecker delta function property. Two types of following Dirichlet boundary conditions are used as follows
Simply supported (rectangular plate): u0 = w0 = θ x = β x = 0 at y = 0, b and v0 = w0 = θ y = β y = 0 at x = 0, a




Fully clamped: u0 = v0 = w0 = θ x = θ y = β x = β y =0 at boundary

5. Numerical results and discussions In this section, we illustrate several examples to show the performance of the present method. Computations are performed on a laptop of CPU Core i7-6500U CPU@ 2.66GHz with RAM 8.0 GB. Various kind of elastic properties of a lamina are taken as follows: •

Material I: E1 = 25E2 , G12 = G13 = 0.5E2 , G23 = 0.2 E2 ,ν 12 = 0.25, ρ = 1.

•

Material II [48]: E1 = 40 E2 , G12 = G13 = 0.6 E2 , G23 = 0.5E2 ,ν 12 = 0.25, ρ = 1.

15

•

Material III [49]: E1 = 2.45 E2 , G12 = G13 = 0.48 E2 , G23 = 0.2 E2 ,ν 12 = 0.23, ρ = 1.

5.1 Static analysis 5.1.1 Four layer [0/90/90/0] laminated square plate under sinusoidally distributed loading

A four-layer cross-ply [0/90/90/0] simply supported square plate, which material I is used, under a  π x  π y  sinusoidally distributed loading q0 = q0 sin    is exampled. The length-to-thickness ratio is  a  b 

taken equal to 4, 10, 20 and 100, respectively. The normalized displacement, in-plane and shear stresses are given by w=

h2 a b h h2 a b h 100 E2 h3 a b ( , ,0); ( , , ); w σ = σ σ = σ y ( , , ); x x y 2 2 2 2 2 2 2 4 q0 a 4 q0b 2 q0b 2

(43) h2 h h b τ xy = τ xy (0,0, ); τ xz = τ xz (0, ,0) q0 b2 2 q0 b 2 To show the accuracy of the present solution, the convergence of the normalized displacement is studied. The length-to-thickness ratio is taken equal to 4 (a/h=4). The plate is modeled by 132, 203, 384 and 593 nodes, as shown in Figure 3, respectively. Figure 4a depicts the percentage error (%) of the normalized displacement of the present solution based on TSDT and ITSDT. The obtained solutions are compared with an exact solution based on TSDT proposed by Reddy [17] and an exact 3D solution proposed by Pagano [1]. The current solution based on TSDT converges to the exactTSDT solution while the current solution based on ITSDT to the exact 3D solution. Besides, the convergence of the computational times is drawn in Figure 4b. It is clearly observed that the present method can produce very accurate results with low computational cost. For comparison results, various length-to-thickness ratios are investigated. Table 2 gives the comparison of the non-dimensional transverse deflection and stresses of the four-layer simply supported square plate. The obtained results are compared with those reported by Reddy [17], Akhras et al. [50], Ferreira et al. [51], Ferreira [52], Ferreira et al. [53], Roque et al. [54], Wang and Shi [55]

and Pagano [1]. It is clearly observed that a good agreement is found for all transverse deflection and stresses corresponding with various length-to-thickness ratios. The obtained results from the presentITSDT solution are more accurate than the present-TSDT ones for all transverse deflection and stresses of thick and moderately thick plates but the differences between two solutions are not significant for the thin plate case when compared to the exact 3D elasticity solution [1] and other solutions. The distribution of stresses through the plate thickness with a/h = 4 and 10 based on TSDT are plotted in Figure 5, respectively.

16

Table 2: The normalized displacement and stresses of the four-layer [0/90/90/0] laminated square plate under a sinusoidally distributed load. a/h

Method

w

σx

σy

σ xz

σ xy

4

Reddy [17] Akhras et al. [50] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Roque et al. [54] Wang and Shi [55] Pagano [1] Present-TSDT Present-ITSDT Reddy [17] Akhras et al. [50] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Roque et al. [54] Wang and Shi [55] Pagano [1] Present-TSDT Present-ITSDT Reddy [17] Akhras et al. [50] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Roque et al. [54] Wang and Shi [55] Pagano [1] Present-TSDT Present-ITSDT Reddy [17] Akhras et al. [50] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Roque et al. [54] Wang and Shi [55] Pagano [1] Present-TSDT Present-ITSDT

1.8939 1.8937 1.8864 1.9075 1.9091 1.8842 1.9073 1.954 1.8979 1.9221 0.7149 0.7147 0.7153 0.7309 0.7303 0.735 0.7368 0.743 0.7178 0.7254 0.5061 0.506 0.507 0.5121 0.5113 0.5127 0.5138 0.517 0.5054 0.5076 0.4343 0.4343 0.4365 0.4374 0.4348 0.4345 0.4355 0.4347 0.4308 0.4309

0.6806 0.6651 0.6659 0.6432 0.6429 0.756 0.7361 0.72 0.7186 0.7160 0.5589 0.5456 0.5466 0.5496 0.5487 0.5637 0.5609 0.559 0.5726 0.5727 0.5523 0.5393 0.5405 0.5417 0.5407 0.544 0.5433 0.543 0.5549 0.5547 0.5507 0.5387 0.5413 0.542 0.5391 0.5388 0.5387 0.539 0.5458 0.5457

0.6463 0.6322 0.6313 0.6228 0.6265 0.6777 0.6994 0.666 0.6637 0.6715 0.3974 0.3888 0.4383 0.3956 0.3966 0.4055 0.4077 0.403 0.4121 0.4156 0.311 0.3043 0.3648 0.3056 0.3073 0.3094 0.3098 0.309 0.3240 0.3252 0.2769 0.2708 0.3359 0.2697 0.2711 0.271 0.271 0.271 0.2919 0.2920

0.2109 0.2064 0.1352 0.2166 0.2173 0.1885 0.211 0.27 0.2265 0.2494 0.2697 0.264 0.3347 0.2888 0.2993 0.2908 0.3002 0.301 0.3054 0.3428 0.2883 0.2825 0.3818 0.3248 0.3256 0.3203 0.3279 0.328 0.3449 0.3888 0.2948 0.2897 0.4106 0.3232 0.3359 0.3354 0.3389 0.339 0.4272 0.4356

0.045 0.044 0.0433 0.0441 0.0443 0.043 0.0435 0.0467 0.0425 0.0433 0.0273 0.0268 0.0267 0.0273 0.0273 0.0272 0.0274 0.0276 0.0253 0.0255 0.0233 0.0228 0.0228 0.023 0.023 0.0223 0.0231 0.023 0.0219 0.0219 0.0217 0.0213 0.0215 0.0216 0.0214 0.0213 0.0214 0.0214 0.0212 0.0212

10

20

100

17

a) 132 nodes

b) 203 nodes

c) 384 nodes

d) 593 nodes

Figure 3. Distributed nodes of the square plate.

100

7 3D Exact−TSDT Present−TSDT Present−ITSDT

5

80 70

4 3 2 1

60 50 40 30

0

20

−1 −2

Present−TSDT Present−ITSDT

90

Times (second)

Error % of normalized displacement

6

10 150

200

250

300 350 400 450 Number of nodes

500

550

600

18

0

150

200

250

300 350 400 450 Number of nodes

500

550

600

a)

b)

0.5

0.4

0.4

0.3

0.3

0.2 0.1

present−TSDT, a/h=10 present−TSDT, a/h=4 present−ITSDT, a/h=10 present−ITSDT, a/h=4

0 −0.1 −0.2

0.1 0 −0.1 −0.2 −0.3

−0.4

−0.4 −0.6

−0.4

−0.2 0 0.2 Stress, σ x (a/2,b/2,z)

0.4

0.6

−0.5 −0.8

0.8

0.5

0.5

0.4

0.4

0.3

0.3

0.2 0.1

present−TSDT, a/h=10 present−TSDT, a/h=4 present−ITSDT, a/h=10 present−ITSDT, a/h=4

0 −0.1 −0.2

0.1

−0.1

−0.4 0.2

0.3 0.4 Stress, τ xz(0,b/2,z)

0.5

0.6

0.7

19

−0.2 0 0.2 Stress, σ y (a/2,b/2,z)

0.4

0.6

0.8

−0.2

−0.4 0.1

−0.4

present−TSDT, a/h=10 present−TSDT, a/h=4 present−ITSDT, a/h=10 present−ITSDT, a/h=4

0

−0.3

0

−0.6

0.2

−0.3

−0.5

present−TSDT, a/h=10 present−TSDT, a/h=4 present−ITSDT, a/h=10 present−ITSDT, a/h=4

0.2

−0.3

−0.5 −0.8

Nomalized thickness z/h

Nomalized thickness z/h

0.5

Nomalized thickness z/h

Nomalized thickness z/h

Figure 4. Convergence of errors (%) of normalized displacement and computational times.

−0.5

0

0.05

0.1 0.15 0.2 Stress, τ yz(a/2,0,z)

0.25

0.3

0.5 present−TSDT, a/h=10 present−TSDT, a/h=4 present−ITSDT, a/h=10 present−ITSDT, a/h=4

0.4

Nomalized thickness z/h

0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.05

0 Stress, τ xy(0,0,z)

0.05

Figure 5. The distribution of stresses across the thickness of the plate with a/h = 4 and 10. 5.1.2 Three-layer sandwich square plate subjected to a uniform loading

We consider a sandwich square plate with simply supported boundary subjected to a uniform transverse loading q0. The sandwich plate is made of three layers including two face sheets with the same thickness and one core. The side-to-thickness ratio is taken equal to 10 (a/h = 10). The thickness of plate to thickness of face sheet ratio is chosen equal to 10 (h/ hface = 10). A distribution of 593 nodes as in Figure 3d is used. The material properties of core and face layers, respectively, are given by

0 0 0  0.999781 0.231192   0.231192 0.524886  0 0 0    and Q face = RQcore Qcore =  0 0 0.262931 0 0   0 0 0 0.266810 0    0 0 0 0 0.159914  where R is a scale factor.

The normalized displacement and stresses of the plate are defined by w=

0.999781 a b a b h a b 2h a b 2h w( , ,0); σ 1x = σ 1x ( , , ) / qo ;σ x2 = σ 1x ( , , ) / qo ;σ x3 = σ x2 ( , , ) / qo hqo 2 2 2 2 2 2 2 5 2 2 5

The exact solution of this problem is found in Srinivas [56]. Three values of the scale factor as R = 5, 10, 15 are studied, respectively. The normalized displacement and stresses obtained from the present solution are given in Table 3. Numerical results are compared with those obtained by Srinivas [56], Pandya and Kant [57], Ferreira et al. [51], Ferreira [52], Ferreira et al. [53], Roque et al. [54],

20

Mantari et al. [58] and Grover et al. [59]. It is found that the obtained results from the current solution are in good agreement with the reference solutions for both three values of R. It can be seen that the normalized displacement decreases when increasing the scale factor R. Table 3: The normalized displacement and stresses of the square sandwich plate under uniform load. R

Method/Authors

w

σ 1x

σ x2

σ x3

5

Pandya and Kant [57] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Grover et al. [59] Mantari et al. [58] Roque et al. [54] Exact [56] Present-TSDT Present-ITSDT Pandya and Kant [57] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Grover et al. [59] Mantari et al. [58] Roque et al. [54] Exact [56] Present-TSDT Present-ITSDT Pandya and Kant [57] Ferreira et al. [51] Ferreira [52] Ferreira et al. [53] Grover et al. [59] Mantari et al. [58] Roque et al. [54] Exact [56] Present-TSDT Present-ITSDT

256.13 257.11 257.523 258.179 255.644 256.706 259.12 258.97 255.7734 254.7220 152.33 154.658 158.38 158.912 154.55 155.498 159.5 159.38 154.9867 154.1756 110.43 114.644 120.988 121.347 115.82 115.919 121.88 121.72 116.2468 115.8825

62.38 60.366 59.968 60.063 60.675 60.525 60.338 60.353 60.7301 60.6698 64.65 65.381 64.846 64.993 65.741 65.542 65.279 65.332 65.9874 65.9192 66.62 66.919 66.291 66.436 67.272 67.185 66.73 66.787 67.7017 67.6110

46.91 47.003 46.291 46.393 47.055 47.061 46.57 46.623 46.5768 46.7230 51.31 49.973 48.443 48.601 49.798 49.708 48.279 48.857 49.0382 49.2034 51.97 50.323 47.899 48.011 49.813 49.769 48.204 48.299 48.8581 48.9569

9.382 9.401 9.258 9.279 9.411 9.412 9.314 9.34 9.3154 9.3446 5.131 4.997 4.844 4.86 4.979 4.971 4.8766 4.903 4.9038 4.9203 3.465 3.355 3.193 3.201 3.321 3.318 3.2136 3.238 3.2572 3.2638

10

15

5.2 Free vibration analysis 5.2.1. Square plates

A four-layer cross-ply [0/90/90/0] square plate with simply supported boundary conditions is considered.

The

material

II

is

used.

The

non-dimensional

frequency

is

defined

by

ϖ = (ω a 2 / h ) ( ρ / E2 ) , where ρ and E2 are a density mass and an elastic modulus of the material in 1/ 2

second direction, respectively. The square plate is modeled by 588 nodes to calculate natural

21

frequency, as shown in Figure 6. The non-dimensional first frequency of the four-layer plate with various length-to-thickness a/h and elastic modulus E1/E2 ratios are computed in Table 4 and Table 5 corresponding with fixed a/h and E1/E2, respectively. The numerical results are compared with those obtained in Kdheir [60], Reddy [41], Liew et al. [48], Ferreira et al. [15], Ferreira [61], Zhen and Wanji [62], Whu and Chen [63], Matsunaga [64], Cho et al. [65] and Thai et al. [66]. From Table 4 and Table 5, a good agreement is found for two cases of various length-to-thickness and elastic modulus ratios, respectively. Table 4: A first non-dimensional frequency of a four-layer simply supported square plate (a/h=5).

E1 /E2

Method/Authors

10 8.2526 8.2794 8.2924 8.2982 8.2982 8.2792 8.2857 8.2929

Ferreira [61] Ferreira et al. [15] Liew et al. [48] Kdheir [60] Reddy [41] Thai et al. [66] Present-TSDT Present-ITSDT

20 9.4974 9.5375 9.5613 9.5671 9.5671 9.5454 9.5291 9.5354

30 10.2308 10.2889 10.3200 10.3260 10.3260 10.3028 10.2599 10.2619

40 10.7329 10.8117 10.8490 10.8540 10.8540 10.8298 10.7594 10.7553

Table 5: A first non-dimensional frequency of a four-layer [0/90/90/0] simply supported square plate (E1/E2 = 40). Method/authors

Zhen and Wanji [62] Whu and Chen [63] Matsunaga [64] Cho et al. [65] Thai et al. [66] Present-TSDT Present-ITSDT

a/h

4 9.2406 9.193 9.1988 9.3757 9.2786 9.2661

5 10.7294 10.682 10.6876 10.673 10.8298 10.7594 10.7553

10 15.1658 15.069 15.0721 15.066 15.1249 15.1275 15.1394

20 17.8035 17.636 17.6369 17.535 17.6521 17.7008 17.7080

25 18.2404 18.055 18.0557 18.054 18.0657 18.1235 18.1287

50 18.9022 18.67 18.6702 18.67 18.6728 18.7487 18.7502

100 19.1566 18.835 18.8352 18.835 18.8359 18.9237 18.9241

A further study of a square plate with three-layer cross-ply [0/90/0/] under fully clamped boundary conditions is examples. Material II is also used. Non-dimensional frequency is defined by

ϖ = (ω b 2 / π 2 ) ( ρ h / D0 ) , in which D0 = E2 h3 / 12 (1 −ν 12ν 21 ) . The first six non-dimensional 1/ 2

frequencies of a three-layer square plate with different values of side-to-thickness ratios are given in Table 6. The obtained results from the present solution are compared with those given by Liew et al. [48], Zhen and Wanji [62], Ferreira et al. [15] and Thai et al. [66]. It can be seen that the present

22

solution is very good for all side-to-thickness ratios when compared to those solutions. The first six mode shapes of a three-layer fully clamped square plate (b/h=10) are plotted in Figure 7. Table 6: A comparison of non-dimensional frequencies of a three-layer clamped square plate. b/h

Authors

5

Liew et al. [48] Zhen and Wanji [62] Ferreira et al. [15] Present-TSDT Present-ITSDT Liew et al. [48] Zhen and Wanji [62] Ferreira et al. [15] Thai et al. [66] Present-TSDT Present-ITSDT Liew et al. [48] Zhen and Wanji [62] Ferreira et al. [15] Thai et al. [66] Present-TSDT Present-ITSDT Liew et al. [48] Zhen and Wanji [62] Ferreira et al. [15] Thai et al. [66] Present-TSDT Present-ITSDT

10

20

100

1 4.447 4.450 4.4466 4.5243 4.5198 7.411 7.484 7.4106 7.7082 7.4814 7.4853 10.953 11.003 10.9528 11.1029 10.9623 10.9670 14.666 14.601 14.4455 14.4386 14.5924 14.5927

2 6.642 6.524 6.6419 6.5987 6.55487 10.393 10.207 10.3944 10.5616 10.4334 10.3803 14.028 14.064 14.036 14.1296 14.2146 14.1922 17.614 17.812 17.5426 17.3929 17.8627 17.8616

Modes 3 4 7.700 9.185 8.178 9.473 7.6996 9.1852 8.0628 9.4532 8.07966 9.4408 13.913 15.429 14.34 14.863 13.9128 15.4403 15.0084 15.4540 14.3144 15.4198 14.3511 15.2710 20.388 23.196 20.321 23.498 20.4533 23.1974 20.3947 23.9449 20.8994 23.5333 20.8008 23.5640 24.511 35.532 25.236 37.168 25.1868 37.8851 24.2701 37.8223 25.5546 38.0299 25.5486 38.0121

Figure 6. A distributed node of a square plate.

23

5 9.738 9.492 9.7379 9.6114 9.5350 15.806 16.07 15.8061 16.7515 16.2078 16.2054 24.978 25.35 24.9827 25.6783 25.5061 25.5200 39.157 38.528 39.5489 39.4330 38.6212 38.6232

6 11.399 11.769 11.3992 11.8192 11.7739 19.572 19.508 19.5797 20.2844 19.9455 19.8598 29.237 29.118 29.2795 29.2604 30.1075 30.0596 40.768 40.668 39.6519 43.4645 40.6383 40.6396

Mode 1, Omega = 7.4814

Mode 2, Omega = 10.4334

Mode 3, Omega = 14.3144

1

1

1

0.5

0.5

0.5

0

0

0

−0.5

−0.5

−0.5

−1 1

−1 1

0.8

−1 1 0.8

0 0.2

0.6

0.2

0.6

0.4

0.4

0.4

0.4

0.4

0.6 0.2

0.8 0

0 0.2

0.6

0.4

0.6 0.2

0.8

0

0

1

Mode 4, Omega = 15.4198

0.6 0.2

0.8

Mode 5, Omega = 16.2078

1

1

0.5

0.5

0.8 0

1

1

Mode 6, Omega = 19.9455

1 0.5 0

0

0

−0.5

−0.5

−1 1

−1 1

−0.5

0.8

0 0.2

0.6 0.4

0.4

−1 1 0.8

0.8

0 0.2

0.6

0.8 0

1

0.4

0.4 0.6

0.4

0.4

0.2

0.8

0.6

0.6 0.2

0 0.2

0.6

0.2

0.8 0

0

1

1

Figure 7. The first six modes shape of a three layer [0/90/0] clamped square plate with a/h= 10 using the present-TSDT. 5.2.2. Elliptical plates A three-layer [0/90/0] fully clamped elliptical plate of two radii being equal a = 5 and b = 2.5 is studied, as shown in Figure 8a. A set of nodes (566) of the elliptical plate is drawn in Figure 8b. Material III is used. The non-dimensional frequency is given by ϖ = (ω a 2 ) ( ρ h / D0 ) , where 1/ 2

D0 = E1h3 / 12 (1 − ν 12ν 21 ) .

In order to validate the accuracy of the present solution for the elliptical plate, a comparison has been carried out with previously published results by Chen et al. [49] for the thin plate and Thai et al. [67, 66] for both thin and thick plates based on the layer-wise plate theories, as reported in Table 7. It is observed that a good agreement with those results is found for both thin and thick plates. For moderately thick and thick plates, the present result is slightly smaller than the one from [67] and [66]. For the thin plate, the present result is slightly larger than those published results. Figure 9 shows the first six mode shapes of three-layer fully clamped laminated elliptical plate (a/h=10).

24

b

a

a)

b)

Figure 8. Geometry and distributed nodes of an ellipse plate. Table 7: A non-dimensional frequencies of the three-layer fully clamped laminated ellipse plate. a/h

Method

5

Layer-wise [67] Layer-wise [66] Present-TSDT Present-ITSDT Layer-wise [67] Layer-wise [66] Present-TSDT Present-ITSDT Layer-wise [67] Layer-wise [66] Present-TSDT Present-ITSDT Layer-wise [67] Layer-wise [66] CLPT [49] Present-TSDT Present-ITSDT

10

20

100

1 14.157 14.1353 13.9944 13.9972 17.184 17.188 16.9830 16.9851 18.329 18.3666 18.1466 18.1472 18.755 18.8085 18.8100 18.8608 18.8607

2 19.969 20.0216 19.6383 19.5740 25.714 25.7979 25.1320 25.0896 28.280 28.4097 27.6288 27.6130 29.332 29.4663 29.5800 29.8159 29.8147

25

Modes 3 4 27.114 28.855 27.2208 28.8527 26.7469 28.4135 26.5900 28.4553 36.982 39.196 37.0987 39.0942 36.2588 38.6334 36.1245 38.6652 42.255 44.321 42.4237 44.2872 41.3751 43.7577 41.3184 43.7696 44.792 46.508 44.8105 46.5297 44.9900 46.7200 46.5424 46.7256 46.5397 46.7250

5 34.943 34.9609 33.8107 33.7833 49.148 49.1092 47.1906 47.1658 57.090 57.1596 54.5920 54.5799 60.792 60.9055 61.3400 59.9272 59.9264

6 35.062 35.2453 34.6643 34.4096 50.259 50.3576 49.5354 49.2638 59.827 59.9436 58.9077 58.7758 65.6230 64.763 65.1400 69.2497 69.2409

Mode 1, Omega = 16.983

Mode 2, Omega = 25.132

Mode 3, Omega = 36.2588

3 2

2

1

1

0

0

−1

−1

−2

−2

5

2 1 0 −1 −2 5

5

0

0

0 −2 0 −5

−2

−2 0

0

2

−5

Mode 4, Omega = 38.6334

2

−5

Mode 5, Omega = 47.1906

Mode 6, Omega = 49.5354

2

2

2

1

1

1

0

0

0

−1

−1

−1

−2

−2

−2

5

5

5

0

0

0

−2

−2

−2

0

0 −5

2

2

−5

2

0 −5

2

Figure 9. The first six modes shape of the three-layer fully clamped ellipse plate with a/h= 10 (using the present-TSDT). 5.2.3. Square plate with a complicated hole

A complex geometry is studied to calculate the natural frequency of the laminated plates. Let us a simply supported square plate with a complicated hole, as shown in Figure 10a. A node distribution (693 nodes) of the plate is plotted in Figure 10b. In this example, the three-layer plate with two different kinds of angle plies as [0/90/0] and [45/-45/45] are presented. The material properties and non-dimensional frequency parameter are taken similar to the previous example (subsection 5.2.2). Numerical solutions for thin plates of this example are available. In addition, obtained results from the present solution for both moderately thick and thick plates are also studied. The first six nondimensional frequencies are listed in Table 8 together with the available results by Yin et al. [69], obtained by using isogeometric analysis (IGA) based on CLPT, by Shojaee et al. [68], obtained by using IGA-CLPT and by Bui et al. [70], obtained by using a meshfree (MF) method based on CLPT.

26

For comparison, a good agreement between the present solution and those solutions is found. The first six mode shapes of the three-layer simply supported square plate with a complicated hole are drawn in Figure 11. Table 8: A non-dimensional frequency parameter ϖ of the three simply supported laminated square plate a complicated hole. Angle plies

a/h

Method

(0/90/0)

166.67

IGA-CPLT [68] IGA-CLPT [69] MF-CLPT [70] Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT IGA-CPLT [68] IGA-CLPT [69] MF-CLPT [70] Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT

100 20 10 5 (45/-45/45)

166.67

100 20 10 5

Modes 1 18.201 18.190 18.278 18.260 18.260 18.190 18.190 17.739 17.732 16.900 16.884 14.764 14.730 21.313 20.982 21.736 21.403 21.403 21.319 21.318 20.744 20.736 19.668 19.645 14.362 14.362

27

2 31.082 31.087 32.264 31.687 31.685 31.360 31.358 29.350 29.327 26.380 26.332 15.565 15.565 34.801 34.848 36.079 35.388 35.386 35.075 35.072 32.810 32.789 28.725 28.725 16.918 16.874

3 36.096 35.655 36.134 36.892 36.889 36.452 36.447 33.891 33.851 30.184 30.112 20.807 20.740 38.289 37.559 39.975 39.098 39.095 38.592 38.586 35.804 35.758 29.433 29.392 23.129 23.089

4 56.473 55.452 57.151 58.160 58.157 57.484 57.477 53.326 53.251 31.129 31.129 23.361 23.279 60.897 59.325 63.897 62.708 62.703 61.972 61.964 57.231 57.163 31.808 31.730 24.395 24.302

5 62.523 62.582 65.853 64.642 64.638 63.623 63.615 57.775 57.704 46.572 46.415 24.649 24.649 66.885 67.518 68.525 68.716 68.713 67.715 67.707 57.449 57.449 49.685 49.556 27.423 27.423

6 83.66 82.383 90.678 86.649 86.647 85.413 85.408 62.258 62.258 49.259 49.127 25.406 25.406 91.601 91.220 96.767 93.127 93.124 92.342 92.335 61.343 61.259 52.061 51.910 29.109 29.109

2

2

4

4

2

10

2

10

a)

b)

Figure 10. Geometry and distributed nodes of a square plate with a complicated hole.

Mode 1, Omega = 21.4029

Mode 2, Omega = 35.3877

Mode 3, Omega = 39.0981

2

2

2

1

1

1

0

0

0

−1

−1

−1

−2 10

−2 10

−2 10

8

8

0 2

6 4

10

Mode 4, Omega = 62.7078

2

2

1

1

1

0

0

0

−1

−1

−1

−2 10

−2 10

−2 10

0 2 4

4 6

2

8 0

10

8

0 2

6

8

0

6

2

8 0

2

6

4

4

10

Mode 6, Omega = 93.1267

2

6

8 0

10

Mode 5, Omega = 68.7164

8

6

2

8 0

4

4

6

2

8 0

0 2

6

4

4

6

2

8

0 2

6

4

10

4

4 6

2

8 0

10

Figure 11. The first six mode shapes of the three-layer [45/-45/45] simply supported square plate with a complicated hole with a/h = 166.67 (the present-TSDT).

28

5.3 Buckling analysis 5.3.1. Uni-axial compression loading

Let us in this problem consider a four-layer cross-ply [0/90/90/0] laminated square plate with a simply supported boundary under axial compression loading. The buckling load factor of the plate is given by λ = λcr a 2 / ( E2 h3 ) , where λcr , a, E2 and h are the critical buckling load, the length, one of the elastic modulus of the material in second direction and the thickness of the plate, respectively. The material II is used. The plate is modeled by 588 nodes. The different values of length-to-thickness a/h and elastic modulus E1/E2 ratios are studied. Firstly, the length-to-thickness ratio is fixed equal to 10 (a/h = 10) with the modulus ratios E1/E2. A comparison of the present solution with the 3D elasticity solution [4] and several published results as a mesh-free solution based on HSDT [71] and a FEM solution based on HSDT [72, 73] are reported in Table 9. It is observed that the obtained numerical results show a good agreement when compared to those solutions. The normalized critical buckling load factor increases when increased to the modulus ratio. Next, the modulus ratio E1/E2 is fixed equal to 40 (E1/E2 = 40) and the length-to-thickness ratio is changed. Table 10 lists the normalized critical buckling load made of various length-to-thickness ratio a/h. It is found that the obtained results are in very close agreement with the values of FEM based on FSDT [74] and HSDT [75] and IGA based on layerwise [66]. Once again, the normalized critical buckling load factor increases when increased to length-to-thickness ratio. From this example, it can be concluded that obtained results from the present solution are very well compared to the 3D elasticity solution and other comparing solutions for the axial compression loading case. Table 9: Normalized critical buckling load λ of the four-layer [0/90/90/0] simply supported square plate with various E1/E2 ratios. Method Meshfree-HSDT [71] FEM-HSDT [72] FEM-HSDT [73] Elasticity 3D [4] IGA-layerwise [66] Present-TSDT Present-ITSDT

E1 /E2 3 5.412 5.114 5.442 5.294 5.3931 5.4040 5.4054

10 10.013 9.774 10.026 9.762 9.9453 9.9628 9.9699

20 15.309 15.298 15.418 15.019 15.3158 15.3226 15.3411

29

30 19.778 19.957 19.813 19.304 19.7089 19.6901 19.7195

40 23.412 23.34 23.489 22.881 23.3949 23.3405 23.3779

Table 10: Normalized critical buckling load λ of the four-layer [0/90/90/0] simply supported square plate with various ratios a/h. a/h Method/Authors 10 20 50 100 FEM-FSDT [74] 23.409 31.625 35.254 35.851 FEM-FSDT [75] 23.471 31.707 35.356 35.955 FEM-HSDT [75] 23.349 31.637 35.419 35.971 IGA-layerwise [66] 23.3949 31.6796 35.3505 35.9537 Present-TSDT 23.3405 31.7333 35.4808 36.1246 Present-ITSDT 23.3779 31.7590 35.4866 36.1261 To show further effectiveness of the current solution for the linear buckling analysis of the complex geometry, a three-layer simply supported square plates with a complicated cutout under the uni-axial compressive loading is studied. The material III is used. The buckling load factor is defined

(

)

by λˆ = λcr a 2 / π 2 Dˆ 0 , in which Dˆ 0 = E1h3 / 12 (1 − ν 12ν 21 ) . The geometry and distributed nodes are shown in Figure 12, respectively. Different values of angle ply and side-to-thickness ratio are investigated. For comparison, numerical results from the IGA based on CLPT for thin plates are given by Yin et al. [69]. Table 11 gives the normalized uni-axial buckling loading of the present method and the IGA one [69]. An excellent agreement is found for both side-to-thickness ratios and angle ply values.

b)

a)

Figure 12. Geometry and distributed nodes of a square plate with a complicated hole. Table 11: Normalized critical buckling load λ of the three-layer simply supported square plate with a complicated hole. a/h

Method

Angle plies

30

166.67

100 10 5

(0/0/0) 1.360 1.3559 1.3560 1.3461 1.3461 1.2092 1.2098 0.9877 0.9894

IGA-CLPT [69] Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT Present-TSDT Present-ITSDT

(15/-15/15) 1.47 9 1.4773 1.4773 1.4661 1.4661 1.3059 1.3059 1.0544 1.0548

(30/-30/30) 1.715 1.7181 1.7180 1.7044 1.7043 1.5074 1.5059 1.1951 1.1927

(45/-45/45) 1.827 1.8328 1.8328 1.8187 1.8186 1.6154 1.6130 1.2700 1.2656

(0/90/0) 1.359 1.3547 1.3547 1.3453 1.3452 1.2222 1.2200 1.0133 1.0084

5.3.2. Bi-axial compression loading

A three-layer [0/90/0] simply supported square plate under the biaxial buckling loading is finally considered. Similar to a previous example, various length-to-thickness and elastic modulus ratios are also studied to verify the buckling load factor. Table 12 and Table 13 give the normalized critical buckling loading corresponding with fixed the value of the length-to-thickness ratio and the elastic modulus ratio, respectively. Numerical results are compared with other solutions by Khdeir and Librescu [73], Fares and Zenkour [76], Thai et al. [66] and Liu et al. [71]. Again, it can be seen that a good agreement is found for both length-to-thickness and elastic

modulus ratios. Table 12: Biaxial critical buckling load λ of three-layer [0/90/0] simply supported square plate with various modulus ratios (a/h=10) Method/Authors Fares and Zenkour [76] Khdeir and Librescu [73] Thai et al. [66] Present-TSDT Present-ITSDT

E1 /E2 10 4.963 4.963 4.9145 4.9275 4.9278

20 7.588 7.516 7.4418 7.458199 7.45960

30 8.575 9.056 8.8473 8.8582 8.81694

40 10.202 10.259 10.0180 9.9977 9.93146

Table 13: Biaxial critical buckling load λ of three-layer [0/90/0] simply supported square plate with various ratios a/h (E1/E2=40) Method/Authors

2 Liu et al. [71] 1.457 Khdeir and Librescu [73] 1.465 Thai et al. [66] 1.5372 Present-TSDT 1.4528 Present-ITSDT 1.4540

a/h 10 10.251 10.259 10.0180 9.9977 9.93146

5 5.519 5.526 5.4813 5.4410 5.3630

31

15 12.239 12.226 12.0725 12.1217 12.08169

20 13.164 13.185 13.0684 13.15629 13.1309

6. Conclusions The consistent weak formulation with naturally stabilized nodal integration (NSNI) in the meshfree method based on HSDT have been studied for static, free vibration and buckling analyses of laminated composite and sandwich plates. Explicit formulations based on HSDT were achieved. The method avoided the instability in the direct nodal integration and reduces the computation cost as compared to traditional high-order Gauss quadrature technique. In comparison with the stabilized conforming nodal integration (SCNI) technique, the NSNI is directly calculated at nodes while the SCNI is indirectly computed through vertices (nodes) of Voronoi diagrams. Essential boundary conditions are directly imposed by the same way as the traditional finite element method due to satisfying the Kronecker delta function property of moving Kriging integration shape functions. The present plate formulations using the NSNI are easily extended for other meshfree methods by changing different approximation functions as the moving least-squares approximation, the reproducing kernel approximation, the radial point interpolation, etc. Through numerical examples, it was found that the present solution is very stable and accurate for analysis of the laminated composite and sandwich plates. Acknowledgement This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2016.19. References 1.

Pagano NJ. Exact solutions for rectangular bidirectional composites and sandwich plates. Journal of Composite Materials 1970; 4:20–34.

2.

Pagano NJ, Hatfield HJ. Elastic Behavior of Multilayered Bidirectional Composites. AIAA Journal 1972; 10:931–933.

3.

Noor AK. Free vibrations of multilayered composite plates. AIAA Journal 1973; 11:1038–1039.

4.

Noor AK, Mathers MD. Shear-flexible finite element method of laminated composite plate. Technical report, NASA; 1975.

5.

Srinivas S. A refined analysis of composite laminates. Journal of Sound and Vibration 1973; 30:495–507.

32

6.

Murakami H. Laminated composite plate theory with improved in-plane responses. Journal of Applied Mechanics 1986; 53:661–666.

7.

Reddy JN. A generalization of two-dimensional theories of laminated composite plates. Communications in Applied Numerical Methods 1987; 3:173–180.

8.

Reddy JN, Arciniega RA. Shear deformation plate and shell theories: from Stavsky to present. Mechanics of Advanced Materials and Structures 2004; 11:535–582.

9.

Reddy JN. An evaluation of equivalent-single-layer and layerwise theories of composite laminates. Composite Structures 1993; 25:21–35.

10. Carrera E. An assessment of mixed and classical theories on global and local response of multilayered orthotropic plates. Composite Structures 2000; 50:183–198. 11. Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Archives of Computational Methods in Engineering 2002; 9:87–140. 12. Reissner E. The effect of transverse shear deformation on the bending of elastic plates. Journal of Applied Mechanics, Transactions ASME 1945; 12:69–77. 13. Mindlin RD. Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates. Journal of Applied Mechanics, Transactions ASME 1951; 18:31–38. 14. Thai HC, Nguyen-Xuan H, Nguyen-Thanh N, Le T-H, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBSbased isogeometric approach. International Journal for Numerical Methods in Engineering 2012; 91:571–603. 15. Ferreira AJM, Castro LMS, Bertoluzza S. A high order collocation method for the static and vibration analysis of composite plates using a first-order theory. Composite Structures 2009; 89: 424–432. 16. Ambartsumian SA. On the theory of bending plates. Izv Otd Tech Nauk ANSSSR 1958; 5:269– 277. 17. Reddy JN. A simple higher-order theory for laminated composite plates. Journal of Applied Mechanics 1984; 51:745–752. 18. Thai HC, Nguyen-Xuan H, Bordas SPA, Nguyen-Thanh N, Rabczuk T. Isogeometric analysis of laminated composite plates using the higher-order shear deformation theory. Mechanics of Advanced Materials and Structures 2015; 22:451–469.

33

19. Thai HC, Zenkour AM, Wahab MA, Nguyen-Xuan H. A simple four unknown shear and normal deformations theory for functionally graded isotropic and sandwich plates based on isogeometric analysis. Composite Structures 2016; 139:77–95. 20. Nguyen-Xuan H, Thai HC, Nguyen-Thoi T. Isogeometric finite element analysis of composite sandwich plates using a higher order shear deformation theory. Composite Part B 2013; 55:558574. 21. Tuan N. Nguyen, D. Hui, J. Lee, H. Nguyen-Xuan. An efficient computational approach for sizedependent analysis of functionally graded nanoplates. Computer Methods in Applied Mechanics and Engineering 2015; 297:191-218. 22. Nguyen NT, Thai HC, Nguyen-Xuan H. On the general framework of high order shear deformation theories for laminated composite plate structures: A novel unified approach. International Journal of Mechanical Sciences 2016; 110:242–255. 23. Soldatos KP. A transverse shear deformation theory for homogenous monoclinic plates. Acta Mechanica 1992; 94:195–220. 24. Touratier M. An efficient standard plate theory. International Journal of Engineering Science 1991; 29:745–752. 25. Arya H, Shimpi RP, Naik NK. A zigzag model for laminated composite beams. Composite Structures 2002; 56:21-24. 26. Thai HC, Ferreira AJM, Rabczuk T, Bordas SPA, Nguyen-Xuan H, Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. European Journal of Mechanics - A/Solids 2014; 43: 89-108. 27. Thai HC, Kulasegaram S, Tran VL, Nguyen-Xuan H. Generalized shear deformation theory for functionally graded isotropic and sandwich plates based on isogeometric approach. Computers & Structures 2014; 141:94-112. 28. Karama M, Afaq KS, Mistou S. Mechanical behavior of laminated composite beam by new multi-layered laminated composite structures model with transverse shear stress continuity. International Journal of Solids and Structures 2003; 40:1525-1546. 29. Aydogdu M. A new shear deformation theory for laminated composite plates. Composite Structures 2009; 89:94-101. 30. Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin meshfree methods. International Journal for Numerical Methods in Engineering 2001; 50:435–466.

34

31. Puso M, Chen JS, Zywicz E, Elmer W. Meshfree and finite element nodal integration methods. International Journal for Numerical Methods in Engineering 2008; 74:416–446. 32. Hillman M, Chen JS, Chi SW. Stabilized and variationally consistent nodal integration for meshfree modeling of impact problems. Computational Particle Mechanics 2014; 1:245–256. 33. Hillman M, Chen JS. An accelerated, convergent, and stable nodal integration in Galerkin meshfree methods for linear and nonlinear mechanics. International Journal for Numerical Methods in Engineering 2016; 107: 603–630. 34. Beissel S, Belytschko T. Nodal integration of the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering 1996; 139:49–74. 35. Duan Q, Belytschko T. Gradient and dilatational stabilizations for stress-point integration in the element-free Galerkin method. International Journal for Numerical Methods in Engineerin 2009; 77: 776–798. 36. Nagashima T. Node-by-node meshless approach and its applications to structural analyses. International Journal for Numerical Methods in Engineering 1999; 46:341–385. 37. Liu GR, Zhang GY, Wang YY, Zhong ZH, Li GY, Han X. A nodal integration technique for meshfree radial point interpolation method (NI-RPIM). International Journal of Solids and Structures 2007; 44:3840–3890. 38. Wu CT, Chi SW, Koishi M, Wu Y. Strain gradient stabilization with dual stress points for the meshfree nodal integration method in inelastic analyses. International Journal for Numerical Methods in Engineering 2016; 107:3–30. 39. Wu CT, Koishi M, Hu W. A displacement smoothing induced strain gradient stabilization for the meshfree Galerkin nodal integration method. Computational Mechanics 2015; 56:19–37. 40. Chen JS, Hillman M, Chi SW. Meshfree Methods: Progress Made after 20 Years. Journal of Engineering Mechanics 2017; 143:1-38. 41. Reddy JN. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis, Second Edition, New York: CRC Press; 2003. 42. Thai HC, Do NVV, Nguyen-Xuan H. An improved moving Kriging meshfree method for analysis of isotropic and sandwich functionally graded material plates using higher-order shear deformation theory. Engineering Analysis with Boundary Elements 2016; 64:122-136. 43. Thai HC, Nguyen NT, Rabczuk T, Nguyen-Xuan H. An improved moving Kriging meshfree method for plate analysis using a refined plate theory. Computer & Structures 2016; 176:34-49.

35

44. Nguyen NT, Thai HC, Nguyen-Xuan H. A novel computational approach for functionally graded isotropic and sandwich plate structures based on a rotation-free meshfree method. Thin Walled Structures 2016; 107:473–488. 45. Liu GR. Meshfree Methods: Moving Beyond the Finite Element Method, CRC Press: U.S.A; 2003. 46. Liu WK, ONG JS, Uras RA. Finite element stabilization matrices - a unification approach. Computer Methods in Applied Mechanics and Engineering 1985; 53:13-46. 47. Koko J. A Matlab mesh generator for the two-dimensional finite element method. Applied Mathematics and Computation 2015; 250:650–664. 48. Liew KM, Huang YQ, and Reddy JN. Vibration analysis of symmetrically laminated plates based on FSDT using the moving least squares differential quadrature method. Computer Methods in Applied Mechanics and Engineering 2003; 192:2203–2222. 49. Chen XL, Liu GR, Lim SP. An element free Galerkin method for the free vibration analysis of composite laminates of complicated shape. Composite Structures 2003; 59:279–289. 50. Akhras G, Cheung MS, Li W. Finite strip analysis for anisotropic laminated composite plates using higher-order deformation theory. Computers and Structures 1994; 52:471–477. 51. Ferreira AJM, Roque CMC, Martins PALS. Analysis of composite plates using higher-order shear deformation theory and a finite point formulation based on the multiquadric radial basis function method. Composites: Part B 2003; 34:627–636. 52. Ferreira AJM. Analysis of composite plates using a layerwise theory and multiquadrics discretization. Mechanics of Advanced Materials and Structures 2005; 12:99–112. 53. Ferreira AJM, Fasshauer GE, Batra RC, Rodrigues JD. Static deformations and vibration analysis of composite and sandwich plates using a layerwise theory and RBF-PS discretizations with optimal shape parameter. Composite Structures 2008; 86: 328–343. 54. Roque CMC, Ferreira AJM, Jorge RMN. Modelling of composite and sandwich plates by a trigonometric layerwise deformation theory and radial basis functions. Composites: Part B 2005; 36:559–572. 55. Wang X, Shi G. A refined laminated plate theory accounting for the third-order shear deformation and inter-laminar transverse stress continuity. Applied Mathematical Modelling 2015; 39:5659– 5680. 56. Srinivas S. A refined analysis of composite laminates. Journal of Sound and Vibration 1973; 30:495–507.

36

57. Pandya BN, Kant T. Higher-order shear deformable theories for flexure of sandwich plates-finite element evaluations. International Journal of Solids and Structures 1988; 24: 419–451. 58. Mantari JL, Oktem AS, Soares CG. A new trigonometric shear deformation theory for isotropic, laminated composite and sandwich plates. International Journal of Solids and Structures 2012; 49: 43–53. 59. Grover N, Maiti DK, Singh BN. A new inverse hyperbolic shear deformation theory for static and buckling analysis of laminated composite and sandwich plates. Composite Structures 2013; 95: 667–675. 60. Kdheir A. Analysis of symmetric cross-ply elastic plates using a higher-order theory, part II: buckling and free vibration. Composite Structures 1988; 9:259–277. 61. Ferreira AJM. A formulation of the multiquadric radial basis function method for the analysis of laminated composite plates. Composite Structures 2003; 59:385–392. 62. Zhen W, Wanji C. Free vibration of laminated composite and sandwich plates using global-local higher-order theory. Journal of Sound and Vibration 2006; 298: 333–349. 63. Wu CP, Chen WY. Vibration and stability of laminated plates based on a local higher-order plate theory. Journal of Sound and Vibration 1994; 177:503–520. 64. Matsunaga H. Vibration and stability of cross-ply laminated composite plates according to a global higher-order plate theory. Composite Structures 2000; 48:231–244. 65. Cho KN, Bert CW, Striz AG. Free vibration of laminated rectangular plates analyzed by higherorder individual-layer theory. Journal of Sound and Vibration 1991; 145:429–442. 66. Thai HC, Ferreira AJM, Wahab MA, Nguyen-Xuan H. A generalized layerwise higher order shear deformation theory for laminated composite and sandwich plates based on isogeometric analysis. Acta Mechanica 2016; 227:1225-1250. 67. Thai HC, Ferreira AJM, Carrera E, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a layerwise deformation theory. Composite Structures 2013; 104:196–214. 68. Shojaee S, Valizadeh N, Izadpanah E, Bui TQ, Vu TV. Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method. Composite Structures 2012; 94:1677–1693. 69 Yin S, Yu T, Bui TQ, Xia S, Hirose S. A cutout isogeometric analysis for thin laminated composite plates using level sets. Composite Structures 2015; 127:152–164.

37

70. Bui TQ, Nguyen NM, Zhang C. An efficient meshfree method for vibration analysis of laminated composite plates. Computational Mechanics 2011; 48:175–193. 71. Liu L, Chua LP, Ghista DN. Mesh-free radial basis function method for static, free vibration and buckling analysis of shear deformable composite laminates. Composite Structures 2007; 78:58– 69. 72. Phan ND, Reddy JN. Analysis of laminated composite plates using a higher-order shear deformation theory. International Journal for Numerical Methods in Engineering 1985; 21:2201– 2219. 73. Khdeir AA, Librescu L. Analysis of symmetric cross-ply elastic plates using a higher order theory: Part II: buckling and free vibration. Composite Structures 1988; 9:259–277. 74. Chakrabarti A, Sheikh AH. Buckling of laminated composite plates by a new element based on higher order shear deformation theory. Mechanics of Composite Materials and Structures 2003; 10:303–317. 75.

Reddy JN, Phan ND. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher order shear deformation theory. Journal of Sound and Vibration 1985; 89: 157–170.

76. Fares ME, Zenkour AM. Buckling and free vibration of non-homogeneous composite cross-ply laminated plates with various plate theories. Composite Structures 1999; 44:279–287.

38