Free vibration analysis of composite plates by higher-order 1D dynamic stiffness elements and experiments

Composite Structures 118 (2014) 654–663

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Free vibration analysis of composite plates by higher-order 1D dynamic stiffness elements and experiments A. Pagani a,b, E. Carrera a,⇑, J.R. Banerjee b, P.H. Cabral c, G. Caprio c, A. Prado c a

Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy School of Engineering and Mathematical Sciences, City University London, Northampton Square, EC1V 0HB London, United Kingdom c Embraer S.A., 12227-901 São José dos Campos, Brazil b

a r t i c l e

i n f o

Article history: Available online 19 August 2014 Keywords: Layered structures Free vibration analysis Carrera Unified Formulation (CUF) Dynamic Stiffness Method (DSM) Experimental tests

a b s t r a c t A novel approach for free vibration analysis of composite plate-like structures is introduced in this paper. Refined beam theories are formulated by making use of the Carrera Unified Formulation (CUF). By exploiting the hierarchical characteristics of CUF, the differential equations of motions and the natural boundary conditions are written in a compact and concise form in terms of fundamental nuclei, whose formal mathematical expressions do not depend on the order of the theory N. After the closed form solution of the N-th order beam model is sought, a general procedure to derive the exact Dynamic Stiffness (DS) matrix is devised by relating the amplitudes of the harmonically varying loads to those of the responses. The global DS matrices of composite laminated plates are then used with reference to the algorithm of Wittrick and Williams to carry out free vibration analyses. The accuracy of the proposed methodology is verified through published literature, finite element solutions from the commercial code MSC/ Nastran and experimental tests. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Laminated composite structures are widely employed in different industrial fields, such as aerospace and automotive engineering. Their success is mainly due to their high strength-to-weight ratio, high stiffness-to-weight ratio, ease of formability, wide range of operating temperatures, and their capability to be tailored according to a given requirement (see the book by Tsai [1]). However, with their anisotropy and non-homogeneity, appropriate theories are required to model multi-layer structures. Several approaches have been proposed in the literature, especially for plates/shells and beams, and a brief although not exhaustive review is given hereinafter. The former and probably the most known theory for composite plate modeling is the Classical Lamination Theory (CLT) [2], which is an extension of the classical plate theory based on the wellknown Love-Kirchoff hypotheses [3]. CLT neglects the effect of out-of-plane strains. However, one of the main issues related to the proper modeling of a composite structure is due to its low ⇑ Corresponding author. E-mail addresses: [email protected] (A. Pagani), [email protected] (E. Carrera), [email protected] (J.R. Banerjee), [email protected] (P.H. Cabral), [email protected] (G. Caprio), alex.prado@embraer. com.br (A. Prado). http://dx.doi.org/10.1016/j.compstruct.2014.08.020 0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.

transverse shear moduli compared to the axial tensile moduli. CLT is therefore inadequate for most of the practical cases. The First-order Shear Deformation Theory (FSDT) is based on the works by Reissner [4] and Mindlin [5] and it accounts for the shear deformation effects by linear variation of in-plane displacements. In order to overcome the limitations of classical theories, several refined models have been developed over the years. These include, for example, the Higher-order Shear Deformation Theories (HSDTs) such as the one developed by Reddy [6,7], the so-called zig-zag theories [8], the theories based on the Reissner’s Mixed Variational Theorem (RMVT) [9], and layer-wise models [10]. Further details about refined theories for composite plates modeling can be found in the review by Khandan et al. [11]. Plate models are suitable when one dimension of the structure is adequately smaller then the others. On the other hand, when slender structures with one dimension higher than the others are considered, beam models are preferable. The classical beam theories by Euler–Bernoulli [12] (hereinafter referred to as EBBM) and Timoshenko [13] (hereinafter referred to as TBM) are not suitable for the analysis of composite laminates as, similarly to CLT and FSDT, they do not correctly foresee shear effects. A comprehensive overview on composite beam works can be found in the excellent review of Kapania and Raciti [14,15]. Some recent noteworthy contributions about refined composite beam analysis are mentioned in the following. The attention is focused on

655

A. Pagani et al. / Composite Structures 118 (2014) 654–663

free vibration analysis, which is the main topic of the present work. A family of sinus models was presented by Vidal and Polit [16] for the vibration analysis of laminated beams. In [17], a trigonometric shear deformation theory was developed and the closed form solution was provided. Subramanian [18] presented two different one-dimensional (1D) Finite Elements (FEs) for laminated beams, in which a 5th order expansion was used to expand the axial displacement and a 4th order power series was used for the transverse displacement. In the work by Marur and Kant [19], Taylor’s series expansions were used for axial displacement in order to describe the warping of cross-sections of sandwich and composite beams. Interesting mixed formulations were presented in [20], where through-the-thickness continuity of transverse stress and displacement fields was enforced. The beam model presented by Heyliger and Reddy [21] accounted for the stress-free conditions on the upper and lower surfaces of the beam while retaining a parabolic shear strain distribution. Other noteworthy contributions are those by Chen et al. [22], Hodges et al. [23], Stemple and Lee [24], Mitra et al. [25], Chandrashekhara et al. [26] and Chandrashekhara and Bangera [27]. The present paper deals with 1D higher-order theories able to capture the mechanical behaviour of laminated orthotropic platelike structures. Refined beam models are developed within the framework of Carrera Unified Formulation (CUF) which was formerly developed for plate and shell models [28–30] and then recently extended to the analysis of beams [31–33]. CUF is a hierarchical formulation that considers the order of the model, N, as a free-parameter (i.e. as an input) of the analysis or in other words, it allows the governing equations (or equivalently, FE matrices) to be written in terms of fundamental nuclei, which do not depend on the beam theory N. In the present work, beam theories using CUF are obtained on the basis of Taylor-type expansions (TE) of the three-dimensional (3D) displacement field. EBBM and TBM can be obtained as particular or special cases. Several works about CUF TE 1D models can be found in the literature, see for example [34–38]. In majority of the papers about 1D CUF models, the Finite Element Method (FEM) has been used to handle arbitrary geometries and loading conditions. Recently, Giunta et al. [39] provided Navier-type closed solutions of 1D refined TE CUF models for modal analysis of composite structures. However, the use of the Navier solution limits the analysis to simply supported beams with cross-ply lamination schemes. In recent works, Pagani et al. [40,41] provided a more powerful approach for CUF TE theories through the application of the Dynamic Stiffness Method (DSM) to carry out the free vibration analysis of isotropic and composite beam structures. The DSM is appealing in dynamics because unlike the FEM, it provides exact solution of the equations of motion of a structure once the initial assumptions on the displacements field have been made and no restrictions applies to the boundary conditions, geometry and stacking sequence. DSM has been quite extensively developed for beam elements by Banerjee [42–46], Banerjee et al. [47], Eisenberger et al. [48], Abramovic et al. [49] and Williams and Wittrick [50]. Plate elements based on DSM were originally formulated by Wittrick [51] and Wittrick and Williams [52]. Recently, DSM has been applied to Mindlin [53,54] and HSDT plate assemblies [55]. In these latter works a more comprehensive review on the use of DSM can be found. In this work, 1D higher-order Dynamic Stiffness (DS) elements based on CUF are extended and applied to the free vibration analysis of laminated plates. In the next section, CUF is introduced and higher-order models are formulated. The principle of virtual displacements is then used to derive the equations of motion and the natural boundary conditions, which are subsequently expressed in the frequency domain by assuming a harmonic

