Spectral element modeling and analysis of the transverse vibration of a laminated composite plate

Accepted Manuscript Spectral element modeling and analysis of the transverse vibration of a laminated composite plate Ilwook Park, Usik Lee PII: DOI: Reference:

S0263-8223(15)00807-7 http://dx.doi.org/10.1016/j.compstruct.2015.08.111 COST 6816

To appear in:

Composite Structures

Please cite this article as: Park, I., Lee, U., Spectral element modeling and analysis of the transverse vibration of a laminated composite plate, Composite Structures (2015), doi: http://dx.doi.org/10.1016/j.compstruct.2015.08.111

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.

Spectral element modeling and analysis of the transverse vibration of a laminated composite plate Ilwook Park and Usik Lee∗ Department of Mechanical Engineering, Inha University 100 Inha-ro, Nam-gu, Incheon 402-751, Republic of Korea

ABSTRACT This paper presents a frequency-domain spectral element model for the symmetric laminated composite plates which have finite dimensions in two orthogonal directions, i.e., in the x- and y-directions. The spectral element model is developed by using two methods in combination: the splitting of original boundary conditions and the socalled super spectral element method in which both the Kantorovich method-based finite strip element method and the frequency-domain waveguide method are utilized. The present spectral element model has nodes (or degrees of freedom (DOFs)) only on four edges of a finite element, i.e., no nodes inside the finite element. Accordingly the total number of DOFs used in the dynamic analysis can be drastically decreased to lead to a significant decrease of the computation cost, when compared with the standard finite element method (FEM). The high accuracy of the present spectral element model is verified in due course by the comparison with the results by two solution methods: the exact theory available in the literature and the standard finite element model developed in this study.

Keywords: Laminated composite plates; Spectral element method; Finite element method; Vibration; Waves

∗

Corresponding author: Tel.: +82 32 860 7318; fax: +82 32 866 1434 E-mail address: [email protected] (U. Lee).

1. Introduction It has been well-recognized that composite materials have many advantages over metallic materials due to their high strength-to-density ratios. Thus composite materials have been increasingly applied in the various fields of engineering including aerospace, mechanical engineering, and civil engineering over the past decades [1]. However, the laminated composite structures are susceptible to various modes of damage including delamination. Therefore, it is very important to accurately estimate the dynamic characteristics of the laminated composite structures during the early design phase. Although the subject of the dynamic analysis of laminated composite plates (simply composite plates) has a quite long history, most existing analytical solutions to this classical problem are considered to be at best approximate because it is difficult to obtain exact closed-form solutions which simultaneously satisfy the governing differential equations and the associated boundary conditions, except for very specific types of plate such as the Levy-type plates [2]. The finite element method (FEM) is certainly one of the most powerful computational methods to analyze the dynamics of a structure with very complex geometry, material distribution, and boundary conditions. The accuracy of FEM depends on the size of finite elements (or meshes) used in the analysis because the displacement fields in a finite element are normally represented by simple polynomial functions which are not dependent of vibrating frequency. Accordingly, as a drawback of FEM, very fine meshing is required to improve the solution accuracy especially at high frequency: this may result in a significant increase of computation cost. Thus, the spectral element method (SEM) can be considered as an alternative to FEM because it can provide extremely accurate solutions even at very high frequencies by using the so-called spectral element matrix (or dynamic stiffness matrix) formulated from frequencydependent (dynamic) shape functions. In the literature, there are two completely different methods which have been called the same name ‘SEM’. The first one is the frequency-domain SEM [3, 4], where accurately formulated dynamic stiffness matrix is used as the finite element stiffness matrix. The second one is the time-domain SEM proposed by Patera [5]. The Patera’s

time-domain SEM is of course formulated in the time domain by using the Legendre or Chebyshev orthogonal polynomials as the shape functions in conjunction with the use of the Gauss-Lobatto-Legendre integration rule. Thus, the word ‘spectral’ is timewise for the frequency-domain SEM, while it is space-wise for the time-domain SEM. The SEM considered in this study is the frequency-domain SEM and it will be called ‘SEM’ throughout the paper. As mentioned previously, the spectral element matrix (or dynamic stiffness matrix) for SEM is formulated from free wave solutions that satisfy the governing differential equations transformed in the frequency domain [4]. Thus the SEM can provide extremely accurate solutions by representing a whole uniform structure member as a single finite element, regardless of its size, and it has been often recognized as an exact element method. The assembly of finite elements to a global structure system can be conducted in an exactly analogous way as that used in the standard FEM. Despite of outstanding features of SEM, however its applications have been limited mostly to one-dimensional (1D) structures (e.g., [3,4]) or the plates with very specific geometries and boundary conditions (e.g., [6-18]). In the literature, there are very few spectral element models for two-dimensional (2D) structures which have finite dimensions in both x and y directions and being subjected to arbitrary boundary conditions. This is because, for such finite 2D structures, it is not easy to obtain exact free wave solutions in analytical forms which are essentially required to formulate spectral element matrix. Only a few researchers [19-20] have presented approximate dynamic stiffness methods or spectral element methods for the in-plane or transverse vibrations of the finite plates. Recently, Park et al. [21] presented a spectral element model for a rectangular membrane by using two methods in combination to obtain the displacement field in a finite membrane element: (1) the splitting of boundary conditions; (2) the super spectral element method (SSEM) in which the Kantorovich method-based finite strip element method and the frequency-domain waveguide method are utilized. As an extension of the author’s previous work [21], this paper proposes a spectral element model for a symmetric composite plate. The proposed spectral element model is formulated by using the splitting of boundary conditions and the SSEM in combination. Though the spectral element modeling method newly proposed in this paper can be applied to various laminated plate theories including the classical laminated plate

theory (CLPT) and the improved theory such as the first-order shear deformation plate theory (FSDT), the discussion in this paper is limited to the CLPT. The performance of the proposed spectral element model is then evaluated in due course by the comparison with exact solutions and the solutions by the standard finite element model developed in this study.

2. Governing equations 2.1. Governing equations in the time domain Consider a symmetric composite plate whose layup is symmetric about the midplane of the plate. The composite plate is rectangular and its dimensions in the xand y-directions are lx and ly, respectively. The origin at coordinates (x, y) is placed at the middle of the composite plate, as shown in Fig. 1. The governing differential equation of motion for the transverse vibration of the symmetric composite plate is given by [1]

D11

∂4w ∂ 4w ∂4w ∂4w ∂4w + 2 ( D + 2 D ) + 4 D + 4 D + D 12 66 16 26 22 ∂x 4 ∂x 2∂y 2 ∂x3∂y ∂x∂y 3 ∂y 4

∂2w + ρ 2 = f ( x, y, t ) ∂t

(1)

where w(x,y,t) is the transverse displacement, f(x,y,t) is the external force applied normal to the surface of the plate, ρ is the mass per unit area of the plate, and Dij are the bending stiffnesses defined by L

Dij =

1 3

∑Q

(k ) ij

( zk3+1 − zk3 )

(i, j = 1, 2, 6,

k = 1, 2,
, L)

(2)

k =1

where Qij( k ) are the transformed plane stress-reduced stiffness coefficients of the kth orthotropic layer, (zk, zk+1) denote the thickness coordinates of the bottom and top, respectively, of the kth layer, and L denotes the total number of layers. The boundary conditions associated with Eq. (1) are given by

M y ( x,−l y / 2, t ) = − M y1 ( x, t )

or

θ y ( x,−l y / 2, t ) = θ y1 ( x, t )

Vy ( x,−l y / 2, t ) = −Vy1 ( x, t )

or

w( x,−l y / 2, t ) = w1 ( x, t )

M y ( x, l y / 2, t ) = M y 2 ( x, t )

or

θ y ( x, l y / 2, t ) = θ y 2 ( x, t )

Vy ( x, l y / 2, t ) = V y 2 ( x, t )

or

w( x, l y / 2, t ) = w2 ( x, t )

M x (−l x / 2, y, t ) = − M x1 ( x, t )

or

θ x (−lx / 2, y, t ) = θ x1 ( y, t )

Vx (−l x / 2, y, t ) = −Vx1 ( x, t )

or

w(−l x / 2, y, t ) = w3 ( y, t )

M x (l x / 2, y, t ) = M x 2 ( x, t )

or

θ x (l x / 2, y, t ) = θ x 2 ( y, t )

Vx (l x / 2, y, t ) = Vx 2 ( x, t )

or

w(l / 2 x , y, t ) = w4 ( y, t )

(3)

where (Mx, My) are the resultant bending moments and (Vx, Vy ) are the resultant transverse shear forces defined by

M x ( x, y, t ) = −D11

∂2w ∂2w ∂2w − D12 2 − 2 D16 , 2 ∂x ∂y ∂x∂y

M y ( x, y, t ) = − D12

