Analysis of composite plates using a layerwise theory and a differential quadrature finite element method

Accepted Manuscript Analysis of composite plates using a layerwise theory and a differential quadrature finite element method Bo Liu, A.J.M. Ferreira, Y.F. Xing, A.M.A. Neves PII: DOI: Reference:

S0263-8223(15)00647-9 http://dx.doi.org/10.1016/j.compstruct.2015.07.101 COST 6682

To appear in:

Composite Structures

Received Date: Accepted Date:

16 July 2015 17 July 2015

Please cite this article as: Liu, B., Ferreira, A.J.M., Xing, Y.F., Neves, A.M.A., Analysis of composite plates using a layerwise theory and a differential quadrature finite element method, Composite Structures (2015), doi: http:// dx.doi.org/10.1016/j.compstruct.2015.07.101

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.

Analysis of composite plates using a layerwise theory and a differential quadrature finite element method Bo Liu a, b, *, A. J. M. Ferreira b, Y. F. Xing a, A. M. A. Neves b, a

b

The Solid Mechanics Research Centre, Beihang University (BUAA), Beijing 100191, China. Departamento de Engenharia Mecânica, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal.

Abstract: A layerwise shear deformation theory for composite laminated plates is discretized using a differential quadrature finite element method (DQFEM). The DQFEM is a weak-form differential quadrature method that can provide highly accurate results using only a few sampling points. The layerwise theory proposed by Ferreira is based on an expansion of Mindlin’s first-order shear deformation theory in each layer and presented results for a laminated plate with three layers as example in the original paper. This work generalized the layerwise theory to plates with any number of layers. The combination of the DQFEM with Ferreira’s layerwise theory allows a very accurate prediction of the field variables. Laminated composite and sandwich plates are analyzed. The DQFEM solutions were compared with various models in literature and especially showed very good agreements with the exact solutions in literature that was based on a similar layerwise theory. The analysis of composite plates based on Ferreira’s layerwise theory indicates that the DQFEM is an effective method for high accuracy analysis of large-scale problems. Keywords: composite plate; layerwise theory; differential quadrature finite element method; bending; high accuracy

1. Introduction Composite and sandwich plates are one of the most significant applications of composite materials in industry. Layers are stacked together to form thin or thick laminates. When the main emphasis of the analysis is to determine the global response of the laminated component, relatively single-layer laminate theories (ESL theories) [1, 2] is accurate enough. The classical laminate plate theory, the first-order shear deformation theory [3-7], and higher-order theories [8] are commonly used examples of simple ESL theories. In some cases, particularly in sandwich applications, the difference between material properties makes it difficult for such theories to fully accommodate the bending behavior. Another set of theories that was introduced back in the 1980s are the layerwise theories, which consider independent degrees of freedom for each layer [9-16]. Layerwise displacement fields provide a much more correct representation of the moderate to severe cross-sectional warping associated with the deformation of thick laminates [2]. A very recent and comprehensive review of such theories in the analysis of multilayered plates and shells has been presented by Carrera [17]. This work adopts Ferreira’s layerwise theory [1] that is based on an expansion of Mindlin’s first-order shear deformation theory in each layer. The displacement continuity at layer’s interface is guaranteed. Also the theory directly produces very accurate transverse shear stress, although constant, in each layer.

*

Corresponding author: [email protected]