solution. After the resulting system of ordinary differential equations of second order with constant coefficients is solved, the frequency dependent DS matrix of the system is derived. Finally, the algorithm of Wittrick and Williams [56] is applied to extrapolate the free vibration characteristics of laminated composite plate and the results are compared to those from the literature, commercial FE tools and experimental tests. 2. Carrera unified formulation 2.1. Preliminaries The adopted rectangular cartesian coordinate system is shown in Fig. 1, together with the geometry of a multi-layered structure. The cross-section of the beam lies on the xz-plane and it is denoted by X, whereas the boundaries over y are 0 6 y 6 L. Superscript k, which is usually used to denote variables and parameters related to the kth layer, is neglected in the following for the sake of simplicity. Let us introduce the transposed displacement vector,

 uðx; y; z; tÞ ¼ ux

uy

uz

T

ð1Þ

The stress, r, and strain, , components are expressed in transposed forms as follows:





r ¼ ryy rxx rzz rxz ryz rxy T ; 

 ¼ yy xx zz xz yz xy

T

ð2Þ

In the case of small displacements with respect to a characteristic dimension in the plane of X, the strain–displacement relations are

 ¼ Du

ð3Þ

where D is the following linear differential operator matrix

2

0

6@ 6 @x 6 60 6 D¼6@ 6 @z 6 60 4 @ @y

@ @y

0 0 0 @ @z @ @x

0

3

7 07 7 @ 7 @z 7 @ 7 7 @x 7 @ 7 @y 5 0

ð4Þ

Constitutive laws are now exploited to obtain stress components to give

r ¼ Ce 

ð5Þ

e is In the case of orthotropic material the matrix C

2

e 33 C 6e 6 C 23 6 6e C 13 e¼6 C 6 6 0 6 6 4 0 e 36 C

e 23 C e 22 C e 12 C

e 13 C e 12 C e 11 C

0

0

0 e 26 C

0 e 16 C

0

0

0

0

0 e 44 C e 45 C

0 e 45 C e 55 C

0

0

3 e 36 C e 26 7 7 C 7 7 e C 16 7 7 0 7 7 7 0 5 e 66 C

Fig. 1. Geometry and reference system for the laminated plate-like beam.

ð6Þ

656

A. Pagani et al. / Composite Structures 118 (2014) 654–663

e ij depend on Young’s and Poisson’s moduli as well as Coefficients C on the fiber orientation angle, h, that is graphically defined in Fig. 2 where ‘‘1’’, ‘‘2’’, and ‘‘3’’ represent the cartesian axes of the material. For the sake of brevity, the expressions for the coefficients e ij are not reported here, but can be found in standard texts, see for C example Tsai [1] and Reddy [7]. Within the framework of CUF, the displacement field uðx; y; z; tÞ can be expressed as

uðx; y; z; tÞ ¼ F s ðx; zÞus ðy; tÞ;

s ¼ 1; 2; . . . :; M

ð7Þ

where F s are the functions of the coordinates x and z on the crosssection. us is the vector of the generalized displacements, M stands for the number of the terms used in the expansion, and the repeated subscript, s, indicates summation. The choice of F s determines the class of the 1D CUF model. For example, TE (Taylor expansion) 1D CUF models consist of a MacLaurin series that uses the two-dimensional (2D) polynomials xi zj as F s functions, where i and j are positive integers. For instance, the displacement field of the secondorder ðN ¼ 2Þ TE model can be expressed as

ux ¼ ux1 þ x ux2 þ z ux3 þ x2 ux4 þ xz ux5 þ z2 ux6 uy ¼ uy1 þ x uy2 þ z uy3 þ x2 uy4 þ xz uy5 þ z2 uy6

ð8Þ

uz ¼ uz1 þ x uz2 þ z uz3 þ x2 uz4 þ xz uz5 þ z2 uz6 The order N of the expansion is set as an input option of the analysis; the integer N is arbitrary and it defines the order of the beam theory. Classical EBBM and TBM can be realised by using suitable F s expansions as shown in [32]. They are, in fact, degenerated cases of the linear ðN ¼ 1Þ TE model.

depend either on the order of the beam theory or on the geometry of the problem. The virtual variation of the inertial work is given by

Z

dLine ¼

Z dus

L

qF s F s dX u€ s dy ¼

Z

€ s dy dus Mss u

ð11Þ

L

X

where q denotes the material density and double over dots stand as second derivative with respect to time (t). The explicit form of the equations of motion is therefore

    26 26 22 44 36 duxs : E66 ss uxs;yy þ Es;x s  Ess;x uxs;y þ Es;x s;x þ Es;z s;z uxs  Ess uys;yy       66 26 45 45 16 þ E23 s;x s  Ess;x uys;y þ Es;x s;x þ Es;z s;z uys þ Es;z s  Ess;z uzs;y   12 q€ þ E44 xs s;z s;x þ Es;x s;z uzs ¼ Ess u     66 23 26 45 33 duys : E36 ss uxs;yy þ Es;x s  Ess;x uxs;y þ Es;x s;x þ Es;z s;z uxs  Ess uys;yy       36 66 55 55 13 þ E36 s;x s  Ess;x uys;y þ Es;x s;x þ Es;z s;z uys þ Es;z s  Ess;z uzs;y   45 q€ þ E16 ys s;x s;z þ Es;z s;x uzs ¼ Ess u       16 45 44 12 55 duzs : Es;z s  Ess;z uxs;y þ Es;x s;z þ Es;z s;x uxs þ E13 s;z s  Ess;z uys;y     16 55 45 45 þ E45 s;x s;z þ Es;z s;x uys  Ess uzs;yy þ Es;x s  Ess;x uzs;y   11 q€ þ E44 zs s;x s;x þ Es;z s;z uzs ¼ Ess u ð12Þ The generic

Eas;bh s;f ¼

Z

2.2. Equations of motions and natural boundary conditions

X

term Eas;bh s;f