∂2w ∂2w ∂2w − D − 2 D 22 26 ∂x 2 ∂y 2 ∂x∂y

Vx ( x, y, t ) = −

∂ ∂2w ∂2w ∂2w   D11 2 + D12 2 + 2 D16  ∂x  ∂x ∂y ∂x∂y 

∂ ∂2w ∂2w ∂2w   −  D16 2 + D26 2 + 2 D66 ∂y  ∂x ∂y ∂x∂y  Vy ( x, y, t ) = −

∂ ∂2w ∂2w ∂2w   D12 2 + D22 2 + 2 D26  ∂y  ∂x ∂y ∂x∂y 

−

∂ ∂2w ∂ 2w ∂2w   D16 2 + D26 2 + 2D66  ∂x  ∂x ∂y ∂x∂y 

(4)

and θx and θy are the slopes defined by

θx =

∂w ∂w , θy = ∂x ∂y

(5)

2.2. Governing equations in the frequency domain

To formulate spectral element model for a composite plate by following the general procedure introduced in Ref. [4], firstly all the time domain quantities in the governing differential equations of motion (Eq. (1)) and the boundary conditions (Eq. (3)) are transformed into the frequency domain quantities by using the discrete Fourier

transform (DFT) theory [22]. For instance, the transverse displacement w(x,y,t) and the external force f(x,y,t) can be represented in the spectral forms as { w( x, y, t ), f ( x, y , t )} =

1 N

N −1

∑{ w ( x, y), n

f n ( x, y )} eiωnt

(6)

n=0

where i = − 1 is the imaginary unit, N is the number of samples for the fast Fourier

transforms (FFT), and ωn = 2πn/T are the discrete frequencies up to the Nyquist frequency. The over-barred quantities wn ( x, y ) and f n ( x, y ) represent spectral or Fourier

components of the corresponding time domain quantifies w(x,y,t) and f(x,y,t), respectively. In the following derivations, the over-bars and the subscripts n will be omitted for the sake of brevity. By using the aforementioned DFT theory, the governing differential equation of motion (Eq. (1)) can be transformed into the frequency domain as

D11

∂ 4w ∂ 4w ∂4w ∂4w ∂ 4w + 2 ( D + 2 D ) + 4 D + 4 D + D 12 66 16 26 22 ∂x 4 ∂x 2∂y 2 ∂x 3∂y ∂x∂y 3 ∂y 4

(7)

2

− ρω w = f ( x, y ) Similarly, the boundary conditions (Eq. (3)) can be also transformed into the frequency domain as M y ( x,−l y / 2) = − M y1 ( x)

or

θ y ( x,−l y / 2) = θ y1 ( x)

Vy ( x,−l y / 2) = −Vy1 ( x)

or

w( x,−l y / 2) = w1 ( x)

M y ( x, l y / 2 ) = M y 2 ( x )

or

θ y ( x, l y / 2) = θ y 2 ( x)

Vy ( x, l y / 2) = V y 2 ( x)

or

w( x, l y / 2) = w2 ( x)

M x (−l x / 2, y ) = − M x1 ( x)

or

θ x (−l x / 2, y) = θ x1 ( y )

Vx (−l x / 2, y ) = −Vx1 ( x)

or

w(−l x / 2, y ) = w3 ( y)

M x (l x / 2, y) = M x 2 ( x)

or

θ x (l x / 2, y ) = θ x 2 ( y )

Vx (l x / 2, y ) = Vx 2 ( x)

or

w(l / 2 x , y ) = w4 ( y )

(8)

3. Free wave solution in the frequency domain To formulate the spectral element model for a finite composite plate element, we need to obtain the frequency domain free wave solution that satisfies the homogeneous

governing differential equation of motion that can be reduced from Eq. (7), by removing external force f(x, y), as follows:

D11

∂ 4w ∂ 4w ∂ 4w ∂ 4w ∂ 4w + 2 ( D + 2 D ) + 4 D + 4 D + D 12 66 16 26 22 ∂x 4 ∂x 2∂y 2 ∂x 3∂y ∂x∂y 3 ∂y 4

(9)

2

− ρω w = 0

The weak form of Eq. (9) can be obtained in the form as  ∂ 2w ∂2w ∂ 2w   ∂ 2w   δ   D + D + 2 D 16 ∫x∫y   11 ∂x2 12 ∂y 2 ∂x∂y   ∂x 2   ∂2w ∂2w ∂2 w   ∂ 2w  δ   +  D12 2 + D22 2 + 2 D26 ∂x ∂y ∂x∂y   ∂y 2  

(10)

  ∂2w ∂2w ∂ 2w   ∂ 2w   δ   − ω 2 ρwδw dxdy = 0 + 2  D16 2 + D26 2 + 2 D66 ∂x ∂y ∂x∂y   ∂x∂y    In this paper, the free wave solution will be obtained approximately so as to satisfy the weak form given by Eq. (10) by using two methods in combination: (1) the splitting of boundary conditions, and (2) the SSEM in which the Kantorovich method-based finite strip element method and the frequency-domain waveguide method are used. The concept of the splitting of boundary conditions is illustrated in Fig. 2. The original problem shown in Fig. 2(a) is represented by the sum of two partial problems, Problem A and Problem B. The geometric boundary conditions for the original problem, Problem A, and Problem B are presented in the simplified forms by using the symbols defined by w x ( x, y) = {w( x,y), θ x ( x,y) = ∂w/∂x}T w y ( x, y ) = {w( x,y), θ y ( x,y) = ∂w/∂y}T w xA ( x, y ) = {wA ( x,y), θ xA ( x,y) = ∂w A /∂x}T w yA ( x, y ) = {wA ( x,y), θ yA ( x,y) = ∂wA /∂y}T

(11)

w xB ( x, y ) = {wB ( x,y), θ xB ( x,y) = ∂wB /∂x}T w yB ( x, y) = {wB ( x,y), θ yB ( x,y) = ∂wB /∂y}T Notice that Problem A has the fixed (null) boundary conditions on two parallel

edges at y = – ly/2 and ly/2 and its free wave solution is represented by wA(x,y), while Problem B has the fixed boundary conditions on two parallel edges at x = – lx/2 and lx/2 and its free wave solution is represented by wB(x,y). Accordingly, once wA(x,y), and wB(x,y) are obtained for Problem A and Problem B, respectively, the free wave solution w(x, y) for the original problem can be readily obtained from w( x, y) = wA ( x, y ) + wB ( x, y )

(12)

3.1. Derivation of wA (x, y) The free wave solution wA(x,y) for Problem A will be obtained by using the SSEM. To that end, as shown in Fig. 3(a), a finite composite plate element that has the dimensions lx and ly in the x- and y-directions, respectively, is divided into Ny finite strip elements in the y–direction. The eth finite strip element shown in Fig. 3(b) has the width of l y( e ) = ye − ye−1 in the y–direction. The displacement field in the eth finite strip element can be represented by w(Ae ) ( x, y) = Y A( e ) ( y )u(Ae ) ( x)

( y e −1 ≤ y ≤ y e )

(13)

where Y A(e ) ( y ) is the one-by-four interpolation function matrix and u(Ae ) ( x) are the nodal line degrees of freedom (DOFs). They are defined by YA( e ) ( y ) = [ YA(1e ) ( y ), YA( e2) ( y) ] w (yAe −1) ( x) u ( x) =  ( e )   w yA ( x)  (e) A

(14)

where

YA(1e ) ( y ) = [ l y( e )−3 ( y − ye ) 2 (2 y + ye − 3 ye−1 ), l y( e )−2 ( y − ye )2 ( y − ye−1 )] YA( e2) ( y ) = [ l y( e )−3 ( y − ye−1 ) 2 (2 y + ye−1 − 3 ye ), l y( e )−2 ( y − ye−1 ) 2 ( y − ye )]

(15)

and (e) w (yAe ) ( x) = w yA ( x, ye ) = {w(Ae ) ( x, ye ), θ yA ( x, ye )}T (e) w (yAe−1) ( x) = w yA ( x, ye−1 ) = {w(Ae ) ( x, ye−1 ), θ yA ( x, ye−1 )}T

(16)

By using Eq. (13) and applying the fixed boundary conditions at y = - ly/2 and y = ly/2, the displacement field wA(x,y) in the whole area of the finite composite plate element can be expressed in the following form: w A ( x, y ) = Y A ( y ) u A ( x )

(−l y / 2 ≤ y ≤ l y / 2)

(17)

where ( N −2 )

YA ( y ) = [ L(A1) ( y ), L(A2) ( y ),
, L(Ae ) ( y ),
, LA y

( N −2 )

u A ( x) = { w (yA1) ( x), w (yA2) ( x),
, w (yAe ) ( x),
, w yA y

( N −1)

( y ), LA y ( y)] ( N −1)

( x), w yA y ( x)}T