Most of the spacial discretization techniques thus far have been based on finite differences (FDM) and finite elements (FEM). Such low order schemes typically use low order basis functions and the accuracy is improved through mesh refinement. High order schemes like the hierarchical finite element method (HFEM) [18, 19], the radial basis functions (RBFs) [20, 21], the mesh free methods [22, 23], the differential quadrature method (DQM) [24], and more recently, the iso-geometric analysis (IGA) [25] and the differential quadrature method finite element method (DQFEM) [26-28], successively emerged as highly accurate numerical methods. All works via high order methods yield excellent results due to the use of the high-order or global basis functions. High order methods tend to give accurate results with far fewer degrees of freedom than low order schemes and have made noticeable success. This work adopts the DQFEM [26, 27] that is a weak-form differential quadrature method in essence. The DQFEM used the differential quadrature (DQ) rule and the Gauss-Lobatto quadrature rule to directly discretize the potential functional of structures to obtain the stiffness and mass matrices that are the same as those in the finite element method. The DQFEM has overcome the limitations of the DQM pointed out by Bert and Malik [24], and was hoped to be a competitive method with FEM for analysis of large-scale problems. This paper focuses for the first time on the analysis of composite laminated plates by the differential quadrature finite element method (DQFEM) [26, 27] and using a layerwise theory. This combination allows the accurate analysis of isotropic, composite, and sandwich plates of arbitrary shape and boundary conditions.

2. A layerwise theory The layerwise theory proposed by Ferreira [1] is used in this work. The theory is based on the assumption of a first-order shear deformation theory in each layer and the imposition of displacement continuity at the layer’s interfaces. Reference [1] considered a laminated plate with three layer for simplicity. This work generalized the theory to plates with any number of layers. The displacement of the i-th layer, according to first order shear deformation theory, can be written as

u i   x, y , z   u0i   x, y   z i  xi  v i   x, y , z   v0i   x, y   z i  yi  w

i 

(1)

 x, y , z   w  x, y 

The continuous of displacement u and v at the layer’s interfaces requires that

u0

i 1

h h  x, y   u0i   x, y   i  xi   i 1  xi 1

2 2 (2) h h i 1 i  i  i 1 i i 1 v0  x, y   v0  x, y    y  y 2 2 where hi are the i-th layer thickness and z(i)  [-hi/2, hi/2] are the i-th layer z coordinates. Using the recursion formula Eq. (2), one only needs to consider one layer. For example, in finite element method, after the stiffness and mass matrices of each layer are derived, the global stiffness and mass matrices of all layers can be obtained by using Eq. (2). Since laminated plates always have identical layers, the stiffness and mass matrices of identical layers only need to compute once.

 x i 1

z i 1 z i 

i 

x

z i 1

 i 1

x

Fig. 1. One-dimensional representation of the layerwise kinematics. For simplicity, in the following the superscript (i) is omitted, since the formulations for all layers are the same. The strain-displacement relations for i-th layer are given by

 xx   Dx    0  yy    xy    D y    0  yz    zx   0

0 Dy Dx 0 0

zDx 0 zD y 0 1

0 zD y zDx 1 0

0   u0  0   v0     0   x   D y   y    Dx   w 

(3a)

or (3b) ε = Du where Dx = /x and Dy = /y are differential operators. Note that D = D(z), namely, the differential operator matrix is a function of z. Neglecting z for each orthotropic layer, the stress–strain relations in the fiber local coordinate system can be expressed as

0 0   1    1  Q11 Q12 0   Q Q22 0 0 0   2   2   12    0 Q33 0 0   12   12    0     0 0 0 Q44 0   23   23     0 0 0 Q55   31   31   0

(4)

where subscripts 1 and 2 are the directions of the fiber and in-plane normal to fiber, respectively; subscript 3 indicates the direction normal to the plate; and the reduced stiffness components, Qij, are given by

Q11 

E1

1  1221

, Q22 

E2

1  1221

, Q12  12Q11

(5)

Q33  G12 , Q44  G23 , Q55  G31 , 21E1  12 E2 in which E1, E2, 12, 21, G12, G23, and G31 are material properties of lamina i. By performing adequate coordinate transformation, the stress–strain relations in the global x–y–z coordinate system can be obtained as

0   xx   xx  Q11 Q12 Q16 0    Q 0   yy   yy   12 Q22 Q26 0   0   xy    xy   Q16 Q26 Q66 0     0 0 0 Q44 Q45   yz  yz      0 0 Q45 Q55   zx    zx   0