above is a cross-sectional moment parameter

e ab F s; F s; dX C h f

ð13Þ

The suffix after the comma denotes the derivatives. Moreover, The principle of virtual displacements is used in this paper to derive the equations of motion.

dLint ¼

Z

dT rdV ¼ dLine

ð9Þ

V

where Lint stands for the strain energy and dLine is the work done by the inertial loadings. d stands for the virtual variation operator. The virtual variation of the strain energy can be written by using Eqs. (3), (5) and (7). After integration by part, it reads

dLint ¼

Z L

 y¼L dus K us dy þ duTs Pss us y¼0 T

ss

ð10Þ

where Kss is the differential linear stiffness matrix and Pss is the matrix of the natural boundary conditions in the form of 3  3 fundamental nuclei. They are not given here for the sake of brevity, but they can be found in [40,41]. The main property of the fundamental nuclei is that their formal mathematical expression does not

Eqss ¼

Z

qF s F s dX

ð14Þ

X

 T Letting Ps ¼ Pxs Pys Pzs to be the vector of the generalized forces, the natural boundary conditions are 26 36 66 16 duxs : Pxs ¼ E66 ss uxs;y þ Ess;x uxs þ Ess uys;y þ Ess;x uys þ Ess;z uzs 23 33 36 13 duys : Pys ¼ E36 ss uxs;y þ Ess;x uxs þ Ess uys;y þ Ess;x uys þ Ess;z uzs

duzs : Pzs ¼

E45 ss;z uxs

þ

E55 ss;z uys

þ

E55 ss uzs;y

þ

ð15Þ

E45 ss;x uzs

For a fixed approximation order N, Eqs. (12) and (15) have to be expanded using the indices s and s from 1 to M ¼ ðN þ 1ÞðN þ 2Þ=2 in order to obtain the governing differential equations and the natural boundary conditions of the desired model. In the case of harmonic motion, the solution of Eq. (12) is sought in the form

us ðy; tÞ ¼ Us ðyÞ eixt

ð16Þ

where Us ðyÞ is the amplitude function of the motion, x is an arbitrary circular or angular frequency, and i is the imaginary unit. Eq. (16) allows the formulation of the equilibrium equations and the natural boundary conditions in the frequency domain. Substituting Eq. (16) into Eq. (12), a set of three coupled ordinary differential equations is obtained which can be written in a matrix form as follows:

~s ¼ 0 dUs : Lss U

ð17Þ

where

n ~ s ¼ U xs U Fig. 2. Fiber orientation angle.

U xs;y

U xs;yy

U ys

U ys;y

U ys;yy

U zs

U zs;y

U zs;yy

oT ð18Þ

657

A. Pagani et al. / Composite Structures 118 (2014) 654–663

and Lss is the 3  9 fundamental nucleus of the matrix containing the coefficients of the ordinary differential equations. The components of matrix Lss are provided below and they are referred to as s LsðijÞ , where i is the row number ði ¼ 1; 2; 3Þ and j is the column number ðj ¼ 1; 2; . . . ; 9Þ s 44 Lsð11Þ ¼ x2 Eqss þ E22 s;x s;x þ Es;z s;z ;

ss

Lð14Þ ¼

E26 s;x s;x

þ

ss

E45 s;z s;z ;

Lð15Þ ¼

s 26 Lsð12Þ ¼ E26 s;x s  Ess;x ;

E23 ss;x



E66 ss;x ;

ss

Lð16Þ ¼

s Lsð13Þ ¼ E66 ss

E36 ss

s 44 Lsð17Þ ¼ E12 s;x s;z þ Es;z s;x ;

s 16 Lsð18Þ ¼ E45 s;z s  Ess;z ;

s Lsð19Þ ¼0

s 45 Lsð21Þ ¼ E26 s;x s;x þ Es;z s;z ;

s 23 Lsð22Þ ¼ E66 s;x s  Ess;x ;

s Lsð23Þ ¼ E36 ss

s 55 Lsð24Þ ¼ x2 Eqss þ E66 s;x s;x þ Es;z s;z ; s 45 Lsð27Þ ¼ E16 s;x s;z þ Es;z s;x ;

ss

E44 s;x s;z

ss

E45 s;x s;z

Lð31Þ ¼ Lð34Þ ¼

þ

E12 s;z s;x ;

þ

E16 s;z s;x ;

s 36 Lsð25Þ ¼ E36 s;x s  Ess;x ;

s 13 Lsð28Þ ¼ E55 s;z s  Ess;z ;

ss

E16 s;z s;x

ss

E13 s;z s

Lð32Þ ¼ Lð35Þ ¼

s 11 Lsð37Þ ¼ x2 Eqss þ E44 s;x s;x þ Es;z s;z ;





E45 ss;z ;

E55 ss;z ;

s Lsð26Þ ¼ E33 ss

ss

Lð33Þ ¼ 0 s Lsð39Þ ¼ E55 ss

For a given expansion order, N, the equations of motion can be obtained in the form of Eq. (20) as given below by expanding Lss for s ¼ 1; 2; . . . ; ðN þ 1ÞðN þ 2Þ=2 and s ¼ 1; 2; . . . ; ðN þ 1ÞðN þ 2Þ=2. It reads:

~ ¼0 LU

ð20Þ

Note that the matrix L is hierarchical. This means that the L matrix related to the N-order model actually contains the L matrix related to the ðN  1Þ-order beam model. Similarly as above, the boundary conditions of Eq. (15) can be written in a matrix form as

bs Ps ¼ Bss U

ð21Þ

where

n b s ¼ U xs U

U xs;y

U ys

U ys;y

U zs

U zs;y

oT

ð22Þ

and Bss is the 3  6 fundamental nucleus which contains the coefficients of the natural boundary conditions

2

E26 ss;x

6 23 Bss ¼ 6 4 Ess;x E45 ss;z

E66 ss

E66 ss;x

E36 ss

E16 ss;z

E36 ss

E36 ss;x

E33 ss

E13 ss;z

0

E55 ss;z

0

E45 ss;x

0

3 2 Z1 d11 6Z 7 6d 2 6 7 6 12 6 . 7¼6 . 6 . 7 6 . 4 . 5 4 . 2

Zn

s Lsð36Þ ¼0

ð19Þ

dUs :

Z;y ðyÞ ¼ SZðyÞ

3

7 0 7 5 E55 ss

ð23Þ

For a given expansion order, N, the natural boundary conditions can be obtained in the form of Eq. (24) by expanding Bss in the same way as Lss to finally give

b P¼B U

ð26Þ

Once the differential problem is described in terms of Eq. (26), the solution can be written as follows:

s Lsð29Þ ¼0