(18)

For Eq. (18), the following definitions are used: L(Ae ) ( y ) = hA( e ) ( y )YA( e2) ( y ) + hA( e+1) ( y )YA(1e+1) ( y ) (e = 1, 2, 3, … , Ny -1)

(19)

and h A( e ) ( y) = H ( y − ye −1 ) − H ( y − ye )

(20)

where H(y) is the Heaviside unit step function. By substituting Eq. (17) into the weak form given by Eq. (12), we can obtain a set of differential equations in the form as ∂ 4 uA ( x) ∂ 3uA ( x) ∂ 2uA ( x) ∂u ( x) + A + A + AA1 A A 3 A 2 4 3 2 ∂x ∂x ∂x ∂x 2 + ( AA0 − ω M A )uA ( x) = 0 AA4

(21)

where AA0 = D22 Λy 6 ,

AA1 = 2 D26 ( Λy 5 − ΛyT5 )

AA2 = D12 ( Λy 3 + ΛTy 3 ) − 4 D66 Λy 4 , AA3 = 2 D16 ( Λy 2 − ΛTy 2 ) AA4 = D11 Λy1 , with

M A = ρΛy1

(22)

ly / 2

Λy1 = ∫

−l y / 2

Λy 3 = ∫

ly / 2

−l y / 2 ly / 2

Λy 5 = ∫

−l y / 2

Y TYdy, YT

∂ 2Y dy, ∂y 2

Λy 2 = ∫

ly / 2

−l y / 2

Λy 4 = ∫

ly / 2

−l y / 2

YT

∂Y dy ∂y

∂Y T ∂Y dy ∂y ∂y

(23)

2 T ly / 2 ∂ Y ∂ 2Y T ∂Y ∂ 2Y dy , Λ = y6 ∫−ly / 2 ∂y 2 ∂y 2 dy ∂y 2 ∂y

The constant matrices Λyj (j = 1, 2, 3, …, 6) are provided in Appendix A. We can assume the solution to Eq. (21) in the form as  1   r   2  k x x − 12 k x l x k x− 1 k l u A ( x) =  ae = ra e x 2 x x     r2( N y −1)   

(24)

where a is constant, kx is the wavenumber in the x-direction, and + kx  k x = − k x  0 

if Re(k x ) > 0 if Re(k x ) < 0

(25)

if Re(k x ) = 0

By substituting Eq. (24) into Eq. (21), we can obtain an eigenvalue problem as follows: [ AA4 k x4 + AA3k x3 + AA2 k x2 + AA1k x + AA0 − ω 2 M A ] rA = 0

(26)

The dispersion relation or the frequency-wavenumber relationship can be derived from det [ AA4 k x4 + AA3k x3 + AA2 k x2 + AA1k x + AA0 − ω 2 M A ] = 0

(27)

By using the wavenumbers kx(j) (j = 1, 2, 3,…, 8×(Ny-1)) computed from Eq. (27), we can write the solution to Eq. (21) in the following form: u A ( x ) = R A E A ( x )a A

where

(28)

a A = { a1 a2 a3 … ai … a8×( N y −1) }T R A = [rA(1) , rA( 2) , rA(3) , …, rA( j ) , …, rA(8×( N y −1)) ] E A ( x; ω ) = diagonal [ e

k x ( j ) x − 12 k x ( j )lx

]

(29)

( j = 1, 2, 3,
, 8( N y − 1))

In Eq. (29), rA(j) represents the jth wave mode vector (or eigenvector) that can be readily computed from Eq. (26) by putting kx = kx(j). The nodal DOFs at the nodes defined on the edges at x = – lx/2 and lx/2 can be written in the vector form as  u A (− l2x )   ∂u (− lx )   A 2   d A(1)   ∂x  dA =  =   d A( 2)   u A (+ l2x )   ∂u A (+ l2x )     ∂x 

(30)

where  w (A1()1)   ( 2)   w A(1)     u A (− l2x ) =  ( e ) ,  w A(1)      ( N y −1)   w A(1)   w A(1()2)   ( 2)   w A( 2 )     u A (+ l2x ) =  (e ) ,  w A( 2 )      ( N y −1)  w A( 2)  and where

(1)  θ xA  (1)  ( 2)   θ xA(1)  lx ∂u A (− 2 )    =  (e )  ∂x  θ xA(1)      ( N y −1)  θ xA(1)  (1)  θ xA  (2)  ( 2)   θ xA( 2 )  ∂u A (+ l2x )    =  (e)  ∂x  θ xA( 2 )      ( N y −1)  θ xA( 2 ) 

(31)

w (e ) (∓ lx , y ) w (Ae()1, 2) = w (yAe ) (∓ l2x ) =  (Ae ) l2x e  θ yA (∓ 2 , ye )  (e) ∂w (yAe ) (∓ l2x ) θ xA (∓ l2x , ye )  (e) θ xA(1, 2) = =  ( e ) lx  ∂x θ xyA (∓ 2 , ye )

(32)

By substituting Eq. (28) into Eq. (30), the nodal DOFs vector d A can be written in terms of constant vector aA as follows:

d A = H A (ω )a A

(33)

where  R A E A (− l2x )    R K E ( − lx ) H A (ω ) =  A A A lx 2   RA E A ( + 2 )   l   R A K A E A (+ 2x )

(34)

The constant vector aA can be removed from Eq. (28) by using Eq. (33) to obtain the expression as u A ( x ) = X A ( x ; ω )d A

(35)

where X A ( x; ω) = R A E A ( x) H A (ω) −1

(36)

Finally, by substituting Eq. (35) into Eq. (17), we can finally obtain the free wave solution wA(x,y) in the form as follows: wA ( x, y) = N A ( x, y; ω)d A

(37)

N A ( x, y; ω) = YA ( y ) X A ( x; ω)

(38)

where

3.2. Derivation of wB(x,y) As shown in Fig. 4, Problem B has fixed boundary conditions at two opposite edges x = – lx/2 and lx/2, and the finite composite plate element is now divided into Nx finite strip elements in the x-direction. As Problem B can be obtained from Problem A by simply rotating the coordinate system (x, y) 90o degrees clock-wise, the solution procedure to obtain wB(x,y) for Problem B is identical to that used previously to obtain the solution wA(x,y) for Problem A, except for the change of coordinate system: x to y and y to x. Thus, without repeating the same solution procedure, the solution wB(x,y) can be written as follows: wB ( x, y ) = X B ( x)YB ( y; ω)d B = N B ( x, y; ω)d B

(39)

 d  d B =  B (1)   d B ( 2) 

(40)

where

and X B ( x) = [ L(B1) ( x), L(B2) ( x),
, L(Be) ( x),
, L(BN x −2) ( x), L(BN x −1) ( x)] YB ( y; ω ) = RB E B ( y) H B (ω ) −1

(41)

The symbols used in Eq. (40) are defined as follows: d B ( k )= { w B(1()k ) , w B( 2()k ) ,
, w B( e()k ) ,
, w B( N( kx −) 1) , ( N x −1) T (1) ( 2) (e) θ yB ( k ) , θ yB ( k ) ,
, θ yB ( k ) ,
, θ yB ( k ) }

(k = 1, 2)

(42)

where wB( e ) ( xe ,∓ l2y ) w = w ( ∓ ) =  (e ) l  θ xB ( xe ,∓ 2y )  l l (e) ∂w (xBe ) (∓ 2y ) θ yB ( xe ,∓ 2y )  (e) θ yB = =  (e) (1, 2 ) l  ∂x θ xyB ( xe ,∓ 2y ) (e) B (1, 2 )

(e ) xB

ly 2

Similarly the symbols used in Eq. (41) are defined as follows:

(43)

L(Be ) ( x) = hB( e ) ( x) X B(e2) ( x) + hB(e+1) ( x) X B( e1+1) ( x) RB = [rB(1) , rB( 2) , rB (3) , …, rB( j ) , …, rB(8( N x −1)) ] E B ( y; ω ) = diagonal [e

k y ( j ) y − 12 k y ( j )l y

] ( j = 1, 2, 3,
, 8( N x − 1))

(44)

 RB E B (− l2y )   l  RB K B E B (− 2y )  H B (ω ) =  R E (+ l y )   B B 2 ly   RB K B E B (+ 2 )

where X B( e1) = [ lx( e )−3 ( x − xe )2 (2 x − 3 xe−1 + xe ), l x(e )−2 ( x − xe−1 )( x − xe ) 2 ] X B( e2) = [ lx( e )−3 ( x − xe−1 ) 2 (2 x + xe−1 − 3xe ), l x( e )−2 ( x − xe )( x − xe−1 ) 2 ]

(45)

(e) B

h ( x) = H ( x − xe−1 ) − H ( x − xe )