(6a)

σ = Qε

(6b)

or

By considering θ as the angle between the x and 1 axes, where the 1 axis is the first principal material axis, usually connected with fiber direction, the components Qij can be calculated by adequate coordinate transformation (as in [2]). As in higher-order theories, this layerwise theory does not require the use of shear-correction factors. The principle of virtual displacements are required to derive the finite element matrices in this work. The virtual strain energy (δU) and the virtual work done by applied forces (δV) are given by

U   σ T εd   



 Du

T

Q  D u d

(7)

and

V    qT ud 

(8)

where  denotes the domain of the laminate, and q is the external distributed load. Usually Eqs. (7) and (8) are expanded by substituting Eqs. (3) and (6). This process will becomes complex with the increase of the number of displacements. Therefore, an alternative but equivalent way will be presented in next section. Discretize the virtual strain energy and virtual work using finite element method, one obtains the finite element formulation for i-th layer as

K i  ui   qi 

(9)

 u0i    i    v0    i  u  ψ xi    i   ψ y   w   

(10)

The displacement vector u(i) can be written as

where u0(i), v0(i),  x(i),  y(i) and w are the discrete values of u0(i), v0(i),  x(i), y(i) and w, respectively. In the present work, only symmetric laminates are considered; therefore, u0 and v0 of layer 2 can be discarded. For a three-layered laminated plate as considered in [1], the displacement field u0(i) and v0(i) of the first and third layers can be expressed as

h 1 h 2 1 u0   x, y    1  x   2  x  2 2 h h v01  x, y    1  y1  2  y 2  2 2

(11)

h2  2  h3  3 x  x 2 2 (12) h h  3 2  3 2 3 v0  x, y    y   y 2 2 Therefore, the displacement field of the three layers can be expressed by the global independent displacement fields as u0   x, y   3

 u01    1 h  1   2 1  v0   0       u0 3    0   3   0  v0         0  w  

 12 h2 0

0  12 h1

1 2

0 0 0

0 1  2 h2

h2 0

0 1 2 h2

0

0

0 0 1 2

0 0

h3 0

1 2

0

0 h3 0

 1  0  x  1 0  y   2    x     0  y 2   0   3   x     3  y  1    w 

(13)

or

u = Hu Therefore, the final stiffness and mass matrices of the three-layered laminated plates can be written as

 K 1  K  HT  0   0

0 K 2 0

 q1  0     0  H , q  H T q  2      3  K  3   q 

(14)

(15)

where HT and H are required by the principle of virtual displacements shown in Eqs. (7) and (8). It is clear that for layers of identical property, one only needs to compute the stiffness matrices of one layer.

3. The differential quadrature finite element method The differential quadrature finite element method (DQFEM) [26, 27] is a weak-form differential quadrature method (DQM) [24]. Namely, the strain energy and work potential of structures are directly discretized by the DQM together with the Gauss-Lobatto integration to obtain the global stiffness and mass matrices of finite element method. The DQM approximates the nth derivatives of a field variable f(x) at point xi by a weighted linear sum as N

f i  n    Aij n  f j

i  1,2,

,N

(16)

j 1

or

f  n   A n  f

(17)

where Aij(n) are the weighting coefficients of the nth order derivatives, and N the number of grid points in the x-direction. For the computation of weighting coefficients and more details about the DQM, one may refer to the survey paper [24]. The DQM for two dimensional problem can be expressed in similar way as Eq. (16) or (17) [24, 26, 27]. The Gauss-Lobatto quadrature is the Gauss integration with two endpoints fixed, which can be found in mathematics handbooks or in [19]. Here a simple introduction of it is presented to make the paper self-contained. The Gauss-Lobatto quadrature rule with precision degree (2n-3) for function f(x) defined at [-1, 1] is



1

1

f   dx   C j f  j  n

(18)

j 1

where the weights Cj of Gauss-Lobatto integration are given by