s 45 Lsð38Þ ¼ E45 s;x s  Ess;x ;

b is the expansion of U b s for a given theory order, M is the where U number of expansion terms for the given N-order beam theory, and n ¼ 6  M is the dimension of the unknown vector as well as the number of differential equations. In [40], an automatic algorithm to transform the L matrix of Eq. (20) into the matrix S of the following linear differential system was described:

ð24Þ

Matrices L and B are evaluated for each layer of the laminated beam; the matrices related to the overall laminate are then obtained by summing the contribution of each layer.

d1n

d21 d22 .. . d2n

3 C 1 e k1 y 7 6 k2 y 7 . . . dn2 76 C 2 e 7 7 6 .. .. 7 76 . 7 . . 54 .. 5 . . . dnn C n e kn y . . . dn1

32

ð27Þ

where ki is the ith eigenvalue of the S matrix, dij is the jth element of the ith eigenvector of the S matrix and C i are the integration constants which need to be determined by using the boundary conditions. The above equation can be written in matrix form as follows:

Z ¼ dCeky

ð28Þ

Vector Z does not only contain the displacements but also their first derivatives. If only the displacements are needed, according to Eq. (25), only the lines 1; 3; 5; . . . ; n  1 should be taken into account, giving a solution in the following form:

U x1 ðyÞ ¼ C 1 d11 ek1 y þ C 2 d21 ek2 y þ . . . þ C n dn1 ekn y U y1 ðyÞ ¼ C 1 d13 ek1 y þ C 2 d23 ek2 y þ . . . þ C n dn3 ekn y U z1 ðyÞ ¼ C 1 d15 ek1 y þ C 2 d25 ek2 y þ . . . þ C n dn5 ekn y .. .

ð29Þ

U zM ðyÞ ¼ C 1 d1ðn1Þ ek1 y þ C 2 d2ðn1Þ ek2 y þ . . . þ C n dnðn1Þ ekn y Once the displacements are known, the boundary conditions are obtained by substituting the solution of Eq. (28) into the boundary b is equal to Z conditions (Eq. (24)). In fact, it should be noted that U (Eq. (25)). It reads

P ¼ BdCeky ¼ KCeky

ð30Þ

where K ¼ Bd. The boundary conditions can be written in explicit form as follows:

Px1 ðyÞ ¼ C 1 K11 ek1 y þ C 2 K12 ek2 y þ . . . þ C n K1n ekn y Py1 ðyÞ ¼ C 1 K21 ek1 y þ C 2 K22 ek2 y þ . . . þ C n K2n ekn y Pz1 ðyÞ ¼ C 1 K31 ek1 y þ C 2 K32 ek2 y þ . . . þ C n K3n ekn y .. .

ð31Þ

PzM ðyÞ ¼ C 1 Kn1 ek1 y þ C 2 Kn2 ek2 y þ . . . þ C n Knn ekn y

3. Closed form solution of the equations of motion 4. Dynamic stiffness formulation Eq. (20) is a system of ordinary differential equations (ODEs) of second order in y with constant coefficients. A change of variables is used to reduce the second order system of ODEs to a first order system, T

Z ¼ fZ 1 Z 2 . . . Z n g b ¼U n ¼ U x1 U x1;y U y1 U yM;y

U zM

U zM;y

U y1;y oT

U z1

U z1;y

...

U xM

U xM;y

U yM ð25Þ

4.1. Dynamic stiffness matrix Once the closed form analytical solution of the differential equations of motion of the structural element in free vibration has been sought, a number of general boundary conditions - which are usually the nodal displacements and forces - equal to twice the number of integration constants in algebraic form needs to be applied (see Fig. 3). Starting from the displacements, the boundary conditions can be written as

658

A. Pagani et al. / Composite Structures 118 (2014) 654–663 P2z

U2z

P2y

U2y U2x

P2x

U1y P1y

By evaluating Eq. (31) at y ¼ 0 and y ¼ L and applying the BCs of Eqs. (36) and (37), the following matrix relation for the nodal forces is obtained:

y

P1x x

Fig. 3. Boundary conditions of the beam element and sign conventions.

8 P1x1 > > > > > P1 > > > y1 > > > P1z1 > > > > .. > > > . > > > < P1

At y ¼ 0 : U x1 ð0Þ ¼ U1x1 U y1 ð0Þ ¼ U1y1 ð32Þ

U z1 ð0Þ ¼ U1z1 .. . U zM ð0Þ ¼ U1zM At y ¼ L :

ð33Þ

U z1 ðLÞ ¼ U2z1 .. . U zM ðLÞ ¼ U2zM

By evaluating Eq. (29) at y ¼ 0 and y ¼ L and applying the boundary conditions of Eqs. (32) and (33), the following matrix relation for the nodal displacements is obtained:

2

d11

6 d13 6 6 6 d15 6 6 .. 6 . 6 6 6 d 1ðn1Þ zM ¼6 6 d11 ek1 L > U2x1 > > > 6 > > > > 6 > > > d e k1 L U2y1 > > > 6 > > 6 13 > > > > 6 k L > > U2z1 > > 6 d15 e 1 > > > > 6 > > . > > 6 > .. > 4 .. > > > > . > > > > : ; d1ðn1Þ ek1 L U2zM

d21

...

d23

...

d25 .. .

... .. .

d2ðn1Þ

...

d21 ek2 L d23 ek2 L

... ...

d25 ek2 L .. .

... .. .

d2ðn1Þ ek2 L

...

3 8C 9 1 > > > > > > > > C2 > > dn3 7 > > 7> > > 7> > > C3 > > dn5 7 > > > > 7> > > > 7> > .. > .. > 7 > > . . > > > 7> > > > 7> < dnðn1Þ 7 C n2 = 7 kn L 7 > C n dn1 e þ1 > > 7> > 2 > > 7> nþ2 > > > C dn3 ekn L 7 > > 2 > > 7> > > > > nþ3 > > C dn5ekn L 7 > 7> 2 > > > > 7 > .. > > 7> .. > > > > 5> . > . > > > > : ; kn L dnðn1Þ e Cn dn1

ð34Þ The above equation can be written in a more compact form as

U ¼ AC

ð35Þ

Similarly, boundary conditions for generalized nodal forces are as follows:

At y ¼ 0 : Px1 ð0Þ ¼ P1x1 Py1 ð0Þ ¼ P1y1 Pz1 ð0Þ ¼ P1z1 .. . PzM ð0Þ ¼ P1zM

2

K11

6 K 21 6 6 6 K31 6 6 .. 6 . 6 6 6 K  n1 M1 ¼6 6 K ek1 L > P2x1 > > > 6 11 > > > > 6 > > > > 6 K21 ek1 L > > P2y1 > > 6 > > > > 6 K ek1 L > > > 6 31 P2z1 > > > > > >. > 6 .. > > 6 > > > > . . > > >. > 4 > > : ; k1 L K n1 e P2M1