and

k y( j)

 + k y( j )  = − k y ( j )  0 

if Re(k y ( j ) ) > 0 if Re( k y ( j ) ) < 0

(46)

if Re(k y ( j ) ) = 0

The jth wavenumber ky(j) and the corresponding wave mode vector rB(j) are determined from the eigenvalue problem given by [ AB 4 k y4( j ) + AB3 k y3( j ) + AB 2 k y2( j ) + AB1k 1y ( j ) + AB 0 − ω 2 M B ] rB ( j ) = 0

(47)

where AB1 = 2 D16 ( Λx 5 − ΛxT5 )

AB 0 = D11 Λx 6 ,

AB 2 = D12 ( Λx 3 + ΛxT3 ) − 4 D66 Λx 4 , AB 3 = 2 D26 ( Λx 2 − ΛxT2 )

(48)

M B = ρΛx1

AB 4 = D22 Λx1 , with Λx1 = ∫

Lx / 2

− Lx / 2

Λx3 = ∫

Lx / 2

Λx5 = ∫

Lx / 2

− Lx / 2

− Lx / 2

X T Xdx, 2

Λx 2 = ∫

Lx / 2

− Lx / 2

Lx / 2 ∂ X dx, Λx 4 = ∫ 2 − Lx / 2 ∂x 2 T Lx / 2 ∂ X ∂X dx, Λx 6 = ∫ 2 − Lx / 2 ∂x ∂x

XT

∂X dx ∂x ∂X T ∂X dx ∂x ∂x ∂2 X T ∂ 2 X dx ∂x 2 ∂x 2

XT

(49)

The constant matrices Λxj (j = 1, 2, 3, …, 6) are provided in Appendix A.

3.3. Derivation of w ( x, y )

The free wave solution w( x, y ) for the original problem shown in Fig. 2(a) can be obtained by substituting Eq. (37) and Eq. (39) into Eq. (12) as follows: w( x, y ) = N A ( x, y; ω )d A + N B ( x, y; ω )d B = N ( x, y; ω ) d

(50)

where d(ω) is the 8(Nx+Ny-2)-by-one spectral nodal DOFs vector defined by d  d =  A d B 

(51)

and N(x, y; ω) is the one-by-8(Nx+Ny-2) dynamic shape function matrix defined by N ( x , y; ω ) = [ N A ( x, y; ω ) N B ( x, y ; ω ) ]

(52)

4. Formulation of spectral element equation The spectral element equation can be formulated from the weak form of the governing differential equation of motion (Eq. (7)) which is given by  ∂ 2w ∂2w ∂ 2w  ∂ 2w δ ∫x ∫y  D11 2 + D12 2 + 2 D16 ∂x ∂y ∂x∂y  ∂x 2   ∂ 2w ∂2w ∂ 2w  ∂ 2w δ +  D12 2 + D22 2 + 2 D26 ∂x ∂y ∂x∂y  ∂y 2    ∂ 2w ∂ 2w ∂2 w  ∂ 2w δ + 2 D16 2 + D26 2 + 2 D66 − ω 2 ρwδw dxdy ∂x ∂y ∂x∂y  ∂x∂y   = ∫x ∫y f ( x, y)δw dxdy + ∫y Vx1δw(− Lx 2 , y)dy + ∫y Vx 2δw( Lx 2 , y )dy + ∫x V y1δw( x, − Ly 2)dx + ∫x V y 2δw( x, L y 2)dx + ∫y M x1δ (∂w(− Lx 2 , y ) / ∂x)dy + ∫y M x 2δ (∂w( Lx 2 , y ) / ∂x)dy + ∫x M y1δ (∂w( x, − Ly 2) / ∂y)dx + ∫x M y 2δ (∂w( x, Ly 2) / ∂y )dx

(53)

where Vx1(y), Vx2(y), Vy1(x), and Vy2(x) are the transverse shear forces applied on four boundary edges, and Mx1(y), Mx2(y), My1(x), and My2(x) are the bending moments applied on four boundary edges. By substituting Eq. (50) into Eq. (53), we can obtain the spectral element equation as follows:

S (ω )d (ω ) = f1 (ω ) + f 2 (ω )

(54)

where

S (ω ) =

∫

x

  ∂2 N T ∂2 N   ∂2 N T ∂2 N ∂2 N T ∂2 N      D + D +  11 12 ∫y   ∂x 2 ∂x 2  2 2 2 2  ∂ x ∂ y ∂ y ∂ x    2 T 2 2 T 2 2 T 2 ∂ N ∂ N ∂ N ∂ N  ∂ N ∂ N   + D22   + 2 D16  + 2 2  2 ∂y 2   ∂x ∂x∂y ∂x∂y ∂x   ∂y  ∂2 N T ∂2 N ∂2 N T ∂2 N   ∂2 N T ∂2 N     + 2 D26  + + 4 D 66 2 2   ∂y ∂x∂y ∂x∂y ∂y   ∂x∂y ∂x∂y 

(55)

]

− ω 2 ρN T N dxdy

and f1 (ω ) = ∫ ∫ f ( x, y) N T ( x, y)dxdy x y

f 2 (ω ) = ∫ Vx1 N T (− 12 l x , y )dy + ∫ Vx 2 N T ( 12 l x , y )dy y

y

+ ∫ Vy1 N T ( x,− 12 l y )dx + ∫ Vy 2 N T ( x, 12 l y )dx x

x

(56)

+ ∫ M x1 (∂N T (− 12 l x , y ) / ∂x)dy + ∫ M x 2 (∂N T ( 12 l x , y ) / ∂x)dy y

y

+ ∫ M y1 (∂N T ( x,− 12 l y ) / ∂y )dx + ∫ M y 2 (∂N T ( x, 12 l y ) / ∂y)dx x

x

Substituting the dynamic shape function matrix N(x, y; ω) given by Eq. (52) into Eq. (55) and taking some mathematical manipulations gives the 8(Nx+Ny-2)-by8(Nx+Ny-2) dynamic stiffness matrix S(ω), which is often called spectral element matrix, in the following form:

S (ω ) = H − T (ω )[ D1 (ω ) + D2 (ω ) + D2T (ω ) + D3 (ω ) + D3T (ω ) + D4 (ω ) + D5 (ω ) + D5T (ω ) + D6 (ω ) − ω 2 Dρ (ω )] H −1 (ω )

(57)

where 0   H ( ω) H ( ω) =  A H B (ω)  0

(58)

and ( K 2 R T Λ R K 2 ). * Γ 1 ( Γ 4 RB ). * ( K A2 RAT Γ 6T ) D1 = D11  A A y1 A A  symmetric ( RBT Λx 6 RB ). * Γ 5   ( K 2 R T Λ R ). * Γ ( Γ 2 RB K B2 ). * ( K A2 R AT Γ 6T ) D2 = D12  A A y 3 A T T 1  ( RBT ΛxT3 RB K B2 ). * Γ 5   ( Γ 8 R A ). * ( RB Γ 4 ) ( K 2 R T Λ R K ). * Γ ( Γ 3 RB K B ). * ( K A2 R AT Γ 6T ) D3 = 2 D16  A A y 2 A A T T 1  ( RBT Λx 5 RB K B ). * Γ 5   ( Γ 7 R A K A ). * ( RB Γ 4 )

( R AT Λy 6 R A ). * Γ 1 ( Γ 2 RB K B2 ). * ( R AT Γ 8T )  D4 = D22   ( K B2 RBT Λx1 RB K B2 ). * Γ 5   symmetric

(59)

 ( R AT Λy 5 R A K A ). * Γ 1 ( Γ 2 RB K B ). * ( R AT Γ 8T )  D5 = 2 D26   2 T T 2 T ( Γ 7 R A K A ). * ( K B RB Γ 2 ) ( K B RB Λx 2 RB K B ). * Γ 5  ( K A R AT Λy 4 RA K A ). * Γ 1 ( Γ 3 RB K B ). * ( K A R AT Γ 7T ) D6 = 4 D66   symmetric ( K B RBT Λx 4 RB K B ). * Γ 5   ( R T Λ R ). * Γ 1 ( Γ 2 RB ). * ( RAT Γ 6T ) Dρ = ρ  A y1 A  ( RBT Λx1 RB ). * Γ 5   symmetric

where the symbol (.*) denotes the element-wise matrix multiplication operator defined in MATLAB [23] (see Appendix B), and the constant matrices Γα (α = 1, 2, …, 8) are defined by Γ1 = ∫ vect( E A ) vect( E A )T dx, x

∂X B dx, ∂x vect( E B ) vect( E B )T dy,

Γ 2 = ∫ vect( E A ) X B dx x

∂2 XB dx ∂x 2 vect( E B )YAdy

Γ 3 = ∫ vect( E A )