C1  Cn 

2 2 , Cj  n(n  1) n(n  1)[ Pn 1 ( j )]2

 j  1, n 

(19)

where  j is the ( j-1)th zero of Pn1 ( ) . Reference [19] presented the detail of computing roots for Legendre polynomials using the recursion formula of Legendre polynomials if more than 40 roots are required. Denote A(1) and B(1) as the two-dimensional DQM weighting coefficient matrix for first order derivative with respect to x and y, respectively. Using the differential quadrature method, the differential operator matrix D in Eq. (3) can be discretized as

 A1   0  D   B1  0   0

0

zA1

0

B 

0

zB  

A 

zB  

zA 

0

0

I

0

I

0

1

1

1

1

1

0   0   0  1 B    1 A  

(20)

where I is a unit matrix. Formulating the Gauss-Lobatto quadrature weight in the material constant matrix Q and denoting it as Q , then the stiffness matrix of i-th layer can be expressed as

K i   

hi /2

 hi /2

DTQDdz  

hi /2

 hi /2

K i   z  dz

(21)

Using the Gauss quadrature rule one only needs two Gauss points to do the integral in Eq. (21). Clearly, simpler than expanding Eq. (7). The obtaining of load vector q(i) is simple, since only Gauss-Lobatto quadrature is needed [27].

4. Numerical examples 4.1 Three-layer square sandwich plate under uniform load A simply supported square sandwich plate under a uniform transverse load is considered. This is a classical sandwich example of Srinivas [29]. The material properties of the sandwich core expressed in the stiffness matrix, Qcore , are expressed as

Qcore

0 0 0  0.999781 0.231192  0.231192 0.524886  0 0 0     0 0 0.262931 0 0   0 0 0 0.266810 0    0 0 0 0 0.159914 

(22)

Skins material properties are related to core properties by a factor R as follows:

Qskin  RQcore

(23)

Transverse displacement and stresses are normalized through the following factors:

w  w  a / 2, a / 2,0 

 1x   x3   y2   xz1 

0.999781 hq

 x1  a / 2, a / 2,  h / 2  q 2

x

,  x2 

 x1  a / 2, a / 2, 2h / 5  q

 a / 2, a / 2, 2h / 5 ,  1   y  a / 2, a / 2, h / 2  y 1

q

q

 y1  a / 2, a / 2, 2h / 5 q

 xz 2   0, a / 2,0  q

,  yz1 

,  y3 

(24)

 y 2   a / 2, a / 2, 2h / 5  q

 yz 2   a / 2,0,0  q

In order to clearly show how well the boundary conditions are satisfied, first the deformation modes of three of the seven independent variables in the right side of Eq. (13) are shown in Fig. 2. The stress resultants at z = ‒2h/5 calculated from layer 1 at z(1) = h1/2 are shown in Fig. 3. It is clear that the boundary conditions are satisfied very well.

 y1

 x1

w

Fig. 2. The deformation shape  x(1),  y(1) and w of the square plate with R = 5.

 y1

 x1

 yz1

 xy1

 xz1

Fig. 3. The stress resultants at z = ‒2h/5 calculated from layer 1 at z(1) = h1/2, R = 5. Transverse displacement and stresses for a sandwich plate are indicated in Tables 1-3 and compared with various formulations, where N is the number of sampling points used in each method for convergence studies. One can see the good agreement between the present results and the exact results. The present results agree with the deflection of exact results for 3 significant digits, agree with the normal stress of exact results for about 2 significant digits in

general, and agree with the shear stress of exact results for one significant digit. One can also see the fast convergence rate of the DQFEM from Tables 1-3. Even with a minimum number of sampling points, N =5, the present method still has good accuracy. The present method has about 3 to 5 significant digits converged when N = 15, while those in literature have only 1 to 2 significant digits converged. Table 1 Square laminated plate under uniform load (R = 5). Method

N