K12

...

K22

...

K32 .. .

... .. .

Kn2 K12 ek2 L K22 ek2 L K32 ek2 L .. . Kn2 ek2 L

... ... ... ... .. . ...

3 8C 9 1 > > > > > > > > 7 > C K2n 7 > > 2 > > > > 7> > > > C3 > > K3n 7 > > > 7> > > 7> > > .. > ... > 7 > > > > . > 7> > > > 7> < Knn 7 C n2 = 7 K1n ekn L 7 > C n2þ1 > > 7> > > > 7> kn L 7 > > C n2þ2 > K2n e 7 > > > > > > > > kn L 7 > n > C K3n e 7 > þ3 > > 2 > > > > 7 > > .. > > 7> . > > .. > 5> . > > > > > : ; kn L Knn e Cn K1n

P ¼ RC

U y1 ðLÞ ¼ U2y1

9 > > > > > > > > > > > > > > > > > > > > > =

9 > > > > > > > > > > > > > > > > > > > > > =

ð38Þ

The above equation can be written in a more compact form as

U x1 ðLÞ ¼ U2x1

8 U1x1 > > > > > U1 > > > y1 > > > > > U1z1 > > .. > > > . > > > < U1

ð37Þ

P zM ðLÞ ¼ P2zM

z

U1x

P y1 ðLÞ ¼ P2y1 P z1 ðLÞ ¼ P2z1 .. .

U1z P1z

At y ¼ L : P x1 ðLÞ ¼ P2x1

ð36Þ

ð39Þ

The integration constants vector C from Eqs. (35) and (39) can now be eliminated by relating the harmonically varying amplitudes of the generalized nodal forces to the corresponding generalized displacements to give the DS matrix of one beam element as follows:

P ¼ KU

ð40Þ

where

K ¼ RA1

ð41Þ

is the required DS matrix. The DS matrix given by Eq. (41) is the basic building block to compute the exact natural frequencies of a higher-order beam. The global DS matrix can be obtained by assembling elemental matrices as in the classical way of FEM. As far as the boundary conditions are concerned, they can be applied by using the wellknown penalty method (often used in FEM) or by simply removing rows and columns of the stiffness matrix corresponding to the degrees of freedom which are zeroes. In this paper, the penalty method is used and more details can be found in [40]. As far as the extension of the present methodology to plate-like structures is concerned, it is clear that the edges x ¼ 0 and x ¼ a have to necessarily be free, since boundary conditions can only be applied at ending nodes (edges y ¼ 0 and y ¼ b) in the case of beam models (see Fig. 1). However, this limitation can be acceptable for most of practical problems. 4.2. Application of the Wittrick–Williams algorithm For free vibration analysis of structures, FEM generally leads to a linear eigenvalue problem. By contrast, the DSM leads to a transcendental (non-linear) eigenvalue problem for which the Wittrick– Williams algorithm [57] can be used. The basic working principle of the algorithm can be briefly summarised in the following steps: (i) A trial frequency x is chosen to compute the dynamic stiffness matrix K of the final structure;

659

A. Pagani et al. / Composite Structures 118 (2014) 654–663

(ii) K is reduced to its upper triangular form by the usual form M of Gauss elimination to obtain K and the number of negaM tive terms on the leading diagonal of K is counted; this is  known as the sign count sðK Þ of the algorithm; (iii) The number, j, of natural frequencies ðxÞ of the structure which lie below the trial frequency ðx Þ is given by:

j ¼ j0 þ sðK Þ

ð42Þ

where j0 is the number of natural frequencies of all individual elements with clamped–clamped (CC) boundary conditions on their opposite sides which still lie below the trial frequency x . Note that j0 is required because the DSM allows for an infinite number of natural frequencies to be accounted for when all the nodes of the structure are fully clamped so that one or more individual elements of the structure can still vibrate on their own between the nodes. j0 corresponds to U ¼ 0 modes of Eq. (40) when P ¼ 0. Assuming that j0 is known, and sðK Þ can be obtained by M

counting the number of negative terms in K , a suitable procedure can be devised, for example the bi-section method, to bracket any natural frequency between an upper and lower bound of the trial frequency x to any desired accuracy. Once the natural frequency has been computed and the related global DS matrix evaluated, the corresponding nodal generalized displacements can be obtained by solving the associated homogeneous system of Eq. (40). By utilizing the nodal generalized displacements U, the integration constants C of the element can be computed with the help of Eq. (35). In this way, using Eq. (29), the unknown generalized displacements can be computed as a function of y. Finally, by using Eqs. (16) and (7), the complete displacement field can be generated as a function of x; y; z and the time t (if an animated plot is needed). Clearly, the plot of the required mode and required element can be visualised on a fictitious 3D mesh. By following this procedure it is possible to compute the exact mode shapes using just one element which is impossible in FEM. 5. Numerical results The capability of the present higher-order beam formulation to deal with free vibration analysis of plate-like structures is demonstrated in this section. First, the analyses of homogeneous isotropic structures are carried out for different boundary conditions. Crossply and arbitrarily laminated plates are subsequently considered. The results from the present CUF-DSM approach are compared to those from the literature and from experimental tests.

A homogeneous isotropic square plate is considered for the first assessment. Table 1 shows the first six natural frequencies for different boundary conditions. For purpose of description, in this work the same notation as used in [58] has been adopted. Therefore, the symbolism C-F-C-F identifies a plate with the edges y ¼ 0; x ¼ a; y ¼ b, and x ¼ 0 having clamped, free, clamped, free boundary conditions, respectively (see Fig. 1). Similarly SS-F-SS-F denotes a plate with two opposite simply-supported edges and two opposite free edges. Finally one single side is clamped in the C-F-F-F configuration. In Table 1 the natural frequencies are given in a non-dimensional form:

pffiffiffiffiffiffiffiffiffi

5.2. Cross-ply laminated plates Analyses of simply-supported (SS-F-SS-F) square cross-ply laminates is discussed hereinafter, even though the present formulation can deal with arbitrary stacking sequences and geometries. The lamination scheme is ð90 =0 =90 =0 =90 Þs . Each lamina has the same thickness and it is made of an orthotropic material that satisfies the following relations:

E1 G12 G13 G23 ¼ 40; ¼ ¼ 0:6; ¼ 0:5; E2 E2 E2 E2

m12 ¼ m13 ¼ m23 ¼ 0:25

In Table 2, the natural frequencies for plates with different side-to-thickness ratio (a=h) are given in form of non-dimensional frequency parameter

5.1. Homogeneous isotropic plate

