A linear smoothed quadratic finite element for the analysis of laminated composite Reissner–Mindlin plates

Accepted Manuscript A linear smoothed quadratic finite element for the analysis of laminated composite Reissner–Mindlin plates Detao Wan, Dean Hu, Sundararajan Natarajan, Stéphane P.A. Bordas, Ting Long PII: DOI: Reference:

S0263-8223(17)31853-6 http://dx.doi.org/10.1016/j.compstruct.2017.07.092 COST 8748

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

13 June 2017 24 July 2017 28 July 2017

Please cite this article as: Wan, D., Hu, D., Natarajan, S., Bordas, S.P.A., Long, T., A linear smoothed quadratic finite element for the analysis of laminated composite Reissner–Mindlin plates, Composite Structures (2017), doi: http://dx.doi.org/10.1016/j.compstruct.2017.07.092

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

A linear smoothed quadratic finite element for the analysis of laminated composite Reissner–Mindlin plates Detao Wan1, 2, Dean Hu1, 2*, Sundararajan Natarajan3, Stéphane P.A. Bordas4 and Ting Long1, 2 1

2

State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, P. R. China Key Laboratory of Advanced Design and Simulation Technology for Special Equipments, Ministry of Education, Hunan University, Changsha, 410082, P. R. China 3

4

Department of Mechanical Engineering, Indian Institute of Technology, Madras, Chennai 600036, India Research Unit in Engineering Science, Luxembourg University, 6 rue Richard Coudenhove-Kalergi, L-1359 Luxembourg, Luxembourg *

E-mail address: [email protected] (Dean Hu)

Abstract: It is well known that the high-order elements have significantly improved the accuracy of solutions in the traditional finite element analysis, but the performance of high-order elements is restricted by the shear-locking and distorted meshes for the plate problems. In this paper, a linear smoothed eight-node Reissner-Mindlin plate element (Q8 plate element) based on the first order shear deformation theory is developed for the static and free vibration analysis of laminated composite plates, the computation of the interior derivatives of shape function and isoparametric mapping can be removed. The strain matrices are modified with a linear smoothing technique by using the divergence theorem between the nodal shape functions and their derivatives in Taylor’s expansion. Moreover, the first order Taylor’s expansion is also employed for the construction of stiffness matrix to satisfy the linear strain distribution. Several numerical examples indicate that the novel Q8 plate element has good performance to alleviate the shear-locking phenomenon and improve the quality of the solutions with distorted meshes.

Keywords: laminated composite plate; linear smoothing technique; Taylor’s expansion; shear-locking; distorted meshes.

1. Introduction Owing to the remarkable weight-to-stiffness, strength-to-stiffness characteristics and the advantage of designable, the laminated composite plates are widely used in engineering structures as diverse as aerospace, aircrafts, automotive structural parts, civil engineering structures, etc. Many plate theories have been successfully applied to analyze the laminated composite plates, such as the classical laminated plate theory [1-3], first-order shear deformation theory [4-11] and high-order shear deformation theory [12-15]. Although many analytical approaches based on different plate theories are able to give deep physical insights and highly accurate results, but they are restricted to the problems with simple geometries and boundary conditions. In the field of computational mechanics, various advanced numerical methods have been developed to analyze the composite plates, such as the finite element method (FEM) [16-18], boundary element method (BEM) [19, 20], mesh-free methods [21-24] and isogeometric analysis (IGA) [9, 25-29]. Among these numerical methods, FEM is considered to be a powerful and effective numerical approach for the plate analysis. The quadratic isoparametric Mindlin plate elements in the FEM are widely used for the plate analysis owning to their higher rate of convergence, highly accurate solutions and their ability to model the complex boundaries. On the other hand, the quadratic isoparametric Mindlin plate element also has one significant deficiency [30-32], viz., the performance of this type elements deteriorates rapidly in thin plate structure with shear locking phenomenon because of the excessive stiffness on account of the transverse displacement constrains. Many remarkable remedial schemes have been developed for relieving the shear-locking phenomenon, such as the reduced or selective integration technique [32, 33], the addition of nonconforming displacement modes [34, 35] and the assumed shear strain fields [36-38]. Choi et al. [34, 35, 39] developed the eight-node Mindlin plate elements by coupling the addition of nonconforming displacement modes, the assumed shear strain fields and the selectively reduced integration to eliminate the shear-locking phenomenon. However, all the proposed schemes are based on the framework of conventional FEM and the computation of interior derivatives of shape function is still required. In addition, the isoparametric mapping is an unavoidable procedure for the conventional FEM. Thus, their quality of the solutions can be broken down for the extremely distorted meshes. For the purpose of improving the quality of finite element solutions and relieving the effect of mesh distortion, Liu et al. [40] proposed a smoothed finite element method (S-FEM), which is based on the stabilized conforming nodal integration (SCNI) of mesh-free method [41]. Then, the node-based smoothed finite element (NS-FEM) [42], the edge-based smoothed finite element (ES-FEM) [43] and the face-based

smoothed finite element (FS-FEM) [44] have been proposed. The S-FEM has been successfully introduced for the analysis of plate and shell [45-52], nearly incompressible solids [53, 54] and other problems. The S-FEM avoids the evaluation of the Jacobian matrix, removes the isoparametric mapping and reduces the mesh sensitivity to some extent. It is a very robust, accurate and computation inexpensive numerical method. It should be noted that all the types of S-FEM use the finite element with linear interpolation, and the strain smoothing technique over the higher order elements exhibits poor performance [55]. The strain smoothing technique used in the conventional S-FEM is generally referred as the constant smoothing technique. Recently, Francis et al. [56] proposed a linear strain smoothing technique in the framework of CS-FEM to improve accuracy of arbitrary convex polytopes with linear interpolation or quadratic interpolation based on the recent work of Duan et al. [57-59], which also has been coupled with the extended finite element method [60]. It was also reported that the linear smoothing technique yields an improved accuracy compared to the constant smoothing technique, and reproduces a linear strain field with quadratic interpolation. The application of the linear smoothing technique for the plate elements has not been reported in any literature. With the aim of eliminating the shear locking phenomenon, relieving the effect of mesh distortion on the quality of finite element solutions and keeping the higher rate of convergence, a linear smoothed quadratic finite elements based on the first-order shear deformation theory is developed and applied for the static and free vibration analysis of laminated composite Reissner–Mindlin plates in this paper, which is the further application of the linear smoothing technique. This also eliminates the need for isoparametric mapping and the computation of interior derivatives of shape function, in addition, all the computation are based on the global Cartesian coordinates. The paper is constructed as follows. The basic formulations of laminated composite plates are given in Section 2. The detailed integration scheme with the linear smoothing technique is presented in Section 3. In Section 4, several numerical examples are presented to demonstrate the performance of the linear smoothed quadratic finite element. The concluding remarks are given in the last section. 2. Basic formulations of laminated composite plates Consider a laminated composite plate with thickness h , as shown in Fig. 1, in which nl number of layers are included and regarded as the orthotropic materials, each layer of laminate can be arbitrarily oriented at an angle θ in reference to x-axis of the coordinate system. According to the first order shear deformation theory, the displacement field of the laminated composite plates can be expressed as

u ( x, y, z ) = u0 ( x, y ) − z β x ( x, y)  v( x, y, z ) = v0 ( x, y ) − z β y ( x, y)  w( x, y, z ) = w ( x, y )  0

(1)

where u0 , v0 and w0 are the displacements of the mid-plane of the plate in x , y and z directions,

β x and β y donate the rotations in y -z and x -z planes.

(a) Coordinate system

(b) Layer details Fig. 1. Laminated composite plate.

The strain vector ε can be written as

ε xx   u0, x   β x , x        b ε yy  =  v0, x  − z  β y , y  = ε 0 − zκ = ε γ  u + v   β + β  y,x   xy   0, y 0, x   x , y

(2)

γ xz   β y − w0, x  s  =  = γ0 = ε γ yz   β x − w0, y 

(3)

and

The strain energy U ε is given by

Uε =

1 ε ΤSεdΩ ∫ Ω 2