w

 1x

 x2

 x3

 1y

 y2

 y3

 xz1

 yz1

HSDT [30]

256.13

62.38

46.91

9.382

38.93

30.33

6.065

3.089

2.566

FSDT [30]

236.10

61.87

49.50

9.899

36.65

29.32

5.864

3.313

2.444

CLT

216.94

61.141

48.623

9.783

36.622

29.297

5.860

4.5899

3.386

Ferreira [31]

258.74

59.21

45.61

9.122

37.88

29.59

5.918

3.593

3.593

Ferreira [32]

15

257.38

58.725

46.980

9.396

37.643

27.714

4.906

3.848

2.839

HSDT [33]

11

253.671

59.6447

46.4292

9.2858

38.0694

29.9313

5.9863

3.8449

1.9650

HSDT [33]

15

256.239

60.1834

46.8581

9.3716

38.3592

30.1642

6.0328

4.2768

2.2227

HSDT [33]

21

257.110

60.3660

47.0028

9.4006

38.4563

30.2420

6.0484

4.5481

2.3910

Ferreira [1]

11

252.084

58.8628

45.4232

9.8846

37.6901

29.4765

5.8953

3.8311

2.5319

Ferreira [1]

15

255.920

59.6503

46.0366

9.2073

38.1408

29.8296

5.9659

3.9773

2.5375

Ferreira [1]

21

257.523

59.9675

46.2906

9.2581

38.3209

29.9740

5.9948

4.0463

2.3901

Present

5

260.791

62.3697

48.1774

9.6355

40.4410

31.6633

6.3327

4.2994

3.5372

Present

7

258.799

60.2134

46.4949

9.2990

38.4624

30.0898

6.0180

4.1297

3.4196

Present

11

258.828

60.2396

46.5201

9.3040

38.4895

30.1130

6.0226

4.1153

3.4118

Present

15

258.833

60.2466

46.5151

9.3030

38.4914

30.1118

6.0224

4.1080

3.4028

Present

25

258.835

60.2540

46.5098

9.3020

38.4951

30.1091

6.0218

4.1076

3.3994

258.97

60.353

46.623

9.340

38.491

30.097

6.161

4.3641

3.2675

Exact [29]

Table 2 Square laminated plate under uniform load (R = 10). Method

N

w

 1x

 x2

 x3

 1y

 y2

 y3

 xz1

 yz1

HSDT [30]

152.33

64.65

51.31

5.131

42.83

33.97

3.397

3.147

2.587

FSDT [30]

131.095

67.80

54.24

4.424

40.10

32.08

3.208

3.152

2.676

CLT

118.87

65.332

48.857

5.356

40.099

32.079

3.208

4.3666

3.7075

Ferreira [31]

159.402

64.16

47.72

4.772

42.970

42.900

3.290

3.518

3.518

Ferreira [32]

15

158.55

62.723

50.16

5.01

42.565

34.052

3.400

3.596

3.053

HSDT [33]

11

153.008

64.7415

49.4716

4.9472

42.8860

33.3524

3.3352

2.7780

1.8207

HSDT [33]

15

154.249

65.2223

49.8488

4.9849

43.1521

33.5663

3.3566

3.1925

2.1360

HSDT [33]

21

154.658

65.3809

49.9729

4.9973

43.2401

33.6366

3.3637

3.5280

2.3984

Ferreira [1]

11

155.037

63.5984

47.4765

4.7476

42.6696

32.7369

3.2737

3.7016

3.3051

Ferreira [1]

15

157.374

64.4828

48.1544

4.8154

43.1887

33.1392

3.3139

3.8447

3.2183

Ferreira [1]

21

158.380

64.8462

48.4434

4.8443

43.3989

33.3062

3.3306

3.9237

2.8809

Present

5

160.497

67.4660

50.4716

5.0472

45.7919

35.2056

3.5206

4.1910

3.5244

Present

7

159.371

65.1601