x ¼ xa2 q=D

and m ¼ 0:3 is the Poisson ratio. Results from both classical EBBM and TBM and up to the sixth-order refined beam theories are presented. The results by the linear beam model ðN ¼ 1Þ are not given since they are equal to those by TBM. The solutions from the present CUF-DSM approach are compared to the results by Leissa [58] and Shu and Du [59]. The former proposed exact solutions of the characteristic equation of plates with two opposite simply-supported sides, whereas the Ritz method was employed to numerically solve problems involving arbitrary boundary conditions. On the other hand, the generalized differential quadrature method was used to solve the governing equation for thin plates in [59]. By adopting the same notation as used in [58], mode shapes are described by the number of half-waves in each direction in the present analysis case. Thus, for example, (m = 2, n = 1) denotes a mode with two half-waves in the x-direction and one half-wave in the y-direction. Fig. 4 shows the trend of the first two fundamental modes according to the classical beam models (m = 1, n = 1 and 2, respectively) versus the order of expansion of higher-order 1D CUF theories for the various boundary conditions considered. It is clear that classical and linear beam models can only detect one half-wave in the x-direction, whereas no limitations apply in the y-direction since DSM closed form solutions are adopted in this paper to solve the equations of motion along the beam axis. On the other hand, Fig. 4 shows that modes with more than one half-wave in the x-direction appear as higherorder beam models are employed. In particular, it is clear that the number of half-waves n is directly related to the expansion order of the beam displacement field N, thus a mode characterized by n = 2 requires at least a second-order ðN ¼ 2Þ beam model and so on. The first analysis case clearly demonstrates that the present CUF-DSM beam modeling approach can successfully deal with free vibration analysis of plate-like structure if sufficiently enriched displacement field are employed.

ð43Þ

where a is the side of the square plate, q is the material density, 3 D ¼ Eh =12ð1  m2 Þ; E is the Young modulus, h is the plate thickness

x ¼ x

a h

rffiffiffiffiffi

q

E2

ð44Þ

The results from both TBM and higher-order beam models are given and compared to the results by Fazzolari et al. [55], where both first order shear deformation theory (FSDT) and higher order shear deformation theory (HSDT) were used in conjunction with the DSM to carry out free vibration analysis of composite plate assemblies. The results by EBBM are not shown in Table 2, since imprecise results were produced by this theory. Figs. 5 and 6 show the first four mode shapes of the cross-ply laminae by the present sixth-order (N ¼ 6) CUF-DSM model for both a=h ¼ 5 and a=h ¼ 10. The proposed analysis show the enhanced capabilities of the present beam formulation that is able to deal with 2D-like analysis of laminate structures.

660

A. Pagani et al. / Composite Structures 118 (2014) 654–663

Table 1 pffiffiffiffiffiffiffiffiffi Non-dimensional natural frequencies, x ¼ xa2 q=D, for the isotropic plate. BCs

Model

Mode 1

Mode 2

Mode 3

Mode 4

Mode 5

Mode 6

C-F-C-F

Leissa [58] Shu and Du [59] N¼6 N¼5 N¼4 N¼3 N¼2 TBM EBBM

22.272 22.237 22.184 22.190 22.237 22.250 22.371 21.330 21.341

26.529 26.594 26.423 26.514 26.542 26.546 28.579 – –

43.664 43.871 43.807 44.020 46.852 48.305 – – –

61.466 61.407 61.153 61.168 61.346 61.397 61.614 58.752 58.821

67.549 67.659 67.153 67.482 67.537 67.555 73.554 – –

79.904 – 82.016 108.095 100.455 – – – –

SS-F-SS-F

Leissa [58] Shu and Du [59] N¼6 N¼5 N¼4 N¼3 N¼2 TBM EBBM

9.631 9.631 9.629 9.630 9.656 9.662 9.867 9.414 9.415

16.135 16.135 16.127 16.185 16.217 16.218 16.881 – –

36.726 36.726 36.916 37.140 39.123 39.160 39.452 37.638 37.654

38.726 38.945 38.926 38.933 40.319 40.868 – – –

46.738 46.739 46.704 47.033 47.090 47.096 50.677 – –

70.740 – 71.439 71.720 74.831 77.158 – – –

C-F-F-F

Leissa [58] Shu and Du [59] N¼6 N¼5 N¼4 N¼3 N¼2 TBM EBBM

3.420 3.898 3.473 3.474 3.487 3.489 3.516 3.354 3.354

8.525 9.459 8.503 8.523 8.537 8.587 8.837 – –

21.429 20.206 21.293 21.303 21.376 21.400 22.026 21.009 21.017

27.331 26.150 27.243 27.437 30.196 31.361 – – –

31.111 26.500 30.906 31.089 31.130 31.612 33.369 – –

54.443 – 54.628 54.878 58.127 59.224 – – –

Fig. 4. Order of appearance of the mode shapes versus beam theory order, isotropic plate.

Table 2 qffiffiffiffi Non-dimensional natural frequencies, x ¼ x ha Eq2 , for the cross-ply SS-F-SS-F plate. a/h

Model

Mode 1

Mode 2

Mode 3

Mode 4

5

HSDT [55] FSDT [55] N¼6 N¼5 N¼4 N¼3 N¼2 TBM

7.263 7.489 7.312 7.312 7.314 7.315 7.905 7.902

7.909 8.073 7.919 7.935 7.937 8.531 8.544 –

18.113 18.916 18.134 18.134 18.143 18.162 20.355 20.351

18.113 19.330 18.642 18.665 18.686 20.762 20.778 –

10

HSDT [55] FSDT [55] N¼6 N¼5 N¼4 N¼3 N¼2 TBM

9.939 9.516 9.544 9.545 9.547 9.547 9.861 9.855

10.248 10.352 10.376 10.394 10.395 10.733 10.755 –

29.017 28.638 28.829 29.142 29.527 34.673 – –

29.054 29.956 29.250 29.250 29.256 29.261 31.621 31.607

5.3. Symmetric 32-layer plate This analysis case consists in a 32-layer rectangular thin plate. The plate has a thickness t ¼ 6:17 mm, a width a equal to 296.5 mm and a length b ¼ 599 mm. The lamination scheme is ½ð0 =45 =  45 =90 Þ4 s and each lamina is made by an orthotropic material with the following data: E1 ¼ 157430 MPa, E2 ¼ E3 ¼ 9430 MPa, m12 ¼ m13 ¼ m13 ¼ 0:4; G12 ¼ G13 ¼ G23 ¼ 4520 MPa, q ¼ 1644:18 Kg m3 . Free edges boundary conditions are considered (F-F-F-F). The first five natural frequencies are given in Table 3, where classical and refined 1D DSM-CUF TE models are compared to the reference solutions. In particular, both a plate FEM solution from Ò MSC=NASTRAN meshed with 30  56 CQUAD4 elements and results from experimental tests are shown in Table 3 for comparison purpose. The experimental setup is shown in Fig. 7. The modal behaviour of the proposed plate is characterized by coupled torsional/bending modes on yz-plane (Modes 1 and 4),