(4)

in which

ε   zκ   0  ε =  0− +   0   0  γ 0   A11 A  12  A16  B S =  11  B12   B16  0   0

A12 A22

A16 A26

B11 B12

B12 B22

B16 B26

0 0

A26 B12 B22

A66 B16 B26

B16 D11 D12

B26 D12 D22

B66 D16 D26

0 0 0

B26 0

B66 0

D16 0

D26 0

D66 0

0 D44s

0

0

0

0

0

s D45

0  0  0   A B 0  0   D 0  = B D 0 =    0  0 Ds    0 0 Ds   0  D45s   D55s 

(5)

(6)

where Aij , Bij , Dij and Dijs represent the extensional, coupling, bending and the transverse shear stiffness, respectively, and those can be given as nl zk  2 A , B , D = ( ij ij ij ) ∑ ∫zk −1 Qij k (1, z , z ) dz , i, j = 1, 2, 6  k =1  nl zk ( D s ) = ∑ ∫z ζ Qij k dz, i, j = 4,5  ij k =1 k −1

( )

(7)

( )

in which

ζ = 5 6 denotes the transverse shear correction coefficient.

stiffness matrix for the k

th

(Q )

ij k

is the transformed reduced

layer of orthotropic laminate and given by

(Q )

ij k

 m2  =  n2  mn 

n2 m

2

−mn

 m2 n2 −2mn    m2 2mn  ( Qij ) k  n 2  − mn mn m 2 − n 2  

2mn   −2mn  m2 − n 2 

(8)

and

m = cos θ k , n = sin θ k where θ k in Eq. (9) is the angle of k

th

layer in reference to x -axis.

(9)

(Qij )k

is the reduced stiffness

matrix for k th layer of orthotropic laminate and the definition can be found in Ref. [61]. For the static analysis of laminated composite plates, the total potential energy can be expressed as

Π s = U ε − ∫ u T bdΩ − ∫ u T tdΓ Ω

Γt

(10)

where t is the external traction and b is the body force. For the free vibration analysis of laminated composite plates, the total energy is given by

Π v = Uε − where

1 ρ h ( u& 2 + v&2 + w& 2 ) dΩ 2 ∫Ω

(11)

ρ is the density.

For the static analysis with Q8 plate elements, the approximation of displacements and rotations of the mid-plane of the plate can be given as

 u0h   uI   h    v0  8  vI  8 h h u 0 =  w0  = ∑ N I ( x )  wI  = ∑ N I ( x ) u I  h  I =1   I =1  βx   β Ix  βh   β Iy     y

(12)

where u I is the vector of degrees of freedom, N I ( x) is the non-mapped Lagrange shape function based on the global Cartesian coordinates and expressed as

N(x) = P −1p ( x, y )

(13)

and

xy 2 ]

(14)

P = p ( x1 , y1 ) , p ( x2 , y2 ) , K , p ( x8 , y8 ) 

(15)