48.7223

4.8722

43.6173

33.4913

3.3491

4.0360

3.4160

Present

11

159.396

65.2015

48.7577

4.8758

43.6378

33.5098

3.3510

4.0040

3.3982

Present

15

159.403

65.2178

48.744

4.8744

43.6425

33.5061

3.3506

3.9957

3.3862

Present

25

159.406

65.2306

48.733

4.8733

43.6494

33.5002

3.3500

3.9957

3.3829

Exact [29]

159.38

65.332

48.857

4.903

43.566

33.413

3.500

4.0959

3.5154

Table 3 Square laminated plate under uniform load (R = 15). Method

N

w

 1x

 x2

 x3

 1y

 y2

 y3

 xz1

 yz1

HSDT [30]

110.43

66.62

51.97

3.465

44.92

35.41

2.361

3.035

2.691

FSDT [30]

90.85

70.04

56.03

3.753

41.39

33.11

2.208

3.091

2.764

CLT

81.768

69.135

55.308

3.687

41.410

33.128

2.209

4.2825

3.8287

Ferreira [31]

121.821

65.650

47.09

3.140

45.850

34.420

2.294

3.466

3.466

Ferreira [32]

15

121.184

63.214

50.571

3.371

45.055

36.044

2.400

3.466

3.099

HSDT [33]

11

113.594

66.3646

49.8957

3.3264

45.2979

34.9096

2.3273

2.1686

1.5578

HSDT [33]

15

114.387

66.7830

50.2175

3.3478

45.5427

35.1057

2.3404

2.6115

1.9271

HSDT [33]

21

114.644

66.9196

50.3230

3.3549

45.6229

35.1696

2.3446

3.0213

2.2750

Ferreira [1]

11

118.298

64.9159

46.8241

3.1216

45.4432

34.2237

2.2816

3.6123

3.8412

Ferreira [1]

15

120.077

65.8418

47.5260

3.1684

46.0049

34.6566

2.3104

3.7556

3.6695

Ferreira [1]

21

120.988

66.2911

47.8992

3.1933

46.2924

34.8898

2.3260

3.8311

3.2562

Present

5

122.549

68.9406

49.8964

3.3264

48.7625

36.8343

2.4556

4.1139

3.5367

Present

7

121.737

66.5786

48.1702

3.2113

46.5131

35.0824

2.3388

3.9683

3.4308

Present

11

121.764

66.6355

48.2074

3.2138

46.5281

35.0953

2.3397

3.9264

3.4062

Present

15

121.774

66.6633

48.1818

3.2121

46.5367

35.0878

2.3392

3.9190

3.3932

Present

25

121.777

66.6800

48.1664

3.2111

46.5463

35.0790

2.3386

3.9192

3.3580

121.72

66.787

48.299

3.238

46.424

34.955

2.494

3.9638

3.5768

Exact [29]

4.2 Three-layer (0/90/0) square cross-ply laminated plate under sinusoidal load A square laminate of length a and thickness h is composed of three equally thick layers oriented at (0/90/0). It is simply supported on all edges and subjected to a sinusoidal vertical pressure of the form

 x   y  pz  P sin   sin    a   a 

(25)

where the origin of the coordinate system is located at the lower-left corner on the mid-plane. For this example, there is a three dimensional exact solution by Pagano [34]. Here we compare the present solution by DQFEM for a/h = 10 with various models, particularly full mixed and hybrid finite element method (FEM) analysis, classical FEM analysis, etc. The material properties are

E1  25.0E2 , G12  G13  0.5E2 G23  0.2 E2 , 12  0.25

(26)

The numerical results are presented in Table 4, in a normalized form, as indicated by the following expressions:

102 w  a / 2, a / 2,0  h 3 E2  xx  a / 2, a / 2, h / 2  h 2 ,   xx Pa 4 Pa 2   a / 2, a / 2, h / 6  h 2   0, a / 2,0  h  yy  yy ,  zx  zx 2 Pa Pa   a / 2,0,0  h  yz  yz Pa w