A. Pagani et al. / Composite Structures 118 (2014) 654–663

661

Fig. 5. First four modes of the SS-F-SS-F cross-ply plate with a=h ¼ 5; N ¼ 6.

Fig. 6. First four modes of the SS-F-SS-F cross-ply plate with a=h ¼ 10; N ¼ 6.

Table 3 Natural frequencies (Hz) of the F-F-F-F symmetric 32-layer composite thin plate. Mode 1

Mode 2

Mode 3

Reference solutions 138.8 302.2 138.6 299.6

Plate FE model Experimental

112.3 110.8

N¼6 N¼5 N¼4 N¼3 N¼2 N¼1 TBM EBBM

Present 1D CUF - DSM solutions 112.4 139.7 303.6 112.5 140.1 303.8 112.5 140.5 305.6 112.6 141.3 307.7 113.4 141.6 311.2 99.8 1484.3 2968.6 99.9 – – 99.9 – –

Mode 4

Mode 5

313.4 309.9

479.3 475.3

314.3 314.4 315.3 316.9 326.9 274.7 274.7 275.3

480.4 484.5 484.8 560.4 564.8 – – –

torsional modes (Modes 2 and 3), and a coupled torsional/bending mode on xy-plane (Mode 5). The first five mode shapes from the sixth-order ðN ¼ 6Þ TE model are shown in Fig. 8. It should be noted that classical as well as linear TE ðN ¼ 1Þ beam models underestimate the natural frequencies of the 32-layer plate. In fact, these models do not foresee any torsional

Fig. 7. Experimental setup for the measurement of the natural frequencies, symmetric 32-layer plate.

and coupling effect, therefore the natural frequencies by lowerorder models that are shown in Table 3 are related to pure bending and not coupled modes. The results actually show that at least a fourth-order ðN ¼ 4Þ TE beam model is necessary to correctly detect coupled bending-torsional modes in accordance to the reference solutions.

662

A. Pagani et al. / Composite Structures 118 (2014) 654–663

Fig. 8. First five modes of the F-F-F-F symmetric 32-layer thin plate, N ¼ 6.

6. Conclusions Higher-order beam kinematics have been obtained in this paper by using the Carrera Unified Formulation (CUF), which is a hierarchical formulation that allows the automatic developments of refined beam theories by using fundamental nuclei. Closed form analytical solution has been sought and the Dynamic Stiffness Method (DSM) has been employed to produce numerical applications. Both homogeneous isotropic and composite plate-like structures have been analysed and the results of the present methodology have been successfully compared to those from the literature and from experimental tests. References [1] Tsai SW. Composites Design. 4th ed. Dayton: Think Composites; 1988. [2] Reissner E, Stavsky Y. Bending and stretching of certain types of heterogeneous aelotropic elastic plates. J Appl Mech 1961;28:402–8. [3] Timoshenko S, Woinowsky–Krieger S. Theory of plates and shells. New York: Mc-Graw Hill; 1959. [4] Reissner E. The effect of transverse shear deformation on the bending of elastic plates. J Appl Mech 1945;12(2):69–77. [5] Mindlin R. Influence of rotary inertia and shear flexural motion of isotropic, elastic plates. J Appl Mech 1951;18:31–8. [6] Reddy J. A simple higher–order theory for laminated composites. J Appl Mech 1984;51:745–52. [7] Reddy JN. Mechanics of laminated composite plates and shells. Theory and Analysis. CRC Press; 2004. [8] Carrera E. Historical review of zig-zag theories for multilayered plates and shells. Appl Mech Rev 2003;56(3):287–308. [9] Carrera E. Developments, ideas, and evaluations based upon Reissner’s mixed variational theorem in the modeling of multilayered plates and shells. Appl Mech Rev 2001;54(4):301–29. [10] Carrera E. Evaluation of layerwise mixed theories for laminated plates analysis. AIAA J 1998;36(5):830–9. [11] Khandan R, Noroozi S, Sewell P, Vinney J. The development of laminated composite plate theories: a review. J Mater Sci 2012;47(16):5901–10. [12] Euler L. De Curvis Elasticis. Lausanne and Geneva: Bousquet; 1744 (English translation: W.A. Oldfather, C.A. Elvis, D.M. Brown, Leonhard Euler’s elastic curves, Isis 20 (1933) 72–160). [13] Timoshenko SP. On the corrections for shear of the differential equation for transverse vibrations of prismatic bars. Philos Mag 1922;41:744–6. [14] Kapania K, Raciti S. Recent advances in analysis of laminated beams and plates, part I: shear effects and buckling. AIAA J 1989;27(7):923–35. [15] Kapania K, Raciti S. Recent advances in analysis of laminated beams and plates, part II: vibrations and wave propagation. AIAA J 1989;27(7):935–46. [16] Vidal P, Polit O. Vibration of multilayered beams using sinus finite elements with transverse normal stress. Compos Struct 2010;92(6):1524–34. [17] Jun L, Hongxing H. Dynamic stiffness analysis of laminated composite beams using trigonometric shear deformation theory. Compos Struct 2009;89(3):433–42. [18] Subramanian P. Dynamic analysis of laminated composite beams using higher order theories and finite elements. Compos Struct 2006;73:342–53. [19] Marur S, Kant T. Free vibration analysis of fiber reinforced composite beams using higher order theories and finite element modelling. J Sound Vibr 1996;194:337–51. [20] Kameswara R, Desai Y, Chistnis M. Free vibrations of laminated beams using mixed theory. Compos Struct 2001;52:149–60. [21] Heyliger P, Reddy J. A higher-order beam finite element for bending and vibration problems. J Sound Vibr 1988;126(2):309–26. [22] Chen WQ, Lv CF, Bian ZG. Free vibration analysis of generally laminated beams via state-space-based differential quadrature. Compos Struct 2004;63:417–25. [23] Hodges DH, Atilgan AR, Fulton MV, Rehfield LW. Free-vibration analysis of composite beams. J Am Helicopter Soc 1991;36(3):36–47.