p( x, y ) = [1 x

x2

y

y2

xy

x2 y

in which N I (x) is composed of polynomial and can be written as

N I (x ) = a1I + a2 I x + a3 I y + a4 I x 2 + a5 I xy + a6 I y 2 + a7 I x 2 y + a8 I xy 2

(16)

here xI , yI ( I = 1L 8 ) are the global Cartesian coordinates of nodes of each element. a1I , a2 I ,…,

a8 I are the corresponding elements of P −1 in Eq. (13). Substituting Eq. (12) into Eq. (10) and taking variation to the energy function, the equilibrium equation is obtained

Ku = f

(17)

K = Kb + K s

(18)

f = ∫ N ΤbdΩ + ∫ N Τ tdΓ

(19)

K bIJ = ∫ BbI T DBbJ dΩ

(20)

K sIJ = ∫ B sI T Ds B sJ dΩ

(21)

where

Ω

Γ

and Ω

Ω

b

s

in which B I and B I are respectively given by

 N I ,x  0  NI,y B bI =   0  0   0

0

0

0

NI ,y NI,x

0 0

0 0

0 0

0 0

NI ,x 0

0

0 NI ,y

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

−N I 0

0  0  0   0  NI,y   N I , x  0  − N I 

(22)

(23)

For the free vibration analysis, the approximation of displacements and rotations of the mid-plane of the plate can be expressed as

8

u 0h == ∑ N I ( x ) u I eiωt

(24)

I =1

Substituting Eq. (24) into Eq. (11) and taking variation to the energy function, the Eigen equation can be obtained by

( K − ω 2M ) u = 0

(25)

where K has the same expression as that of Eq. (18). The mass matrix is express as

M = ∫ N ΤmNdΩ

(26)

Ω

and

1 0  m = ρ h 0  0 0

     2 0 0 h 12 0  0 0 0 h 2 12  0 0 1 0 0 1

0 0 0

0 0 0

(27)

3. Integration scheme with linear smoothing technique

% b and The linear smoothing technique is employed in the process of computing the modified strain B I B% sI to construct the smoothed stiffness matrix, respectively. According to the framework of conventional CS-FEM, the Q8 plate element is divided into eight smoothing cells, as shown in Fig. 2, the eight smoothing cells are respectively constructed by points 1-5-9, 5-2-9, 2-6-9, 6-3-9, 3-7-9, 7-4-9, 4-8-9, 8-1-9. The boundary integral along the curve edge is approximated by dash line, as shown in Fig. 2(b).

(a) With straight edges (b) With a curve edge Fig. 2. Division of a Q8 plate element into eight smoothing cells. In order to obtain the modified strain field ε%

q(x) = [1 x

b

in each smoothing cell, the linear polynomial

y ]Τ is used, then the modified strain is given by

ε% b ( x k ) = ∫ ε b ( x)q( x)dΩ Ωs

(28)

In the process of computing the modified strain matrix B bI , the consistency form should be met by the terms related to the shape function N I ( x) and the derivative of shape function N I ,i ( x) (i = x, y) . According to Eq. (28) and implementing the divergence theorem, we obtain

∫

Ωs

N I ,iq ( x ) dΩ = ∫ N I q ( x ) ni dΓ − ∫ N I q ,i ( x ) dΩ Γs

Ωs

(29)

where ni ( i = x, y ) is the outward normal. The derivatives of q ( x ) are

q, x ( x ) = [ 0 1 0 ] and q , y ( x ) = [ 0 0 1]

(30)

The evaluation of the modified derivatives N% I , x (x ) is first presented in this section. To begin with, according to Ref. [59], N% I , x ( x) can be computed by only one evaluation point with high order derivatives, viz.,

{N%

I,x

}

(x c ), N% I , xx (x c ), N% I , xy (x c ) , where x c = ( xc , yc ) is the interior evaluation point

(central point) of the smoothing cell. The Taylor’s expansion is employed, and the expanded form of

N% I , x ( x) and q ( x ) are, respectively, expressed as N% I , x (x ) = N% I , x ( x c ) + ( x − xc ) N% I , xx ( x c ) + ( y − yc ) N% I , xy ( x c ) + H.O.T

(31)

q( x) = q( x c ) + ( x − xc ) q, x ( x c ) + ( y − yc ) q, y ( x c ) + H.O.T

(32)

where H.O.T is the high order term that can be ignored. N I , x (x ) is replaced by N% I , x (x ) , then, substituting Eqs. (31) and (32) into Eq. (29), we obtain q ( x c ) A N% I , x ( xc ) + q, x ( x c ) I cxx + q, y ( x c ) I cxy  N% I , xx ( xc ) + q, x ( x c ) I cxy + q, y ( xc ) I cyy  N% I , xy ( xc ) = ∫ N I ( x ) q ( x ) nx dΓ − ∫ Γs

Γs

(∫ N

I

( x ) q, x ( x ) dx ) nx dΓ

(33)

where A is the area of the smoothing cell, and

 I cxx = ( x − xc ) 2 dΩ ∫Ωs   I xy =  c ∫Ωs ( x − xc )( y − yc ) dΩ  2  I cxy = ∫ ( y − yc ) dΩ Ωs 

(34)

x

Eq. (34) is the second order area moments of the smoothing cell Ω s , and the first order area moments ( I c ,

I cy ) with the center xc of smoothing cells are zero [59]. In addition, the domain integral of the second term in the right hand side of Eq. (29) has been transformed into line integral and expressed in Eq. (33), the detail theory about the transformation can be found in Refs. [62-64]. The indefinite integral within the second term of Eq. (33) can be calculated easily by analytical integral, because the matrix of shape function

N I ( x ) is composed by polynomial and q,x ( x ) is always a constant. Note that in such a way, the domain integral in Eq. (34) also can be transformed into boundary integral.

Substituting Eq. (30) into Eq. (33) and two Gauss points on each edge are used for evaluation of the boundary integral, the expanded form is expressed as

 A  Ax  c  Ayc

0

I cxx I

xy c

0   N% I , x ( x c )   F11x      I cxy   N% I , xx ( x c )  =  F12x  I cyy   N% I , xy ( x c )   F13x 

(35)

where N eg 2   N I ( x kG ) nxkWG   ∑∑ k =1 G =1   x  F11  N   x   eg 2 k k k k  F12  = ∑∑ N I ( xG ) xG − ∫ N I dx ( x G ) nx WG   x   k =1 G =1   F13   N eg 2    N I ( x kG ) yGk nxkWG ∑∑   i =1 G =1

(

(

)

)

(36)

and N eg is the number of edges of the smoothing cell, WG is the weight corresponding to Gauss integral point x Gk on each boundary segment Γ ks , nxk is the outward normal on Γ ks . It can be observed from Eq. (35) that the linear smoothing technique essentially contains three equations corresponding to the smoothing function q( x) = [1

x

y ]Τ , then the modified derivatives N% I , x (xc ) , N% I , xx (xc ) and

N% I , xy ( xc ) can be obtained by analytical solutions of Eq. (35), we have

% F11x  N I , x (x c ) = A   I cyy ( F12x − F11x xc ) − I cxy ( F13x − F11x yc )  N% I , xx ( x c ) = 2  I cxx I cyy − ( I cxy )   I xx ( F x − F11x yc ) − I cxy ( F12x − F11x xc )  N% I , xy (x c ) = c 13 2  I cxx I cyy − ( I cxy )  here A > 0 and I cxx I cyy > ( I cxy )

2

(37)

ensure that the Eq. (35) can always has the unique solution [59].

In such a way, the modified derivatives N% I , y ( x c ) , N% I , yx ( x c ) and N% I , yy (x c ) can also be given as

 A  Ax  c  Ayc where

0 I cxx

I

xy c

0   N% I , y ( x c )   F11y      I cxy   N% I , yx ( x c )  =  F12y   I cyy   N% I , yy ( x c )   F13y 

(38)

N eg 2   N I ( x Gk ) n ykWG   ∑∑ k =1 G =1 y    F11  N eg 2   y   N I ( x kG ) xGk n kyWG  F12  =   ∑∑ k =1 G =1  y    F13   Neg 2   ∑∑ N I ( x kG ) yGk − ∫ N I dy ( xGk ) n ykWG   k =1 G =1 

(

(

)

(39)

)

the analytical solutions of Eq. (38) are

% F11y  N I , y (x c ) = A   I cyy ( F12y − F11y xc ) − I cxy ( F13y − F11y yc )  N% I , yx ( x c ) = 2  I cxx I cyy − ( I cxy )   I xx ( F13y − F11y yc ) − I cxy ( F12y − F11y xc )  N% I , yy (x c ) = c xx yy xy 2  I I − I ( ) c c c 

(40)

% b and the modified derivatives of the strain matrix ∂B% b ∂x , After that, the modified strain matrix B I I ∂B% bI ∂x can be respectively expressed as

 N% I , x   0  N% B% bI =  I , y  0  0   0

0 N% I , y N%

 N% I , xx   0 b ∂B% I  N% I , yx = ∂x  0  0   0

0 % N I , yx N%

 N% I , xy   0 ∂B% bI  N% I , yy = ∂y  0  0   0

0 0

0 0

0

0 0 % 0 NI,x

0 0

0 0 0 N% I , y

I,x

0 0

0 0

0

0 0 0 N% I , xx

0 0

0 0 % 0 N I , yx

0

0 0

I , xx

N% I , yy N%

0 0

0

0 0 % 0 N I , xy

0 0

0 0 0 N% I , yy

I , xy

0   0  0  , 0  N% I , y   N% I , x 

I = 1L8

(41)

0   0  0  , 0  N% I , yx   N% I , xx 

I = 1L8

(42)

0   0  0  , 0  N% I , yy   N% I , xy 

I = 1L8

(43)

% s and the modified derivatives of the shear strain In the same way, the modified shear strain matrix B I % s ∂x , ∂B% s ∂y are given as matrix ∂B I I

0 0 N% I , x B% sI =  0 0 N% I , y

− N% I 0

0  , − N% I 

I = 1L8

(44)

∂B% sI 0 0 N% I , xx = ∂x 0 0 N% I , yx

− N% I , x 0

0  , − N% I , x 

I = 1L 8

(45)

∂B% sI 0 0 N% I , xy = ∂y 0 0 N% I , yy

− N% I , y 0

0  , − N% I , y 

I = 1L 8

(46)

where N% I is computed at the central of the smoothing cell, N% I , x , N% I , y , N% I , xx , N% I , yy , N% I , xy and

N% I , yx are given in Eqs. (37) and (40). Based on the framework of CS-FEM, the smoothed element stiffness matrices are computed by the sum of the contributions of the smoothing cells and given by

) D B% dΩ )