Γ 4 = ∫ vect( E A )

Γ5 = ∫

Γ6 = ∫

x

y

Γ 7 = ∫ vect( E B ) y

∂YA dy, ∂y

x

y

Γ 8 = ∫ vect( E B ) y

(60)

∂ 2YA dy ∂y 2

The constant matrices Γα (α = 1, 2,…, 8) are provided in explicit forms in Appendix B.

5. Numerical results and discussion For the numerical studies in this study, we considered simply supported symmetric square composite plates as shown in Fig. 5(a). The square composite plates have the same dimensions 1 m×1 m× 0.001m, and they all consist of eight layers. Each layer was assumed to have the same material properties as follows: E1 = 144.48 GPa, E2 = 9.632 GPa, E3 = 9.412 GPa, G12 = 4.128 GPa, G13 = 7.457 GPa, G23 = 6.516 GPa, ν12 = 0.3, ν13 = 0.3, ν23 = 0.49, and ρb = 1389.2 kg/m3. Figure 6 shows five types of lay-ups considered in this study: two cross-ply lay-ups (Figs. 6(a) and 6(d)) and three angle-ply lay-ups (Figs. 6(b), 6(c), and 6(e)). The bending stiffnesses Dij (i, j = 1, 2, 6) computed for each type of lay-up are displayed in Table 1. At first, to evaluate the accuracy of the present spectral element model, the natural frequencies of a simply supported square composite plate whose lay-up is [02/902]S were obtained by the present spectral element model and the standard finite element model introduced in Appendix D and these results were compared in Table 2 with exact analytical results given by [1]

ωmn =

π 2l

1 2

ρ

D11m 4 + 2( D12 + 2 D66 )m 2 n 2 + D22n 4 (Hz)

(61)

where (m, n) denote the mode numbers, and l is the dimension of a square composite plate. For the results by standard FEM (denoted by ‘FEM’ in Table 2), the square composite plate was represented by the finite element model as shown in Fig. 5(b). Similarly, for the results by SEM (denoted by ‘SEM’ in Table 2), the square composite plate was represented by the spectral element model as shown in Fig. 5(c). Notice that the present spectral element model has nodes (i.e., nodal DOFs) only on four boundary edges while the standard finite element model has a huge number of additional nodes inside the composite plate. Table 2 shows that the results obtained by both FEM and SEM converge to the exact results as we increase the total number of nodes used in the analysis by decreasing the size of grids (or meshes) shown in Figs. 5(b) and 5(c). In addition, Table 2 also shows that the present spectral element model can provide nearly exact results by using only an extremely small number of nodes (i.e., nodal DOFs) when compared to that used for the FEM results. Table 3 shows the natural frequencies obtained by the present spectral element

model for the simply supported square composite plates with five types of lay-ups. For the results shown Table 3, total of 200 nodes were used. The number (m, n) in the parenthesis beside each natural frequency denotes the corresponding mode number. Table 3 shows that, as we can readily expect, the cross-ply lay-ups [0 2/902]S and [90 2/0 2]S have the same natural frequencies in the order of magnitude, but with switching the mode number (m, n) into (n, m). Form Table 3, we can find that the natural frequencies of a composite plate in general vary depending on the lay-up of composite lamina. To investigate the dynamic responses and wave propagation in the simply supported square composite plates with various types of lay-ups as shown in Fig. 6, an impulse point force was applied at the middle of each composite plate. The dynamic responses were then predicted in the various direction (θ), but at equal radial distance r = 0.25 m from the middle of the composite plates as illustrated in Fig. 7. Figure 8 compares the dynamic responses predicted at the radial distance r = 0.25 m in various directions represented by the direction angle θ (degrees) for four types of lay-ups: [0 2/902]S, [302/-602]S , [452/-452]S, and [0/90/45/-45]S. The dynamic responses were normalized with respect to the maximum value of the dynamic response predicted in the direction θ = 0 (degrees) for the lay-up [02/902]S. It is obvious from Fig. 8 that the dynamic responses in a specific direction strongly depend on the types of layups and this observation can be further clearly investigated from Fig. 9. It is also interesting to investigate from Fig. 8 that the very early parts of the dynamic responses in the directions θ = 0, 30, and 45 (degrees) for the lay-ups [0 2/902]S, [302/-602]S , and [452/-45 2]S, respectively, are nearly identical to each other. However, the dynamic responses become different from each other after some times (about 0.32 ms) at which the backward propagating waves reflected from boundary edges start being superposed with the original forward propagating waves. Finally, the waves that are propagating in the radial direction of the simply supported square composite plates with five types of lay-ups ([0 2/90 2]S, [302/-602]S , [452/45 2]S, [902/0 2]S, and [0/90/45/-45]S) are shown in Fig. 10. The propagating waves were predicted at the times t = 0.1 ms, 0.5 ms, 1.0 ms, and 2.0 ms. Figure 10 shows that the shapes of wave fronts strongly depend on the types of lay-ups. In general, the wave velocity will be faster in the direction of the higher bending stiffness than in the direction of lower bending stiffness. Accordingly, by investigating the bending stiffnesses

provided in Table 1, we can find that the wave front for the lay-up [02/902]S certainly moves faster in the x-direction than in the y-direction, and vice versa for the lay-up [902/02]S.

6. Conclusions In this study, the followings have been conducted: (1) A frequency-domain spectral element model is developed for the symmetric composite plates which have finite dimensions in two orthogonal directions, i.e., in the x- and y-directions. (2) To obtain the displacement field in a rectangular finite element that is required for the spectral element formulation, the splitting of original boundary conditions is used in combination with the so-called super spectral element method in which both the Kantorovich method-based finite strip element method and the frequency-domain waveguide method are used. (3) The spectral element model developed in this study has nodal points (i.e., nodal DOFs) only on four edges of a rectangular finite element, but no nodes inside the finite element. As the result, the total number of DOFs to be used in SEM will be drastically decreased to lead to a significant decrease of the computation cost, when compared with standard FEM. (4) Numerical studies were conducted to verify the high accuracy of the present spectral element model by the comparison with the results by exact theory and the standard FEM developed separately in this study.

(5) By using the present spectral element model, the effects of the lay-ups on the natural frequencies, dynamic responses, and wave propagations are investigated.

Acknowledgement This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (Grant No: NRF-2015R1A2A2A01003653)

Appendix A: Constant Matrices Defined in Eqs. (23) and (49)

 Λq(11)  (2)T  Λq1 (1)  lq 0 Λq1 =  420    0   0  Λ  − Λ 1  0 Λq 2 =  60    0   0

(1) q2 (2)T q2

Λ  Λ  0     0   0

(1) q3 ( 2)T q3

Λq 3 =

1 30lq(1)

Λq( 12 ) Λq(11) Λq( 12 ) T  0 0

0 Λq(12 ) Λq(11)  0 0

( 2) q2 (1) q2 (2)T q2


0
0
0  
Λq(11)
Λq( 12 ) T

( 2) q3 (1) q3 (2)T q3

0
Λq( 23)
Λq(13)


Λ Λ Λ

 0 0

 0 0

2 ( N q −1)×2 ( N q −1)

0   0  0     Λq( 22)   Λq(12) 


0
0
0  
Λq(12)
− Λq( 22) T

0 Λq( 22) Λq(12)  0 0

Λ Λ −Λ  0 0

0   0  0     Λq(12 )   Λq(11) 

2 ( N q −1)×2 ( N q −1)

       Λq( 23)   Λq(13)  2 ( N q −1)×2 ( N q −1)

0 0 0

0 0 0

 
Λq(13)
Λq( 23) T

Λq 4 = − Λq 3  Λq(15)  ( 2)T − Λq 5 1  0 Λq 5 = (1)  2lq    0   0  Λq(16)  (2)T  Λq 6 2  0 Λq 6 = (1) 3  lq    0   0

Λq( 25) Λq(15) − Λq( 25) T 

0 Λq( 25) Λq(15) 

0 0

0 0

Λq( 26) Λq(16) Λq( 26) T  0 0

0 Λq( 26) Λq(16)  0 0






0 0 0 


Λq(15)
− Λq( 25) T


0
0
0  
Λq(16)
Λq( 26) T

0   0  0     Λq( 26)   Λq(16) 

      Λq( 25)   Λq(15)  2 ( N q −1)×2 ( N q −1) 0 0 0 

2 ( N q −1)×2 ( N q −1)

where q = x and y, and the following sub-matrices are defined:

(A1)

 54 Λq(12 ) =  (1) 13l q

0  312 Λq(11) =  (1) 2  ,  0 8lq 