(27)

As can be seen that, the present methodology converges to very good results, especially agree with the layerwise results of Ferreira [1] very well, since this work are based on the layerwise formulation of [1]. The convergence of the DQFEM in this case is even better than the case of section 4.1, as can be seen that the DQFEM already converged very well when the sampling points number N = 7.

Table 4 Laminated square plate (0/90/0) under sinusoidal load

w

 xx

 yy

 zx

 yz

Three-dimensional (Pagano [34])

0.7530

0.590

0.285

0.357

0.1228

Liou and Sun [35]

0.7546

0.580

0.285

0.367

0.127

Layerwise linear LD1 Carrera [36]

0.7371

0.5608

0.2740

0.3726

0.1338

Mixed layerwise LM4 Carrera [36]

0.7528

0.5801

0.2796

0.3626

0.1249

Reddy [8]

0.7125

0.5684

0.1033

Ferreira [1], layerwise (N = 11)

0.7399

0.5711

 0.2801

0.3560

 0.0872

Ferreira [1], layerwise (N = 15)

0.7420

0.5731

0.2808

0.3582

0.0931

Ferreira [1], layerwise (N = 21)

0.7427

0.5738

0.2810

0.3590

0.0953

Present (N = 5)

0.7391

0.5615

0.2758

0.3526

0.0923

Present (N = 7)

0.7402

0.5715

0.2806

0.3583

0.0958

Present (N = 11)

0.7402

0.5717

0.2807

0.3582

0.0958

Present (N = 15)

0.7402

0.5717

0.2807

0.3582

0.0958

Method

5. Conclusions The authors have used the differential quadrature finite element method (DQFEM) for vibration of isotropic beams, plates and three-dimensional elasticity [26, 27]. In all cases, excellent results that closely agree with exact results were obtained. In this paper, the analysis of composite laminated plates by the use of the DQFEM [26, 27] and using a layerwise theory [1] with independent rotations in each layer is performed for the first time. The layerwise theory [1] was also generalized to plates with any number of layers. The stiffness and mass matrices were schematically formulated. Composite laminated plates and sandwich plates were considered for testing of the present methodology and the results showed excellent accuracy for all cases. This layerwise theory combined with DQFEM discretization is a simple yet very effective and accurate numerical technique for the analysis of thick or thin composite or sandwich laminates and their structures.

Acknowledgments The financial support of National Natural Science Foundation of China (Grant Nos. 11402015, 11372021, 11172028), the specialized research fund for the doctoral program of higher education (20131102110039), FCT-Fundação para a Ciência e a Tecnologia to Project PTDC/EMS-PRO/2044/2012 and to grant SFRH/BPD/99591/2014 is gratefully acknowledged.

Reference [1] Ferreira AJM. Analysis of Composite Plates Using a Layerwise Theory and Multiquadrics Discretization. Mechanics of Advanced Materials and Structures 2005;12:99-112. [2] Reddy JN. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. Second ed: CRC Press 2003.