( K% ) = ∑ ( ∫

B% bI T DB% bJ dΩ

(47)

( K% ) = ∑ ( ∫

% sT B I

(48)

nc

b IJ

Ωs

C =1 nc

s IJ

C =1

Ωs

s

s J

where nc is the number of smoothing cells within an element （e.g., as shown in Fig.2, nc = 8 ）. In order

%b to guarantee the linear field of the strain matrix in the process of construction the stiffness matrices K IJ % s , the modified derivatives of strain matrix of Eqs. (42) and (43) should be introduced into the final and K IJ discretized equation. Therefore, the first-order Taylor’s expansion with respect to the center x c of each smoothing cell is applied to the Eq. (47), leads to

( K% ) = ∑ ( ∫ nc

b IJ

C =1

Ωs

B% b Τ DB% b dΩ

)

nc  %b  % b Τ ∂B% b Τ   % b ∂B% b   ∂B% b Τ ∂B = ∑  ∫  B + + x − xc ) + y − yc )  D  B x − xc ) + y − yc )  dΩ  (49) ( ( ( ( Ωs ∂x ∂y ∂x ∂y C =1       b Τ b b Τ b b Τ b b Τ nc  % % % % % % % %b  % b Τ DB % b Τ + I xx ∂B D ∂B + I yy ∂B D ∂B + I xy  ∂B D ∂B + ∂B D ∂B   = ∑  AB c c c ∂x ∂x ∂y ∂y ∂y ∂y ∂x   C =1   ∂x

where

B% b = [B% 1b

B% b2 L B% 8b ]

% b , ∂B% b ∂x , ∂B% b ∂y in Eq. (49) can be determined by Eqs. (41)-(43). and B I I I Similarly, Eq. (48) can be rewritten as

(50)

nc

( K% ) = ∑( ∫ s IJ

C =1

Ωs

)

% s ΤDs B % s dΩ B

nc  % sΤ % sΤ %s %s  % sΤ ∂B   % s ∂B   ∂B ∂B = ∑ ∫ B + x − xc ) + y − yc )  Ds B + x − xc ) + ( ( ( ( y − yc )  dΩ  (51) Ωs ∂x ∂y ∂x ∂y C =1       sΤ s sΤ s sΤ s sΤ nc  % % % % % % % %s  % sΤDs B % sΤ + I xx ∂B Ds ∂B + I yy ∂B Ds ∂B + I xy  ∂B Ds ∂B + ∂B Ds ∂B   = ∑ AB c c c ∂x ∂x ∂y ∂y ∂y ∂y ∂x   C =1   ∂x

where

B% s = [B% 1s

B% 2s L B% 8s ]

(52)

% s , ∂B% s ∂x , ∂B% s ∂y are given in Eqs. (44)-(46). It can be observed from Eq. (51) that the and B I I I material matrix D s related to the shear term has been replaced by the modified material matrix D s . The modified material matrix D s is computed based on the stabilization technique proposed in Ref. [65], and expressed as

D s = Ds

h2 h2 + α l 2

(53)

where l is the longest length of the edges of the element, α is a correction coefficient and generally given in the interval of 0.05 ≤ α ≤ 0.15 [9]. It is found from numerical experiments of isotropic plates and laminated composite Reissner-Mindlin plates that the present method can provide reasonable accuracy in solutions.

% is rewritten as Finally, the smoothed stiffness matrix K IJ % =K % b +K %s K IJ IJ IJ

(54)

The computation of the force vector and mass matrix is not discussed in details in this paper, which can also be computed over the smoothing cells as those of Refs. [59] and [66], respectively.

4. Numerical examples In this section, several numerical examples are presented to demonstrate the performance of the proposed approach for static and free vibration analysis of isotropic plates and laminated composite Reissner-Mindlin plates. 4.1 Numerical examples for the analysis of isotropic plates 4.1.1 Shear locking and convergence test As shown in Figs. 3(a) and 3(b), the clamped and simply supported square plates subjected to transverse uniform load P = 1 are first analyzed, the length of the square plate is L = 80 and length to thickness

ratio is L h . The material properties of the isotropic plates are taken as: Young’s modulus E = 3.0 × 106 , Poisson ratio v = 0.3 . Firstly, the shear locking phenomenon is investigated with the L h ranging from

101 to 107 , and the upper right quadrant of the plate is discretized with 16 × 16 number of Q8 plate elements. The central deflection is given as w = wEh3 /12(1 − v 2 ) PL4 . In this numerical example, the influence of α within the modified material matrix D s of Eq. (53) is also investigated.

(a) Clamped square plate

NEL = 1

(b) Simply supported square plate

NEL = 4

NEL = 4

NEL = 1

NEL = 16

NEL = 16

(c) Discretization of upper right quadrant of the plate with regular and distorted Q8 elements ( NEL : Number of elements) Fig. 3. Clamped and simply supported square plates and its discretization with Q8 elements. -3

4.3

x 10

1.7

x 10

-3

Present method (α=0.05) Present method (α=0.10) Present method (α=0.15) QSN QSR Exact thin plate solution

1.6

4.2 Central deflection

Central deflection

1.5 4.1 4 Present method (α=0.05) Present method (α=0.10) Present method (α=0.15) QSN QSR Exact thin plate solution

3.9 3.8 3.7 1 10

10

2

10

3

4

10 L/h

5

10

1.4 1.3 1.2 1.1 1

7

10

0.9 1 10

10

2

10

3

4

10 L/h

10

5

10

7

(a) Simply supported square plate (b) Clamped square plate Fig. 4. Shear locking test for simply-supported and clamped square plates. QSN [33]: Q8 plate element with normal integration scheme. QSR [33]: Q8 plate element with reduced integration scheme. Figs. 4(a) and 4(b) show the shear locking test results with simply supported and clamped boundary conditions, respectively. Compared with the results of Ref. [33], it can be observed from Fig. 4 that the

present method yields a locking-free solution and agrees well with the exact thin plate solution. However, a significant shear locking phenomenon is observed with the isoparametric Q8 plate elements with normal integration (QSN) and reduced integration (QSR). In addition, it can be seen that the correction coefficient

α has little effect on the results of the present method. By careful comparison, we can see that α = 0.15 can provide more accurate solutions, thus α = 0.15 is employed for the rest numerical examples unless other specified. The convergence of the present method with regular and distorted meshes are also studied. Fig. 3(c) illustrates the discretization of upper right quadrant of the plate with regular and distorted Q8 plate elements, in which NEL=1, 4,16，64 are used respectively (note that only the discretization of

NEL=1，4， 16 are present in Fig. 3(c)). The distorted elements can be constructed by the following equations

 x′ = x + rcγ∆x   y′ = y + rcγ∆y

(55)

where x , y are the initial coordinates of the internal nodes of regular meshes, x′ , y′ are the new coordinates of node of the distorted meshes. rc is a random number between -1~1, and γ is a prescribed distortion ratio whose value is chosen between 0~0.5. ∆x and ∆y are initial nodal spacing in x- and y-directions of the initial meshes, respectively. The convergence of the present method with simply supported and clamped boundary conditions are given in Fig. 5. For the purpose of comparison, the results obtain from QSR [31] and NC8-DS [34] are also given in Fig. 5. It is concluded that the present method shows excellent convergence performance, especially, the mesh distortion has little effect on the convergence performance.

1.05

1

Normalized central deflection

Normalized central deflection

1.05

0.95 0.9 Present method (γ=0.0) Present method (γ=0.3) QSR NC8-DS Exact solution

0.85 0.8 0.75 1

4

16 NEL

64

1

0.95

0.9 Present method (γ=0.0) Present method (γ=0.3) QSR NC8-DS Exact solution

0.85

0.8 1

4

16 NEL

(a) Simply supported square plate (b) Clamped square plate Fig. 5. Results of the convergence test for square plates with different method ( L h = 100 ).

64

4.1.2 Sensitivity test to mesh distortion The effect of mesh distortion on computational accuracy is investigated. A circular plate subjected to a transverse uniform load ( P = 1 ) with clamped boundary condition is analyzed in this numerical example. The geometry and material properties are given as: radius R = 5 , thickness to radius ratio h / R = 0.02 and 0.2, Young’s modulus E = 10.92 , Poisson’s ratio v = 0.3 . The value of irregularity factor γ is chosen as 0.0, 0.05, 0.1, 0.2 and 0.3, respectively. Fig. 6 shows the distribution of Q8 elements with

γ = 0.0 , 0.2, where γ = 0.0 indicates the initial meshes. The analytical solutions of the circular isotropic plates can be found in Ref. [67], the central deflection with h R = 0.02 , 0.2 are respectively given as wcanalytical =9783.48 and 11.5513. The normalized central deflections of the present method with various meshes are given in Table 1. It can be seen that the results of the present method agree well with the analytical solutions, the present method is quite insensitive to the mesh distortion.

(a) Initial meshes ( γ =0.0 ) (b) Distorted meshes ( γ =0.2 ) Fig. 6. Discretization of circular isotropic plate with Q8 elements. Table 1. Normalized central deflections of clamped circular plate subjected to transverse uniform load. Distortion ratio γ =0.0 γ =0.05 γ =0.1 γ =0.2 γ =0.3

Normalized central deflection ( wc / wcanalytical ) h R = 0.02 ( h = 0.1)

h R = 0.2 ( h = 1)

1.0007 1.0008 1.0011 1.0029 1.0051

1.0016 1.0018 1.0024 1.0044 1.0099

4.1.3 Free vibration analysis for L-shaped isotropic plate A clamped L-shaped isotropic plate with thickness to length ratio of h L = 0.1 and L = 2 is considered for the free vibration analysis. The regular and distorted meshes of the L-shaped isotropic plate with Q8 elements are shown in Fig. 7. The material properties are given as: Young’s modulus E = 10.92 , Poisson’s ratio v = 0.3 , density ρ = 1 . A non-dimensional natural frequency is defined by ω =ω L ρ G ,

G = E 2 (1 + v ) . The first four non-dimensional natural frequencies of the clamped L-shaped isotropic plate

are listed in Table 2. The present solutions are compared with those obtained from isogeometric approach with quartic elements and mixed interpolated tensorial component (MITC16) elements [9]. It can be seen that the present method yields reasonable solutions for both regular and distorted meshes. Moreover, the first six mode shapes of the L-shaped plate are plotted in Fig. 8.

(b) Distorted meshes （ γ =0.3 ）

(a) Regular meshes

Fig. 7. Clamped L-shaped isotropic plate and its discretization with Q8 elements. Table 2. Non-dimensional natural frequencies of clamped L-shaped isotropic plate. Non-dimensional natural frequency ( ω =ω L ρ G )

Mode

Regular meshes 1.8422 2.3752 2.7460 3.5862

1 2 3 4

Distorted meshes 1.8369 2.3672 2.7372 3.5743

IGA [9] 1.8147 2.3741 2.7570 3.6091

MITC16 [9] 1.8397 2.3712 2.7419 3.5827

0

0.02

0. 02

- 0. 02

0

0

- 0. 04 0

-0.02 0

0

0

1

1

1

1

2 2

2 2

(b) Mode 2

(c) Mode 3

0. 02

0. 02

0. 02

0

0

0

- 0. 02 0

- 0. 02 0

0 1 2 2

(d) Mode 4

0 1

1 2 2

1

1

2 2

(a) Mode 1

1

0

- 0. 02 0

0

- 0. 02 0

1

1 2 2

(e) Mode 5 (f) Mode 6 Fig. 8. First six mode shapes of the L-shaped isotropic plate.

4.2. Numerical examples for the analysis of laminated composite plates 4.2.1 Static analysis of square laminated composite plate under sinusoidal load A further application of the proposed method for the analysis of laminated composite plates is presented in this section. As shown in Fig. 9, we consider a simply supported square laminated composite plate

[0° / 90° / 90° / 0° ] subjected to a sinusoidal load q0 . The geometry and material properties are given as: L = 10 ,

L / h = 4,10, 20, 100 ,

E1 = 25 E2 ,

G12 = G13 = 0.5 E2 ,

G23 = 0.2 E2

and

v12 = 0.25 ,

respectively. The present solutions are compared with those obtained from isogeometric approach with quartic elements [9], FSDT [11] and high-order collocation method [5]. It can be observed from Table 3 that the present method also shows excellent performance for the analysis of both thick and thin square laminated composite plates. Meanwhile, it indicates that the present method can yield reasonable solutions even the distorted meshes are used.

q0

Fig. 9. Simply supported square laminated composite plate subjected to a sinusoidal load. Table 3. Non-dimensional central deflections of simply supported square laminated plate [0° / 90° / 90° / 0° ] under sinusoidal load. Non-dimensional central deflection w = (100 E2 h 3 ) w( L / 2, L / 2, 0) qL4 L h Number of elements Present method IGA [9] FSTD [11] Wavelets [5] Distorted ( γ =0.3 ) Regular 4
4 1.7023 1.7059 8
8 1.7089 1.7246 4 1.7164 1.71 1.7095 16
16 1.7095 1.7194 20
20 1.7095 1.7195

10

4
4 8
8 16
16 20
20

0.6641 0.6641 0.6632 0.6630

0.6683 0.6698 0.6698 0.6680

0.6654

0.6628

0.6627

20

4
4 8
8 16
16 20
20

0.4911 0.4925 0.4917 0.4915

0.4934 0.4959 0.4946 0.4956

0.4931

0.4912

0.4912

100

4
4 8
8 16
16 20
20

0.4047 0.4338 0.4337 0.4338

0.3973 0.4374 0.4362 0.4366

0.4354

0.4337

0.4335

4.2.2 Free vibration analysis of square laminated composite plate A symmetric cross-ply [0° / 90° / 0° ] square laminated composite plate is employed in the free vibration analysis. The geometry and material properties are given as: L = 10 , L / h = 100 , E2 = 7.6 × 109 ,

E1 E2 =40 , G12 = G13 = 0.6 E2 , G23 = 0.5 E2 and v12 = 0.25 , ρ = 1643 . Table 4. Non-dimensional fundamental frequency ( ωnum =ω = (ω L2 h ) ρ E2 ) of simply supported and clamped symmetric cross-ply [0° / 90° / 0° ] square composite plates with regular meshes. Boundary condition Number of elements Simply supported (SS) Clamped (CC) 4
4 19.1089 45.7528 6
6 18.9297 41.3924 8
8 18.8890 40.8685 10
10 18.8813 40.7511 Exact [68] 18.891 40.743 -0.5

log10 (|w num/w exact -1|)

-1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5

Simply supported square plate Convergence rate =-2.13 (SS) Clamped square plate Convergence rate =-3.35 (CC) 1.3

1.4

1.5

1.6

1.7

1.8

1.9

2

log10 (Number of element)

Fig. 10. Convergence rate for a fully simply supported and clamped laminated composite square plate

(a) γ =0.0 (b) γ =0.2 (c) γ =0.4 Fig. 11. Discretization schematic of a square laminated composite plate with Q8 elements. Firstly, the convergence of the present method for the free vibration analysis of square laminated composite plates with simply supported (SS) and clamped (CC) boundary conditions is investigated. The present results with regular meshes are given in Table 4 compared with the exact non-dimensional fundamental frequency ( ωexact ) obtained from Khdeir and Librescu [68]. It can be seen that the present

results with 10
10 elements agree very well with the exact solutions. In addition, the convergence of the non-dimensional fundamental frequency is also presented in Fig. 10, it can be seen that the good convergence rates are obtained by the present method. Secondly, the sensitivity to mesh distortion is also investigated with the meshes of 10
10 in this numerical example. Fig .11 only illustrates the discretization with regular ( γ =0.0 ) and distorted meshes ( γ =0.2，0.4 ). It can be seen from Fig. 12 that the present method is quite insensitive to the mesh distortion, even the extremely distorted meshes can yield reasonable solutions. 1.05 Simply supported square plate Clamped square plate Exact

wnum/wexact

1.03

1.01 1 0.99

0.97

0.95

0

0.1

0.2 0.3 Distortion ratio γ

0.4

0.5

Fig. 12. Effect of mesh distortion on the accuracy of the fundamental frequency with simply supported and clamped boundary conditions.

4.2.3 Free vibration analysis of skew laminated composite plate Five-layer symmetric cross-ply [90° / 0° / 90° / 0° / 90° ] and angle-ply [ 45° / −45° / 45° / −45° / 45° ] skew laminated plates with clamped boundary condition are considered in this section, and the skew angle ϕ is chosen as 0° , 15° , 30° , 45° and 60° , respectively. The geometry and material properties are given as:

L = 10 , L / h = 10 , E2 = 7.6 × 109 , E1 = 40 E2 , G12 = G13 = 0.5 E2 , G23 = 0.2 E2 , and v12 = 0.25 ,

ρ = 1643 . The schematic of discretization of the skew laminated composite plate with Q8 elements is shown in Fig. 13. The present numerical results with regular and distorted meshes are listed in Tables 5 and 6, respectively. For the purpose of comparison, the results obtained from other published paper are also listed in Tables 5 and 6. We can find that a good agreement between the present results and the reference results for both cases of cross-ply and angle-ply skew laminated composite plates can be observed from Tables 5 and 6. In addition, the first six mode shapes of cross-ply and angle-ply skew laminated composite °

plates ( ϕ = 15 ) are plotted in Figs. 14 and 15.

ϕ

(a) Mesh A (b) Mesh B ( γ =0.3 ) Fig. 13. Schematic of discretization of a skew laminated composite plate with Q8 elements. Table 5. Non-dimensional natural frequencies of clamped symmetric cross-ply [90° / 0° / 90° / 0° / 90° ] skew composite plates. Non-dimensional natural frequency ω = (ω L2 ρ E2 ) (π 2 h )

ϕ

Present method Mesh A Mesh B 2.3710 2.3619 2.4553 2.4525

°

0 15° 30° 45° 60°

IGA [9]

RBF [6]

MISQ20 [69]

MLSDQ [70]

2.3724 2.4648

2.4021 2.4932

2.3869 2.4803

2.3790 2.4725

2.7656

2.7613

2.7799

2.8005

2.7998

2.7927

3.4323 4.8802

3.4235 4.8540

3.4573 4.9282

3.4923 4.9541

3.4893 4.9989

3.4723 4.9430

Table 6. Non-dimensional natural frequencies of clamped symmetric angle-ply [ 45° / −45° / 45° / −45° / 45° ] skew composite plates. Non-dimensional natural frequency ω = (ω L2 ρ E2 ) (π 2 h )

ϕ 0° 15° 30° 45° 60°

Present method Mesh A Mesh B 2.2820 2.2795

IGA [9]

RBF[6]

MISQ20 [69]

MLSDQ [70]

2.2776

2.3324

2.2908

2.2787

2.3477 2.6578

2.3451 2.6533

2.3432 2.6527

2.3962 2.6981

2.3570 2.6708

2.3504 2.6636

3.3447

3.3359

3.3383

3.3747

3.3683

3.3594

4.8372

4.8137

4.8281

4.8548

4.8982

4.8566

0

0. 1

0. 1

- 0. 05

0

0

- 0. 1 10

20

5

10

- 0. 1 10

20

5

10

0 0

0 0

(a) Mode 1

(b) Mode 2

- 0. 1 10

20

5

10 0 0

(c) Mode 3

0. 1

0. 1

0. 1

0

0

0

- 0. 1 10

20

5

- 0. 1 10

20

5

10 0 0

- 0. 1 10

20

5

10 0 0

10 0 0

(d) Mode 4 (e) Mode 5 (f) Mode 6 Fig. 14. First six mode shapes of clamped symmetric cross-ply [90° / 0° / 90° / 0° / 90° ] skew composite plates ( ϕ = 15° ).

0

0. 1

0. 1

- 0. 05

0

0

- 0. 1 10

20

5

- 0. 1 10

20

5

10

- 0. 1 10

20

5

10

10

0 0

0 0

0 0

(a) Mode 1

(b) Mode 2

(c) Mode 3

0. 1

0. 1

0. 1

0

0

0

- 0. 1 10

20

5

- 0. 1 10

10 0 0

20

5

10 0 0

- 0. 1 10

20

5

10 0 0

(d) Mode 4 (e) Mode 5 (f) Mode 6 Fig. 15. First six mode shapes of clamped symmetric angle-ply [ 45° / −45° / 45° / −45° / 45° ] skew composite plates ( ϕ = 15° ).

4.2.4 Free vibration analysis of circular laminated composite plate In this numerical example, a clamped symmetric four-layer [θ / −θ / −θ / θ ] laminated composite plate with various ply angles is analyzed. The material properties are same as those used in Section 4.2.3, and the radius to thickness ratio is R h = 5, ( R = L 2 ) , the fiber orientation angle is given as θ = 15° , 30° , 45° , respectively. As shown in Fig. 16, two different types of meshes are used in this example. The effect of the ply angle θ on the non-dimensional natural frequency of the clamped laminated composite plate is presented in Table 7. Good agreement between the present results and the reference results for both meshes can also be observed from Table 7. The first six mode shapes of circular laminated composite plate

[ 45° / −45° / −45° / 45° ] are plotted in Fig. 17.

Mesh A Mesh B ( γ =0.3 ) Fig. 16. Schematic of discretization of circular laminated composite plate with Q8 elements. Table 7. Non-dimensional natural frequencies of circular four-layer [θ / −θ / −θ / θ ] clamped laminated composite plate. Non-dimensional natural frequency ω = (ω L2 h ) ρ E2 θ Mode Present method IGA [9] MISQ20 [69] MLSDQ [7] Mesh A Mesh B 22.7630 1 22.7314 22.6751 22.698 22.774 2 30.8107 31.3359 31.568 31.455 30.8596 3 43.2648 43.3550 43.635 43.350 43.3301 15° 4 43.6983 43.4165 44.318 43.469 43.7583 5 52.0427 52.9486 53.468 52.872 52.4372 6 59.1587 57.8349 60.012 57.386 59.2438

30°

1 2 3 4 5 6

23.6699 36.0653 43.0849 51.6651 56.1719 64.8238

23.6270 35.9966 43.0152 51.5722 56.0548 64.7233

23.9703 36.0298 43.8390 51.0024 56.7337 66.1316

24.046 36.399 44.189 52.028 57.478 67.099

24.071 36.153 43.968 51.074 56.315 66.220

45°

1 2 3 4 5 6

24.7233 39.1438 43.6284 57.2138 57.2302 65.5721

24.6770 39.0515 43.5588 57.0599 57.1275 65.4522

24.6622 38.9814 43.4559 56.9205 56.9844 65.3320

24.766 39.441 43.817 57.907 57.945 66.297

24.752 39.181 43.607 56.759 56.967 65.571

(a) Mode 1

(b) Mode 2

(c) Mode 3

(e) Mode 5 (f) Mode 6 (d) Mode 4 Fig. 17. First six mode shapes of clamped circular four-layer [ 45° / −45° / −45° / 45° ] laminated composite plate with Mesh A.

4.2.5 Free vibration analysis of square plate with a complicated cutout We adopt the present method to analyze a three-layer [θ / −θ / θ ] symmetric laminated composite square plate with complicated cutout, in which the fiber orientation angle is given as θ = 0° , 15° , 30° , 45° , respectively. The geometry and the discretization are shown in Fig. 18. The material and geometrical parameters are same as those used in Ref. [27] and given as: E1 = 2.45 E2 , G12 = G13 = 0.48 E2 ,

G23 = 0.2 E2 , v12 = 0.23 , ρ = 8000 , L = 10 and h = 0.06 . The non-dimensional natural frequency is 12

defined by ω = (ω 2 L4 ρ h D0.1 )

with D0.1 = E1h 3 12 (1 − v12 v21 ) . The non-dimensional first six frequencies

of the full simply supported three-ply laminated composite plate with various ply angles are presented in Table. 8. For the purpose of comparison, the solutions obtained by IGA with level set method [27], EFG and MKI methods [23] are also listed in Table 8. It can be clearly observed that the present solutions agree well with the solutions obtained from other published paper for different fiber orientations, and the present method can yield reasonable solutions even the distorted meshes are involved. In addition, the non-dimensional frequency parameters of the plate with clamped boundary condition are also listed in Table 9. It also shows that the present results obtained by regular and distorted meshes are in good agreement with the reference results. Furthermore, the first six mode shapes of the three-layer

[ 45° / −45° / 45° ] laminated composite square plate with complicated cutout under full simply supported and clamped boundary conditions are plotted in Figs. 19 and 20.

Mesh A Mesh B ( γ =0.3 ) Fig. 18. Schematic of discretization of circular laminated composite plate with Q8 elements. Table 8. Non-dimensional natural frequencies of fully simply supported three-ply [θ / −θ / θ ] laminate with complicated shaped cutout for various orientations. 12

0°

1 2 3 4 5 6

Non-dimensional natural frequency ω = (ω 2 L4 ρ h D0.1 ) Present method EFG [23] Mesh A Mesh B 18.3531 18.226 18.3780 31.1012 31.0681 31.127 36.0741 36.237 36.1096 56.7603 56.874 56.8200 61.9892 62.390 62.0291 83.3851 83.565 83.4181

15°

1 2 3 4 5 6

19.2960 32.3156 36.4212 57.9682 63.2451 85.2110

19.2655 32.2843 36.3693 57.9156 63.1915 85.1946

19.177 32.445 37.238 58.716 63.994 86.500

18.323 31.472 37.617 63.077 66.538 86.486

19.097 32.141 36.449 57.541 63.320 84.702

30°

1 2 3 4 5 6

20.8182 34.1579 37.5490 60.2124 65.4363 89.3077

20.7943 34.1229 37.4962 60.1444 65.3588 89.2360

20.926 34.915 39.101 62.222 67.054 92.715

20.310 33.987 39.898 58.111 69.699 92.099

20.603 33.987 37.600 59.762 65.646 88.731

45°

1 2 3 4 5 6

21.5324 34.9588 38.2188 61.3236 66.5572 92.4156

21.5028 34.9199 38.1953 61.2598 66.5390 92.3969

21.736 36.079 39.975 63.897 68.525 96.767

20.987 34.897 39.269 63.375 69.017 96.588

21.310 34.791 38.278 60.859 66.841 91.527

θ

Mode

MKI [23]

IGA [27]

18.169 30.303 36.581 57.429 64.145 85.656

18.190 30.928 36.073 56.390 61.986 82.894

(a) Mode 1

(b) Mode 2

(c) Mode 3

(d) Mode 4

(e) Mode 5

(f) Mode 6

Fig. 19. First six mode shapes of three-layer [ 45° / −45° / 45° ] laminated composite square plate with complicated cutout for full simply supported boundary condition (with Mesh A). Table 9. Non-dimensional natural frequencies of clamped three-ply [θ / −θ / θ ] laminate with complicated shaped cutout for various orientations. 12

θ

Mode

0°

1 2 3 4 5 6

Non-dimensional natural frequency ω = (ω 2 L4 ρ h D0.1 ) Present method IGA [27] Mesh A Mesh B 44.6274 44.239 44.6460 72.1580 71.348 72.1461 82.5299 81.520 82.5274 97.8056 96.969 97.8061 102.4001 101.449 102.3943 119.3124 119.3769 119.014

15°

1 2 3 4 5 6

45.4698 72.5055 82.1270 97.5403 102.6507 119.5709

45.4179 72.5278 82.0781 97.5559 102.6381 119.5840

45.052 71.720 81.148 96.742 101.823 119.294

30°

1 2 3 4 5 6

46.1429 74.1136 80.1125 97.2169 102.2393 121.6724

46.0894 74.0854 80.1038 97.2280 102.2638 121.7483

45.704 73.315 79.179 96.413 101.720 121.293

45°

1 2 3 4 5 6

46.4141 76.9313 76.9572 97.0625 101.9821 124.8695

46.4169 76.8906 77.0115 97.0688 102.0280 124.9527

45.965 75.993 76.171 96.262 101.674 124.967

(a) Mode 1

(b) Mode 2

(c) Mode 3

(d) Mode 4 (e) Mode 5 (f) Mode 6 ° ° ° Fig. 20. First six mode shapes of three-layer [ 45 / −45 / 45 ] laminated composite square plate with complicated cutout for clamped boundary condition (with Mesh A).

5. Conclusions The conventional isoparametric eight-node Reissner-Mindlin plate element using the selectively reduced integration still produces shear-locking phenomenon below a certain length to thickness ratio and also yields poor solutions and sub-optimal convergence rates with distorted meshes. In this paper, a linear smoothed quadratic eight-node finite element based on the first order shear deformation theory has been developed for the static and free vibration analysis of isotropic and laminated composite plates without isoparametric mapping. The modified strain matrices are computed by using the divergence theorem and the first order Taylor’s expansion over the smoothing cells. In addition, the first order Taylor’s expansion is also used for the construction of stiffness matrix. All the numerical examples indicate that the present method can effectively eliminate shear-locking and yield more reasonable results for the distorted meshes. The excellent performance of the distorted meshes also promotes the applicative competence of the present method for the problems with complicated domains.

Acknowledgments The supports from National Natural Science Foundation of China (11672105, 11372106) and Natural Science Foundation of Hunan Province (2016JJ1009) are gratefully acknowledged.

References [1] Reissner E, Stavsky Y. Bending and stretching of certain types of heterogeneous aeolotropic elastic plates. J Appl Mech 1961;28(3): 402-408. [2] Yang PC, Norris CH, Stavsky Y. Elastic wave propagation in heterogeneous plates. Int J Solids Struct 1966;2(4): 665-684. [3] Whitney JM, Leissa AW. Analysis of heterogeneous anisotropic plates. J Appl Mech 1969;36(2): 261. [4] Whitney JM. The effect of transverse shear deformation on the bending of laminated plates. J Compos Mater 1969;3(3): 534-547. [5] Ferreira A, Castro LM, Bertoluzza S. A high order collocation method for the static and vibration analysis of composite plates using a first-order theory. Compos Struct 2009;89(3): 424-432. [6] Ferreira A, Roque C, Jorge R. Free vibration analysis of symmetric laminated composite plates by FSDT and radial basis functions. Comput Methods Appl Mech Eng 2005;194(39): 4265-4278. [7] Liew K, Huang Y, Reddy J. Vibration analysis of symmetrically laminated plates based on FSDT using the moving least squares differential quadrature method. Comput Methods Appl Mech Eng 2003;192(19): 2203-2222. [8] Liew K, Huang Y, Reddy J. Moving least squares differential quadrature method and its application to the analysis of shear deformable plates. Int J Numer Methods Eng 2003;56(15): 2331-2351. [9] Thai CH, Nguyen-Xuan H, Nguyen-Thanh N, Le TH, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite Reissner-Mindlin plates using NURBS-based isogeometric approach. Int J Numer Methods Eng 2012;91(6): 571-603. [10] Liew KM. Solving the vibration of thick symmetric laminates by Reissner/Mindlin plate theory and the p-Ritz method. J Sound Vib 1996;198(198): 343-360. [11] Akhras G. Static and vibration analysis of anisotropic composite laminates by finite strip method. Int J Solids Struct 1993;30(22): 3129-3137. [12] Reddy J, Phan N. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher-order shear deformation theory. J Sound Vib 1985;98(2): 157-170. [13] Phan ND, Reddy JN. Analysis of laminated composite plates using a higher-order shear deformation theory. Int J Numer Methods Eng 1985;21(12): 2201-2219. [14] Reddy JN, Khdeir A. Buckling and vibration of laminated composite plates using various plate theories. Aiaa J 2012;27(12): 1808-1817. [15] Khdeir AA, Librescu L. Analysis of symmetric cross-ply laminated elastic plates using a higher-order theory: Part II—Buckling and free vibration. Compos Struct 1988;9(4): 259-277. [16] Panda SC, Natarajan R. Finite element analysis of laminated composite plates. Int J Numer Methods Eng 1979;14(1): 69-79. [17] Reddy J. Free vibration of antisymmetric, angle-ply laminated plates including transverse shear deformation by the finite element method. J Sound Vib 1979;66(4): 565-576. [18] Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Arch Comput Method Eng 2002;9(2): 87-140. [19] Sollero P, Aliabadi M. Fracture mechanics analysis of anisotropic plates by the boundary element method. Int J Fract 1993;64(4): 269-284. [20] Paiva WP, Sollero P, Albuquerque EL. Modal analysis of anisotropic plates using the boundary element method. Eng Anal Bound Elemen 2011;35(12): 1248-1255. [21] Liu GR, Chen X. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. J Sound Vib 2001;241(5): 839-855. [22] Dai K, Liu G, Lim K, Chen X. A mesh-free method for static and free vibration analysis of shear deformable laminated composite plates. J Sound Vib 2004;269(3): 633-652. [23] Bui TQ, Nguyen MN, Zhang C. An efficient meshfree method for vibration analysis of laminated composite plates. Comput Mech 2011;48(2): 175-193. [24] Liu GR, Zhao X, Dai KY, Zhong ZH, Li GY, Han X. Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method. Compos Sci Technol 2008;68(2): 354-366. [25] Thai CH, Ferreira A, Carrera E, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a layerwise deformation theory. Compos Struct 2013;104: 196-214. [26] Yin S, Yu T, Bui TQ, Xia S, Hirose S. A cutout isogeometric analysis for thin laminated composite plates using level sets. Compos Struct 2015;127: 152-164. [27] Yu T, Yin S, Bui TQ, Xia S, Tanaka S, Hirose S. NURBS-based isogeometric analysis of buckling and free vibration problems for laminated composites plates with complicated cutouts using a new simple FSDT theory and level set method. Thin-Walled Struct 2016;101: 141-156.

[28] Thai CH, Nguyen-Xuan H, Bordas SPA, Nguyen-Thanh N, Rabczuk T. Isogeometric analysis of laminated composite plates using the higher-order shear deformation theory. Mech Adv Mater Struct 2015;22(6): 451-469. [29] Thai CH, Ferreira A, Bordas SPA, Rabczuk T, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. Eur J Mech A-Solids 2014;43: 89-108. [30] Bathe KJ, Dvorkin EN. A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. Int J Numer Methods Eng 1985;21(2): 367-383. [31] Dvorkin EN, Bathe KJ. A continuum mechanics based four‐node shell element for general non‐linear analysis. Eng Comput 1984;1(1): 77-88. [32] Zienkiewicz OC, Taylor RL, Too JM. Reduced integration technique in general analysis of plates and shells. Int J Numer Methods Eng 1971;3(2): 275-290. [33] Pugh EDL, Hinton E, Zienkiewicz OC. A study of quadrilateral plate bending elements with ‘reduced’ integration. Int J Numer Methods Eng 1978;12(7): 1059-1079. [34] Choi CK, Kim SH. Coupled use of reduced integration and non‐conforming modes in quadratic Mindlin plate element. Int J Numer Methods Eng 1989;28(8): 1909-1928. [35] Kim SH, Choi CK. Improvement of quadratic finite element for Mindlin plate bending. Int J Numer Methods Eng 1992;34(1): 197-208. [36] Katili I. A new discrete Kirchhoff-Mindlin element based on Mindlin‐Reissner plate theory and assumed shear strain fields-part I: An extended DKT element for thick-plate bending analysis. Int J Numer Methods Eng 1993;36(11): 1859-1883. [37] Katili I. A new discrete Kirchhoff-Mindlin element based on Mindlin-Reissner plate theory and assumed shear strain fields-part II: An extended DKQ element for thick-plate bending analysis. Int J Numer Methods Eng 1993;36(11): 1885-1908. [38] Polit O, Touratier M, Lory P. A new eight-node quadrilateral shear-bending plate finite element. Int J Numer Methods Eng 1994;37(3): 387-411. [39] Choi CK, Park YM. Quadratic NMS Mindlin-plate-bending element. Int J Numer Methods Eng 1999;46(8): 1273-1289. [40] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Comput Mech 2007;39(6): 859-877. [41] Chen JS, Yoon S, Wu CT. A stabilized confirming nodal integration for Galerkin meshfree methods. Int J Numer Methods Eng 2002;53(12): 2587-2615. [42] Liu GR, Nguyen-Thoi T, Nguyen-Xuan H, Lam KY. A node-based smoothed finite element method (NS-FEM) for upper bound solutions to solid mechanics problems. Comput Struct 2009;87(1-2): 14-26. [43] Liu GR, Nguyen-Thoi T, Lam KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. J Sound Vib 2009;320(4-5): 1100-1130. [44] Nguyen-Thoi T, Liu GR, Lam KY, Zhang GY. A face-based smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements. Int J Numer Methods Eng 2009;78(3): 324-353. [45] Thai CH, Tran LV, Tran DT, Nguyen-Thoi T, Nguyen-Xuan H. Analysis of laminated composite plates using higher-order shear deformation plate theory and node-based smoothed discrete shear gap method. Appl Math Model 2012;36(11): 5657-5677. [46] Phan-Dao H, Nguyen-Xuan H, Thai-Hoang C, Nguyen-Thoi T, Rabczuk T. An edge-based smoothed finite element method for analysis of laminated composite plates. Int J Comput Methods 2013;10(01): 1340005. [47] Cui XY, Liu GR, Li GY, Zhao X, Nguyen-Thoi T, Sun G. A smoothed finite element method (SFEM) for linear and geometrically nonlinear analysis of plates and shells. Comp Model Eng 2008;28(2): 109-125. [48] Nguyen-Van H, Mai-Duy N, Karunasena W, Tran-Cong T. Buckling and vibration analysis of laminated composite plate/shell structures via a smoothed quadrilateral flat shell element with in-plane rotations. Comput Struct 2011;89(7-8): 612-625. [49] Phung-Van P, Thai CH, Nguyen-Thoi T, Nguyen-Xuan H. Static and free vibration analyses of composite and sandwich plates by an edge-based smoothed discrete shear gap method (ES-DSG3) using triangular elements based on layerwise theory. Compos Pt B-Eng 2013;60(4): 227-238. [50] Nguyen-Thoi T, Phung-Van P, Nguyen-Xuan H, Thai-Hoang C. A cell-based smoothed discrete shear gap method using triangular elements for static and free vibration analyses of ReissnerMindlin plates. Int J Numer Methods Eng 2012;91(7): 705 -741. [51] Phung-Van P, Nguyen-Thoi T, Dang-Trung H, Nguyen-Minh N. A cell-based smoothed discrete shear gap method (CS-FEM-DSG3) using layerwise theory based on the C0-HSDT for analyses of composite plates. Compos Struct 2014;111(1): 553-565. [52] Nguyen-Thanh N, Rabczuk T, Nguyen-Xuan H, Bordas SP. A smoothed finite element method for shell analysis. Comput Methods Appl Mech Eng 2008;198(2): 165-177.

[53] Jiang C, Zhang ZQ, Han X, Liu GR. Selective smoothed finite element methods for extremely large deformation of anisotropic incompressible bio-tissues. Int J Numer Methods Eng 2014;99(8): 587-610. [54] Jiang C, Liu GR, Han X, Zhang ZQ, Zeng W. A smoothed finite element method for analysis of anisotropic large deformation of passive rabbit ventricles in diastole. Int J Numer Meth Biomed 2015;31(1): 1-25. [55] Bordas SPA, Natarajan S, Kerfriden P, Augarde CE, Mahapatra DR, Rabczuk T, Pont SD. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). Int J Numer Methods Eng 2011;86(4-5): 637-666. [56] Francis A, Ortiz-Bernardin A, Bordas SPA, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. Int J Numer Methods Eng 2016;109: 1263-1288. [57] Duan Q, Gao X, Wang B, Li X, Zhang H, Belytschko T, Shao Y. Consistent element-free Galerkin method. Int J Numer Methods Eng 2014;99(2): 79-101. [58] Duan Q, Li X, Zhang H, Belytschko T. Second-order accurate derivatives and integration schemes for meshfree methods. Int J Numer Methods Eng 2012;92(4): 399-424. [59] Duan Q, Li X, Zhang H, Wang B, Gao X. Quadratically consistent one-point (QC1) quadrature for meshfree Galerkin methods. Comput Methods Appl Mech Eng 2012;s 245–246(3): 256-272. [60] Surendran M, Natarajan S, Bordas SPA, Palani GS. Linear smoothed extended finite element method. Int J Numer Methods Eng 2017; DOI: 10.1002/nme.5579 [61] Reddy JN. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. Crc Press 2004. [62] Wan D, Hu D, Natarajan S, Bordas SPA, Yang G. A fully smoothed XFEM for analysis of axisymmetric problems with weak discontinuities. Int J Numer Methods Eng 2016; 110(3):203-226. [63] Wan D, Hu D, Yang G, Long T. A fully smoothed finite element method for analysis of axisymmetric problems. Eng Anal Bound Elemen 2016;72: 78-88. [64] Dasgupta G. Integration within polygonal finite elements. J Aerospace Eng 2003;16(1): 9-18. [65] Lyly M, Stenberg R, Vihinen T. A stable bilinear element for the Reissner-Mindlin plate model. Comput Methods Appl Mech Eng 1993;110(3-4): 343-357. [66] Yang G, Hu D, Ma G, Wan D. A novel integration scheme for solution of consistent mass matrix in free and forced vibration analysis. Meccanica 2016;51(8): 1897-1911. [67] Rezak A, Alain R. An improved four-node hybrid-mixed element based upon Mindlin's plate theory. Int J Numer Methods Eng 2002;55(6): 705-731. [68] Khdeir A, Librescu L. Analysis of symmetric cross-ply laminated elastic plates using a higher-order theory: Part II-Buckling and free vibration. Compos Struct 1988;9(4): 259-277. [69] Nguyen-Van H, Mai-Duy N, Tran-Cong T. Free vibration analysis of laminated plate/shell structures based on FSDT with a stabilized nodal-integrated quadrilateral element. J Sound Vib 2008;313(1-2): 205-223. [70] Liew KM, Huang YQ, Reddy JN. Vibration analysis of symmetrically laminated plates based on FSDT using the moving least squares differential quadrature method. Comput Methods Appl Mech Eng 2003;192(19): 2203-2222.