− 13lq(1)   − 3lq(1) 2 

 0  30 − 6lq(1)  12lq(1)  ( 2) Λq(12) =  , Λ =   (1)  q2 (1) 0  − lq(1) 2  6lq  − 12lq 0   36 − 3lq(1)   − 72 (2) Λq(13) =  , Λ =  (1)  q3 − 8lq(1) 2  lq(1) 2   0 3lq 0   24 Λq(14) =  (1) 2  ,  0 8l q 

 − 12 Λq( 24) =  (1) − 6lq

0 − 4 Λq(15) =  , 4 0 

2 0 Λq( 25) =  (1)  − 2 l q   − 6 3lq(1)  Λq( 26) =   (1) lq(1) 2  − 3lq

0  12 Λq(16) =  (1) 2  ,  0 4l q 

(A2)

6lq(1)   2lq(1) 

Appendix B: Constant Matrices Defined in Eq. (60) (i, j = 1, 2, 3,…, 8( N y − 1))

Γ1 = [ Γ1ij ]

Γ m = [ Γ mij ] (m = 2, 3, 4; i = 1, 2, 3,…, 8( N y − 1); j = 1, 2, 3,…, ( N x − 1)) Γ 5 = [ Γ5ij ]

(i, j = 1, 2, 3,…, 8( N x − 1))

Γ n = [ Γ nij ] (n = 6, 7, 8; i = 1, 2, 3,…, 8( N x − 1); j = 1, 2, 3,…, ( N y − 1)) where − (k +k )l  l e 2 x(i) x( j) x  − 1 ( k x ( i ) + k x ( j ) ) l x x 1 ( k x ( i ) + k x ( j ) )l x − 1 ( k x ( i ) + k x ( j ) )l x =  (e 2 − e2 )e 2  k x (i ) + k x ( j )  1

Γ1ij

Γ 2ij = k x−(4i ) [( p x( 1j +1) − p x( 3j ) )e

k x ( i ) x j −1 − 12 k x ( i ) l x

( p x( 5j +1) + p x( 7j ) )e

Γ 4ij = k x−(2i ) [( p x( 3j +1) − p x( 3j ) )e

,

k x ( i ) x j − 12 k x ( i ) l x

k x ( i ) x j − 12 k x ( i ) l x

+ (− p x( 8j +1) − p x( 10j ) )e

,

](1×2)

k x ( i ) x j − 12 k x ( i ) l x

+ ( p x( 8j +1) + p x( 6j ) )e

+ (− p x( 4j +1) + p x( 4j ) )e

k x ( i ) x j −1 − 12 k x ( i ) l x

k x ( i ) x j − 12 k x ( i ) l x

k x ( i ) x j − 12 k x ( i ) l x

+ ( p x( 4j +1) − p x( 4j ) )e

k x ( i ) x j −1 − 12 k x ( i ) l x

k x ( i ) x j −1 − 12 k x ( i ) l x

( p x( 9j +1) + p x( 7j ) )e

− ( p x( 8j +1) + p x( 6j ) )e

k x ( i ) x j −1 − 12 k x ( i ) l x

(− p x( 5j +1) − p x( 7j ) )e

if k x ( i ) + k x ( j ) ≠ 0

+ (− p x( 4j +1) + p x( 2j ) )e

k x ( i ) x j −1 − 12 k x ( i ) l x

Γ 3ij = k x−(3i ) [(− p x( 3j +1) + p x( 3j ) )e

if k x (i ) + k x ( j ) = 0

](1× 2)

,

k x ( i ) x j − 12 k x ( i ) l x

](1×2 )

(B1)

− (k +k )l  l e 2 y (i ) y ( j ) y  − 1 ( k y ( i ) + k y ( j ) )l y y 1 ( k y ( i ) + k y ( j ) )l y − 1 ( k y ( i ) + k y ( j ) )l y =  (e 2 − e2 )e 2  k y (i) + k y ( j )  1

Γ 5ij

Γ 6ij = k y−(4i ) [( p (y1j +1) − p (y 3j ) )e

k y ( i ) y j −1 − 12 k y ( i ) l y

( p (y 5j +1) + p (y 7j ) )e

(− p (y 5j +1) − p (y 7j ) )e Γ 8ij = k y−(2i ) [( p (y 3j +1) − p (y 3j ) )e

( p (y 9j +1) + p (y 7j ) )e

+ ( p (y 4j +1) − p (y 4j ) )e

k y ( i ) y j −1 − 12 k y ( i ) l y

k y ( i ) y j −1 − 12 k y ( i ) l y

,

k y ( i ) y j − 12 k y ( i ) l y

k y ( i ) y j − 12 k y ( i )l y

+ ( p (y8j +1) + p (y 6j ) )e

+ (− p (y 4j +1) + p (y 4j ) )e

k y ( i ) y j −1 − 12 k y ( i ) l y

k y ( i ) y j − 12 k y ( i ) l y

+ (− p (y 8j +1) − p (y 6j ) )e

k y ( i ) y j −1 − 12 k y ( i ) l y

(B2)

if k y (i ) + k y ( j ) ≠ 0

+ (− p (y 4j +1) + p (y 2j ) )e

k y ( i ) y j −1 − 12 k y ( i ) l y

Γ 7ij = k y−(3i ) [(− p (y 3j +1) + p (y 3j ) )e

if k y ( i ) + k y ( j ) = 0

j) + (− p (y 8j +1) − p (y10 )e

,

k y ( i ) y j − 12 k y ( i ) l y

k y ( i ) y j − 12 k y ( i ) l y

](1×2)

](1× 2)

,

k y ( i ) y j − 12 k y ( i ) l y

](1×2 )

with

p

( j) x1

p

( j) x3

p

( j) x5

p

( j) x7

= = = =

p x( 9j ) =

p (y1j ) = p (y 3j ) = p (y 5j ) = p

( j) y7

p

( j) y9

= =

12 + 6k x (i )l x( j ) − k x3( j )l x( j )3 l x( j )3 12 + 6k x (i )l x( j ) l x( j )3

,

6 + 4k x ( i )l x( j ) + k x2(i )l x( j ) 2 l x( j ) 2 6 + 2k x ( i )l x( j ) l x( j ) 2 6 + 4k x ( i )l x( j ) l x( j ) 2

p , p

l x( j ) 3 6 − 4k x (i ) l x( j ) + k x2(i ) l x( j ) 2 l x( j ) 2 6 − 2k x ( i )l x( j )

=

p x( 10j ) =

l y( j )3

,

6 + 4k y (i )l y( j ) + k y2(i )l y( j ) 2 l y( j ) 2

l y( j ) 2

( j) x8

=

12 − 6k x ( i )l x( j )

,

12 + 6k y (i ) l y( j )

6 + 4k y (i )l y( j )

( j) x6

=

l x( j ) 3

p

l y( j )3

l y( j ) 2

( j) x4

=

12 − 6k x (i ) l x( j ) − k x3( i )l x( j ) 3

,

12 + 6k y (i ) l y( j ) − k y3( i )l y( j ) 3

6 + 2k y (i )l y( j )

, p

( j) x2

, ,

, p (ye2) = p (y 4j ) =

, p (y 6j ) = p

( j) y8

p

( j) y10

= =

l x( j ) 2 6 − 4k x (i )l x( j ) l x( j ) 2

12 − 6k y ( i )l y( j ) − k y3( i )l y( j ) 3 l y( j ) 3 12 − 6k y (i )l y( j ) l y( j )3 6 − 4k y ( i )l y( j ) + k y2( i )l y( j ) 2 l y( j ) 2 6 − 2k y (i )l y( j ) l y( j ) 2 6 − 4k y ( i )l y( j ) l y( j ) 2

(B3)

Appendix C: Element-wise matrix multiplication [23] The symbol (.*) used in Eq. (59) denotes the element-wise matrix multiplication that is introduced in MATLAB as follows: A. * B = C

(C1)

where the components of the matrix C are defined by Cij = Aij × Bij

(C2)

Appendix D: Finite element model Figure D-1 shows a rectangular symmetric composite plate element with four-node points, one at each corner. Assuming there are three degrees of freedom (DOFs) at each node (i.e., w, θx, θy), the finite element model is derived in the following form: Md + Kd = f (t )

(D1)

where d(t) and f(t) are the nodal DOFs vector and nodal forces vectors, respectively, defined by

d (t ) = { w1 , θ x1, θ y1 , w2 , θ x2 , θ y 2 , w3 , θ x3 , θ y 3 , w4 , θ x4 , θ y 4}T f (t ) = {V1 , M x1, M y1 , V2 , M x2 , M y 2 , V3 , M x3 , M y 3 , V4 , M x4 , M y 4 }T

(D2)

and M and K are the finite element mass matrix and finite element stiffness matrix, respectively, defined by M = ρ ∫ N T Ndxdy A

K = D11 K 1 + D22 K 2 + D12 ( K 3 + K 3T ) + 2 D16 ( K 4 + K 4T ) + 2 D26 ( K 5 + K 5T ) + 4 D66 K 6