[24] Stemple AD, Lee SW, Large deflection static and dynamic analyses of composite beams with arbitrary crosssectional warping. In: Proceedings of the AIAA/ASME/ASCE/AHS/ASC 30th Structures, Structural Dynamics, and Materials Conference, AIAA, Washington, DC; 1989, p. 1788–1798. [25] Mitra M, Gopalakrishnan S, Bhat MS. A new super convergent thin walled composite beam element for analysis of box beam structures. Int J Solids Struct 2004;41:1491–518. [26] Chandrashekhara K, Krishnamurthy K, Roy S. Free vibration of composite beams including rotary inertia and shear deformation. Compos Struct 1990;14:269–79. [27] Chandrashekhara K, Bangera KM. Free vibration of composite beams using a refined shear flexible beam element. Comput Struct 1992;43(4):719–27. [28] Carrera E. A class of two dimensional theories for multilayered plates analysis. Atti Accademia delle Scienze di Torino, Memorie Scienze Fisiche 1995;1920:49–87. [29] Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Arch Comput Methods Eng 2002;9(2): 87–140. [30] Carrera E. Theories and finite elements for multilayered plates and shells: a unified compact formulation with numerical assessment and benchmarking. Arch Comput Methods Eng 2003;10(3):216–96. [31] Carrera E, Giunta G. Refined beam theories based on a unified formulation. Int J Appl Mech 2010;2(1):117–43. http://dx.doi.org/10.1142/S1758825110000500. [32] Carrera E, Giunta G, Petrolo M. Beam Structures: Classical and Advanced Theories. John Wiley & Sons; 2011. [33] Carrera E, Petrolo M. Refined beam elements with only displacement variables and plate/shell capabilities. Meccanica 2012;47(3):537–56. http://dx.doi.org/ 10.1007/s11012-011-9466-5. [34] Carrera E, Giunta G, Nali P, Petrolo M. Refined beam elements with arbitrary cross-section geometries. Comput Struct 2010;88(5–6):283–93. http:// dx.doi.org/10.1016/j.compstruc.2009.11.002. [35] Carrera E, Petrolo M, Zappino E. Performance of CUF approach to analyze the structural behavior of slender bodies. J Struct Eng 2012;138(2):285–97. doi:10.1061/(ASCE)ST.1943-541X.0000402. [36] Carrera E, Petrolo M, Nali P. Unified formulation applied to free vibrations finite element analysis of beams with arbitrary section. Shock Vibr 2011;18(3):485–502. http://dx.doi.org/10.3233/SAV-2010-0528. [37] Carrera E, Pagani A. Multi-line enhanced beam model for the analysis of laminated composite structures. Compos Part B: Eng 2014;57:112–9. http:// dx.doi.org/10.1016/j.compositesb.2013.09.046. [38] Carrera E, Pagani A. Analysis of reinforced and thin-walled structures by multiline refined 1D/beam models. Int J Mech Sci 2013;75:278–87. http:// dx.doi.org/10.1016/j.ijmecsci.2013.07.010. [39] Giunta G, Biscani F, Belouettar S, Ferreira AJM, Carrera E. Free vibration analysis of composite beams via refined theories. Compos Part B: Eng 2013;44(1):540–52. [40] Pagani A, Boscolo M, Banerjee J, Carrera E. Exact dynamic stiffness elements based on one-dimensional higher-order theories for free vibration analysis of solid and thin-walled structures. J Sound Vibr 2013;332(23):6104–27. http:// dx.doi.org/10.1016/j.jsv.2013.06.023. [41] Pagani A, Carrera E, Boscolo M, Banerjee J. Refined dynamic stiffness elements applied to free vibration analysis of generally laminated composite beams with arbitrary boundary conditions. Compos Struct 2014;110(23):305–16. http://dx.doi.org/10.1016/j.compstruct.2013.12.010. [42] Banerjee JR. Dynamic stiffness formulation for structural elements: a general approach. Comput Struct 1997;63(1):101–3. [43] Banerjee JR. Coupled bending-torsional dynamic stiffness matrix for beam elements. Int J Numer Methods Eng 1989;28:1283–98. [44] Banerjee JR. Free vibration analysis of a twisted beam using the dynamic stiffness method. Int J Solids Struct 2001;38(38–39):6703–22. [45] Banerjee JR. Free vibration of sandwich beams using the dynamic stiffness method. Comput Struct 2003;81(18–19):1915–22. [46] Banerjee JR. Development of an exact dynamic stiffness matrix for free vibration analysis of a twisted timoshenko beam. J Sound Vibr 2004;270(1– 2):379–401. [47] Banerjee JR, Cheung CW, Morishima R, Perera M, Njuguna J. Free vibration of a three-layered sandwich beam using the dynamic stiffness method and experiment. Int J Solids Struct 2007;44(22–23):7543–63.

A. Pagani et al. / Composite Structures 118 (2014) 654–663 [48] Eisenberger M, Abramovich H, Shulepov O. Dynamic stiffness analysis of laminated beams using a first order shear deformation theory. Compos Struct 1995;31:265–71. [49] Eisenberger M, Abramovich H, Shulepov O. Vibrations of multi-span nonsymmetric composite beams. Compos Eng 1995;5(4):397–404. [50] Williams FW, Wittrick WH. An automatic computational procedure for calculating natural frequencies of skeletal structures. Int J Mech Sci 1970;12(9):781–91. [51] Wittrick WH. A unified approach to initial buckling of stiffened panels in compression. Int J Numer Methods Eng 1968;11:1067–81. [52] Wittrick WH, Williams FW. Buckling and vibration of anisotropic or isotropic plate assemblies under combined loadings. Int J Mech Sci 1974;16(4):209–39. [53] Boscolo M, Banerjee J. Dynamic stiffness formulation for composite mindlin plates for exact modal analysis of structures. part I: theory. Comput Struct 2012;96-97(0):61–73. http://dx.doi.org/10.1016/j.compstruc.2012.01.002.

663

[54] Boscolo M, Banerjee J. Dynamic stiffness formulation for composite mindlin plates for exact modal analysis of structures. part II: results and applications. Comput Struct 2012;96-97(0):74–83. http://dx.doi.org/10.1016/j.compstruc. 2012.01.003. [55] Fazzolari F, Boscolo M, Banerjee J. An exact dynamic stiffness element using a higher order shear deformation theory for free vibration analysis of composite plate assemblies. Compos Struct 2013;96:262–78. [56] Wittrick WH, Williams FW. A general algorithm for computing natural frequencies of elastic structures. Q J Mech Appl Math 1971;24(3):263–84. [57] Wittrick WH, Williams FW. A general algorithm for computing natural frequencies of elastic structures. Q J Mech Appl Sci 1970;24(3):263–84. [58] Leissa A. The free vibration of rectangular plates. J Sound Vibr 1973;31(3):257–93. [59] Shu C, Du H. A generalized approach for implementing general boundary conditions in the GDQ free vibration analysis of plates. Thin-Walled Struct 1997;34(7):837–46.