[3] Whitney JM, Pagano NJ. Shear Deformation in Heterogeneous Anisotropic Plates. Journal of Applied Mechanics 1970;37:1031-1036. [4] Whitney JM. The Effect of Transverse Shear Deformation on the Bending of Laminated Plates. Journal of Composite Materials 1969;3:534-547. [5] Reissner E. On transverse bending of plates, including the effect of transverse shear deformation. International Journal of Solids and Structures 1975;11:569-573. [6] Reissner E. A Consistent Treatment of Transverse Shear Deformations in Laminated Anisotropic Plates. AIAA Journal 1972;10:716-718. [7] Mindlin RD. Influence rotatory inertia and shear in flexural motion of isotropic, elastic plates. ASME Journal of Applied Mechanics 1951;18:31–38. [8] Reddy JN. A Simple Higher-Order Theory for Laminated Composite Plates. Journal of Applied Mechanics 1984;51:745-752. [9] Mau ST. A Refined Laminated Plate Theory. Journal of Applied Mechanics 1973;40:606-607. [10] Chou PC, Carleone J. Transverse shear in laminated plate theories. AIAA Journal 1973;11:1333-1336. [11] Di Sciuva M. An Improved Shear-Deformation Theory for Moderately Thick Multilayered Anisotropic Shells and Plates. Journal of Applied Mechanics 1987;54:589-596. [12] Murakami H. Laminated Composite Plate Theory With Improved In-Plane Responses. Journal of Applied Mechanics 1986;53:661-666. [13] Ren JG. A new theory of laminated plate. Composites Science and Technology 1986;26:225-239. [14] Carrera E. C0 reissner–mindlin multilayered plate elements including zig-zag and interlaminar stress continuity. International Journal for Numerical Methods in Engineering 1996;39:1797-1820. [15] Carrera E, Kröplin B. Zigzag and interlaminar equilibria effects in large-deflection and postbuckling analysis of multilayered plates. Mechanics of Composite Materials and Structures 1997;4:69-94. [16] Carrera E. Evaluation of Layerwise Mixed Theories for Laminated Plates Analysis. AIAA Journal 1998;36:830-839. [17] Carrera E. Historical review of Zig-Zag theories for multilayered plates and shells. Applied Mechanics Reviews 2003;56:287-308. [18] Zienkiewicz OC, De S.R. Gago JP, Kelly DW. The hierarchical concept in finite element analysis. Computers & Structures 1983;16:53-65. [19] Liu B, Xing Y, Wang W, Yu W. Thickness-shear vibration analysis of circular quartz crystal plates by a differential quadrature hierarchical finite element method (in press). Composite Structures. [20] Kansa EJ. Multiquadrics—A scattered data approximation scheme with applications to computational fluid-dynamics—I surface approximations and partial derivative estimates. Computers & Mathematics with Applications 1990;19:127-145. [21] Kansa EJ. Multiquadrics—A scattered data approximation scheme with applications to computational fluid-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications 1990;19:147-161. [22] Liu GR. Meshfree Methods: Moving Beyond the Finite Element Method, Second Edition: CRC Press 2009. [23] Liu GR, Trung NT. Smoothed Finite Element Methods: CRC Press 2010. [24] Bert CW, Malik M. Differential Quadrature Method in Computational Mechanics: A Review. Applied Mechanics Reviews 1996;49:1-28. [25] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 2005;194:4135-4195. [26] Xing Y, Liu B. High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain. International Journal for Numerical Methods in Engineering 2009;80:1718-1742.

[27] Xing Y, Liu B, Liu G. A differential quadrature finite element method. International Journal of Applied Mechanics 2010;2:207-227. [28] Zhong H, Yu T. A weak form quadrature element method for plane elasticity problems. Applied Mathematical Modelling 2009;33:3801-3814. [29] Srinivas S. A refined analysis of composite laminates. Journal of Sound and Vibration 1973;30:495-507. [30] 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:1267-1286. [31] Ferreira AJM, Barbosa JT. Buckling behaviour of composite shells. Composite Structures 2000;50:93-98. [32] Ferreira AJM. A formulation of the multiquadric radial basis function method for the analysis of laminated composite plates. Composite Structures 2003;59:385-392. [33] 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: Engineering 2003;34:627-636. [34] Pagano NJ. Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates. Journal of Composite Materials 1970;4:20-34. [35] Liou W-J, Sun CT. A three-dimensional hybrid stress isoparametric element for the analysis of laminated composite plates. Computers & Structures 1987;25:241-249. [36] Carrera E. Developments, ideas, and evaluations based upon Reissner’s Mixed Variational Theorem in the modeling of multilayered plates and shells. Applied Mechanics Reviews 2001;54:301-329.