where Ki (i = 1, 2, 3, 4, 5, 6) are defined by

(D3)

K1 = ∫

∂2 N T ∂2 N ∂ 2 N T ∂2 N dxdy , K = 2 ∫A ∂y 2 ∂y 2 dxdy ∂x 2 ∂x 2

K3 = ∫

∂2 N T ∂ 2 N ∂2 N T ∂2 N dxdy , K = 4 ∫A ∂x 2 ∂x∂y dxdy ∂x 2 ∂y 2

K5 = ∫

∂2 N T ∂2 N ∂2 N T ∂2 N dxdy , K = 6 ∫A ∂x∂y ∂x∂y dxdy ∂y 2 ∂x∂y

A

A

A

(D4)

The one-by-twelve shape function matrix N(x,y) used in Eqs. (D3) and (D4) is defined by

N ( x, y ) = [ N1 , N 2 , N 3 , N 4 , N 5 , N 6 , N 7 , N 8 , N 9 , N10 , N11 , N12 ]

(D5)

where N 1 ( x, y ) = α(l x − 2 x)(l y − 2 y)(l x2l y2 − l x l y2 x − l x2l y y − 2l y2 x 2 − 2l x2 y 2 ) N 2 ( x, y) = β (l x + 2 x)(l y − 2 y )(l x − 2 x) 2 N 3 ( x, y ) = γ (l x − 2 x)(l y + 2 y )(l y − 2 y ) 2 N 4 ( x, y) = α(l x + 2 x)(l y − 2 y)(l x2l y2 + l x l y2 x − l x2 l y y − 2l y2 x 2 − 2l x2 y 2 ) N 5 ( x, y ) = − β (l x − 2 x)(l y − 2 y )(l x + 2 x) 2 N 6 ( x, y ) = γ (l x + 2 x)(l y + 2 y )(l y − 2 y) 2 N 7 ( x, y) = α(l x + 2 x)(l y + 2 y )(l x2l y2 + l x l y2 x + l x2l y y − 2l y2 x 2 − 2l x2 y 2 )

(D6)

N 8 ( x, y ) = − β (l x − 2 x)(l y + 2 y )(l x + 2 x) 2 N 9 ( x, y) = −γ (l x + 2 x)(l y − 2 y )(l y + 2 y) 2 N 10 ( x, y ) = α (l x − 2 x)(l y + 2 y )(l x2l y2 − l x l y2 x + l x2l y y − 2l y2 x 2 − 2l x2 y 2 ) N 11 ( x, y ) = β (l x + 2 x)(l y + 2 y )(l x − 2 x) 2 N 12 ( x, y ) = −γ (l x − 2 x)(l y − 2 y )(l y + 2 y ) 2

with α=

1 1 1 , β= , γ= 3 3 2 4l x l y 16l x l y 16l x l y2

(D7)

References [1] Reddy JN. Mechanics of laminated composite plates, theory and analysis. New York: CRC Press; 1997. [2] Gorman DJ. Free vibration analysis of rectangular plates. New York: Elsevier; 1982. [3] Doyle JF. Wave propagation in structures. New York: Springer; 1999. [4] Lee U. Spectral element method in structural dynamics. Singapore: John Wiley & Sons; 2009. [5] Patera AT. A spectral element method for fluid dynamics: laminar flow in a channel expansion. J Comput Phys 1984;54(3):468-88. [6] Berçin AN. Analysis of orthotropic plate structures by the direct-dynamic s Mmethod. Mech Res Commun 1995;22(5):461-6. [7] Leung AYT, Zhou WE. Dynamic stiffness analysis of laminated composite plates. Thin Wall Struct 1996;25(2):109-33. [8] Orrenius U, Finnveden S. Calculation of wave propagation in rib-stiffened plate structures. J Sound Vib 1996;198(2):203-24. [9] Lee U, Lee J. Spectral element method for Levy-type plates subject to dynamic loads. J Eng Mech 1999;125(2):243-7. [10] Krawczuk M, Palacz M, Ostachowicz W. Wave propagation in plate structures for crack detection. Finite Elem Anal Des 2004;40(9-10):991-1004. [11] Casimir JB, Kevorkian S, Vinh T. The dynamic stiffness matrix of twodimensional elements: Application to Kirchhoff’s plate continuous elements. J Sound Vib 2005;287(3):571-89. [12] Birgersson F, Finnveden S, Nilsson CM. A spectral super element for modeling of plate vibration. Part 1: General theory. J Sound Vib 2005;287(1-2):297-314. [13] Chakraborty A, Gopalakrishnan S. A spectral finite element model for wave propagation analysis in laminated composite plate. J Vib Acoust 2006;128 (4):477-88. [14] Ostachowicz WM. Damage detection of structures using spectral finite element method. Comput Struct 2008;86(3-5):454-462.

[15] Barbieri E, Cammarano A, De Rosa S, Franco F. Waveguides of a composite plate by using the spectral finite element approach. J Vib Control 2009;15 (3):347-67. [16] Boscolo M, Banerjee JR. Dynamic stiffness elements and their appplications for plates using first order shear deformation theory. Comput Struct 2011;89(34):395-410. [17] Hajheidari H, Mirdamadi HR. Free and transient vibration analysis of an unsymmetric cross-ply laminated plate by spectral finite elements. Acta Mech 2012;223(11):2477-92. [18] Hajheidari H, Mirdamadi HR. Frequency-dependent vibration analysis of symmetric cross-ply laminated plate of Levy-type by spectral element and finite strip procedures. Appl Math Model 2013;37(12-13):7193-205. [19] Kulla PH. High precision finite elements. Finite Elem Anal Des 1997;26(2):97114. [20] Gupta KK, Meek JL. Finite element multidisciplinary analysis. Reston, Virginia: American Institute of Aeronautics and Astronautics; 2000. [21] Park J, Park I, Lee U. Transverse vibration and waves in a membrane: frequency domain spectral element modeling and analysis, Math Probl Eng 2014;2014:114. Article ID 642782. [22] Newland DE. Random vibrations, spectral and wavelet analysis. New York: Longman; 1993. [23] MATLAB User’s Guide, The MathWorks Inc., MA, USA, 1993.

Table Captions Table 1. The bending stiffnesses (GPa·mm3) of the symmetric composite plates with

various types of lay-ups. Table 2. Comparison of the natural frequencies (Hz) obtained by exact theory, FEM,

and SEM for the simply supported square composite plate with the lay-up [02/902]S. Table 3. Comparison of the natural frequencies (Hz) of the simply supported square

composite plates with various types of lay-ups.

Figure Captions Fig. 1. The geometry of a symmetric laminated composite plate. Fig. 2. Splitting of boundary conditions for a rectangular composite plate element

subjected to arbitrary boundary conditions. Fig. 3. For Problem A: (a) the finite strip element representation of a finite composite

plate element, and (b) a representative finite strip element: where ● = nodes. Fig. 4. For Problem B: (a) the finite strip element representation of a finite composite

plate element, and (b) a representative finite strip element: where ● = nodes. Fig. 5 . Example problem: (a) a simply supported square composite plate, (b) finite

element model, (c) spectral element model; where ● = nodes.. Fig. 6. Five types of lay-ups for the symmetric composite plates considered in this

study. Fig. 7 . The location (r,θ) in the simply supported square composite plates with vari-

ous lay-ups at which the dynamic responses are predicted, where r = 0.25 m.. Fig. 8 . The impulse-induced transient dynamic responses in the various directions θ

(degrees) of the simply supported square composite plates with four types of lay-ups: (a) [02/902]S; (b) [302/-602]S (c) [45 2/-452]S; and (d) [0/90/45/-45]S. Fig. 9 . The impulse-induced transient dynamic responses of the simply supported

square composite plates with three types of lay-ups [02/90 2]S, [30 2/-602]S , and [452/-452]S predicted in the directions: (a) θ = 0 (degrees), (b) θ = 30 (degrees), and (c) θ = 45 (degrees). Fig. 10 . The impulse-induced waves predicted at the times t = 0.1 ms, 0.5 ms, 1.0 ms,

and 2.0 ms for the simply supported square composite plates with five types of lay-ups: (a) [02/90 2]S, (b) [302/-60 2]S, (c) [452/-45 2]S, (d) [02/902]S , and (e) [0/90/45/-45]S. Fig. D-1. Four-node 12-DOF rectangular plate bending element.

Table 1

The bending stiffnesses (GPa·mm3) of the symmetric composite plates with various types of lay-ups. Lay-ups

Bending stiffnesses

[02/902]S

[302/-60 2]S

[452/-452]S

[902/02]S

[0/90/45/-45]S

D11 D22 D66 D12 D16 D26

10.70 2.221 0.344 0.242 0 0

6.506 2.257 2.418 2.316 3.033 0.639

3.695 3.695 3.109 3.007 2.120 2.120

2.221 10.70 0.344 0.242 0 0

7.705 4.525 0.690 0.588 0.265 0.265

Table 2

Comparison of the natural frequencies (Hz) obtained by exact theory, FEM, and SEM for the simply supported square composite plate with the lay-up [02/902]S. Modes (m, n)

Exact [1]

(1, 1) (1, 2) (2, 1) (1, 3) (2, 2) (2, 3) (1, 4) (2, 4) (3, 1) (3, 2)

5.124 9.764 17.92 19.19 20.49 27.25 32.89 39.05 39.66 41.48

FEM

SEM

nn=2,501 n n=10,201 n n=14,641

n n=40 n n=120 n n=200

5.122 9.759 17.92 19.18 20.47 27.20 32.87 38.98 39.65 41.43

5.123 9.762 17.92 19.19 20.49 27.24 32.88 39.03 39.66 41.48

Note: nn = total number of nodes used in the analysis

5.124 9.764 17.92 19.19 20.49 27.25 32.89 39.04 39.66 41.48

5.133 9.784 17.93 19.21 20.52 27.32 32.91 39.13 39.66 41.53

5.125 9.766 17.92 19.19 20.50 27.26 32.89 39.06 39.66 41.49

5.124 9.764 17.92 19.19 20.49 27.25 32.89 39.05 39.66 41.49

Table 3

Comparison of the natural frequencies (Hz) of the simply supported square composite plates with various types of lay-ups. Lay-ups [02/902]S

[302/-60 2]S

[452/-45 2]S

[902/0 2]S

[0/90/45/-45]S

5.124 (1,1) 9.764 (1,2) 17.92 (2,1) 19.19 (1,3) 20.49 (2,2) 27.25 (2,3) 32.89 (1,4) 39.05 (2,4) 39.66 (3,1) 41.49 (3,2)

5.827 (1,1) 11.76 (1,2) 16.33 (2,1) 20.25 (1,3) 24.17 (2,2) 30.46 (1,4) 33.28 (3,1) 36.69 (2,4) 41.41 (3,2) 42.79 (1,5)

6.083 (1,1) 12.30 (1,2) 15.98 (2,1) 20.61 (1,3) 26.25 (2,2) 30.46 (3,1) 30.83 (1,4) 38.28 (2,3) 43.09 (1,5) 44.98 (3,2)

5.124 (1,1) 9.764 (2,1) 17.92 (1,2) 19.19 (3,1) 20.49 (2,2) 27.25 (3,2) 32.89 (4,1) 39.05 (4,2) 39.66 (1,3) 41.49 (2,3)

5.343 (1,1) 13.00 (1,2) 15.96 (2,1) 21.32 (2,2) 26.97 (1,3) 33.23 (2,3) 34.33 (3,1) 38.60 (3,2) 46.70 (3,3) 47.67 (1,4)

Note: (m, n) denote the mode numbers.

z

y

z …

n-th lamina

ly

x

zk+1 zk

… …

k-th lamina

2nd lamina 1st lamina

lx

Fig. 1. The geometry of a symmetric laminated composite plate.

x

(a) Original problem

wyA(x,-ly/2)=0 0 (b) Problem A

w (-L /2,y)

wB(x,y)

wxB(lx/2,y)=0 0

w (-L /2,y)

w (-L /2,y)

+

wxB(-lx/2,y)=0 0

w (-L /2,y)

wA(x,y)

wyB(x,ly/2)=w wy2(x) wxA(lx/2,y)=w wx2(y)

wy(x,-ly/2)=w wy1 (x)

=

wxA(-lx/2,y)=w w x1 (y)

x

w (-L /2,y)

w(x,y)

wyA(x,ly/2)=0 0 w x(lx/2,y)=w w x2 (y)

y

w (-L /2,y)

w x(-lx /2,y)=w wx1(y)

wy(x,ly/2)=w wy2(x)

wyB(x,-ly/2)=w wy1 (x) (c) Problem B

Fig. 2. Splitting of boundary conditions for a rectangular composite plate element subjected to arbitrary boundary conditions.

y (Ny )

w y A (x) = w y 2 (x) = 0 (Ny −1)

wyA

(Ny)

 (e) ye ly dA1



(x)

w (e) y A (x)

w (e) y A (x) = w yA (x, y e )

w (ey A−1) (x)

ye-1

x

(3)

w (2) y A (x)

(2)

w (1) y A (x)

(1)

w (0) y A (x) = w y1(x) = 0

dA2

(e)

w (Ae ) (x, y)

l y(e)

−1) w (e y A (x) = w yA (x, y e −1 )

lx

Fig. 3. For Problem A: (a) the finite strip element representation of a finite composite plate element, and (b) a representative finite strip element: where ● = nodes.

y dB2

(e)

x

ly

xe xe−1

(1) (2) (3)

…

(e)

… (Nx) l x(e)

dB1

lx

Fig. 4. For Problem B: (a) the finite strip element representation of a finite composite plate element, and (b) a representative finite strip element: where ● = nodes.

lx = 1 m

Simply-supported (a) Example problem

ly = 1 m

Simply-supported

Simply-supported

Simply-supported

(b) Finite element model

(c) Spectral element model

Fig. 5. Example problem: (a) a simply supported square composite plate, (b) finite element model, (c) spectral element model; where ● = nodes.

0˚ 90˚ 90˚ 0˚

30˚ -60˚ -60˚ 30˚

45˚ -45˚ -45˚ 45˚

90˚ 0˚ 0˚ 90˚

(a) [02/902]S

(b) [302/-602] S

(c) [452 /-452] S

(d) [902/02]S

(e) [0/90/45/-45]S

Fig. 6. Five types of lay-ups for the symmetric composite plates considered in this study.

y r = 0.25m

r 1m

θ

x

1m

Fig. 7. The location (r,θ) in the simply supported square composite plates with various lay-ups at which the dynamic responses are predicted, where r = 0.25 m.

180

165

165

150

150

135

135

120

120

Direction θ (degrees)

Direction θ (degrees)

180

105 90 75

105 90 75

60

60

45

45

30

30

15

15

0

0

(a) [02/902]S

(b) [302/-602]S

Fig. 8. The impulse-induced transient dynamic responses in the various directions θ (degrees) of the simply supported square composite plates with four types of lay-ups: (a) [02/902]S, (b) [302/-602]S, (c) [452/-452]S, and (d) [0/90/45/-45]S.

180

165

165

150

150

135

135

120

120

Direction θ (degrees)

Direction θ (degrees)

180

105 90

75

105 90

75

60

60

45

45

30

30

15

15

0

0

(c) [452/-452]S

(d) [0/90/45/-45]S

Fig. 8 (Continued). The impulse-induced transient dynamic responses in the various directions θ (degrees) of the simply supported square composite plates with four types of lay-ups: (a) [0 2/90 2]S, (b) [30 2/-602]S, (c) [452/-452]S, and (d) [0/90/45/-45]S.

[02/902]S [302/-602] S [452/-452] S

(a) θ = 0 (degrees) [02/902]S [302/-602] S [452/-452] S

(b) θ = 30 (degrees) [02/902]S [302/-602] S [452/-452] S

(c) θ = 45 (degrees)

Fig. 9. The impulse-induced transient dynamic responses of the simply supported square composite plates with three types of lay-ups [02/902]S, [302/-602]S , and [452/45 2]S predicted in the directions: (a) θ = 0 (degrees), (b) θ = 30 (degrees), and (c) θ = 45 (degrees).

t = 0.1 ms

t = 0.5 ms

t = 1.0 ms

t = 2.0 ms

(a) [02/902]S t = 0.1 ms

t = 0.5 ms

t = 1.0 ms

t = 2.0 ms

(b) [302/-602]S t = 0.1 ms

t = 0.5 ms

t = 1.0 ms

t = 2.0 ms

(c) [452/-452]S t = 0.1 ms

t = 0.5 ms

t = 1.0 ms

t = 2.0 ms

(d) [902/02]S t = 0.1 ms

t = 0.5 ms

t = 1.0 ms

t = 2.0 ms

(e) [0/90/45/-45]S

Fig. 10. The impulse-induced waves predicted at the times t = 0.1 ms, 0.5 ms, 1.0 ms, and 2.0 ms for the simply supported square composite plates with five types of layups: (a) [02/90 2]S, (b) [30 2/-60 2]S, (c) [452/-452]S, (d) [02/902]S , and (e) [0/90/45/-45]S.

(w4, θx4, θy4)

(w3, θx3, θy3)

4

3

1

2

(w1, θx1, θy1)

(w2, θx2, θy2)

Fig. D-1. Four-node 12-DOF rectangular plate bending element.